Modeling Coupled Water and Heat Transport in the Root Zone of Winter Wheat under Non-Isothermal Conditions

Temperature is an integral part of soil quality in terms of moisture content; coupling between water and heat can render a soil fertile, and plays a role in water conservation. Although it is widely recognized that both water and heat transport are fundamental factors in the quantification of soil mass and energy balance, their computation is still limited in most models or practical applications in the root zone under non-isothermal conditions. This research was conducted to: (a) implement a fully coupled mathematical model that contains the full coupled process of soil water and heat transport with plants focused on the influence of temperature gradient on soil water redistribution and on the influence of change in soil water movement on soil heat flux transport; (b) verify the mathematical model with detailed field monitoring data; and (c) analyze the accuracy of the model. Results show the high accuracy of the model in predicting the actual changes in soil water content and temperature as a function of time and soil depth. Moreover, the model can accurately reflect changes in soil moisture and heat transfer in different periods. With only a few empirical parameters, the proposed model will serve as guide in the field of surface irrigation.


Introduction
Soil temperature is not stable in the field, especially in northern semi-arid areas where drastic changes in season and temperature occur [1,2].Sensible heat flux, latent heat flux, and net radiation are all affected by constant changes in surface soil temperature, thereby affecting soil heat flux, resulting in changes in soil internal temperature [3,4].The presence of a soil temperature gradient affects not only the transfer of soil water by affecting the physical properties of soil water, but also affects the redistribution of soil water [5,6].The thermodynamic characteristics of soil are therefore altered, consequently affecting heat flow in soil [7].This phenomenon shows that a coupling relationship exists between water and heat transfer in soil [8].Therefore, if a model of soil water movement established under normal isothermal conditions is used to simulate soil water movement under non-isothermal conditions, the results will show a certain degree of deviation.In studying soil water movement, the effect of temperature gradient on soil water movement parameters and of heat flux on soil water movement cannot be ignored [9].
A mathematical description of soil water movement under both pressure head of soil water (isothermal) and soil temperature (thermal) gradients was provided [10].The theory used consists of classical simulation models, applied by many investigators:

The Cropland
The experiments were conducted in the training base of Shanxi Conservancy Technical College (110°41′23″ E, 34°48′27″ N) in Yanhu District, Yuncheng City, Shanxi Province in China, in October 2015-June 2016.Figure 1 is the location of study sites.The average annual temperature of the soil is 13.6 °C, and the area receives an average annual rainfall of 559.3 mm.This area experiences a warm temperate monsoon climate, characterized by strong semi-arid continental monsoon climate.The soil in the study site is loam.Texture composition and physical parameters of the soil are in Tables 1  and 2, respectively.Depth (cm) (g cm −3 ) (cm 3 cm −3 ) (cm 3 cm −3 ) (cm 3 cm −3 ) (m −1 ) 0-20

Measurements
This study used the polyvinyl chloride pipe soil column method and winter wheat.The outer and inner diameters were 20 cm and 18.6 cm, respectively, and the length of the pipe was 300 cm.The bottom of the column was sealed.Soil columns were built to represent the real conditions in the field, as follows: Before the formal experiment, holes which could just accommodate the soil columns were be dug in a specific position of the test field according to the experiment plan.Additionally, the soil dug out from the holes was marked by the depth of soil layers according to the soil bulk density of the field.The soil moisture testing pipe was arranged along the axis of the soil column vessel.Then, the vessel was filled every 5 cm, according to the corresponding density, with the marked soil.The soil was compacted by the tailor-made tamping device, and soil column could be filled layer by layer from the bottom to the top using the same method.It should be pointed out that the surface of each soil layer was burred to ensure the continuity between the two soil layers before filling the upper soil layer.Next, small holes were drilled in the soil column every 20 cm, in a vertical direction.Then, the temperature probes were inserted into the small holes.The others were manufactured by the same method.Finally, the completed soil columns were placed into the corresponding holes.Figure 2 shows the test equipment.Data were collected during the green, jointing, heading, and filling and maturation stages of winter wheat.The land area of the study site is 61 m 2 , in which several wheat plants were grown in the soil column to simulate the climatic conditions of the field.

Measurements
This study used the polyvinyl chloride pipe soil column method and winter wheat.The outer and inner diameters were 20 cm and 18.6 cm, respectively, and the length of the pipe was 300 cm.The bottom of the column was sealed.Soil columns were built to represent the real conditions in the field, as follows: Before the formal experiment, holes which could just accommodate the soil columns were be dug in a specific position of the test field according to the experiment plan.Additionally, the soil dug out from the holes was marked by the depth of soil layers according to the soil bulk density of the field.The soil moisture testing pipe was arranged along the axis of the soil column vessel.Then, the vessel was filled every 5 cm, according to the corresponding density, with the marked soil.The soil was compacted by the tailor-made tamping device, and soil column could be filled layer by layer from the bottom to the top using the same method.It should be pointed out that the surface of each soil layer was burred to ensure the continuity between the two soil layers before filling the upper soil layer.Next, small holes were drilled in the soil column every 20 cm, in a vertical direction.Then, the temperature probes were inserted into the small holes.The others were manufactured by the same method.Finally, the completed soil columns were placed into the corresponding holes.Figure 2 shows the test equipment.Data were collected during the green, jointing, heading, and filling and maturation stages of winter wheat.The land area of the study site is 61 m 2 , in which several wheat plants were grown in the soil column to simulate the climatic conditions of the field.A capacitance probe (Diviner 2000, Sentek (Pty) Ltd., Melbourne, Australia) was used to measure the soil moisture content once every week; if it rained, more observations were conducted.The temperature of each soil layer was measured by a multi-way soil temperature tester.Data were recorded once every 2 h during the observation period.The depth profiles for measuring soil temperature and water content were at 0, 20, 40, 60, 80, 100,120, 140, 160, 180, 200, 220, 240, 260, 280, and 300 cm, and 0, 20, 40, 60, 80, 100, 120, 140, 160, 180, 200, 220, 240, 260, 280, and 300 cm, respectively.Meteorological data were collected by the ADCON Wireless Automatic Weather Monitoring Station, which was located in the field.Meteorological data, including rainfall, temperature, and humidity, were automatically collected every 10 min.The root system was scanned by a J131A Epson Scanner (Seiko Epson Corporation, Nagano, Japan), and a WinRHIZO version 5.0 (Regent Ins., Quebec, QC, Canada) was used to analyze the scans.As a result, root distribution in each soil layer was obtained.Figure 3 shows precipitation and irrigation during the growing period of winter wheat.
Water 2017, 9, 290 5 of 21 distribution in each soil layer was obtained.Figure 3 shows precipitation and irrigation during the growing period of winter wheat.

