Ippei Iiyama
School of Agriculture, Utsunomiya University, Utsunomiya, Japan.
DOI: 10.4236/gep.2026.144007   PDF    HTML   XML   44 Downloads   291 Views   Citations

Abstract

Soil temperature is a key factor for growth and development of plants, and its regime should be quantified spatially and temporally through numerical simulations with concise boundary conditions. Aiming at evaluating the effectiveness of applying either of the time series of daily-averaged air temperatures and that of daily-averaged soil surface temperature as the surface boundary condition (air-temperature BC and surface-temperature BC), this study numerically simulated seasonal change in soil temperature profile with the time-dependent temperature boundary conditions. The resultant numerical solutions reproduced well the general trends of day-to-day and seasonal changes in measured soil temperatures in both the simulations with the air-temperature BC and that with the surface-temperature BC, while the latter gave the slightly better solution than the former. Based on the concept of thermal time and literature values about thermal time for germination of various herbaceous plants, the numerical solutions also indicated that the sizes of simulation errors can cause the estimation error of less than 1 [d] for predicting the number of dates required for germination under commonly-found situations. Therefore, it was concluded that the time series of air temperatures can become a concise surface boundary condition for predicting or reproducing the seasonal change in soil temperature profile when the temporal resolution of the simulation is set as one day, while the time series of surface temperatures is superior to those of air temperatures as the surface boundary condition.

Keywords

Bulk Density, Heat Balance, Heat Capacity, Newton-Raphson Method, Nonlinear Equations, Particle Density, Relative Humidity, Thermal Conductivity, Three Phases, Water Content, Water Potential

Share and Cite:

1. Introduction

Soil temperature is one of deterministic factors for plants in the geosphere since it affects development and growth of plants in various stages. As the place of germination for plants is in soil, optimum ranges in soil temperature have been clarified (Haj Sghaier et al., 2022; Khaeim et al., 2022), as well as mechanisms with which soil temperature can regulate dormancy and germination of plant seeds (Batlla & Benech-Arnold, 2015; Kuroda & Sawada, 2021). Root vegetables are typical plants whose growth and quality are directly affected by root zone temperature (Sakamoto & Suzuki, 2015). Fruit vegetables also have optimum ranges of root zone temperature for the rates of transpiration/respiration (Shishido & Kumakura, 1994) and of nutrient uptake (Nakano, 2007).

Soil temperature is often more influential than air temperature for growth and development of plants. For instance, optimizing soil temperature is effective in improving plants’ physiological conditions even when they are exposed to mal air temperature conditions (Xu & Huang, 2001; Kawasaki et al., 2013, 2014; Wang et al., 2016). High root temperature also has more adverse effects than high shoot temperature in terms of root growth and nutrient uptake (Huang & Xu, 2000).

In addition, degrees to which soil temperature affects plants’ physiology vary depending on development stages of plants, as Arai-Sanoh et al. (2010) reported that high soil temperature before heading stage is more influential on yield, grain quality, and growth in rice plants than that during ripening stage. At the same time, optimum root zone temperatures differ when air temperature changes (Wang et al., 2022; Yamori et al., 2022). These facts suggest that better knowledge of soil temperature regime allows to achieve better crop yields and, thus, soil thermal environment should be quantified as both spatial distribution and temporal change in soil temperature.

Predicting and reproducing temporal changes in soil temperature profile are realized by numerical analysis. The analysis needs to formulate and solve an initial-boundary problem based on the laws of heat transport and heat balance in soil. Among the arguments for the system of analysis, boundary conditions are utilized to model the ways of thermal interaction between the outside and inside of the soil layer in question. Especially, the boundary condition for the soil surface strongly defines the temporal behavior of soil temperature profile, since warming and cooling of a soil layer result mainly from heat exchange between the soil layer and the atmosphere, while temperature in deep part of the soil layer is relatively stable particularly when there is no specific belowground heat source.

A surface boundary condition for a problem of heat flow in soil is primarily formulated as a heat flux flowing into or out of the soil layer at any time. One of typical ways to derive this flux-type boundary condition is to develop and solve the heat balance relation at the land surface under a given set of weather conditions. To develop the relation requires many known inputs, while to solve the relation requires knowledge of some iterative scheme for solving non-linear equations. However, it is often difficult to obtain all the environmental variables necessary for the development of the surface heat balance relation. In such situations, more simplified manners should be favorable for formulating the surface boundary conditions.

A time series of temperature, instead of a time series of heat flux, may be an option of the simplifications of the surface boundary condition. The straightforward candidate for the temperature boundary condition should be soil surface temperature. However, since temperatures tend to be almost uniform vertically when averaged throughout a cycle of day, daily averages in air temperature and in soil surface temperature in a day may be close to each other, even though surface temperature can markedly differ from air temperature at any time in the day. Some past studies suggested that, while heat flux boundary condition should be used for physically based studies, air temperature boundary condition may also be used for predicting soil temperature (Thunholm, 1990), and it may perform with low error (Naranjo-Mendoza et al., 2018). Therefore, a time series of daily averages in air temperature, as well as that in surface temperature, may be used for the surface boundary condition in predicting year-round change of soil temperature profile with the time resolution of one day.

This study aimed at numerically simulating seasonal change in soil temperature profile by applying time-dependent temperature boundary conditions. The surface boundary condition was formulated by the time series of surface temperature and that of air temperature, and each of the numerical results was validated with measured soil temperature profiles for the comparisons of the model performances.

2. Materials and Methods

2.1. Field Data Sets

The data sets for calculating and validating the time series of soil temperature profile were obtained in a meadow field (350 × 80 m; 36˚29'23''N, 139˚59'14''E) located in the Utsunomiya University Farm in Moka city, Tochigi, Japan. The study period was from May 2018 to March 2019. The soil layers below the field were of volcanic ash soils which were classified into the topsoil layer of 0.25 [m] in thick, the transitional layer with 0.3 [m] in thick below the topsoil layer, and the subsoil layer below the two layers. The details of the meadow field in this study were described in Iiyama (2023).

The data sets were comprised of the above-ground variables and the depth profiles of state variables and physical properties of the soil layers in the harvesting area in the meadow field. The variables included air temperatures T_a [K], land surface temperatures T_sur [K], soil temperatures T [K], and volumetric water contents θ [m³∙m⁻³]. All the data sets were rearranged to be the series of daily averages.

