An Empirical Method of Estimating Soil Thermal Inertia

Amethod of estimating soil thermal inertia (STI), which usesmidday soil heat flux (G m ) and diurnal surface temperature amplitude as the inputs, is presented in the paper. G m is achieved from an empirical relationship between net radiation (Rn) and soil heat flux (G). To validate the STI method, a method proposed by Verhoef, which requires STI and a Fourier series analysis on surface temperature, is used to estimate diurnal G. By comparing diurnal G estimates and diurnal G measurements, the STI method is evaluated indirectly. The results show that the diurnal curve of G estimates can coincide with that of G measurements for bare soil, with the correlation coefficient (R) of 0.64, bias of 10.1W⋅m, and root mean squared errors (RMSE) of 40.9W⋅m. For the vegetated surface, R2 is 0.56, bias is −11.9W⋅m, and RMSE is 49.2W⋅m. The large uncertainty in the estimation of G m resulting from the wider variation of the empirical relationship between Rn andG and the difference betweenmixed surface temperature and soil surface temperature may be the two primary factors for the larger deviation of the diurnal shape and the magnitude between G estimates and Gmeasurements for the vegetated surface.


Introduction
Soil thermal inertia (STI, J m −2 K −1 s −1/2 ) formulated as (1) is an important surface property and is a function of specific heat capacity (, J kg −1 K −1 ), soil bulk density (, kg m −3 ), and thermal conductivity (, W m −1 K −1 ), defined in STI = √. ( It determines the resistance of a material or surface to surface temperature variations and is used to determine surface characteristics since the 1970s, such as geological interpretation, soil moisture, and soil heat flux () [1][2][3][4].
It is difficult to achieve spatial distributions of STI based on (1) because it is almost impossible to estimate , , and  at a regional scale.At present, the achievement of the spatial distributions of STI is mainly dependent on remote sensing methods.The method proposed by Price [5], shown in (2), is a milestone in the estimation of STI with remote sensing method.The primary contribution of this method is that it builds the connection between STI and the remotely sensed parameters of surface albedo (, %) and day-night surface temperature difference ( day −  night , K), which makes the retrieval of STI possible at a regional scale: where  is the solar constant (J m −2 s −1 ),  (%) is the atmospheric transmittance in the visible spectrum,  1 is the first item of the Fourier series analysis on surface temperature (),  is (rad s −1 ) the frequency of diurnal variation, and  1 (rad) is solar declination.In this method, solar radiation absorbed by the Earth's surface is taken as the driving force of the variation of .Following this method, thermal inertia has been approximated by land surface temperature differences as acquired by infrared thermometers for several decades [6].The most widely used method is the apparent thermal inertia (ATI, K −1 ) named by Short and Stuart [7].ATI method is 2 Advances in Meteorology a simplification of (2) and requires only surface albedo and diurnal surface temperature amplitude (Δ, K): where  is the solar correction factor that can be calculated from latitude and solar declination.Similar to (2), the ATI method also takes solar radiation absorbed by the Earth's surface as the driving force of the variation of .However, ATI has been found to be of very limited use in areas with strong evapotranspiration [8] as evapotranspiration can reduce surface temperature through evaporative cooling.
Additionally, there is large uncertainty when ATI is applied to vegetated surfaces because the canopy obscures the soil surface and the remote determination of soil temperature for soils under vegetation is still problematic [9].As a result, there are still some constraints in the applications of ATI.
To reduce the effects of evapotranspiration and sensible heat flux () on surface temperature, Brunt [10] indicates that it is better to use surface temperature difference between sunset and sunrise (Δ  ) than to use diurnal surface temperature when estimating STI because evapotranspiration and  are close to zero at night and have few effects on surface temperature.Equation (4) gives the function, which uses Δ  and the average net radiation (Rn  , J m −2 s −1 ) between sunset and sunrise to estimate real thermal inertia.It is only suitable for bare soil because of the difficulty in acquiring Rn  and Δ  (K) for soils under vegetation.Consider where Δ  (s) is the time between sunset and sunrise and Rn  is the average net radiation between sunset and sunrise.
Obviously, Rn  is regarded as the driving force of the decrease of  at night.There are two limitations for this method.One is that this method can only be applied to the clear and calm nights because only in this condition are evapotranspiration and  close to zero.The other limitation is that it is difficult to get Rn  and Δ  from remotely sensed data.Idso et al. [11,12] and Menenti [13] proposed a simple model to derive remote estimates of STI with Δ and the variation of  during the time between the maximum and minimum surface daily temperature.On the basis of this method, Minacapilli et al. [14] mapped STI from the airborne data and then used it to retrieve surface soil water content.So far, developing STI estimation methods is still an ongoing progress.Many methods have been proposed.Xue and Cracknell [15] proposed a real thermal inertia analytical model by taking into consideration the phase information of surface temperature changes and Δ.Murray and Verhoef [16] proposed a method to calculate STI based upon the normalized theory of soil thermal conductivity.On the basis of this method, Lu et al. [17] developed this method further by improving the Kersten function and the estimation of STI at saturation.Ma et al. [18] derived a STI method based on the heat conduction equation and an approximated energy budget equation.However, the above methods require the information on soil properties, which limit their applications of remote sensing methods to a certain degree.Using multitemporal land surface temperatures retrieved from satellite data, some other STI estimation methods were presented.Sobrino and El Kharraz [19] presented a method which uses satellite data at four moments (2:30, 7:30, 14:30, and 20:30 in local time).Cai [20] modified Sobrino's method for MODIS (Moderate Resolution Imaging Spectroradiometer) data.Nearing et al. [21] developed a method with two daily measurements of surface temperature.The limitation of these methods is that they are restricted to the availability of satellite overpasses or field measured data [22].Geostationary satellite data may provide a good alternative input for these methods.
According to the soil conductive heat transfer, soil heat flux is the real driving force of  [23,24].However, because  is difficult to acquire accurately at a regional scale by remote sensing method,  is usually replaced by  or Rn when estimating STI, as used in (3) and (4).When the alternative driving energy is close to , the estimated STI would be close to the real STI [25].Otherwise, the estimated STI would deviate greatly.
As is well known,  has a close relationship with Rn.Many studies have formulated the empirical function between them ((, Rn)) [26][27][28].Equations ( 5) and ( 6) are the two typical methods widely used to estimate  for bare soil (see (5)) and for vegetated surface (see (6)) at a regional scale, respectively: where VI is the vegetation index; , , and  are the empirical coefficients, which are mainly determined by soil properties, soil moisture, , and LE [29,30].The objective of this paper is to present a STI estimation method which uses midday  (  ) and Δ as the inputs.Here,   is acquired from the empirical relationship between  and Rn.Another purpose of this paper is to investigate whether the empirical function about , respectively, for bare soil and for vegetated surfaces can be used to estimate soil thermal inertia in combination with Δ.Due to the large uncertainties in the estimation of STI for vegetated surface, STI method is limited in a variety of applications, such as soil moisture retrieval.If the method proposed in the paper is applicable, it will provide a feasible way to estimate STI for vegetated surface and thereby expand the applications of STI.The major difference between the method proposed in this paper and the previous methods is that the combination of   is used to estimate STI, which could be better than Rn or solar radiation according to the soil conductive heat transfer.The relationship between  and Rn is a simple method for acquiring   .For an area with relatively uniform soil properties and meteorological conditions (namely, there is a reliable relationship between  and Rn), STI could be estimated conveniently using the proposed method.
In the paper, a method of estimating soil thermal inertia, which uses   and Δ as the inputs, is described in Section 2. Study site, data description, and the evaluation method are

Site Description.
Measurement data used in the study are from the Yucheng (YC) Agro-Ecosystem Station, Chinese Academy of Sciences (CAS).The station is located in an irrigated agricultural field (36.8570N, 116.8360E) in the North China Plain characterized by a semihumid and monsoon climate.Mean annual precipitation, temperature, and global solar radiation at the station are 528 mm, 13.1 ∘ C, and 5225 MJ m −2 , respectively.The dominant soil type is sandy loam with an average bulk density of 1.43 (g cm −3 ), soil wilting point of 0.07 (v/v), and saturated soil moisture of 0.462 (v/v).The groundwater table varies from 1.5 to 3.5 m with an average of 2.5 m.This site is a typical representative of the agricultural area of winter wheat and summer maize production in the North China Plain.According to the traditional tillage practice, winter wheat is sown in middle October and harvested in early or mid-June of the next year, and summer maize is planted in mid-to late June and harvested at the end of September.It is almost a bare soil surface in the nongrowing season (Jan, Feb, mid-and late Oct, Nov, and Dec).In the growing season, winter wheat and summer maize are the dominant vegetation types, respectively, in the periods from March to June and from July to September.

Soil Thermal Inertia Method.
Similar to the concepts of STI  and ATI, the method of estimating STI used in this study is given in where   (J m −2 s −1 ) is the midday soil heat flux.Δ (s) is the time difference between daily maximum and minimum surface temperature.The difference between (7), STI  , and ATI is that (7) uses   rather than Rn or solar radiation absorbed by the surface as the driving force of the variation of .Because it is easier to get midday soil heat flux from remotely sensed data, ( 7) can be applied to remote sensing data.
In the study,   is determined by the empirical function between Rn and , like ( 5) and (6).The relationship between Rn and  is established on the basis of in situ measurements of Rn and  in 2003 at YC Station.Under the condition that there is only a small annual variation of surface conditions at YC Station in a few years, the relationship between Rn and  can be assumed to be almost constant; namely, the empirical coefficients of , , and  in ( 5) and ( 6) are constant.In this case,   in other years can be estimated on the basis of the established function of Rn and .With the midday Rn (Rn  ) as the input,   of 2004 was estimated and used to estimate STI of 2004 in the study.Note that Rn  , Δ, and Δ all can be estimated from remotely sensed data, for example, MODIS.Therefore, as long as the empirical function between Rn and  is determined, STI can be estimated by (7).

Diurnal Soil Heat Flux Estimation.
Because there are no data of in situ soil thermal inertia measurements, an indirect validation for the STI method was used.Firstly, diurnal  in 2004 at YC Station was estimated based on the method proposed by Verhoef [9].By comparing diurnal  estimates with  measurements, the STI method was then evaluated indirectly.
Based on the method of Horton [4], Verhoef [25] presented a remote sensing method of estimating diurnal  for bare soil, as shown in (8).This approach requires a Fourier series analysis on surface temperature first: where  (s) is the time during the 24 h cycle,  is the total number of harmonics,  is the harmonic number, and  (rad s −1 ) is the radial frequency used in the harmonic analysis.  ( ∘ C) is the overall amplitude of the th harmonic and   (rad) is the phase shift of the th harmonic, which are calculated from where   ( ∘ C) and   ( ∘ C) are the amplitudes in the harmonic representation of soil surface temperature and can be obtained from the Fourier series analysis on .In this way,  for each 24 h period was represented by a sum of sine and cosine terms (harmonics) according to the following equation: where the overbar denotes the average of  over 24 h ( ∘ C).
Other parameters have the same meanings as in (8).It can be seen from ( 8) that the summation of harmonic terms determines the diurnal shape of  estimates and STI determines the magnitude of  estimates.With ( 8) and ( 10), diurnal  in 2004 is estimated based on the measurements of  and STI estimated from (7).Diurnal  measurements in 2004 are then compared with the measurements of  of 2004 at YC Station.
where  is the surface emissivity, taken as 0.95 for bare soil and 0.98 for vegetation surface, and  (W m −2 K −4 ) is the Stefan-Boltzmann constant.

Advances in Meteorology
Soil heat flux at a depth of 0.05 m was measured by a soil heat flow plate (HFP01SC, Hukseflux).Soil heat flux at the surface was estimated using the following formula: where   (J m −2 s −1 ) is the plate heat flow measurement, Δ  (K) is the change in mean soil temperature during the measurement period,   (J m −3 K) is the volumetric heat capacity of the soil,  (m) is the depth of the plate, and  (s) is the length of the measurement period.The standard equation given by De Vries [31] was used to estimate the volumetric heat capacity,   : where   and   are the volumetric heat capacity of water and solid soil minerals, respectively (i.e., 4.2 and 2.0 MJ m −3 K −1 ) and  (m 3 m −3 ) is the actual soil volumetric water content and was measured with a water content reflectometer (Model CS615-L, Compbell Scientific) at a depth of 20 cm. * is the saturated soil moisture content.Data in the rainy days were not considered in the study due to the large variation of soil moisture, which would have great effect on the determination of the empirical coefficients between Rn and .Rainfall was monitored with a rain gauge (52203, RM Young, Inc.).
Three steps were used to perform the data quality control: (1) checking for null values for ,   , and Rn data for the same day, which is mainly caused by the broken electricity problem (if one null value was found, the data for the day were eliminated); (2) eliminating data on rainy and snowy days based on the rain and snow observations; (3) checking whether there were diurnal variations for ,   , and Rn for the same day.Only data with obvious diurnal variations for all the three parameters were kept.After this screening, 96 days and 128 days of observations for bare soil and vegetated surfaces, respectively, in 2003 were used to establish the regression function between Rn and .In 2004, 100 days and 115 days of observations, respectively, for bare soil and for vegetated surface were used to validate the method.
The data of VI at YC Station are from the MODIS 16-day vegetation indices product (MOD13Q1) acquired from the Land Process Distributed Active Archive Center (LPDAAC) (https://lpdaac.usgs.gov/).Its spatial resolution is 250 m.A linear interpolation method is used to retrieve the Normalized Difference Vegetation Index (NDVI) of each day.

The Determination of the Coefficients of 𝑎, 𝑏, and 𝑐.
In the study, the measurements of Rn and  in 2003 at YC Station in nongrowing season (Jan, Feb, mid-and late Oct, Nov, and Dec) were used to determine the empirical coefficient of  in (5), which is realized in two steps: First,  for each day is determined by the linear regression between the measurements of Rn and , which is the slope of the regression.Then, the average of  is calculated by averaging  of each day and is used to calculate  of bare soil in 2004.
A regression function between Normalized Difference Vegetation Index (NDVI) and the ratio of  and Rn was established for the growing season using the field measurements of  and Rn in 2003, from which  and  in (6) are determined.Figure 1(a) shows the scatter plot of Rn and  in the nongrowing season.Figure 1(b) shows the scatter plot of NDVI and the ratio of  and Rn in the growing season.
Using the above methods, the regression functions between Rn and  were established, shown in (14).The coefficients of , , and  are −0.413,0.457, and 0.472, respectively.They were used to estimate   at YC Station in 2004, with Rn  as the input:

Results and Discussions
Figure 2 shows the scatter plot between   estimates calculated from ( 14) and   measurements at YC Station in 2004.It shows that their values are close to 1 : 1 line as a whole, with a bias of −6.5 W⋅m −2 and a RMSE (root mean squared error) of 30.7 W⋅m −2 .Integrating   estimates and Δ determined from  measurements into (7), STI was estimated.Figure 3 shows the annual variations of STI in 2004.It can be seen that STI ranges from 500 to 2500 (J m −2 K −1 s −1/2 ) and higher STI was found in the growing season, with average STI of 1061.2 (J m −2 K −1 s −1/2 ) in the nongrowing season and of 1370.6 (J m −2 K −1 s −1/2 ) in the growing season.Rainfall and frequent irrigation may be the primary reason for the higher STI in the growing season because high soil water content can usually produce high soil thermal inertia [17].Integrating STI obtained above into (8) and on the basis of the Fourier series analysis on the surface temperature, diurnal  in 2004 was estimated.The comparisons between diurnal  estimates and  measurements during the growing season and during the nongrowing season are given in Figure 4.The scatter plots are given in Figure 5.
As a whole,  estimates can fit with  measurements, with correlation coefficient ( 2 ) of 0.64 and 0.56 for the nongrowing season and the growing season, respectively.The biases and the RMSE for the nongrowing season and for the growing season are 10.1 W⋅m −2 , 40.9 W⋅m −2 and −11.9 W⋅m −2 , 49.2 W⋅m −2 , respectively. estimates for the bare soil are slightly better than that for the vegetated surface.The biases show that  is slightly overestimated in the nongrowing season and is underestimated in the growing season.It was also seen from Figure 4 that the discrepancies between  estimates and  measurements in the nighttime are larger than that in the daytime, which suggests that   estimates used in (8) deviate more greatly from the nighttime soil heat flux in comparison with the daytime soil heat flux.
As mentioned in Section 2.3, the diurnal variation of  determines the shape of the diurnal curve of  estimates by determining the summation of harmonic terms in (8) and STI determines the magnitude.Therefore, the larger errors in the STI estimates for the vegetated surface may be the main reason for the larger errors in  estimates.As described    in Section 2,   for vegetated surface is obtained from the empirical function between NDVI and the ratio of  and Rn (see ( 14)).Vegetation conditions, irrigation activity, radiation conditions, and soil properties all influence the relationship between NDVI and the ratio of  and Rn.Therefore, there are wide variations in the empirical coefficients of the function, which leads to the uncertainties in the   estimation for the vegetated surface and thereby produces the errors in the STI estimation.Comparatively, the relationship between  and Rn for bare soil is relatively stable because soil properties play a dominant role in their relationship, which is relatively invariable in the same place.It seldom rains at YC Station and there is no irrigation activity in the nongrowing season, so soil moisture is relatively stable, which reduces the effects of the variation of soil moisture on the relationship between  and Rn.
To see the comparisons of  estimates and  measurements clearly, Figure 6 displays the 5-day results.
It is obvious that the diurnal curve of  estimates coincides with the diurnal shape of , which illustrates that the diurnal variation of  determines the shape of the diurnal curve of  estimates by determining the summation of harmonic terms in (8).Compared with the in situ ,  estimates have a similar diurnal shape in the nongrowing season, while in the growing season, the fitness of the diurnal shape is worse.As mentioned above, the errors in the estimation of   are a primary reason for the discrepancy between  estimates and  measurements.For the vegetated surface, there is another reason.Surface temperature used to calculate diurnal  at vegetated surface is a mixed temperature which is composed of soil surface temperature and canopy temperature.However, the fact is that  only has relationship with soil temperature.For vegetated surface, soil temperature for soils under vegetation is usually not available from the conventional observation and remotely sensed data.Therefore, the difference between mixed surface temperature and soil temperature would produce the difference in the diurnal shape of  estimates and  measurements in the growing season.An underestimation of  was found during nighttime periods ( estimates are lower than  measurements at nighttime).This underestimation arises from differences between canopy and soil temperatures, where temperatures at the top of the canopy fall below those of      the soil due to the effect of the sheltering of the soil surface by the canopy coverage.In the study, canopy temperature calculated from the radiative measurements by a four-component radiometer was used to compute the summation of harmonic terms in (10).Thus, the influence of canopy temperature on the radiation measurements induces a lower estimate of the summation of harmonic terms at nighttime.The same phenomenon was found in the study of Murray and Verhoef [9].They attributed this to the differences between canopy and soil temperatures.The above results of the comparison between diurnal  and in situ  illustrate that (8) can be used to estimate soil thermal inertia correctly.However, there are larger uncertainties in the estimation of STI at vegetated surfaces than that at bare soil due to the larger uncertainties in the estimation of   caused by the wider variation of the relationship between Rn and .Using mixed surface temperature to replace soil temperature is another uncertainty in the estimation of STI for the vegetated surface.

Advances in Meteorology
The method can be applied to map STI with remotely sensed data.As is known, Rn can be retrieved from remotely sensed albedo, surface temperature, and some auxiliary meteorological data, such as solar radiation and air temperature.On the basis of the regression function between Rn and  in the study area,   can be calculated according to Rn. Δ in (7) can also be obtained from satellite data, such as MODIS.Therefore, STI at a regional scale can be estimated.By using the regression function between Rn and  for vegetated surfaces, STI at vegetated surface can be retrieved, which provide an opportunity to overcome the difficulty in acquiring STI at vegetated surface.The method is simpler than the physically based STI methods that usually require soil properties.It is difficult to acquire soil properties by remote sensing; as a result, STI at a regional scale usually is not dependent on the physically based STI methods.
The limitation of the method presented in the paper is that the estimation of   depends on the empirical relationship between Rn and .Actually, the relationship between Rn and  is not constant.It varies with soil properties, meteorological condition, and vegetation condition.Therefore, if the surface condition changes, errors will be produced in the estimation of   and thereby in the estimation of STI.In the application of the method, it is recommended that the regression function between Rn and  for the study area is established first, which would make   estimation more accurate and thereby make the STI estimation more reliable.The investigation of the accuracy of the  estimation method at a regional scale will be one of the future researches, which is expected to reduce the limitation of the empirical method of the determination of   .Applying the method to retrieve soil water content for both bare soil surface and vegetated surface is another implication of this research.

Conclusions
On the basis of the fact that soil thermal inertia (STI) can be estimated by the combinations of the driving force of the variation of surface temperature and the diurnal amplitude of surface temperature, a STI method, using midday soil heat flux (  ) and Δ as the inputs, is presented in the paper.  is achieved from an empirical method.
In the study, the measurements of net radiation and  in 2003 at the Yucheng (YC) Agro-Ecosystem Station, Chinese Academy of Sciences, were used to establish the empirical functions between Rn and , respectively, for bare soil (nongrowing season) and for vegetated surface (growing season).According to the empirical functions,   of 2004 at YC Station is estimated and thereby STI of 2004 is estimated.
To validate the STI method, diurnal  of 2004 at YC Station is estimated with the method proposed by Verhoef [25], which requires STI and a Fourier series analysis on surface temperature.The comparisons of diurnal  estimates and the diurnal  measurements show that (1)  estimates agree well with  measurements, with correlation coefficient ( 2 ) of 0.64 and 0.56 for the nongrowing season and the growing season, respectively.The biases and the root mean squared errors (RMSE) for the nongrowing season and for the growing season are 10.1 W⋅m −2 , 40.9 W⋅m −2 and −11.9 W⋅m −2 , 49.2 W⋅m −2 , respectively.(2) The diurnal curve of  estimates coincides well with the diurnal curve of  measurements for bare soil.For the vegetated surface, the consistence between them is poorer than that for bare soil.The large uncertainties in the estimation of   resulting from the wider variation of the empirical relationship between Rn and  and the difference between mixed surface temperature and soil temperature may be the two primary reasons for the larger discrepancies in the diurnal shape and in the magnitude of  estimates and  measurements in the growing season.(3) The performance of  estimates in the daytime is better than that in the nighttime.
The comparison results between  estimates and in situ  illustrated that the STI method presented in the paper can be used to estimate soil thermal inertia reliably and the method can estimate STI more accurately for bare soil than for vegetated surface.

Figure 1 :Figure 2 :
Figure 1: The linear relationship between the ratio of  to Rn and NDVI at YC Station in 2003.

Figure 4 :
Figure 4: The diurnal course of  estimates and  measurements for the nongrowing season (a) and for the growing season (b).

Figure 5 :
Figure 5: The scatter plots between  estimates and  measurements for the nongrowing season (a) and for the growing season (b).

Figure 6 :
Figure 6: The diurnal course of  estimates,  measurements, and surface temperatures of 5 days.