Description of the Model
Coupled transfer of water and heat in soil mainly occurs in the vertical direction.As such, it is considered a one-dimensional problem.Given that the soil texture is silty clay loam, the sand content is relatively small, and the soil moisture content is relatively high (the minimum moisture content is 0.14 cm 3 cm −3 ) in the cropland, the effect of the gas phase on the soil was ignored.Therefore, local thermodynamic equilibrium conditions were assumed [35,36].

Water Movement Equation
Soil is a typical porous medium.In porous media, the motion of fluid satisfies Darcy's Law.Movement of soil moisture is driven by soil water potential; that is, soil moisture spontaneously moves from a region of high soil water potential to a region of low soil water potential.Moreover, the movement and changes in matter follows the principle of mass conservation.The mass conservation principle is applied as a continuity equation for fluid motion in porous media.In view of the effect of soil temperature gradient on the physical properties of soil moisture, the equation for water movement under the influence of a temperature gradient can be established with Richards equation and the continuity equation.
where h is the soil water potential (cm); z is the spatial coordinates (in a downward direction) (cm); t is the time (min); T is the temperature (change with time t and spatial coordinate z) (°C); C(h) is the specific water capacity; K(h) is the unsaturated soil water conductivity (cm min −1 ); DT is the moisture diffusivity under the influence of temperature gradient; Sr(z,t) is a function of winter wheat root water uptake rate.

Description of the Model
Coupled transfer of water and heat in soil mainly occurs in the vertical direction.As such, it is considered a one-dimensional problem.Given that the soil texture is silty clay loam, the sand content is relatively small, and the soil moisture content is relatively high (the minimum moisture content is 0.14 cm 3 cm −3 ) in the cropland, the effect of the gas phase on the soil was ignored.Therefore, local thermodynamic equilibrium conditions were assumed [35,36].

Water Movement Equation
Soil is a typical porous medium.In porous media, the motion of fluid satisfies Darcy's Law.Movement of soil moisture is driven by soil water potential; that is, soil moisture spontaneously moves from a region of high soil water potential to a region of low soil water potential.Moreover, the movement and changes in matter follows the principle of mass conservation.The mass conservation principle is applied as a continuity equation for fluid motion in porous media.In view of the effect of soil temperature gradient on the physical properties of soil moisture, the equation for water movement under the influence of a temperature gradient can be established with Richards equation and the continuity equation.

C(h)
∂h ∂t Water 2017, 9, 290 where h is the soil water potential (cm); z is the spatial coordinates (in a downward direction) (cm); t is the time (min); T is the temperature (change with time t and spatial coordinate z) ( • C); C(h) is the specific water capacity; K(h) is the unsaturated soil water conductivity (cm min −1 ); D T is the moisture diffusivity under the influence of temperature gradient; S r (z,t) is a function of winter wheat root water uptake rate.

Heat Flux Equation of Soil
Soil water movement affects the TRANSFORMATION of soil heat, and heat flow in the soil follows Fourier's Law.The system follows the law of conservation of energy, and an equation for heat flux influenced by flow can be established.
where C v is the soil volumetric heat capacity (J cm −3 • C −1 ); K h is the soil thermal conductivity (J cm −1 • C −1 min −1 ); C w is the volumetric heat capacity of soil water (J cm −3 • C −1 ); and q is the soil moisture flux (in a downward direction), calculated by the following formula: where D q is the soil thermal diffusivity (10 −3 cm 2 min −1 ); ρ w is the soil water density (g cm −3 ); v is the Darcy velocity of soil water (cm min −1 ), calculated by the following formula: Equations for coupled soil water and heat transfer consist of Formulas (1)-( 4).The mathematical model for the coupled heat transfer of water can be established when the initial and boundary conditions are given.