The time series of T_a [K] had been measured at 0.7 [m] and 1.8 [m] in height by using temperature loggers (HOBO U23 Pro v2 External temperature/relative humidity data logger “U23-002A”; Onset Computer Corp; Bourne, MA, USA) enveloped in radiation shields (Solar radiation shield “RS1”; Onset Computer Corp; Bourne, MA, USA). This study denoted the two data sets as T_a070 [K] and T_a180 [K]. The time series of T_sur [K] had been evaluated by using the measured T_a070 [K] and T_a180 [K] with the Monin-Obukhov similarity theory. The measurement details about T_a070 [K] and T_a180 [K] and the method for deriving T_sur [K] were explained in Iiyama (2025).

The time series of T [K] had been measured at 0.05, 0.10, 0.20, 0.40, 0.70, and 1.00 [m] below the land surface by using micro-loggers (HOBO 8K pendant temperature/alarm (waterproof) data logger (UA-001-08); Onset Computer Corp; Bourne, MA, USA), and used for validating the numerically determined soil temperature profiles. The time series of θ [m³∙m⁻³] had been measured every 0.1 [m] along depth by using a capacitance-type soil moisture sensing system (Diviner 2000, Sentek Pty Ltd., Stepney, South Australia). The details of the measurements about T [K] and θ [m³∙m⁻³] were described in Iiyama (2023).

The volumetric water content at the surface θ_sur [m³∙m⁻³] was evaluated by assuming the equilibrium between water vapor potential in the air ψ_va [J∙kg⁻¹] and soil water potential at land surface ψ_sur [J∙kg⁻¹].

${\psi }_{sur}={\psi }_{va}$ (1)

And according to the thermodynamics, the water vapor potential ψ_va [J∙kg⁻¹] can be related to the relative humidity of the air h_r [−] as below:

${\psi }_{va}=\frac{R{T}_a}{M_w}\ln h_r$ (2)

where R [J∙mol⁻¹∙K⁻¹] is the ideal gas constant (=8.314 [J∙mol⁻¹∙K⁻¹]), M_w [kg∙mol⁻¹] is the molar weight of water (=0.018 [kg∙mol⁻¹]). The daily averages of h_r [−] were calculated by assuming that a daily average water vapor pressure in air e_a [Pa] is nearly equal to a saturated water vapor pressure e_as [Pa] for the daily lowest temperature T_{a min} [K]:

$h_r=\frac{e_a}{e_{as}}$ (3)

$e_a=e_s\left(T_{amin}\right)$ (4)

$e_{as}=e_s\left(T_a\right)$ (5)

$e_s\left(T_a\right)={\epsilon }_0\exp \frac{{\epsilon }_1\left(T_a-273.15\right)}{{\epsilon }_2+T_a-273.15}$ (6)

where ε₀ = 611 [Pa], ε₁ = 17.502 [K⁻¹], and ε₂ = 240.97 [K] (Buck, 1981).

Then, the time series of ψ_sur [J∙kg⁻¹] calculated from Equations (2)-(6) was converted to the series of θ_sur [m³∙m⁻³] by using the soil water retention curve for the topsoil layer of the study site. The soil water retention curve was expressed as below:

$\theta ={\theta }_c+{\theta }_f$ (7)

${\theta }_c={\theta }_{cmin}+\frac{{\theta }_{cmax}-{\theta }_{cmin}}{1+{\left|{\gamma }_c\psi \right|}^{{\mu }_c}}$ (8)

${\theta }_f={\theta }_{fmin}+\frac{{\theta }_{fmax}-{\theta }_{fmin}}{1+{\left|{\gamma }_f\psi \right|}^{{\mu }_f}}$ (9)

where θ_c [m³∙m⁻³] and θ_f [m³∙m⁻³] are of coarse pore and fine pore systems in the soil, respectively. The values of the parameters were γ_c = 0.2843 [kg∙J⁻¹], γ_f = 0.001400 [kg∙J⁻¹], μ_c = 1.0000 [−], μ_f = 1.0486 [−], θ_{c max} = 0.1686 [m³∙m⁻³], θ_{f max} = 0.5874 [m³∙m⁻³], θ_{c min} = 0.0000 [m³∙m⁻³], and θ_{f min} = 0.0307 [m³∙m⁻³], obtained by fitting this model to the laboratory-measured data sets presented in Iiyama & Hirai (2014). Then, the data set of θ_sur [m³∙m⁻³] was used with the data set of θ [m³∙m⁻³] to form the time series of the entire depth profile of water content.

The soil physical properties including soil bulk density ρ_d [kg∙m⁻³] and soil particle density ρ_s [kg∙m⁻³] were measured. The measurement details were described in Iiyama (2023). The saturated volumetric water content θ_s [m³∙m⁻³] was evaluated as 1 - ρ_d/ρ_s. These soil physical properties were tabulated on Table 1.

Table 1. The soil bulk densities, particle densities, and saturated volumetric water contents of the soil layer in this study. The measurement details are found in Iiyama (2023). Each value is an average of triplicate measurements with a standard deviation in the parentheses. The soil textures of all the soil layers were classified as clay loam soils based on the soil-texture classification defined by the International Union of Soil Science.

Depth range [m]	0 - 0.25	0.25 - 0.55	0.55 - 1.00
Bulk density ρd [kg∙m−3]	601 (62)	600 (23)	602 (23)
Particle density ρs [kg∙m−3]	2463 (287)	2808 (38)	2784 (23)
Saturated water content θs [m3∙m−3]	0.756 (0.156)	0.786 (0.041)	0.784 (0.039)

2.2. Heat Balance in Soil

The heat balance in the soil layer was formulated for numerically evaluating the seasonal change in soil temperature profile. Firstly, a soil layer to be analyzed was separated by a series of calculation nodes arranged along depth. The calculation nodes were numbered from 0 to N, and named as z = z_i (0 ≤ i ≤ N). z₀ [m] was set at the bottom of the soil layer while z_N [m] to the land surface. Soil temperatures at z = z_i (0 ≤ i ≤ N) were denoted as T_i [K].

