Two-dimensional numerical simulations of soil-water and heat flow in a rainfed soybean field under plastic mulching

Numerical simulation can help understanding waterand heat-flow systems through plastic-mulched soils. An effective simulation approach is crucial to know the role of plastic mulch in a soil ecosystem, which can save water in agriculture. A field experiment was conducted at Gifu University in a rainfed soybean cultivation under plastic mulch and bare soil treatments to clarify the soil water and heat flow mechanism. Furthermore, the two-dimensional numerical software HYDRSUS-2D model with different boundaries at the soil surface was used to simulate water and heat flows. Firstly, soil hydraulic parameters were estimated by inverse solution using laboratory-measured data and then coupled soil-water and heat flows were simulated by optimizing soil thermal parameters by inverse solution. The HYDRUS-2D model simulated water and heat flow through the root zone depths satisfactorily. The root-mean square error (RMSE) was 0.015–0.030, and 0.046–0.055 cm cm 3 for the plastic mulch, and bare soil, respectively, in estimating soil moisture and 0.66–1.28, and 0.70–1.54 C, respectively in estimating soil temperature. Water infiltration was 61% lower in the plastic-mulched soil, which reduced soil evaporation as well as soil-moisture storage changes compared to bare soil. This study can be applied to design and manage different plastic mulching patterns in rain-fed crop cultivation.