Determination of Model Variables
The initial conditions of temperature and matric potential at each discretization node in the study area at the beginning of the simulation are expressed as follows: where z is the depth of soil (cm); t is the time (min); h 0 (z) and T 0 (z) are obtained by field observation, whilst values at the nodes without field data are obtained by quadratic interpolation method.The top and bottom boundary conditions are Neumann conditions and Zero-Neumann conditions, respectively.The top boundary conditions are described as follows: The bottom conditions are described as follows: Water 2017, 9, 290 7 of 22 ∂T ∂z = 0z = 300 cm t ≥ 0 (10) where R(t) as rainfall intensity (cm min −1 ); E s (t) is the surface soil evaporation intensity (cm min −1 ); G n (t) is surface heat flux (J cm −2 min −1 ).The relationship between the soil water content and pressure head proposed by Van Genuchten [37] is adopted as follows: The Van Genhten (VG) equation is applied to determine the unsaturated hydraulic conductivity of soil.The unsaturated hydraulic conductivity can be determined from the characteristic curve of soil water and from the saturated hydraulic conductivity, as follows: where θ s and θ r are saturated and residual values of the soil water content (cm 3 cm −3 ), respectively; m and n 0 are shape parameters following the relationship m = 1 − 1/n 0 , and α (m −1 ) is a scaling parameter of the matric potential; K s is the saturated hydraulic conductivity of soil (cm min −1 ); S e , θ s , and θ r follow the relationship: In this study, the characteristic curve of soil water and the saturated hydraulic conductivity were measured with a pressure membrane and a permeability meter, respectively.Parameters α, n 0 , θ s , and θ r were obtained as detailed in the model calibration below.Table 2 shows the character parameters of the soil.
The presence of a temperature gradient alters the moisture diffusion rate D T .An empirical formula of temperature gradient was used to describe the effect of temperature on water content, as follows [38]: The explicit differential equation method was used to separate the equation.According to the estimation obtained from observed data in the laboratory test (i.e., of the water content and temperature of different sections at different times), the discrete equation was solved with programming.
Therefore, D T determination shows that D T is a function of temperature T: The water balance method was used to calculate crop evapotranspiration, according to the irrigation experiment standard [39].Crop evapotranspiration can be calculated by the following equation: where ET is the farmland evapotranspiration during the calculation period (cm); i is the number of soil layer; N is the total number of soil layers; H i is the soil thickness of layer i (cm); θ i1 is the initial soil moisture content of layer i (cm 3 cm −3 ); θ i2 is the end soil moisture content of layer i (cm 3 cm −3 ); M is the irrigation amount during the calculation period (cm); P is the rainfall during the calculation period (cm); W is the groundwater recharge during the calculation period (cm); and D is the water discharge during the calculation period (cm).
Given that one end of the soil column was sealed with a plastic cloth, then W = 0. Given that the water holding capacity of root layer was larger, then D = 0.
Soil evaporation depends on meteorological conditions, basic soil factors, and the characteristics of the crop growth period [40]: where E s is the surface evaporation; θ is the volumetric soil water content (cm 3 cm −3 ); θ p is the crop wilting moisture (cm 3 cm −3 ); θ f is the field capacity (cm 3 cm −3 ); LAI is the leaf area index; and K w is the soil water supply coefficient [41], as shown in the following equation: where A w is the available water in soil, and it can be calculated as follows: E g is the atmospheric evaporation, and it can be calculated as follows [42]: where T d is the daily mean temperature; n is the sunshine hours; N is the duration of possible sunshine; V 10 is the wind speed at 10 cm high altitude.
In this study, ET and E s can be computed after having the measured data.Then, the crop transpiration intensity T r (t) can be calculated as follows: where ET(t) is the evapotranspiration intensity; E s (T) is the surface soil evaporation intensity.The root water uptake rate S r (z,t) is determined through reverse seeking based on the soil water content [43]; S r (z,t) is expressed as follows [44]: where T r (t) is the transpiration intensity of winter wheat (m s −1 ); L(t) is root depth at time t (cm); A, B, and C are the parameters based on the relationship between root density (g cm −3 ) and soil depth (m).
In the whole trial period, the values of A, B, and C are shown in Table 3.Moreover, they were obtained as detailed in the model calibration below.Soil thermal properties mainly include soil heat capacity (C v ), soil thermal conductivity (K h ), and soil thermal diffusivity (D q ).The relationship among the three is as follows: Soil heat capacity C v can be calculated as follows: where θ is the volumetric liquid water content (cm 3 cm −3 ); θ a is the volumetric air content (cm 3 cm −3 ); n is the porosity of the soil; C s , C w , and C a are the heat capacity of solid matter, water, and gas, respectively.The unsteady state method was used to calculate soil thermal diffusivity by conducting the horizontal soil column experiment in the laboratory [45].The results show that the soil thermal diffusivity is a function of water content: K h represents the thermal conductivity of the soil, and is mainly determined from the soil moisture, soil texture, and gaps in the soil.The thermal conductivity can be obtained after having determined soil heat capacity and diffusivity by Equation (23).
In this study, soil surface heat flux G n can be calculated after having the meteorological data according to the energy balance as follows [41]: where R n is the surface net radiation according to the measured data (W cm −2 ), and it is expressed as follows: The short wave radiation, R s , is the arrival of the ground of the total solar radiation (W cm −2 ); r is the surface reflectance; R L is the effective radiation of the ground (W cm −2 ); R A is the sky radiation (W cm −2 ); empirical coefficients a and b are 0.1668 and 0.5105, respectively [46].δ is the Stefan-Boltzmann constant, and the value is 5.670 × 10 −12 (W cm −2 K −4 ); T K is the thermodynamic temperature, T K = 273 + T ( • C); e is the local water vapor pressure, e = 100 pa; H is the sensible heat flux, and it is expressed as follows: where ρ is the air density (kg cm −3 ); c p is the air specific heat capacity at constant pressure (KJ kg −1 K −1 ); T 0 is the surface temperature (K); T z is the temperature of a high degree of Z (generally take 2 m) (K); r ah is the impedance of air to sensible heat flux (s cm −1 ); Z 0 is the rough height according to experience (cm); K a is the Karman constant, and the value is 0.41; u z is the wind speed of a high Water 2017, 9, 290 10 of 22 degree of Z (generally take 2 m) (cm s −1 ); LE t is the potential latent heat transfer flux from the plant layer to atmosphere layer (W cm −2 ).

Calibration
The implicit difference method was used to discretize the governing equations.The water movement equation and the heat flux equation were solved with an iterative Newton-Raphson technique.Before the model could be used to simulate, the relevant parameters had to be determined.The parameters α, n 0 , θ r , θ s , A, B, and C, which were used to compute the unsaturated soil water conductivity of the soil and the root water uptake rate were determined in the calibration process.The least square method was used to establish the objective functions, and the objective functions are expressed as follows: where θ i is the measured water content in the sample i in the experiment of the characteristic curve for soil water (cm 3 cm −3 ); h i is the soil water potential corresponding to the measured water content (cm); θ(h i , X) is the soil water content calculated by Equation ( 11); X 1 is the parameter vector (θ r , θ s , α, n), which will be determined; N 1 is the number of measured point samples.S r i is the root water uptake rate, which was determined through reverse seeking based on the measured water content in the field experiment; S r (T r (t), L(t), z i , X 2 ) is the root water uptake rate calculated by Equation ( 22); T r (t) is the transpiration rate of winter wheat (m s −1 ); L(t) is the root depth at time t along z coordinates (cm); z i is the depth of the sample i; N 2 is the number of measured point samples in the field experiment; X 2 is the parameter vector (A, B, C), which will be determined.
A hybrid genetic algorithm [47] was employed to optimize these parameters in this study because determining them manually would have been a tedious job.The optimized parameters are listed in Table 2.
The incremental time step ∆t and spatial step ∆z adopted were 1 min and 0.01 m, respectively.The tolerant differential errors of h and T between any two consecutive computational times were 0.01 cm H 2 O and 0.01 • C at the steady state, respectively.Moreover, tolerances ε h and ε T of h and T were 0.1 cm H 2 O and 0.1 • C, respectively, to prevent the occurrence of any undesirable numerical instability during the modeling process.
To further test the accuracy of the simulation results obtained using the proposed model, the simulated soil moisture content and the simulated soil temperature were compared with the measured values by using the following formula: where N is the number of measured point samples for variable; y i is the simulated value of a variable (the water content and temperature) in the sample i; x i is the measured value of a variable (the water content and temperature) in the sample i; MRE is the average relative error; MAE is the mean absolute error; MRE max is the maximum relative error; MAE max is the maximum absolute error.

