Comparison of six algorithms to determine the soil thermal diffusivity at a site in the Loess Plateau of China

Soil thermal diffusivity is a crucial physical parameter that affects soil temperature. Six prevalent algorithms to calculate soil thermal diffusivity are inter-compared by using soil temperature data collected at the depths of 0.05 m and 0.10 m at a bare site in the China Loess Plateau from DOY 201 through DOY 207 in 2005. Five of the 5 six algorithms (i.e., Amplitude, Phase, Arctangent, Logarithm, and Harmonic or HM algorithms) are developed from the traditional one-dimensional heat conduction equation. The other algorithm is based on the one-dimensional heat conduction-convection equation which considers the vertical heterogeneity of thermal diffusivity in soil and couples thermal conduction and convection processes (hereinafter referred to as the 10 Conduction-convection algorithm). To assess these six algorithms, we (1) calculate the soil thermal diffusivities by using each of the algorithms, (2) use the soil thermal diffusivities to predict soil temperature at the 0.10 m depth, and (3) compare the estimated soil temperature against direct measurements. Results show that (1) HM algorithm gives the most reliable estimates among the traditional five algorithms; and (2) gener15 ally, the Conduction-convection algorithm provides the second best estimates. Among all of the algorithms, the HM algorithm has the best description of the upper boundary temperature with time, but it only includes conduction heat transfer in the soil. Compared to the HM algorithm, the Conduction-convection algorithm has a less accurate description of the upper boundary temperature, but by accounting for the vertical gra20 dient of soil diffusivity and the water flux density it includes more physics in the soil heat transfer process. The Conduction-convection algorithm has potential application within land surface models, but future effort should be made to combine the HM and Conduction-convection algorithms in order to make use of the advantages of each.


Introduction
Soil temperature plays an important role in land surface processes, and it is critical in energy balance applications such as land surface modeling, numerical weather forecasting, and climate prediction (Holmes et al., 2008).It is especially true for the soil surface.Accurate prediction of soil surface temperature requires a realistic understanding of the soil thermal properties, i.e., thermal conductivity (λ), thermal diffusivity (k), and volumetric heat capacity (C g ).The volumetric heat capacity C g can be estimated as follows (De Vries, 1963), where η is the volumetric water content, η s is the saturated value of η, and C s and C w are the volumetric heat capacities of dry soil and water respectively.If C g is known, only thermal conductivity, λ, or thermal diffusivity, k, must be determined to characterize the thermal properties of a soil (Passerat et al., 1996).k is of primary importance in determining soil temperature propagation (Zhang and Osterkamp, 1995).
Several algorithms have been proposed to estimate soil thermal diffusivity.Most of the algorithms are based on solutions of the one-dimensional conduction heat transfer equation.Lettau (1971) calculated the thermal diffusivity as a function of depth below the soil surface.In order to utilize this algorithm, measurements of soil temperature with time are required at the soil surface and at several subsurface depths.However, the lack of soil temperature data at several subsurface depths often limits the utility of this algorithm (Horton et al., 1983).Assuming that the thermal diffusivity is independent of depth, and considering that temperature at the upper boundary is well described by a sinusoidal function, the analytical solution of the one-dimensional heat conduction equation can be used to estimate k.Based on this solution, the thermal diffusivity can be estimated by the Amplitude algorithm and Phase algorithm.Errors due to the assumption of single sinusoidal temperature wave at the soil surface can be reduced by using a Fourier series to accurately describe the diurnal variation in surface soil temperature (Van Wijk, 1963).In this way, the thermal diffusivity is estimated by 2249 the Arctangent algorithm (Nerpin and Chudnovskii, 1967) with two harmonics.It was shown by Seemann (1979) that, in analogy to the Arctangent algorithm, the thermal diffusivity can also be calculated by the Logarithmic algorithm.These two algorithms are analogous to the Amplitude algorithm and Phase algorithm but take advantage of greater number of temperature observations to approximate a potentially nonsinusoidal behavior (Horton et al., 1983).However, a two-harmonic function cannot describe the surface temperature very well, and a series of harmonics for the upper boundary offers advantages.Based on this boundary, the solution of the one-dimensional heat conduction equation is developed from the assumption of a sinusoidal function.According to the solution, the thermal diffusivity can be selected to minimize the sum of squared differences between the calculated and measured soil temperature values (Horton et al., 1983).And it can also be estimated from an iteration process by fitting the amplitude and phase of soil temperature at one depth (Heusinkveld et al., 2004).All algorithms mentioned above are based on solutions of the one-dimensional conduction heat equation and constant diffusivity, and thus apply to uniform soils only.In fact, soil heat transfer is caused by a combination of conduction and intra-porous convection (Passerat et al., 1996).Gao et al. (2003) pointed out that soil temperature changes in response to both conduction and convection processes, where convection was understood as "vertical heat transfer caused by the vertical movement of liquid water in the soil".They solved analytically the equation for one-dimensional conduction-convection, and derived a simple algorithm to accurately estimate soil thermal diffusivity.Few efforts have been made to quantitatively test the various k algorithms by using an identical soil temperature data set.The land-air interaction over the Loess Plateau located in mid-western China affects the weather and climate in northwest China.A realistic description of soil temperature helps better understanding the land-air interaction over the Loess Plateau however few attempts have been made to determine the Loess Plateau soil thermal parameters for soil temperature algorithms.The measurements of soil temperature and soil water content during the LOess Plateau land surface process field EXperiment (LOPEX) in 2005 provided us an opportunity to evaluate the k algo-rithms for use on Loess Plateau soil.In order to improve the accurate knowledge of the soil thermal diffusivity in this area, the objective of this paper is to compare six algorithms for determination of the soil thermal diffusivity by using the direct measurements of soil temperature restricted to the upper 0.1 m of soil.
2.1 Classical thermal conduction equation for soil temperature Conduction heat transfer in a one-dimensional isotropic medium is described by where T is the soil temperature (K), t the time (s), z the depth (m), C g the volumetric heat capacity , and λ the thermal conductivity (W m −1 K −1 ).Assuming that a soil is vertically homogeneous, in this case that C g and λ are independent of depth, provides where the thermal diffusivity k=λ/C g (m 2 s −1 ).The following five algorithms based on the solution of Eq. ( 3) have been used to estimate k.