INTRODUCTION
Global warming and irregular rainfall patterns have been limiting water resources in arid and semi-arid regions, which need effective utilization of water for enhancing crop production (Bradford et al. ). Plastic film mulching, considered as water-saving technology, is a well-known farming practice in dryland areas (Qin et al. ) for regulating soil temperature and conserving soil moisture by reducing soil evaporation (Kader et al. a). So far, China has become the largest user of plastic film mulching (Daryanto et al. ), which is mainly popular in the Loess Plateau area where frequent water shortage occurs due to splash rainfall (Wang et al. ; He et al. ). Plastic mulch improves soil-moisture retention and availability, and exhibits soilwarming effects by regulating radiation and thermal properties (Tian et al. ; Kader et al. ). The use of plastic film increases the resistance to water vapour transport from the soil to atmosphere; this resistance is the major controller of evaporation (Davarzani et al. ). In addition, plastic mulch, by reducing soil evaporation and altering soil temperature in crop fields, has demonstrated great potential for improving water-use efficiency as well as crop yield (Li et  In all these application methods, the farm land is covered with plastic film, which limits the entry of rainwater and air into the soil surface. Only a little amount of water can enter into the soil through planting holes and infiltration at the edge of the mulched-bed (film-side infiltration) (Chen et al. ). Therefore, new application patterns of plastic mulching need to be developed in humid subtropical climates like the central Japan area for increasing water use efficiency in order to ensure water sustainability in irrigated farming.
Knowledge of soil thermal and optical properties is essential when assessing heat transport in mulching soil.
Planting holes in plastic mulching support soil aeration for plant emergence and rainwater distribution into the soil that would increase water availability in the root zone. It is important to manage the critical parameters of plastic film (e.g. diameter of the hole and spacing of the holes) for effective conservation of water with mulches (Li et al. ). The percent open area on the plastic film may have a significant impact on soil evaporation and infiltration into soil; a greater hole ratio provides a smaller vapour transfer resistance so that water vapour can easily cross through plastic mulch to the atmosphere (Qi et al. ).
Until now, soil-water and heat movement and their distributions by plastic mulch is not completely understood due to close interactions of microclimate, soil environment and plant growth (Yang et al. ; Steinmetz et al. ). Therefore, quantification of the effects of plastic mulching by numerical simulations of water and heat flows along with root-water uptake may increase the effectiveness and popularity of the system. Numerical modelling is an effective way to explore water-and heat-flow distributions into soil ecosystems (Qi et al. ). HYDRUS is a popular software that can simulate water, heat and nutrient dynamics in variably-saturated porous media (Šimunek et al. ). The HYDRUS-2D version has superiority over its 1D  both water and heat flow for the plastic-covered area, assuming that film mulch protects water flow as well as heat flow into the soil. But, those studies did not quantify the effects of water fluxes on different boundaries and the soil-water balance. In the case of coupled flow of water and heat, the noflux heat-flow boundary is inappropriate for plastic mulching since heat flow to the soil profile occurs through the planting holes. In actual conditions, the plastic covers prevent water flow but can influence (increase or reduce) heat flow. It is therefore important to select the correct boundary conditions in HYDRUS-2D for plastic mulching. A variable flux boundary in a mulching area where water flux is to be zero and heat flux is to be governed by air and soil surface temperatures may be the correct boundary for HYDRUS-2D simulation with plastic mulching. In this study, the crop planting holes and of the plastic mulching have been taken into the modeling approach by considering the atmospheric boundary, and the plastic-mulched area has been taken as a variable flux boundary where water flow is restricted to 0 cm d À1 and heat flow occurs from the plastic film to soil. The study thus simulates fully coupled water and heat flows in a rainfed soybean field.
Simulations with HYDRUS-2D under drip irrigation with ridge-furrow/raised-bed mulching were done in several regions of the world (Saglam et al. ; He et al. ; Parihar et al. ). However, in order to understand the effects of film mulching comprehensively, both plastic mulching having a potential of considerable open-hole area in plastic mulching and bare soil need to be modeled. Therefore, a comprehensive research must investigate the complete coupling of water-and heat-flow regimes with plastic mulching. With this view, the coupled-flow regimes of water and heat as well as water-and heat-distribution patterns under plastic mulching was simulated by the HYDRUS-2D model and the results were compared to that of bare soil in a rain-fed soybean field of central Japan. Furthermore, the effects of plastic mulching on soil-water balance were also quantified by using the established model.

Description of field experiment
A field experiment was done at Gifu University farm in Japan (35 27 0 N latitude and 136 44 0 E longitude, 12 m above sea level) by growing a rainfed soybean cultivar (Glycine max cv. Meguro) under two treatments of plastic mulching, and no mulching (bare soil). The mean air temperature at the study site is 16.1 C and annual rainfall is 1,849 mm. Figure 1 illustrates the layout of the treatments and setting of various devices in the experimental field (Kader et al. b). The research field was divided into two plots of size 12.5 m 2 (each 5 m long and 2.5 m wide) with a buffer zone of 0.5 m surrounding each plot. A raised bed was prepared in the plots and each treatment was employed in one plot. A silver-colour plastic film (one layer of 20 μm thickness) was selected for the plastic mulching and spread manually over the entire raised bed. Soybean was grown in four rows in each plot, with row to row and plant to plant spacing of 40 and 30 cm, respectively. The diameter of the planting hole for each soybean plant was 3 cm (Figure 1). The experiment was conducted under natural rainfall (rainfed) without any irrigation; a total of 874 mm rainfall occurred during the soybean-growing period. The average pan evaporation was 3.9 mm d À1 (except the rainy days) during the study period. The distribution of rainfall over the crop period ( Figure 2) was used as input information for defining the boundary condition in the simulation.

Measurements
Relevant meteorological data like air temperature at 2 m height, rainfall, solar radiation, relative humidity and wind speed were measured at the experimental field during the soybean cultivation period. The hourly soil moisture and temperature were measured at 5, 15 and 25 cm depths from each treatment throughout the experimental period. Necessary sensors (5TM of Decagon Devices, Inc., USA) were installed horizontally at the centre, between the rows and between the columns of each treatment, and the entire raised bed of each plot was covered with a plastic film by ensuring maximum contact of the plastic film with the soil surface. Each depth had one sensor and the sensors were connected to Em50 data loggers (Decagon Devices, Inc., USA).
where θ is volumetric moisture content (cm 3 cm À3 ), h is pressure head (cm), K(h) is unsaturated hydraulic conductivity of the liquid phase (cm d À1 ), t is time (d), x and z are horizontal and vertical coordinates (cm), respectively and S(h) is a sink term referring to root-water uptake (d À1 ).
The unsaturated soil hydraulic properties were modeled for the entire flow domain by using the van Genuchten-

Mualem model (van Genuchten ) expressed by
and where θ s is soil-moisture content at saturation (cm 3 cm À3 ), θ r is residual soil-moisture content (cm 3 cm À3 ), K s is saturated hydraulic conductivity (cm d À1 ), α (cm À1 ) and n (À) are shape parameters of the soil-water retention curve, S e (À) is effective saturation, and l (À) is a pore connectivity parameter that is assumed to be 0.50 for all cases (Mualem ). Here, m is given by (1-1/n). Temperature dependence of soil hydraulic parameters was taken into account in the study.

Root-uptake functions
The root zone of soybean is usually shallow; the maximum depth with roots (50 cm where α(h, x, z) is water-stress response function (-) due to rootwater uptake that is expressed as a function of soil-water pressure head, b(x, z) is water-uptake distribution function (cm À2 ), T p is potential transpiration rate (cm d À1 ) and L t is surface length associated with transpiration (cm). The root-water uptake parameters were explained based on the Feddes model for bean (soybean crop) from HYDRUS internal database (Šimunek et al. ) and are given in Table 1.

Soil-heat flow
Neglecting the effects of water-vapour diffusion, two-dimensional heat-transport function is given by (Sophocleous ) is apparent thermal conductivity of the soil (W m À1 C À1 ), and C p and C w are volumetric heat capacities of moist soil and liquid water (J m À3 C À1 ), respectively. The thermal conductivity as a function of soil-moisture content was expressed by (Chung & Horton ) where b 1 , b 2 and b 3 are empirical parameters (W m À1 C À1 ).
In Equation (6), the volumetric heat capacity, originally given by de Vries (), was further modified by Šimunek where subscripts n, o and w represent solid phase, organic matter and liquid phase, respectively in the soil.

Model parameterizations
Geometry of HYDRUS-2D domain HYDRUS-2D with general type of geometry uses the Galerkin finite element method in the two-dimensional  spacing of 2 cm was generated by HYDRUS-2D, with a finer mesh size close to the soil surface.

Setting of boundary
The boundary and initial conditions used in the HYDRUS- where k is a constant related to the leaf extinction coefficient of soybean, which was assumed to be 0.49 (Adeboye et al.

)
and LAI is Leaf Area Index, which varied linearly from 0 to 3 during the development stages of soybean (Ahmad et al. ). The variation of daily E p and T p over the soybean-growing period is illustrated in Figure 2. Initial conditions at the beginning of the simulation period were set in terms of field-measured soil moisture and temperature at 5, 15 and 25 cm depths for each treatment. Three observation nodes located at 5, 15 and 25 cm depths represented soil moisture and temperature sensors (Figure 3(a)).

Heat-transport boundary
First-type (Dirichlet-type) boundary, specified by temporal changes in temperature of each treatment (Figure 3(b)) at the upper surface, was adopted for heat-transport simulation against the atmospheric boundary in water-transport simulation. Heat flows from the bare soil, and planting areas occurred to the atmosphere. So, the time-dependent temperature boundaries at the upper surface were the first-type, with a pointer vector boundary of 1 and 2 (Figure 3(b)).
The condition of point vector 1 was defined by air temperature, which was specified for bare soil, and planting areas.
The point vector 2 was defined by soil-surface temperature just below the mulch that was specified for the plastic mulch-covered area (Figure 3(b)). The soil surface temperature of the silver colour plastic mulched soil was estimated from another study reported by Kader et al. (). Both vertical sides of the flow domain were 'no heat flux' boundary and the bottom boundary was third-type (Cauchy-type), which is of zero gradient soil temperature (Figure 3(b)).

Parameter optimization
Numerical simulations of soil-water and heat flow under the two treatments were done for the period from 12 June to 27 August (77 days (Table 2). Then, water flow of each treatment was simulated by optimizing van Genuchten's parameters using the Levenberg-Marquardt non-linear minimization method (Šimunek et al. ). Soil hydraulic parameters from the laboratorymeasured data were used as initial soil-moisture retention functions and then θ r , θ s , α, n and K s (Equations (2) and (3)) were optimized to obtain good fitting of soil-moisture contents. The optimized soil hydraulic parameters are listed in Table 2. Afterward, the coupled water and heat flow simulations were run for each treatment, where the water-flow parameters were fixed and heat-flow parameters were estimated by inverse solutions. The thermal parameters in Equation (6) for the relationship between soil thermal conductivity and volumetric moisture content for each soil layer were finally optimized by inverse solutions (Table 3).
The longitudinal dispersivity (D L ) of the soil was optimized between 4 and 10 cm by inverse simulations under condition of fixed transverse dispersivity (D T ) of 1 cm for each soil layer.

Model performance
Coefficient of determination (R 2 ), as a statistical indicator, was used to reflect the degree of correlation between the simulated and observed values of soil moisture and soil temperature of each treatment. The performance of the model was evaluated by root-mean-square error, RMSE, as expressed by where N is number of observations, and P i and O i are the simulated and observed values, respectively, of soil moisture or soil temperature.

Soil-moisture retention curves
Soil-moisture retention curves from the inverse optimization and laboratory measurements at 0-10, 10-20 and 20-30 cm Values shown as mean ± SE (Standard Error). θ n , b 1 , b 2 , b 3 , D L and C n were optimized by inverse solution. depths are depicted in Figure 4. The inverse-optimized hydraulic parameters varied slightly from those obtained from the laboratory-measured data, while the overall trend of the soil-moisture retention curves was similar for all three depths (Figure 4). It is noted that the soil hydraulic parameters obtained in laboratory tests would not be appropriate since the soil was heterogeneous over the year due to cultivation and the parameters might vary over time. The optimized soil hydraulic parameters were largely different among the three treatments to reproduce the observed soil-moisture contents. The optimized soil-moisture content at saturation, θ s (Equation (2)), showed a decreasing pattern with increasing depth; θ s for 0-10 cm soil layer was the highest in the plastic mulch treatment and lowest in the bare soil. The water-retention characteristics of the soils were important since they influenced water flow during the simulation process. The laboratorymeasured soil hydraulic parameters (θ r , θ s , α, n; Equation (2)) are not always used for correct fitting in simulation study; they need to be updated in most studies (Zhao et al.

Soil moisture
The observed and simulated soil-moisture contents at 5, 15 and 25 cm depths in two treatments are compared in      (Table 4) (7) and (8)) were optimized by inverse simulations for each treatment (Table 3). Moreover, the differences between the measured and simulated soil temperatures could also have occurred due to the effect of specified surface and bottom heat transport boundaries.
Overall, the predicted results from the HYDRUS-2D model were fairly close to the measurements, and the model was able to accurately simulate the dynamics of soil temperature in the soybean field.

Spatial distributions of soil-moisture and temperature
Simulated distributions of soil moisture and temperature in the vertical sections of the soil profile for the two treatments are illustrated in Figures 7 and 8, respectively, for rainfall and no-rainfall (drying) days. Soil-moisture distribution greatly varied when rainwater entered the soil profile. After a rainfall event, soil moisture moved in the vertical (from surface to deeper soil layers) and horizontal directions. In the dry days, soil moisture constantly moved, mainly in the vertical direction ( Figure 7). The bare soil treatment provided higher soil moisture (wettest) at the surface during rainfall (20th day of simulation) followed by the plastic mulch, which provided the lowest soil moisture (driest) (Figure 7).
On the 31st day of simulation (dry period), lower soil moisture was obtained for all three treatments compared to the wetting days (Figure 7). During the dry period, soil moisture evaporated from the bare soil and also the soil moisture moved downward due to free drainage condition, while the plastic mulch restricted evaporation from the soil.

Soil-water balance
Water balance in the modeling with cumulative water fluxes in and out of the simulated flow domain is described in Table 5. The water balance components include water fluxes like cumulative infiltration, evaporation and transpiration, change in water storage and deep drainage from the nodes at the bottom of the soil profile during simulation.
The negative values of soil-moisture storage indicated water loss from the soil domain during the simulation period.
Rainwater infiltration was 61% higher in bare soil than in plastic mulching. Similarly, cumulative evaporation from the bare soil was 62% higher compared to plastic mulching.
The amount of cumulative evaporation was in the order: bare soil (126 mm) > plastic (48 mm). The bare soil caused greater infiltration of water and showed higher evaporation and bottom fluxes, indicating downward movement of water below the root zone. The change in soil-moisture storage was also the greatest in bare soil followed by plastic mulching. The water balance analysis revealed that the bare soil increased rainwater infiltration compared to plastic mulching but it reduced soil evaporation compared to the bare soil during the soybean-growing period. The plastic mulch provided direct pathways between soil and air, from which soil-moisture exited during evaporation; thus, size and density of holes might control the evaporation that occurred from the soil surface.

Modeling implications
The two-dimensional numerical simulations of soil-water and heat flows provided better understanding of infiltration and evaporation processes under plastic mulching.
Although difficult to choose a correct boundary for plastic-  In our research, we analyzed the effects of soil evaporation from plastic mulch by quantifying the soil-water balance.

Water flow
Generally, plastic film mulching restricts water flow into the soil profile. The two-dimensional HYDRUS model quantified the effect of plastic mulching on rainwater infiltration into soil by describing water flow. The temporal variations of cumulative total seasonal infiltration and evaporation of the two treatments are illustrated in Figure 9. The cumulative water infiltration into the bare soil was greater than that into the plastic mulch. In the plastic mulch treatment, rainwater infiltrated through the planting holes. The bare soil caused more infiltration than the plastic mulch treatments since rainwater fell directly on the soil surface (Kader et al. b). Cumulative evaporation was also greater in bare soil treatment than in plastic mulch treatment ( Figure 9). In Figure 9, the fluctuations of soilmoisture storage (ΔS) are higher in bare soil than in plastic  (Parihar et al. ). In our rainfed crop cultivation system, rainfall was the only input of water balance where water loss included evaporation from soil, transpiration from the rooting zone, deep percolation below the root zone and lateral water movement due to surface runoff. A large part of infiltrated water in the bare soil moved to the deeper soil layers by free drainage (Table 5). Thus, the plastic mulch greatly reduced evaporation from the soil surface compared to bare soil.

CONCLUSIONS
The two-dimensional HYDRUS model first time successfully mulching. This study confirmed that plastic mulch treatments reduced the changes in soil-moisture storage and soil evaporation compared to bare soil treatment, indicating greater crop-water availability by mulching. Therefore, the HYDRUS-2D model may be useful to optimize the best density of holes in the plastic mulching by scenario analysis for enhancing soil-moisture conservation under different climatic conditions; this topic can be investigated in future studies.