Validation/Prediction
Profiles of soil water content and soil temperature can be simulated based on the calibrated model.Simulations were performed in four important growth periods of winter wheat (greening stage, jointing stage, heading stage, and filling stage).These simulations were performed as follows: (a) During the greening stage, simulation was conducted from 25  The model simulated the soil water content and soil temperature at the end of forecast period; the indications at the start of forecast period were used as the initial values in the model.
Moreover, simulations of the whole growth period were carried out under three conditions; namely, condition one, condition two, and condition three.Under condition one, D T and S r were taken into account in the proposed model; under condition two, S r was taken into account in the heat flux equation, but D T was neglected in the proposed model; and under condition three, D T was taken into account, but S r was neglected in the heat flux equation.

Simulation Results and Analysis
Figures 4 and 5 show the changes in initial, simulated, and observed values for soil moisture content and temperature, respectively, as a function of soil depth in each growing period of winter wheat.The simulated values for a variable (the water content and temperature) matched the measured values well.
Figure 4a-d show that the trend of simulated values for soil water content was the first to increase (peaked at a depth of nearly 50 cm), and then decreased (a valley existed at approximately 100 cm deep), and then gradually increased again (the maximum value was observed at the bottom of the soil column).The trend for simulated value was consistent with the measured soil water content, demonstrating that the model was effective, and that the results were reliable.Moreover, the simulated and measured values for the degree of agreement increased with soil depth.A large difference was observed between the simulated value of water content in soil surface (0-20 cm) and the measured value, mainly because the soil surface was the layer that was most affected by weather.Weather changes affected evapotranspiration in soil surface and crop transpiration, as well as the water content of the surface soil.Moreover, the structure of the soil surface was less stable than that of the deep soil.When the rainfall was frequent or when the amount of rainfall was high, either the surface soil was compacted, or wetting was produced.These phenomena would alter the soil structure [48], resulting in deviation between the simulated and measured surface soil water contents.Furthermore, changes in temperature of the surface soil would affect soil moisture transfer.Compared with the surface soil, the deep soil was less affected by weather conditions, and its structure was relatively stable.Therefore, the water content of the fitting effect was better at demonstrating that the model produces accurate long-term simulation results.
Figure 4a shows that the soil water content decreased gradually with time during the simulated period characterized by the absence of rain and irrigation.This result was observed mainly because some soil moisture evaporated from the surface, while some of the moisture was released from the soil through transpiration.Soil water content decreased with soil depth, mainly because evapotranspiration of soil moisture decreased gradually with soil depth.Moreover, the root length density of winter wheat showed exponential decline with soil depth (the root length accounted for 81% of the total root length at 0-50 cm, and the proportion decreased to 19% at 50-120 cm).The average root water uptake rates at 0-50 cm and 50-120 cm were 0.004 (cm 3 cm −3 day −1 ) and 0.0015 (cm 3 cm −3 day −1 ), respectively.Thus, the water uptake of roots decreased significantly with soil depth during this period.
Figure 4b shows that the simulated value decreased first and then increased slightly compared with the initial value with soil depth when wheat was irrigated only during the simulated period.It was mainly caused by the loss of soil water through evaporation from the surface and through crop transpiration (the root length accounted for 79% of the total root length and the average root water uptake rate of the layer was 0.0075 (cm 3 cm −3 day −1 )) in the shallow soil layer, which is characterized as being <50 cm deep.A greater amount of soil water was transported into the lower soil layer because of the action of gravity.At >50 cm deep, soil water content increased gradually relative to the initial value as a result of transportation of soil moisture from the upper soil layer.However, the increase in soil water content decreased with soil depth, indicating that the lower soil layers were less affected by rainfall compared with the upper layer.
Figure 4d shows that during the simulated period characterized by the presence of rainfall and irrigation, the simulated soil water content increased obviously with time, mainly because considerable soil water exists, and was transported after surface evaporation and crop transpiration and transporting took place.Moreover, the increase in soil water content progressively decreased with soil depth, mainly because the maximum root depth of winter wheat was 240 cm during this period.The percentage of root length to total root length decreased to 65% in the shallow soil layer, which is characterized as being <50 cm deep.However, the root length density of winter wheat showed exponential decline with increase of soil depth.Moreover, the wheat roots in shallow soil became old, and there was little difference between the shallow and deep soil layer on the root vigor; the root activity decreased to 34.15 (µg g −1 h −1 ) in the shallow soil layer, and the root activity increased to 26.56 (µg g −1 h −1 ) in the deep soil layer.Therefore, the water absorption capacity of the deeper roots increased gradually.
Figure 4c shows that the simulated value of soil water content increased with time relative to the initial value.The increase in soil water content progressively decreased with soil depth.At >220 cm deep, the predicted value was nearly similar to the initial value.This was because the maximum root depth of winter wheat was 220 cm during this period, and the lower soil layers were less affected by crop root water uptake and rainfall than the upper layer.
Figure 5a-d show that the simulated soil temperature in each soil layer increased as growth progressed.This trend was consistent with the reality, demonstrating that the model was effective, and that the results were reliable.The difference between the simulated and observed values decreased with soil depth.This was mainly attributed to factors such as changes in climate, day and night alternation, rainfall, or accumulation of calculation error.
Figure 5a shows that the simulated value increased gradually with soil depth during the green period.The temperature of the topsoil was lowest at the greening stage; the soil temperature gradually increased with soil depth, reaching the maximum at the bottom of the soil column.This was because the average atmospheric temperature was 1.8 • C lower than the topsoil temperature, and the heat eventually transferred from the bottom to the upper soil layer.Then, the air temperature gradually increased with time.Thus, the simulated soil temperature was higher than the initial value.
period.The percentage of root length to total root length decreased to 65% in the shallow soil layer, which is characterized as being <50 cm deep.However, the root length density of winter wheat showed exponential decline with increase of soil depth.Moreover, the wheat roots in shallow soil became old, and there was little difference between the shallow and deep soil layer on the root vigor; the root activity decreased to 34.15 (µg g −1 h −1 ) in the shallow soil layer, and the root activity increased to 26.56 (µg g −1 h −1 ) in the deep soil layer.Therefore, the water absorption capacity of the deeper roots increased gradually.Figure 5a-d show that the simulated soil temperature in each soil layer increased as growth progressed.This trend was consistent with the reality, demonstrating that the model was effective, and that the results were reliable.The difference between the simulated and observed values decreased with soil depth.This was mainly attributed to factors such as changes in climate, day and night alternation, rainfall, or accumulation of calculation error.
Figure 5a shows that the simulated value increased gradually with soil depth during the green period.The temperature of the topsoil was lowest at the greening stage; the soil temperature gradually increased with soil depth, reaching the maximum at the bottom of the soil column.This was because the average atmospheric temperature was 1.8 °C lower than the topsoil temperature, and the heat eventually transferred from the bottom to the upper soil layer.Then, the air temperature gradually increased with time.Thus, the simulated soil temperature was higher than the initial value.Figure 5b shows that the temperature of the topsoil and the bottom of the soil column did not considerably differ (only 3.3 • C difference).The soil temperature decreased first and then increased gradually with soil depth, forming a C-shaped curve.The result is mainly caused by the average atmospheric temperature, which was 3.3 • C higher than the topsoil temperature.The topsoil temperature rose because of solar radiation.The heat eventually transferred to the lower soil layer and thus the temperature of the topsoil decreased gradually.During heat transfer, the moisture-which displays a heat capacity higher than that of the soil-permeated the lower soil.Thus, at <120 cm deep, the temperature of the lower soil remained lower than that of the upper soil layer, although the temperature of the lower soil layer increased.At >120 cm deep, the soil temperature rose gradually because of the heat coming from the bottom of the soil column.Eventually, the atmospheric temperature rose gradually.Thus, the simulated value was higher than the initial value.
Figure 5c shows that the soil temperature decreased slowly and the reduction rate gradually decreased with soil depth during the wheat heading stage.This was because the average atmospheric temperature was 7.8 • C higher than the topsoil temperature.Moreover, the topsoil temperature was increased by solar radiation.The heat eventually transferred to the lower soil layer.However, the lower soil layer obtained less energy than the upper layer because of the energy dissipation in the process of heat transfer.Conversely, the soil moisture increased with soil depth, which is shown in Figure 4.The upper soil temperature increased rapidly because the heat capacity of water was larger than that of the soil.At >220 cm deep, the simulated and initial values nearly formed a curve, indicating that the temperature of the soil layer was stable during the period.
Figure 5d shows that in the filling stage, the soil temperature decreased gradually with soil depth, and the reduction rate also gradually decreased with soil depth.This phenomenon was similar to that in the heading stage.However, the temperature difference between the topsoil and the bottom of the soil column was smaller in the filling stage than in heading stage, as shown in Figure 5c.This phenomenon was observed because of the higher rainfall and irrigation during this period (the total rainfall was 45.5 mm, and the irrigation amount was 141.62 mm), resulting in increased soil moisture.Moreover, soil water transfer directly affected the soil temperature.Thus, the increase in temperature was slightly lower in the surface soil than in the bottom of the soil column.
and the reduction rate also gradually decreased with soil depth.This phenomenon was similar to that in the heading stage.However, the temperature difference between the topsoil and the bottom of the soil column was smaller in the filling stage than in heading stage, as shown in Figure 5c.This phenomenon was observed because of the higher rainfall and irrigation during this period (the total rainfall was 45.5 mm, and the irrigation amount was 141.62 mm), resulting in increased soil moisture.Moreover, soil water transfer directly affected the soil temperature.Thus, the increase in temperature was slightly lower in the surface soil than in the bottom of the soil column.

