Soil & Tillage Research 244 (2024) 106259

Contents lists available at ScienceDirect

Soil & Tillage Research
journal homepage: www.elsevier.com/locate/still

Evaluating the effects of different irrigation water sources on soil
temperature using HYDRUS (2D/3D) and considering the coupled
movement of water and heat
Yuehong Zhang a, c , Xianyue Li a, c, * , Jiří Šimůnek b , Ning Chen a, ** , Qi Hu a , Haibin Shi aaCollege of Water Conservancy and Civil Engineering, Inner Mongolia Agricultural University, Huhhot 010018, China
Department of Environmental Sciences, University of California Riverside, Riverside, CA 92521, USA
Collaborative Innovation Center for Integrated Management of Water Resources and Water Environment in the Inner Mongolia Reaches of the Yellow River, Hohhot,
China
b
c

A R T I C L E I N F O

A B S T R A C T

Keywords:
Alternate use of surface water and groundwater
Daily soil temperature fluctuations
Soil temperature distribution
Irrigation schedule
The HYDRUS (2D/3D) software

Drip irrigation with alternate use of surface water and groundwater (ADI) has been widely applied in arid regions
to relieve the effects of heat stress on crop growth. However, the heat dynamics under ADI are still unclear,especially concerning the impacts of ADI on daily and seasonal fluctuations of soil temperature (Ts). Thus, a twoyear experiment was carried out during 2019–2020 in the Jiuzhuang comprehensive saving-water experimental
station to continuously monitor soil water content (SWC) and Ts variations. Moreover, the HYDRUS (2D/3D)
software was used to simulate Ts fluctuations under various evaluated scenarios involving a) surface water
irrigation (SW), b) groundwater water irrigation (GW), c) alternate use of groundwater and surface water irrigation (1G1S), d) two groundwater irrigations and one surface water irrigation (2G1S), and e) three groundwater
irrigations and one surface water irrigation (3G1S). The result showed that the HYDRUS (2D/3D) software could
precisely simulate soil water content and Ts dynamics under all irrigation treatments, with the root mean squareerror of 0.01–0.06 cm3 cm− 3 and 1.25–1.57 ◦ C for SWC and Ts in the verification period, respectively. Apparent
spatial-temporal differences in diurnal Ts fluctuations under different ADI treatments were found, especially inthe 5 cm soil depth. In general, Ts decreased in response to an increase in the frequency of groundwater irrigation. The lowest Ts occurred in the 3G1S treatment under different ADI treatments. The average Ts in both
years under 3G1S was 12.2 % lower than under SW and 4.4 % higher than under GW. However, the highest Ts
occurred in the 1G1S treatment under different ADI treatments. Average Ts in both years under 1G1S increasedby 8.6 % and 15.4 % compared to 2G1S and 3G1S, respectively. Meanwhile, the difference in Ts fluctuations
under different ADI treatments during daytime was substantially higher than during nighttime. The largest area
(1271.8 cm2) of “moderate Ts” (20–22 ℃) occurred in the 2G1S treatment. Moreover, the longest “optimal Ts”
duration occurred for the 22.5 mm irrigation depth under 2G1S. Therefore, the irrigation depth of 22.5 mm andthe 2G1S treatment is recommended as the optimal irrigation strategy in this region.

1. Introduction
The technology combining drip irrigation and plastic film mulchinghas been widely applied in agricultural production worldwide due to its
obvious advantages, such as improving soil temperature (Ts) (Lee andPark, 2020), reducing soil evaporation (Ibrahim et al., 2020), andincreasing soil water content (Lamptey et al., 2019). Moreover, thetechnology could effectively restrain weed growth and increase crop

yield (Sahai et al., 2023; Scarascia-Mugnozza et al., 2011; Menossi et al.,2021). However, Ts under plastic film mulching could exceed the variation range of “optimal Ts” during the late crop growth stage, restricting
root water uptake due to heat stress (Moreno and Moreno, 2008).
Therefore, to relieve the effects of heat stress on crop yield,
high-frequency irrigation was often adopted to decrease Ts during highair temperature periods (Garg and Cheema. 2011).
The reasonable irrigation strategy not only maintained the “optimal* Corresponding author at: College of Water Conservancy and Civil Engineering, Inner Mongolia Agricultural University, Huhhot 010018, China.
** Corresponding author.
E-mail addresses: lixianyue80@126.com (X. Li), chenning0801@126.com (N. Chen).
https://doi.org/10.1016/j.still.2024.106259
Received 21 September 2023; Received in revised form 16 July 2024; Accepted 29 July 2024
Available online 31 July 2024
0167-1987/© 2024 Elsevier B.V. All rights are reserved, including those for text and data mining, AI training, and similar technologies.

Y. Zhang et al.

Soil & Tillage Research 244 (2024) 106259

Ts” (Soureshjani et al., 2019), but also promoted root growth and crop
yield (Ierna and Mauromicale, 2022). Ts fluctuations under different
irrigation depths have been thoroughly studied in the previous literature(e.g., Kishore et al., 2022; Filipović et al., 2016; Liao et al., 2016). For
example, Karandish and Shahnazari (2016) compared Ts differences
under full irrigation, partial root-zone drying (PRD) (i.e., moderate
irrigation), and deficit irrigation. They found that the PRD treatment
kept Ts in a favorable range the longest. Jha et al. (2017) reported that Ts
decreased by 4.3 ℃ when the irrigation depth was increased by 20 %.
Additionally, the probability of a reduction in crop yield caused byextremely high temperatures can be decreased by irrigation (Wechsung
et al., 2021). For example, Deb et al. (2022) illustrated that the potential
corn yield increased from 10 % to 259 % under different irrigation water
application rates relative to rainfed crop yield in the drought year.
It is thus vital to control Ts fluctuations by formulating a reasonableirrigation schedule. However, few studies have focused on Ts fluctuations under drip irrigation using alternate surface water and groundwater (ADI, alternate drip irrigation) sources. ADI is an innovative
irrigation technology that can be used to decrease Ts by increasing theirrigation frequency with groundwater due to its lower water temperature. Nevertheless, root activity was still restrained when Ts overly
decreased with an increase in the frequency of groundwater irrigation.
Therefore, it is urgent to quantify Ts fluctuations to determine theoptimal irrigation schedule since an unreasonable irrigation strategy
could increase Ts variability and further limit crop root activity.
Ts fluctuations were in the past monitored mainly using field experiments, which were highly time-consuming and costly. On the other
hand, mathematical models could be alternatively used to save time andother resources. Moreover, mathematical models could provide additional details about Ts fluctuations in the soil (Bristow et al., 1986;
Chung and Horton, 1987; Mehdizadeh et al., 2020). Among various
mathematical models, the HYDRUS software has been widely used tosimulate various processes in soils, including soil temperature dynamics,
due to its flexibility in handling various boundary conditions, friendly
graphical user interface, and ability to accurately reproduce field results.
For instance, Rai et al. (2019) used the HYDRUS software to evaluate theimpacts of different agricultural conservation practices on Ts in the soil
profile and found that the “permanent broad bed with crop residue”
strategy could increase Ts. Kodešová et al. (2014) simulated soil heat
dynamics under various surface covers using the HYDRUS software andillustrated that the highest Ts during the cold period occurred for thebasalt gravel mulching treatment. Kassaye et al. (2021) reported thatsoil heat was transported to a depth of about 70 cm through an Andosol
vadose zone under a temperate monsoon climate.
However, previous studies focused only on simulating Ts under
different irrigation and mulching scenarios. For instance, Daniel et al.
(2022) demonstrated that Ts at 10, 20, 30, and 40 cm soil depths in an
irrigated field were 1.3, 0.26, 0.59, and 0.4 ℃ lower than in a rainfed
field. Matsue et al. (2021) indicated that Ts under the full irrigation
treatment decreased by 0.8–2.0 ℃ compared to intermittent irrigation.
Chen et al. (2021) found that Ts under plastic film mulching increased by
2.2–4.7 ℃ compared to no mulching. However, the effects of ADI on Ts
fluctuations have not yet been quantified using mathematical models.
Therefore, capturing Ts fluctuation under different ADI treatments andformulating an optimal irrigation schedule is necessary.
Therefore, the objectives of this study are 1) to calibrate and validate
the HYDRUS (2D/3D) software using the observed SWC and Ts data, 2)
to evaluate differences in Ts fluctuations under different ADI treatments,
3) to analyze the impacts of ADI on daily and seasonal Ts fluctuations,and 4) to reveal the two-dimensional distribution of Ts under different
ADI treatments.

Table 1
Different irrigation strategies (i.e., the surface water irrigation (SW), one
groundwater irrigation after one surface water irrigation (1G1S), two groundwater irrigations after one surface water irrigation (2G1S), three groundwater
irrigations after one surface water irrigation (3G1S), and groundwater water
irrigation (GW)) during the growing seasons in 2019 and 2020.
NumberDate

Remarks

1

May 9, 2019, andMay 14, 20202

May 16, 2019, andMay 23, 20203

May 23, 2019, andJune 2, 20204

May 30, 2019, andJunu 15, 2020

5

June 5, 2019, andJune 23, 2020

6

June 12, 2019, andJuly 2, 2020

7

June 20, 2019, andJuly 9, 2020

8

