Simulation of Daily Mean Soil Temperatures for Agricultural Land Use Considering Limited Input Data

Abstract

A one-dimensional simulation model that simulates daily mean soil temperature on a daily time-step basis, named AGRISOTES (AGRIcultural SOil TEmperature Simulation), is described. It considers ground coverage by biomass or a snow layer and accounts for the freeze/thaw effect of soil water. The model is designed for use on agricultural land with limited (and mostly easily available) input data, for estimating soil temperature spatial patterns, for single sites (as a stand-alone version), or in context with agrometeorological and agronomic models. The calibration and validation of the model are carried out on measured soil temperatures in experimental fields and other measurement sites with various climates, agricultural land uses and soil conditions in Europe. The model validation shows good results, but they are determined strongly by the quality and representativeness of the measured or estimated input parameters to which the model is most sensitive, particularly soil cover dynamics (biomass and snow cover), soil pore volume, soil texture and water content over the soil column.

1. Introduction

Soil temperature plays an important role in many soil processes and is related to atmospheric, soil and surface conditions. Temperature as one of the driving factors of soil genesis was recognized in the late 19th century by Dokuchaev [1] and independently by Hilgard [2], who listed climate, plants and organisms, the parent material and time as key soil-forming factors. Later, for example, Ellenberg [3] and Gray et al. [4] showed there is a dominant influence of climate and parent material on numerous soil properties. The soil properties, in turn, influence the water and energy transfer, which directly affect the soil temperature, then directly determine the activity levels and survival of soil fauna and flora [5]. Soil temperature measurements are critical for calibrating many soil temperature response functions in simulation models [6], since soil temperature and water content affect physical, chemical and soil biological processes [7,8]. A good performance of soil temperature simulation is therefore very important, especially for crop models, which include nutrient balance or models simulating leaching processes and gas emissions from soils. For example, the soil carbon (C) and nitrogen (N) balance is determined by soil temperature and water content [9,10,11], because they regulate the rate of N-mineralization and, further, the emission of gases from soils [12,13]; C/N balance models are a critical part of most crop growth models, e.g., [14,15]. Many related phenomena are related closely to soil temperatures, such as the effects on soil properties [16], the soil carbon balance under climate change [17], etc. As soil properties and surface conditions vary significantly over space, this creates a big challenge for representative site-specific soil temperature simulations [18,19] as well as for spatial applications of soil temperature models, as shown by, e.g., [6]. Most of the data on near-ground atmospheric conditions are available from in situ weather stations, which are irregularly distributed, creating substantial uncertainties about the spatial distribution, especially for precipitation. However, soil properties and important surface conditions, such as snow cover and plant canopy characteristics [20,21] that influence soil temperatures, are much more difficult to estimate due to their even stronger spatial variability.

Regarding agricultural land use, management factors such as irrigation, soil cultivation and crop management (i.e., harvest, cutting, and mulching) influence soil physical properties and/or soil surface characteristics, affecting soil temperature and wetness [22,23,24] and can add significant spatial and temporal changes compared to natural or undisturbed land. In conclusion, agricultural management strongly affects the dynamics of soil heat flux and needs to be considered in soil temperature modeling to get the most accurate results for soil processes such as soil respiration [25,26], decomposition rates [27,28], and the microbiological or biotic activity [29,30] of arable lands—all affecting final crop productivity [31,32].

Impact of seasonal snow cover and soil freezing on soil temperature is well documented [33,34,35] but not sufficiently captured by soil-vegetation-atmosphere-transfer (SVAT) models and most crop models. Neglecting snow cover and freezing effects can lead to significant error in the assessment of the soil temperature’s impact on winter crops’ development in crop models or the calculation of canopy micrometeorological conditions (canopy air temperature and humidity, leaf temperature and wetness, for example), commonly performed by SVAT models. Currently, numerical weather prediction (NWP) models and climate models [36] implement treatment of snow presence based on the energy balance of the surface (surface albedo) [37], but horizontal spatial grid spacing of these models (particularly the vegetation/land use maps) cannot satisfy the needs for simulations at the field or farm level [38]. However, NWP model outputs can be used as input weather data for more sophisticated soil temperature models, which can operate as standalones or as part of SVAT, crop, or other agrometeorological models.

Inventory of different soil temperature (sub)models of SVAT models offer comprehensive modeling of physical processes by describing soil surface-atmosphere energy transfer and simulating the soil temperature on different levels [39,40,41]. Due to the time scales of the driving processes (surface exchange of short-wave and long-wave radiation), hourly meteorological data are used. On the other hand, for the process-oriented models used in crop production (such as crop growth and yield models, pest models, irrigation models, agrometeorological models) and for site-specific or spatial applications, only limited sets of hourly input data are available [42]. However, these are not necessary because a daily simulation time step is sufficient for most of these applications [43]. For example, the daily course of temperatures can be assessed from daily extremes of air temperatures. In regard to soil conditions, short-term variations in precipitation, air temperature, or radiation and their impact on soil temperature and wetness are buffered in the soil column, averaging out short-term variations that could affect simulated ecological processes.

Therefore, the main objective of our study is the development and validation of a robust and tailored soil temperature model for applications with limited input data requirements (including daily time step of weather data and surface conditions as well as easily available or measurable soil physical data). In this context, reduced uncertainty can be achieved by an optimized balance between model complexity and input data requirements. Furthermore, the model is designed for low computational requirements, allowing for applications for decision-making in crop production as well as applications over larger scales (e.g., using GIS).

The AGRIcultural SOil TEmperature Simulation (AGRISOTES) model developed in our study, is based on model approaches with limited complexity [44,45,46] in order to calculate daily soil temperatures on different flexible depths using daily weather data and temporal changes of soil surface characteristics as well as soil pore volume, which allows for a consideration of the soil cultivation effects. Considering the specific needs of agricultural applications (such as the living conditions for soil-borne pests), particular attention is further devoted to the parameterization of processes in the presence of soil cover (snow and vegetation) and ice in the soil. AGRISOTES is designed as a tool for a wide range of possible combinations with other ecosystem models or algorithms, measured datasets as well as model implementations.

2. Materials and Methods

2.1. Model Description

In this section, the developed AGRISOTES model will be described. The main model inputs addressing variable temporal conditions (Figure 1) on a daily time step are meteorological elements (mean, maximum and minimum dry-bulb air temperature, solar radiation, and actual evapotranspiration), soil characteristics in order to address soil cultivation events (soil-layer-specific pore volume and soil water content), and soil cover characteristics (snow cover and soil surface biomass), which will allow us to simulate sudden changes of snow or crop cover. Static input parameters include soil texture and configuration parameters, such as simulated soil depth, number of layers and layer thickness, etc. (Figure 1). In general, all input data can be either based on measured data or estimated through other models or simplified algorithms. However, the final selection of input data sources depends on the accuracy and representativeness for the specific application. Our model input scheme allows for maximum flexibility in model application, but it is up to the decision of the model users if a specific type of dataset is suitable for a relevant application. The notations and units of input/output or processed parameters can be found in the Appendix D.

The model consists of two parts: (i) the soil surface temperature module and (ii) the subsurface soil temperature module, which includes heat transport and freezing and thawing of soil water (Figure 2). These two modules are coupled with the ground heat flow at the soil surface. In the soil surface temperature module of AGRISOTES, with respect to snow cover, three different conditions will be described: without snow cover, with dense snow cover and with partial snow cover (limits to be defined as calibrated input).