Error Analysis of Simulated and Observed Values
The relative and absolute errors of the simulated and observed values for soil water content and temperature in different time points were calculated, and the results are shown in Tables 4 and 5.
Tables 4 and 5 demonstrate the accuracy of the simulation result for the dynamic changes in soil moisture and temperature, respectively, in the soil column based on root water uptake by using the proposed mathematical model for soil water heat transfer.The average MRE max of the simulated values for soil water content in the soil column is 5.57%, the average MRE is 1.84%, the average MAE max is 0.0108 (cm 3 cm −3 ), and the average MAE is 0.0039 (cm 3 cm −3 ).The maximum relative error is mainly caused by the unstable soil water content in the surface soil.The relative error decreased with soil depth.The average MRE max of the simulated values for soil temperature in the soil column is 8.42%, the average MAE is 2.48%, the average MAE max is 0.9320 ( • C), and the average MAE is 0.1941 ( • C).The changes in climate, weather, and day and night alternation affected the topsoil temperature, resulting in the presence of maximum relative error.The relative error decreased with soil depth.

Sensitivity Analysis
Given the effect of temperature gradient on the water movement and the effect of root water uptake on the heat transfer, the simulated curves of the soil water content and temperature were obtained under the three conditions.Table 6 shows the error between simulated and observed values for soil water content and temperature under different conditions.The results are shown as follows: Figure 6a-c show that whether or not the effect of temperature gradient on the water movement was taken into consideration, the simulated curves showed roughly the same trend.The simulated curves show that the soil water content decreased with time, but increased when there was rainfall or irrigation.Compared with the measured values at the depths of 20 cm, 140 cm, and 300 cm, the maximum relative error of the soil water content under condition two was 1.716 times, 1.851 times, and 1.543 times, respectively, all of which were greater than that of condition one; the average relative error of the soil water content under condition two was 1.614 times, 1.702 times, and 1.730 times, respectively, all of which were greater than that of condition one.These indicate that the simulated accuracy for the water content under condition one was higher than condition two.
or irrigation.Compared with the measured values at the depths of 20 cm, 140 cm, and 300 cm, the maximum relative error of the soil water content under condition two was 1.716 times, 1.851 times, and 1.543 times, respectively, all of which were greater than that of condition one; the average relative error of the soil water content under condition two was 1.614 times, 1.702 times, and 1.730 times, respectively, all of which were greater than that of condition one.These indicate that the simulated accuracy for the water content under condition one was higher than condition two.Moreover, at depths of 20 cm, 140 cm, and 300 cm, the maximum relative error of the soil water content under condition one was 4.467%, 2.508%, and 1.357%, respectively, all of which were lower than that of condition two; and the average relative error of the soil water content under condition one was 2.427%, 1.262%, and 0.528%, respectively, all of which were lower than that of condition two.These show that the difference of the simulated curves under the two conditions and the influence of temperature gradient on the water movement decreased with soil depth.This was because at depths of 20 cm, 140 cm, and 300 cm, the highest soil temperature was higher than the minimum soil temperature, measuring at 17.13 °C, 9.5 °C, and 4.3 °C, respectively, throughout the whole simulated period.In the topsoil, the surface soil temperature gradient was large because of the diurnal and seasonal changes.Therefore, the effect of the temperature gradient on water movement was greater than that of the lower soil depth.
Figure 7a-c show that the two simulated curves of the soil temperature under the two conditions coincide during the simulated period.This shows that there was little difference on the simulated curve regardless of whether the effect of the temperature gradient on the flow was taken into account.At depths of 20 cm, 140 cm, and 300 cm, the maximum relative error of the soil temperature under condition two was 1.014 times, 1.062 times and 1.062 times, respectively, all of which were greater than that of condition one; and the average relative error of the soil temperature under Moreover, at depths of 20 cm, 140 cm, and 300 cm, the maximum relative error of the soil water content under condition one was 4.467%, 2.508%, and 1.357%, respectively, all of which were lower than that of condition two; and the average relative error of the soil water content under condition one was 2.427%, 1.262%, and 0.528%, respectively, all of which were lower than that of condition two.These show that the difference of the simulated curves under the two conditions and the influence of temperature gradient on the water movement decreased with soil depth.This was because at depths of 20 cm, 140 cm, and 300 cm, the highest soil temperature was higher than the minimum soil temperature, measuring at 17.13 • C, 9.5 • C, and 4.3 • C, respectively, throughout the whole simulated period.In the topsoil, the surface soil temperature gradient was large because of the diurnal and seasonal changes.Therefore, the effect of the temperature gradient on water movement was greater than that of the lower soil depth.
Figure 7a-c show that the two simulated curves of the soil temperature under the two conditions coincide during the simulated period.This shows that there was little difference on the simulated curve regardless of whether the effect of the temperature gradient on the flow was taken into account.At depths of 20 cm, 140 cm, and 300 cm, the maximum relative error of the soil temperature under condition two was 1.014 times, 1.062 times and 1.062 times, respectively, all of which were greater than that of condition one; and the average relative error of the soil temperature under condition two was 1.020 times, 1.080 times and 1.095 times, respectively, all of which were greater than that of condition one.These indicate that the simulated accuracy for the soil temperature under condition one was higher than condition two.However, this advantage was not obvious.This was mainly because thermal parameters can be changed by the change of water content, which resulted in the change of soil temperature.However, the simulated values for the water content of the maximum difference was 0.0352 (cm 3 cm −3 ) under both conditions throughout the whole simulated period.With this difference, the thermal conductivity and diffusivity would be changed by 1.5 × 10 −3 (J cm −1 • C −1 min −1 ) and 0.66 × 10 −5 (cm 2 min −1 ), respectively.This small difference did not cause significant changes in the soil temperature.
condition two was 1.020 times, 1.080 times and 1.095 times, respectively, all of which were greater than that of condition one.These indicate that the simulated accuracy for the soil temperature under condition one was higher than condition two.However, this advantage was not obvious.This was mainly because thermal parameters can be changed by the change of water content, which resulted in the change of soil temperature.However, the simulated values for the water content of the maximum difference was 0.0352 (cm 3 cm −3 ) under both conditions throughout the whole simulated period.With this difference, the thermal conductivity and diffusivity would be changed by 1.5 × 10 −3 (J cm −1 °C−1 min −1 ) and 0.66 × 10 −5 (cm 2 min −1 ), respectively.This small difference did not cause significant changes in the soil temperature.Figure 8a-c show that the two simulated curves of the soil water content under both conditions (conditions one and three) coincide during the simulated period.This shows that there was little difference on the simulated curve regardless of whether the effect of heat conduction in the root water uptake on the heat transfer equation was considered.At depths of 20 cm, 140 cm, and 300 cm, the maximum relative error of the soil water content under condition three was 1.015 times, 1.026 times, and 1.073 times, respectively, all of which were greater than that of condition one; the average relative error of the soil water content under condition three was 1.018 times, 1.048 times, and 1.047 times, respectively, all of which were greater than that of condition one.These indicate that the simulated accuracy for the soil water content under condition one was higher than that of condition three.However, this advantage was not obvious.This was mainly because the soil water movement could be affected by the change of soil temperature.However, the simulated values for the soil temperature of the maximum difference was 0.19 °C under these two conditions throughout the whole simulated period.The soil water diffusivity would be changed by 3.710 × 10 −6 °C min −1 with this difference, which did not cause significant changes in the soil water content.Figure 8a-c show that the two simulated curves of the soil water content under both conditions (conditions one and three) coincide during the simulated period.This shows that there was little difference on the simulated curve regardless of whether the effect of heat conduction in the root water uptake on the heat transfer equation was considered.At depths of 20 cm, 140 cm, and 300 cm, the maximum relative error of the soil water content under condition three was 1.015 times, 1.026 times, and 1.073 times, respectively, all of which were greater than that of condition one; the average relative error of the soil water content under condition three was 1.018 times, 1.048 times, and 1.047 times, respectively, all of which were greater than that of condition one.These indicate that the simulated accuracy for the soil water content under condition one was higher than that of condition three.However, this advantage was not obvious.This was mainly because the soil water movement could be affected by the change of soil temperature.However, the simulated values for the soil temperature of the maximum difference was 0.19 • C under these two conditions throughout the whole simulated period.The soil water diffusivity would be changed by 3.710 × 10 −6 • C min −1 with this difference, which did not cause significant changes in the soil water content.
Figure 9a-c show that whether or not the effect of heat conduction in the root water uptake on the heat transfer is taken into account, the simulated curves showed roughly the same trend.The simulated curves show that the soil temperature increased with the fluctuation of the external climate trend with time.Compared with the measured values, at depths of 20 cm, 140 cm, and 300 cm, the maximum relative error of the soil temperature under condition three was 1.242 times, 1.685 times, and 1.103 times, respectively, all of which were greater than that of condition one; the average relative error of the soil temperature under condition three was 1.531 times, 1.945 times, and 1.075 times, respectively, all of which were greater than that of condition one.These indicated that the simulated accuracy for the soil temperature under condition one was higher than that of condition three.Figure 9a-c show that whether or not the effect of heat conduction in the root water uptake on the heat transfer is taken into account, the simulated curves showed roughly the same trend.The simulated curves show that the soil temperature increased with the fluctuation of the external climate trend with time.Compared with the measured values, at depths of 20 cm, 140 cm, and 300 cm, the maximum relative error of the soil temperature under condition three was 1.242 times, 1.685 times, and 1.103 times, respectively, all of which were greater than that of condition one; the average relative error of the soil temperature under condition three was 1.531 times, 1.945 times, and 1.075 times, respectively, all of which were greater than that of condition one.These indicated that the simulated accuracy for the soil temperature under condition one was higher than that of condition three.
Moreover, at depths of 20 cm, 140 cm, and 300 cm, the maximum relative error of the soil temperature under condition one was 5.543%, 4.018%, and 0.552%, respectively, all of which were lower than that of condition three; the average relative error of the soil temperature under condition one was 5.121%, 2.302%, and 0.137%, respectively, all of which were lower than that of condition three.These show that the difference of the simulated curves under these two conditions and the influence of heat conduction in the root water uptake on the heat transfer decreased with soil depth.This was because that the winter wheat root was mainly distributed in the range of 0-30 cm in depth, and the range of root length accounted for 64% of the total root length.Moreover, in the greening stage, jointing stage, heating stage, and grain filling and maturity stage, at a depth of 20 cm, the root length density values were 3.063 (cm cm −3 ), 6.755 (cm cm −3 ), 7.889 (cm cm −3 ), and 5.195 (cm cm −3 ), respectively; the root water uptake rate values were 0.004 (cm 3 cm −3 day −1 ), 0.0075 (cm 3 cm −3 day −1 ), 0.0108 (cm 3 cm −3 day −1 ), and 0.0087 (cm 3 cm −3 day −1 ), respectively; at a depth of 140cm, the root length density values were 0 (cm cm −3 ), 0.2295 (cm cm −3 ), 0.4393 (cm cm −3 ), and 0.4432 (cm cm −3 ), respectively; the root water uptake rate values were 0 (cm 3 cm −3 day −1 ), 0.00015 (cm 3 cm −3 day −1 ), 0.00045 (cm 3 cm −3 day −1 ), and 0.00066 (cm 3 cm −3 day −1 ), respectively.There was no winter wheat root at a depth of 300 cm.Thus, the root water uptake rate was always 0 (cm 3 cm −3 day −1 ).In the topsoil, the effect of heat Moreover, at depths of 20 cm, 140 cm, and 300 cm, the maximum relative error of the soil temperature under condition one was 5.543%, 4.018%, and 0.552%, respectively, all of which were lower than that of condition three; the average relative error of the soil temperature under condition one was 5.121%, 2.302%, and 0.137%, respectively, all of which were lower than that of condition three.These show that the difference of the simulated curves under these two conditions and the influence of heat conduction in the root water uptake on the heat transfer decreased with soil depth.This was because that the winter wheat root was mainly distributed in the range of 0-30 cm in depth, and the range of root length accounted for 64% of the total root length.Moreover, in the greening stage, jointing stage, heating stage, and grain filling and maturity stage, at a depth of 20 cm, the root length density values were 3.063 (cm cm −3 ), 6.755 (cm cm −3 ), 7.889 (cm cm −3 ), and 5.195 (cm cm −3 ), respectively; the root water uptake rate values were 0.004 (cm 3 cm −3 day −1 ), 0.0075 (cm 3 cm −3 day −1 ), 0.0108 (cm 3 cm −3 day −1 ), and 0.0087 (cm 3 cm −3 day −1 ), respectively; at a depth of 140cm, the root length density values were 0 (cm cm −3 ), 0.2295 (cm cm −3 ), 0.4393 (cm cm −3 ), and 0.4432 (cm cm −3 ), respectively; the root water uptake rate values were 0 (cm 3 cm −3 day −1 ), 0.00015 (cm 3 cm −3 day −1 ), 0.00045 (cm 3 cm −3 day −1 ), and 0.00066 (cm 3 cm −3 day −1 ), respectively.There was no winter wheat root at a depth of 300 cm.Thus, the root water uptake rate was always 0 (cm 3 cm −3 day −1 ).In the topsoil, the effect of heat conduction in the root water uptake on the heat transfer was more obvious than in the lower soil.The effect of root water uptake on the heat transfer decreased because of the weakened root absorption with soil depth.conduction in the root water uptake on the heat transfer was more obvious than in the lower soil.The effect of root water uptake on the heat transfer decreased because of the weakened root absorption with soil depth.In summary, the results show that the model under condition one demonstrated higher accuracy in simulating the soil water content and temperature than that of conditions two and three.Regardless of whether the effect of temperature gradient on the water movement was taken into account, there was a great influence on the simulated values of soil water content.Moreover, whether or not the effect of heat conduction in the root water uptake on the heat transfer was considered, there was a great influence on the simulated values of soil temperature.Therefore, both the effect of temperature gradient on the water movement and the effect of heat conduction in the root water uptake on the heat transfer should not be ignored.

