PLANT DENSITY, IRRIGATION AND NITROGEN MANAGEMENT: THREE MAJOR PRACTICES IN CLOSING YIELD GAPS FOR AGRICULTURAL SUSTAINABILITY IN NORTH-WEST CHINA

Highlights

• A relative yield of 70% was obtained under both border and drip irrigation.

• Drip irrigation saved water and lowered yield variability compared to border irrigation.

• Drip irrigation led to accumulation of soil nitrogen and phosphorus in the root zone.

• Relative yield may increase 8% to 10% by optimizing field management.

• Plant density, irrigation and nitrogen are major factors closing yield gap in NW China.

Abstract

Agriculture faces the dual challenges of food security and environmental sustainability. Here, we investigate current maize production at the field scale, analyze the yield gaps and impacting factors, and recommend measures for sustainably closing yield gaps. An experiment was conducted on a 3.9-ha maize seed production field in arid north-western China, managed with border and drip irrigation, respectively, in 2015 and 2016. The relative yield reached 70% in both years. However, drip irrigation saved 227 mm irrigation water during a drier growing season compared with traditional border irrigation, accounting for 44% of the maize evapotranspiration (ET). Yield variability under drip irrigation was 12.1%, lower than the 18.8% under border irrigation. Boundary line analysis indicates that a relative yield increase of 8% to 10% might be obtained by optimizing the yield-limiting factors. Plant density and soil available water content and available nitrogen were the three major factors involved. In conclusion, closing yield gaps with agricultural sustainability may be realized by optimizing agronomic, irrigation and fertilizer management, using water-saving irrigation methods and using site-specific management.

Graphical abstract

Keywords

boundary line analysis / irrigation method / precision agriculture / spatial variability / yield gaps / yield-limiting factors

1 INTRODUCTION

Two main challenges facing agriculture are food security and environmental sustainability^[1,2]. A 60% to 110% increase in global crop production (from the 2005 to 2007 levels) may be needed by the year 2050 because of the increasing population and energy needs^[3,4]. Slowing down of rate of increase in cereal yield or even yield stagnation has also been a global trend in recent years^[5]. Increasing crop production without expanding the land area used involves increasing productivity on the existing arable land. One way is to increase potential yields and another is to close the yield gap between those achieved by farmers and potential yields^[1,6]. However, increasing potential yields is challenging because crop harvest index and radiation interception efficiency have been greatly improved since the first green revolution. Great room for improvement exists in increasing the conversion efficiency of intercepted radiation into biomass energy, but there is a two- to three-decade gap between the conceptual design at the experimental level and availability to farmers^[7]. It is vital to close existing yield gaps for food security.

Crop yields at the field scale are the comprehensive result of multiple factors. For example, climatic conditions and cultivar features define the potential yield, water and nutrients limit the potential yield, and weeds and pests serve as reducing factors^[8]. Irrigation and nutrient management are essential parts of crop production. Global agriculture consumes about 70% of the freshwater withdrawals^[1,9]. Global water consumption is projected to increase by 21% by 2050^[9] and domestic, industrial, and livestock water demands are still growing rapidly. Thus, shrinkage of the share for irrigation may be inevitable in the face of limited water resources. Climate change, ecosystem degradation, unsustainable use of groundwater, and pollution exacerbate water scarcity^[9]. Global fertilizer use also faces serious problems as in some regions insufficient nutrients are a major factor limiting crop yields. However, excessive fertilization has caused environmental problems in other regions, especially in East Asia^[1,2,10]. Increasing crop yields with sustainability requires greater efforts to increase water and fertilizer use efficiency.

Numerous studies have investigated the interrelationships among crop yield, resource use efficiency and multiple variables at different scales. Agronomists have screened minimum data sets (MDS) to assess correlations between soil quality and yield at the farm or regional scales^[11–13]. The MDS may include soil physical variables (e.g., bulk density and clay, silt and sand contents), soil chemical properties (e.g., pH in water and organic carbon, total nitrogen, ammonium nitrogen, exchangeable calcium and available phosphorus contents), and soil microbial variables (e.g., sucrase activity, earthworm numbers, relative amounts of bacterivorous nematodes, microbial biomass C and biomass N). Mueller et al.^[2] examined yield gaps at the global scale and concluded that 75% of attainable yields could be realized by increasing irrigation area and nutrient application on all underachieving areas. Furthermore, minimal changes in total worldwide nitrogen and phosphate use would be possible by reducing nutrient imbalances and inefficiencies. Chen et al.^[14] evaluated factors affecting summer maize yield and nitrogen use efficiency on small farms on the North China Plain and concluded that plant density was more important than other management practices in closing yield gaps. Li et al.^[15,16] analyzed the spatiotemporal distribution of irrigation water productivity and its driving factors in the Hexi Corridor, north-west China. They concluded that increasing irrigation water productivity should rely on the use of water-saving technology, optimum fertilizer rates, agricultural mulching film, and pesticides. Although a large amount of research has been conducted, sustainable agricultural development in China continues to face many challenges such as water shortage, low resource use efficiency, and environmental contamination^[17]. Moreover, few system analyses and interdisciplinary studies have been conducted.

Further information is needed on current crop production conditions to increase crop yield, reduce agricultural resource inputs or increase resource use efficiency, and reduce the impact of agriculture on the environment. For example, the extent of current yield gaps, the major causes of these, and how much they can be closed and by what means are all poorly understood. Boundary line analysis (BLA) has been widely used to determine critical levels of nutrients^[18,19] and analyze yield gaps^[20]. It simplifies a multivariate model into analysis of the relationships between the dependent variable and each independent variable to make interpretation easier^[21]. The objectives of the present study were to investigate current yields and yield-limiting factors in maize production at the field scale, determine the major yield-limiting factors and quantify the yield gaps using BLA, and compare the feasibility of and measures for closing yield gaps and maintaining agricultural sustainability.

2 MATERIALS AND METHODS

2.1 Experimental site and field management

The experiment was conducted on a 3.9-ha maize seed production field during the growing seasons of 2015 and 2016 at Shiyanghe Experimental Station of China Agricultural University (37°52′ N, 102°50′ E, 1581 m), located in Wuwei City, Gansu Province, north-west China. This region is located in an arid inland temperate climatic zone with mean annual pan evaporation of 2000 mm, mean annual precipitation of 164 mm, and an average groundwater table lower than 25 m below the ground surface^[22]. Precipitation during the maize growing season was 149 and 115 mm in 2015 and 2016, respectively.

Field management was operated by a farmer according to local practices. The cultivars (the female cultivar was FN1390 in 2015 and WH01 in 2016) and agronomic recommendations were provided by local seed companies. The field was partly covered by plastic film mulch, with every 1.2 m width covered strip across the whole field was separated by 0.46 m width uncovered strip. In 2015, every four rows of the female parents with one row of the male parents were sown on each mulched strip at a sowing rate of 9.34 × 10⁴ female plants ha⁻¹ with 2.34 × 10⁴ male plants ha⁻¹. In 2016, every seven rows of the female parents with one row of the male parents were planted on two adjacent mulched strips at a sowing rate of 8.17 × 10⁴ female plants ha⁻¹ with 1.17 × 10⁴ male plants ha^-1. Each mulched strip was controlled by two laterals of the drip irrigation system (the inner flat emitter drip tape was used with an emitter spacing of 30 cm and a flow rate of 3 L·h⁻¹) in 2016. The average spacing within rows was 25.8 cm in both years. Female parents were sown on April 13–14, 2015 and April 23–24, 2016. The pollination period was extended by sowing the male parents twice, about 5 and 10 days later than the sowing date of the female parents. Female parents were detasseling on July 7 and harvested on September 16–18 in both years.

Basal fertilizers were spread on the soil surface and plowed evenly into the soil before soil mulching with plastic film. In 2015 the amounts of basal fertilizers were 216 kg·ha⁻¹ N, 105 kg·ha⁻¹ P, 62 kg·ha⁻¹ K and 37.5 kg·ha⁻¹ magnesium-zinc sulfate (Zn+ MgO ≥ 17%). A further 207 kg·ha⁻¹ N as urea was top-dressed and then dissolved and infiltrated into the soil with irrigation water at the first irrigation event. Also, foliar spraying with KH₂PO₄ was conducted at the jointing stage. In 2016, only 205 kg·ha⁻¹ N and 101 kg·ha⁻¹ P were applied as basal fertilizers, and two top-dressings were fertigated with 139 kg·ha⁻¹ N on June 23–24 and 70 kg·ha^-1 N on July 15–16.

In 2015, border irrigation (irrigation controlled or directed by short dikes around the areas treated) was used. Five irrigation events were used with the irrigation units marked with a red dashed line (Fig. 1(a)). Outside the red dashed line, only the first four irrigation events were used. Also, winter irrigation with more than 150 mm was applied in the last winter to control salinity and provide adequate soil moisture for seedling emergence in 2015. Winter irrigation in 2015 was canceled because drip irrigation (surface drip irrigation) was planned in 2016. Instead, irrigation was applied with a drip irrigation system just after the female parents were sown to assure seedling emergence in 2016. Also, another seven irrigation events were applied during the growing season. The total amount of irrigation was 679 mm (winter irrigation was not included) under border irrigation in 2015 and 452 mm under drip irrigation in 2016. Irrigation schedules are listed in Table 1.

Fig.1 Sampling location and elevation map of the study area. (a) Border irrigation in 2015; (b) drip irrigation in 2016. The area with red dashed line demonstrates the irrigation units with five irrigation events; the other irrigation units applied only the first four irrigation events. The red stars mark the location of the eddy covariance system.

Full size|PPT slide

Tab.1 Irrigation schedule summary

No.	Border irrigation in 2015		Drip irrigation in 2016
	Date	Irrigation amount (mm)	Date	Irrigation amount (mm)
