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

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. 2017). Plastic film mulching, considered as water-saving technology, is a well-known farming practice in dryland areas (Qin et al. 2018) for regulating soil temperature and conserving soil moisture by reducing soil evaporation (Kader et al. 2017a). So far, China has become the largest user of plastic film mulching (Daryanto et al. 2017), which is mainly popular in the Loess Plateau area where frequent water shortage occurs due to splash rainfall (Wang et al. 2015; He et al. 2016). Plastic mulch improves soil-moisture retention and availability, and exhibits soil-warming effects by regulating radiation and thermal properties (Tian et al. 2003; Kader et al. 2020). 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. 2014). 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 al. 2018; Zhang et al. 2018). The application methods of plastic mulching in the field have been practiced by ridge furrow (Zhao et al. 2014), raised bed (Alliaume et al. 2017; Zhang et al. 2018), partial ridge and furrows (Chen et al. 2018), double layer mulching (straw incorporated with plastic) (Yin et al. 2016) and traditional flat film mulching (Zhang et al. 2015) with various thicknesses (10–20 μm) and colours. 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. 2018). 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. 2017). 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. 2018). 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. 2015; Steinmetz et al. 2016). 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. 2018). HYDRUS is a popular software that can simulate water, heat and nutrient dynamics in variably-saturated porous media (Šimůnek et al. 2016). The HYDRUS-2D version has superiority over its 1D version in choosing various boundaries and providing water-flow distribution into soils. Simulations with HYDRUS-2D provides quick and accurate estimation of the hydro-thermal environment of field soils (Chen et al. 2018; Rai et al. 2019). However, selection of accurate soil water- and heat-flow boundaries in HYDRUS-2D software under plastic mulch cover is critical. In previous researches (Saglam et al. 2017; Chen et al. 2018), a ‘no flux’ water-flow boundary was considered in the plastic-mulched area. Moreover, several studies (e.g. Saglam et al. 2017; He et al. 2018; Qi et al. 2018), simulating only water flow using HYDRUS-2D under drip emitters with plastic mulching, did not discuss heat-flow distribution, which is influenced by mulch. Previous researches (e.g. Zhang et al. 2018; Zhao et al. 2018) considered a no-flux boundary for 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 no-flux 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⁻¹ 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. 2017; He et al. 2018; Parihar et al. 2019). 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.

A field experiment was done at Gifu University farm in Japan (35° 27′N latitude and 136° 44′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. 2017b). The research field was divided into two plots of size 12.5 m² (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⁻¹ (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.

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). Soil-surface temperature (at 0 cm depth) was calculated by correlating air temperature in 2015 and regression co-efficient (R²) of soil-surface temperature from another experiment done in 2016 (Kader et al. 2019). Soil samples were collected after the experiment from 0–10, 10–20 and 20–30 cm soil layers, with three replications for each treatment by using core samplers of 100 cm³. Texture of the 0–30 cm soil profile was loamy sand (53% sand, 31% silt, 16% clay); the classifications of the soil were described by Kader et al. (2017b). Soil-water retention functions were determined for 0–10, 10–20 and 20–30 cm soil layers with three replications in the laboratory by centrifugation method (Russel & Richards 1938) with a Kokusan H-2000B centrifuge machine at equilibrium matric potentials of –10, –32, –100, –316 and –1,500 kPa. The saturated hydraulic conductivity (K_s) of the field soil at 5, 15 and 25 cm depths with three replications was measured in the laboratory by the falling head method (Klute & Dirksen 1986).

HYDRUS-2D with general type of geometry uses the Galerkin finite element method in the two-dimensional vertical x-z plane. The transport domain was 250 cm deep and 250 cm wide where the length of plastic mulching was 200 cm (Figure 3(a)). Due to the deep groundwater table (much below the simulation domain) at the study site, a free drainage boundary was considered at the bottom (at 250 cm depth). A finite element mesh with radical grid spacing of 2 cm was generated by HYDRUS-2D, with a finer mesh size close to the soil surface.

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. (2019). 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)).

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) in 2015. Firstly, the soil hydraulic parameters were estimated from the laboratory-measured soil-moisture retention data of each treatment using RETC software (van Genuchten et al. 1991). The estimated parameters were used as initial values in inverse solution for estimating the most appropriate soil hydraulic parameters (Table 2). Then, water flow of each treatment was simulated by optimizing van Genuchten's parameters using the Levenberg-Marquardt non-linear minimization method (Šimůnek et al. 2018). Soil hydraulic parameters from the laboratory-measured 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.