Conclusions
Based on the relationship between soil water and heat transfer, a mathematical model for the fully-coupled water and heat transfer in soil under non-isothermal conditions was established, considering the influence of temperature gradient on soil water redistribution and the influence of root water uptake on soil heat flux transport.The soil water content and temperature in the column of soil were measured during surface irrigation.The experimental data were used to evaluate the proposed model.The results show that the model demonstrated high accuracy in predicting changes in soil water content and soil temperature in the winter wheat root area as a function of time and soil depth.Moreover, the model can accurately reflect the changes in soil moisture and heat transfer in winter wheat in different periods.With only a few empirical parameters, the proposed model will serve as guide in the field of surface irrigation, and the use of which is of practical significance.
Given that the surface layer is affected by numerous external factors, the error in simulation of water content and temperature is relatively large.Thus, dealing with the upper boundary In summary, the results show that the model under condition one demonstrated higher accuracy in simulating the soil water content and temperature than that of conditions two and three.Regardless of whether the effect of temperature gradient on the water movement was taken into account, there was a great influence on the simulated values of soil water content.Moreover, whether or not the effect of heat conduction in the root water uptake on the heat transfer was considered, there was a great influence on the simulated values of soil temperature.Therefore, both the effect of temperature gradient on the water movement and the effect of heat conduction in the root water uptake on the heat transfer should not be ignored.