0	Last winter	>150	–	–
1	May 30–June 1	158	April 24–25	50
2	Jun 29–July 2	158	June 10–12	51
3	July 20, 22 and 25	150	June 22–24	51
4	August 9 and 12	158	July 2–3	50
5	September 5–6	95	July 15–16	66
6	–	–	July 27–29	55
7	–	–	August 6–7	68
8	–	–	August 21–22	62
Total		679		452

Note: Irrigation amount is the area-weighted average for the irrigated units. The total is the area-weighted average of total amount of irrigation during the growing season, and winter irrigation is omitted.

2.2 Data collection and preprocessing

In 2015, 185 sampling locations were used, comprising 169 sampling points on an approximate 15-by-15-m grid and another 16 sampling points placed randomly with the distance to the grid points less than 15 m. In 2016, 184 sampling points were used in the same field, comprising 156 sampling points on an approximate 15-by-15-m grid and another 28 sampling points placed randomly with the distance to the grid points less than 15 m. A further 18 sampling points, labeled with ID numbers 185–202 (Fig. 1(b)), were added on an adjacent field under the same agronomic management in 2016. These additional sampling points were not included in the statistical analysis as there are no corresponding samples from 2015; however, they were used in searching the maximum yield, clustering soil texture zones and fitting boundary lines. Sampling points were adjusted slightly when too close to the field ridge. The coordinates of each sampling point were obtained using a GPS tracker (Trimble Recon, Sunnyvale, CA) and then projected into a Gauss-Kruger projected coordinate system (Beijing 1954, 3 Degree GK CM 102E), and the elevation was obtained with an optical level (Fig. 1).

