A Coupled Model for Simulating Water and Heat Transfer in Soil-Plant-Atmosphere Continuum with Crop Growth

Crop growth is influenced by the energy partition and water–heat transfer in the soil and canopy, while crop growth affects the land surface energy distribution and soil water-heat dynamics. In order to simulate the above processes and their interactions, a new model, named CropSPAC, was developed considering both the growth of winter wheat and the water–heat transfer in Soil-Plant-Atmosphere Continuum (SPAC). In CropSPAC, the crop module depicts the dynamic changes of leaf area index (LAI), crop height, and the root distribution and outputs them to the SPAC module, while the latter outputs soil moisture conditions for the crop module. CropSPAC was calibrated and validated by field experiment of winter wheat in Yongledian, Beijing, with five levels of irrigation treatments, namely W0 (0 mm), W1 (60 mm), W2 (110 mm), W3 (170 mm), and W4 (230 mm). Results show that CropSPAC could predict the soil water and temperature distribution, and winter wheat growth with acceptable accuracy. For example, for the 0–1 m soil water storage, the R2 for W0, W1, W2, W3, and W4 is 0.90, 0.88, 0.90, 0.91, and 0.79, and the root mean square error (RMSE) is 17.24 mm, 27.65 mm, 20.47 mm, 22.35 mm, and 12.88 mm, respectively. For soil temperature along the soil profile, the R2 ranges between 0.96 and 0.98, and the RMSE between 1.22 ◦C and 1.94 ◦C. For LAI, the R2 varied from 0.76 to 0.96, and the RMSE from 0.52 to 0.67. We further compared the simulation results by CropSPAC and its two detached modules, i.e., crop and the SPAC modules. Results demonstrate that the coupled model could better reflect the interactions between crop growth and soil moisture condition, more suitable to be used under deficit irrigation conditions.


Introduction
Crop habitat, including the water and heat conditions in the soil root zone and crop canopy, is crucial for crop growth and yield formation.Plants extract water from soil, transport water within xylem vessels, and lose water by transpiration through stomata in plant leaves [1].Transpiration flux provides the driving force to capture nutrients from the soil.By absorbing the energy in solar radiation, crop photosynthesis could be accomplished which produce carbohydrate for biomass accumulation and distribution [2].Additionally, the surrounding temperature is a vital signal for determining the crop phenological period.On the other hand, crop growth conditions also affect crop habitats by influencing the land surface energy distribution, partition of solar energy, and root zone development and water uptake, etc. [3].Understanding their interactions is essential for scientific management of As shown in Figure 1, the output of the crop module provides input data for the SPAC module, mainly the crop information, such as LAI, root distribution, plant height, etc.The output of the SPAC module provides the input data for the crop module, mainly for the soil moisture condition, such as the averaged soil moisture content in the root zone.The calculation flowchart is illustrated in Figure 2. As shown in Figure 1, the output of the crop module provides input data for the SPAC module, mainly the crop information, such as LAI, root distribution, plant height, etc.The output of the SPAC module provides the input data for the crop module, mainly for the soil moisture condition, such as the averaged soil moisture content in the root zone.The calculation flowchart is illustrated in Figure 2.
Depending on the requirement of stability, explicit or implicit methods could be applied for the coupling between Crop module and SPAC module.In the implicit method, first, the crop data from previous/initial time step is put into the SPAC water-heat transfer module to obtain the corresponding predicted soil water/temperature status at the end of the time step.Then, the predicted soil water/temperature data was used in the crop module to derive the predicted crop data at the end of the time step.The predicted crop data is then used as the input for the other round of simulation till the error of predicted data in the latest two iterations is below a certain threshold.Then the simulation of a particular time step is over and the next time step starts till the process reaches the end simulation period, generally the end of the crop season.The crop module and SPAC module are described separately in the following sections.Depending on the requirement of stability, explicit or implicit methods could be applied for the coupling between Crop module and SPAC module.In the implicit method, first, the crop data from previous/initial time step is put into the SPAC water-heat transfer module to obtain the corresponding predicted soil water/temperature status at the end of the time step.Then, the predicted soil water/temperature data was used in the crop module to derive the predicted crop data at the end of the time step.The predicted crop data is then used as the input for the other round of simulation till the error of predicted data in the latest two iterations is below a certain threshold.Then the simulation of a particular time step is over and the next time step starts till the process reaches the end simulation period, generally the end of the crop season.The crop module and SPAC module are described separately in the following sections.The GAI is the green area index; E vp is the potential transpiration; E s is the soil evaporation rate; G is the downward heat flux in the soil surface; u is the wind speed; T is the air temperature; RH is the relative humidity; and n are the sunshine hours.

Crop Module
The crop module consists of three parts, i.e., simulation of crop stage development, simulation of biomass accumulation and simulation of dry biomass distribution and final yield.Here, the crop considered is winter wheat, which is a common cereal planted in North China.The calculation flowchart of the crop module is illustrated in Figure 3.

Simulation of Crop Stage Development
In this paper, the description of crop stage development is based on Cao et al. [27,28].The physiological processes of stage development are generally regulated and influenced by the temperature and photoperiod.The stage development response to temperature is manifested in two forms.One is the thermal effect (TE), which the rate of stage development increases with the increase of thermal effect.The other is vernalization effect (VE), which only functioned in cold season crops such as winter wheat, generally referring to the phenomenon that plants must undergo a period of sustained low temperature then it can move from the vegetative growth stage to the reproductive stage.The sensitivity of stage developmental rate to photoperiod or daytime are called photoperiod effect (PE).For cold season crops, it is generally characterized by the longer sunshine time which will promote the crop growth rate.In addition to the temperature and light response characteristics, the crop stage development has also received species and genetic influences.
Physiological development time (PDT) is a time scale of crops in the most suitable development environment.Under any light and temperature conditions, the time for a specific crop to complete the entire growth period is essentially constant, i.e., the PDT value is basically constant, which is called the principle of constant physiological development time.Therefore, the development stage of specific genotypes crop can be predicted by [27]: where BDF is physiological development factor, a specific genetic parameter related to the variety of the crop, the range of wheat is 0.6-1.0 according to Cao et al. [27]; the PDTt is the physiological development time without consideration of crop variety impact, which was accumulated by the daily physiological effect (DPE): The daily physiological effect (DPE) can be represented by daily thermal effects (DTE) and daily thermal sensitive (DTS) interactions:

Simulation of Crop Stage Development
In this paper, the description of crop stage development is based on Cao et al. [27,28].The physiological processes of stage development are generally regulated and influenced by the temperature and photoperiod.The stage development response to temperature is manifested in two forms.One is the thermal effect (TE), which the rate of stage development increases with the increase of thermal effect.The other is vernalization effect (VE), which only functioned in cold season crops such as winter wheat, generally referring to the phenomenon that plants must undergo a period of sustained low temperature then it can move from the vegetative growth stage to the reproductive stage.The sensitivity of stage developmental rate to photoperiod or daytime are called photoperiod effect (PE).For cold season crops, it is generally characterized by the longer sunshine time which will promote the crop growth rate.In addition to the temperature and light response characteristics, the crop stage development has also received species and genetic influences.
Physiological development time (PDT) is a time scale of crops in the most suitable development environment.Under any light and temperature conditions, the time for a specific crop to complete the entire growth period is essentially constant, i.e., the PDT value is basically constant, which is called the principle of constant physiological development time.Therefore, the development stage of specific genotypes crop can be predicted by [27]: where BDF is physiological development factor, a specific genetic parameter related to the variety of the crop, the range of wheat is 0.6-1.0 according to Cao et al. [27]; the PDT t is the physiological development time without consideration of crop variety impact, which was accumulated by the daily physiological effect (DPE): Water 2019, 11, 47 6 of 32 The daily physiological effect (DPE) can be represented by daily thermal effects (DTE) and daily thermal sensitive (DTS) interactions: DTS can be calculated from [28]: where the VP is vernalization process; PDTTS and PDTHD are PDT requirements for the terminal spikelet stage and heading stages, 10 and 27, respectively.The detail calculation methods of the DTE, VP, and PE can be found in Appendix A, and the values of some parameters above (e.g., BDF) are shown in Appendix C, Table A2.

Simulation of Biomass Accumulation
The formation of biomass is mainly related to the physiological and ecological processes, such as photosynthesis and respiration.The photosynthesis and biomass accumulation is based on the calculation of photosynthetic rate of single leaf, then the daily assimilation of canopy are calculated, subtracting the consumption of respiration, finally the daily net assimilation amount is calculated. (

1) Photosynthesis of crop canopy
The photosynthetic rate of single leaf is calculated as [29]: where FG is the photosynthetic rate of leaves, kgCO 2 •ha −1 •h −1 ; PLMX 0 is the maximum photosynthetic rate of leaves, kgCO 2 •ha −1 •h −1 ; ε is the initial utilization efficiency of absorbed light, kgCO 2 •ha −1 •h −1 •(J•m −2 •s −1 ) −1 ; I L is the photosynthetically active radiation intensity at L (depth of crop canopy), J•m −2 •s −1 ; FH is the water influence factor; FN is the nitrogen influence factor (note that the current CropSPAC model could not simulate the soil nitrogen migration, however, any nitrogen stress could be empirically indicated by this parameter); FT is the temperature influence factor; FC is the CO 2 concentration influence factor.Calculation methods and details are shown in Appendix A.
Based on the single leaf photosynthesis model, the Gauss integral method [30] is used to calculate the daily photosynthesis rate of the canopy, which can greatly reduce the amount of calculation under the premise of ensuring the accuracy of calculation.The Gauss integral method divides the leaf canopy into five layers, calculates the instantaneous assimilation rate of each layer and obtains the instantaneous assimilation rate of the whole canopy.By selecting the three time points from noon to sunset, the canopy assimilation rate at three time points were calculated by weighted summation, and the daily canopy assimilation rate (FDTGA) was obtained.
(2) Biomass formation The net assimilation rate of the population is calculated by: where PND is net daily assimilation of population, kgCO  For daily increment of dry biomass in population: where TDRW is daily increment of dry biomass in population, kgDM•hm −2 •day −1 ; ξ is conversion coefficients between CO 2 and organic compounds.
For accumulation of dry biomass: where W DAY is accumulation of dry biomass, kgDM•hm −2 •day −1 .

Dry Biomass Partitioning and Yield
Based on the dynamic relationship between the distribution index and the growth process [31], the distribution index of wheat organs during the growing period was established as: P STEM = 1 − P LEAF − P EAR (13) where P TOP is the aboveground part distribution index, P ROOT is the underground part distribution index, P LEAF is the leaf distribution index, P EAR is the ear distribution index, P STEM is the stem and sheath distribution index, and HI is the harvest index.
In the distribution of dry biomass, the distribution index of aboveground and underground parts of the soil is also affected by the water status, and the distribution index of the leaf dry biomass is affected by both water and nitrogen.
The approximate linear relationship between winter wheat yield (YIELD) and spike weight (W EAR ) was established by Cong [29]: Based on previous study [32], the empirical relationship between LAI and leaf dry biomass quality (W LEAF ) was given as: where SLA is leaf weight per unit, general range 0.0022~0.0045ha•kg −1 [32], and in this study it is calibrated by field data (value shown in Appendix A).

SPAC Module
In this study, the SPAC module was built on the previous work by Mao et al. [33], who established the water-heat transfer model in SPAC under soil water stress, with one canopy layer, based on soil hydrodynamics, micro-meteorology, plant physiology theory and energy balance principle.The calculation flowchart of the SPAC module is shown in Figure 4.