Conclusions
Based on the relationship between soil water and heat transfer, a mathematical model for the fully-coupled water and heat transfer in soil under non-isothermal conditions was established, considering the influence of temperature gradient on soil water redistribution and the influence of root water uptake on soil heat flux transport.The soil water content and temperature in the column of soil were measured during surface irrigation.The experimental data were used to evaluate the proposed model.The results show that the model demonstrated high accuracy in predicting changes in soil water content and soil temperature in the winter wheat root area as a function of time and soil depth.Moreover, the model can accurately reflect the changes in soil moisture and heat transfer in winter wheat in different periods.With only a few empirical parameters, the proposed model will serve as guide in the field of surface irrigation, and the use of which is of practical significance.
Given that the surface layer is affected by numerous external factors, the error in simulation of water content and temperature is relatively large.Thus, dealing with the upper boundary conditions, climate change, day and night alternation, and redistribution of the rainfall are the concerns that must be considered to improve the accuracy of the proposed model.

Figure 1 .
Figure 1.Location of study sites.

Figure 3 .
Figure 3. Precipitation and irrigation during the growth period of winter wheat.

Figure 3 .
Figure 3. Precipitation and irrigation during the growth period of winter wheat.

3. 3 .
Determination of Model Parameters 3.3.1.Determination of the Parameters for Soil Water Movement