2.1.1. Soil Surface Temperature: Vegetated Soil Surface without Snow Cover

The calculation of soil surface temperature is based on Best’s canopy model [47], which is a simplification of the model proposed by Deardorff (1978), [48]. The plant canopy is treated as a single foliage layer covering a fraction ν_f of the ground surface area. The vegetation cover (or fractional cover) ν_f must be between 0 and 1 and has to be chosen to be representative of the type and amount of aboveground biomass CV (kg ha⁻¹).

The heat balance equation for the soil surface is formulated in the AGRISOTES model in a similar way to that in Herb et al. (2006) [49], as follows:

$$h_{net,g}=\left(1-{\nu }_f\right)\left(\left(1-{\alpha }_g\right)h_s+{\epsilon }_g{h}_{l,ac}\right)+h_{l,fg}-h_{l,ga}-h_{evap,g}-h_{conv,g}$$

(1)

where h_net,g is the net ground heat flow (W m⁻²), α_g is the soil surface albedo, h_s is the observed solar radiation per unit area of the ground surface (W m⁻²), ε_g is the soil emissivity, h_l,ac is the incoming atmospheric long-wave radiation (W m⁻²), h_l,fg is the net long-wave radiation from the foliage to the ground (W m⁻²), h_l,ga is the outgoing long-wave radiation emitted from the bare soil (W m⁻²), h_evap,g is the evaporative heat transfer from the soil surface to the atmosphere (W m⁻²) and h_conv,g is convective heat transfer from the soil surface to the atmosphere (W m⁻²).

The convective and evaporative heat flux components between the ground and the canopy are neglected, as in Best’s model. As indicated by [49], this simplification is only suitable in the presence of a dense canopy. Therefore, partial canopies are treated here by linearly combining the dense canopy model and the model for bare soil using the vegetation cover factor ν_f in a similar manner as that described by [49].

The soil surface albedo α_g and soil emissivity ε_g are calculated from the surface wetness, as proposed by Van Bavel and Hillel (1976) [50]. The incoming atmospheric long-wave radiation, h_l,ac, is according to [51], expressed as follows:

$$h_{l,ac}={\epsilon }_{ac}\sigma {\left(T_a+273.15\right)}^4$$

(2)

where T_a is the air temperature (°C), σ is the Stefan-Boltzmann constant (≈5.67 × 10⁻⁸ W m⁻² K⁻⁴), and ε_ac is the atmospheric emissivity (Appendix A).

For the calculation of h_l,fg, the soil surface and vegetation layer are considered as two parallel plates for which the net radiation exchange is calculated and scaled with vegetation cover ν_f. Following [52], using air temperature in the vicinity of plants as an approximation for canopy temperature leads to an expression for h_l,fg in the form

$$h_{l,fg}={\nu }_f{\mu }_{fg}\sigma \left({\left(T_a+273.15\right)}^4-{\left(T_g+273.15\right)}^4\right)$$

(3)

where T_g is the soil surface temperature (°C) and μ_fg is the degree of radiation exchange defined as [53]:

$${\mu }_{fg}=\frac{1}{1/{\epsilon }_f+1/{\epsilon }_g-1},$$

(4)

where ε_f is the foliage emissivity. The outgoing long-wave radiation emitted from the bare soil h_l,ga is calculated as

$$h_{l,ga}=\left(1-{\nu }_f\right){\epsilon }_g\sigma {\left(T_g+273.15\right)}^4.$$

(5)

The evaporative heat transfer h_evap,g from the ground to the atmosphere is roughly approximated to be a fraction (1 − ν_f) of the total evaporative heat transfer (ground + foliage) h_evap (W m⁻²):

$$h_{evap,g}\approx \left(1-{\nu }_f\right)h_{evap}.$$

(6)

It should be pointed out that Equation (6) is not based on physical principles, but is just a simple approximation using the vegetation cover ν_f to interpolate linearly between two limit cases: the case of bare soil, in which ν_f = 0 and h_evap,g = h_evap, and the case of full vegetation cover, in which ν_f = 1 and the assumption is made that soil evaporation can be neglected (h_evap,g = 0). This approximation allows the AGRISOTES model to perform its calculations when only actual daily evapotranspiration is provided as input, without knowledge of the splitting up between ground and foliage of the total evaporative heat transfer. In contrast to this, in applications where daily ground evaporation is available as an input, the model accuracy could probably be improved by keeping h_evap,g as a direct model input instead of using the approximation in Equation (6). However, this will, in general, require recalibration of the model, so that the calibration results presented afterwards no longer remain valid.

The convective heat flow from the soil surface to the atmosphere is assumed to occur only from the uncovered fraction of the soil surface (for the covered fraction, convective heat flow is assumed to occur only between the foliage and air) and is expressed by a ground convective heat transfer h_conv,g in the form

$$h_{conv,g}=\left(1-{\nu }_f\right)k_{cg}\left(T_g-T_a\right)$$

(7)

where k_cg is the convective heat transfer coefficient between soil surface and air (W m⁻² K⁻¹), which is dependent on the bare soil surface roughness length and the wind speed. In the model, k_cg is considered a direct model input parameter and must either be chosen to be representative of the average wind speed and surface characteristics in the considered region or be determined by model calibration. Combining Equations (1)–(3) and (5)–(7), the heat balance equation for the vegetated soil surface can be written in the form:

$$\begin{array}{ll}{h}_{net,g}=& \left(1-{\nu }_f\right)\left(\left(1-{\alpha }_g\right)h_s+{\epsilon }_g\sigma \left({\epsilon }_{ac}{\left(T_a+273.15\right)}^4-{\left(T_g+273.15\right)}^4\right)-h_{evap}-k_{cg}\left(T_g-T_a\right)\right)+\\ & +{\nu }_f{\mu }_{fg}\sigma \left({\left(T_a+273.15\right)}^4-{\left(T_g+273.15\right)}^4\right).\end{array}$$

(8)

Equation (8) is highly nonlinear with respect to temperature. To calculate T_g from Equation (8), it is linearized using a first-order Taylor series approximation:

$$\sigma {\left(T_g+273.15\right)}^4\approx \sigma {\left(T_a+273.15\right)}^4+4\sigma {\left(T_a+273.15\right)}^3\Delta T_g=h_{l0}+k_{l0}\Delta T_g,$$

(9)

where the following abbreviations are introduced:

$$\Delta T_g=T_g-T_a, h_{l0}=\sigma {\left(T_a+273.15\right)}^4, k_{l0}=4\sigma {\left(T_a+273.15\right)}^3, h_{l,g0}={\epsilon }_g{h}_{l0}, k_{l,g0}={\epsilon }_g{k}_{l0}.$$

This transforms the heat balance equation into the form

$$h_{net,g}=\left(1-{\nu }_f\right)\left(\left(1-{\alpha }_g\right)h_s-\left(1-{\epsilon }_{ac}\right)h_{l,g0}-h_{evap}\right)-\left(\left(1-{\nu }_f\right)\left(k_{l,g0}+k_{cg}\right)+{\nu }_f{\mu }_{fg}{k}_{l0}\right)\Delta T_g.$$

(10)

Equation (10) establishes a relationship between the output variables T_g and h_net,g and the input parameters of the model. It provides a basis to derive a top boundary condition for the simulation of soil temperatures within the soil column during snow-free periods. To aggregate the input model variables, the following abbreviations are introduced:

$$k_1=\left(1-{\nu }_f\right)\left(k_{l,g0}+k_{cg}\right)+{\nu }_f{\mu }_{fg}{k}_{l0},$$

(11)

$$h_{t,1}=\left(1-{\nu }_f\right) \left(\left(1-{\alpha }_g\right)h_s-\left(1-{\epsilon }_{ac}\right)h_{l,g0}-h_{evap}\right).$$

(12)

So, Equation (10) can be rewritten as

$$h_{net,g}=h_{t,1}-k_1\Delta T_g,$$

(13)

or in the form

$$\Delta T_g=\Delta T_1-R_1{h}_{net,g},$$

(14)

where $R_1=1/k_1$ and $\Delta T_1=h_{t,1}/k_1$.

To calculate the daily mean soil surface temperature T_g,mean, a daily time averaging procedure is applied to Equation (13) (see in Appendix B) in order to obtain the difference between the daily mean soil surface temperature T_g,mean and the daily mean air temperature T_a,mean in the form

$$\Delta T_{g,mean}=\Delta T_{1,mean}-R_{1,mean}{h}_{net,g,mean},$$

(15)

which will be used as the top boundary condition for the simulation of soil temperature during snow-free periods.

2.1.2. Soil Surface Temperature: Soil Surface with Dense Snow Cover

Dense snow cover affects temperature and soil moisture on a daily basis. In the AGRISOTES model, the calculation of the daily mean soil surface temperature T_g,mean with dense snow cover, is implemented in a similar way to that of the soil surface temperature module of the DAYCENT model [54] as