Calculation of Net Radiation and Water-Heat Transport in Crop Canopy
The formula for the net radiation of the lower cushion surface can be expressed as: ( ) where Rn is the net radiation, W•m −2 ; Rg is total short-wave radiation, W•m −2 ; A is the land surface albedo; F is net long-wave radiation, W•m −2 .The Rg, and F can be calculated by theoretical or empirical methods shown in the Appendix B [34], the value of A can be seen in Appendix C, Table A2.Net radiation of underlying surface Rn can be divided into two parts: where Rv is the net radiation absorbed by crop canopy, W•m −2 ; Rs is the net radiation absorbed by ground surface, W•m −2 .
The Rs is calculated from the empirical formula summarized from field observations [35]: where A and B is the empirical coefficient: for winter wheat, A = 0.10364, B = 0.3937; t is the local time; and SN is the time of solar noon.
Neglecting the energy consumed by photosynthesis and the amount of moisture stored in the canopy air, the canopy energy distribution can be described by [36]:

Calculation of Net Radiation and Water-Heat Transport in Crop Canopy
The formula for the net radiation of the lower cushion surface can be expressed as: where R n is the net radiation, W•m −2 ; R g is total short-wave radiation, W•m −2 ; A is the land surface albedo; F is net long-wave radiation, W•m −2 .The R g , and F can be calculated by theoretical or empirical methods shown in the Appendix B [34], the value of A can be seen in Appendix C, Table A2.Net radiation of underlying surface R n can be divided into two parts: where R v is the net radiation absorbed by crop canopy, W•m −2 ; R s is the net radiation absorbed by ground surface, W•m −2 .The R s is calculated from the empirical formula summarized from field observations [35]: where A and B is the empirical coefficient: for winter wheat, A = 0.10364, B = 0.3937; t is the local time; and SN is the time of solar noon.
Water 2019, 11, 47 Neglecting the energy consumed by photosynthesis and the amount of moisture stored in the canopy air, the canopy energy distribution can be described by [36]: where

Calculation of Crop Transpiration under Water Stress in Root Zone
When the root zone is in low soil water condition, the water absorption rate by the crop root would decrease, resulting in the reduction of stomatal conductance, the decrease of leaf water potential, and the increase of stomatal resistance (r v ), which induce the reduction of canopy transpiration.Plant physiologists have studied the effects of soil moisture changes on leaf water potential and stomatal resistance [37], but due to the complex mechanism involved, it is difficult to quantify directly the solve leaf water potential and stomatal resistance under soil water stress.Neglecting the water storage variation in the crop, we assume the root water uptake equals the canopy transpiration.Therefore, the response of canopy transpiration to soil water stress can be reflected by the actual root water uptake described by the soil water stress coefficient of root water uptake k(ψ).The actual transpiration of crops under the condition of water stress, E v , can be calculated by: Water 2019, 11, 47 10 of 32 where k(ψ) is the soil water stress coefficient of root layer, depends on the distribution of root density and the corresponding function of water stress, which is derived in next section and shown in Equation ( 34); E vp is the potential transpiration without water stress in root zone.
Hence the latent heat flux between crop leaves and canopy air (LE v , see Equation (A42)) can be replaced by Equation ( 26): where r v0 is the canopy total stomatal resistance under no soil water stress, s•m −1 .

Water and Heat Transfer in Soil with Root Water Uptake
Richard's equation under root water absorption conditions and heat conduction equation is used to describe soil water-heat transfer: where θ is the soil water content, cm 3 •cm −3 ; T is the soil temperature, • C; t is the time, s; z is the depth from the surface, m; D w is the coefficient of soil water diffusion, m 2 •s −1 ; D wh are the diffusion coefficients of temperature gradients to water flow, There are three types of root density distribution function commonly used in the research, including the evenly distribution, the linear variation distribution with root depth, and the exponential distribution with root depth.Here we used the linear variation function, assuming the ratio of upper part root density to lower part root density is 0.7/0.3,which is in accordance with previous research [38].
We assume the maximum potential transpiration E vp , which is not limited by soil water stress, equal to the potential total root water uptake rate, i.e., the integration of potential root water uptake rate along the entire root zone considering the root density distribution.The actual root water uptake rate s(z, t) is equal to the potential root water uptake rate s 0 (z, t) multiplied by water stress response function α (ψ, π) [39]: where ξ'(z, t) is relative root density function; ξ(z, t) is root density function; ψ(z), π(z) is soil matric potential and solute potential at the depth z below ground surface respectively, P is an empirical constant, ψ 50 is the specific soil water potential, when E vp decreases by 50%(ψ 50 ≈ −0.43 MPa); Z r is the depth of the root zone.Through Equations ( 32)-( 34), the expressions of the root water uptake rate distribution, the crop transpiration and the corresponding soil water stress coefficient of root water uptake at a certain time can be derived as: For the two non-linear partial differential equations (see Equation ( 27)), they are high non-linear, could be solved by implicit finite difference methods in order to reduce the numerical oscillations.Iterations between the above ground and soil processes might be required considering the calculation stability.