A small segment was defined for every calculation node so that T_i [K] represents the value of soil temperature for a segment around z = z_i [m]. Any of the segment boundaries was placed in the middle of two adjacent calculation nodes so that a segment around z = z_i covers the domain (z_i + z_i−1)/2 ≤ z [m] ≤ (z_i+1 + z_i)/2 for 1 ≤ i ≤ N − 1. The segments for z = z₀ [m] and z = z_N [m] cover the domains z₀ ≤ z [m] ≤ (z₁ + z₀)/2 and (z_N + z_N−1)/2 ≤ z [m] ≤ z_N, respectively. In this study, 20 was assigned to N while z₀ [m] and z_N [m] were set at z = −1 [m] and z = 0 [m]. Other calculation nodes were located every 0.05 [m] interval between z₀ [m] and z_N [m].

The heat flow in the soil layer was formulated between any of two adjacent soil segments. When heat flow in soil owes almost solely on heat conduction driven by temperature gradient along depth, a value of heat flux f_{h i} [J∙m⁻²∙s⁻¹] can be evaluated on a segment boundary between z = z_i and z = z_i+1 as follows:

$f_{h i}=-{\lambda }_i\frac{T_{i+1}-T_i}{z_{i+1}-z_i}$ (10)

where λ_i [W∙m⁻¹∙K⁻¹] is the soil thermal conductivity assigned to the domain z_i ≤ z [m] ≤ z_i+1 (1 ≤ i ≤ N − 1). The formulation of λ_i [W∙m⁻¹∙K⁻¹] was described in the next section.

By using heat flux terms defined above, the heat balance equations for the segments i (1 ≤ i ≤ N − 1) can be expressed as below:

$\frac{C_{hi}^{j+1}{T}_i^{j+1}-C_{hi}^j{T}_i^j}{t_{j+1}-t_j}\frac{z_{i+1}-z_{i-1}}{2}=\frac{f_{hi-1}^{j+1}+f_{hi-1}^j}{2}-\frac{f_{hi}^{j+1}+f_{hi}^j}{2}\left(1\le i\le N-1\right)$ (11)

where $C_{hi}^j$ [J∙m⁻³∙K⁻¹] is the heat capacity of soil in the segment i at a time t_j [s]. The superscript j denotes a time-step number (j ≥ 0) and can be used for indicating that variables with the superscripts j and j + 1 have values occurring at times t_j [s] and t_j+1 [s], respectively. The time step Δt = t_j+1 − t_j [s] of calculation was set as 86400 [s]. This temporal resolution can be practical for prediction purposes because even a modern 5-day weather forecast is as accurate as a 1-day forecast was in 1980 (Alley et al., 2019), though the skills in medium-range weather forecasting have been rapidly improved. The formulation of $C_{hi}^j$ [J∙m⁻³∙K⁻¹] was explained in the section after the next section. The expressions in the right side of the equations for the segments i (1 ≤ i ≤ N − 1) were defined as their time-averaged values between t = t_j [s] and t = t_j+1 [s] and, thus, the Crank-Nicholson scheme was applied.

The two boundary conditions were expressed as below:

$T_i^{j+1}=T_{bottom}^{j+1}\left(i=0\right)$ (12)

$T_i^{j+1}=T_{top}^{j+1}\left(i=N\right)$ (13)

In these equations, T_bottom [K] is the soil temperature at the bottom of the soil layer of interest, and the time series of the daily averages of soil temperature measured at z = −1 [m] was used to define this value-type boundary condition. Neither of a constant value boundary condition nor the zero flux boundary condition can be applicable to setting up this bottom boundary condition, because, according to the possible sizes of heat capacity and thermal conductivity found in the soil layer in this study, the damping depth for the yearly soil temperature fluctuation was likely to be deeper than 1 [m]. T_top [K] is the temperature at the top boundary of the region of analysis, and either of the time series of the daily averages of surface temperature T_sur [K], air temperature measured at 0.7 [m] in height T_a070 [K], and air temperature measured at 1.8 [m] in height T_a180 [K].

The initial condition was set by interpolating the soil temperature profile measured at the first date of the study period.

As described in the later sections, both the soil thermal conductivity and the soil heat capacity were formulated as functions of soil temperature. Therefore, the equations were of non-linear and, the Newton-Raphson iterative method was applied to solve them every time step. The condition of convergence in the Newton-Raphson iteration was defined in such a way that the sum of squared residuals in Equation (11) becomes lower than 10⁻¹², and the iteration numbers found in all of the time steps in every numerical analysis did not go beyond 3.

2.3. Soil Thermal Conductivity

The soil thermal conductivity λ [W∙m⁻¹∙K⁻¹] was calculated by using the following equation:

$\lambda =\frac{{\lambda }_{max}}{1+\left(\frac{{\lambda }_{max}}{{\lambda }_{min}}-1\right)\exp \left(-\kappa S_r\right)}$ (14)

where λ_max [W∙m⁻¹∙K⁻¹] is the ideal maximum value of λ [W∙m⁻¹∙K⁻¹] for the moisture saturated condition, λ_min [W∙m⁻¹∙K⁻¹] is the ideal minimum value of λ [W∙m⁻¹∙K⁻¹] for the driest condition, S_r [−] is the degree of saturation of the soil in question, defined as the ratio of the volumetric water content θ [m³∙m⁻³] to the saturated volumetric water content θ_s [m³∙m⁻³], κ [−] is a parameter that regulates the slope of the S_r − λ relation. This sigmoidal model of λ [W∙m⁻¹∙K⁻¹] can be derived from the following ordinary differential equation:

$\frac{\text{d}\lambda }{\text{d}{S}_r}=\kappa \lambda \left({\lambda }_{max}-\lambda \right)$ (15)

The form of this ODE is based on the following empirical characteristics about the S_r − λ relation: 1) λ [W∙m⁻¹∙K⁻¹] increases monotonically with the increase in S_r [−], 2) λ [W∙m⁻¹∙K⁻¹] has some upper limit λ_max [W∙m⁻¹∙K⁻¹] to which λ [W∙m⁻¹∙K⁻¹] asymptotically gets close when S_r [−] gets close to unity, 3) dλ/dS_r becomes very small when S_r converges on 0.

The ideal maximum value of λ [W∙m⁻¹∙K⁻¹], λ_max [W∙m⁻¹∙K⁻¹], was determined by using the geometric mean model (Woodside & Messmer, 1961; Sass et al., 1971; Johansen, 1975) as follows:

${\lambda }_{max}={\lambda }_w^{1-\sigma }{\lambda }_s^{\sigma }$ (16)