June 29, 2019, andJuly 16, 2020

9July 13, 2019, andJuly 26, 2020

10

July 23, 2019, andAugust 4, 2020

11

August 1, 2019, andAugust 11, 2020

12

Auguset 15, 2019,
August 18, 2020

Surface water was applied in the SW treatment,
and groundwater was applied in the 1G1S,
2G1S, 3G1S, and GW treatments.
Surface water was applied in the SW and 1G1S
treatments, and groundwater was applied in the
2G1S, 3G1S, and GW treatments.
Surface water was applied in the SW and 2G1S
treatments, and groundwater was applied in the
1G1S, 3G1S, and GW treatments.
Surface water was applied in the SW, 1G1S, and
3G1S treatments, and groundwater was appliedin the 2G1S and GW treatments.
Surface water was applied in the SW treatment,
and groundwater was applied in the 1G1S,
2G1S, 3G1S, and GW treatments.
Surface water was applied in the SW, 1G1S, and
2G1S treatments, and groundwater was appliedin the 3G1S and GW treatments.
Surface water was applied in the SW treatment,
and groundwater was applied in the 1G1S,
2G1S, 3G1S, and GW treatments.
Surface water was applied in the SW, 1G1S, and
3G1S treatments, and groundwater was appliedin the 2G1S and GW treatments.
Surface water was applied in the SW and 2G1S
treatments, and groundwater was applied in the
1G1S, 3G1S, and GW treatments.
Surface water was applied in the SW and 1G1S
treatments, and groundwater was applied in the
2G1S, 3G1S, and GW treatments.
Surface water was applied in the SW treatment,
and groundwater was applied in the 1G1S,
2G1S, 3G1S, and GW treatments.
Surface water was applied in the SW, 1G1S,
2G1S, and 3G1S treatments, and groundwater
was applied in the GW treatment.

Note that the irrigation depth (for each irrigation event) and the total irrigation
depth (for the entire season) were 22.5 and 270 mm, respectively.

2. Materials and methods
2.1. Site description
Experimental plots are located at the Jiuzhuang Comprehensive
Saving-Water Experimental Station (40◦41’N, 107◦18’E) in Linhe, Inner
Mongolia, China. The region’s climate is typically meso temperate, withdaily average minimum and maximum air temperatures of − 4.0 ℃ and36.8 ℃ during the crop growth season. The average annual rainfall is116.3 mm, with monthly coefficients of variation above 50 %. The soil
type at this experiment site is sandy loam soil (belonging to Plaggeptsoil) based on the USDA classification system, with moderate fertility
and low organic matter content (pH: 8.0; TOC: 8 g kg− 1; N: 0.8 g kg− 1; P:
10 mg kg− 1; K: 200 mg kg− 1). The average bulk density, permanent
wilting point water content, and field capacity in the 0–100 cm soil layer
were 1.42 g cm− 3, 0.12 cm3 cm− 3, and 0.36 cm3 cm− 3, respectively.
2.2. Experimental layout and treatments
The experimental area of 3332 m2 was divided into three blocks with
3 m spacing. Each block was divided into five plots, and the treatments
were allocated randomly. This resulted in 15 experimental plots, i.e.,
three replicates on five plots. The irrigation treatment included irrigation with surface water (SW) or groundwater (GW) during the entire
crop growth season as control treatments. Moreover, three treatments
2

Y. Zhang et al.

Soil & Tillage Research 244 (2024) 106259

Fig. 1. Potential evapotranspiration, cumulative rainfall, and irrigation strategy (small blue graph) in 2019 and 2020. Five rows in the small blue graph represent
five irrigation strategies. Note that black lines represent surface water irrigation, red lines groundwater irrigation, and pink tics rainfall time.

Fig. 2. A schematic showing the modeling domain, boundary conditions, sensor locations (a), and root distribution (b).

with alternate use of surface water and groundwater were used: one
groundwater irrigation and one surface water irrigation (1G1S), twogroundwater irrigations and one surface water irrigation (2G1S), andthree groundwater irrigations and one surface water irrigation (3G1S).
Irrigations were applied when the soil water content was lower than
0.23 cm3 cm− 3, i.e., 65 % of field capacity. The management allowable
depletion (with MAD=0.525) method (Allen et al., 1998) was used todetermine corn’s irrigation lower limit. Moreover, corn’s recommended
irrigation upper limit was set to 95 % of field capacity. Thus, the irrigation depth was 22.5 mm. The irrigation times and total irrigation
depths in 2019 and 2020 were the same, i.e., 12 times and 270 mm(Table 1). One drip line with an emitter spacing of 30 cm and an irrigation rate of 2.4 L h− 1 was placed in the middle of two cornrows using
an automatic sowing-mulching machine. The average groundwater andsurface water temperatures were 17.1 and 21.6 ℃, respectively.
Groundwater and surface water salinities were 1.9 and 0.2 dS m− 1,
respectively.
Corn seeds were sown on May 1 and 4 and harvested on September
17 and 25 in 2019 and 2020, respectively. The planting depth, density,and row spacing were 3.0 cm, 4.5 seeds m− 2, and 40 cm, respectively. A
60 cm wide and 0.1 cm thick plastic film, produced from polyethylene
using a polymerization reaction, was applied above two cornrows. The
plastic film’s percentage of breaking elongation, maximum load, andtensile strength were 363.8 %, 0.97 N, and 24.1 MPa, respectively. Field
management and fertilization strategy were identical for all plots. Before
sowing, based on the soil analysis recommendation, a basal fertilizer
with 140 kg ha− 1 of diammonium phosphate ((NH4)2HPO4, N≥18 %),140 kg ha− 1 of potassium sulfate (K2SO4), and 50 kg ha− 1 of carbamide(CO(NH2)2, N≥32 %) was applied. The solubilized topdressing fertilizer
was carbamide (CO(NH2)2, N≥32 %), which was applied during theelongation stage (30 % of the total N-fertilizer) and the grain-fillingstage (30 % of the total N-fertilizer) along with irrigation.

2.3. Measurements and methods
2.3.1. Meteorological data
An automatic meteorological station (Onset Computer Inc.; U30,
Hobo, USA) was installed about 50 m from the experimental field.
Rainfall, solar radiation, wind speed, and air temperature were monitored once every hour and averaged every 24 hours (Fig. 1). The highest
solar radiation and air temperatures occurred in July and August, while
the largest wind speed was recorded in May and September. The prevailing winds come from the southeast during May-September. Additionally, intensive rainfalls occurred mainly in August and September.
2.3.2. Soil water content
The soil water content (SWC) was measured every 20 cm between
0 and 100 cm depths using TDR probes (IMKO GmbH Inc.; IPH, TRIMEPICO, Germany) connected to a data logger. TDR probes were installed
between two cornrows and in a bare area (Fig. 2). Measurements were
performed once every 5–7 days. Meanwhile, gravimetric measurements
were adopted to verify the accuracy of TDR probes at regular intervals.
2.3.3. Soil water balance components
The calculation of the soil water balance was based on the principleof mass conservation:
SWS = I + P + R − ET − K(1)where I is irrigation (mm), P is precipitation (mm), R is capillary rise(mm), K is surface runoff (mm), and SWS is the change in soil water
storage within the 0–100 cm soil layer. Surface runoff K could be
neglected because of the negligible flow rate of drip irrigation. Thus, thesoil water balance equation was simplified as follows:
ET = I + P + R − SWS

3(2)

Y. Zhang et al.

Soil & Tillage Research 244 (2024) 106259

SWS = (θt2 − θt1 ) × H

2.4. HYDRUS (2D/3D) software(3)

2.4.1. Governing equations
The HYDRUS (2D/3D) software simulates the movement of soil
water and heat in variably saturated porous media under various
irregular and complex boundary conditions (Šimůnek et al., 2016;2024). Soil water flow is simulated by numerically solving the modified
form of the two-dimensional Richards equation (Šimůnek et al., 2016):
[
]
(
)
∂θ
∂
∂h
Kij (h)(9)+ Kiz − S h, hφ=∂t ∂xi∂xj

where H is the soil depth (cm) and θt1 and θt2 are the soil water contents
at times t1 and t2, respectively (cm3cm− 3).
Capillary rise R was calculated as follows:(Ψ − Ψ
)m2
m1
R = K(θ)+1(4)
Z
where ψ m1 and ψ m1 are the soil water potentials at depths of 90 and
110 cm (hPa), respectively, Z is the vertical distance of 20 cm, and K(θ)
is the saturated hydraulic conductivity (cm day− 1).