Each sampling site was a 3.3 m (two mulched strips) by 3.1 m (12 plants per row) area. Emergence was counted to calculate plant density. Averaged emergence was 74% under border irrigation and 81% under drip irrigation. All maize ears were harvested in the representative area then weighed for fresh mass (M_fresh), ten ears were oven-dried at 80°C to determine ear water content (q', according to Eq. (1)) and the mass ratio of grain to ear (p) was calculated. The grain yield was calculated on dry matter basis (Eq. (2)). Then the relative yield (RY, expressed as actual crop yield divided by the observed maximum yield in each year. Equation (3) was used in the analysis to avoid the influences of meteorological and cultivar differences between the two production years.

(1)$\mathrm{\theta }\text{'}=({\mathrm{m}}_{\mathrm{fresh}}-{\mathrm{m}}_{\mathrm{dry}})/{\mathrm{m}}_{\mathrm{dry}}$

(2)$\mathrm{Yield}={\displaystyle \frac{{\mathrm{M}}_{\mathrm{dry}}\times p}{A}}\times \mathrm{10}={\displaystyle \frac{p\times {\mathrm{M}}_{\mathrm{fresh}}/(\mathrm{1}+{\theta }^{\prime })}{A}}\times \mathrm{10}$

(3)$\mathrm{RY}=\mathrm{Yield}/{\mathrm{Yield}}_{\max }$

where q' is the ear water content. m_fresh and m_dry are the fresh and dry mass of ten ears, respectively, g. Yield is the grain yield in dry mass, Mg·ha^-1. M_dry and M_fresh are the dry and fresh mass of all the ears harvested on the representative area, kg. p is the mass ratio of grain to ear. A is the harvested area, m². 10 is the coefficient converting grain yield from kg·m⁻² to Mg·ha^-1. RY is the relative yield. Yield_max is the maximum of all the observed grain yield in each year.

Soil samples were collected four times each year to a depth of 100 cm at 20 cm depth intervals in the representative area at the start and end of the growing season, and on June 8 and July 10, 2015 (about one week after irrigation), and June 30 and August 6, 2016 (one day before irrigation). Four soil samples at the same soil depth were mixed and processed to determine soil properties. All samples were used to determine soil gravimetric water content (q') and nitrate- and ammonium-nitrogen (NO₃^–-N and NH₄⁺-N). All samples collected at the start of each growing season were used to determine soil texture (clay, silt and sand fractions) and the volume-weighted mean diameter (D_vol). Also, samples from 0 to 40 cm depth at the start and the end of each growing season were used to determine pH in water, electrical conductivity (EC), organic matter (OM) content, total nitrogen (TN), total phosphorus (TP), soil available phosphorus (AP), and soil available potassium (AK) contents. Bulk density (BD) in the top 100 cm of the soil profile was determined on the 130 subsample sites with drip irrigation. These sampling sites included ID numbers 1–155 (odd numbers only), all ID numbers 157–202, and other ID numbers 2, 32, 34, 48, 76 and 110.

Soil pH was determined in a 1:2.5 soil-water suspension with a pH meter after shaking for 2 min and standing for 30 min. Soil electrical conductivity was determined in a 1:5 soil-water extract using a conductivity meter after shaking for 3 min. Soil OM was determined by the Walkley-Black rapid dichromate oxidation method^[23]. TN was determined by Kjeldahl digestion and continuous flow analysis (AutoAnalyzer 3, Bran+ Luebbe, SEAL Analytical GmbH, Norderstedt, Germany). TP was determined by digesting soil samples with an acid mixture (H₂SO₄-HClO₄) and continuous flow analysis. NO₃^–-N and NH₄⁺-N were determined by extracting samples with 1 mol·L⁻¹ KCL solution and continuous flow analysis. AP was determined by extracting samples with 0.5 mol·L⁻¹ NaHCO₃ solution and the molybdate colorimetric method. AK was determined by extracting samples with 1 mol·L⁻¹ NH₄OAc solution and flame photometry. Soil gravimetric water content (q') and bulk density (BD) were determined by the oven drying method. Soil texture was determined with a laser particle size analyzer (Mastersizer2000, Malvern Instruments Ltd., Malvern, UK).

Soil available nitrogen (AN_{0–100 cm}) is expressed as the average of available nitrogen (nitrate- and ammonium-nitrogen summed) in top 100 cm of the soil profile. Other soil properties (pH, EC, OM, TN, TP, AP, AK) were the averages in the top 40 cm of the soil profile.

Soil available water content in the top 100 cm of the soil profile (PropAWC) is expressed as the fraction of the total available water content (the moisture held between field capacity and permanent wilting point). PropAWC was calculated using Eq. (4) to Eq. (8). The permanent wilting point (q_WP) was estimated using the pedotransfer functions (PTFs) described by van den Berg et al.^[24] (Eq. (6)). Field capacity (q_FC) and BD were estimated using the PTFs when no observed data were available (Eq. (7) and Eq. (8)).

(4)$\mathrm{PropAWC}={\displaystyle \frac{{\mathrm{SWC}}_{\mathrm{Actual}}-{\mathrm{SWC}}_{\mathrm{WP}}}{{\mathrm{SWC}}_{\mathrm{FC}}-{\mathrm{SWC}}_{\mathrm{WP}}}}={\displaystyle \frac{{\displaystyle {\sum }_{i=\mathrm{1}}^{\mathrm{5}}{\mathrm{200}\theta }_{\mathrm{Actual}[i]}-{\displaystyle {\sum }_{i=\mathrm{1}}^{\mathrm{5}}{\mathrm{200}\theta }_{\mathrm{WP}[i]}}}}{{\displaystyle {\sum }_{i=\mathrm{1}}^{\mathrm{5}}{\mathrm{200}\theta }_{\mathrm{FC}[i]}}-{\displaystyle {\sum }_{i=\mathrm{1}}^{\mathrm{5}}{\mathrm{200}\theta }_{\mathrm{WP}[i]}}}},(i=\mathrm{1},\mathrm{2},\dots ,\mathrm{5})$

(5)${\mathrm{\theta }}_{\mathrm{Actual}}=\mathrm{\theta }{\text{'}}_{\mathrm{Actual}}\times \mathrm{BD}$

(6)${\mathrm{\theta }}_{\mathrm{WP}}=\mathrm{}\mathrm{0.00}\mathrm{334}\times \mathrm{Clay}\times \mathrm{BD}+\mathrm{}\mathrm{0.00}\mathrm{1}\mathrm{0}\mathrm{4}\times \mathrm{Silt}\times \mathrm{BD}$

(7)${\mathrm{\theta }}_{\mathrm{FC}}=\mathrm{pr1}+\mathrm{pr2}\times \mathrm{Clay}+\mathrm{pr3}\times {\mathrm{Clay}}^{\mathrm{2}}+\mathrm{pr4}\times \mathrm{Silt}+\mathrm{pr5}\times {\mathrm{Silt}}^{\mathrm{2}}+\mathrm{pr6}\times \mathrm{Clay}\times \mathrm{Silt}$

(8)$\mathrm{BD}=\mathrm{pr1}+\mathrm{pr2}\times \mathrm{Clay}+\mathrm{pr3}\times {\mathrm{Clay}}^{\mathrm{2}}+\mathrm{pr4}\times \mathrm{Silt}+\mathrm{pr5}\times {\mathrm{Silt}}^{\mathrm{2}}+\mathrm{pr6}\times \mathrm{Clay}\times \mathrm{Silt}$

where SWC_Actual, SWC_WP, and SWC_FC are the soil water storage in the top 100 cm of the soil profile when soil was sampled, soil potential reached -1.5 MPa, and field capacity reached, respectively, mm. i (i = 1, 2, ⋯, 5) is the number of the soil layer 0–20, 20–40, 40–60, 60–80 and 80–100 cm. q'_Actual and q_Actual are the soil gravimetric and volumetric water contents when the soil was sampled. q_FC and q_WP are the volumetric water contents when field capacity reached, and soil potential reached -1.5 MPa, respectively. Clay and silt are the fractions of soil particle size ˂ 2 mm and between 2 and 50 mm (when clay was 30%, then 30 would be used in Eq. (6) to Eq. (8)), respectively. BD is the soil bulk density, Mg·m⁻³. pr1 to pr6 are the parameters of the PTFs estimated by nonlinear regression with XLSTAT version 2014.5.03 (Table 2). When estimating the parameters of q_FC (Eq. (7)), data was collected from other research at the same study area (sample size: n = 26; the measured q_FC ranged from 18.4% to 38.2%). For the parameters of BD at different soil depths (Eq. (8)), both the measured data set in this study (bulk density of the 0–100 cm depth soil profile with 20 cm depth intervals from 130 sampling sites) and collected data set from other studies (5 sampling sites) were used (sample size: n = 135).

Tab.2 Nonlinear regression parameters qFC and BD

Parameter	θFC	BD0–20cm	BD20–40cm	BD40–60cm	BD60–80cm	BD80–100cm
pr1	2.281E-01	1.568E+00	1.642E+00	1.460E+00	1.600E+00	1.555E+00
pr2	–2.288E-02	6.582E-03	6.495E-03	–7.208E-03	–2.454E-02	3.084E-02
pr3	1.202E-03	3.729E-04	–1.786E-04	1.821E-03	2.280E-03	–1.260E-05
pr4	3.103E-03	–1.175E-03	–4.515E-04	9.433E-03	4.317E-03	–3.196E-03
pr5	1.727E-05	–2.200E-05	–2.170E-05	–1.049E-04	–2.940E-05	4.800E-05
pr6	–1.410E-04	–1.533E-04	–1.654E-04	–4.962E-04	–4.701E-04	–6.000E-04

Note: The parameters were estimated using the nonlinear regression method in XLSTAT version 2014.5.03. For qFC, data were collected from other studies in the same study area (sample size, n = 26; measured qFC ranged from 18.4% to 38.2%). For BD at different soil depths, both the measured data set from the present study (bulk density of 0–100 cm depth of soil profile with 20–cm depth intervals from 130 sampling sites) and a data set from other studies (5 sampling sites) were used (sample size, n = 135).

2.3 Data analysis

2.3.1 Boundary line analysis and quantification of yield gaps

The yield gap is commonly defined as the difference between yield potential and actual yield achieved in the field. Yield potential, also called potential yield, is defined as the yield of a cultivar when grown in an environment to which it is adapted, with nutrients and water non-limiting and other stresses effectively controlled^[8]. However, it is difficult to measure the yield potential as it can be difficult to provide an entirely stress-free testing environment. Yield potential is therefore usually estimated by the maximum yield of sizable samples under optimum management. As the main concern of this study was closing yield gap as determined by optimizing field management, and RY (each season) was used to avoid the influence of meteorological difference between the two growing seasons, we define the yield gap (D_Total) as the difference between the maximum RY (RY_max) and the RY achieved:

(9)${\mathrm{\Delta }}_{\mathrm{Total}}={\mathrm{RY}}_{\max }-\mathrm{RY}$

BLA has been described in previous studies^[18–20]. The boundary line is usually determined by fitting a line through the upper boundary points of data scatter (yield as the dependent variable and yield-limiting factors as the independent variables). It is assumed that the boundary line represents the dependent variable influenced only by the single independent variable in the absence of any other limiting factors if the data set is sufficient^[21]. Here, a large area was affected by plant density deficiency. Thus the boundary lines of other yield-limiting factors (except for plant density) were inevitably influenced by two or more yield-limiting factors in most cases (at least one was the independent variable, and another was plant density). This conflicts with the assumption of the BLA. We therefore gave priority to determining the boundary line of plant density. We used RY as the dependent variable for the boundary line of plant density, but adjusted RY as the dependent variable for boundary lines of other yield-limiting factors by eliminating the reduction effect of all the yield-limiting factor(s) that were analyzed (Fig. 2). The procedure is described below.

Fig.2 Boundary line analysis (BLA) of different yield-limiting factors. (a) Plant density of female parents, 104 plants ha−1; (b) soil available potassium at 0–40 cm depth at the end of the growing season, mg·kg−1 K; (c) soil available nitrogen at 0–100 cm depth at the end of the growing season, mg·kg−1 N; (d) soil available water content at 0–100 cm depth on June 8 under border irrigation in 2015; (e,f) soil available water content at 0-100 cm depth on 6 August and 13 September 2016 under drip irrigation; (g) soil electrical conductivity at 0–40 cm depth when maize was sown, mS·cm−1; (h) soil pH in water at 0–40 cm depth when maize was sown. White circles (o) represent data serial in 2015; black circles (•) represent data serial in 2016; red crosses (×) represent the upper points eliminated as outliers; black crosses ( + ) represent the upper points used in fitting the upper boundary line.

Full size|PPT slide

Supposing crop yield (we used RY) was the outcome that potential yield (as the potential yield was difficult to measure, we used the observed maximum yield instead, RY_max = 1) under specific climatic conditions constrained by yield-limiting factors (X_i, i = 1, 2, ..., n) and unknown factors (e, including the reduction effect of others unknown yield-limiting factors, random error and interaction effects between factors X_i if they exist). The reduction effect of each yield-limiting factor was expressed as a yield correction coefficient, which was a function of each factor f(X₁), f(X₂), ..., f(X_n), f(e). This may be expressed as follows.

(10)$\mathrm{RY}={\mathrm{RY}}_{\max }\times \mathrm{f}({\mathrm{X}}_{\mathrm{1}})\mathrm{}\times \mathrm{f}({\mathrm{X}}_{\mathrm{2}})\mathrm{}\times \cdots \times \mathrm{f}({\mathrm{X}}_n)\mathrm{}\times \mathrm{f}(\mathrm{\epsilon })$

As the assumption of BLA, the boundary line for the factor X₁ represents the dependent variable influenced only by a single independent variable X₁ in the absence of any other limiting factors when the data are sufficient. Then Eq. (10) can be simplified as follows, as all the reduction coefficient function of factors (X₂, X₃, ..., X_n, e) are equal to 1:

(11)$\mathrm{BLA}({\mathrm{X}}_{\mathrm{1}})\mathrm{}={\mathrm{RY}}_{\max }\times \mathrm{f}({\mathrm{X}}_{\mathrm{1}}),\mathrm{i}.\mathrm{e}.,\mathrm{f}({\mathrm{X}}_{\mathrm{1}})\mathrm{}=\mathrm{BLA}({\mathrm{X}}_{\mathrm{1}})/{\mathrm{RY}}_{\max }$

where RY_max is the maximum RY, BLA(X₁) is the boundary line function of factor X₁. According to Eq. (10) and Eq. (11), the reduction effect of other factors (X₂, X₃, ..., X_n, e) can be expressed as:

(12)${\mathrm{f(X}}_{\mathrm{2}})\times \dots \times {\mathrm{f(X}}_n)\times \mathrm{f(\epsilon })={\displaystyle \frac{\mathrm{RY}}{{\mathrm{RY}}_{\max }\times {\mathrm{f(X}}_{\mathrm{1}})}}\stackrel{\mathrm{Eq}\mathrm{.11}}{\iff }{\displaystyle \frac{\mathrm{RY}}{{\mathrm{BLA(X}}_{\mathrm{1}})}}$

When factor X₁ is managed to the optimum level, the reduction effect of X₁ can be eliminated, i.e., the yield correction coefficient function f(X'₁) ≡1 (the symbol f(X'_i) is used to represent the reduction function after optimizing factor X_i distinguished from the yield correction coefficient function f(X_i) before optimizing factor X_i; RY'_Xi is used to represent the predicted RY when the first i yield-limiting factor(s) X₁, X₂, ..., X_i is/are managed to the optimum level(s). Then the predicted RY (${{\mathrm{R}\mathrm{Y}}^{\prime }\mathrm{}}_{{\mathrm{X}}_{\mathrm{1}}}$) after optimizing factor X₁ can be expressed as:

(13)$\mathrm{R}{{\mathrm{Y}}^{\prime }}_{{\mathrm{X}}_{\mathrm{1}}}={\mathrm{RY}}_{\max }\times {\mathrm{f(}{\mathrm{X}}^{\prime }}_{\mathrm{1}})\times {\mathrm{f(X}}_{\mathrm{2}})\times \dots \times {\mathrm{f(X}}_n)\times \mathrm{f(\epsilon })\stackrel{{\mathrm{f(}{\mathrm{X}}^{\prime }}_{\mathrm{1}})\equiv \mathrm{1}}{\iff }{\mathrm{RY}}_{\max }\times {\mathrm{f(X}}_{\mathrm{2}})\times \dots \times {\mathrm{f(X}}_n)\times \mathrm{f(\epsilon })={\mathrm{RY}}_{\max }\times {\displaystyle \frac{\mathrm{RY}}{{\mathrm{BLA(X}}_{\mathrm{1}})}}$

The difference ($\mathrm{\Delta }{\mathrm{RY}}_{{\mathrm{X}}_{\mathrm{1}}}$) between RY and ${{\mathrm{R}\mathrm{Y}}^{\prime }\mathrm{}}_{{\mathrm{X}}_{\mathrm{1}}}$ can be explained as the yield gap determined by single factor X₁ when factor X₁ was managed to the optimum level:

(14)$\mathrm{\Delta }{\mathrm{RY}}_{{\mathrm{X}}_{\mathrm{1}}}={\mathrm{R}{\mathrm{Y}}^{\prime }}_{{\mathrm{X}}_{\mathrm{1}}}-\mathrm{RY}={\mathrm{RY}}_{\max }\times {\displaystyle \frac{\mathrm{RY}}{{\mathrm{BLA(X}}_{\mathrm{1}})}}-\mathrm{RY}$

Then the ${{\mathrm{R}\mathrm{Y}}^{\prime }\mathrm{}}_{{\mathrm{X}}_{\mathrm{1}}}$ can be used as the new dependent variable to determine the boundary line of the following factor X₂ free from the yield reduction effect of the previous factor X₁. Similarly, the boundary lines of other factors can be determined in sequence. The general expressions for factors X_i (i= 2, 3, ..., n) are as follows:

(15)${\mathrm{R}{\mathrm{Y}}^{\prime }}_{{\mathrm{X}}_i}{\mathrm{=RY}}_{\max }\times {\displaystyle \frac{{\mathrm{R}{\mathrm{Y}}^{\prime }}_{{\mathrm{X}}_{i\mathrm{-1}}}}{{\mathrm{BLA(X}}_i)}},i=\mathrm{2},\mathrm{3},\mathrm{}\mathrm{...},n$

(16)$\mathrm{\Delta }{\mathrm{RY}}_{{\mathrm{X}}_i}={\mathrm{R}{\mathrm{Y}}^{\prime }}_{{\mathrm{X}}_i}-{\mathrm{R}{\mathrm{Y}}^{\prime }}_{{\mathrm{X}}_{i\mathrm{-1}}},i=\mathrm{2},\mathrm{3},\mathrm{}\mathrm{...},n$

where ${{\mathrm{R}\mathrm{Y}}^{\prime }\mathrm{}}_{{\mathrm{X}}_{i-\mathrm{1}}}$ is the predicted RY with the reduction effect of previous factors (X₁, X₂, ..., X_i–1) eliminated by managing them to the optimum levels, which was also the dependent variable when fitting the boundary line of factors X_i; ${{\mathrm{R}\mathrm{Y}}^{\prime }\mathrm{}}_{{\mathrm{X}}_i}$is the predicted RY with the reduction effect of the first i factors (X₁, X₂, ..., X_i–1, X_i) eliminated by managing them to the optimum levels; BLA(X_i) s the boundary line function of factor X_i; $\mathrm{\Delta }{\mathrm{RY}}_{{\mathrm{X}}_i}$ is the yield gap determined by a single factor X_i.

When all the boundary lines of the yield-limiting factors considered had been determined, the yield gap of unknown factors e was estimated by Eq. (17):

(17)$\mathrm{\Delta }{\mathrm{RY}}_{\mathrm{\epsilon }}={\mathrm{RY}}_{\max }-(\mathrm{RY}+\mathrm{\Delta }{\mathrm{RY}}_{{\mathrm{X}}_{\mathrm{1}}}+\mathrm{\Delta }{\mathrm{RY}}_{{\mathrm{X}}_{\mathrm{2}}}+\dots +\mathrm{\Delta }{\mathrm{RY}}_{{\mathrm{X}}_n})$

where RY was the relative yield achieved in the field, RY_max was the maximum RY; $\mathrm{\Delta }{\mathrm{RY}}_{{\mathrm{X}}_{\mathrm{1}}}$, $\mathrm{\Delta }{\mathrm{RY}}_{{\mathrm{X}}_{\mathrm{2}}}$, ..., $\mathrm{\Delta }{\mathrm{RY}}_{{\mathrm{X}}_n}$ and $\mathrm{\Delta }{\mathrm{RY}}_{\mathrm{\epsilon }}$ were the yield gaps determined by single yield-limiting factors X_i (i = 1, 2, ..., n) and unknown factors (e), respectively.

The main steps for quantifying the yield gaps by multiple yield-limiting factors are as follows:

(1) Determining the sequence of fitting the boundary lines for different yield-limiting factors. The factor that was more likely to be the dominant yield-limiting factor and with less uncertainty should be given priority (plant density in our study). However, the order should be regulated if the upper boundary points, when fitting the boundary line of the current factor, are commonly outside the optimum levels of the latter factors (RY is constrained by the latter factors).

(2) Fitting the boundary line of every single factor (X_i). We used yield-limiting factor X_i as the independent variable and the RY as the dependent variable for the first factor X₁, but we adjusted relative yield (${{\mathrm{R}\mathrm{Y}}^{\prime }\mathrm{}}_{{\mathrm{X}}_{\mathrm{1}}}$) as the dependent variable for the subsequent factors X_{i + 1} (i= 1, 2, ..., n). For the convenience of description, the dependent variable is expressed as Y_i. First, searching the point when the maximum y-value is reached (if two or more points are identified, using the point with the largest x-value), expressed as (X_i,k, Y_i,max). Secondly, sorting (X_j, Y_j) by the independent variables X_i in ascending order. For (X_j, Y_j) with the same x-value, using the dependent variable as the second sorting variable in descending order when X_i,j<X_i,k (j = 1, 2, ..., k, ..., m), and in ascending order when X_i,j>X_i,k. Thirdly, searching the upper boundary points. The three points (X_i,1, Y_i,1), (X_i,m, Y_i,m), and (X_i,k, Y_i,max) are saved as the upper boundary points. Using (X_i,1, Y_i,1) as the current identified upper boundary point (X_i,temp = X_i,1, Y_i,temp = Y_i,1) to search other points from X_i,1 to X_i,k (searching from the left side). Another point (X_i,j, Y_i,j) would be identified as the upper boundary point if Y_i,j>Y_i,temp, and it is used as the new current identified upper boundary point for further searching. Similarly, using the (X_i,m, Y_i,m) as the current identified upper boundary point, searching other points from X_i,m to X_i,k (searching from the right side). Fourthly, excluding outliers of the upper boundary points data set. Points deviating from the main trends are treated as outliers. Also, if a point of the x-value of any yield-limiting factors (X_{i + 1}, ..., X_n) (after the current analysis factor X_i) is not within its optimum range, it is also deleted as an outlier. Finally, fitting an appropriate model through the upper boundary points for every single factor. More details about the searching procedure are given by Schnug et al.^[19] and we conducted the procedure using the R version 3.5.2 statistical package^[25].

(3) Calculating the yield gaps for different yield-limiting factors (X₁, X₂, ..., X_n) (Eq. (16)) and unknown factors (e) (Eq. (17)).

To describe the reduction intensity of different yield-limiting factors we defined the optimum range and acceptable range. The optimum range is defined as the independent variable interval with RY of the boundary line (maximum attainable RY) equal to 100%, and the acceptable range as the independent variable interval with RY of the boundary line ˃ 95%. The independent variable between the acceptable range was supposed to be acceptable stress and outside the acceptable range was regarded as non-ignorable stress.

The results of the BLA are shown in Fig. 2 and Table 3. Soil properties, comprising soil organic matter, total nitrogen, total phosphorus and available phosphorus, were not included in the BLA, as the boundary lines of these properties were approaching the horizontal line RY= 1.

Tab.3 Parameters of the BLA method for different yield-limiting factors and areal percentage under different yield reduction levels

ID	Factor		Optimum range		Acceptable range		Function form				Parameters (significance codes and sample size)
I	Num (104 plants ha−1)		7.61–7.84		6.57–8.88		min (y=ax2+bx+c, 1)				a = –3.78E-02 (***); b = 0.584 (***); c = –1.26 (***); (F-test: ***, n = 5)
II	AK (mg·kg−1 K) §		>80.5		>70.8		y = a(x–t1)+1; x<t1 y = 1; x>t1				a = 5.17E-03 (***); t1 = 80.5 (※); (n = 6)
III	AN (mg·kg−1 N) §		6.04–26.1		3.14–66.2		y = a(x–t1)+1; x<t1 y = 1; t1<x<t2 y = b(x–t2)+1; x>t2				a = 1.73E−02 (***); b= –1.25E−03 (***); t1 = 6.04 (※); t2 = 26.1 (※); (n = 7)
IV-1	PropAWC 20150608 ¶		>0.739		>0.679		y = a(x–t1)+1; x< t1 y = 1; x> t1				a = 0.835 (***); t1 = 0.739 (※); (n = 5)
IV-2	PropAWC 20160806 ¶		>0.371		>0.282		y = a(x–t1)+1; x< t1 y = 1; x>t1				a = 0.559 (**); t1 = 0.371 (**); (n = 3)
V	PropAWC 20160913 ¶		>0.507		>0.422		y = a(x–t1)+1; x<t1 y = 1; x> t1				a = 0.583 (***); t1 = 0.507 (***); (n = 6)
VI	EC (μS·cm−1)		<224		<323		y = 1; x<t2 y = b(x–t2)+1; x>t2				b = –5.06E-04 (***); t2 = 224 (※); (n = 9)
VII	pH		<8.41		<8.49		y = 1; t1<x<t2 y = b(x–t2)+1; x> t2				b = –0.616 (***); t2 = 8.41 (***); (n = 14)
ID	Percentage to the whole area under different reduction level (border irrigation in 2015) (%)							Percentage to the whole area under different reduction level (drip irrigation in 2016) (%)
	Deficiency			Optimum	Excessive			Deficiency		Optimum		Excessive
	<0.95†	[0.95, 1]		1	[0.95, 1]	<0.95		<0.95	[0.95, 1]	1		[0.95, 1]	<0.95
I	38	43		6	12	1		44	56	0		0	0
II	2	6		92	–	–		3	9	87		–	–
III	0	18		68	12	2		0	0	56		39	5
IV-1	5	5		89	–	–		–	–	–		–	–
IV-2	–	–		–	–	–		11	16	73		–	–
V	–	–		–	–	–		2	7	91		–	–
VI	–	–		86	12	2		–	–	78		19	3
VII	–	–		87	12	1		–	–	100		0	0

Note: ID, order of fitting the boundary line; Num, female plant density; AK, soil available potassium (at 0–40 cm) at the end of the growing season; AN, soil available nitrogen (at 0–100 cm) at the end of the growing season; PropAWC, soil available water content (in the fraction of total available water content); EC, soil electrical conductivity at the start of the growing season; pH, soil pH at the start of the growing season; optimum range, the independent variable interval that relative yield of the boundary line equaled 100%; acceptable range, the independent variable interval that relative yield of the boundary line was ˃ 95%. When determining the parameters of boundary lines with R version 3.5.2[25], lm() function was used for the factor Num (ID: I), and both the t-test and F-test were conducted; nls() function was used for other factors (ID: II-VII), and only the t-test was conducted. Significance of t-test and F-test: ***, P<0.001; **, P<0.01; *, P<0.05. ※ The value of t1 or t2 was set manually to the minimum or maximum of the independent variable with the adjusted RY approaching 1.000 (>0.999) if the fitting program failed for too many parameters or too small a sample size. ¶ Soil sampling dates were June 8, 2015, August 6, 2016 and September 13, 2016. § Nutrient convertion coefficient: 1 mg·kg−1 = 6.02 kg·ha−1 (at 0–40 cm soil), 1 mg·kg−1 = 14.1 kg·ha−1 (at 0–100 cm soil); † Yield reduction level is expressed as the maximum attainable relative yield (relative yield of the boundary line) that was reached, and 0.95 or 1 represents the relative yield of the boundary line that reached 0.95 or 1. This is consistent with the lower or upper limits of the acceptable range and optimum range.

2.3.2 Hierarchical cluster analysis

To analyze the relationship among soil moisture, nutrients and crop yield at the same soil zones, hierarchical cluster analysis was conducted using XLSTAT version 2014.5.03. Soil D_vol in the top 100 cm of the soil profile (20-cm-depth intervals) was selected as the clustering variable. Ward's method was chosen as the agglomeration method, and Euclidian distance was used to measure the dissimilarity coefficient. As a result, the whole field was divided into four soil texture zones (Fig. 3). Also, the box plots of the remaining soil NPK in different soil texture zones are shown in Fig. 4.

Fig.3 Hierarchical cluster analysis results. (a) Dendrogram of the hierarchical cluster analysis; (b) the volume-weighted mean diameter of the class centroids in the soil profile at 0–100 cm depth (20 cm depth intervals) in each soil zone; (c) box plots of the volume-weighted mean diameter at 0–60 cm depth (orange lines) and 60–100 cm depth (purple lines) in each soil zone.

Full size|PPT slide

Fig.4 Box plots of remaining soil available nitrogen (at 0–100 cm depth), available phosphorus (0–40 cm depth) and available potassium (0–40 cm) at the end of the growing seasons. (a) Soil available nitrogen; (b) soil available phosphorus; (c) soil available potassium. The green line represents the 2015 results and the black line the 2016 result. Significance by independent-samples t-test using the SPSS 20.0 software package: ***, P<0.001; **, P<0.01; *, P<0.05; n.s., not significant (P>0.05).

Full size|PPT slide

2.3.3 Interpolation method

The inverse distance weighted method was used for interpolation when plotting the spatial distribution of yield and yield-limiting factors^[26].

(18)$Z(x_{\mathrm{0}})=\left[{\displaystyle {\sum }_{i=\mathrm{1}}^m{\displaystyle \frac{\mathrm{1}}{d_i^p}}Z(x_i)}\right]/\left[{\displaystyle {\sum }_{i=\mathrm{1}}^m\frac{\mathrm{1}}{d_i^p}}\right]$

where $Z(x_{\mathrm{0}})$ is the interpolated value at point $x_{\mathrm{0}}$, $Z(x_i)$ is the measured value at sampling points $x_i$ (i = 1, 2, ..., m), d_i is the distance between $x_{\mathrm{0}}$ and $x_i$, p (ranging from 1 to 100) is the power value of the distance d_i. For the parameters settings when interpolating RY and yield-limiting factors, both the searching major and minor semiaxes were 15 m, minimum and maximum neighbors of each sector were 1 and 2, four-sector with 45° offset was selected, and an optimum power value p was used according to cross-validation by minimizing the root mean square prediction error. The value of p was 1.75, 1.96, 3.94, 1.58, 1.00, 1.00, 1.00 and 1.00 for RY, RY', Num, AK, AN, PropAWC_20150608, pH and EC in 2015, respectively, and p was 1.00, 1.00, 1.00, 1.00, 1.24, 1.00, 1.00, 1.00 and 1.40 for RY, RY', Num, AK, AN, PropAWC_20160806, PropAWC_20160913, pH and EC in 2016. However, to maintain the soil texture of the predicted point consistent with the nearest observed point when interpolating soil texture zone, one sector and one neighbor for each sector were selected, and p could be any value. Interpolation was conducted using ArcGIS 10.4 and the interception plots are shown in Fig. 5 and Fig. 6.

Fig.5 Spatial interpolation plots of remaining soil available nitrogen (at 0–100 cm depth), available phosphorus (0–40 cm depth) and available potassium (0–40 cm) at the end of the growing seasons. (a–c) Soil available nitrogen, available phosphorus and available potassium in 2015; (d–f) soil available nitrogen, available phosphorus and available potassium in 2016. Different soil zones are also shown.

Full size|PPT slide

Fig.6 Interpolation maps of the relative yield, locations of different soil zones and locations where yield was affected by different yield-limiting factors. (a,f) Relative yield in 2015 (a) and 2016 (f); (b,g) locations where plant density deficit or available potassium deficit occurred in 2015 (b) and 2016 (g); (c,h) locations where water deficit, available nitrogen deficit or surplus occurred in 2015 (c) and 2016 (h); (d,i) locations where high soil salinity or high pH occurred in 2015 (d) and 2016 (i); (e,g) predicted relative yield when all the considered yield-limiting factors (not including unknown factors) were at the optimum ranges in 2015 (e) and 2016 (g). The relative yield map in 2015 is also shown in (b–d) and the relative yield map in 2016 is shown in (g–i). Soil texture zones are also depicted in (a), (e), (f), and (j).

Full size|PPT slide

3 RESULTS

3.1 Hierarchical cluster analysis and zonal statistics of soil NPK

The hierarchical cluster analysis results are shown in Fig. 3. The whole field was divided into four soil zones. The soil texture according to the USDA soil texture classification system of the cluster centroid in each soil zone is listed in Table 4. The area percentages of Zones 1, 2, 3 and 4 were 41%, 30%, 16% and 13%, respectively.

Tab.4 Summary of the soil texture zones and residual soil NPK

Zone (year)	Soil texture (0–60 cm/60–100 cm)	Percentage of the whole area (samples count) (%)	AN0–100 cm (mean±SD) (mg·kg−1 N)	AP0–40 cm (mean±SD) (mg·kg−1 P)	AK0–40 cm (mean±SD) (mg·kg−1 K)
All Zones (2015)	–	100 (185)	16.3±17.3	12.1±3.3	104.3±21.2
Zone 1 (2015)	Silt loam/Silt loam	41 (71)	23.1±21.5 b¶	13.2±3.3 c¶	116.8±22.6 a¶
Zone 2 (2015)	Loam/Silt loam	30 (53)	14.1±14.0 cd	11.6±2.9 d	103.5±13.5 b
Zone 3 (2015)	Sandy loam/Silt loam	16 (36)	12.0±11.4 cd	10.2±3.3 e	87.5±15.7d
Zone 4 (2015)	Sandy loam/Loamy sand	13 (25)	8.1±10.0 d	13.0±2.9 cd	94.7±14.7 cd
All Zones (2016)	–	100 (184)	30.4±25.0	14.8±3.4	101.9±22.6
Zone 1 (2016)	Silt loam/Silt loam	41 (82)	39.3±31.1 a	14.9±3.7 ab	112.9±22.4 a
Zone 2 (2016)	Loam/Silt loam	30 (64)	25.2±18.0 b	15.0±3.2 ab	95.4±19.1 c
Zone 3 (2016)	Sandy loam/Silt loam	16 (23)	21.4±9.4 b	13.5±2.9 bc	87.8±17.6 cd
Zone 4 (2016)	Sandy loam/Loamy sand	13 (15)	17.8±8.3 bc	15.8±2.6 a	90.5±17.2 cd

Note: Soil texture, the soil texture of the class centroid, calculated by averaging the sand, silt and clay contents at 0–60 cm or 60–100 cm depth in each soil zone, according to the USDA soil texture classification system; AN0–100 cm, AP0–40 cm and AK0–40 cm, the remaining soil available nitrogen at 0–100 cm, available phosphorus at 0–40 cm and available potassium at 0–40 cm depths at the end of the growing season. ¶ Significance test of difference by independent-samples t-test using the SPSS 20.0 software package (P<0.05).

Soil nutrients in the four different zones under border and drip irrigation at the end of the growing season are summarized in Table 4. Throughout the field the remaining soil available nitrogen at 0–100 cm depth increased from 16.3 mg·kg^-1 N under border irrigation to 30.4 mg·kg⁻¹ N under drip irrigation, and available phosphorus at 0–40 cm depth increased from 12.1 mg·kg⁻¹ P under border irrigation to 14.8 mg·kg^-1 P under drip irrigation. Both the soil available nitrogen and available phosphorus under drip irrigation were significantly higher than in the same soil texture zones under border irrigation (P<0.05). Accumulation of soil available nitrogen and available phosphorus is also shown in box plots (Fig. 4) and spatial interception plots (Fig. 5). However, throughout the field the difference in remaining soil available potassium at 0–40 cm depth was not significant under the border and drip irrigation, with an average of 104.3 mg·kg^-1 K under border irrigation and 101.9 mg·kg⁻¹ K under drip irrigation. Only in Zone 2 was the difference in remaining available potassium significant (P<0.05) with average values of 103.5 mg·kg⁻¹ K under border irrigation and 95.4 mg·kg^-1 K under drip irrigation (Table 4).

3.2 Grain yield and water consumption under both irrigation methods

Maize grain yield and water consumption are summarized in Table 5. The maximum yields were 6.25 and 7.02 Mg·ha^-1 in 2015 and 2016, respectively. RY ranged from 0.39 to 1.0 with an average of 0.70 under border irrigation in 2015. RY ranged from 0.48 to 0.97 with an average of 0.71 under drip irrigation in 2016 (Table 5). The difference in RY between the two irrigation methods was not significant in the same soil zone (P<0.05). In the same growing seasons Zone 4 was the zone with the lowest RY under both irrigation methods. However, the only significant difference was between Zones 2 and 4 under border irrigation (P<0.05). Drip irrigation showed a marked advantage in increasing RY at area A (mainly belonging to Zones 3 and 4 with coarser-textured soil) and area B (belonging to Zone 1 with finer-textured subsoil) where soil nitrogen and salinity were both at higher levels. However, drip irrigation decreased RY at area C (Fig. 6(a)). Drip irrigation decreased the yield variability with a coefficient of variation of 12.1% throughout the field compared with 18.8% under border irrigation.

Tab.5 Crop yield and water consumption

Zone (year)	RY (mean±SD)	CV (%)	Yieldmax (Mg·ha−1)	P mm	I mm	ET† mm	ET/(P+I) (%)
All Zones (2015)	0.70±0.13	18.8	6.25	149	679	508	61
Zone 1 (2015)	0.70±0.14 ab¶	19.7
Zone 2 (2015)	0.73±0.11 a	15.6
Zone 3 (2015)	0.71±0.14 ab	19.1
Zone 4 (2015)	0.65±0.14 b	21.2
All Zones (2016)	0.71±0.09	12.1	7.02	115	452	518	91
Zone 1 (2016)	0.71±0.08 a	11.3
Zone 2 (2016)	0.71±0.10 ab	13.5
Zone 3 (2016)	0.73±0.08 a	10.5
Zone 4 (2016)	0.69±0.09 ab	12.7

Note: RY, the relative yield; CV, the coefficient of variation; Yieldmax, the observed maximum yield (dry grain); P and I, the total precipitation and irrigation water during the growing season, respectively; ET, crop evapotranspiration by the eddy covariance method. ¶ The significance test of difference by independent-samples t-test using the SPSS 20.0 software package (P<0.05). † ET was the eddy covariance result in both years, according to other literature on the same field. ET from April 15 to September 16, 2015 was the eddy covariance results from Qin et al.[27]; ET from April 25 to September 16, 2016 was adjusted from the result of He et al.[28] by plus or minus the ET of missing or extra days.

Irrigation water input during the growing season was 679 and 452 mm under border and drip irrigation, respectively. Drip irrigation saved 227 mm irrigation water during a drier growing season compared with border irrigation and this was 44% of the average ET throughout the growing season (ET, the eddy covariance measurements of the same field by Qin et al.^[27] and He et al.^[28]). Also, 91% of the total water input (precipitation plus irrigation) was converted into crop ET under drip irrigation but the corresponding percentage was only 61% under border irrigation (Table 5).

3.3 Yield-limiting factor analysis and yield gap quantification

The parameters of the BLA method, the optimum range and the acceptable range, are listed in Table 3. We recommend the optimum range as the thresholds for field management and yield-limiting factors outside the acceptable range should be given priority in closing yield gaps. Also, the area proportion of different stress intensity by different yield-limiting factors is also summarized in Table 3. The area in which yield was constrained by each yield-limiting factor (deficiency or excessive) is delineated in the interpolation maps (Fig. 6).

Yield gaps determined by yield-limiting factors are summarized for different soil texture zones in Table 6. Throughout the field the yield gap determined by all the currently considered yield-limiting factors (D) was 7.9% under border irrigation and 9.7% under drip irrigation (Table 6). However, yield gaps differed with soil zone and irrigation method. Zone 1 with finer-textured soil had the minimum yield gaps determined by currently considered factors and Zone 4 with coarser-textured soil had the maximum yield gaps under both irrigation methods. The yield gaps determined by currently considered factors (D) ranged from 5.9% at Zone 1 to 12.4% at Zone 4 under border irrigation, and 7.0% at Zone 1 to 13.8% at Zone 4 under drip irrigation. The RY was expected to increase from ~70% at present to ~80% (range 76% to 84% in different soil zones under both irrigation methods) if the currently considered yield-limiting factors were managed to the optimum ranges. The yield gaps determined by unknown factors (D_ε) throughout the field were 21.6% under border irrigation and 19.2% under drip irrigation. Zone 1 had the maximum unexplained yield gap with values of 24.2% under border irrigation and 21.8% under drip irrigation.

Tab.6 Yield gaps determined by different yield-limiting factors in different soil zones

Zone (year)	Count		RY (mean±SD)			RY' (mean±SD)		Δ (%)		Δε (%)
All Zones (2015)	185		0.70±0.13			0.78±0.13		7.9		21.6
Zone 1 (2015)	71		0.70±0.14			0.76±0.14		5.9		24.2
Zone 2 (2015)	53		0.73±0.11			0.82±0.12		8.5		18.0
Zone 3 (2015)	36		0.71±0.14			0.79±0.13		7.9		21.3
Zone 4 (2015)	25		0.65±0.14			0.77±0.11		12.4		22.5
All Zones (2016)	184		0.71±0.09			0.81±0.09		9.7		19.2
Zone 1 (2016)	82		0.71±0.08			0.78±0.08		7.0		21.8
Zone 2 (2016)	64		0.71±0.10			0.82±0.09		11.8		17.6
Zone 3 (2016)	23		0.73±0.08			0.84±0.09		11.1		15.9
Zone 4 (2016)	15		0.69±0.09			0.83±0.10		13.8		16.9
Zone (Year)		ΔNum (%)		ΔPropAWC (%)	ΔAN (%)		ΔEC (%)		ΔAK (%)		ΔpH (%)
All Zones (2015)		4.1		1.4	1.1		0.4		0.4		0.6
Zone 1 (2015)		3.9		0.1	0.8		0.5		0.1		0.7
Zone 2 (2015)		5.3		1.0	1.0		0.6		0.0		0.5
Zone 3 (2015)		3.5		1.2	1.2		0.1		1.5		0.5
Zone 4 (2015)		3.3		6.0	1.9		0.1		0.4		0.8
All Zones (2016)		4.9		2.5	1.0		0.6		0.7		0.0
Zone 1 (2016)		3.4		0.6	1.6		0.9		0.4		0.0
Zone 2 (2016)		7.0		3.1	0.6		0.4		0.7		0.0
Zone 3 (2016)		4.4		4.5	0.2		0.0		1.8		0.1
Zone 4 (2016)		4.8		7.5	0.1		0.1		1.4		0.0

Note: RY, the observed relative yield; RY', the expected relative yield when all the factors considered (unknown factors ε not included) were managed to the optimum ranges; RY' = RY+ D; D, the total yield gap of the current factors considered (unknown factors ε not included); D = DNum + DPropAWC + DAN + DEC + DAK + DpH; DNum, DPropAWC, DAN, DEC, DAK, DpH, Dε, yield gap determined by female plant density, soil available water content, soil available nitrogen, soil electrical conductivity, soil available potassium, soil pH, and unknown factors, respectively; RY' + Dε = 1.

Irrespective of irrigation method the first three yield-limiting factors were plant density (D_Num), soil available water content (D_PropAWC), and soil available nitrogen (D_AN) throughout the field (Table 6). The yield gaps determined by plant density (D_Num) were 4.1% and 4.9% in 2015 and 2016, respectively. The slightly larger yield gap determined by plant density in 2016 may be explained by the larger areas constrained by female plant deficiency that occurred almost throughout the study area compared with 81% of the study area in 2015 (Table 3). Throughout the field the yield gaps determined by soil available water content (D_PropAWC) were 1.4% under border irrigation and 2.5% under drip irrigation (Table 6). The larger yield gap determined by soil available water content under drip irrigation was the result of the larger area affected by water stress. Yield gap determined by soil available water content was considerable in all soil textures except Zone 1 under drip irrigation with values of 3.1%, 4.5% and 7.5% in Zones 2, 3 and 4, respectively. However, only Zone 4 had a large yield gap of 6.0% determined by soil available water content under border irrigation (Table 6). In Zone 4, soil available water content had become the predominant yield-limiting factor under both irrigation methods (Table 6). The yield gap determined by soil available nitrogen (D_AN) was similar throughout the field with values of 1.1% and 1.0% under the border and drip irrigation, respectively. However, the cause of the yield gap was different. Both soil nitrogen deficit and excess occurred under border irrigation, with soil nitrogen deficit occurring mainly in the coarser-textured Zones 3 and 4, and soil nitrogen surplus occurred mainly in finer-textured Zone 1 (Fig. 6(c)). However, only soil nitrogen surplus occurred under drip irrigation, mainly in Zone 1 (Fig. 6(h)). Accumulation of soil available nitrogen led to the area proportion affected by excessive nitrogen increasing from 14% of the total area under border irrigation to 44% under drip irrigation (Table 3).

Although other yield-limiting factors considered had little influence on yield gaps throughout the field, they need to be considered in specific soil texture zones. For example, a large area at Zone 1 was affected by soil salinity at the early growth stages, accompanied by soil nitrogen surplus (Fig. 6(c,d,h,i), and Zones 3 and 4 with coarser-textured topsoil were also affected by available potassium deficit (Table 6, Fig. 6(b,g)).

4 DISCUSSION

4.1 Thresholds of yield-limiting factors for field management

Studies have been conducted to recommend the thresholds for field management with different methods such as BLA^[18,19], the Cate-Nelson graphical procedure^[29], or least square regression between crop yield and yield-limiting factors. The results may differ slightly from each other. As these ranges are the results of analyzing the relationship between the yield of a specific crop and yield-limiting factors at specific sampling time, they would be misleading if the sampling time or crop species were ignored.

Here, the optimum range of female parents was (7.6–7.8) × 10⁴ plants ha⁻¹, a lower value than the results of Jiang et al.^[22] who found that maximum yield was achieved at a female plant density of (9.4–10.6) × 10⁴ plants ha⁻¹ (a further 20% male parents were not included). The optimum range was reached earlier in our study because the BLA method is concerned with the upper boundary points of scatter, and factors other than plant density were at optimum levels. Also, the narrow optimum range was due to few data being distributed at higher plant density in our study. However, this range is consistent with the optimum plant density of (9.7–16.4) × 10⁴ plants ha⁻¹ for hybrids released between 1965 and 2010^[30] if the male parents are taken into account. A slightly higher plant density is recommended for increased yield stability as modern hybrids are tolerant to high plant density^[30].

According to the BLA of soil available water content in 2016 we recommend 37% to 51% of soil available water content as the lower limit for irrigation management. This is consistent with the recommendation of 35% to 55% (the maximum allowable water depletion is 45% to 65% of total soil available water) to meet full water requirements^[31,32]. Soil available water content exceeded 100% at some sampling sites (Fig. 2), especially under border irrigation, so that soil water content exceeded field capacity. This is related to the common occurrence of finer-textured subsoil in the field, and this may lead to a longer time required when downward flux becomes negligible after large rates of irrigation^[33]. The linear-plateau model was used in fitting the relationship of yield response to soil available water content, as there was no apparent yield reduction by excessive irrigation in our study. However, large rates of irrigation under border irrigation should be avoided because of the risk of hypoxia^[34].

Soil nutrients at the end of the growing season were used to indicate thresholds for fertilizer management, as soil nutrients during this period tended to be lower than that of the earlier growth stages. Another consideration was that residual soil nitrogen was more stable than at other growth stages, as nitrogen was applied more than once and complex movement with soil water. The optimum range of average soil available nitrogen in the 100 cm was 6.04 to 26.1 mg·kg⁻¹ (Table 3). However, high nitrogen inputs should be avoided because of the risk of nitrate leaching^[35] and the negative influence of salinity on germination and early seedling establishment^[36]. The lower limit of the optimum range in our study is comparable to the critical level of 8 mg NO₃⁻-N kg⁻¹ (bulk density assumed to be 1.49 Mg·m⁻³)^[29]. Soil phosphorus limitation did not occur according to the BLA method, with soil available phosphorus ranging from 5 (at the end of the growing season) to 51 mg·kg⁻¹ P (at the start of the growing season when maize was sown). However, soil available phosphorus outside this range should be considered for phosphorus deficiency or excess, which may lead to deficits of micronutrients such as Zn or Fe^[37]. We recommend the lower limit of 81 mg·kg⁻¹ K for soil available potassium which is comparable to the lower limit of 90 mg·kg⁻¹ K according to other studies^[38].

4.2 The current situation and field management problems

Here, irrespective of whether border or drip irrigation was used, current field management achieved a RY of 70% for the whole field level, which was at a high level given the yield plateau 75% to 85% of the potential yield in the literature^[8]. However, it was achieved at the cost of high resource consumption.

Under border irrigation only 61% of the total water input (precipitation plus irrigation) was converted to crop ET (Tab.5), so that the remaining water was wasted in the form of deep percolation or water storage in the root zone at the end of the growing season. Drip irrigation showed a substantial advantage over border irrigation in saving irrigation water. Other studies also report that drip irrigation required 30% to 56% less water with an equal or higher yield compared with the sprinkler or border irrigation^[39,40]. The poor performance of border irrigation is well known, as irrigation efficiency is significantly affected by field microtopography even with laser leveling^[41]. Also, application efficiency (Ea, the percentage of irrigation water contributing to target irrigation depth) was only 49% to 82% under border irrigation^[42,43], but reached ~ 90% under both sprinkle and drip irrigation^[42]. Expanding drip irrigation may ease the water scarcity situation in north-west China.

Using drip irrigation calls for corresponding changes in fertilizer management because current high fertilizer inputs greatly threaten resources and the environment. Fertilizer overuse in China has been widely reported in recent years^[1,2,10,17]. Here, despite the application of ˃ 400 kg·ha^-1 N, nitrogen deficit occurred at some locations in the coarser-textured Zones 3 and 4 under border irrigation (Fig. 6(c)). However, drip irrigation succeeded in maintaining soil nitrogen within the root zone, as top-dressing nitrogen was fertigated with greater frequency but lower quantity at each event and there was less nitrate leaching due to reduced irrigation water volume. As a result, nitrogen deficit was eliminated. However, only one growing season after using drip irrigation the area in which remaining available nitrogen at 0–100 cm depth exceeded the upper limit of the optimum range (26.1 mg·kg⁻¹ or 369 kg·ha⁻¹ N) increased to 44% under drip irrigation compared to 14% under border irrigation (Table 3). This upper limit reached 1.8 times the maize demand of 200 kg·ha⁻¹ N throughout the growing season^[31], and crop yield began to be constrained by nitrogen surplus (Fig. 2(c)). In these areas, reducing nitrogen application would save nitrogen, reduce the energy used by the fertilizer industry^[1], decrease agricultural impact on drinking water, atmosphere or ecosystem^[44], increase the partial factor productivity of nitrogen^[14,45], and also increase crop yields.

4.3 Measures and suggestions for closing yield gaps and sustainable agriculture

We recommend the optimum range as the threshold for field management as discussed in Section 4.1. Any practices that can maintain the yield-limiting factors within the optimum ranges should be taken into consideration. Here, plant density was the predominant factor limiting yield throughout the field irrespective of border or drip irrigation. This is consistent with the conclusion of Chen et al. that plant density makes a slightly higher contribution to the yield gap than other management practices^[14]. Both the emergence and sowing rate determined plant density. Emergence was influenced by numerous variables such as soil temperature, soil moisture, soil air quality, soil burial, compaction, soil fertility and salinity^[36,46], which can be manipulated indirectly or directly through field management^[46]. Any practices that can guarantee emergence should be considered, including delaying sowing date to avoid low temperatures at the germination-emergence stage, ensuring suitable soil moisture content, lowering soil salinity and even increasing the sowing rate if reduced emergence is unavoidable.

Soil available water content and available nitrogen were the other two important yield-limiting factors following plant density, and they can be manipulated through irrigation and nitrogen management. Modern irrigation methods such as drip irrigation should be given special priority, because of the advantages in irrigation management even in soils with small ranges of non-limiting water^[47]. Other advantages include lower consumption of irrigation water^[40], improved soil physical properties^[48], higher soil aeration^[34], and reduced weed growth^[40] compared to border irrigation. Also, balanced fertilizer management was needed to maintain a high yield level with lower fertilizer input. However, it seems to be difficult or impossible to achieve the optimum range for the whole field with current uniform irrigation and fertilizer management, as these factors are closely related to soil texture and spatial variability commonly occurs in soils. Hence, site-specific management is needed. Frequent irrigation with a smaller volume each time would be favorable in the coarser-textured soil zones as irrigation frequency had a large influence on grain yield even using the same amount of total irrigation^[49]. In the finer-textured soil zones in which nitrogen and salt accumulation occurred, low fertilizer rates are recommended. Utilizing site-specific management can improve nitrogen management and increase nitrogen use efficiency, and reduce nitrogen losses to the environment^[50], which was also true for irrigation management. This supports the use of site-specific management on larger farms and at regional scale.

According to the yield gap analysis of different yield-limiting factors, multiple goals of 8% to 10% RY increase at the whole field scale, reducing agricultural inputs and promoting agricultural sustainability may be realized by managing the important factors within the optimum ranges. However, this relies on optimizing agronomic management^[45], balanced fertilizer management, modern irrigation methods, site-specific management^[50], and other technical support of precision agriculture, such as the timely acquisition of field information, variable-rate fertilization and irrigation techniques^[51]. However, further yield increases appear not to be economical or even possible when a RY of 78% to 81% was attained throughout the field but a further yield gap of 19% to 22% RY occurred due to unknown factors (Table 6). Unknown factors may include weeds and pests, variation in seed vigor and plant health status, spatiotemporal variability of soil water and nutrients, and random error. Sustainable agriculture must be practiced to ensure future food, water resources and environmental security.

5 CONCLUSIONS

Using a field investigation and yield-limiting factor analysis, the main conclusions of this study are as follows. (1) Whatever the irrigation method, current field management can achieve a RY of 70%. However, drip irrigation consumed less irrigation water and decreased yield variability compared to traditional border irrigation. Drip irrigation is recommended to minimize water scarcity in arid north-west China. (2) The optimum ranges of different yield-limiting factors can be used as the thresholds for field management. About 8% to10% RY can be increased by managing the considered yield-limiting factors within the optimum ranges by relying on the theoretical and technical support of precision agriculture. When a RY of 78% to 81% is reached, further yield increase appears to be uneconomic or impossible. (3) Plant density, followed by soil available water content and remaining available nitrogen, explained the major causes of yield gaps by considered yield-limiting factors (not including the unknown factors). Also, in coarser-textured soil zones soil available water content explained a greater yield gap than the other two factors. Hence, plant density, irrigation and nitrogen management must be given special emphasis in closing yield gaps and developing sustainable agricultural systems.

References

[1] Foley J A, Ramankutty N, Brauman K A, Cassidy E S, Gerber J S, Johnston M, Mueller N D, O’Connell C, Ray D K, West P C, Balzer C, Bennett E M, Carpenter S R, Hill J, Monfreda C, Polasky S, Rockström J, Sheehan J, Siebert S, Tilman D, Zaks D P M. Solutions for a cultivated planet. Nature, 2011, 478(7369): 337–342

[2] Mueller N D, Gerber J S, Johnston M, Ray D K, Ramankutty N, Foley J A. Closing yield gaps through nutrient and water management. Nature, 2012, 490(7419): 254–257

[3] Tilman D, Balzer C, Hill J, Befort B L. Global food demand and the sustainable intensification of agriculture. Proceedings of the National Academy of Sciences of the United States of America, 2011, 108(50): 20260–20264

[4] Alexandratos N, Bruinsma J. World agriculture towards 2030/2050: the 2012 revision. Rome: FAO, 2012

[5] Food and Agriculture Organization of the United Nations (FAO). FAOSTAT statistics database. 2017. Avaible at FAO website on May 1, 2019

[6] Cassman K G, Dobermann A, Walters D T, Yang H S. Meeting cereal demand while protecting natural resources and improving environmental quality. Annual Review of Environment and Resources, 2003, 28(1): 315–358

[7] Long S P, Marshall-Colon A, Zhu X G. Meeting the global food demand of the future by engineering crop photosynthesis and yield potential. Cell, 2015, 161(1): 56–66

[8] van Ittersum M K, Cassman K G, Grassini P, Wolf J, Tittonell P, Hochman Z. Yield gap analysis with local to global relevance—a review. Field Crops Research, 2013, 143: 4–17

[9] Rosegrant M W, Ringler C, Zhu T J. Water for agriculture: maintaining food security under growing scarcity. Annual Review of Environment and Resources, 2009, 34(1): 205–222

[10] Lu C Q, Tian H Q. Global nitrogen and phosphorus fertilizer use for agriculture production in the past half century: shifted hot spots and nutrient imbalance. Earth System Science Data, 2017, 9(1): 181–192

[11] Yemefack M, Jetten V G, Rossiter D G. Developing a minimum data set for characterizing soil dynamics in shifting cultivation systems. Soil & Tillage Research, 2006, 86(1): 84–98

[12] D’Hose T, Cougnon M, De Vliegher A, Vandecasteele B, Viaene N, Cornelis W, Van Bockstaele E, Reheul D. The positive relationship between soil quality and crop production: a case study on the effect of farm compost application. Applied Soil Ecology, 2014, 75: 189–198

[13] Yang W G, Zheng F L, Han Y, Wang Z L, Yi Y, Feng Z Z. Investigating spatial distribution of soil quality index and its impacts on corn yield in a cultivated catchment of the Chinese mollisol region. Soil Science Society of America Journal, 2016, 80(2): 317–327

[14] Chen G F, Cao H Z, Liang J, Ma W Q, Guo L F, Zhang S H, Jiang R F, Zhang H Y, Goulding K W T, Zhang F S. Factors affecting nitrogen use efficiency and grain yield of summer maize on smallholder farms in the North China Plain. Sustainability, 2018, 10(2): 363

[15] Li X, Zhang X, Niu J, Tong L, Kang S, Du T, Li S, Ding R. Irrigation water productivity is more influenced by agronomic practice factors than by climatic factors in Hexi Corridor, Northwest China. Scientific Reports, 2016, 6(1): 37971

[16] Li X L, Tong L, Niu J, Kang S Z, Du T S, Li S E, Ding R S. Spatio-temporal distribution of irrigation water productivity and its driving factors for cereal crops in Hexi Corridor, Northwest China. Agricultural Water Management, 2017, 179: 55–63

[17] Shen J B, Zhu Q C, Jiao X Q, Ying H, Wang H L, Wen X, Xu W, Li T Y, Cong W F, Liu X J, Hou Y, Cui Z L, Oenema O, Davies W J, Zhang F S. Agriculture Green Development: a model for China and the world. Frontiers of Agricultural Science and Engineering, 2020, 7(1): 5–13

[18] Evanylo G K, Sumner M E. Utilization of the boundary line approach in the development of soil nutrient norms for soybean production. Communications in Soil Science and Plant Analysis, 1987, 18(12): 1379–1401

[19] Schnug E, Heym J, Murphy D P. Boundary line determination technique (BOLIDES). In: Robert P C, Rust R H, Larson W E, eds. Site-specific management for agricultural systems. American Society of Agronomy, Crop Science Society of America, Soil Science Society of America, 1995, 899–908

[20] Hajjarpoor A, Soltani A, Zeinali E, Kashiri H, Aynehband A, Vadez V. Using boundary line analysis to assess the on-farm crop yield gap of wheat. Field Crops Research, 2018, 225: 64–73

[21] Shatar T M, McBratney A B. Boundary-line analysis of field-scale yield response to soil properties. Journal of Agricultural Science, 2004, 142(5): 553–560

[22] Jiang X L, Tong L, Kang S Z, Li F S, Li D H, Qin Y H, Shi R C, Li J B. Planting density affected biomass and grain yield of maize for seed production in an arid region of Northwest China. Journal of Arid Land, 2018, 10(2): 292–303

[23] Walkley A. A critical examination of a rapid method for determining organic carbon in soils—effect of variations in digestion conditions and of inorganic soil constituents. Soil Science, 1947, 63(4): 251–264

[24] van den Berg M, Klamt E, van Reeuwijk L P, Sombroek W G. Pedotransfer functions for the estimation of moisture retention characteristics of Ferralsols and related soils. Geoderma, 1997, 78(3–4): 161–180

[25] R Core Team. R: A Language and Environment for Statistical Computing (Version 3.5.2). Vienna, Austria: R Foundation for Statistical Computing,2018

[26] Bartier P M, Keller C P. Multivariate interpolation to incorporate thematic surface data using inverse distance weighting (IDW). Computers & Geosciences, 1996, 22(7): 795–799

[27] Qin S J, Li S E, Kang S Z, Du T S, Tong L, Ding R S. Can the drip irrigation under film mulch reduce crop evapotranspiration and save water under the sufficient irrigation condition? Agricultural Water Management, 2016, 177: 128–137

[28] He Q S, Li S E, Kang S Z, Yang H B, Qin S J. Simulation of water balance in a maize field under film-mulching drip irrigation. Agricultural Water Management, 2018, 210: 252–260

[29] Magdoff F R, Ross D, Amadon J. A soil test for nitrogen availability to corn. Soil Science Society of America Journal, 1984, 48(6): 1301–1304

[30] Di Matteo J A, Ferreyra J M, Cerrudo A A, Echarte L, Andrade F H. Yield potential and yield stability of Argentine maize hybrids over 45 years of breeding. Field Crops Research, 2016, 197: 107–116

[31] Doorenbos J, Kassam A H. FAO Irrigation and Drainage Paper 33: Yield response to water. Rome: FAO, 1979

[32] Panda R K, Behera S K, Kashyap P S. Effective management of irrigation water for maize under stressed conditions. Agricultural Water Management, 2004, 66(3): 181–203

[33] Cassel D K, Nielsen D R. Field capacity and available water capacity. In: Klute A, eds. Methods of soil analysis: Part 1—physical and mineralogical methods. Madison, WI: Soil Science Society of America, American Society of Agronomy, 1986, 901–926

[34] Bhattarai S P, Su N H, Midmore D J. Oxygation unlocks yield potentials of crops in oxygen-limited soil environments. Advances in Agronomy, 2005, 88: 313–377

[35] Perego A, Basile A, Bonfante A, De Mascellis R, Terribile F, Brenna S, Acutis M. Nitrate leaching under maize cropping systems in Po Valley (Italy). Agriculture, Ecosystems & Environment, 2012, 147: 57–65

[36] Farooq M, Hussain M, Wakeel A, Siddique K H M. Salt stress in maize: effects, resistance mechanisms, and management. A review. Agronomy for Sustainable Development, 2015, 35(2): 461–481

[37] Jones J B. Plant nutrition and soil fertility manual. 2nd eds. New York: CRC Press, 2012

[38] Mallarino A P, Oyarzabal E S, Hinz P N. Interpreting within-field relationships between crop yields and soil and plant variables using factor analysis. Precision Agriculture, 1999, 1(1): 15–25

[39] Adamsen F J. Irrigation method and water quality effect on peanut yield and grade. Agronomy Journal, 1989, 81(4): 589–593

[40] Camp C R, 0. Subsurface drip irrigation: a review. Transactions of the ASAE. American Society of Agricultural Engineers, 1998, 41(5): 1353–1367

[41] Playán E, Faci J M, Serreta A. Modeling microtopography in basin irrigation. Journal of Irrigation and Drainage Engineering, 1996, 122(6): 339–347

[42] Battikhi A M, Abu-Hammad A H. Comparison between the efficiencies of surface and pressurized irrigation systems in Jordan. Irrigation and Drainage Systems, 1994, 8(2): 109–121

[43] Lecina S, Playán E, Isidoro D, Dechmi F, Causapé J, Faci J M. Irrigation evaluation and simulation at the Irrigation District V of Bardenas (Spain). Agricultural Water Management, 2005, 73(3): 223–245

[44] Socolow R H. Nitrogen management and the future of food: lessons from the management of energy and carbon. Proceedings of the National Academy of Sciences of the United States of America, 1999, 96(11): 6001–6008

[45] Song Z W, Feng X M, Lal R, Fan M M, Ren J, Qi H, Qian C R, Guo J R, Cai H G, Cao T H, Yu Y, Hao Y B, Huang X M, Deng A X, Zheng C Y, Zhang J, Zhang W J. Optimized agronomic management as a double-win option for higher maize productivity and less global warming intensity: a case study of Northeastern China. Advances in Agronomy, 2019, 157: 251–292

[46] Forcella F, Benech Arnold R L, Sanchez R, Ghersa C M. Modeling seedling emergence. Field Crops Research, 2000, 67(2): 123–139

[47] Letey J. Relationship between soil physical properties and crop production. Advances in Soil Science, 1958, 1: 277–294

[48] Nelson S D, Terry R E. The effects of soil physical properties and irrigation method on denitrification. Soil Science, 1996, 161(4): 242–249

[49] El-Hendawy S E, Schmidhalter U. Optimum coupling combinations between irrigation frequency and rate for drip-irrigated maize grown on sandy soil. Agricultural Water Management, 2010, 97(3): 439–448

[50] Khosla R, Fleming K, Delgado J A, Shaver T M, Westfall D G. Use of site-specific management zones to improve nitrogen management for precision agriculture. Journal of Soil and Water Conservation, 2002, 57(6): 513–518

[51] Balafoutis A, Beck B, Fountas S, Vangeyte J, van der Wal T, Soto I, Gómez-Barbero M, Barnes A, Eory V. Precision agriculture technologies positively contributing to GHG emissions mitigation, farm productivity and economics. Sustainability, 2017, 9(8): 1339

Acknowledgements

This research was funded by the National Natural Science Foundation of China (51725904, 51621061, 51861125103) and the Discipline Innovative Engineering Plan (111 Program, B14002). We thank Mr. Quan Lu (Shiyanghe Experimental Station, China Agricultural University, China) for field management. We are grateful to the anonymous reviewers and the editors for their work on this manuscript.

Compliance with ethics guidelines

Xiuwei Guo, Manoj Kumar Shukla, Di Wu, Shichao Chen, Donghao Li, and Taisheng Du declare that they have no conflicts of interest or financial conflicts to disclose.

This article does not contain any studies with human or animal subjects performed by any of the authors.

RIGHTS & PERMISSIONS

The Author(s) 2020. Published by Higher Education Press. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0)