where λ_w [W∙m⁻¹∙K⁻¹] and λ_s [W∙m⁻¹∙K⁻¹] are the thermal conductivities of soil water and soil solids, respectively, σ [m³∙m⁻³] is the volume fraction of solid phase in the soil.

The thermal conductivity of water λ_w [W∙m⁻¹∙K⁻¹] was expressed as the polynomial function of soil temperature T [K] as below:

${\lambda }_w=l_{w0}+l_{w1}T+l_{w2}{T}^2$ (17)

where l_w0 = −9.003748 × 10⁻¹, l_w1 = 8.387698 × 10⁻³, and l_w2 = −1.118205 × 10⁻⁵. The coefficient of determination R² was more than 0.999998. The reference of this function was Ramires et al. (1995).

The thermal conductivity of soil solid materials λ_s [W∙m⁻¹∙K⁻¹] was modeled as below:

${\lambda }_s={\lambda }_o^{1-q}{\lambda }_q^q$ (18)

where λ_q [W∙m⁻¹∙K⁻¹] and λ_o [W∙m⁻¹∙K⁻¹] are the thermal conductivities of quartz and other soil minerals, respectively, q [kg∙kg⁻¹] is the mass fraction of quartz in the soil solid phase. In this equation, 7.7 [W∙m⁻¹∙K⁻¹] was assigned to λ_q [W∙m⁻¹∙K⁻¹] (Johansen, 1975), while 2.1 was given to λ_o [W∙m⁻¹∙K⁻¹] as a typical value for primary parent materials for soils (Haynes & Lide, 2010). The quartz contents q [kg∙kg⁻¹] was assumed to be equivalent to sand content of a soil and 0.36 [kg∙kg⁻¹] (Iiyama, 2024) was assigned to q [kg∙kg⁻¹].

The ideal minimum value of λ [W∙m⁻¹∙K⁻¹], λ_min [W∙m⁻¹∙K⁻¹], was determined by the linear regression analysis for literature values (Lu et al., 2007, 2014; Tarnawski et al., 2013, 2015), expressed as follows:

${\lambda }_{min}=l_{min0}+l_{min1}\left(1-\sigma \right)$ (19)

where l_{min 0} = 0.5126 and l_{min 1} = −0.5894. The coefficient of determination R² was more than 0.6950. This relation had been obtained from a variety of soil classes and covered the wide range of total porosity n = θ_s = 1 − σ [m³∙m⁻³], including both the soil class and the range of total porosity found in this study.

The parameter κ [−] can be determined by using the condition such as:

${\lambda }_{max}=\underset{S_r\to 1}{\lim }\lambda$ (20)

Considering Equation (14), this condition can be transformed as below:

$\delta =\underset{S_r\to 1}{\lim }\left(\frac{{\lambda }_{max}}{{\lambda }_{min}}-1\right)\exp \left(-\kappa S_r\right)\ll 1$ (21)

And, solving this equation for κ [−] gives:

$\kappa =\ln \frac{\frac{{\lambda }_{max}}{{\lambda }_{min}}-1}{\delta }$ (22)

In this study, 0.005 was assigned to δ [−] as a small value.

As a result, λ [W∙m⁻¹∙K⁻¹] was modeled as the function of q [kg∙kg⁻¹], T [K], σ [m³∙m⁻³], and θ [m³∙m⁻³]. In the numerical scheme, λ_i [W∙m⁻¹∙K⁻¹] was defined in each of the domain z_i ≤ z [m] ≤ z_i+1 so that each of the arguments for the mathematical function λ [W∙m⁻¹∙K⁻¹] was the average of the values on the two adjacent calculation nodes at z = z_i [m] and z = z_i+1 [m].

2.4. Heat Capacity of Soil

The heat capacity of soil C_h [kg∙m⁻³] was expressed as the weighted mean of the heat capacities for solid-, liquid-, and air-phases. The weighting factors for the expression were the volume fractions of the soil three phases such as:

$C_h=\sigma C_{hs}+\theta C_{hw}+a{C}_{ha}$ (23)

where a [m³∙m⁻³] is the air-filled porosity of the soil, identified as a = 1 − θ − σ [m³∙m⁻³], C_hs [J∙m⁻³∙K⁻¹], C_hw [J∙m⁻³∙K⁻¹], and C_ha [J∙m⁻³∙K⁻¹] are the volumetric heat capacity of solid-, liquid-, and air-phases of the soil, respectively. The value of C_hs [J∙m⁻³∙K⁻¹] was approximated as 2.4 × 10⁶ [J∙m⁻³∙K⁻¹], since typical parent materials of soils such as basalt and granite have specific heat ranging from 0.8 to 1.0 [J∙g⁻¹∙K⁻¹] (National Astronomical Observatory of Japan, 2016) with particle density around 2.65 [Mg∙m⁻³]. C_hw [J m⁻³∙K⁻¹] was expressed by such a polynomial function as:

$C_{hw}=c_{hw0}+c_{hw1}{T}_c+c_{hw2}{T}_c^2+c_{hw3}{T}_c^3+c_{hw4}{T}_c^4+c_{hw5}{T}_c^5$ (24)

where T_c [˚C] is the Celsius temperature of the system of interest, equivalent to T – 273.15 [K], c_hw0 = 4.216947 × 10⁶ [J∙m⁻³∙K⁻¹], c_hw1 = −3.252814 × 10³ [J∙m⁻³∙K⁻²], c_hw2 = 8.710579 × 10¹ [J m⁻³ K⁻³], c_hw3 = −1.819027 [J∙m⁻³∙K⁻⁴], c_hw4 = 1.666381 × 10⁻² [J m⁻³∙K⁻⁵], and c_hw5 = −5.707307 × 10⁻⁵ [J∙m⁻³∙K⁻⁶]. The coefficient of determination R² was more than 0.999950. The reference of this function was National Astronomical Observatory of Japan (2016). C_ha [J∙m⁻³∙K⁻¹] is defined as the specific heat of air c_p [J∙kg⁻¹∙K⁻¹] multiplied by the density of air ρ_a [kg∙m⁻³]:

$C_{ha}=c_p{\rho }_a$ (25)

c_p [J∙kg^-1∙K^-1] was expressed as a function of temperature T [K] and air pressure p_a [Pa] derived from the multivariate regression for the data set tabulated in Lemmon (2010):

$c_p=c_{p0}+c_{p1}T+c_{p2}{p}_a$ (26)