Study Area and Experiment
The experiment data were collected from a farmland planted with winter wheat (Triticum aestivum L.) in Yongledian Experiment Station, Beijing, China from October 1998 to June 1999.The experiment station is located at 116 • 40 E and 39 • 47 N with an average elevation of 12 m above sea-level.The average annual temperature in the region is about 11.5 • C. The average annual rainfall is about 570 mm, and most of the rain falls in the period from June to August.The soil is sandy loam with the bulk density of 1.4 g•cm −3 , and the depth of water table is mostly larger than 5 m in this area.
Daily meteorological data were recorded by a weather station at the experimental site, including solar radiation, sunshine hours, air temperature, humidity, atmospheric pressure, wind speed, and precipitation.Soil water content was measured with neutron hydroprobe and tensiometer (503DR, Campbell Pacific Nuclear Corp., Martinez, San Francisco, USA), which were measured every 20 cm soil layer till the depth of 100 cm.For the soil temperature measurement, thermistors manufactured by the Chinese Academy of Agriculture Sciences were installed and were monitored in soil depth of 10 cm, 20 cm, 30 cm, 40 cm, 60 cm, 80 cm, and 100 cm.Crop characteristics, e.g., crop height, was measured with a meter ruler, and leaf area index was calculated as the product of plant density and average green leaf area of individual plant, with the leaf area estimated from the product of length, maximum width and leaf form factor (0.83 for winter wheat) [40].Dry biomass was calculated from the product of plant density and average dry weight of individual plant, with the dry weight determined after drying for 24 h at 80 • C in a fan-forced convection oven.
The experiment field was divided into 24 experiment plots, with the area of each plot 50 m 2 (5 m × 10 m), as shown in Figure 6.It includes six rows for five irrigation treatments (two rows for W2 treatment), and four columns for three fertilization treatments (two columns for medium level fertilization treatment).Since only column FM-1 has crop growth index monitoring, the medium fertilization with five levels of water treatments were used for model calibration and validation.In those treatments, the level of fertilization is the same, which was reflected in the factor FN (e.g., FN is assumed to be 0.95) in CropSPAC.The details of fertilization volume, irrigation volume, and time were shown in Table 1.
Usually, the winter wheat grows in dry seasons from October to early June of the next year.During the soil frozen period from late November to next February, wheat grows very slowly or even stops for about 100 days.Note that our current model for soil water-heat transfer does not include the soil water freezing and thawing period.Under such condition, more complex mechanism should be involved and could only be solved by some specific soil water-heat transfer model, like THAW [41].Since our model has difficulty to simulate the soil frozen period, we set the CropSPAC model to start simulation on March 20, and assume the initial value of the crop information at this time based on empirical experience (the initial value of the crop information can be found in Table 2 in Section 3.2).
× 10 m), as shown in Figure 6.It includes six rows for five irrigation treatments (two rows for W2 treatment), and four columns for three fertilization treatments (two columns for medium level fertilization treatment).Since only column FM-1 has crop growth index monitoring, the medium fertilization with five levels of water treatments were used for model calibration and validation.In those treatments, the level of fertilization is the same, which was reflected in the factor FN (e.g., FN is assumed to be 0.95) in CropSPAC.The details of fertilization volume, irrigation volume, and time were shown in Table 1.

Model Input
The input for the simulation includes meteorological conditions, crop parameters, initial crop physiological index, soil water and thermal characteristics parameters, and initial and boundary conditions of soil temperature and moisture.
Meteorological parameters of the model input include daily average wind speed, sunshine hours, daily maximum and minimum temperature, daily maximum and minimum water pressure, irrigation and rainfall, etc., which were provided by weather station.The crop parameters related with winter wheat were illustrated in the Appendix C, Table A2, as mentioned in Section 2.1.The simulation period starts from March 20.The initial crop physiological index on this date was determined based on the empirical experience with calibration, as shown in Table 2.
The relationships among soil water content, matric potential and hydraulic conductivity can be described by van Genuchten-Mualem (VGM) model [42,43], and the related soil hydraulic parameters for the station were calibrated with field data and shown in Table 3. Soil thermal parameters, K h and C v , mainly vary with the soil water content, can be described by Johansen model [44].For the soil water and heat movement, the equations were solved numerically by finite difference method.The time step usually was in hourly level which could be automatically adjusted according to the numerical stability.Note that the weather information and the other boundary conditions were usually obtained on daily bases, which needed to be refined accordingly to fit the requirement of this numerical simulation, e.g., for the hourly temperature, it is obtained based on the sine function and the daily monitored maximum and minimum temperature.For the more detailed temperature, it is obtained by linear interpolation based on the hourly temperatures.At the beginning of the simulation, the initial soil moisture and temperature along the soil profile were given by the measured values.
When there was no precipitation and irrigation, the upper boundary of soil moisture movement was the evaporation boundary, and the evaporation rate Es was calculated based on the above ground water-heat transfer and energy participation.When there was irrigation or precipitation, the upper boundary condition was switched to the Dirichlet boundary.The upper boundary of soil temperature was Neumann boundary, with the flux, i.e., the downward heat flux to the soil, obtained from the above ground calculations.The lower boundaries of soil water and soil temperature were set to be the Dirichlet boundary, based on the measured values.

Simulation Results
The established CropSPAC model was manually calibrated and validated by field data in Yongledian.W0, W2, and W4 treatments were used for calibration, then the W1 and W3 treatments were used for validation to the calibrated model.

Soil Water Content
The calibration and validation results of soil water storage in the 0-1 m soil layer are shown in Figures 7 and 8.All the treatments show that the soil water storage in the root zone was affected strongly by the surface irrigation/precipitation and the evapotranspiration process.After each event of irrigation or large rainfall, the soil moisture increased sharply (such as, on March 26, April 21, May 4, and May 9 in W4), and then decreased gradually due to the spring wheat roots water uptake, soil surface evaporation and deep percolation, similar to the previous studies, e.g., Shang et al. [45].In the early stage of crop growth, the soil water storage decreased more slowly because of the less root water uptake and small evapotranspiration in the early stage of crop growth.In the middle and late growth stages, with the increase of LAI and the increase of potential atmospheric evaporation, the evapotranspiration of crops increased, resulting in the rapid decline of soil water storage.The R 2 in W0, W2, W4, W1, and W3 is 0.90, 0.90, 0.79, 0.88, and 0.91, and is RMSE 17.24 mm, 20.47 mm, 12.88 mm, 27.65 mm, and 22.35 mm, respectively.

