Soil heat transfer is complex, and conduction-alone models may not always perform well in estimating soil apparent thermal diffusivity. Soil apparent thermal diffusivity is related to soil temperature change propagation rates. Soil temperature data collected at the Tazhong station in China were used to examine the characteristics of soil apparent thermal diffusivity determined by three different algorithms and the sum of vertical gradient of soil apparent thermal diffusivity and apparent water flux density . The results showed that 1) soil apparent thermal diffusivity obtained with a conduction–convection algorithm had a better agreement with soil apparent thermal diffusivity obtained with a phase algorithm than with soil apparent thermal diffusivity obtained with an amplitude algorithm except for the case of = 0; 2) when > 0, , and when < 0, ; 3) for a given soil temperature phase shift, increased (decreased) with increasing logarithmic amplitude attenuation when the phase shift was larger (smaller) than the logarithmic amplitude attenuation, reached a maximum value when the phase shift equaled the logarithmic amplitude attenuation, and increased with increasing logarithmic amplitude attenuation; and 4) for a given logarithmic amplitude attenuation, decreased with increasing phase shift and increased (decreased) with increasing phase shift when the phase shift was larger (smaller) than times the logarithmic amplitude attenuation. These mathematical conclusions were also confirmed with field data.
Soil apparent thermal diffusivity was reported to be associated with transient processes of heat conduction and intraporous convection (Zhang and Osterkamp 1995; de Silans et al. 1996), and it had importance for estimating soil temperature in land surface models (Dickinson et al. 1993; Sellers et al. 1996; Chen and Dudhia 2001; Dai et al. 2003). Gao et al. (2008) presented three soil apparent thermal diffusivity algorithms and tested them using data collected in 2005 during the Loess Plateau Land Surface Process Field Experiment (LOPEX). Two of the algorithms, the phase algorithm and the amplitude algorithm, included conduction heat transfer only and assumed that the soil was vertically homogenous. The conduction–convection algorithm included both conduction and convection heat transfer and allowed for vertical heterogeneity of soil apparent thermal diffusivity. Considering soil convection heat transfer, Gao et al. (2010) characterized the phase difference between soil surface heat flux and temperature and investigated its influences on soil surface energy balance closure by using theoretical analysis and experimental evaluations. Wang et al. (2010) compared six soil thermal diffusivity algorithms (amplitude, phase, arctangent, logarithmic, harmonic, and conduction–convection) using data collected at a site in the Loess Plateau of China, and the results showed that the harmonic algorithm gave the most reliable estimates among the six algorithms. Wang (2014) examined the wave phase differences between soil temperature and soil heat fluxes and found that phase lags exist among all land surface energy budgets and the land surface temperature. Differences in phase and time evolution of the land surface temperature and the ground heat flux manifested in diurnal and in annual cycles. An et al. (2016) evaluated four apparent thermal diffusivity algorithms (amplitude, phase, harmonic, and conduction–convection) under different climatic conditions and found that the performances of different algorithms varied with different climatic conditions. Hu et al. (2016) replaced the soil upper boundary condition described with a sinusoidal function in the conduction–convection algorithm of Gao et al. (2008) by that described with a Fourier series, presented an analytical solution to the conduction–convection equation, and derived the mean soil apparent thermal diffusivity and apparent water flux density using a genetic algorithm method. To date, the relative magnitudes of soil apparent thermal diffusivity (obtained with a conduction–convection algorithm , a phase algorithm , and an amplitude algorithm ) derived from the three algorithms have not been evaluated. The comprehensive influences of soil temperature phase shift and logarithmic amplitude attenuation on and (the sum of vertical gradient of soil apparent thermal diffusivity and apparent water flux density), and the relationship between and have not yet been fully quantified. Therefore, to understand and values quantitatively, the objectives of this paper were to 1) compare three soil apparent thermal diffusivity algorithms, that is, , , and ; 2) analyze the influence of phase shift and logarithmic amplitude attenuation on and ; and 3) validate the mathematical findings with field data.
2. Materials and methods
Data were collected at the Taklimakan Desert Atmosphere and Environment Observation Experiment Station located at Tazhong, China (hereinafter referred to as the Tazhong station), from day of year (DOY) 98 to 104 in 2011. This site was located at 38.98°N, 83.64°E, with an altitude of 1103 m. The ground surface was bare and relatively flat. The soil at the site was predominantly fine sand. The site had an arid climate zone with maximum and minimum air temperatures of 319 and 240 K, respectively. The mean annual air temperature and precipitation were 285 K and 24 mm, respectively. The site had 2690 h of sunshine and 263 frost-free days per year (averaged from 1996 to 2010).
Soil temperatures were measured at the surface and at 0.10-, 0.20-, and 0.40-m depths by temperature probes (model 109L, Campbell Scientific Inc.), with an accuracy of ±0.2 K. Soil volumetric water contents (VWCs) were measured at 0.025-, 0.10-, 0.20-, and 0.40-m depths by soil moisture sensors (model CS616, Campbell Scientific Inc.), with an accuracy of ±2.5%. The soil temperature and moisture sensors were sampled each second, and all of the sensor outputs were averaged over 30-min time periods and recorded. Soil sensor information is provided in Table 1. Standard micrometeorological measurements were also made at the site, including four radiation components, wind speed, wind direction, air temperature, relative humidity, air pressure, and precipitation.
Gao et al. (2008) reviewed the classical thermal conduction equation for soil temperature and expanded the equation by adding a term for convection:
where . Variable consists of two parts: 1) , the vertical gradient of soil apparent thermal diffusivity; and 2) , which is associated with the apparent water flux density (positive downward). Variable is the heat capacity of water, is the heat capacity of soil, represents the apparent liquid water velocity (positive downward), and is the volumetric soil water content.
For a boundary condition at depth of , where is the arithmetical average of the daytime maximum soil temperature and the nighttime minimum soil temperature, is half the difference between the daytime maximum value and the nighttime minimum value for soil depth , is the angular velocity of Earth’s rotation, and is the initial phase of soil temperature at depth . The soil temperature at a depth is calculated as (Gao et al. 2008)
where , , and , where is the damping depth of the diurnal temperature wave. It must be noted that Eq. (2) is not a solution to the equation with variable , but it approximates the solution of a soil layer for which varies slightly with depth.
1) The mathematical interrelationships among the three algorithms for calculating apparent thermal diffusivity
If the soil phase shift and the logarithmic amplitude attenuation for two depths are known, the ratio of the soil apparent thermal diffusivities can be determined.
Relationships among the three thermal diffusivity algorithms were derived mathematically. Based on the derived relationships, it can be concluded that when > 0, , and when < 0, . No matter whether > 0 or < 0, is always closer to than to .
2) The mathematical relationships of and to and
Partial derivatives were calculated to obtain the variations of and with phase shift or with logarithmic amplitude attenuation. In summary,
is closer to than to , except for the case of = 0;
if > 0, ;
if < 0, ;
when the phase shift is constant, increases (decreases) with increasing if and reaches a maximum value at and increases with increasing ;
when is constant, decreases with increasing , and increases (decreases) with increasing when .
3) Data processing
Sine functions, , were used to approximate the temporal variation curves of soil temperature collected at the surface and at 0.05-, 0.10-, and 0.20-m depths. Values of and were determined at each depth. The temporal variations of and were determined with Eqs. (3), (4), (6), and (7) using and , where i = 1, 2, and 3.
The interrelationships among , , and obtained in section 2b were evaluated by using soil temperature values and volumetric water contents measured at the Tazhong station from DOY 98 to 104 in 2011. These 7 days were selected for this analysis because of cloud-free conditions. The weather conditions of the 7-day period are given in Table 2.
Figure 1 shows temporal variations in the soil temperature values measured at the surface and at 0.10-, 0.20-, and 0.40-m depths and volumetric water contents measured at 0.025-, 0.10-, 0.20-, and 0.40-m depths. The soil temperature values fluctuated diurnally. Soil temperature amplitudes decreased with depth, and soil temperatures shifted phase as soil depth increased. The maximum (minimum) surface soil temperature during this 7-day period was 330.8 K (268.5 K), recorded on DOY 104. It can be seen in Fig. 1a that the vertical temperature difference between the surface and the 0.10-m depth was 35.4 K at 1300 local time on DOY 102.
Figure 1b shows the temporal variation in the volumetric soil water content. The volumetric water content also changed diurnally, and the diurnal variations at shallow depths were larger than those at deeper depths. The volumetric water content was greatest at the 0.025-m depth, and the largest diurnal variation occurred at the surface. The magnitudes of the volumetric water content were similar at the 0.10- and 0.20-m depths.
To verify the mathematical conclusions for the relative magnitudes of the values from the three algorithms given in section 2b(1), and for two soil layers (0.00–0.10 and 0.10–0.20 m) determined by Eqs. (3), (4), (6), and (7) using the data in Fig. 1 are shown in Fig. 2. Here, and were determined by using the amplitude and phase values of the observed soil temperatures (Gao et al. 2008). Because of small diurnal variation in soil temperature at the 0.40-m depth, and for the 0.20–0.40-m soil layer were not calculated. It can be seen that the value is more similar to the value than it is to the value, when > 0 (i.e., for the shallower soil layer), and when < 0 (i.e., for the deeper soil layer). The phase shift and logarithmic amplitude attenuation were also calculated: and for the shallower soil layer and and for the deeper layer. Substituting the phase shift and logarithmic amplitude attenuation values into Eqs. (A1) and (A2) gave , , and [the ratio of Eqs. (A1) and (A2)] for the shallower soil layer and , , and [the ratio of Eqs. (A1) and (A2)] for the deeper soil layer. These results also supported the mathematical conclusions regarding the relative magnitudes of , , and given in section 2b(1), that is, was always closer to than to , when > 0, and when < 0.
The diurnal and average values of and are presented in Table 3. The , , and values were stable for each soil layer over the 7-day time period. The values of and differed for each soil layer. The values varied from layer to layer, and they varied over time. The average values of , , and were 2.24 × 10−7 (1.96 × 10−7) m2 s−1, 1.79 × 10−7 (2.35 × 10−7) m2 s−1, and 2.22 × 10−7 (1.96 × 10−7) m2 s−1 at the shallower (deeper) layer, respectively, showing that was vertically heterogeneous. The average value of for the shallower soil layer (0.00–0.10 m) was positive (6.41 × 10−7 m s−1), and it was negative (−4.82 × 10−7 m s−1) for the deeper layer (0.10–0.20 m).
To illustrate the mathematical conclusions of Eqs. (3) and (4) given in section 2b(2), Figs. 3 and 4 present the variations of and against phase shift and logarithmic amplitude attenuation . It is obvious that systematically increased with decreasing , reached its maximum value at point A, where was equal to , and reached its minimum value at point B, where was maximum and was minimum. When the phase shift was constant, increased (decreased) with increasing when , and reached its maximum value when ; when was constant, increased with decreasing . Figure 3 comprehensively illustrates the influences of phase shift and logarithmic amplitude attenuation on . Figure 4 is the same as Fig. 3, except that it shows the variation of against phase shift and logarithmic amplitude attenuation , as shown in Eq. (4). Variable reached its maximum value at point A when was minimum and was maximum and reached its minimum value at point B when both variables were at the minimum. These results again supported the mathematical conclusions of Eq. (4) in section 2b(2), that is, when was constant, increased with increasing , and when was constant, increased (decreased) with increasing when .
Both and were functions of and . To illustrate the relationship between and , Fig. 5 shows the variation of with respect to obtained by combining Figs. 3 and 4. Each point expresses the magnitude of and at a certain coordinate pair . It can be seen that 1) when is constant (along the curved arrow), first increases then decreases with increasing and increases with increasing ; and 2) when is constant (along the oblique arrow), decreases with increasing , and increases (decreases) with at smaller (larger) .
Gao et al. (2008) presented two-dimensional figures of and against phase shift and against logarithmic amplitude attenuation. They reported that and varied with phase shift and with logarithmic amplitude attenuation, and they showed the relative magnitudes of obtained with the three algorithms (the conduction–convection, the phase, and the amplitude) when had different signs. The mathematical conclusions based on the relative magnitudes of the values from the three algorithms given in section 2b(1) supported the conclusions reported by Gao et al. (2008). The relative magnitudes of , , and were proven correct theoretically for the condition of with opposite signs in soil layers. This discovery improves our comprehensive understanding of the three soil apparent thermal diffusivity values.
To further understand the variations of and against phase shift and logarithmic amplitude attenuation, theoretical three-dimensional plots of and are shown in Figs. 3 and 4. The three-dimensional plots provide guidance for deriving several mathematical conclusions based on how and vary with phase shifting and logarithmic amplitude attenuation.
Equations (3) and (4) demonstrated that and depended significantly on phase shift and logarithmic amplitude attenuation of soil temperature collected at two depths; thus, accurate determination of the soil temperature phase and amplitude was essential. Figures 3 and 4 also showed that small changes in phase shift and/or logarithmic amplitude attenuation could lead to significant changes in and . Thus, these findings provided an explanation on why temporal variations of field soil temperature yielded different and .
Variable had minimum values when phase shifts were maximum, namely, when the soil phase difference between two depths was largest. Physically speaking, the smaller the soil apparent thermal diffusivity, the larger the phase shift with depth and the greater the decrease in temperature amplitude with depth. Variable reached its maximum values when phase shifts equaled logarithmic amplitude attenuations, providing equal values of and .
The relationships of and are shown in Fig. 5 with different coordinate pairs of phase shift and logarithmic amplitude attenuation. The figure represents a new finding visually by simultaneously capturing the influences of phase shift and logarithmic amplitude attenuation on and .
The and and and values were derived based on the one-dimensional heat conduction and conduction–convection equations when the soil surface temperature boundary condition was assumed to be a periodic sinusoidal curve. In reality, the in situ diurnal change in soil temperature did not strictly follow a sinusoidal curve. If a Fourier series had been used as the surface boundary condition, the apparent thermal diffusivity could be selected to minimize the sum of squared differences between the calculated and measured soil temperature values (Horton et al. 1983), or it could have been estimated from an iteration process by fitting the amplitude and phase of soil temperature at one depth (Heusinkveld et al. 2004). Using a genetic algorithm, Hu et al. (2016) determined and for a Fourier series surface boundary condition. Thus far, no analytical expressions are available to derive relationships between and for a Fourier series surface boundary condition.
This study did not account for water phase changes in soil, although such phase changes (i.e., water vaporization, vapor condensation, ice freezing, and ice thawing) could at times play an important role in the surface energy budget. Future studies should further investigate the role of soil water phase changes on estimation of soil thermal properties.
This study extends earlier findings on soil apparent thermal diffusivity and the sum of vertical gradient of soil apparent thermal diffusivity and apparent water flux density. Important mathematical relationships are derived, including the relative magnitudes of determined by the three algorithms (the conduction–convection, the phase, and the amplitude), and the variations of and with phase shift and/or logarithmic amplitude attenuation in theoretical three-dimensional plots. The derived mathematical relationships are validated with field data collected at the Tazhong station. The relationships derived and validated in this study provide increased understanding of and and of soil thermal property determination.
This study was supported by the Ministry of Science and Technology of China (JFYS2016ZY01002213), National Natural Science Foundation of China (Grant 41275022), the State Scholarship Fund of China Scholarship Council (201608320197), Jiangsu Province graduate education innovation project (KYLX16_0943), and International S&T Cooperation Program of China (ISTCP; Grant 2011DFG23210). We are very grateful to Dr. Aili Maimaitiming (Institute of Desert Meteorology, China Meteorological Administration) for providing data from the Tazhong station. We appreciate the valuable comments of Profs. Zhihua Wang (Arizona State University), Yubin Li (Nanjing University of Information Science and Technology), and Beyong Lee (Nanjing University of Information Science and Technology). We are very grateful to two anonymous reviewers for their careful review and valuable comments, which led to substantial improvement of this manuscript. The code (in MATLAB format) used in this paper can be obtained from the corresponding author and first author.
Details of the Derivation Process
a. The mathematical interrelationships among the three algorithms for calculating apparent thermal diffusivity
If the soil phase shift and the logarithmic amplitude attenuation for two depths are known, the ratio of the soil apparent thermal diffusivities can be determined.
When > 0, the following relationship holds:
Therefore, it can be concluded that when > 0, .
When < 0, the following relationship holds:
Therefore, when < 0, .
We assume that the logarithmic amplitude attenuation and the phase shift . Note that and are both positive real numbers.
Therefore, combining the results in (i) and (ii), we conclude that, except for the case of = 0, is always smaller than . In other words, is always closer to than to .
b. The mathematical relationships of and to and
1) Variation of with and
- When the phase shift is constant, and are both positive real numbers, , and the sign of depends on the relative magnitudes of and , thus,
when , that is, , increases with increasing ;
when , that is, , decreases with increasing ; and
when , that is, , is maximized.
- When the logarithmic amplitude attenuation is constant, the following can be shown from Eq. (A12):
Note that and , so ; thus, decreases with increasing .
2) Variation of with and
- When the logarithmic amplitude attenuation is constant, the following can be obtained from Eq. (A15): , , and the sign of depends on the relative magnitudes of and , thus,
when , that is, , increases with increasing , and
when , that is, , decreases with increasing .