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

: 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 ﬁelds 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), pore the column.


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.

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).

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 −1 ). 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: where h net,g is the net ground heat flow (W m −2 ), α g is the soil surface albedo, h s is the observed solar radiation per unit area of the ground surface (W m −2 ), ε g is the soil emissivity, h l,ac is the incoming atmospheric long-wave radiation (W m −2 ), h l,fg is the net long-wave radiation from the foliage to the ground (W m −2 ), h l,ga is the outgoing long-wave radiation emitted from the bare soil (W m −2 ), h evap,g is the evaporative heat transfer from the soil surface to the atmosphere (W m −2 ) and h conv,g is convective heat transfer from the soil surface to the atmosphere (W m −2 ). 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: where T a is the air temperature ( • C), σ is the Stefan-Boltzmann constant (≈5.67 × 10 −8 W m −2 K −4 ), 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 where T g is the soil surface temperature ( • C) and µ fg is the degree of radiation exchange defined as [53]: where ε f is the foliage emissivity. The outgoing long-wave radiation emitted from the bare soil h l,ga is calculated as 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 −2 ): 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 where k cg is the convective heat transfer coefficient between soil surface and air (W m −2 K −1 ), 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: 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: where the following abbreviations are introduced: This transforms the heat balance equation into the form 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: So, Equation (10) can be rewritten as or in the form where R 1 = 1/k 1 and ∆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 ∆T g,mean = ∆T 1,mean − R 1,mean h net,g,mean , which will be used as the top boundary condition for the simulation of soil temperature during snow-free periods.

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 where k snow is defined as 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): where and For a given SNO value, which is designated as an input to the AGRISOTES model, the snow coverage factor ν sn is calculated as for SNO ≤ SNO limit,1 , SNO−SNO limit,1 SNO limit,2 −SNO limit,1 for SNO limit,1 < SNO < SNO limit,2 , 1 for 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 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:

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: Here, h is the soil heat flow (W m −2 ), T s is the soil temperature ( • C), λ is the thermal conductance of soil (W m −1 K −1 ), C is the volumetric heat capacity of soil (J m −3 K −1 ), F H2O is the latent heat of ice fusion related to the volume of liquid water (≈3.34 × 10 8 J m −3 ), θ ice is the volumetric ice content (cm 3 cm −3 ) 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Ĉ j (z), 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 ∂ z h do not vary significantly during a 24-h period and ∂ z h 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:ĥ In Equations (27) and (28), valuesĈ whereζ is the approximate daily mean ice fraction, defined aŝ with θ denoting the volumetric water content (cm 3 cm −3 ), 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.

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 For the bottom boundary condition, three different modes are available for the soil temperature model: (i) To set the approximated daily mean heat fluxĥ j (z max ) 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 −1 ) and the approximated damping depth at the bottom of the soil profile d j a (z max ) (m) for annual frequency, this leads to the following bottom boundary condition (detailed derivation can be found in Appendix C): .
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 −3 ). 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 T s and ice content θ 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 T s > 0 • C, but during some night hours, freezing may occur so that, nevertheless, θ 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: 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ζ on 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.

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). 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). Meteorological input parameters: Daily mean, maximum and minimum temperature, and solar radiation (global radiation), 2 Soil: Sand and clay content (VP sand , VP clay ) and total soil pore volume (measurement samples from up to four soil depths), 3 θ: Volumetric soil water content (measurements based on 1-2 soil depths), 4 CV: Total above ground dry biomass, 5 Soil texture estimated from soil map (mean value over soil profile), 6 Estimated for calibration of snow cover effect, but measured for validation.

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) as well as the convective heat transfer coefficient (k cg ) for bare soil (Table 3). Taking the form of Equation (2) in [62], a simple exponential curve was fitted for the data in Table 2. The resulting relationship was 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 −1 , 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 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 Figures 4 and 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 −1 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 −1 (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 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 −2 d −1 vs. > 20 MJ m −2 d −1 ), 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.
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.
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). 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  1 All sensitivity test input factors modified in respect to the reference. 2 SBM: Surface dry biomass change, reference in July = 4000 kg ha −1 , in October = 1000 kg ha −1 (straw-soil mulch). Table 7. Model validation results based on available data sets (see Table 1 for site characteristics).   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 Figures 6 and 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.     -d). Comparison of measured and simulated daily mean soil temperatures at three climatologically different sites in Belgium, the Czech Republic, and Turkey over the full measured temperature ranges (statistics shown in Table 7). 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 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 −1 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 wellmonitored 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 −2 , yielding RMSE values of 1.24-2.15 • C, depending on the site.