where c_p0, c_p1, and c_p2 are 1120.42 [J∙kg⁻¹∙K⁻¹], −0.36905 [J∙kg⁻¹∙K⁻²], and 1.55494 × 10⁻⁵ [J∙kg⁻¹∙K⁻¹∙Pa⁻¹], respectively. The coefficient of determination R² of this regression was 0.90953. ρ_a [kg∙m⁻³] was calculated by assuming that the soil air behaves as an ideal gas:

${\rho }_a=\frac{p_a{M}_a}{RT}$ (27)

where M_a [kg∙mol⁻¹] is the molar weight of soil air, assumed to be 0.028966 [kg∙mol⁻¹], R [J∙mol⁻¹∙K⁻¹] is the ideal gas constant (=8.314 [J∙mol⁻¹∙K⁻¹]). The air pressure p_a [Pa] in the soil layer was assumed to be equal to atmospheric pressure.

As described above, C_h [J∙m⁻³∙K⁻¹] was formulated as the function of p_a [Pa], T [K], σ [m³∙m⁻³], and θ [m³∙m⁻³].

3. Results and Discussion

Figure 1 shows the time series of T_a180 [K] and T [K], T_a180 - T_a070 [K], T_a180 - T_sur [K], and θ [m³∙m⁻³]. T_a180 - T_a070 [K] and T_a180 - T_sur [K] were plotted, instead of each of T_a070 [K], T_a180 [K], and T_sur [K], because it was difficult to distinguish the curves of T_a070 [K], T_a180 [K], and T_sur [K] when any two of the three were plotted on the same graph area.

The time series of air and soil temperatures (Figure 1(i)) indicated that the diurnal changes in air temperature had affected mainly those in the soil temperatures in the topsoil layer while the soil temperatures below 0.4 [m] in depth had not had clear daily fluctuations but showed only seasonal changes.

The time series of the difference among the daily-averaged air temperatures (Figure 1(ii)) showed that the difference between the two heights was negligible or very small at any time in the year. When they had evaluated every 30 minutes, the two temperatures significantly differed from each other (Iiyama, 2025). Therefore, it was suggested that the difference in air temperature between the two heights can be cancelled when it is summed up for each day.

Figure 1. The time-series of (i) the air temperature measured at 1.8 [m] in height T_a180 [K] and the soil temperatures T [K] measured at −0.05, −0.10, −0.20, −0.40, −0.70, and −1.00 [m] in height from the land surface, (ii) the difference between T_a180 [K] and the air temperature measured at 0.7 [m] in height T_a070 [K], (iii) the difference between T_a180 [K] and the surface temperature T_sur [K] evaluated in Iiyama (2025), and (iv) the volumetric water content θ [m³∙m⁻³] measured at every 0.1 [m] along depth. Since the difference in θ [m³∙m⁻³] between any two adjacent measurement depths had been small below 0.2 [m] in depth, only the data sets from 6 measurement heights of −0.1, −0.2, −0.4, −0.5, −0.9, −1.0 [m] from the land surface were plotted.

The time series of the difference between the surface temperature and the air temperature at z = 1.8 [m] had fluctuated from day to day (Figure 1(iii)). The sizes in the difference were within 2 [K] in many of the dates and rarely went beyond 3 [K]. The signs of T_a180 - T_sur [K] tended to be negative in the spring and summer, and to be positive in the winter.

The time series of the volumetric water contents θ [m³∙m⁻³] (Figure 1(iv)) showed that the topsoil layer had experienced severely-dried conditions several times in the summer. On the other hand, the moisture condition below 0.5 [m] in depth was stable with more than 0.55 [m³∙m⁻³] of water content throughout the year. Even the values of θ [m³∙m⁻³] at z = −0.2 [m] and −0.4 [m] had fluctuated more moderately than at z = −0.1 [m]. These moisture conditions suggested that the soil thermal conductivity and the soil heat capacity were temporally stable in most of the soil layer in the study site.

Figure 2 shows the comparisons between the time series of the measured and simulated soil temperatures. As the results in Figure 1(ii) meant, both of the time series of the two air temperatures T_a070 [K] and T_a180 [K] gave almost the same numerical solutions of the soil temperature profiles when they had been used as the top boundary condition in Equation (13). Thus, only the simulation results for the case in which T_a180 [K] was adopted were displayed in Figure 2(a1) through Figure 2(a5). And the results obtained with T_sur [K] as the top boundary condition were plotted in Figure 2(s1) through Figure 2(s5).

The general trends of day-to-day and seasonal fluctuations in the measured soil temperatures were well followed by both of the two simulations in all of the five depths of measurement.

The simulated curves tended to underestimate the measured curves from May to June in 2018 in many of the depths of measurement. In addition, in z = −0.05 [m], the portions of the simulated curves from January to February also overestimated the measured curves. And the sizes of these overestimations and underestimations were larger in the simulation results for which T_a180 [K] was used as the top boundary condition.

The features of these overestimations and underestimations implied that although the time series of T_a180 [K] and T_sur [K] shared the same characteristics in day-to-day and seasonal fluctuations, the amplitude of the yearly temperature cycle was very slightly narrower in T_a180 [K] than in T_sur [K], and thus, the simulation results for which T_sur [K] was adopted reproduced better the actual soil temperature regimes.

For quantifying the simulation performances between the three simulation cases in which the surface boundary conditions were given by (i) T_a180 [K], (ii) T_a070 [K], and (iii) T_sur [K], the absolute errors of the simulated soil temperatures compared to the measured soil temperatures were tabulated (Table 2). A value in each component in Table 2 is calculated by the following expression:

$\text{d}{T}_{ave}=\frac{1}{k_{max}}{\displaystyle \sum }_{k=1}^{k_{max}}\left|T_{meask}-T_{simulk}\right|$ (28)

where T_meas [K] and T_simul [K] are the measured and simulated soil temperatures at a vertical location for a date, and k_max is the number of data samples equal to 333 as the total dates in the study period from 2018/5/3 to 2019/3/31.

The sizes of estimation error dT_ave [K] took the maximum in the case (i) among the three kinds of simulations at every vertical location while the minima were found in the case (iii). Almost the same trend was found for the values of dT_sd [K], the standard deviation for dT_ave [K]. Thus, the surface temperature T_sur [K] gave the best result as the surface boundary condition of the simulation, and the simulation performance certainly lowered as the location of air temperature measurement was set further away from the land surface.