where θ is the soil water content (cm3cm− 3), h is the pressure head (cm),
K is the unsaturated hydraulic conductivity (cm day− 1), t is time (day), S
is the root water uptake term (cm3 cm− 3 day− 1), and xi are spatial coordinates (cm) (x and z are the horizontal and vertical coordinates (cm),respectively).
Root water uptake was computed using the Feddes model (Feddes
et al., 1978) that considers separately water and salinity stresses:{
Tpωω > ωc
S(h, hφ , x, z) = α1 (h)α2 (hφ )b(x, z) , ω =(10)ωc ω < ωcω

2.3.4. Soil temperature
Two methods were used to monitor Ts, i.e., thermocouples and anangle thermometer. The thermocouples (Beijing Huize Agriculture
Technologies Inc., HZR-8 T, CHN) have an accuracy of ±0.2–0.5 ℃. The
instruments were installed in each plot between two cornrows and in abare area to continuously monitor Ts at soil depths of 5, 10, 15, and20 cm throughout the day. Data was recorded every two hours andaveraged every 24 hours. An angle thermometer (Beijing Beixin Keyuan
Instrument Inc., HG04-SYQX-2, CHN) was used to calibrate andcompare the results of the thermocouples at regular intervals.
Two indicators, i.e., the average increase (Tu) and decrease (Td) indaily temperatures, were used to evaluate increasing and decreasing
trends in Ts in the tillage layer (Chen et al., 2021):
)
n− 1
1 ∑
Tu =
max(Ti+1 − Ti , 0(5)m1 i=1
n− 1
1 ∑
Td =
min(Ti+1 − Ti , 0m2 i=1

where α1(h) and α2(hφ) are the stress response functions representing
water and osmotic stresses (-), respectively, b(x,z) is a normalized distribution of root density in the two-dimensional soil domain (cm− 2), andω is a dimensionless water stress index. The water stress reduction
function of Feddes et al. (1978) and the threshold and slope function ofMaas (1990) were adopted to represent water and osmotic stresses,
respectively. The root water uptake parameters for corn proposed by
Wesseling and Brandyk (1985) were used in this study. Moreover, corn’ssalinity threshold and slope values were set to 1.7 dS m− 1 and 12 %,respectively (Karandish and Šimůnek, 2019).
Soil heat movement is simulated by numerically solving the twodimensional heat transport equation as follows:[]∂C(θ)T
∂
∂T
∂T
(11)
=
λij (θ)− Cw qi
∂xi
∂t
∂xj
∂xi

)

n − 1 = m1 + m2

(6)
(7)

where Ti is the average daily soil temperature (◦ C), n is the total numberof data points (d), m1 is the total number of days with increasing temperature (d), and m2 is the total number of days with decreasing temperatures (d).

where C(θ) and Cw are the volumetric heat capacities of soil and water (Jcm− 3 ◦ C− 1), respectively, λij(θ) is the soil apparent thermal conductivity(W cm− 1 ◦ C− 1), T is the temperature (◦ C), and qi is the water flux (cmday− 1). C(θ) was determined using the method proposed by de Vries(1963):

2.3.5. Groundwater level and groundwater temperature
The groundwater level and temperature at the experimental field
were monitored using a pressure transducer (Onset Computer Inc., U20,
Hobo, USA). The transducer was installed in the observation well 3 m
deep with a 5 cm internal diameter. Data were averaged and recorded at
one-day intervals and stored in a datalogger.

C(θ) = Cn θn + Co θo + Cw θ

where Cn, Co, and Cw are the volumetric heat capacities of the solid
phase, organic matter, and liquid phase (J m− 3 ◦ C− 1), respectively, andθn and θo are volumetric fractions of the solid phase and organic matter(cm3 cm− 3), respectively.
(
)
qj qi
λij (θ) = λT Cw |q|δij + λL − λT Cw(13)+ λ0 (θ)δij|q|

2.3.6. Leaf area index
Five plants of each species from each plot were randomly sampled todetermine leaf areas using a leaf-area meter (LI-COR Inc. Li- 3000 C, LI,
USA). The leaf area index (LAI) was calculated using the FAO method(Allen et al., 1998). LAI was 1.8, 4.1, 3.8, and 3.3 cm2cm− 2 in theelongation, tasseling, filling, and maturity stages in 2019, respectively,
while it was 2.7, 4.5, 4.2, and 3.5 cm2 cm− 2 in 2020.

where λL and λT are the longitudinal and transverse thermal dispersivities (cm), respectively, δij is the Kronecker delta, and λ0(θ) is thethermal conductivity (W cm− 1 ◦ C− 1) that can be described according toChung and Horton (1987) as follows:

2.3.7. Crop yield and water use efficiency
Ten corn plants from each plot were sampled to measure crop yield(Y) at physiological maturity. Water use efficiency (WUE) was calculated as follows:
WUE =Y
× 100%ET(12)

λ0 (θ) = b1 + b2 + b3 θ0.5

(14)

where b1, b2, and b3 are empirical parameters (W cm− 1 ◦ C− 1).

(8)

2.4.2. Modeling parameters
The soil hydraulic properties were described using the standard
equations of van Genuchten (1980). The soil hydraulic parameters fordifferent soil layers were estimated using the Rosetta Lite functions of
Schaap et al. (2001) and the class pedotransfer function of Carsel andwhere Y is crop yield (kg ha− 1) and ET is evapotranspiration (mm).

4

Y. Zhang et al.

Soil & Tillage Research 244 (2024) 106259

and bare areas (W m− 2), respectively.
Rnm was determined by considering the effects of crop canopy andmulching on net radiation received by soil (Liakatas et al., 1986):)(
Rn = 1 − ρs τS + Lʹn(18)

Table 2
Soil hydraulic parameters (the residual water content θr, the saturated watercontent θs, the shape parameters (α, n, and l), and the saturated hydraulic conductivity Ks) of the 0–20, 20–40, and 40–100 cm soil depths in the experimentfield.
Soil depth(cm)

θr(cm3
cm− 3)

θs(cm3
cm− 3)α(cm− 1)

n(-)

Ks(cm
day− 1)

l(-)

0–20
20–40
40–100

0.03
0.04
0.05

0.41
0.42
0.41

0.02
0.01
0.01

1.45
1.63
1.63

55.4
52.0
37.9

0.5
0.5
0.5

where ρs is the reflection coefficient of the soil (-), τ is the transmission
coefficient of the mulch (-), S is the flux density of solar radiation on ahorizontal surface (W m− 2), andLʹn is net thermal radiation (W m− 2)
received by the soil beneath a mulch calculated using a thermodynamical function proposed by Liakatas et al. (1986).
Since the energy balance equation did not continuously determine
the surface heat flux densities due to discontinuity in evaporation ratedata, the temperature gradient method proposed by Rai et al. (2019) was
adopted to calculate the continuous surface heat flux densities:

Parrish (1988) as implemented in HYDRUS (2D/3D). Estimated soil
hydraulic parameters were further calibrated while considering theirtemperature dependence using HYDRUS “inverse solution” and manualadjustments (Table 2). The thermal parameters were simultaneously
optimized using observed SWCs and Ts (Table 3).
Potential evapotranspiration (ETp) was calculated using the modified
Penman-Monteith equation (Allen et al., 1998) (Fig. 1). Potential transpiration (Tp) and evaporation (Ep) were required in the simulations, andthey were obtained by partitioning of ETp using the dual Kc method(Rosa et al., 2012).

Gb = λvolumeofsolutionapplied/day
Aw

Rnb − Hb − LEb − Gb = 0

(17)

(20)

The performance of HYDRUS (2D/3D) was evaluated using the root
mean square error (RMSE), the mean relative error (M RE ), and themodeling efficiency (EF) to determine the correlation and discrepancy ofsimulated values around observed ones for both the calibration andvalidation stages (Amato, Giménez, 2024; Ebrahimi et al., 2022; Gong
et al., 2023):
[
]1/2
2
1∑n
(S
i − Mi )n i=1

RMSE =

where Aw is the surface wetted area below drippers considered in a timevariable flux boundary condition (cm2).
The time-variable pressure head boundary condition was specified at
the bottom of the flow domain to represent measured groundwater table
depths. In addition, the zero-flux boundary condition was prescribed on
the left and right sides of the domain.
The third-type boundary condition was applied on the upperboundary to represent the actual heat flux in the mulched and bareareas. Heat fluxes in the mulched and bare areas for heat transport
simulations were obtained from the energy balance equations:
(16)

dTb
dZ

2.5. The performance of the HYDRUS (2D/3D)(15)

Rnm − Hm − LEm − Gm = 0

(19)

dTb
m
where dT
dZ and dZ are the temperature gradients in the mulched region
and bare area, respectively, and λ is the thermal conductivity determined by jointly solving Eqs. (19–20) (W cm− 1 ◦ C− 1).

2.4.3. Initial and boundary conditions
The SWCs and Ts measured at the beginning of the simulation period
were taken as the initial conditions. The transport domain was set as arectangle with a width of 50 cm and a depth of 100 cm, and there thus
were four boundaries. The upper boundary was divided into three sections with the time-variable flux, zero-flux, and atmospheric boundary
conditions representing the wetted area below the dripper, the mulchedregion, and the bare area, respectively (Fig. 2). The time-variable flux
with a flux q was defined as (Skaggs et al., 2004):
q=

dTm
dZ

Gm = λ1

MRE = n∑ni=1 |Si −

Mi

Mi |

(21)

× 100%(22)∑n(Si − Mi )2
EF = 1 − ∑ni=12
i=1 (Mi − M)(23)

where n is the number of data points, and Si and Mi are simulated andmeasured values, respectively.
3. Results
3.1. Model evaluation

where Rnm and Rnb are net radiations received by the soil in the mulched
and bare areas (W m− 2), respectively Hm and Hb are sensible heat flux
densities in the mulched and bare areas (W m− 2), respectively, L is thelatent heat of vaporization (equal to 2.45 MJ kg− 1), Em and Eb are theevaporation rates in the mulched and bare areas (kg m− 2 s− 1), respectively, and Gm and Gb are the surface heat flux densities in the mulched

Observed data obtained from 2019 was used to calibrate soil hydraulic and thermal parameters of the HYDRUS (2D/3D) software, andthe corresponding data obtained from 2020 was used to validate theoptimized parameters. The results showed that the HYDRUS (2D/3D)
software could precisely simulate soil water and heat dynamics. In the

Table 3
Soil thermal parameters (the longitudinal thermal dispersivity λL, transverse thermal dispersivity λT, volumetric heat capacity of the solid phase Cn, volumetric heat
capacity of organic matter C0, volumetric heat capacity of the liquid phase Cw, and the empirical parameters b1, b2, b3) of the 0–20, 20–40, and 40–100 cm soil depths in
the experiment field.
Soil depth (cm)
0–20
20–40
40–100

λL

λT

b1

b2

b3

Cn

C0

Cw

(cm)

(cm)

(J m− 1 s− 1 ◦ C− 1)(J m− 1 s− 1 ◦ C− 1)(J m− 1 s− 1 ◦ C− 1)(J m− 3 ◦ C− 1)(J m− 3 ◦ C− 1)(J m− 3 ◦ C− 1)

5.0
6.3
4.4

1.0
1.3
0.9

0.57
1.24
0.55

2.53
3.19
1.69

1.89
1.54
3.75

1.92×106
1.92×106
1.92×106

2.51×106
2.51×106
2.51×106

4.18×106
4.18×106
4.18×106

5

0.86
0.96
0.99
0.80
0.92
0.90
0.85
0.82
0.86
0.77
0.87
0.71
0.88
0.86
0.83
0.80
0.03
0.01
0.04
0.04
0.99
0.94
1.05
0.97
0.05
0.06
0.05
0.04
1.37
1.25
1.38
1.38
0.89
0.95
0.87
0.84
0.91
0.89
0.88
0.85
0.93
0.83
0.80
0.88
0.85
0.84
0.84
0.79

7.94
7.38
8.02
7.30
5.06
4.64
5.23
4.92
8.28
8.04
7.80
8.12
5.39
4.61
5.24
5.11

RMSE
(cm3 cm− 3/℃)
EF
(-)

7.68
7.25
6.75
7.11
5.40
4.74
5.16
4.89
7.66
8.33
9.00
7.58
5.26
4.96
4.74
5.16
8.05
8.27
7.29
7.52
4.81
5.09
4.74
5.28
8.15
8.80
8.39
7.82
5.41
4.84
5.05
4.59

0.99
0.84
0.83
0.86
0.96
0.94
0.95
0.93
0.82
0.81
0.79
0.80
0.86
0.86
0.83
0.84

0.05
0.03
0.06
0.06
1.19
1.16
1.21
1.25
0.04
0.02
0.06
0.04
1.34
1.28
1.35
1.31

7.91
7.32
7.39
7.91
5.10
5.04
5.19
5.32
9.74
9.74
9.02
8.84
5.37
4.94
5.13
4.86

1.00
0.90
0.95
0.81
0.94
0.93
0.92
0.90
0.87
0.83
0.70
0.75
0.86
0.84
0.79
0.75

0.03
0.02
0.05
0.02
1.10
1.03
1.09
1.04
0.03
0.06
0.02
0.01
1.41
1.38
1.35
1.43

Apparent spatial-temporal differences in daily Ts fluctuations under
different ADI treatments were found, especially in the 5 cm soil layer(Fig. 3). In general, Ts decreased in response to increased irrigation
frequency with groundwater. The lowest Ts occurred in the 3G1S
treatment among different ADI treatments. The average Ts in both years
under 3G1S was 12.2 % lower than under SW and 4.4 % higher thanunder GW. However, the highest Ts occurred in the 1G1S treatment anddecreased only 3.9 % compared with SW. Meanwhile, the average Ts inboth years under 2G1S was 10.7 % lower than under SW. Additionally,
differences in Ts under different ADI treatments decreased withincreasing soil depth. The average Ts at the 20 cm soil depth under
1G1S, 2G1S, 3G1S, GW, and SW was 23.9, 23.0, 21.1, 20.7, and 24.0 ℃,
respectively.
Three consecutive days during the early (DAS 22–24 in 2019 and
DAS 29–31 in 2020), middle (DAS 73–75 in 2019 and 2020), and late
(DAS 106–108 in 2019 and 2020) crop growth stages were chosen tofurther evaluate the impact of ADI on diurnal Ts fluctuations (Fig. 4). In
the early crop growth stage, there was an apparent difference in diurnal
Ts fluctuations under different ADI treatments. The average Ts at 5 cm
soil depth during daytime under 1G1S, 2G1S, and 3G1S was 1.2 %,1.5 %, and 2.0 % lower than under SW, respectively, while it was 1.6 %,1.0 %, and 1.8 % higher than under GW. Moreover, the average Ts at
5 cm soil depth during nighttime under 1G1S, 2G1S, and 3G1S was
1.6 %, 1.6 %, and 1.9 % lower than under SW, respectively, while it was
0.8 %, 1.8 %, and 1.9 % higher than under GW. In the middle crop
growth stage, the maximally diurnal Ts fluctuations occurred in the
3G1S treatment among different ADI treatments (Fig. 4). The average Ts
during daytime and nighttime under 3G1S decreased by 16.6 % and12.3 % compared with SW, respectively. However, 1G1S has thesmallest diurnal Ts fluctuations among different ADI treatments. The
average Ts during daytime and nighttime under 1G1S decreased by
8.5 % and 4.1 % compared with SW, respectively. The diurnal Ts differences under different ADI treatments evidently decreased in the latecrop growth stage. Similar to the daily Ts fluctuations, differences in thediurnal Ts fluctuations under different ADI treatments also decreased in
response to increased soil depth.
3.3. Effect of irrigation event on soil temperature distributions

0.96
0.85
0.87
0.97
0.95
0.96
0.95
0.95
0.84
0.90
0.80
0.73
0.91
0.89
0.89
0.86
0.05
0.07
0.06
0.02
1.45
1.29
1.29
1.29
0.03
0.05
0.03
0.01
1.49
1.56
1.44
1.57

8.12
8.69
7.72
8.94
5.65
4.99
4.86
4.75
9.43
9.16
9.80
9.80
5.21
5.29
4.66
5.38

0.04
0.02
0.05
0.05
1.20
1.29
1.16
1.33
0.02
0.03
0.02
0.01
1.45
1.31
1.37
1.27

RMSE
(cm3 cm− 3/℃)
RMSE
(cm3 cm− 3/℃)
RMSE
(cm3 cm− 3/℃)
EF
(-)
RMSE
(cm3 cm− 3/℃)

One day after irrigation in the early (DAS 31) and middle (DAS 75)
crop growth stages and the harvest day (DAS 140) in 2020 were chosen
to evaluate the effect of irrigation events on soil temperature distributions (Fig. 5).
A noticeable difference in Ts distributions on DAS 31 under different
ADI treatments was found. On DAS 31, the “moderate Ts” (20–22 ℃)
areas under the 1G1S, 2G1S, and SW treatments were 1332.3, 1644.4,
and 1974.2 cm2, respectively. However, there were no “moderate Ts”
areas under the 3G1S and GW treatments. The largest difference in Ts
distributions under different ADI treatments occurred in the middle crop
growth stage (i.e., DAS 75) (Fig. 5). The Ts distribution under 1G1S was

Ts

SWC
2020

Ts

0–20
20–40
40–60
60–100
5
10
15
20
0–20
20–40
40–60
60–100
5
10
15
20
2019

SWC

Soil depth(cm)

calibration period (2019), the range of the RMSE, MRE , and EF values
for SWC under different irrigation treatments were 0.01–0.07 cm3 cm− 3,
6.8 %-8.9 %, and 0.8–1.0, respectively, while they were 0.94–1.45 ◦ C,
4.6 %-5.7 %, and 0.8–1.0 for Ts (Table 4). In the validation period(2020), the range of the RMSE, MRE , and EF values for SWC under
different irrigation treatments were 0.01–0.06 cm3 cm− 3, 7.6 %-9.8 %,and 0.7–0.9, respectively, while they were 1.25–1.57 ◦ C, 4.6 %-5.4 %,and 0.8–0.9 for Ts. In general, the simulation accuracy of the HYDRUS(2D/3D) software was high in terms of SWC and Ts. Additionally, there
was a high consistency between simulated and observed values based on
the result of visual validation (Fig. 2). Therefore, the HYDRUS (2D/3D)
software can be recommended to evaluate the differences in SWC and Ts
under different irrigation treatments.
3.2. Daily and diurnal fluctuations of soil temperature

MRE
(%)

1G1S
SW

MRE
(%)

EF
(-)

2G1S

MRE
(%)

EF
(-)

3G1S

MRE
(%)

GW

MRE
(%)

EF
(-)

Soil & Tillage Research 244 (2024) 106259

Year

Table 4
Statistical performance of HYDRUS (2D/3D) for soil water contents (SWC) in the 0–10, 10–40, 40–60, and 60–100 cm soil layers and soil temperature (Ts) at 5, 10, 15, and 20 cm soil depths during calibration (2019) andvalidation (2020) for the surface water irrigation (SW), one groundwater irrigation after one surface water irrigation (1G1S), two groundwater irrigations after one surface water irrigation (2G1S), three groundwater
irrigations after one surface water irrigation (3G1S), and groundwater water irrigation (GW) treatments. RMSE, MRE , and EF are the root mean square error, the mean relative error, and the modeling efficiency,
respectively.

Y. Zhang et al.

6

Y. Zhang et al.

Soil & Tillage Research 244 (2024) 106259

Fig. 3. The daily fluctuations of the soil temperature at 5, 10, 15, and 20 cm soil depths for the surface water irrigation (SW), one groundwater irrigation after one
surface water irrigation (1G1S), two groundwater irrigations after one surface water irrigation (2G1S), three groundwater irrigations after one surface water irrigation (3G1S), and groundwater water irrigation (GW) treatments during the whole crop growth season in 2019 and 2020.

similar to that under SW, but the Ts distributions under the 2G1S, 3G1S,
and GW treatments were better than under SW, especially under 2G1S.
The “moderate Ts” area under 2G1S reached 1271.8 cm2, representing
an increase of 23.2 % compared with SW. Meanwhile, the “moderate Ts”
area under 2G1S increased by 97.2 %, 96.7 %, and 89.4 % comparedwith 1G1S, 3G1S, and GW, respectively. There was no apparent difference in Ts distributions under different ADI treatments on DAS 140. The“moderate Ts” areas on DAS 140 under the 1G1S, 2G1S, 3G1S, GW, andSW treatments were 20.3, 19.8, 18.9, 17.6, and 17.5 cm2, respectively.
In general, the 2G1S treatment was the most optimal irrigation
treatment.

differences under different ADI treatments occurred. The largest Ts
fluctuations occurred in the 3G1S treatment. The average Tu in both
years under 3G1S increased by 0.53 and 0.23 ℃ compared with 1G1S
and 2G1S, respectively, and the average Td under 3G1S decreased by
0.41 and 0.27 ℃ compared with 1G1S and 2G1S. However, Tu and Td
under 3G1S decreased by 0.44 and 0.16 ℃ compared with GW,
respectively. In the late crop growth stage (the maturation stage), thedifferences in Ts fluctuations under different ADI treatments decreased
compared with the middle crop growth stage. Tu at 5 cm soil depth
under 1G1S, 2G1S, 3G1S, GW, and SW was 1.07, 1.10, 1.14, 1.17, and0.98 ℃, respectively, while it was − 1.15, − 1.22, − 1.24, − 1.25, and− 1.09 ℃ for Td. Additionally, the differences in Ts fluctuations under
different ADI treatments significantly decreased with increased soil
depth.

3.4. Seasonal fluctuations of soil temperature under different irrigation
treatments
Significant seasonal Ts fluctuations under different ADI treatments
were found. There were no apparent differences in seasonal Ts fluctuations under different ADI treatments in the early crop growth stage (theseeding and elongation stages). Tu at 5 cm soil depth under 3G1S only
increased by 0.24 and 0.09 ℃ compared with 2G1S and 1G1S, respectively, while Td at 5 cm soil depth under 3G1S decreased by 0.10 and
0.02 ℃ compared with 2G1S and 1G1S (Table 5). However, in the
middle crop growth stage (the tasseling and filling stages), the largest Ts

3.5. Soil water balances, crop yields, and water use efficiency for different
irrigation treatments
There are apparent differences in soil water balance components indifferent irrigation treatments (Table 6). In general, evapotranspiration(ET) in different ADI treatments was lower than in SW and higher than in
GW. The highest ET in different ADI treatments occurred in the 1G1S,
which increased by 0.4 % and 3.1 % compared with 2G1S and 3G1S,
7

Y. Zhang et al.

Soil & Tillage Research 244 (2024) 106259

Fig. 4. The diurnal fluctuations of the soil temperature at 5, 10, 15, and 20 cm soil depths under the surface water irrigation (SW), one groundwater after one surface
water irrigation (1G1S), two groundwater after one surface water irrigation (2G1S), three groundwater after one surface water irrigation (3G1S), and groundwater
water irrigation (GW) treatments during the early (DAS 22–24 in 2019; DAS 29–31 in 2020), middle (DAS 73–75 in 2019 and 2020), and late crop growth stages(DAS 106–108 in 2019 and 2020) in 2019 and 2020.

Fig. 5. Soil temperature distributions in the 0–40 cm soil layer one day after irrigation during the early (DAS 31) and middle (DAS 75) crop growth stages and theharvest day (DAS 140) in 2020 under the surface water irrigation (SW), one groundwater irrigation after one surface water irrigation (1G1S), two groundwater
irrigations after one surface water irrigation (2G1S), three groundwater irrigations after one surface water irrigation (3G1S), and groundwater water irrigation(GW) treatments.

respectively. Average ET in both years in different ADI treatments
decreased by 2.2 % compared to SW, while it increased by 2.4 %compared to GW. However, capillary water rise (R) decreased in
response to increased irrigation frequency. The highest R occurred in the
1G1S treatment among different ADI treatments. Average R in both
years under 1G1S was 4.3 %, 7.2 %, and 10.3 % higher than under
2G1S, 3G1S, and GW, respectively, while it was 2.3 % lower than underSW. However, there was an opposite trend in the soil water storagechange in different treatments. SWS in different ADI treatments washigher than in SW and lower than in GW. 3G1S had the lowest SWS,
12.2 % and 5.9 % lower than 1G1S and 2G1S, respectively. Meanwhile,
average crop yield (Y) in both years under 1G1S increased by 3.6 % and
8.8 % compared to 3G1S and GW, respectively, and decreased by 2.9 %compared to SW. However, there was no apparent difference in crop
yield between 1G1S and 2G1S. The highest water use efficiency (WUE)of 24.9 kg ha− 1 mm− 1 occurred in the 2G1S treatment among different

ADI treatments.
4. Discussion
4.1. The effects of different irrigation water sources on soil temperatures
Apparent differences in Ts fluctuations under different irrigationwater sources were found. Generally, Ts fluctuations under GW irrigation were higher than under SW. For instance, the daily minimum Ts
under GW decreased by 27.9 % compared with SW (Fig. 3 and Fig. 4).
This is caused by root water uptake under SW being higher than under
GW, which resulted in decreased root zone pressure head (hroot) under
SW. The average hroot in both years under SW was 14.5 % lower thanunder GW (Fig. 6). Meanwhile, Ts increased in response to a decrease in
SWC under SW. The main reason was a decreased soil thermal conductivity for reduced hroot (SWCs) (Obia et al., 2020; Nyawade et al.,8

− 1.64− 1.35− 1.07− 0.64− 1.7− 1.56− 1.18− 0.73− 2.25− 1.66− 1.39− 1.07− 2.38− 1.4− 0.98− 0.59− 1.41− 0.91− 0.8− 0.54− 1.63− 1.33− 1.05− 0.62− 1.51− 1.18− 1.11− 0.69− 2.14− 1.53− 1.34− 0.8− 2.07− 1.25− 0.93− 0.56− 1.39− 0.89− 0.79− 0.54

3G1S
2G1S− 1.63− 1.33− 1.05− 0.62− 1.45− 1.18− 1.1− 0.54− 1.76− 1.44− 1.22− 0.72− 1.73− 1.2− 0.88− 0.52− 1.38− 0.89− 0.7− 0.53− 1.52− 1.05− 0.94− 0.6− 1.33− 1.13− 1.01− 0.58− 1.72− 1.38− 1.03− 0.62− 1.64− 1.1− 0.81− 0.5− 1.29− 0.85− 0.7− 0.491G1S
SW− 1.48− 0.79− 0.68− 0.57− 1.04− 1.06− 0.88− 0.51− 1.26− 1.32− 0.99− 0.44− 1.55− 1.03− 0.7− 0.44− 1.25− 0.84− 0.64− 0.49
1.92
1.13
0.89
0.54
2.19
1.7
1.16
0.73
2.73
1.76
1.26
0.57
2.04
1.32
1.25
0.63
1.25
1.28
0.93
0.54
1.92
1.13
0.89
0.54
1.97
1.44
0.98
0.73
1.99
1.52
1.15
0.55
1.74
1.29
1.17
0.62
1.23
1.2
0.86
0.44
1.92
1.13
0.89
0.54
1.68
1.11
0.98
0.63
1.91
1.34
0.93
0.49
1.63
1.04
0.77
0.58
1.18
1.17
0.86
0.41
1.59
1.07
0.82
0.54
1.49
1.07
0.92
0.62
1.38
1.27
0.86
0.44
1.44
1.02
0.76
0.48
1.15
1.15
0.82
0.37
1.11
0.86
0.72
0.54
1.45
1.05
0.88
0.5
1.26
0.77
0.71
0.41
1.09
0.75
0.74
0.29
1.02
0.96
0.81
0.36
− 3.39− 3.22− 2.65− 2.2− 2.91− 2.16− 1.75− 0.98− 2.17− 1.68− 1.17− 0.92− 2.74− 1.76− 1.56− 0.95− 1.08− 0.88− 0.83− 0.56− 3.37− 3.18− 2.65− 2.2− 1.97− 2.02− 1.49− 0.9− 2.05− 1.49− 1.14− 0.89− 2.64− 1.67− 1.48− 0.8− 1.08− 0.88− 0.61− 0.43− 3.39− 3.18− 2.5− 2.2− 1.95− 1.64− 1.47− 0.82− 1.93− 1.46− 1.04− 0.83− 2.42− 1.57− 1.39− 0.76− 1.05− 0.82− 0.61− 0.42− 3.28− 2.69− 2.4− 2.18− 1.94− 1.6− 1.39− 0.79− 1.7− 1.43− 1.01− 0.67− 2.22− 1.56− 1.32− 0.72− 1.01− 0.82− 0.61− 0.42

There is an apparent difference in Ts variations under different ADI
treatments. In general, Ts decreased in response to increased irrigation
frequency (Kincaid et al., 1993; Seo, 1989; Wang et al., 2008). For
example, this study found that the average Ts under 3G1S was 0.6 and
0.4 ℃ lower than under 1G1S and 2G1S, respectively. Sun et al. (2017)
similarly reported that Ts under high irrigation frequency (17 irrigation
times) of groundwater significantly decreased by 1.5 ℃ compared with
low irrigation frequency (9 irrigation times). The main reason could be
attributed to the temperature difference between soil water and soil
solids. The temperature of infiltrating soil water was lower than that ofsoil, which caused heat to be transferred from the soil solids to infiltrating soil water along the energy gradient after irrigation. Higher
SWCs prolonged the duration of increased Ts (Tomoyuki et al., 1984;
Mitsuno et al., 1983; Votrubova et al., 2012). Thus, low-temperaturestress may restrain crop growth under high groundwater irrigation

3.85
2.81
2.59
1.50
3.19
2.6
2.38
1.54
2.8
1.87
1.6
0.64
3.37
2.52
1.65
0.72
1.08
0.62
0.41
0.50
3.76
2.80
2.58
1.49
2.86
2.52
2.32
1.51
2.23
1.63
1.51
0.63
3.24
2.52
1.60
0.72
1.04
0.62
0.41
0.42
3.70
2.46
2.35
1.48
2.79
2.2
1.92
1.45
1.44
1.16
1.14
0.59
2.83
2.18
1.40
0.69
0.98
0.45
0.41
0.4

3.74
2.51
2.43
1.48
2.81
2.38
2.02
1.50
1.62
1.33
1.23
0.61
3.13
2.32
1.56
0.69
1.02
0.62
0.41
0.42

4.2. Soil temperature variations under different alternate irrigations withsurface water and groundwater

Elongation stage(DAS 36–71 in
2019; DAS 31–68
in 2020)
Tasseling stage(DAS
72–115 in 2019;
DAS 69–107 in
2020)
Filling stage(DAS
116–127 in 2019;
DAS 108–120 in2020)
Maturation stage(DAS 127–140 in
2019; DAS
121–140 in 2020)

5
10
15
20
5
10
15
20
5
10
15
20
5
10
15
20
5
10
15
20
Seedling stage(DAS
0–35 in 2019; DAS
0–30 in 2020)

3.69
2.45
2.28
1.47
2.67
1.71
1.65
1.37
1.38
1.05
1.02
0.56
2.5
2.07
1.36
0.65
0.94
0.44
0.38
0.4

Td (℃)

GW
3G1S
2G1S
1G1S
Tu (℃)

SW
SW
GW
3G1S
2G1S
SW

Tu (℃)1G1S

Td (℃)

1G1S

2G1S

3G1S

GW

2020

Soil
depths(cm)20192018). Thus, Ts fluctuations under SW were lower than under GW.
Meanwhile, ET in SW was higher than in the ADI and GW treatments
due to higher soil evaporation caused by higher soil temperatures.
Moreover, SWS decreased in response to an increase of ET in SW, while R
increased with an increase in the soil water potential gradient between
surface soil water and groundwater (Table 6). Additionally, the differences in Ts fluctuations under irrigations with different water sources
mainly occurred in the 5 cm soil depth during the middle crop growthstage. There were two main reasons for that. On the one hand, Ts fluctuations increased due to a higher irrigation frequency to meet cropwater requirements and avoid the effects of soil heat stress on root water
uptake during the middle crop growth stage. In general, Ts decreased
after irrigation due to increasing soil thermal conductivity and soil
thermal capacity. However, Ts rapidly increased once hroot decreased in
response to increased crop water requirement (i.e., root water uptake)under high air temperature conditions, particularly in the SW treatment.
For instance, the average hroot in both years under SW during the middle
crop growth stage decreased by 2.7 %, 7.4 %, 11.5 %, and 14.4 %compared to 1G1S, 2G1S, 3G1S, and GW, respectively. On the other
hand, Ts in the 5 cm soil depth was affected by irrigation since thewetting front under drip irrigation was mainly in the 0–10 cm soil layer
due to small drip flow. Similar locations of the wetting front have beenreported in the literature (Chen et al., 2021; Sampathkumar et al., 2012;
Surendran and Chandran, 2022). Additionally, Ts fluctuations during
the entire crop growth season in 2019 were higher than in 2020 due tolower SWCs (Fig. 3 and Fig. 6). The maximum and minimum hroot in2019 decreased by 32.4 % and 48.1 % compared with 2020, respectively (Fig. 6).
The average Ts under GW was 2.0 ℃ lower than under SW. The
reason was mainly attributed to the temperature difference between
groundwater and the soil, which was higher than the difference between
surface water and the soil. Ts decreased after groundwater irrigation
since the heat was transferred from soil particles (with higher temperature) to infiltrating soil water (with lower temperature) based on theenergy balance law. Similar results were found in the literature(Wierenga et al., 1970; Shi et al., 2023). For example, Wierenga et al.
(1970) reported that the difference in surface Ts (in a 1 cm soil depth)
between two plots irrigated with cool (4.1 ℃) and warm (21.6 ℃) waterreached 3.0 ℃ during nighttime. Shi et al. (2023) demonstrated that Ts
under GW was 0.2–1.4 ℃ lower than under SW. Additionally, areasonable irrigation strategy could effectively eliminate the adverse
effects of heat stress on crop growth (Becker et al., 2023; Birthal et al.,2021; Meena et al., 2015). For example, this study found that the“moderate Ts” area under GW increased by 96.3 % compared with SW
on DAS 75 (Fig. 5). Therefore, groundwater irrigation should beconsidered during high air temperature periods to overcome the heatand drought stresses.

− 2.98− 2.66− 1.95− 1.66− 1.85− 1.49− 1.22− 0.7− 1.08− 1.38− 0.94− 0.64− 2.02− 1.52− 1.32− 0.72− 0.92− 0.78− 0.61− 0.42GW

Soil & Tillage Research 244 (2024) 106259

Growth stage

Table 5The average increase (Tu) and decrease (Td) in daily temperatures at 5, 10, 15, and 20 cm soil depths under the surface water irrigation (SW), one groundwater irrigation after one surface water irrigation (1G1S), twogroundwater irrigations after one surface water irrigation (2G1S), three groundwater irrigations after one surface water irrigation (3G1S), and groundwater water irrigation (GW) treatments during different crop growthstages in 2019 and 2020. DAS - days after seedling.

Y. Zhang et al.

9

Y. Zhang et al.

Soil & Tillage Research 244 (2024) 106259

Table 6
Soil water balance components (i.e., I, P, ET, R, and SWS), crop yields (Y), and water use efficiency (WUE) for different irrigation treatments (i.e., surface water
irrigation (SW), one groundwater irrigation after one surface water irrigation (1G1S), two groundwater irrigations after one surface water irrigation (2G1S), threegroundwater irrigations after one surface water irrigation (3G1S), and groundwater water irrigation (GW) in 2019 and 2020. Note that I, P, ET, R, and SWS represent
irrigation, rainfall, evapotranspiration, capillary rise, and the change in soil water storage, respectively.
Years
2019

2020

Treatments
SW
1G1S
2G1S
3G1S
GW
SW
1G1S
2G1S
3G1S
GW

Soil water balances (mm)
I

P

ET

R

SWS

270
270
270
270
270
270
270
270
270
270

87.9
87.9
87.9
87.9
87.9
144.7
144.7
144.7
144.7
144.7

505.6
498.1
492.4
486.3
477.8
598.0
595.9
585.0
579.2
575.8

100.1
96.3
93.1
90.0
85.4
120.9
119.6
113.6
110.3
108.4

− 47.6− 43.9− 41.4− 38.4− 34.5− 62.4− 61.6− 56.7− 54.2− 52.7Y (t ha− 1)

WUE (kg ha− 1 mm− 1)

13.7b
13.3c
13.3c
13.0d
12.2e
14.0a
13.5bc13.4bc
13.0d
12.3e

27.0a
26.8ab
27.0a
26.6b
25.6c
23.4d
22.7de22.8de22.4e21.3 fFig. 6. Time series of root zone pressure heads for different irrigation treatments (i.e., surface water irrigation (SW), one groundwater irrigation after one surface
water irrigation (1G1S), two groundwater irrigations after one surface water irrigation (2G1S), three groundwater irrigations after one surface water irrigation(3G1S), and groundwater water irrigation (GW)) in 2019 (left) and 2020 (right).

Fig. 7. Simulated soil temperature variations at 5, 10, 15, and 20 cm soil depths in 2019 and 2020 under surface water irrigation (SW), one groundwater irrigation
after one surface water irrigation (1G1S), two groundwater irrigations after one surface water irrigation (2G1S), three groundwater irrigations after one surface water
irrigation (3G1S), four groundwater irrigations after one surface water irrigation (4G1S), five groundwater irrigations after one surface water irrigation (5G1S), andgroundwater water irrigation (GW) treatments.

frequency. For instance, the “relatively low Ts” area under 3G1S
increased by 7.2 % and 4.4 % compared with the 1G1S and 2G1S
treatments. Wanjura et al. (1996) similarly demonstrated that higher
irrigation frequency quickly resulted in low-temperature stress despite
effectively avoiding water stress in the rapid crop growth stage. Additionally, although 1G1S effectively decreased the “relatively low Ts”
area compared with 3G1S, the “relatively high Ts” area under 1G1S
noticably increased. Crop growth also may be restrained byhigh-temperature stress under 1G1S. Overall, the 2G1S treatment was
the most optimal irrigation treatment among different ADI treatments,
with a “moderate Ts” area of 1271.8 cm2. Therefore, 2G1S could berecommended as an optimal irrigation strategy in this region.
Five different irrigation frequency treatments (i.e., 1G1S, 2G1S,
3G1S, 4G1S, and 5G1S) were further simulated to determine theresponse characteristic of Ts to irrigation frequency for the same

irrigation depth (22.5 mm) (Fig. 7). Additionally, five different irrigation depth treatments, i.e., 15, 22.5, 30, 37.5, and 45 mm, werecompared for the same irrigation frequency (2G1S). The simulated results for these scenarios showed that the sensitivity of Ts to the irrigation
depth variation was higher than that to the irrigation frequency. Ts only
decreased by 0.9 ℃ for the same irrigation depth and increased irrigation frequency by one. However, Ts decreased by 1.2 ℃ when the irrigation depth increased by 7.5 mm. The longest duration of the“moderate Ts” (88 d) occurred for the 22.5 mm irrigation depth under
2G1S. Therefore, 22.5 mm irrigation depth and the 2G1S treatment
could be recommended as an optimal irrigation strategy.
5. Conclusions
The HYDRUS (2D/3D) software could precisely simulate soil water
10

Soil & Tillage Research 244 (2024) 106259

Y. Zhang et al.

content and heat dynamics under different irrigation treatments, with
the RMSE, M RE , and EF values of 0.01–0.06 cm cm− 3, 7.6 %-9.8 %,and 0.7–0.9 for SWC, respectively, and 1.25–1.57 ◦ C, 4.6 %-5.4 %, and
0.8–0.9 for Ts in the verification period. Significant spatial-temporal
differences in daily Ts fluctuations under different ADI treatments
were found, especially at the 5 cm soil depth. In general, Ts decreased in
response to an increase in the irrigation frequency of groundwater. The
lowest Ts occurred in the 3G1W treatment among different ADI treatments. The average Ts in both years under 3G1W was 12.2 % lower thanunder SW and 4.4 % higher than under GW. However, the highest Ts
occurred in the 1G1S treatment and decreased by only 3.9 % comparedwith SW. The largest “moderate Ts” area of 1271.8 cm2 occurred in the
2G1S treatment. Meanwhile, the longest duration of the “moderate Ts”
occurred for the 22.5 mm irrigation depth under 2G1S. Therefore, theirrigation depth of 22.5 mm and the 2G1S treatment could be recommended as the optimal irrigation strategy in this region.

Becker, R., Schüth, C., Merz, R., Khaliq, T., Usman, M., Beek, T.A.D., Kumar, R.,
Schulz, S., 2023. Increased heat stress reduces future yields of three major crops inPakistan’s Punjab region despite intensification of irrigation. Agric. Water Manag.
281, 108243.
Birthal, P.S., Hazrana, J., Negi, D.S., Pandey, G., 2021. Benefits of irrigation against heatstress in agriculture: evidence from wheat crop in India. Agric. Water Manag. 255,
106950.
Bristow, K.L., Campbell, G.S., Papendick, R.I., Elliott, L.F., 1986. Simulation of heat andmoisture transfer through a surface residue-soil system. Agric. . Meteorol. 36,
193–214.
Carsel, R.F., Parrish, R.S., 1988. Developing joint probability distributions of soil water
retention characteristics. Water Resour. Res. 24, 755–769.
Chen, N., Li, X.Y., Šimůnek, J., Shi, H.B., Hu, Q., Zhang, Y.H., 2021. Evaluating theeffects of biodegradable and plastic film mulching on soil temperature in a dripirrigated field. Soil Tillage Res. 213, 105116.
Chung, S.O., Horton, R., 1987. Soil heat and water flow with a partial surface mulch.
Water Resour. Res. 23, 2175–2186.
Daniel, D.F., Dallacort, R., Barbieri, J.D., Modolo, A.J., Carvalho, M.A.C.D., Tieppo, R.C.,
Yamashita, O.M., 2022. Crotalaria sowing times intercropped with off-season maizein the variability of soil temperature and moisture. J. Agric. Stud. 10 (3), 126–151.
Deb, P., Moradkhani, H., Han, X.L., Abbaszadeh, P., Xu, L., 2022. Assessing irrigation
mitigating drought impacts on crop yields with an integrated modeling framework.
J. Hydrol. 609, 127760.
Ebrahimi, E., Tekeste, M.Z., Huth, N.I., Antille, D.L., Archontoulis, S.V., Horton, R., 2022.
Measured and modeled maize and so water use on pipeline disturbed land. Soil
Tillage Res 220, 105340.
Filipović, V., Romić, D., Romić, M., Borošić, J., Filipović, L., Mallmann, F.J.K.,
Robinson, D.A., 2016. Plastic mulch and nitrogen fertigation in growing vegetables
modify soil temperature, water and nitrate dynamics: experimental results and amodeling study. Agric. Water Manag. 176, 100–110.
Garg, N., Cheema, D.S., 2011. Assessment of fruit quality attributes of tomato hybrids
involving ripening mutants under high temperature conditions. Sci. Hortic. 131,
29–38.
van Genuchten, M.T., 1980. A closed-form equation for predicting the hydraulicconductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44, 892–898.
Gong, X.W., Li, X.M., Li, Y., Bo, G.K., Qiu, R.J., Huang, Z.D., Gao, S.K., Wang, S.S., 2023.
An improved model to simulate soil water and heat: A case study for drip-irrigatedtomato grown in a greenhouse. Agric. Water Manag. 277, 108121.
Ibrahim, M., Khan, A., Anjum Ali, W., Akbar, H., 2020. Mulching techniques: An
approach for offsetting soil moisture deficit and enhancing manure mineralization
during maize cultivation. Soil Tillage Res 200, 104631.
Ierna, A., Mauromicale, G., 2022. How irrigation water saving strategy can affect tuber
growth and nutritional composition of potato. Sci. Hortic. 299, 111034.
Jha, S.K., Gao, Y., Liu, H., Huang, Z.D., Wang, G.S., Liang, Y.P., Duan, A.W., 2017. Rootdevelopment and water uptake in winter wheat under different irrigation methods
and scheduling for north China. Agric. Water Manag. 182, 139–150.
Karandish, F., Shahnazari, A., 2016. Soil temperature and maize nitrogen uptake
improvement under partial root-zone drying irrigation. Pedosphere 26 (6), 872–886.
Karandish, F., Šimůnek, J., 2019. A comparison of the HYDRUS (2D/3D) and SALTMED
models to investigate the influence of various water-saving irrigation strategies on
the maize water footprint. Agric. Water Manag. 213, 809–820.
Kassaye, K.T., Boulange, J., Tu, L.H., Saito, H., Watanabe, H., 2021. Soil water content
and soil temperature modeling in a vadose zone of Andosol under temperatemonsoon climate. Geoderma 384, 114797.
Kincaid, D.C., Westermann, D.T., Trout, T.J., 1993. Irrigation and soil temperatureeffects on Russet Burbank quality. Am. Potato J. 70 (10), 711–723.
Kishore, G., Babu, B.M., Mattaparti, L.D., 2022. Influence of plastic mulching andirrigation levels on soil temperature, moisture and water use efficiency of tomato
crop (Solanum lycopersicum). Int. J. Plant Soil Sci. 34 (20), 277–282.
Kodešová, R., Fér, M., Klement, A., Nikodem, A., Teplá, D., Neuberger, P., Bureš, P.,
2014. Impact of various surface covers on water and thermal regime of Technosol.
J. Hydrol. 519, 2272–2288.
Lamptey, S., Li, L., Xie, J., Coulter, J.A., 2019. Tillage system affects soil water andphotosynthesis of plastic-mulched maize on the semiarid loess plateau of China. Soil
Tillage Res 196, 104479.
Lee, O.N., Park, H.Y., 2020. Effects of different colored film mulches on the growth andbolting time of radish (Raphanus sativus l.). Sci. Hortic. 266, 109271.
Liakatas, A., Clark, J.A., Monteith, J.L., 1986. Measurements of the heat balance underplastic mulches. Part I. Radiation balance and soil heat flux. Agric. . Meteorol. 36 (3),
227–239.
Liao, X.L., Su, Z.H., Liu, G.D., Zotarelli, L., Cui, Y.Q., Snodgrass, C., 2016. Impact of soil
moisture and temperature on potato production using seepage and center pivot
irrigation. Agric. Water Manag. 165, 230–236.
Matsue, Y., Takasaki, K., Abe, J., 2021. Water management for improvement of rice
yield: appearance quality and palatability with high temperature during ripening
period. Rice Sci. 28 (4), 409–416.
Meena, J.R., Behera, U.K., Chakraborty, D., Sharma, A.R., 2015. Tillage and residue
management effect on soil properties, crop performance and energy relations ingreen gram (Vigna radiata L.) under maize-based cropping systems. Int. Soil Water
Conserv. Res 3 (4), 261–272.
Mehdizadeh, S., Fathian, F., Safari, M.J.S., Khosravi, A., 2020. Developing novel hybrid
models for estimation of daily soil temperature at various depths. Soil Tillage Res
197, 104513.
Menossi, M., Cisneros, M., Alvarez, V.A., Casalongue, C., 2021. Current and emerging
biodegradable mulch films based on polysaccharide bio-composites. A review.
Agron. Sustain. Dev. 41 (53), 1–27.

CRediT authorship contribution statement
Jiří Šimůnek: Writing – review & editing, Software, Methodology.
Ning Chen: Writing – review & editing, Writing – original draft,
Methodology, Investigation. Yuehong Zhang: Writing – review &editing, Writing – original draft, Visualization, Methodology, Investigation, Formal analysis, Data curation. Xianyue Li: Writing – review &editing, Validation, Methodology, Funding acquisition. Qi Hu: Investigation. Haibin Shi: Writing – review & editing, Validation,
Conceptualization.
Declaration of Competing Interest
The authors declare that they have no known competing financial
interests or personal relationships that could have appeared to influence
the work reported in this paper: No conflict of interest exits in thesubmission of this manuscript, and manuscript is approved by all authors for publication. I would like to declare on behalf of my co-authorsthat the work described was original research that has not been published previously, and not under consideration for publication elsewhere, in whole or in part. All the authors listed have approved themanuscript that is enclosed.
Data Availability
Data will be made available on request.
Acknowledgments
This research was jointly supported by the National Key Researchand Development Program of China (2021YFC3201202), the Inner
Mongolia Natural Science Foundation (2024QN05039), the Scientificand Technological Transformative Project of Inner Mongolia(2021EEDSCXSFQZD011, 2021CG0022), the Inner Mongolia Scienceand Technology Program (2022YFHH0039), the Scientific Research
Projects of High-Level and Excellent Doctoral Talents for Inner Mongolia
Agricultural University (NDYB2023-20), the First-Class Academic Subjects Special Research Project of the Education Department of Inner
Mongolia Autonomous Region (YLXKZX-NND-022 and YLXKZX-NND023), and the Inner Mongolia Autonomous Region’s Talent Introduction
and Research Support Project.
References
Allen, R.G., Pereira, L.S., Raes, D., Smith, M., 1998. Crop evapotranspiration: guidelines
for computing crop water requirements, irrigation and drainage Paper 56. United
Nations FAO, Rome, p. 300.
Amato, M.T., Giménez, D., 2024. Predicting monthly near-surface soil temperature from
air temperature and the leaf area index. Agric. . Meteorol. 345, 109838.

11

Y. Zhang et al.

Soil & Tillage Research 244 (2024) 106259
Šimůnek, J., van Genuchten, M.T., Šejna, M., 2016. Recent developments andapplications of the HYDRUS computer software packages. Vadose Zone J. 15, 25.
Šimůnek, J., Brunetti, G., Jacques, D., van Genuchten, M.Th, Šejna, M., 2024. Recent
developments and applications of the HYDRUS computer software packages since2016. Vadose Zone J. 23 (1), 29 e20310.
Skaggs, T.H., Trout, T.J., Šimůnek, J., Shouse, P.J., 2004. Comparison of HYDRUS-2D
simulations of drip irrigation with experimental observations. J. Irrig. Drain. Div.
Am. Soc. Civ. Eng. 130 (4), 304–310.
Soureshjani, H.K., Dehkordi, A.G., Bahador, M., 2019. Temperature effect on yield ofwinter and spring irrigated crops. Agric. . Meteorol. 279, 107664.
Sun, G.F., Qu, Z.Y., Du, B., Ren, Z.S., Liu, A.Q., 2017. Water heat salt effects of mulcheddrip irrigation maize with different irrigation scheduling in Hetao Irrigation District.
Trans. Chin. Soc. Agric. Eng. 33 (12), 144–152.
Surendran, U., Chandran, K.M., 2022. Development and evaluation of drip irrigation andfertigation scheduling to improve water productivity and sustainable crop
production using HYDRUS. Agric. Water Manag. 269, 107668.
Tomoyuki, Ishida, Tooru, Mitsuno, Toshisuke, Maruyama, 1984. Heat conduction of soil
under relatively higher water content condition. Trans. Jpn. Soc. Irrig., Drain. Rec.
Eng. 1984 (110), 67–74.
Votrubova, J., Dohnal, M., Vogel, T., Tesa, M., 2012. On parameterization of heat
conduction in coupled soil water and heat flow modelling. Soil Water Res 7 (4),
125–137.
de Vries, D.A., 1963. The thermal properties of soils, In Physics of Plant Environment,edited by R. W. van Wijk, pp. 210-235, North Holland, Amsterdam.
Wang, J.D., Gong, S.H., Yu, Y.D., Liang, Y.P., 2008. Effect of surface irrigation frequency
on soil water and temperature as well as growth of spring maize. J. Hydraul. Eng. 39(4) 500–505.
Wanjura, D.F., Mahan, J.R., Upchurch, D.R., 1996. Irrigation starting time effects on
cotton under high-frequency irrigation. Agron. J. 88 (4), 561.
Wechsung, F., Ritter, M., Wall, G.W., 2021. The upper homeostatic range for thetemperature–yield response of irrigated US wheat down revised from a theoreticaland experimental perspective. Agric. . Meteorol. 307, 108478.
Wierenga, P.J., Hagan, R.M., Nielsen, D.R., 1970. Soil temperature profiles during
infiltration and redistribution of cool and warm irrigation water. Water Resour. Res.
6 (1), 230–238.

Mitsuno, T., Ishida, T., Maruyama, T., 1983. Model for soil heat conduction under athree-phase condition and calculation of soil heat conductivity. Trans. Jpn. Soc.
Irrig. Drain. Rec. Eng. 1983, 35–43.
Moreno, M.M., Moreno, A., 2008. Effect of different biodegradable and polyethylene
mulches on soil properties and production in a tomato crop. Sci. Hortic. 116 (3),256–263.
Nyawade, S.O., Karanja, N.N., Gachene, C.K.K., Geldermann, E.S., Parker, M., 2018.
Effect of potato hilling on soil temperature, soil moisture distribution and sedimentyield on a sloping terrain. Soil Tillage Res 184, 24–36.
Obia, A., Cornelissen, G., Martinsen, V., Smebye, A.B., Mulder, J., 2020. Conservation
tillage and biochar improve soil water content and moderate soil temperature in a
tropical Acrisol. Soil Tillage Res 197, 104521.
Rai, V., Pramanik, P., Das, T.K., Aggarwal, P., Bhattacharyya, R., Krishnan, P., Sehgal, V.
K., 2019. Modelling soil hydrothermal regimes in pigeon pea under conservation
agriculture using Hydrus-2D. Soil Tillage Res 190, 92–108.
Rosa, R.D., Paredes, P., Rodrigues, G.C., Alves, I., Fernando, R.M., Pereira, L.S., Allen, R.
G., 2012. Implementing the dual crop coefficient approach in interactive software. 1.
Background and computational strategy. Agric. Water Manag. 103, 8–24.
Sahai, H., Valverde, M.G., Morales, M.M., Hernando, M.D., Real, A.M.A.D., FernándezAlba, A.R., 2023. Exploring sorption of pesticides and PAHs in microplastics derived
from plastic mulch films used in modern agriculture. Chemosphere 333, 138959.
Sampathkumar, T., Pandian, B.J., Mahimairaja, S., 2012. Soil moisture distribution androot characters as influenced by deficit irrigation through drip system incotton–maize cropping sequence. Agric. Water Manag. 103, 43–53.
Scarascia-Mugnozza, G., Sica, C., Russo, G., 2011. Plastic materials in European
agriculture: actual use and perspectives. J. Agric. Eng. 42 (3).
Schaap, M., Leij, F., van Genuchten, M.Th, 2001. Rosetta: A computer program forestimating soil hydraulic parameters with hierarchical pedotransfer functions.
J. Hydrol. 251, 163–176.
Seo, S., 1989. Studies on the pasture irrigation. 5; Effect of irrigation and nitrogen
fertilization during mid-summer season on growth, yield and nitrate nitrogenconcentration of grasses, and botanical composition in pasture. Korean J. Anim. Sci.
11, 569–576.
Shi, W.J., Zhang, X.J., Yang, J., Impa, S.M., Wang, D., Lai, Y.S., Yang, Z.J., Xu, H., Wu, J.
S., Zhang, J.H., Jagadish, S.V.K., 2023. Irrigating with cooler water does not reverse
high temperature impact on grain yield and quality in hybrid rice. Crop J. 11 (3),904–913.

12