Location and Soil Surface Not Directly Measured Settings/Inputs
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 ( Figures A1-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, 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, 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.

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.  where ∆T g,mean is defined as the difference between the daily mean soil surface temperature T g,mean and the daily mean air temperature T a,mean ∆T g,mean = T g,mean − T a,mean . (A7) The value of k 1 , as defined in Equation (11), changes during the day because of k l0 air temperature dependence. The daily mean value of k 1 is calculated as follows: using the approximation k l0,mean = 4σ(T a + 273.15) 3 ≈ 4σ(T a,mean + 273.15) 3 (A9) and k l,g0,mean = ε g k l0,mean , where T a,mean is the daily mean air temperature ( • C). Now, Equation (A3) can be written, approximately, in the form h net,g,mean = h t,1,mean − k 1,mean ∆T g,mean , or in the form of Equation (15), which is more convenient for further calculations, where ∆T 1,mean = h t,1,mean k 1,mean . (A13)

Appendix C. Derivation of Bottom Boundary Equations
Combining Equations (25) and (26) (neglecting θ ice change, because the bottom layer of the simulated soil profile is assumed to be deep enough that θ ice does not change there) and assuming constant thermal conductance λ and heat capacity C yields The solution to Equation (A14) for sinusoidal variation with annual angular frequency ω a and amplitude T ampl around a mean value of T AA is known as [70]: where the angular frequency of a year ω a (rad s −1 ) is defined as ω a = 2π 365 · 86400 (A16) and the damping depth d a for annual frequency is given as [70]: Using Equations (25) and (A15) and the trigonometric addition theorems it can be verified that the following equation for the heat flow is true: Time-averaging yields (∂ z T s is assumed to be continuous) From Equation (A19) now follows the bottom boundary condition Equation (33), where the approximated damping depth d j a (z) (m) for annual frequency has been defined as follows [70]: (A20) vol % volume percentage of total solid soil matter that is clay VP sand vol % volume percentage of total solid soil matter that is sand w J m −3 heat energy density z m depth in soil profile (z = 0 at soil surface, z > 0 inside the soil) z max m maximum simulation depth (bottom of soil profile) Part II α g 1 soil surface albedo (0-1) β ha kg −1 parameter in ν f (CV)-relation ∆T 1 • C difference between soil surface and air temperature during snow-free periods for zero ground heat flow ∆T 1,mean • C difference between daily mean soil surface and air temperature during snow-free periods for zero ground heat flow ∆T 2,mean • C difference between daily mean soil surface and air temperature for dense snow cover and zero ground heat flow ∆T g • C difference between soil surface and air temperature ∆T g,mean • C difference between daily mean soil surface and air temperature

Appendix D. Notations and Units
∆T mean • C difference between daily mean soil surface and air temperature for partial snow cover and zero ground heat flow ∆t s simulation time step (= 86,400 s) ε ac 1 atmospheric emissivity (0-1) ε f 1 foliage emissivity (0-1) ε g 1 ground (soil surface) emissivity (0-1) ζ 1 ice fraction (0-1) θ cm 3 cm −3 volumetric total soil water content (liquid + ice, but neglecting volume expansion of frozen water) θ ice cm 3 cm −3 volumetric ice content in soil, but neglecting volume expansion of frozen water θ sat cm 3 cm −3 saturated volumetric water content or total soil pore volume λ W m −1 K −1 soil thermal conductivity µ fg 1 radiation exchange degree between foliage and ground ν f 1 vegetation cover (0-1) ν sn 1 snow coverage factor (0-1) σ W m −2 K −4 Stefan-Boltzmann constant (≈5.67 · 10 −8 W m −2 K −4 ) ω a rad s −1 angular frequency of the year Appendix E. 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).  Table A1) of Czech Republic (graphical representation).