Figure 2. The time-series of the measured and simulated soil temperatures. The dotted and solid lines in each sub-graph indicate the measured and simulated curves, respectively. The simulated curves in the sub-graphs (a1) through (a5) were obtained when the top boundary condition was defined by the time series of T_a180 [K], while those in the sub-graphs (s1) through (s5) were given when the time series of T_sur [K] was used to formulate the top boundary condition. The vertical locations of evaluation were denoted by z = −0.05, −0.10, −0.20, −0.40, and −0.70 [m], as labelled in each of the sub-graphs.

Table 2. The absolute errors of the simulated soil temperatures compared to the measured soil temperatures. The rows are assigned to the three cases of simulation in which the surface boundary conditions were given by the time series of (i) air temperatures measured at 1.8 [m] in height, (ii) air temperatures measured at 0.7 [m] in height, and (iii) surface temperatures evaluated in Iiyama (2025). The columns are assigned to the vertical locations of measurements z = −0.7, −0.4, −0.2, −0.1, and −0.05 [m]. A value in each component is defined as dT_ave [K], which indicates the average of |T_meas - T_simul| [K] for totally 333 dates in the study period from 2018/5/3 to 2019/3/31, where T_meas [K] and T_simul [K] are the measured and simulated soil temperatures at a vertical location for any date. The values in the parentheses indicate the standard deviations dT_sd [K] for dT_ave [K].

vertical locations z [m]	−0.7	−0.4	−0.2	−0.1	−0.05
(i) surface BC with Ta180	0.399 (0.220)	0.807 (0.455)	0.985 (0.615)	0.874 (0.659)	1.002 (0.780)
(ii) surface BC with Ta070	0.386 (0.225)	0.779 (0.455)	0.978 (0.601)	0.852 (0.647)	0.973 (0.759)
(iii) surface BC with Tsur	0.330 (0.204)	0.657 (0.374)	0.825 (0.517)	0.711 (0.566)	0.844 (0.676)

The sizes of dT_ave [K] became smaller at z = −0.7 [m] and −0.4 [m] than those in the topsoil layer in any of the three simulation cases, while the maxima of dT_ave [K] were found for z = −0.05 [m]. At the same time, dT_sd [K] were obtained in the ascending size order from the lower to higher vertical locations. Thus, although the orders of the sizes of dT_ave [K] did not completely agree with the order of the vertical location z [m], it was suggested that soil temperatures can be reproduced or predicted more correctly in deeper layers than in shallower layers, presumably because of relatively stable soil temperature regime in deeper part of the entire soil layer.

Although the superiority of using T_sur [K] as the top boundary condition was decisive, this result can also support that air temperatures can become a concise top boundary condition for predicting or reproducing a time series of soil temperature profile. In quantity, the sizes of dT_ave [K] were about 1 [K] at maximum, and even the sizes of the “three-sigma” 3 dT_sd [K] were less than 2.34 [K]. These results implied that even a time series of soil temperatures simulated with air temperatures enables to draw a curve that seems so well imitating a measured time series that the two curves on a chart are similar to each other when looked briefly, as shown in Figure 2(a).

To evaluate whether these sizes of estimation error dT_ave [K] may be practically small, it was considered how an estimation error for soil temperature can affect the prediction of some phenology of plants. A straightforward example is seed germination that is often predicted by applying the concept of thermal time for germination with temperature regime in a topsoil layer. According to the concept of thermal time, when a base temperature and a thermal time for germination of a plant species is expressed as T_b [K] and S [K d], the days for germination n_g [d] under a temporal mean of soil temperature T_m [K] can be estimated as below:

$n_g=\frac{S}{T_m-T_b}$ (29)

Then, if the soil temperature is erroneously evaluated as T_me [K] with an error dT [K] such as:

$T_{me}=T_m+\text{d}T$ (30)

The days required for germination should also be wrongly estimated as n_ge [d] in such a way as:

$n_{ge}=\frac{S}{T_{me}-T_b}$ (31)

Therefore, subtracting Equation (29) from Equation (31) gives the error in the estimation of the time period for germination dn_g = |n_ge − n_g| [d] as below:

$\text{d}{n}_g=\left|\frac{-\text{d}T\cdot S}{\left(T_m+\text{d}T-T_b\right)\left(T_m-T_b\right)}\right|$ (32)

Pairs of T_b [K] and S [K d] had been studied for many herbaceous plants with wide ranges in soil temperature (Trudgill et al., 2000; Lonati et al., 2009; Wu et al., 2023; Sampayo-Maldonado et al., 2025). By applying Equation (32) to the data sets of T_b [K] and S [K∙d] in these past studies, it was evaluated that dT of 1 [K] can cause dn_g [d] of less than 1 [d] with dn_g/n_g [−] of less than 0.1 in most cases, particularly when T_m [K] is sufficiently higher than T_b [K]. In this sense, it was suggested that dT ≤ 1 [K] can be practically small, and that the time series of air temperatures can be concisely adopted for the formulation of the surface boundary condition of a soil temperature simulation, though surface temperatures can give better simulation results than air temperatures.

4. Conclusion

Based on the hypothesis that a time series of daily averages in air temperature, as well as that in surface temperature, may be available as the surface boundary condition for simulating seasonal change in soil temperature profile, this study numerically reproduced seasonal change in soil temperature profile by using time-dependent temperature boundary conditions. The surface boundary condition was formulated by the time series of daily averaged air temperatures measured at 0.7 and 1.8 [m] in height (T_a070 [K] and T_a180 [K]) and by that of surface temperature which had evaluated with the measured air temperatures, wind speed, and turbulent theory (T_sur [K]). The simulation process used in this study was based on the heat balance concept and the law of heat conduction, and all the parameters required for the simulation were soil physical properties of easily obtainable, though some difficulties remain on parameterizing the soil quartz content q [kg∙kg⁻¹], the small value δ [−] for regulating the slope of the S_r − λ relation, and the water potential at the surface ψ_sur [J∙kg⁻¹]. Measured soil temperatures were used for validating each of the numerical results for the comparisons among the simulation performances.

Since T_a070 [K] and T_a180 [K] were very similar to each other, the simulated soil temperature profiles were almost equivalent among the cases in which either of them was used as the top boundary condition, suggesting that the selection of the vertical location of air temperature measurement does not affect so strongly the simulation results with one-day resolution as selecting either of the air temperatures or the surface temperature.