Soil Temperature
The simulated and measured soil temperature at different soil layers were shown in Figure 11.(note that only W0 and W4 have measured temperature data).The net radiation absorbed by the ground surface, Rs, was mainly consumed by the soil latent heat, sensible heat, and the downward heat flux.In this simulation period, with increase in atmospheric temperature and surface downward heat flux (G), the soil temperature in each layer showed a gradually increasing trend.Additionally, soil temperature fluctuated more greatly in upper soil layers, because the upper soil was more exposed to atmospheric temperature and surface radiation.
The simulated results showed good agreement with the measured values.The range of R 2 varied between 0.96 and 0.98, and RMSE ranges from 1.22 °C to 1.94 °C.However, it seems the simulated results systematically underestimated soil temperature generally in the early crop stage.The reason may be, (1) according to the observed air temperature, there is a large drop on March 22.It is captured in the simulation with the soil temperature dropped accordingly, especially in the near surface layers.However, strangely the observed data did not show this drop, possibly because some artificial measure had been taken to keep the soil temperature which was unable to be captured in our model; and (2) due to the uncertainty of local crop parameters, the simulated LAI tended to be higher in the early crop stage.Thus, the larger canopy meant larger interception of radiation and less radiation absorbed by the soil surface.Therefore, the simulated soil temperature tended to be small in the observed data in the early crop stage.With more work on the accuracy of crop data and the detailed investigation of agriculture management, the simulation performance was able to be improved in future research.

Leaf Area Index
Figures 12 and 13 present the comparison between simulated and measured LAI results, which showed the dynamic changes of LAI of winter wheat from turning green to maturity stage.Simulated and measured values showed similar trends, with the R 2 ranging between 0.76 to 0.96, RMSE between 0.52 and 0.67.From early March when winter wheat turned green, leaf area and dry biomass gradually increased with the increase of air temperature.When winter wheat entered the grain filling stage, the LAI reached the maximum value.Then the leaves got senescent and LAI declined, and the Although the trend between observation and simulation are similar, there was still some error between the two sets of data, especially in the late stage where the simulated results showed higher soil water storage than the observed data, e.g., Figure 7b,c, and in some soil profiles simulation results, e.g., Figures 9d and 10b.One of the reason for simulated deviation should be the lack of local crop parameter information which requires further study to determine them more accurately.The other should be caused by the soil heterogeneity.In general, the simulated soil water content values were in reasonable agreement with the measured value under the different water treatments.

Soil Temperature
The simulated and measured soil temperature at different soil layers were shown in Figure 11.(note that only W0 and W4 have measured temperature data).The net radiation absorbed by the ground surface, Rs, was mainly consumed by the soil latent heat, sensible heat, and the downward heat flux.In this simulation period, with increase in atmospheric temperature and surface downward heat flux (G), the soil temperature in each layer showed a gradually increasing trend.Additionally, soil temperature fluctuated more greatly in upper soil layers, because the upper soil was more exposed to atmospheric temperature and surface radiation.
The simulated results showed good agreement with the measured values.The range of R 2 varied between 0.96 and 0.98, and RMSE ranges from 1.22 • C to 1.94 • C.However, it seems the simulated results systematically underestimated soil temperature generally in the early crop stage.The reason may be, (1) according to the observed air temperature, there is a large drop on March 22.It is captured in the simulation with the soil temperature dropped accordingly, especially in the near surface layers.However, strangely the observed data did not show this drop, possibly because some artificial measure had been taken to keep the soil temperature which was unable to be captured in our model; and (2) due to the uncertainty of local crop parameters, the simulated LAI tended to be higher in the early crop stage.Thus, the larger canopy meant larger interception of radiation and less radiation absorbed by the soil surface.Therefore, the simulated soil temperature tended to be small in the observed data in the early crop stage.With more work on the accuracy of crop data and the detailed investigation of agriculture management, the simulation performance was able to be improved in future research.However, the model overestimated the value of LAI in early stage of crop growth.One possible reason is that the LAI in the model was calculated by the relationship between the LAI and the dry