$$T_{g,mean}=\left\{\begin{array}{ll}0{ }^{°}\mathrm{C}\hfill & \mathrm{for} T_{a,mean}\ge 0{ }^{°}\mathrm{C},\hfill \\ 0.3\cdot T_{a,mean}\cdot k_{snow}\hfill & \mathrm{for} T_{a,mean}<0{ }^{°}\mathrm{C},\hfill \end{array}\right.$$

(16)

where k_snow is defined as

$$k_{snow}=\max \left(-0.015\cdot SNO+1.0, 0\right)$$

(17)

and SNO is the snow water equivalent of the snow layer (mm). By subtracting T_a,mean from Equation (16), it may also be written in the form of Equation (15):

$$\Delta T_{g,mean}=\Delta T_{2,mean}-R_{2,mean}\cdot h_{net,g,mean},$$

(18)

where

$$R_{2,mean}=0 {\mathrm{m}}^2{ \mathrm{K} \mathrm{W}}^{-1}$$

(19)

and

$$\Delta T_{2,mean}=\left\{\begin{array}{ll}-T_{a,mean}\hfill & \mathrm{for} T_{a,mean}\ge 0{ }^{°}\mathrm{C},\hfill \\ 0.3\cdot T_{a,mean}\cdot k_{snow}-T_{a,mean}\hfill & \mathrm{for} T_{a,mean}<0{ }^{°}\mathrm{C}.\hfill \end{array}\right.$$

(20)

2.1.3. Soil Surface Temperature: Soil Surface with Partial Snow Cover (General Case)

For a given SNO value, which is designated as an input to the AGRISOTES model, the snow coverage factor ν_sn is calculated as

$${\nu }_{sn}=\left\{\begin{array}{ll}0\hfill & \mathrm{for} SNO\le SN{O}_{limit,1},\hfill \\ \frac{SNO-SN{O}_{limit,1}}{SN{O}_{limit,2}-SN{O}_{limit,1}}\hfill & \mathrm{for} SN{O}_{limit,1}<SNO<SN{O}_{\mathrm{l}imit,2},\hfill \\ 1\hfill & \mathrm{for} SNO\ge SN{O}_{limit,2},\hfill \end{array}\right.$$

(21)

where SNO_limit,1 is the lower limit of SNO below which the snow layer is neglected, and SNO_limit,2 is the upper limit of SNO beyond which snow cover is treated as dense (completely covered soil surface).

For partial snow cover, the soil surface is considered to be a surface consisting of patches of vegetated soil surface and soil with dense snow cover. The ground heat flow h_net,g,mean is assumed to be homogeneous with respect to the horizontal surface, whereas the soil surface temperature difference ΔT_g,mean is assumed to vary according to the local snow cover. Therefore, to obtain an approximation for a horizontally averaged difference in soil surface temperature, a linear combination of Equations (15) and (18) can be used in the form

$$\Delta T_{g,mean}=\left(1-{\nu }_{sn}\right)\left(\Delta T_{1,mean}-R_{1,mean}{h}_{net,g,mean}\right)+{\nu }_{sn}\left(\Delta T_{2,mean}-R_{2,mean}{h}_{net,g,mean}\right)$$

so that the definitions

$$\Delta T_{mean}=\left(1-{\nu }_{sn}\right)\Delta T_{1,mean}+{\nu }_{sn}\Delta T_{2,mean}$$

(22)

$$R_{mean}=\left(1-{\nu }_{sn}\right)R_{1,mean}+{\nu }_{sn}{R}_{2,mean}$$

(23)

can be used to obtain the following equation, which will be used as top boundary condition for soil temperature simulation when snow cover is present:

$$\Delta T_{g,mean}=\Delta T_{mean}-R_{mean}\cdot h_{net,g,mean}.$$

(24)

2.1.4. Subsurface Soil Temperature: Soil Temperature Equation in Finite-Difference Form

The soil temperatures at different depths are calculated by spatial and temporal discretization of the differential equations for a soil heat flow. The simulated layers are aligned to soil horizon boundaries and thin near the soil surface, with gradually increasing thickness towards the bottom boundary [49]. We applied Fourier’s law of heat conduction (Equation (25)) and the energy conservation law in differential form (Equation (26)), including a term accounting for the latent heat of the freezing or thawing of soil water, both tailored for a one-dimensional heat flow:

$$h\left(z,t\right)=-\lambda \left(z,t\right)\frac{\partial T_s\left(z,t\right)}{\partial z},$$

(25)

$$C\left(z,t\right)\frac{\partial T_s\left(z,t\right)}{\partial t}-F_{H2O}\frac{\partial {\theta }_{ice}\left(z,t\right)}{\partial t}=-\frac{\partial h\left(z,t\right)}{\partial z}.$$

(26)

Here, h is the soil heat flow (W m⁻²), T_s is the soil temperature (°C), λ is the thermal conductance of soil (W m⁻¹ K⁻¹), C is the volumetric heat capacity of soil (J m⁻³ K⁻¹), F_H2O is the latent heat of ice fusion related to the volume of liquid water (≈3.34 × 10⁸ J m⁻³), θ_ice is the volumetric ice content (cm³ cm⁻³) in the soil (neglecting volume expansion of frozen water), z is the depth (m) and t is time (s).

The equations above assume negligible water migration, which can lead to short-term deviations, e.g., after heavy precipitation and on strongly draining soils. An implicit, semidiscrete finite difference scheme to solve these partial differential equations approximately is derived. The simulation domain for the depth z in the soil profile is defined as [0, z_max], where z_max is the depth of the bottom of the simulated soil profile (m). The simulation time step Δt is one day (Δt = 86400 s), while every time point t_j is defined as t_j = jΔt. The evaluation of quantities at time point t_j will further be denoted by an upper index j.

In the further derivations, time averaging will be used in the form ${\widehat{\overline{C}}}^j\left(z\right)$, where the bar ( ¯ ) indicates the averaging of variable C over the 24-h period [t, t + Δt], with t = t_j in this example, while the caret (ˆ) indicates approximated quantities. If the soil thermal properties λ and C as well as the time-averaged heat flow gradient $\overline{{\partial }_{z}h}$ do not vary significantly during a 24-h period and ∂_zh is continuous, by applying daily averaging on Equation (25) once and on Equation (26) twice, it can be concluded that the following semidiscrete finite difference equations hold approximately:

$${\widehat{\overline{h}}}^j\left(z\right)=-{\widehat{\overline{\lambda }}}^j\left(z\right)\cdot {\partial }_z{\overline{T}}_s^j\left(z\right),$$

(27)

$${\widehat{\overline{C}}}^j\left(z\right)\frac{{\overline{T}}_s^j\left(z\right)-{\overline{T}}_s^{j-1}\left(z\right)}{\Delta t}-F_{H2O}\frac{{\overline{\theta }}_{ice}^j\left(z\right)-{\overline{\theta }}_{ice}^{j-1}\left(z\right)}{\Delta t}=-{\partial }_z{\widehat{\overline{h}}}^j\left(z\right).$$

(28)

In Equations (27) and (28), values ${\widehat{\overline{C}}}^j\left(z\right)$ and ${\widehat{\overline{\lambda }}}^j\left(z\right)$ are approximations for ${\overline{C}}^j\left(z\right)$ and ${\overline{\lambda }}^j\left(z\right)$, defined as

$${\widehat{\overline{C}}}^j\left(z\right)=C\left(V{P}_{sand}\left(z\right),V{P}_{clay}\left(z\right),{\overline{\theta }}_{sat}^j\left(z\right),{\overline{\theta }}^j\left(z\right),{\widehat{\overline{\zeta }}}^{j-1}\left(z\right)\right)$$

(29)

and

$${\widehat{\overline{\lambda }}}^j\left(z\right)=\lambda \left(V{P}_{sand}\left(z\right),V{P}_{clay}\left(z\right),{\overline{\theta }}_{sat}^j\left(z\right),{\overline{\theta }}^j\left(z\right),{\widehat{\overline{\zeta }}}^{j-1}\left(z\right)\right),$$

(30)

where $\widehat{\overline{\zeta }}$ is the approximate daily mean ice fraction, defined as

$${\widehat{\overline{\zeta }}}^j\left(z\right)=\frac{{\overline{\theta }}_{ice}^j\left(z\right)}{{\overline{\theta }}^j\left(z\right)},$$

(31)

with θ denoting the volumetric water content (cm³ cm⁻³), while the five-parametric C(VP_sand, VP_clay, θ_sat, θ, ζ) and λ(VP_sand, VP_clay, θ_sat, θ, ζ) functions denote the soil thermal property calculation procedures, which are additionally dependent on the volume percentages of clay VP_clay (%) and sand VP_sand (%) relative to the total solid matter and the saturated volumetric water content θ_sat [55]. To calculate the thermal properties of the soil, the method described for the model SWAP version 3.2 [55], which is based on [56,57,58,59,60] was used, together with some enhancements to treat the case of partially or totally frozen soil.

The application of the finite difference method to solve partial differential equations generally requires sufficient differentiability of the involved space- and time-dependent quantities. Therefore, regularization techniques have been used for the numerical solution of Equations (27) and (28). In particular, the ideally discontinuous dependence of ζ on T_s (see Section 2.1.6) has been replaced by a smoothed function that does not differ significantly from the idealized behavior but is sufficiently often differentiable to allow the application of the finite difference method for the spatial discretization of Equations (27) and (28). Spatial discontinuities of soil texture properties at the boundary between soil horizons, as long as they are not moving and their depths are known, are allowed in our model and have been appropriately accounted for in the numerical implementation.

2.1.5. Subsurface Soil Temperature: Top and Bottom Boundary Conditions

To formulate the top boundary condition for the finite difference scheme approximation of soil heat flow, the soil surface heat balance (Equation (24)) must be incorporated, which is approximated for the finite-difference scheme using the formulation

$${\overline{T}}_s^j\left(0\right)=T_{a,mean}^j+\Delta T_{mean}^j-R_{mean}^j\cdot {\widehat{\overline{h}}}^j\left(0\right).$$

(32)

For the bottom boundary condition, three different modes are available for the soil temperature model:

(i)

To set the approximated daily mean heat flux ${\widehat{\overline{h}}}^j\left(z_{max}\right)$ at the bottom of the simulated soil profile to zero;

(ii)

To fix the temperature at the bottom of the soil profile to the annual mean air temperature T_AA (°C), which is an input to the model; and

(iii)

To assume that the soil thermal properties λ and C and the ice content θ_ice near the bottom of the soil profile are constant over z and t for depths z ≥ z_max and that the temperature at the bottom of the soil profile varies sinusoidally with a one-year period around a mean value of T_AA. With the angular frequency of a year ω_a (rad s⁻¹) and the approximated damping depth at the bottom of the soil profile $d_a^j\left(z_{max}\right)$ (m) for annual frequency, this leads to the following bottom boundary condition (detailed derivation can be found in Appendix C):

$$-{\partial }_z{\overline{T}}_s^j\left(z_{max}\right)=\frac{1}{d_a^j\left(z_{max}\right)}\left(\left(1-\tan \left(\frac{{\omega }_a\Delta t}{2}\right)\right)\left({\overline{T}}_s^j\left(z_{max}\right)-T_{AA}\right)+\frac{{\overline{T}}_s^j\left(z_{max}\right)-{\overline{T}}_s^{j-1}\left(z_{max}\right)}{\sin \left({\omega }_a\Delta t\right)}\right).$$

(33)

Simulations at site Groß-Enzersdorf, which is located northeast of Vienna, Austria (see Section 2.2), representing a Central European semihumid climate, have shown that the first two modes for bottom boundary conditions were only applicable to depths greater than approximately 10 m. The third mode for the bottom boundary condition was applicable even at lower depths of approximately 2.5 m. Therefore, when simulation results of the soil temperature at depths between 2.5 m and 10 m are not required, the third bottom boundary condition mode can be used to reduce the computing time of the soil temperature simulation.

2.1.6. Subsurface Soil Temperature: Treatment of Freezing/Thawing of Soil Water

Figure 3 shows an idealized graph for some fixed depth z (neglecting freezing-point depression) of the dependence of soil temperature and ice content on heat energy density w (J m⁻³). The zero-point of w is defined as the state where T_s = 0 °C and θ_ice = 0. It can be seen what occurs if the soil at first has a temperature T_s,1 and ice content θ_ice,1 and the heat energy density changes by Δw, so that afterwards, the soil has temperature T_s,2 and ice content θ_ice,2.

As can be seen, during time periods in which the soil is partially frozen (0 < θ_ice < θ), the temperature stays constant at 0 °C. For T_s > 0 °C, the soil is unfrozen, so θ_ice = 0 and the derivative of T_s with respect to w is given as 1/C_unfrz, where C_unfrz is the volumetric heat capacity for unfrozen soil. For T_s < 0 °C, the soil is completely frozen, so θ_ice = θ and the derivative of T_s with respect to w is given as 1/C_frz, where C_frz is the volumetric heat capacity for completely frozen soil. It should be noted that the shape of the curves in Figure 3 depends on the soil water content θ, as well as on soil composition and pore volume.

Considering the daily mean soil temperature ${\overline{T}}_s$ and ice content ${\overline{\theta }}_{ice}$, their behavior is comparable to, but not the same as that shown in Figure 3 for T_s and θ_ice. Due to the diurnal variation of T_s and θ_ice, on some days ${\overline{T}}_s$ > 0 °C, but during some night hours, freezing may occur so that, nevertheless, ${\overline{\theta }}_{ice}$ > 0. Since it is not possible to take into account precisely such effects in a model that is only dealing with daily averaged quantities, the following assumptions are introduced:

$${\overline{T}}_s^j\left(z\right)>0{ }^{°}\mathrm{C}\Rightarrow {\overline{\theta }}_{ice}^j\left(z\right)=0,$$

(34)

$${\overline{T}}_s^j\left(z\right)<0{ }^{°}\mathrm{C}\Rightarrow {\overline{\theta }}_{ice}^j\left(z\right)={\overline{\theta }}^j\left(z\right).$$

(35)

The above approach is an approximation and, therefore, in the simulation, particularly during those days on which soil water freezing starts or ends (considering some fixed depth z), greater deviations of simulated daily mean soil temperature and ice content from the real values can occur (see validation results in Figure 6). As has already been mentioned in Section 2.1.4, the discontinuous dependence of $\widehat{\overline{\zeta }}$ on ${\overline{T}}_s$ resulting from Equations (34) and (35) was replaced by a similar smoothed function for the numerical implementation of the model. In summary, the applied simplifications as described above can lead to short-term deviations, especially due to sudden changes in surface conditions (e.g., timing of snow or vegetation cover), extreme weather impacts such as heavy precipitation, which can influence additionally heat transport in soils, and deviations of the freezing status of soil layers.

2.2. Calibration and Validation Sites

AGRISOTES was calibrated and validated at different sites with various climates, soils and agricultural land uses in Austria, Belgium, Turkey and the Czech Republic (

Table 1. Main characteristics of site conditions used for generating representative inputs for the AGRIcultural SOil TEmperature Simulation (AGRISOTES) model for the calibration and validation periods (all sites characterized by no or only negligible slope inclination and no shading from larger structures).

Site	Soil and Surface Management	Station Type	Annual Means of Temperature and Precipitation (1981–2010)	Main Soil Type	AGRISOTES Model Inputs M: Measured at Site E: Indirectly Estimated S: Simulated from Measured Weather Data
					Met 1	Soil 2	θ3	CV4
Groß-Enzersdorf Northeast Austria (Calibration)	Arable – Bare soil and straw mulch	Field experiment station	10.3 °C – 516 mm	Para-chernosem	M	E/M 5	M	M
Goggendorf Northeast Austria (Validation)	Arable – Bare soil and fallow after winter wheat	Field experiment station	9.2 °C – 519 mm	Chernosem	M	M	S	E
Doksany Northwest Czech Republic (Validation)	Permanent short grass	Long term field experiment station	9.5 °C – 466 mm	Chernosem	M	M	M	E
Pucking Northwest Austria (Validation)	Arable – Crop rotation	Lysimeter	9.1 °C – 790 mm	Cambisol	M	M	S	E
Kirklareli, Thrace region, Turkey (Validation)	Winter wheat/Bare soil and fallow after winter wheat	Field experiment station	13.3 °C – 570 mm	Vertisol Cambisol	M	M	M	M
Hamme, Belgium (Validation)	Permanent short grass	Field experiment station	10.6 °C – 852 mm	Eutric Fluvic Gleyic Cambisol	M	M	E	E
Purbach East Austria (Validation)	Arable – Vineyard	Field experiment station	10.4 °C – 758 mm	Chernosem	M	M	S	E
Grafendorf Southeast Austria (Calibration of snow cover effect/Validation)	Permanent grassland (3–cut)	Field experiment station	9.4 °C – 716 mm	Cambisol	M	E/M 6	M	E
Obersiebenbrunn Northeast Austria (Validation)	Natural grassland		10.3 °C – 516 mm	Chernosem	M	M	S	E

¹ Meteorological input parameters: Daily mean, maximum and minimum temperature, and solar radiation (global radiation), ² Soil: Sand and clay content (VP_sand, VP_clay) and total soil pore volume (measurement samples from up to four soil depths), ³ θ: Volumetric soil water content (measurements based on 1–2 soil depths), ⁴ CV: Total above ground dry biomass, ⁵ Soil texture estimated from soil map (mean value over soil profile), ⁶ Estimated for calibration of snow cover effect, but measured for validation.

). Table 1 presents an overview and comparison of the primary sites, from which measured daily soil temperature time series were available. Site conditions and measured/estimated/simulated model input parameters are indicated. Further details on site conditions are presented with the results.

For all sites, the measured physical soil properties (texture and soil pore volume) were available from up to four depths over the vertical soil profile, which were used as direct inputs and assumed to represent the full soil column. As a daily input, soil pore volume at all soil layers was assumed to be constant over the validation periods due to lack of temporal data. In some cases, the soil water content was simulated (S) using the data of on-site or nearby weather stations by the FAO model approach [61]. In all cases, daily actual evapotranspiration was simulated for AGRISOTES input using the FAO model approach as well. The daily input of soil surface dry matter dynamics was estimated (E) at most validation sites based on reported biomass at harvest and according to linear interpolation between phenological development stages [61] or by direct measurements during the growth period (M).

3. Results and Discussion

AGRISOTES Model Calibration, Validation and Sensitivity Analysis on Daily Mean Soil Temperatures

A calibration of two important model parameters was performed based on measured soil temperatures in a field trial (Groß-Enzersdorf, Austria). Calibration was carried out for arable (plowed) land with bare soil and with different straw mulch coverage to estimate the dampening effect of total aboveground dry matter biomass, CV (

Table 2. Calibrated cover effect factor (v_f) for total aboveground dry matter biomass (CV; based on straw mulch).

CV (kg ha−1)	vf (Dimensionless)
0	0
2500	0.848
4000	0.933
5000	0.965

) as well as the convective heat transfer coefficient (k_cg) for bare soil (

Table 3. Model parameters based on the calibration sites.

Parameter	Unit	Calibrated Value	First Estimate (Value Used in Subsequent Simulations)
kcg	W m−2 K−1	42.0	41.0
β	ha kg−1	0.000729	0.000663
εf	1	0.95	0.95
SNOlimit,1	mm H2O	0.3	0.4
SNOlimit,2	mm H2O	13.8	13.8

).

Taking the form of Equation (2) in [62], a simple exponential curve was fitted for the data in Table 2. The resulting relationship was

$${\nu }_f\left(CV\right)=1-\exp \left(-\beta \cdot CV\right),$$

(36)

with β given in Table 3. This relationship can be used to interpolate for values of CV other than the calibrated ones, but care needs to be taken for very low values, i.e., less than 2000 kg ha⁻¹, because gaps can occur in canopies under these conditions. For other types of biomass coverage (e.g., living crops), different values of v_f (in relation to CV) may apply depending on the canopy structure; however, the validation showed good results in our case for grassland and wheat canopy (see

Table 4. Soil temperature model calibration/validation on study specific field experiments.

Location	Surface Conditions	Simulation Period	Available Days	Soil Depth (cm)	R	σo (°C)	σs (°C)	RMSE (°C)	RRMSE (%)	IA (0–1)
Groß-Enzersdorf Calibration	4 × 4 m plot bare soil after seedbed preparation	1 April–24 May 2012	54	10	0.98	4.45	4.59	0.89	6.3	0.99
Goggendorf Validation	4 × 4 m plot bare soil in winter wheat stand since 22 June; soil water content simulated	22 June–9 December 2016	116	5	0.999	8.87	8.63	1.21	8.2	0.99
			116	10	0.999	8.79	8.31	1.17	7.9	0.99
			109	20	0.99	8.07	7.75	1.01	7.2	0.99
			117	40	0.999	6.73	6.66	0.59	4.1	0.99
	Winter wheat, 15 August harvest –stubble field – 20 September soil mulching– Soil water cont. simulated	22 June–9 December 2016	123	5	0.999	7.47	7.55	0.69	5.2	0.99
			123	20	0.999	6.71	6.77	0.55	4.1	0.99
			123	40	0.999	5.87	5.71	0.42	3.1	0.99

and Table 6).

The value for k_cg (convective heat transfer coefficient from the soil surface to air) was first optimized for bare soil with v_f set to zero (for the obtained value of k_cg see Table 3). Then, k_cg was fixed, and optimal values for v_f were estimated for the mulch-covered trial.

In Table 3, the calibrated model parameters are given. The value of ε_f was not optimized, but taken from [49]. The values of SNO_limit,1 and SNO_limit,2 were taken from parameter optimizations from a grassland site with regular snow cover occurrence (Grafendorf), since the calibration at the other site (Groß-Enzersdorf) was performed only during snow-free periods.

The goal of the optimizations was to reduce the root-mean-square-difference between simulated and measured daily mean soil temperatures near the soil surface. However, the presented calibration result can be understood as a first guess only, because our available calibration data sets are based on field measurements, which are related to inherent uncertainties, especially due to the nature of inhomogeneous physical properties of soils. Extended model calibrations are therefore recommended for other different sites to ensure the maximum accuracy of simulation results (where the number of measured parameters and measurement quality are considered as suitable, including soil physical parameters, daily LAI, canopy or snow cover gaps, and biomass measurements). In our case, a correction for the value of volumetric heat capacity of air has been made in the soil thermal property calculation routine of the model after a first parameter estimation. This caused a change in the calibrated model parameter values compared with their first estimates (Table 3), which were determined before the model correction. The parameter change, however, showed only very marginal effects on the simulation results. Therefore, all simulation results shown subsequently were generated using the corrected model version but with the first estimates of the model parameters since they were available earlier. The main limitation, however, is the limited availability of calibration data sets of suitable quality.

The quality of model calibration and validation was tested using: (a) R—the coefficient of determination between observed and simulated values; (b) σ—the standard deviation of observed (σ_o) and simulated values (σ_s); (c) RMSE—the root mean square error; (d) RRMSE—the relative RMSE and e) IA—the index of agreement [63]. According to [64], the simulation is performed more realistic if: (a) RMSE is less than the standard deviation of observed values, σ_o, and (b) the standard deviation of simulated values, σ_s, is close to the standard deviation of observed values, σ_o. Also, good model performance is indicated if the index of agreement is close to 1 [63]. Furthermore, the model accuracy is considered excellent when RRMSE < 10%, good if 10% < RRMSE < 20%, fair if 20% < RRMSE < 30%, and poor if RRMSE > 30% [65].

Results of the model for bare soil for Groß-Enzersdorf and Goggendorf are shown in Figure 4 and Figure 5, and related statistics are shown in Table 4.

Calibration site statistics over the whole simulation period for Groß-Enzersdorf (Table 4) indicate good model performance: R and IA are close to 1, the RMSE is far below σ_o and σ_s is just 0.14 °C higher than σ_o. RRMSE is well below 10%, indicating excellent model performance.

Validation using a controlled experiment was carried out for a site with arable soil (Goggendorf), with similar conditions to the site of Groß-Enzersdorf. Two stations were placed within a wheat field, representing bare soil and wheat field surface biomass dynamics. The aboveground biomass of the full canopy of standing wheat was estimated from measured yield (set as 7000 kg ha⁻¹ dry matter biomass), and soil porosity and texture at three soil depths were measured once (kept constant during the simulation period), while the daily soil water content and actual evapotranspiration were simulated using the FAO approach [61]. The FAO approach is based on calculating the actual evapotranspiration of crop canopies from grass reference evapotranspiration by using correction factors, which depend on soil water content and crop development. Stubble after wheat harvest was set to 1000 kg ha⁻¹ (straw was removed from field after harvest). Table 4 shows similar performance statistics for Goggendorf for R and IA as at the Groß-Enzersdorf site. For both variants, σ_o and σ_s decreased with depth, following the temperature trend. For both bare soil and vegetated surface simulations in Goggendorf, the RRMSE was slightly higher for bare soil, decreased with depth, and was always well below 10%.

These results, further supported by Figure 5 for the 20 cm soil depth, show that the model can represent soil temperatures within a deviation range of 2 °C on a daily basis at a 20 cm soil depth, considering the existing uncertainties in inputs on soil surface biomass and soil characteristics (simulated soil water content, constant porosity, and soil texture).

As solar radiation interception into canopies and related energy transfer to the soil surface vary not only by canopy structure but also between diffuse and direct solar radiation,

Table 5. Model performance (measured as RMSE, °C) in response to daily solar radiation level (global radiation of 7–15 MJ m⁻² d⁻¹ vs. > 20 MJ m⁻² d⁻¹) for four simulation sites as described in Table 1 (n = number of affected days).

Site Name	Soil Depth		7–15 MJ m−2 d−1		>20 MJ m−2 d−1
	cm	n	RMSE (°C)	n	RMSE (°C)
Doksany, CZ	5	1366	1.36	1207	1.88
Kirklareli, TR	20	724	1.25	560	1.71
Hamme, BE	15	1000	1.85	837	2.15
Grafendorf, AT	10	1717	1.15	1630	1.24

shows a comparison of AGRISOTES performance for days with high versus low solar radiation. Due to the simplified approach used in AGRISOTES (Equation (1)), some differences in model behavior can be expected. Using our database, including four different sites with different canopy conditions (Table 1), we found that on the simulation days with a higher radiation level (7–15 MJ m⁻² d⁻¹ vs. > 20 MJ m⁻² d⁻¹), the model performance is slightly lower (in the range of 0.09 °C to 0.52 °C in RMSE) for both grassland and arable land. However, in order to better calibrate this function, detailed canopy characteristics and light interception data would be necessary.

For the Goggendorf site, we carried out a model sensitivity analysis of the most important parameters for soil temperature simulation dynamics, including surface biomass, soil layer specific soil water content, soil porosity, and soil texture. We defined an artificial simple stepwise change with the most likely realistic ranges for these parameters in agricultural soils, as shown in

Table 6. Basic model sensitivity analysis of the effect of the main important model input parameters on simulated soil temperature at the Goggendorf validation site (see also Table 4).

Reference Simulation 1	Reference Conditions	Soil Depth (cm)	Soil Water Content +10 vol %	Soil Water Content −10 vol %	Porosity +10 vol %	Porosity −10 vol %	SBM 2 +1000 kg ha−1	SBM 2 −1000 kg ha−1	Sand +55% (Ref: 10%)
Simulation period: 1 January 2016 till 31 December 2016	Fallow after winter wheat–15 August harvest – straw mulch – 20 September soil mulching		Δ °C—July 2016 mean deviation from reference
		5	0.06	−0.09	−0.04	0.05	−0.15	−0.19	−0.07
		10	0.05	−0.09	0.00	0.02	−0.15	−0.17	−0.17
		20	0.01	−0.08	0.04	−0.04	−0.14	−0.12	−0.33
		40	−0.02	−0.08	0.11	−0.11	−0.15	−0.07	−0.57
		60	−0.03	−0.09	0.13	−0.14	−0.18	−0.04	−0.75
		100	−0.02	−0.11	0.10	−0.13	−0.23	0.00	−0.90
			Δ °C—October 2016 mean deviation from reference
		5	−0.10	0.07	0.05	−0.07	−0.35	0.58	−0.02
		10	−0.08	0.06	0.01	−0.03	−0.34	0.57	0.10
		20	−0.04	0.04	−0.05	0.03	−0.35	0.53	0.26
		40	0.03	0.00	−0.17	0.16	−0.35	0.47	0.58
		60	0.08	−0.02	−0.27	0.26	−0.33	0.41	0.91
		100	0.15	−0.06	−0.41	0.42	−0.31	0.29	1.45

¹ All sensitivity test input factors modified in respect to the reference. ² SBM: Surface dry biomass change, reference in July = 4000 kg ha⁻¹, in October = 1000 kg ha⁻¹ (straw-soil mulch).

.

The model sensitivity analysis demonstrates that uncertainties in the soil water content may have less impact on model performance than soil porosity and sand content within the often observed vertical and horizontal variations in the soil column of soils under field conditions. The highest impact on soil temperature and its simulation is from the dynamics of soil surface biomass cover under common field conditions (see also Figure 5). The sensitivity to the changing sand content increases with soil depth, whereas the sensitivity to the soil surface biomass affects all depths to the same extent (Table 6).

Finally, for the model validation on a wider range of soil and surface conditions, independent datasets of several years at different sites were used to demonstrate the model performance under various conditions of input data availability and quality. The selected sites include different climates, agricultural soils, and land use, where for model parameterization the local available data were used (such as annual temperature and soil characteristics). The sites are located in Austria, the Czech Republic, Belgium, and Turkey. In all cases, the model input “actual evapotranspiration” was calculated based on the Penman-Monteith equation according to the FAO approach [61]. If measured data were not available, the soil water content input was also calculated by the FAO approach [61]. Snow cover was calculated by the model Snow Mouse [66] when measured data were not available. The statistics of model performance at all validation sites are shown in

Table 7. Model validation results based on available data sets (see Table 1 for site characteristics).

Location and Soil Surface	Not Directly Measured Settings/Inputs	Simulation Period	Available Days	Soil Depth (cm)	R	σo (°C)	σs (°C)	RMSE (°C)	RRMSE (%)	IA (0–1)
Doksany – Permanent short grass	Constant surface biomass assumed	Full years 2001–2004	1454	5	0.98	8.42	8.41	1.65	15.0	0.99
			1454	50	0.99	6.75	6.70	0.97	9.2	0.99
			1454	100	0.99	5.17	5.24	0.87	8.2	0.99
Kirklareli – Winter wheat/fallow	Constant pore volume assumed	Full years 2010–2011	714	5	0.99	8.23	8.69	1.57	10.7	0.99
			714	20	0.99	7.57	8.11	1.45	10.0	0.99
Hamme – Permanent Short grass	Constant soil water content and surface biomass assumed	Full years 2012–2014	1096	15	0.96	5.60	5.64	1.94	15.5	0.97
Pucking – Cultivated soil on lysimeter	Constant pore volume Crop rotation related surface biomass and soil water content calculated	Sept 1996–Feb 1999	745	10	0.97	6.32	6.36	1.59	21.7	0.98
			854	30	0.99	6.86	6.32	1.25	14.2	0.99
Purbach – Cultivated soil under vineyard	Constant pore volume and surface biomass (3000 kg ha−1) Soil water content calculated	July 1999–May 2002	941	10	0.99	7.67	7.72	2.20	23.8	0.98
			941	20	0.99	7.66	7.30	1.19	11.4	0.99
			941	40	0.99	7.43	6.49	1.30	12.1	0.99
			941	70	0.99	7.22	5.45	2.06	19.0	0.97
Grafendorf – 3–cuts Permanent grassland	Surface biomass calculated according to reported cuts–Snow cover simulated	Years 2004–2007	1714	10	0.98	7.11	6.98	1.45	12.9	0.99
			1744	20	0.98	6.88	6.52	1.24	11.1	0.99
			1744	30	0.99	6.59	6.09	1.13	10.2	0.99
			1739	40	0.99	6.37	5.69	1.13	10.3	0.99
			1744	60	0.99	5.26	5.16	0.91	8.3	0.99
Obersieben-brunn – Permanent natural grassland	Constant surface biomass assumed (3000 kg ha−1) Soil water content calculated	Years 2000–2009	2557	10	0.98	6.88	7.66	2.56	28.9	0.97
			2557	50	0.99	5.80	5.23	1.11	10.2	0.99

.

It can be seen that the various uncertainties and indirect estimates of critical soil temperature model inputs are reflected in the performance of results. As expected, most errors can occur when the dynamics of surface biomass are not measured or are difficult to estimate. For example, in the vineyard and permanent grassland sites in Purbach and Obersiebenbrunn, respectively, an increase in the RRMSE, accompanied by a higher deviation between σ_o and σ_s (e.g., 6.88 and 7.66, respectively at 10 cm soil depth in Obersiebenbrunn), is caused by temporal deviation of simulated near soil surface temperatures from measurements (not shown). It seems that known grassland cutting dates and related calculated biomass dynamics contribute to a better performance, such as at sites Doksany and Grafendorf, where the RRMSE values are mostly well below 15. Even at site Pucking, where there was rotation of various crops, and the biomass dynamics were calculated using a simple linear approach based on [61] according to exact dates of growing periods of the different crops, an acceptable result was achieved. Conversely, soil parameter uncertainty and their variation over the vertical soil profile led to deviations, especially in deeper soil layers, such as those at Purbach and Grafendorf. In fact, the physical soil properties for the deeper soil layer of these two sites were extrapolated from the next shallower soil layer. Furthermore, we found that soil conditions beyond the common range for agricultural soils, such as gravel, are not well represented by the current model parameterization (data not shown).

From the weather inputs, the existence of snow cover is the most critical component for a representative simulation of soil temperatures during the winter period. This is especially clear for Grafendorf and Doksany, in Figure 6 and Figure 7, respectively. The site Grafendorf is characterized by a semihumid climate, with regular snow cover during the winter. Figure 6 shows the model performance for the selected soil depth of 10 cm over a four-year period. The overall performance is good, but there are temporal underestimations of simulated soil temperatures during the growing periods, which result from the fact that the grassland surface biomass was calculated according to the reported cutting dates. The largest short-term deviations occur during the winter period and are caused by inadequate (simulated) snow cover duration, which is an input to the soil temperature model. The majority of the daily variation range is below 1 °C.

Figure 6 demonstrates the performance of the simulations over a four-year period for the Grafendorf site, including information on snow cover duration. A strong effect of snow cover on the underlying soil temperatures and a freeze/thaw effect (soil temperatures stay close to zero during ongoing freezing or thawing) can be seen; however, it is difficult to meet the conditions of a not fully frozen or fully frozen state, so that deviations occur more frequently. The impact of snow cover and soil freezing on simulated soil temperature is, as expected, much more pronounced at the upper layer than in deeper soil layers, leading to better model performance in deeper soil layers, e.g., decreasing RMSE and RRMSE (see Table 7 for statistical results). Furthermore, the greatest deviations between measured and simulated soil temperatures are related to sites with periods of snow cover occurrence and/or soil freezing events (see Figure 7a—a site with regular soil frost versus Figure 7c—a site without significant soil frost). Especially at the freezing point (around 0 °C), deviations occur (see also Figure A1, Figure A2, Figure A3, Figure A4, Figure A5, Figure A6, Figure A7, Figure A8, Figure A9, Figure A10, Figure A11, Figure A12, Figure A13, Figure A14, Figure A15, Figure A16, Figure A17 and Figure A18 of nine validation sites in Appendix E).

In general, also the representation of surface dry biomass (SBM) can contribute to deviations and uncertainties (see sensitivity analysis, Table 6), where a deviation of about 1000 kg ha⁻¹ SBM leads to a change of up to 0.6 °C in mean monthly soil temperature. Especially during the crop growing period, this parameter is continuously changing and deviations increase when it is assessed and/or kept constant in the simulation (as at Doksany and Hamme, Figure 7a,b, and other sites, as statistically presented in Table 7, when, e.g., RRMSE values in the top soil layers rise above 15%). The effect of solar radiation levels on model performance is presented statistically in Table 5 and visually for the well-monitored site in Kirklareli, Turkey in Figure 7c,d. It can be seen that, with increasing radiation, the model performance declines, e.g., for levels above 20 MJ m⁻², yielding RMSE values of 1.24–2.15 °C, depending on the site.

Additional validation results for nine independent Czech soil temperature measurement sites for 10 cm and 50 cm soil depth, respectively, are presented in Appendix E (Figure A1, Figure A2, Figure A3, Figure A4, Figure A5, Figure A6, Figure A7, Figure A8, Figure A9, Figure A10, Figure A11, Figure A12, Figure A13, Figure A14, Figure A15, Figure A16, Figure A17 and Figure A18). These sites are characterized by grass cover that is kept at or below 10 cm. At these sites, the soil physical characteristics (especially soil pore volume, wilting point and field capacity) were derived from FAO soil map and soil type description (

Table A1. Description of nine validation sites in the Czech Republic (related to Figure A1, Figure A2, Figure A3, Figure A4, Figure A5, Figure A6, Figure A7, Figure A8, Figure A9, Figure A10, Figure A11, Figure A12, Figure A13, Figure A14, Figure A15, Figure A16, Figure A17 and Figure A18).

Station Name	Station Code	Long	Lat	Altitude	Soil Description
		°	°	m
KROMĚŘÍŽ	B1KROM01	17.3653	49.2847	233	0–100 cm loamy soil, chernozem
LEDNICE	B2LEDN01	16.7989	48.79262	177	0–100 cm loamy soil, chernozem
KOCELOVICE	C1KOCE01	13.83861	49.46722	519	0–25 cm sandy loam, 25–50 cm loam 50–100 cm sandy, Cambisol
UPICE	H1UPIC01	16.01162	50.50642	413	0–70 cm loamy-sand
BĚLOTÍN	O1BELO01	17.8042	49.5869	306	0–100 cm loamy soil, luvisol
LUČINA	O1LUCI01	18.4425	49.7308	300	0–100 cm loam. Pseudogleyic luvisol
LUKÁ	O2LUKA01	16.95333	49.65222	510	0–30 cm sandy-loam 30–100 cm loam, Cambisol
USTÍ NAD LABEM	U1ULKO01	14.04111	50.68333	375	0–30 cm sandy loam 30–100 cm loam, typical Cambisol with strong human influence
LIBEREC	U2LIBC01	15.02389	50.76972	398	0–22 cm sandy-loam 22–100 cm loam, gleyic soil
Surface conditions: All measurement sites are characterized by no or only negligible slope inclination (flat area) and grass cover with regular mowing; AGRISOTES simulations were carried out with constant surface biomass of 1000 kg ha−1 (representing short grass).

, Appendix E). Therefore, larger biases can be expected for these sites due to larger input data uncertainty as other AGRISOTES inputs (e.g., soil water content and evapotranspiration) were simulated by the FAO approach based on these data as well. These additional uncertainties lead to higher RMSE values (overall in the range of 1–2 °C;

Table A2. Statistical validation results of AGRISOTES (measured vs. simulated soil temperatures) for nine validation sites (see Table A1) in the Czech Republic (related to Figure A1, Figure A2, Figure A3, Figure A4, Figure A5, Figure A6, Figure A7, Figure A8, Figure A9, Figure A10, Figure A11, Figure A12, Figure A13, Figure A14, Figure A15, Figure A16, Figure A17 and Figure A18).

Station Name	Station Code	Soil Depth 10 cm			Soil Depth 50 cm
		R2	RMSE	MBE	R2	RMSE	MBE
			°C	°C		°C	°C
KROMĚŘÍŽ	B1KROM01	0.98	1.26	0.52	0.99	1.10	0.72
LEDNICE	B2LEDN01	0.97	1.43	−0.08	0.98	1.16	0.44
KOCELOVICE	C1KOCE01	0.96	1.67	0.45	0.96	1.57	0.65
UPICE	H1UPIC01	0.95	1.83	0.78	0.95	1.83	0.82
BĚLOTÍN	O1BELO01	0.95	1.82	−0.06	0.95	1.65	0.04
LUČINA	O1LUCI01	0.94	1.95	0.20	0.94	1.96	0.36
LUKÁ	O2LUKA01	0.97	1.60	0.70	0.97	1.64	0.72
USTÍ NAD LABEM	U1ULKO01	0.96	1.58	0.37	0.94	1.74	0.31
LIBEREC	U2LIBC01	0.95	1.98	1.08	0.95	1.84	1.14

, Appendix E) compared to our other validation sites (Table 7) where mostly in situ-based soil characteristics were available (Table 1). However, the overall performance of the simulation results are satisfactory over the whole temperature range in view of these uncertainties, and useable for many applications in crop management (e.g., assessing optimum sowing time, soil born pest activities, N-mineralization rates, etc.). Accepted applicability for farmers, however, depends on the specific use cases, which need to be tested before implementation.

4. Conclusions

According to the results of the AGRISOTES model validation performed for different European climatic conditions, land use, soil type, biomass development, snow presence and freezing conditions, it can be concluded that the model:

(a)

is able to produce a smooth transition between bare soil and vegetated surface at the beginning of the vegetation period, as well as after sudden surface condition changes such as harvest;

(b)

represents the influence of surface biomass dynamics on soil temperatures with acceptable accuracy for many agricultural applications;

(c)

RMSE is higher for bare soil than for vegetated surfaces and intensively decreases with depth, indicating that the model is more sensitive to atmospheric forcing (a stronger effect on bare soil) than to the given soil characteristics and the parameterization of the energy transport over the whole profile;

(d)

is highly sensitive to the presence of snow cover; and

(e)

reflects well the most determining soil physical factors such as porosity, soil texture and soil water content. Soil surface conditions, however, remain the greatest influence under agricultural land use.

However, currently, there are still potential limitations on the application of AGRISOTES in different climate zones and under various types of surface conditions. For example, under tropical and other climates and surface conditions such as in wetlands, deserts or forests, AGRISOTES should be evaluated on extended suitable data sets. If the model applicability seems to be given, at least a validation should be carried out first. Based on these validation results, if necessary, a recalibration of the model should be considered, especially of the ν_f(CV)-relationship. In general, we recommend testing with highly accurate measured input datasets for calibration that reflect the greater variability of canopy structure and crop types.