The general trends of day-to-day and seasonal fluctuations in the measured soil temperatures were reproduced well by both of the simulations with T_a180 [K] and T_sur [K]. And the case in which T_sur [K] was used as the top boundary condition gave the slightly better solution than that with T_a180 [K], mainly because the amplitude of the yearly temperature cycle was slightly narrower in the air temperature than in the surface temperature.

For quantifying the simulation performances between the simulation cases, the absolute errors of the simulated soil temperatures from the measured soil temperatures were evaluated. According to the temporal averages of the absolute error dT_ave [K], T_sur [K] gave the best result as the surface boundary condition of the simulation, and the simulation performance lowered as the location of air temperature measurement was set further away from the land surface. However, since the sizes of dT_ave [K] were about 1 [K] at maximum, and even the sizes of the “three-sigma” 3 dT_sd [K] were less than 2.34 [K] for all of the three simulation cases, it was implied that even using a time series of air temperatures as the top boundary condition makes the simulation results well imitate a measured time series of soil temperatures.

The sizes of dT_ave [K] obtained in this study were seemingly small from a practical point of view, as exemplified in the evaluation of the effect of dT_ave [K] on the estimation of the time period required for plant germination. Based on the concept of thermal time and literature values about thermal time for germination, it was estimated that dT_ave [K] in this study can cause the estimation error of less than 1 [d] for predicting the date of germination under commonly-found situations.

In conclusion, while the time series of surface temperatures is superior to those of air temperatures as the top boundary condition, the time series of air temperatures can become a concise top boundary condition for predicting or reproducing the seasonal change in soil temperature profile when the temporal resolution of the simulation is set as one day. In this study, however, the effects of neglecting the surface heat balance relation on the simulation results have yet to be evaluated. Therefore, further studies should include the evaluation of possible difference in simulated soil temperature profiles between the case in which the surface heat balance relation is incorporated into the process of analysis and the case in which one of the simplified surface boundary conditions is used as in this study.

Acknowledgements