Amplitude algorithm
Given the surface boundary condition: where T is the mean soil surface temperature, A is the amplitude of the diurnal soil surface temperature wave, and ω is the angular velocity of the Earth's rotation and ω=2π/P (rad s −1 ) with P representing the period of the diurnal cycle.The soil temperature (T ) at a depth z can be calculated via here d = 2k/ω is the damping depth of the diurnal temperature wave.Soil temperature measured at two different depths (z 1 and z 2 ) are often assumed to be approximated by a sinusoidal function when estimating k.The sinusoidal functions are given by and where A 1 (A 2 ), Φ 1 (Φ 2 ) and T 1 (T 2 ) are the amplitude, phase and mean soil temperature at the depth z 1 (z 2 ).T 1 (T 2 ) is the arithmetical average of the daytime maximum soil temperature and the nighttime minimum soil temperature; and A 1 (A 2 ) is half of the difference between the daytime maximum value and the nighttime minimum value for soil depth of z 1 (z 2 ); and Φ 1 (Φ 2 ) is the initial phase of soil temperature at depth z 1 (z 2 ), obtained by using the best fit algorithm (Horton et al., 1983).Then the thermal diffusivity k is determined by the Amplitude algorithm

Phase algorithm
If the time interval between the measured occurrences of maximum soil temperature at the depths of z 1 and z 2 is ∆t=t 2 −t 1 , the Phase algorithm stemming from Eq. ( 4) is (Horton et al., 1983)

Arctangent algorithm
Soil surface temperature can be described by a Fourier series: where n is the number of harmonics, and a i and b i are the amplitudes.With boundary condition n=2, k can be calculated by the Arctangent algorithm where temperatures T j and T j are recorded each 6 h (j =1, 2, 3, and 4) at two different depths z 1 and z 2 , respectively.The first reading is taken at 02:00 (Local time, hereinafter, referred to as LT), then 08:00 (LT), 14:00 (LT), 20:00 (LT).

Logarithmic algorithm
Using the same assumption of the Arctangent algorithm, k is expressed by

