Sustainable irrigation based on co-regulation of soil water supply and atmospheric evaporative demand

Field measurements

We used two sets of field measurements of soil moisture, VPD, and stomatal conductance of maize at the daily scale to illustrate a proof-of-concept for the co-regulation of soil moisture and VPD on stomatal conductance.

The first set was measurements from greenhouse experiments of maize (seed: Dekalb hybrid DKC52-04) at Colorado State University during the 2013 growing season (planted on June 10, 2013)⁴⁹. There were two treatments (well-watered, WW, and water-stressed, WS) with five plants per treatment. The soil of the greenhouse experiments was the air-dried soilless substrate (8.8 kg) consisting of a 1:1.3 by volume ratio of Greens GradeTM, Turface^® Quick Dry^® and Fafard 2SV in 26 L pots⁴⁹. The soil moisture measurements came from soil moisture sensors (Decagon5TM sensors) installed in the middle of the pots (~6 inches from top). The greenhouse measurements of leaf-level stomatal conductance and soil moisture were performed in approximately 2-week intervals beginning in the vegetative stage and continuing until plant senescence (DOY 198–199, 210–211, 217–218, 233–234, 247), with 11 replicates for each plant under two treatments (WW and WS). The environmental variables, such as relative humidity and air temperature, were continuously measured in minutes. Other detailed experimental setups can be found in Miner and Bauerle (2017)⁴⁹.

