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

Abstract

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 R² 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 R² ranges between 0.96 and 0.98, and the RMSE between 1.22 °C and 1.94 °C. For LAI, the R² 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.

1. 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 farmland, especially in conditions with great changes in farmland soil moisture and heat (such as deficit irrigation, greenhouses, coverage, etc.), so as to provide guidance for enhancing the crop yields and water use efficiency (WUE).

Crop models were commonly used for quantitative description of the crop growth process, which can also interpret experimental results [4], and as agronomic research tools for crop production management [5,6], land resource evaluation, and agricultural production impact assessment [7,8]. Additionally, crop models have also been used to assess and develop irrigation scheduling strategies under limited available water for increasing crop water productivity [9,10]. Examples of such models that have been used for crop production systems include EPIC [11], WOFOST [12], DNDC [13], DSSAT [14], and AquaCrop [15], etc.

Current crop models normally did not include the influence of energy partitioning on the crop growth and water consumption. However, under the circumstances of deficit irrigation or film mulching, which had been extensively applied in North and Northwest China, it has profound influence on the latent heat consumption through soil surface evaporation, thus affecting the energy partition and the advancement of crop growth stage and the crop water consumption [16]. Therefore, a model that could describe not only the crop growth, but also the energy partition, and the water and heat transfer in the soil and canopy, should be more suitable under such circumstances.

Philips [17] proposed the concept of the Soil-Plant-Atmosphere Continuum (SPAC), considering the water movement in the SPAC system as a continuous process. Based on this concept, the water–heat transfer in SPAC describes the process of the water and heat movement in the soil, moisture, and heat transfer in the canopy, as well as the energy partition and consumption in the soil surface and canopy [18], which is essential for crop growth, formation of yield, and a key subject in agriculture, hydrology, and atmospheric science [19]. SPAC theory laid a theoretical foundation for the study of modern farmland water, which has prompted many researchers to conduct in-depth studies in this field [20,21,22]. Software for simulating the water and heat exchanges near land surface were also developed under this concept, e.g., SWAP [23], for simulating the vertical transport of water, solutes, and heat in variably saturated, cultivated soils, SWAT [24], for assessing the impact of management on water supplies and nonpoint source pollution in watersheds and large river basins, and land surface model (LSMs), for simulating surface energy, water fluxes, and surface energy and water budgets in response to near-surface atmospheric forcing [25,26].

In SPAC model the plant variables was generally considered as static parameters, such as leaf area index (LAI), which was just an input data to drive SPAC model simulation. Those input data were usually obtained from field experiment or literature and did not change with the calculated water and heat condition in SPAC. Actually, the crop growth indices, as the output of the crop model, should be expected to have a noticeable impact on the water–heat transfer in the SPAC system, such as root uptake from the soil and water lost to the atmosphere through transpiration. On the other hand, the crop growth was affected by radiation, temperature, and humidity of the atmosphere, soil water, etc., which was the output of SPAC water–heat transfer model.

Given the limitations of previous crop models and SPAC models as outlined above, it is necessary to develop a new model in which both the growth of crop, the water–heat transfer in SPAC and the interactions between them should be considered. Crop index from crop models provides input data for SPAC, while SPAC provides water-heat condition for crop model. The objectives of this paper are to develop such a coupled model (called CropSPAC), that can better characterize not only the crop growth but also the water–heat transfer in SPAC and their interactions. The model will be calibrated and validated based on the field experiment data and its performance will be evaluated by comparing with the results of the detached crop module and SPAC module, respectively.

2. Model Conceptualization and Formulation

2.1. Coupling of Crop Growth and SPAC Water–Heat Transport

CropSPAC includes two main modules, i.e., crop module for simulating crop growth and SPAC module for simulating water-heat transfer in the SPAC system. The coupling between them is shown in Figure 1.

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.

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

2.2.1. 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]:

$$PDT=PD{T}_t\times BDF$$

(1)

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

$$PD{T}_t=PD{T}_{t-1}+DPE$$

(2)

The daily physiological effect (DPE) can be represented by daily thermal effects (DTE) and daily thermal sensitive (DTS) interactions:

$$DPE=DTE\times DTS$$

(3)

DTS can be calculated from [28]:

$$DTS=\{\begin{array}{ll}PE\times VP\hfill & \left(PDT<PDTTS\right)\\ PE\times VP+\left(1-PE\times VP\right)\times \frac{PDT-PDTTS}{PDTHD-PDTTS}\hfill & \left(PDTTS<PDT<PDTHD\right)\\ 1\hfill & \left(PDT>PDTHD\right)\end{array}$$

(4)

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. Parameter values of the CropSPAC.

Sort	Symbol	Parameter Name	Unit	Value	Source
Genetic parameter	ts	Temperature sensitivity		0.9
	PVT	Physiological vernalization time	day	50
	PS	Photoperiod sensitivity		0.005	[27]
	BDF	Basic developmental factor		0.99
	HI	Harvest index		0.5
	SLA	leaf weight per unit	ha/kg	0.0022
Parameters in water influence factor	θOH	Upper limit of optimum soil water content	cm3/cm3	0.35	Calibrated with the observed data
	θOL	Lower limit of optimum soil water content	cm3/cm3	0.18
	θWP	Wilting point soil content	cm3/cm3	0.05
Influence factor	FC	CO2 concentration influence factor		0.95	[27]
	FN	Nitrogen influence factor		0.95
Photosynthesis	κ	Extinction coefficient		0.6	[27]
	PLMX0	Maximum photosynthetic rate of leaves	kgCO2·ha−1·h−1	40
	ε	Initial utilization efficiency of absorbed light	kgCO2·ha−1·h−1/(J·m−2·s−1)	0.5
	σ	Single leaf dissipation coefficient		0.2
Respiration	T0	Optimum temperature of respiration	°C	25	[27]
	Q10	Temperature coefficient of respiration		2	[27]
	RM(T0)	Sustained respiration coefficient at T0	gCO2/gCO2	0.010	Calibrated with the observed data
	Rg	Growth respiration coefficient	gCO2/gCO2	0.26
	Rp(T0)	Photorespiration coefficient at T0	gCO2/gCO2	0.20
Radiation	A	Surface albedo		0.25	[34]
	μ	Radiation ratio		0.97	[34]
	A	Net radiation distribution coefficient		0.3973	[35]
	B	Net radiation distribution coefficient		1.036	[35]

.

2.2.2. 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]:

$$FG=\min \left(FH,FN\right)·FT·FC·PLM{X}_0·\left[1-e^{\left(-\epsilon \times I_L/PLM{X}_0\right)}\right]$$

(5)

where FG is the photosynthetic rate of leaves, kgCO₂·ha⁻¹·h⁻¹; PLMX₀ is the maximum photosynthetic rate of leaves, kgCO₂·ha⁻¹·h⁻¹; ε is the initial utilization efficiency of absorbed light, kgCO₂·ha⁻¹·h⁻¹·(J·m⁻²·s⁻¹) ⁻¹; I_L is the photosynthetically active radiation intensity at L (depth of crop canopy), J·m⁻²·s⁻¹; 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₂ 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:

$$PND=FDTGA-RM-RG-RP$$

(6)

where PND is net daily assimilation of population, kgCO₂·hm⁻²·day⁻¹; RM is consumption of maintain respiration, kgCO₂·hm⁻²·day⁻¹; RG is consumption of growth respiration, kgCO₂·hm⁻²·day⁻¹; RP is consumption of photorespiration, kgCO₂·hm⁻²·day⁻¹, with details shown in Appendix A.

For daily increment of dry biomass in population:

$$\begin{array}{l}TDRW=\xi \times 0.95\times PND/(1-0.05)\\ \xi =\frac{(C{H}_{2}O)}{\left(C{O}_2\right)}=\frac{30}{44}=0.682\end{array}$$

(7)

where TDRW is daily increment of dry biomass in population, kgDM·hm⁻²·day⁻¹; ξ is conversion coefficients between CO₂ and organic compounds.

For accumulation of dry biomass:

$$W_{DAY,t}=W_{DAY,t-1}+TDRW$$

(8)