Leaf Area Index
Figures 12 and 13 present the comparison between simulated and measured LAI results, which showed the dynamic changes of LAI of winter wheat from turning green to maturity stage.Simulated and measured values showed similar trends, with the R 2 ranging between 0.76 to 0.96, RMSE between 0.52 and 0.67.From early March when winter wheat turned green, leaf area and dry biomass gradually Water 2019, 11, 47 18 of 32 increased with the increase of air temperature.When winter wheat entered the grain filling stage, the LAI reached the maximum value.Then the leaves got senescent and LAI declined, and the dry biomass accumulated by the plant transfers from the organs to the ear.The peak and time when the LAI reached this maximum value in the three treatments were different as shown in Figure 11.For example, in W4 treatment with the adequate irrigation water supply, it is in the beginning of May that LAI reached its peak value of 5.6, while in W2 and W0 treatments the time is in early April, with the peak values also smaller at between 4 and 5.This illustrated that deficit irrigation limited the growth of crop and the extension of the foliage, and led to the crop early mature.This is also a kind of physiological protection measures for plants in the shortage of water supply, mature in advance to avoid large-scale production loss, as demonstrated by other research [46,47].
However, the model overestimated the value of LAI in early stage of crop growth.One possible reason is that the LAI in the model was calculated by the relationship between the LAI and the dry biomass of leaf, and the simulation of dry biomass by the model was slightly larger in early stage (seen in Figures 13 and 14), so that the simulated value of the LAI also appeared to be large in the previous period.The coefficient between LAI and leaf dry biomass was empirical and may differ from this specific local crop, thus contributing to the inaccurate simulation value.The gap between measured and simulated value can also be induced by the measurement errors.In general, the gradient of LAI between different water treatments is obvious, and the model fitting results were basically acceptable.
biomass of leaf, and the simulation of dry biomass by the model was slightly larger in early stage (seen in Figures 13 and 14), so that the simulated value of the LAI also appeared to be large in the previous period.The coefficient between LAI and leaf dry biomass was empirical and may differ from this specific local crop, thus contributing to the inaccurate simulation value.The gap between measured and simulated value can also be induced by the measurement errors.In general, the gradient of LAI between different water treatments is obvious, and the model fitting results were basically acceptable.4 shows the measured and simulated yields.With the growth of spring wheat, the dry aboveground biomass gradually increased.The CropSPAC model simulated the accumulation of D-AGBin spring wheat with the R 2 varied from 0.93 to 0.99 and RMSE from 993 kg/ha to 1729 kg/ha in different treatments.The increase in irrigation contributed to the accumulation of dry biomass.The reduction in water directly led to a reduction in dry biomass accumulation and then affected the final yield.As a whole, the simulation of the CropSPAC on the accumulation of winter wheat biomass and the final yield is reliable.4 shows the measured and simulated yields.With the growth of spring wheat, the dry aboveground biomass gradually increased.The CropSPAC model simulated the accumulation of D-AGBin spring wheat with the R 2 varied from 0.93 to 0.99 and RMSE from 993 kg/ha to 1729 kg/ha in different treatments.The increase in irrigation contributed to the accumulation of dry biomass.The reduction in water directly led to a reduction in dry biomass accumulation and then affected the final yield.As a whole, the simulation of the CropSPAC on the accumulation of winter wheat biomass and the final yield is reliable.

Comparison between CropSPAC and the Detached Crop Module
To illustrate the advantage of CropSPAC compared with the crop growth model without considering the influence of SPAC water-heat transfer, the crop module was detached from CropSPAC for simulation of crop growth.Since the soil moisture was regarded as specified input, the detached crop module may not transfer soil water deficit information to crop growth under deficit irrigation, therefore the crop growth indexes such as LAI and dry aboveground biomass could not dynamically respond to soil moisture changes, resulting in higher simulated values as shown in Figure 16 for simulating W0 treatment.While in CropSPAC the crop growth status could be directly

Comparison between CropSPAC and the Detached Crop Module
To illustrate the advantage of CropSPAC compared with the crop growth model without considering the influence of SPAC water-heat transfer, the crop module was detached from CropSPAC for simulation of crop growth.Since the soil moisture was regarded as specified input, the detached crop module may not transfer soil water deficit information to crop growth under deficit irrigation, therefore the crop growth indexes such as LAI and dry aboveground biomass could not dynamically respond to soil moisture changes, resulting in higher simulated values as shown in Figure 16 for simulating W0 treatment.While in CropSPAC the crop growth status could be directly influenced by the updated soil moisture status, making the coupled model, CropSPAC, more accurate for simulating crop growth under water deficit conditions.
influenced by the updated soil moisture status, making the coupled model, CropSPAC, more accurate for simulating crop growth under water deficit conditions.

Comparison between CropSPAC and SPAC
CropSPAC model could also outcompete the detached SPAC module for water-heat transfer that ignores the influence of crop growth.SPAC model only accounts for the water-heat transfer from soils to plants, and then to atmosphere, while the crop growth and canopy development is not taken into account.Therefore the detached module of SPAC is inadequate because there is a close relationship between crop growth and water-heat transfer, where the leaf extension and root development of crops will actively affect energy distribution in the canopy and root water absorption.
The simulation results of the CropSPAC model and the SPAC hydrothermal transport model under water deficit conditions (W0 treatment) were, thus, compared and analysed (see Figures 17  and 18).The crop evapotranspiration simulated by the coupled CropSPAC model was less than the SPAC module, while the simulated value of soil water storage in 0-1 m was larger than the SPAC module.It means under deficit irrigation, the crop module of the CropSPAC model could predict the reduction of the LAI of the crop and feed it back to the SPAC module.While in the detached SPAC model, crop information, such as LAI, plant height, root growth, etc., were specified as input data, which does not change when water deficit occurs.Therefore, in case the actual LAI decreased, the detached SPAC module still used the given value (this value generally is greater than the actual LAI), which led to the larger crop transpiration in detached SPAC module.Accordingly, the root water absorption was also greater than the actual value, and then the soil water storage capacity was less than the actual value.Therefore, CropSPAC model that takes into account the dynamic changes in crop growth can better describe and simulate the dynamic changes of soil moisture and evapotranspiration under water stress.

Comparison between CropSPAC and SPAC
CropSPAC model could also outcompete the detached SPAC module for water-heat transfer that ignores the influence of crop growth.SPAC model only accounts for the water-heat transfer from soils to plants, and then to atmosphere, while the crop growth and canopy development is not taken into account.Therefore the detached module of SPAC is inadequate because there is a close relationship between crop growth and water-heat transfer, where the leaf extension and root development of crops will actively affect energy distribution in the canopy and root water absorption.
The simulation results of the CropSPAC model and the SPAC hydrothermal transport model under water deficit conditions (W0 treatment) were, thus, compared and analysed (see Figures 17  and 18).The crop evapotranspiration simulated by the coupled CropSPAC model was less than the SPAC module, while the simulated value of soil water storage in 0-1 m was larger than the SPAC module.It means under deficit irrigation, the crop module of the CropSPAC model could predict the reduction of the LAI of the crop and feed it back to the SPAC module.While in the detached SPAC model, crop information, such as LAI, plant height, root growth, etc., were specified as input data, which does not change when water deficit occurs.Therefore, in case the actual LAI decreased, the detached SPAC module still used the given value (this value generally is greater than the actual LAI), which led to the larger crop transpiration in detached SPAC module.Accordingly, the root water absorption was also greater than the actual value, and then the soil water storage capacity was less than the actual value.Therefore, CropSPAC model that takes into account the dynamic changes in crop growth can better describe and simulate the dynamic changes of soil moisture and evapotranspiration under water stress.