The second set was eddy-covariance measurements of maize cropping systems (seed: Pioneer 33P67/33B51) from 2001 to 2012 at three AmeriFlux sites (US-Ne1, Ne2, and Ne3). US-Ne1 and Ne2 were irrigated sites, with a continuous maize cropping system during 2001–2012 for US-Ne1 and with a maize-soybean rotation cropping system during 2001-2009 and then a continuous maize cropping system during 2010-2012 for US-Ne2. US-Ne3 was rainfed with a maize-soybean rotation cropping system during 2001–2012. The soil at the three AmeriFlux sites was a deep silty clay loam consisting of four soil series: Yutan, Tomek, Filbert, and Filmore. There are three replicates with the soil moisture sensors (theta probes: ML2, Dynamax Inc.) installed horizontally with the profile of soil depth (10, 25, 50, and 100 cm) in the US-Ne1 and US-Ne2, and four replicates with soil moisture sensors (theta probes: ML2, Dynamax Inc.) installed horizontally with the profile of soil depth (10, 25, 50, and 100 cm) in the US-Ne3 (http://csp.unl.edu/public/G_moist.htm). The soil moisture data used here was from the top soil layer (10–25 cm). The canopy-level stomatal conductance (G_s) was derived by inverting the Penman-Monteith equation⁵⁰ (Equations 1 and 2) from the eddy-covariance measurements at the hourly scale^18,24,51, and the averaged value near midday (from 12:00 to 14:00) was applied as the daily canopy-level stomatal conductance to remove the diurnal cycle. This inversion was only conducted during peak growing season (July and August) to avoid the impact of LAI²⁴. The impact of evaporation from canopy interception and of low incoming shortwave radiation was removed by data filtering²⁴, i.e., excluding the data within 2 days following every precipitation and irrigation event, and periods of low incoming shortwave radiation conditions (<500 Wm⁻²).

where Δ is the slope of the water vapor deficit; R_n and G are net radiation and soil heat flux, respectively; (rho ) is the air density; c_p is the specific heat capacity for dry air; g_a is aerodynamic conductance; (gamma ) is the psychometric constant; and (lambda E) is evapotranspiration.

Advanced process-based model (ecosys)

We further used an advanced process-based agroecosystem model, ecosys, to reproduce the co-regulation pattern at the daily scale and to investigate the performances of different irrigation schemes with continuous maize cropping systems across 12 sites in Nebraska with a large rainfall gradient under current climate (2001–2019) and RCP-8.5 scenario (2058–2076) (Fig. 3 and Table S1). Ecosys simulates coupled energy, water, carbon, and nutrient cycles^{52,53,54,55,56}, and has been extensively tested in various agricultural ecosystems^{52,53,54,56,57}. Ecosys could simulate all major agricultural management practices, such as tillage⁵⁸, crop rotation⁵⁸, fertilizer⁵⁹, and irrigation⁶⁰.

Physical processes

Ecosys uses a multilayered soil-root-canopy system to get hourly two-stage convergent solutions for crop carbon assimilation, water uptake, and energy fluxes^{52,53,54,55,60,61}. The first stage focuses on the convergence of canopy temperature for the first-order closure of canopy energy balance in Equation 3, including net radiation R_n, latent heat flux LE (including from evaporation, LE_v in Equation 4a, and transpiration, LE_c in Equation 4b), sensible heat flux H, and soil heat flux G. Canopy latent heat is controlled by aerodynamic resistance (r_a)^52,53 and canopy stomatal resistance (r_c). Canopy stomatal resistance is postulated by two controlling mechanisms: canopy photosynthesis (leaf level driven by rates of carboxylation vs. diffusion, Equation 5) and canopy turgor potential (ψ_t, canopy level constrained by water status, Equation 6). The second stage focuses on the convergence of canopy water potential (ψ_c) and thereby r_c (=1/g_s) at which transpiration (E_c) based on canopy energy balance (left term in Equation 7) equals root water uptake from multiple soil layers (right term in Equation 7). Root water uptake is calculated using water potential differences between soil (ψ_s) and root (ψ_r), and the soil and root hydraulic resistances (Ω_s and Ω_r) in each rooted soil layer.

where r_c(min) is the minimum r_c at ψ_c= 0 MPa; r_c(max) is canopy cuticular resistance to vapor flux; C_b is [CO₂] in canopy air; C_i΄ is [CO₂] in canopy leaves at ψ_c= 0 MPa; V_c΄ is the potential canopy CO₂ fixation rate at ψ_c= 0 MPa; ψ_π is canopy osmotic potential; e_a is atmospheric vapor density at air temperature (T_a) and ambient humidity; e_c is canopy vapor density at canopy temperature (T_c) and ψ_c; Ω_s,r,l is radial resistance to water transport from soil to surface of roots or mycorrhizae; Ω_r,r,l is radial resistance to water transport from the surface to axis of roots or mycorrhizae; Ω_a,r,l,x is axial resistance to water transport along axes of primary (x = 1) or secondary (x = 2) roots or mycorrhizae; l is soil or canopy layer; r is root or mycorrhizae; and X_c is canopy capacitance.

Photosynthesis at the leaf-level is calculated using the Farquhar model for C3 plants and the Farquhar model plus mesophyll-bundle sheath carbon exchange for C4 plants with specific azimuth, leaf inclination, exposure of light conditions (i.e. sunlit and shaded leaves), and canopy height. Canopy photosynthesis is the sum of the photosynthesis of all individual leaves. The carbohydrate is then allocated for maintenance respiration (R_m) in both the shoot and root, and the remainder for growth respiration (R_g), and dry mass (DM) formation. The phenologically-driven plant carbon allocation ratio among shoot and root organs is impacted by the number of phyllochron intervals and the water and nutrient status of the plant. DM of shoots is partitioned to seven organs (leaf, sheath, stalk, soluble reserves, husk, cob, and grain) with dynamic partitioning ratios varying with growing stages. Seed number and kernel mass are set during postanthesis growth stages to determine the yield upon harvest. More details about the biophysical and biochemical processes in ecosys can be referred to the supplement of Grant, et al (2020)⁶¹.

Model setup

Ecosys was validated at three AmeriFlux sites (US-Ne1, Ne2, and Ne3, https://ameriflux.lbl.gov/)⁶⁰ and 12 sites across Nebraska (Table S5, Fig. S12–S15).

Three AmeriFlux sites (US-Ne1, Ne2, and Ne3) had the complete data from 2001 to 2012 to test model performance, including the meteorological variables (such as surface air temperature, downward shortwave radiation, wind speed, and precipitation), eddy-covariance fluxes (such as GPP, net ecosystem exchange-NEE, and LE) from the FLUXNET2015 dataset⁶² (http://fluxnet.fluxdata.org/data/fluxnet2015-dataset/) and detailed ground-based crop growth observations (i.e. planting/harvest date, irrigation/fertilization records, LAI, and yield) from the Carbon Sequestration Program (CSP) at University of Nebraska-Lincoln’s Agricultural Research and Development Center (http://csp.unl.edu/Public/sites.htm)⁶³. The daily and monthly model simulations of GPP, ET, LAI, and yield matched well with the eddy-covariance and ground-based observations (Table S5, Fig. S12–S15).

The hourly meteorological variables, including surface air temperature, humidity, wind speed, precipitation, and downward shortwave radiation, from 1979 to 2019 (1979–2000 set for model spin-up and 2001–2019 set as current climate conditions) across 12 sites in Nebraska were obtained from the North American Land Data Assimilation System (NLDAS-2). For climate change, the ensemble predictions of air temperature, precipitation, and carbon dioxide concentration [CO₂] at the same irrigated sites in Nebraska were ensemble average projections from 15 Coupled Model Intercomparison Project phase 5 (CMIP5) models under Representative Concentration Pathway 8.5 (RCP-8.5)⁶⁴ (Table S4). RCP-8.5 was selected as the future scenario to investigate the high warming conditions with the highest greenhouse gas (GHG) emissions⁶⁵. ecosys model simulation under the current climate (2001–2019) was extended through three additional 19-year cycles from January 1, 2020 to December 31, 2076 under the RCP-8.5 scenario using the incremental change scheme without the acclimation of parameters in ecosys, and the fourth cycle 2058–2076 was selected as the future climate to investigate the climate change effects. It should be noted that the acclimation of parameters in ecosys could be applied to maintain high maize yields with increased air temperature under climate change, like the high maize yields under hot, semiarid to arid climate with intensive irrigation in the Texas Panhandle. As we aimed to investigate the impacts of environmental variables to maize under climate change, no acclimation of parameters in ecosys is suggested under climate change.

The soil information (i.e., bulk density, field capacity, wilting point, soil texture, saturated hydraulic conductivity, soil organic carbon, pH, and cation exchange capacity) was acquired from the Gridded Soil Survey Geographic Database (gSSURGO) dataset⁶⁶. There were 12 soil layers (the depth of the bottom for each layer: 0.01, 0.05, 0.1, 0.15, 0.18, 0.28, 0.35, 0.59, 0.92, 1.32, 1.6, and 2.00 m) with a maximum depth of 2.0 m. The planting date of the continuous maize cropping systems at the 12 sites was obtained from the USDA NASS weekly Crop Progress Reports (2001–2019) with the fertilizer (18 g N m^-2 and 5 g P m^-2 per year) applied two days before planting, and the crops were harvested on October 31. Other land management practices were set as the same across the 12 irrigated sites in Nebraska, including planting density (8.4 plants m^-2), tillage practice (no tillage), and crop type (continuous maize cropping systems). The auto-irrigation scheme in ecosys with the widely used soil-based MAD-50% was applied to determine the irrigation scheduling at 12 sites to test the model performance. The NASS county-level irrigated maize yield (2010–2019) at the counties where the 12 sites were located was used as the observations. The probability density function with Gaussian kernel density estimation of irrigated yields from ecosys model at 12 sites showed good agreement with that of the National Agricultural Statistics Service (NASS) county-level irrigated yields during the period from 2010 to 2019 (Fig. S15h), which further validated the reliability of the ecosys model.

Empirical nonlinear statistical model

There were linear and nonlinear empirical models available for modeling stomatal conductance^20,24,67. As the exponent of VPD was close to one for croplands²⁵ and the linear function between the stomatal conductance and soil water potential led to the logistic function between the stomatal conductance and soil moisture, we used an empirical nonlinear statistical model to describe the co-regulation of soil moisture and VPD on stomatal conductance (Equation 8), including two sub-functions²⁴. The first one denoted an inverse proportional relationship between the VPD and stomatal conductance^25,68, and the other represented the logistic function between soil moisture and stomatal conductance⁶⁹. The impacts of [CO₂], nutrient, radiation, and temperature on stomatal conductance were not considered in this research.

The field measurements of soil moisture, VPD, and stomatal conductance at the daily scale from the greenhouse experiments and three AmeriFlux sites (US-Ne1, Ne2, and Ne3) were applied to investigate the co-regulation of soil moisture and VPD on stomatal conductance (Equation 8) as observational evidences. With the validated ecosys model at three AmeriFlux sites (US-Ne1, Ne2, and Ne3) and 12 sites across Nebraska (Table S5, Fig. S12–15), the model simulations of daily soil moisture, VPD, and canopy-level stomatal conductance during the peak growing season (July and August) were applied to investigate the co-regulation patterns across 12 sites in Nebraska. It needs to be noted that the fitted parameters under the current climate (2001–2019) should be updated under RCP-8.5 scenario, as maize may respond differently under climate change due to increased air temperature and VPD.

In addition, Lindeman, Merenda, and Gold method (LMG)³⁰ was used to identify the relative importance of soil moisture and VPD on the stomatal conductance. LMG decomposed the determination coefficients of a linear regression (R²) to the contributions of soil moisture and VPD, i.e., to quantify the variation of stomatal conductance that can be explained by soil moisture and VPD, while taking the correlation between the soil moisture and VPD into account.

Agricultural irrigation management

We proposed the plant-centric irrigation scheme based on water supply-demand dynamics (SDD). By using ecosys model simulations, we further compared its performance with the widely used soil moisture-based irrigation scheme, i.e. management allowable depletion (MAD) under both current (2001–2019) and future climate conditions (RCP-8.5 scenario, 2058–2076). The universal (under current climate and RCP-8.5 scenario) and site-specific (under current climate) parameters of SDD and MAD irrigation schemes across 12 sites in Nebraska were optimized for maximizing economic profit (Equation 9). For the irrigation setting in the ecosys model, irrigation was triggered when soil moisture was lower than the soil moisture threshold of SDD (varying with VPD) and MAD (constant) at the daily scale during the growing season. Soil moisture was the weighted average of soil water content in different soil layers within the root zone, which was varying with dynamic root growth but no more than the top 9 soil layers with a depth of 0.92 m to reduce the impacts of deep wet soil. For simplification, we ignored the constraints on irrigation amount and duration from irrigation infrastructures. Irrigation amount was determined by water required to fill current soil water content to field capacity. Each irrigation event lasting 24 h was incorporated into the ecosys model at a daily time step in real time. Furthermore, we used two indexes, economic profit and irrigation water productivity, to evaluate the performances of SDD and MAD irrigation schemes. Economic profit was the net revenue based on marketable yields and costs, including irrigation costs and fixed costs of production (Equation 9). Irrigation water productivity (IWP)⁷⁰ was the ratio between the marketable yields and irrigation amount during the growing season (Equation 10). To reduce the impacts of parameters in performance assessment, we used the recorded parameters at Nebraska, United States in 2019 (Table S3). Relative and absolute differences in irrigation amount, yield, profit, and irrigation water productivity between SDD and MAD irrigation schemes under current climate (2001–2019) and RCP-8.5 scenario (2058–2076) were calculated (Equations 11 and 12).

where y is maize yield (t ha^-1); p_maize is the price of maize ($ t^-1); I is the irrigation amount (mm); ({varGamma }_{irrigation}) is the price of irrigation ($ m^-3); λ is the irrigation application efficiency of the center pivots; K_fixed is the fixed costs of production ($ ha^-1), including the costs of seeds, fertilizer, storage, and so on; and IWP is the irrigation water productivity (kg m⁻³).

Source: Ecology - nature.com