where W_DAY is accumulation of dry biomass, kgDM·hm⁻²·day⁻¹.

2.2.3. 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_{TOP}=\frac{0.95}{1+e^{-0.18PDT}}$$

(9)

$$P_{ROOT}=1-P_{TOP}$$

(10)

$$P_{LEAF}=\frac{0.6}{1+78.01e^{-0.118\left(56-PDT\right)}}$$

(11)

$$P_{EAR}=\{\begin{array}{ll}\left(0.11+HI\right){\left[\sin \left(\frac{PDT-7.4}{43.6}\times \frac{\pi }{2}\right)\right]}^{\left(4+HI\right)}\hfill & PDT\ge 7.4\hfill \\ 0\hfill & PDT<7.4\end{array}$$

(12)

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

$$\begin{array}{l}{W}_{TOP,t}=W_{TOP,t-1}+FH\left(P_{TOP,t}\times W_{DAY,t}-W_{TOP,t-1}\right)\\ W_{ROOT,t}=W_{DAY,t}-W_{TOP,t}\\ W_{LEAF,t}=W_{LEAF,t-1}+\left(P_{LEAF,t}\times W_{TOP,t}-W_{LEAF,t-1}\right)\times MIN\left(FH,FN\right)\\ W_{STEM}=P_{STEM}\times W_{TOP}\\ W_{EAR,t}=W_{TOP,t}-W_{LEAF,t}-W{}_{STEM,t}\end{array}$$

(14)

The approximate linear relationship between winter wheat yield (YIELD) and spike weight (W_EAR) was established by Cong [29]:

$$YIELD=HI\times W_{EAR}$$

(15)

Based on previous study [32], the empirical relationship between LAI and leaf dry biomass quality (W_LEAF) was given as:

$$LAI=SLA\times W_{LEAF}$$

(16)

where SLA is leaf weight per unit, general range 0.0022~0.0045 ha·kg⁻¹ [32], and in this study it is calibrated by field data (value shown in Appendix A).

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

2.3.1. 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:

$$R_n=R_g\left(1-A\right)-F$$

(17)

where R_n is the net radiation, W·m⁻²; R_g is total short-wave radiation, W·m⁻²; A is the land surface albedo; F is net long-wave radiation, W·m⁻². 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:

$$R_n=R_v+R_s$$

(18)

where R_v is the net radiation absorbed by crop canopy, W·m⁻²; R_s is the net radiation absorbed by ground surface, W·m⁻².

The R_s is calculated from the empirical formula summarized from field observations [35]:

$$R_s=R_n\exp \left[-B\left(1+A\left|\sin \frac{\pi (t-SN)}{12}\right|\right)LAI\right]$$

(19)

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]:

$$R_v=L{E}_v+H_v$$

(20)

$$R_s=L{E}_s+H_s+G$$

(21)

$$H=H_v+H_s$$

(22)

$$LE=L{E}_v+L{E}_s$$

(23)

where LE_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⁻². 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.

2.3.2. 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:

$$E_v=E_{vp}k\left(\psi \right)$$

(24)

$$E_{vp}=\frac{\rho C_p\left(e_{vs}-e_b\right)}{L\gamma \left(r_{vb}+r_{v0}\right)}$$

(25)

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

$$L{E}_v=\frac{\rho C_p\left(e_{vs}-e_b\right)}{\gamma \left(r_{vb}+r_{v0}\right)}k\left(\psi \right)$$

(26)

where r_v0 is the canopy total stomatal resistance under no soil water stress, s·m⁻¹.

2.3.3. 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:

$$\{\begin{array}{l}\frac{\partial \theta }{\partial t}=\frac{\partial }{\partial z}\left(D_w\frac{\partial \theta }{\partial z}\right)+\frac{\partial }{\partial z}\left(D_{wh}\frac{\partial \theta }{\partial z}\right)-\frac{\partial K_w}{\partial z}-s\left(z,t\right)\\ C_v\frac{\partial T}{\partial z}=\frac{\partial }{\partial z}\left(K_h\frac{\partial T}{\partial z}\right)-S_h\end{array}$$

(27)

where θ is the soil water content, cm³·cm⁻³; 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²·s⁻¹; D_wh are the diffusion coefficients of temperature gradients to water flow, m²·s⁻¹·°C⁻¹; K_w is soil water conductivity, m·s⁻¹; C_v is soil volumetric heat capacity, J·m⁻²·°C⁻¹; K_h is the soil thermal conductivity, J· s⁻¹·°C⁻¹; s(z,t) is the distribution function of the actual water uptake rate of roots, m³·m⁻³·s⁻¹; S_h is the heat source sink term, J·m⁻²·°C⁻¹, generally negligible.

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₀(z, t) multiplied by water stress response function α (ψ, π) [39]:

$$s(z,t)=s_0(z,t)\alpha (\psi ,\pi )$$

(28)

$$s_0(z,t)=E_{vp}{\zeta }^{\prime }(z,t)$$

(29)

$${\zeta }^{’}(z,t)=\frac{\zeta (z,t)}{{\displaystyle {\int }_0^{z_r}\zeta (z,t)dz}}$$

(30)

$$\alpha \left(\psi ,\pi \right)=\frac{1}{1+{\left(\frac{\psi (z)+\pi (z)}{{\psi }_{50}}\right)}^P}$$

(31)

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, ψ₅₀ is the specific soil water potential, when E_vp decreases by 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:

$$s(z,t)=\frac{E_{vp}(t){\zeta }^{\prime }(z,t)}{1+{\left(\frac{\psi (z)+\pi (z)}{{\psi }_{50}}\right)}^P}$$

(32)

$$E_v(t)=\frac{E_{vp}(t){\displaystyle {\int }_0^{z_r}{\zeta }^{\prime }(z,t)}}{\left[1+{\left(\frac{\psi (z)+\pi (z)}{{\psi }_{50}}\right)}^P\right]dz}$$

(33)

$$k(\psi )=\frac{{\displaystyle {\int }_0^{z_r}{\zeta }^{\prime }(z,t)}}{\left[1+{\left(\frac{\psi (z)+\pi (z)}{{\psi }_{50}}\right)}^P\right]dz}$$

(34)

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.

3. Experiment and Model Input

3.1. 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⁻³, 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² (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. Irrigation schedule.

Treatment	Fertilization (kg/ha)	Irrigation Volume (mm)				Total Irrigation (mm)
	21 April 1999	March 26	April 21	May 4	May 19
W0	150	-	-	-	-	0
W1	150	-	60	-	-	60
W2	150	-	60	-	50	110
W3	150	60	60	-	50	170
W4	150	60	60	60	50	230

.

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. Initial conditions of winter wheat.

Initial Conditions	Symbol	Unit	Value
Physiological development time	PDT		8
Vernalization	VP		1
Dry biomass amount	Wday	kg/ha	2585
Dry aboveground biomass amount	Wtop	kg/ha	2226
Dry leaf biomass amount	Wleaf	kg/ha	601

in Section 3.2).

3.2. 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 water characteristic parameters.

Soil Texture	Ks (cm/min)	θs	θr	α (cm−1)	n
Sandy loam	0.02	0.48	0.05	0.02	1.34

Note that 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.

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

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

4.1. Soil Water Content

The calibration and validation results of soil water storage in the 0–1 m soil layer are shown in Figure 7 and Figure 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² 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.

Comparison between simulated and measured soil water content of different soil layer are put in Figure 9 and Figure 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² 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.

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., Figure 9d and Figure 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.

4.2. 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² 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.

4.3. Leaf Area Index

Figure 12 and Figure 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² 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 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 Figure 13 and Figure 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.4. Biomass and Yield

Figure 14 and Figure 15 show the results between the measured and simulated dry aboveground biomass (D-AGB) and

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

Treatments	Measured Yield (kg/ha)	Simulated Yield (kg/ha)	Relative Error
W0	3395	3987	17.4%
W1	4722	4192	11.2%
W2	5378	4518	15.9%
W3	5641	4925	12.6%
W4	5340	4937	7.5%

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

5. Discussion on the Coupling Effect of Model

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

5.2. 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 Figure 17 and Figure 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.

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