Conclusions and Future Work
In this study, we present a new coupled model (CropSPAC) in which the mechanism of interactions between SPAC water-heat transfer and crop growth is considered.Winter wheat, as one of the most important food crops in North China, is affected directly by the soil hydrothermal conditions, and meteorological conditions, such as solar radiation, air temperature, air humidity, wind speed, precipitation, etc., during its growth season and yield formation.Additionally, the crop growth can also affect the water-heat transfer in SPAC.Thus, the CropSPAC model not only reflects the effect of soil moisture on the growth and yield of winter wheat, but also the effect of crop growth on the water-heat transfer process of SPAC.
The model was validated by using field measurement data from the Yongledian Experimental Station (1998)(1999)

Conclusions and Future Work
In this study, we present a new coupled model (CropSPAC) in which the mechanism of interactions between SPAC water-heat transfer and crop growth is considered.Winter wheat, as one of the most important food crops in North China, is affected directly by the soil hydrothermal conditions, and meteorological conditions, such as solar radiation, air temperature, air humidity, wind speed, precipitation, etc., during its growth season and yield formation.Additionally, the crop growth can also affect the water-heat transfer in SPAC.Thus, the CropSPAC model not only reflects the effect of soil moisture on the growth and yield of winter wheat, but also the effect of crop growth on the water-heat transfer process of SPAC.
The model was validated by using field measurement data from the Yongledian Experimental Station (1998-1999) in Beijing.The simulated soil water content, soil temperature at different depths, evapotranspiration, and crop growth index were in reasonable agreement with the measured values.The CropSPAC model can simulate the water-heat transfer process in the soil, and the process of crop growth more accurately, reliably, and effectively compared with the detached crop module and SPAC module.CropSPAC can be used for scenario analysis, whereby varying the irrigation times and irrigation volumes, the best combination of yield, economic benefit, and water use efficiency could be selected.Therefore, the model can provide decision support and a theoretical basis for the water-saving and yield-increasing of winter wheat and evaluation of irrigation water use efficiency.For this purpose, an optimization procedure would be designed for this model in further research.
Note that the current CropSPAC model contains only one crop, i.e., winter wheat.With the deepening of research, various types of crops, such as maize, rice, soybeans, etc., should be taken into account and the model should be validated on more climatic and geographical conditions, e.g" maize, as a C4 crop, is different from C3 crops (e.g., wheat) in phenology development, photosynthesis, and biomass partitioning from wheat (C3 crop) and need to be considered in the future version of CropSPAC.Furthermore, with the advantage of considering both water and heat transfer, as well as energy partitions in the canopy, the coupled model can be potentially used under field mulching conditions, which is widely applied in North China and should be considered in future study.A Net radiation distribution coefficient 0.3973 [35] B Net radiation distribution coefficient 1.036 [35]

Figure 1 .
Figure 1.Coupling structure of the CropSPAC model.θ is the soil water content; Zr is the root distribution; and h is the plant height.

Figure 1 .
Figure 1.Coupling structure of the CropSPAC model.θ is the soil water content; Zr is the root distribution; and h is the plant height.

Figure 2 .
Figure 2. Calculation flowchart of the CropSPAC model.The GAI is the green area index; Evp is the potential transpiration; Es is the soil evaporation rate; G is the downward heat flux in the soil surface; u is the wind speed; T is the air temperature; RH is the relative humidity; and n are the sunshine hours.

Figure 2 .
Figure 2. Calculation flowchart of the CropSPAC model.The GAI is the green area index; E vp is the potential transpiration; E s is the soil evaporation rate; G is the downward heat flux in the soil surface; u is the wind speed; T is the air temperature; RH is the relative humidity; and n are the sunshine hours.

Figure 3 .
Figure 3. Calculation flowchart of the crop module.The EAI is the ear area index; FT, FC, FN, and FH are the temperature influence factor, CO2 concentration influence factor, nitrogen influence factor, and water influence factor, respectively.

Figure 3 .
Figure 3. Calculation flowchart of the crop module.The EAI is the ear area index; FT, FC, FN, and FH are the temperature influence factor, CO 2 concentration influence factor, nitrogen influence factor, and water influence factor, respectively.

Figure 4 .
Figure 4.The calculation flowchart of the SPAC module.

34 where
v is the latent heat flux by leaf transpiration; H v is the sensible heat flux between crop leaves and canopy air; LE s is the latent heat flux by soil surface evaporation; H s is the sensible heat flux between soil and canopy air; G is the downward heat flux in the soil surface; H is the sensible heat flux between canopy air and the above atmosphere (usually considered at the reference height); LE is the latent heat flux between canopy air and the reference height atmosphere.The unit of the above variables is W•m −2 .Schematic diagrams of energy distribution between atmospheric-ground interfaces and resistance of the hydrothermal transmission are shown in Figure 5. Detailed calculation methods are shown in Appendix B. Water 2019, 11, x FOR PEER REVIEW 9 of LEv is the latent heat flux by leaf transpiration; Hv is the sensible heat flux between crop leaves and canopy air; LEs is the latent heat flux by soil surface evaporation; Hs is the sensible heat flux between soil and canopy air; G is the downward heat flux in the soil surface; H is the sensible heat flux between canopy air and the above atmosphere (usually considered at the reference height); LE is the latent heat flux between canopy air and the reference height atmosphere.The unit of the above variables is W•m −2 .Schematic diagrams of energy distribution between atmospheric-ground interfaces and resistance of the hydrothermal transmission are shown in Figure 5. Detailed calculation methods are shown in Appendix B.