The author thanks Mr. T. Shiozawa, Mr. E. Saito and Mr. N. Yamaguchi in the Utsunomiya University Farm for their supports in the managements of the study field.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Alley, R. B., Emanuel, K. A., & Zhang, F. (2019). Advances in Weather Prediction: Better Weather and Environmental Forecasting Will Continue to Improve Well-Being. Science, 363, 342-344. [Google Scholar] [CrossRef] [PubMed]
[2]	Arai-Sanoh, Y., Ishimaru, T., Ohsumi, A., & Kondo, M. (2010). Effects of Soil Temperature on Growth and Root Function in Rice. Plant Production Science, 13, 235-242. [Google Scholar] [CrossRef]
[3]	Batlla, D., & Benech-Arnold, R. L. (2015). A Framework for the Interpretation of Temperature Effects on Dormancy and Germination in Seed Populations Showing Dormancy. Seed Science Research, 25, 147-158. [Google Scholar] [CrossRef]
[4]	Buck, A. L. (1981). New Equations for Computing Vapor Pressure and Enhancement Factor. Journal of Applied Meteorology, 20, 1527-1532. [Google Scholar] [CrossRef]
[5]	Haj Sghaier, A., Tarnawa, Á., Khaeim, H., Kovács, G. P., Gyuricza, C., & Kende, Z. (2022). The Effects of Temperature and Water on the Seed Germination and Seedling Development of Rapeseed (Brassica napus L.). Plants, 11, 2819. [Google Scholar] [CrossRef] [PubMed]
[6]	Haynes, W., M., & Lide, D., R. (2010). CRC Handbook of Chemistry and Physics (91st ed.). CRC Press.
[7]	Huang, B., & Xu, Q. (2000). Root Growth and Nutrient Element Status of Creeping Bentgrass Cultivars Differing in Heat Tolerance as Influenced by Supraoptimal Shoot and Root Temperatures. Journal of Plant Nutrition, 23, 979-990. [Google Scholar] [CrossRef]
[8]	Iiyama, I. (2023). Seasonal Change in CO2 Production Rate along Depth in a Grassland Field. Journal of Geoscience and Environment Protection, 11, 106-124. [Google Scholar] [CrossRef]
[9]	Iiyama, I. (2024). Characteristics of In-Situ Soil Water Hysteresis Observed through Multiple-Years Monitoring. Journal of Geoscience and Environment Protection, 12, 162-175. [Google Scholar] [CrossRef]
[10]	Iiyama, I. (2025). Diurnal and Seasonal Trends in Aerodynamic Conductance over a Grassland Surface. Journal of Geoscience and Environment Protection, 13, 42-58. [Google Scholar] [CrossRef]
[11]	Iiyama, I., & Hirai, T. (2014). Subsoil Water Available in Suction and Poor in Amount under Continuous Grassland Use. Soil Science and Plant Nutrition, 60, 439-447. [Google Scholar] [CrossRef]
[12]	Johansen, O. (1975). Varmeledningsevne av jordarter (Thermal Conductivity of Soils). Trondheim, Group for Thermal Analysis of Frost in the Ground. Institute for Kjoletek-nikk, 231p. (Translated in English)
[13]	Kawasaki, Y., Matsuo, S., Kanayama, Y., & Kanahama, K. (2014). Effect of Root-Zone Heating on Root Growth and Activity, Nutrient Uptake, and Fruit Yield of Tomato at Low Air Temperatures. Journal of the Japanese Society for Horticultural Science, 83, 295-301. [Google Scholar] [CrossRef]
[14]	Kawasaki, Y., Matsuo, S., Suzuki, K., Kanayama, Y., & Kanahama, K. (2013). Root-zone Cooling at High Air Temperatures Enhances Physiological Activities and Internal Structures of Roots in Young Tomato Plants. Journal of the Japanese Society for Horticultural Science, 82, 322-327. [Google Scholar] [CrossRef]
[15]	Khaeim, H., Kende, Z., Jolánkai, M., Kovács, G. P., Gyuricza, C., & Tarnawa, Á. (2022). Impact of Temperature and Water on Seed Germination and Seedling Growth of Maize (Zea mays L.). Agronomy, 12, Article 397. [Google Scholar] [CrossRef]
[16]	Kuroda, A., & Sawada, Y. (2021). Effects of Temperature on Seed Dormancy and Germination of the Coastal Dune Plant Viola grayi: Germination Phenology and Responses to Winter Warming. American Journal of Botany, 109, 237-249. [Google Scholar] [CrossRef] [PubMed]
[17]	Lemmon, E., W. (2010). Thermophysical Properties of Air. In W. M. Haynes, & D. R. Lide (Eds.), CRC Handbook of Chemistry and Physics (pp. 615-637). CRC Press.
[18]	Lonati, M., Moot, D. J., Aceto, P., Cavallero, A., & Lucas, R. J. (2009). Thermal Time Requirements for Germination, Emergence and Seedling Development of Adventive Legume and Grass Species. New Zealand Journal of Agricultural Research, 52, 17-29. [Google Scholar] [CrossRef]
[19]	Lu, S., Ren, T., Gong, Y., & Horton, R. (2007). An Improved Model for Predicting Soil Thermal Conductivity from Water Content at Room Temperature. Soil Science Society of America Journal, 71, 8-14. [Google Scholar] [CrossRef]
[20]	Lu, Y., Lu, S., Horton, R., & Ren, T. (2014). An Empirical Model for Estimating Soil Thermal Conductivity from Texture, Water Content, and Bulk Density. Soil Science Society of America Journal, 78, 1859-1868. [Google Scholar] [CrossRef]
[21]	Nakano, Y. (2007). Response of Tomato Root Systems to Environmental Stress under Soilless Culture. Japan Agricultural Research Quarterly, 41, 7-15. [Google Scholar] [CrossRef]
[22]	Naranjo-Mendoza, C., Wright, A. J., Oyinlola, M. A., & Greenough, R. M. (2018). A Comparison of Analytical and Numerical Model Predictions of Shallow Soil Temperature Variation with Experimental Measurements. Geothermics, 76, 38-49. [Google Scholar] [CrossRef]
[23]	National Astronomical Observatory of Japan (2016). Rika nenpyo (Chronological Scientific Tables 2017). Maruzen Publishing Ltd.
[24]	Ramires, M. L. V., Nieto de Castro, C. A., Nagasaka, Y., Nagashima, A., Assael, M. J., & Wakeham, W. A. (1995). Standard Reference Data for the Thermal Conductivity of Water. Journal of Physical and Chemical Reference Data, 24, 1377-1381. [Google Scholar] [CrossRef]
[25]	Sakamoto, M., & Suzuki, T. (2015). Elevated Root-Zone Temperature Modulates Growth and Quality of Hydroponically Grown Carrots. Agricultural Sciences, 6, 749-757. [Google Scholar] [CrossRef]
[26]	Sampayo-Maldonado, S., Cabrera-Santos, D., Dávila-Aranda, P., Rodríguez-Arévalo, N. I., Orozco-Segovia, A., Gianella, M. et al. (2025). Using the Optimal Seed Germination Temperature Approach to Determine the Potential Distribution of Inga Jinicuil in Mexico under Climate Change Scenarios. Scientific Reports, 15, Article No. 3951. [Google Scholar] [CrossRef] [PubMed]
[27]	Sass, J. H., Lachenbruch, A. H., & Munroe, R. J. (1971). Thermal Conductivity of Rocks from Measurements on Fragments and Its Application to Heat-Flow Determinations. Journal of Geophysical Research, 76, 3391-3401. [Google Scholar] [CrossRef]
[28]	Shishido, Y., & Kumakura, H. (1994). Effects of Root Temperature on Photosynthesis, Transpiration, Translocation and Distribution of 14c-Photoassimilates and Root Respiration in Tomato. Journal of the Japanese Society for Horticultural Science, 63, 81-89. [Google Scholar] [CrossRef]
[29]	Tarnawski, V. R., McCombie, M. L., Momose, T., Sakaguchi, I., & Leong, W. H. (2013). Thermal Conductivity of Standard Sands. Part III. Full Range of Saturation. International Journal of Thermophysics, 34, 1130-1147. [Google Scholar] [CrossRef]
[30]	Tarnawski, V. R., Momose, T., McCombie, M. L., & Leong, W. H. (2015). Canadian Field Soils III. Thermal-Conductivity Data and Modeling. International Journal of Thermophysics, 36, 119-156. [Google Scholar] [CrossRef]
[31]	Thunholm, B. (1990). A Comparison of Measured and Simulated Soil Temperature Using Air Temperature and Soil Surface Energy Balance as Boundary Conditions. Agricultural and Forest Meteorology, 53, 59-72. [Google Scholar] [CrossRef]
[32]	Trudgill, D. L., Squire, G. R., & Thompson, K. (2000). A Thermal Time Basis for Comparing the Germination Requirements of Some British Herbaceous Plants. New Phytologist, 145, 107-114. [Google Scholar] [CrossRef]
[33]	Wang, R., Isozaki, M., Iwasaki, Y., & Muramatsu, Y. (2022). Root-Zone Temperature Effects on Spinach Biomass Production Using a Nutrient Film Technique System. HortScience, 57, 532-540. [Google Scholar] [CrossRef]
[34]	Wang, X., Zhang, W., Miao, Y., & Gao, L. (2016). Root-Zone Warming Differently Benefits Mature and Newly Unfolded Leaves of Cucumis sativus L. Seedlings under Sub-Optimal Temperature Stress. PLOS ONE, 11, e0155298. [Google Scholar] [CrossRef] [PubMed]
[35]	Woodside, W., & Messmer, J. H. (1961). Thermal Conductivity of Porous Media. I. Unconsolidated Sands. Journal of Applied Physics, 32, 1688-1699. [Google Scholar] [CrossRef]
[36]	Wu, Y., Zhang, H., Tian, Y., Song, Y., & Li, Q. (2023). A Thermal Time Basis for Comparing the Germination Requirements of Alfalfa Cultivars with Different Fall Dormancy Ratings. Agronomy, 13, Article 2969. [Google Scholar] [CrossRef]
[37]	Xu, Q., & Huang, B., (2001). Lowering Soil Temperature under Supraoptimal Air Temperature Improved Shoot and Root Growth in Creeping Bent Grass. Crop Science, 41, 1878-1883.
[38]	Yamori, N., Levine, C. P., Mattson, N. S., & Yamori, W. (2022). Optimum Root Zone Temperature of Photosynthesis and Plant Growth Depends on Air Temperature in Lettuce Plants. Plant Molecular Biology, 110, 385-395. [Google Scholar] [CrossRef] [PubMed]