Soil-moisture retention curves from the inverse optimization and laboratory measurements at 0–10, 10–20 and 20–30 cm 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 laboratory-measured 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. 2018).

The observed and simulated soil-moisture contents at 5, 15 and 25 cm depths in two treatments are compared in Figure 5(a) and 5(b), respectively. The fluctuation of soil moisture was the largest in bare soil followed by plastic mulch treatment. The plastic mulch showed moderate fluctuation of soil moisture. Rainfall was the main factor in altering soil moisture in all treatments. The bare soil, due to its whole open surface, received more rainwater than other treatments and showed the greatest fluctuations in soil moisture. The plastic mulch protected rainwater to a large extent and showed moderate fluctuations of soil moisture over the soybean-growing season. In plastic mulch treatment, the fluctuation of soil-moisture content at 5 cm depth was smaller than at 15 and 25 cm depths (Figure 5(a)). The variation of soil-moisture content in the bare soil was obviously larger near the soil surface (at 5 cm depth) since the effects of rainfall on evaporation and root-water uptake was greater at the soil surface than at the deeper layers.

The statistical performance of HYDRUS-2D model for water-flow simulation is summarized in Table 4 in terms of root-mean square error, RMSE, and coefficient of determination, R². The model provided fairly good estimate of soil-moisture contents, with RMSE between the measured and simulated soil-moisture contents of 0.015–0.030, and 0.046–0.055 cm³ cm⁻³ for the plastic mulch and bare soil treatment, respectively (Table 4). Moreover, the best fitting between the measured and simulated soil-moisture contents was obtained for the plastic mulching (R² = 0.62) followed by bare soil (R² = 0.47). There were little differences between the observed and simulated values, possibly due to close proximity of the sensors' positions in the field and observation nodes in the simulation model. Our simulation model performed well, particularly when considering the complexity of the conditions like heterogeneous soil properties, plastic mulch-covered soil, several rainfall events and non-homogeneous root distribution of soybean plants to which the model was applied. Finally, we confirmed that HYDRUS-2D model reflected reliability in simulating soil-moisture contents at the soil root zone (0–30 cm) depth.