Figure 4 .
Figure 4. Comparison between simulated and observed values for soil water content in different stages of wheat growth.

Figure 4 .
Figure 4. Comparison between simulated and observed values for soil water content in different stages of wheat growth.

Figure 5 .
Figure 5.Comparison between simulated and observed values for soil temperature in different stages of wheat growth.

Figure 5 .
Figure 5.Comparison between simulated and observed values for soil temperature in different stages of wheat growth.

Figure 6 .
Figure 6.Comparison between simulated and observed values for soil water content with time at different soil depths under conditions one and two.

Figure 6 .
Figure 6.Comparison between simulated and observed values for soil water content with time at different soil depths under conditions one and two.

Figure 7 .
Figure 7.Comparison between simulated and observed values for soil temperature with time at different soil depths under conditions one and two.

Figure 7 .
Figure 7.Comparison between simulated and observed values for soil temperature with time at different soil depths under conditions one and two.

Figure 8 .
Figure 8.Comparison between simulated and observed values for soil water content with time at different soil depths under conditions one and three.

Figure 8 .
Figure 8.Comparison between simulated and observed values for soil water content with time at different soil depths under conditions one and three.

Figure 9 .
Figure 9.Comparison between simulated and observed values for soil temperature with time at different soil depths under conditions one and three.

Figure 9 .
Figure 9.Comparison between simulated and observed values for soil temperature with time at different soil depths under conditions one and three.

Table 1 .
Texture composition of the soil.

Table 2 .
The physical parameters of the soil.

Table 1 .
Texture composition of the soil.
Figure 1.Location of study sites.

Table 2 .
The physical parameters of the soil.

Table 3 .
Parameters of the root water uptake rate.
3.3.2.Determination of the Parameters for Soil Heat Transfer February to 3 March, which was a period characterized by sunny days without rain and irrigation.(b) During the jointing stage, simulation was conducted from 19 March to 26 March, which was characterized by sunny days without rain.Irrigation (70.81 mm) was implemented on 19 March.(c) During the heading stage, simulation was conducted from 18 April to 27 April, which was a period characterized by rainy days; irrigation was not implemented.The amounts of rainfall received in 19, 21, 23, and 26 April were 0.2, 0.6, 6.9, and 5.4 mm, respectively.(d) To demonstrate the accuracy of the long-term forecast of the model during the filling stage, we conducted the simulation from 18 April to 16 May, which was a period characterized by rains, and irrigation was implemented.The amounts of rainfall received in 19, 21, 23, 26 April, and in 7 and 14 May were 0.2, 0.6, 6.9, 5.4, 9.2, and 19.9 mm, respectively.Irrigation was implemented in 28 April and 9 May, and the amount of irrigation for each time was 70.81 mm.

Table 4 .
Error between simulated and observed values for soil water content.MAE: mean absolute error; MAE max : maximum absolute error; MRE: mean relative error; MRE max maximum relative error.

Table 5 .
Error between simulated and observed values for soil temperature.

Table 6 .
Error between simulated and observed values for soil water content and temperature under different conditions.