Harmonic algorithm
Equation ( 10) can also be changed into another form as where C i is the amplitude of the harmonic i : and Φ i is the phase of the harmonic i : Given the following boundary condition: the solution of Eq. ( 3) is where C 0i and Φ 0i are the amplitude and phase of the harmonic i for the upper depth, respectively, and d i = 2k/(i ω) corresponds to the depth at which the signal is propagated during a period P/i .Based on Eq. ( 17), the thermal diffusivity k can be determined by the Least Squares Algorithm (Horton et al., 1983).On the other hand,C 1i (C 2i ) and Φ 1i (Φ 2i ) at the depth of z 1 (z 2 ) can be obtained by the approximation of the observed data at these two depths with the harmonic curve fit.In addition, according to Eq. ( 17), the amplitude C 2i and initial phase Φ 2i at the depth of z 2 can be predicted from and After an initial guess of k, the predicted results of amplitude and initial phase are compared with the fitted ones, and the parameter is adjusted depending on the differences in amplitude and initial phase (Heusinkveld et al., 2004).
2.2 Soil temperature rate equation with vertical heterogeneity of soil thermal diffusivity coupled with thermal conduction and heat transfer by water flux Equation (3) assumes that k is independent of depth, however, k can vary (increase or decrease) from the surface downward in the shallow surface layer of most soils.Equation (2) can therefore be improved as follows (Gao et al., 2008): Neglecting the vertical heterogeneity of k, Gao et al. (2003) incorporated thermal conduction and convection together as follows: where w is the liquid flow rate (positive downward) and θ is the volumetric water content of the soil.C w is the heat capacity of water.Assuming these three parameters are also independent of z for a thin soil layer in present work.− C W C g wθ was defined as water flux by Gao et al. (2003).Based on Eqs. ( 19) and ( 21), Gao et al. (2008) presented the following equation With the boundary condition Eq. ( 7), the expression of the soil temperature at the depth z 1 is where They also gave the expression of k and We call it as the Conduction-convection algorithm in this paper.Applying W =0 to Eq. ( 25) results in Φ 2 −Φ 1 =− ln(A 2 /A 1 ) or− ln(A 2 /A 1 )=Φ 2 −Φ 1 , then Eq. ( 24) reduces to be Eq.( 8) or Eq. ( 9).

Field experiments
The experiment was conducted on soil in the China Loess Plateau during an intensive observation period from DOY 197 through 241 in LOPEX in 2005.The soil measurements were collected at a bare soil site located at 106.42The ground surface of this site was bare, flat and homogeneous.The soil at the site was predominantly medium loam with a high proportion of silt.The site is located within a semiarid climate zone.The maximum air temperature was 307 K and the lowest was 249 K, the annual air temperature and precipitation were 279 K and 510 mm with 2425 h of sunshine, and 170 frost-free days per year all averaged over the last 50 years (Gao et al., 2008).
Soil temperature was measured with four TCAV averaging soil thermocouple probes (Campbell Scientific Inc., USA) at 0.05 and 0.10 m depths.The volumetric water content of the soil was measured at 0.05 and 0.10 m depths by two soil moisture reflectometers (CS615, Campbell Scientific Inc. USA).All of the sensor outputs were recorded and averaged over 10 min intervals.