Figure 6(a) and 6(b) illustrate the measured and simulated soil temperatures at 5, 15 and 25 cm depths under mulching and bare soil treatments over the soybean-growing period. With the occurrence of rainfall, soil temperature decreased from top to bottom of the soil profile. The differences in soil temperature between the plastic treatment was minimal during the crop period. The simulation of soil temperatures with thermal parameters, which were estimated by heat-transport inverse model, was reasonably good (Table 4, Figure 6). The model performed fairly well, with RMSE of 0.66–1.54 °C and R² of 0.70–0.95 for the two treatments (Table 4). The plastic mulching provided greater modeling efficiency, with RMSE varying from 0.90 to 1.32 °C at 5, 15 and 25 cm soil depths compared to RMSE 0.70–1.54 °C for bare soil (Table 4). Bare soil received rainfall and sunlight directly on the surface of the soil and part of the sunlight was blocked by the leaves, which causes higher degrees of fluctuation of soil moisture and temperature compared to plastic treatment. So, it is difficult to accurately capture the trend of fluctuations of the observed soil moisture and temperature of bare soil by HYDRUS-2D simulations. For this reason, the simulation model showed inferior performance in simulating soil moisture and temperature in bare soil (Table 4) compared to the plastic mulch treatment. Due to the lack of directly measured data of thermal conductivity and volumetric heat capacity of the soil, we optimized a minimum number of heat transport parameters in our simulations. Only the least number of heat-transport parameters (D_L, b₁, b₂, b₃ and C_n in Equations (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.

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.

With increasing rainfall, both the soil moisture and temperature changed significantly at the soil surface and moved in the vertical and horizontal directions from the planting and mulching holes. Rainwater could not infiltrate directly into soil through the plastic mulching as it did through the bare soil. Our simulation revealed that, with 97.5 mm rainfall on 30 June (20th day of simulation), the bare soil was completely saturated; plastic mulching partially wetted the soil, especially the planting-hole area. On the other hand, the drying day on 25 July (31st day of simulation) showed reduced soil moisture at the rooting zone under all three treatments in the order: plastic mulch > bare soil.

Figure 8 focuses on the distribution of soil temperature on day 49 (rainfall) and 56 (no rainfall, dry day) of simulations. On day 49, the variation of soil temperature was very small, especially in bare soil compared to the mulching treatments. Large variation in soil temperature was observed on the drying day 56. The distribution of soil temperature was similar for the plastic mulching, while the bare soil showed steady change in temperature within the soil rooting zone (Figure 8). The surface mulch cover had a greater effect on the soil environment near the surface compared to the deeper soil layers, with the development of strong horizontal gradients in temperature across different mulched-soil boundaries (Zhao et al. 2018). Therefore, the spatial distributions of soil moisture and heat flow are influenced by the shadow between two plants, the presence of plastic film mulch and the root density of soybean plants.

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.

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-mulched soil in the HYDRUS-2D model (Chen et al. 2018; Zhao et al. 2018), our study accurately selected the boundary conditions for the plastic mulched area, including planting holes for simulating coupled flow of water and heat under two different treatments. Therefore, the HYDRUS-2D model provided an efficient way to simulate coupled distribution of water and heat under a plastic-mulched field in a rainfed system. In our HYDRUS-2D modeling, water-flow boundary was followed by heat-flow boundary, and the model performed better in heat-flow simulation than water-flow simulation (Table 4). However, similar to the finding of Xi et al. (2016), infiltration of rainfall and soil evaporation considerably influenced soil-moisture contents in the upper soil layers compared to the deeper soil layers and, consequently, resulted in some deviations in the simulated soil-moisture contents. The simulation results can, however, be substantially improved by calibrating the input parameters (effective soil hydraulic parameters) against a dynamic flow experiment (Kandelous et al. 2011; Rai et al. 2019). Previous researches (Han et al. 2015; Saglam et al. 2017) simulated plastic mulch with drip irrigation, in which soil-water dynamics were controlled by both irrigation and rainfall. This study reveals that the performance of the HYDRUS-2D model strongly depended on rainfall and selection of various surface boundaries in the plastic mulch area; in this study, canopy distribution was not considered in the modeling. Moreover, the simulation of soil-surface temperature under plastic mulch is affected by latent heat, which has received limited attention in the present model. This is because it is a complex system, which involves the optical properties of the film, the thickness of the film and incoming solar radiation (Zhao et al. 2018). In the heat-transport models of previous research (Zhang et al. 2018), soil evaporation was neglected since the soil surface was fully covered with plastic film mulch. In our research, we analyzed the effects of soil evaporation from plastic mulch by quantifying the soil-water balance.

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. 2017b). Cumulative evaporation was also greater in bare soil treatment than in plastic mulch treatment (Figure 9). In Figure 9, the fluctuations of soil-moisture storage (ΔS) are higher in bare soil than in plastic mulching. The smaller changes of soil-moisture storage by plastic mulching indicate minimized soil-moisture fluctuations, which is important for crop growth. Larger soil-moisture fluctuations may be induced frequently under soil wetting and drying conditions that have a high risk of getting good crop growth and yield (Parihar et al. 2019). 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.

The two-dimensional HYDRUS model first time successfully simulated soil-water dynamics and heat transport under plastic mulching in a rainfed soybean field in central Japan. In this study, water and heat distribution characteristics in soils were evaluated by HYDRUS-2D simulations under plastic and bare soil treatments. The results suggest that the measured and simulated soil moistures and temperatures are in good agreement and able to accurately capture the dynamic changes in water and heat at 5, 15 and 25 cm soil depths. The results of the soil-moisture balance of each treatment show that the bare soil can greatly increase rainwater infiltration compared to plastic 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.

The authors gratefully acknowledge the Gifu University of Japan for providing agricultural land to conduct the field experiment. The first author thanks the Japanese Government for providing financial support during his doctoral study through the MEXT scholarship.

None.

All relevant data are included in the paper or its Supplementary Information.