Figure 5 .
Figure 5. Schematic diagram of energy distribution in soil surface and the canopy resistance of the hydrothermal transmission.

Figure 5 .
Figure 5. Schematic diagram of energy distribution in soil surface and the canopy resistance of the hydrothermal transmission.

Figure 6 .
Figure 6.Field experiment plot layout.The open circle (○) denote the position of the neutron probe; the solid circle (•) denote the position of soil temperature thermistor; the plus symbol (+) denote the crop growth index monitoring.

Figure 6 .
Figure 6.Field experiment plot layout.The open circle ( ) denote the position of the neutron probe; the solid circle ( ) denote the position of soil temperature thermistor; the plus symbol (+) denote the crop growth index monitoring.

Figure 7 .
Figure 7.Comparison between simulated and measured soil water storage in 0-1 m soil layer in calibration treatments.The error bars are the standard deviation of field measurements in the two or four medium fertilization plots, as shown in Figure 6 (e.g., the error bars in W0 is calculated from W0FM-1 and W0FM-2; W2 from W2-1FM-1, W2-1FM-2, W2-2FM-1, and W2-2FM-2).The P and I represent precipitation and irrigation.

Figure 7 .
Figure 7.Comparison between simulated and measured soil water storage in 0-1 m soil layer in calibration treatments.The error bars are the standard deviation of field measurements in the two or four medium fertilization plots, as shown in Figure 6 (e.g., the error bars in W0 is calculated from W0FM-1 and W0FM-2; W2 from W2-1FM-1, W2-1FM-2, W2-2FM-1, and W2-2FM-2).The P and I represent precipitation and irrigation.

Figure 8 .
Figure 8.Comparison between simulated and measured soil water storage in 0-1 m soil layer in validation treatments.

Figure 8 . 34 Figure 8 .
Figure 8.Comparison between simulated and measured soil water storage in 0-1 m soil layer in validation treatments.Comparison between simulated and measured soil water content of different soil layer are put in Figures 9 and 10 (only W4 and W1 was shown in this article due to the limits of length).Most of the soil profiles simulation results agreed well with the experiment, with the R 2 on March 25, April 27, May 20, and June 2 in W4 of 0.92, 0.24, 0.74, 0.61, March 25, April 27, May 20, and June 5 in W1 0.89, 0.78, 0.23, 0.91, and the corresponding RMSE of 0.01, 0.03, 0.02, 0.03, 0.01, 0.05, 0.03, and 0.02, respectively.

Figure 9 .
Figure 9.Comparison between simulated and measured soil water content of different soil layer in W4 treatment.(a-d) is March 25, April 27, May 20 and June 2, respectively.

Figure 9 .
Figure 9.Comparison between simulated and measured soil water content of different soil layer in W4 treatment.(a-d) is March 25, April 27, May 20 and June 2, respectively.

Figure 10 .
Figure 10.Comparison between simulated and measured soil water content of different soil layer in W1 treatment.(a-d) is March 25, April 27, May 20 and June 5, respectively.

Figure 10 .
Figure 10.Comparison between simulated and measured soil water content of different soil layer in W1 treatment.(a-d) is March 25, April 27, May 20 and June 5, respectively.

Figure 11 .
Figure 11.Comparison between simulated and measured soil temperature at soil depth 10 cm, 40 cm, 60 cm under different water treatments.

Figure 11 .
Figure 11.Comparison between simulated and measured soil temperature at soil depth 10 cm, 40 cm, 60 cm under different water treatments.

Figure 12 .
Figure 12.Comparison of simulated and measured LAI in calibration treatments.Figure 12.Comparison of simulated and measured LAI in calibration treatments.

Figure 12 .
Figure 12.Comparison of simulated and measured LAI in calibration treatments.Figure 12.Comparison of simulated and measured LAI in calibration treatments.

Figure 13 .
Figure 13.Comparison of simulated and measured LAI in validation treatments.

Figure 14 .
Figure 14.Comparison of measured and simulated above-ground dry biomass in calibration treatments.

Figure 13 . 34 Figure 13 .
Figure 13.Comparison of simulated and measured LAI in validation treatments.

Figure 14 .
Figure 14.Comparison of measured and simulated above-ground dry biomass in calibration treatments.

Figure 14 .Figures 14
Figure 14.Comparison of measured and simulated above-ground dry biomass in calibration treatments.

Figure 15 .
Figure 15.Comparison of measured and simulated above-ground dry biomass in validation treatments.

Figure 15 .
Figure 15.Comparison of measured and simulated above-ground dry biomass in validation treatments.

Figure 16 .
Figure 16.Comparison between CropSPAC and crop module simulated D-AGB and LAI in W0 treatment.

Figure 16 .
Figure 16.Comparison between CropSPAC and crop module simulated D-AGB and LAI in W0 treatment.

Figure 17 .
Comparison between CropSPAC and SPAC model simulated daily ET (mm).

Figure 18 .
Figure 18.Comparison between CropSPAC and SPAC model simulated soil water storage in 0-1 m soil layer in W0 treatment.

Table 2 .
Initial conditions of winter wheat.

Table 3 .
Soil water characteristic parameters.Ks is the saturated hydraulic conductivity; θ s is the saturated water content; θ r is the residual moisture content; and α and n are the parameters in the VGM model.

Table 4 .
Comparison between simulated and measured yield under different water treatments.Figures 14 and 15 show the results between the measured and simulated dry aboveground biomass (D-AGB) and Table

Table 4 .
Comparison between simulated and measured yield under different water treatments.

Measured Yield (kg/ha) Simulated Yield (kg/ha) Relative Error
in Beijing.The simulated soil water content, soil temperature at different depths,

Table A2 .
Parameter values of the CropSPAC.