Results
The data used in this paper were collected in a 7-day period from (DOYs 201 through 207 (i.e.,20  for the soil layer from 0.05 to 0.10 m depths at 14:45 (LT) on DOY 204 at this site.

2257
The temporal variations in volumetric soil water content at 0.05 and 0.10 m depths during the same period are shown in Fig. 1b.Precipitation occurred from DOYs 199 to 201 with an amount of 15.7 mm.Since then, owing to evaporation from the bare soil surface, the soil volumetric water content was decreasing gradually at both depths from DOYs 202 through 207.Gao et al. (2008) pointed out that under evaporation conditions there is a net upward flux of water (liquid and vapor) that responds to the progressively drying surface condition.The net flux of water causes an associated net convective heat flux.The soil physics implies that the heat transfer should incorporate a vertical convective heat transfer component.Figure 1b also shows that the soil volumetric water content changed diurnally and soil volumetric water content at the shallow level was lower than that at the deeper level in rain-free days.
After using a 2-h smoothing technique for soil temperature measured at the depth of 0.05 m, A 1 , A 2 , Φ 1 , Φ 2 , T 1 , T 2 , T , C i and Φ i are obtained by using the approximation of soil temperature collected at the depths of 0.05 m and 0.10 m for each day, respectively.The temporal variations of k are calculated by using the Amplitude, Phase, Conduction-convection and HM algorithms (Horton et al., 1983;Heusinkveld et al., 2004;Gao et al., 2008) for the soil layer from 0.05 m to 0.10 m.Temperatures at each depth for the arbitrary times of 02:00 (LT), 08:00 (LT), 14:00 (LT), and 20:00 (LT) were used for Arctangent algorithm and Logarithmic algorithm.Results are shown in Fig. 2. The Amplitude, Arctangent, Logarithmic and HM algorithms provided relatively low values of k during this period.This is especially true for the Arctangent algorithm.The Phase algorithm and the Conduction-convection algorithm provided k values approximately twice as large as the other algorithms.The two HM algorithms provided similar values of k.The thermal diffusivity estimated by the Logarithmic algorithm changed from day to day in the drying period.The maximum, minimum and mean values of k calculated by these six algorithms are listed in

Discussion
The variations of the thermal diffusivity k obtained by the six algorithms with the volumetric soil water content θ for the 0.05 m to 0.10 m layer are shown in Fig. 3. Estimates of k from five of the six algorithms change in a narrow range with θ during this drying period (θ<25%).The exception is the Logarithmic algorithm.The values of k from the Phase algorithm and the Conduction-convection algorithm have similar trends with θ.The variations of k shown by the other four algorithms show a similar trend with θ.
These results make sense because values of k from the Phase and the Conductionconvection algorithms mainly depend on Φ 1 −Φ 2 , while estimates of k with the other algorithms mainly depend on the amplitude.All of the algorithms indicate that the largest value of k does not occur at the largest soil water content.Earlier researchers have reported that k does not monotonically increase with increasing θ.It tends to increase as dry soil begins to wet, but it approaches a constant value or even decreases as the soil continues to wet.
In the Conduction-convection algorithm, another parameter W is needed and is calculated with Eq. ( 25) (see values in Table 2).
The measured and modeled soil temperatures at the 0.10 m depth from DOY 201 through 207 are presented in Fig. 4. It is obvious that the fitted temperature by using Arctangent, Logarithmic and the two HM algorithms gave poor approximation in the early morning of DOY 203 because the assumption of repeating surface periodic temperature is not valid between DOY 202 and 203.To better show the model outputs, we take DOY 204 as an example (see Fig. 5).Overall, The Phase algorithm reasonably estimated the soil temperature phases but overestimated the amplitudes, and the amplitude algorithm reasonably estimated the soil temperature amplitudes but overestimated 2259 the phase shift.In fact, as shown in Table 3, the values of ln(A 1 A 2 ) are larger than Φ 1 −Φ 2 for the whole 7-day period.Using the Phase algorithm to estimate k implies forcing ln(A 1 A 2 ) to be equal to Φ 1 −Φ 2 , which overestimates the soil temperature amplitude by about 0.74 K on average for DOYs 201 to 207.Similarly, using the Amplitude algorithm to estimate k implies that Φ 1 −Φ 2 is equal to ln(A 1 A 2 ), which overestimates the soil temperature phase shift by ln(A 1 A 2 )−(Φ 1 −Φ 2 )=0.2397 rad (55 min) on average for DOYs 201 to 207.The Logarithmic and Arctangent algorithms require four pairs of soil temperature measurements.The modeled k values are very sensitive to the measurement time of four pairs of soil temperatures, so we have to average the calculated values of k for different selections of four pairs of soil temperatures for each day.
The two HM algorithms generated similar values of k.For most of the study days, the HM algorithm gave realistic estimations of soil temperature at the depth of 0.10 m.However, as mentioned above, it did not give a good estimation in the early morning of DOY 203 because of the invalid assumption of the repeating periodic for soil temperature.Results obtained by two HM algorithms indicate that fitting measurements of soil temperature by using the Least Squares approach directly and simultaneously determining the amplitude and initial phase of the soil temperature may provide realistic values of k.Comparison of the results obtained by the HM algorithms, Phase algorithm, Amplitude algorithm, Logarithmic algorithm and Arctangent algorithm shows that the HM algorithms gave the most accurate values of k.This conclusion agrees with those by both Horton et al. (1983) and Verhoef et al. (1996).
The Conduction-convection algorithm provided realistic daytime soil temperature values.However, it underestimated the soil temperature during the period from 18:00 (LT) to 08:00 (LT), and a noteworthy difference between the measurements and the model output occurred around 00:00 (LT) for all days in this study.A similar underestimation was also encountered by Lin (1980) and Gao et al. (2008).Our explanation is that the model always keeps W constant although the actual W may decrease to zero or even become negative during the nighttime, and also the model does not account for water phase changes that usually happen in shallow soil during the nighttime.
Scatter plots of soil temperature modeled by using the six algorithms against the measured soil temperature at the depth of 0.1 m are given in Fig. 6.The results show that the HM algorithms and Conduction-convection algorithm generated larger correlation coefficients (r), than did the other algorithms.All of the regression lines had slopes of 1.
Statistical analyses are also used to examine the error of model output as follows, Where T m is modeled temperature; Another comparison of the accuracies of the six algorithms is shown in Fig. 7 using the empirical probability distribution functions (PDF) of difference between the modeled and measured soil temperatures at the depth of 0.10 m.The differences between the modeled and measured soil temperatures using the HM algorithm ranged between −1 K and 1 K, and most were near zero.The Conduction-convection algorithm generated the second best results.2261

Conclusions
Six algorithms for calculating soil thermal diffusivity are evaluated with shallow soil measurements collected during LOPEX in 2005.The Phase algorithm and the Amplitude algorithm overestimated the phase and overestimated the amplitude of the soil temperature, respectively.Although the Arctangent algorithm and the Logarithmic algorithm only required four measures of temperature spaced equally in time at two depths, the timing of the four measures of temperature affected the values of soil thermal diffusivity greatly.The HM algorithm gave a reasonable result for most days.However, the assumption of repeating periodicity for soil temperature is invalid on cloudy or rainy days.The algorithms mentioned above are based upon the one-dimensional conduction equation.The Conduction-convection algorithm which is based on the onedimensional conduction-convection equation, provided satisfactory results for daytime temperatures, but it systematically underestimated nighttime soil temperatures.Overall, the Conduction-convection algorithm provided better results than all of the other algorithms except for the HM algorithm.Future efforts should focus on combining the HM and the Conduction-convection algorithms in order to develop an improved method that combines the advantages of each algorithm.The new method should include multiple harmonics to describe the upper boundary temperatures and include conduction and convection heat transfer processes.
the vertical gradient of soil diffusivity ( ∂k ∂z ) and the water flux density (− through 26 July), 2005.The soil temperature measured at 0.05 and 0.10 m depths changed diurnally during DOYs 203-207, as shown in Fig. 1a.The amplitudes of the soil temperature decreased and the phases shifted ahead when the soil depth increased.The soil temperature at 0.05 m changed in response to intermittent cloudiness during DOYs 202 and 203.The maximum soil temperature reached 310.78 K on DOY 204, and the minimum soil temperature was 287.77K at 0.05 m depth on DOY 203. Figure 1a also shows that the soil vertical temperature gradient reached 191.20 K m −1

Figure 3 .
Figure 3. Variation of soil thermal diffusivity k (m 2 s -1 ) with volumetric soil water content θ (%) at a bare soil site over the Loess Plateau from DOYs 201 through 207, 2005.

Figure 4 .
Figure 4. Comparisons of soil temperature modeled by using Phase algorithm, Amplitude algorithm, Arctangent algorithm, Logarithmic algorithm, algorithm by Horton et al.(1983), algorithm by Heusinkveld et al.(2004), and Conduction-convection algorithm, against measurements of soil temperature at 0.10 m depth at a bare soil site over the Loess Plateau from DOYs 201 through 207, 2005.
• E, 35.35 • N at an altitude of 1592 m in Pingliang county of Gansu Province in western China.

Table 1 .
The smallest value of k is 0.6×10 −7 m the Conduction-convection algorithm are larger than the values derived by the other algorithms.The two HM algorithms and the Amplitude algorithm provide similar estimates of the mean values of k.
2 s −1 , and it is obtained from the Arctangent algorithm on DOY 202.The maximum value is 5.47×10 −7 m 2 s −1 , and it is obtained from the Phase algorithm on DOY 207.The maximum, minimum and mean values of k calculated by the Phase al-gorithm and by T is measured temperature; n is the total number of data points; SEE is the standard error of the estimate; and NSEE is a normalized SEE which denotes an estimate of relative uncertainty.The statistical indices SEE and NSEE are presented in Table 4 for the modeled period.It is obvious that the HM algorithm has the lowest values both of SEE and NSEE.The Conduction-convection algorithm has the second lowest values of SEE and NSEE.

Table 2 .
The values of W calculated by Conduction-convection algorithm for the layer of 0.05-0.10mon the Loess Plateau fromDOY 201 to DOY 207, 2005.

Table 3 .
The values of the phase shift and the logarithm of amplitude ratio of soil temperature obtained by using one sine function approximation algorithm at the 0.05 m and 0.10 m depths on the Loess Plateau from DOY 201 toDOY 207, 2005.

Table 4 .
Computed Standard Error of the Estimate (SEE) and Normalized Standard Error of the Estimate (NSEE) of soil temperature at 0.10 m depth on the Loess Plateau from DOY 201 to DOY 207, 2005.