Nitrogen and phosphorus leaching losses under cropping and zone-specific variable-rate irrigation

Overview of site and soils

The study was conducted at Massey University’s No. 1 Farm, near Palmerston North, New Zealand. The farm is located at 40.22°S and 175.36°E. The climate of the region is temperate. Median annual rainfall at Palmerston North (1981–2010) was 900 mm; mean annual temperature range was 9.1°C (Chappell 2015).

The study area (1.2 ha) is on a mixed-cropping site, typically cropped with peas, beans, and spring wheat, with sheep grazing between cropping events. The site has a centre-pivot irrigator with variable rate control that irrigates two management zones. The two management zones are based on soil type. These two zones were delineated using information from a soil survey (Pollok et al. 2003), an electromagnetic survey (El-Naggar et al. 2021), and measurements of a range of soil physical properties. Zone 1 soil is a Manawatū fine sandy loam, a deep, free-draining soil classified as a Fluvial Recent Soil (Hewitt 2010). Zone 2 is a Manawatū silt loam classified as a Fluvial Recent Soil and is poorly drained. Both soils are classified as a Fluvisol in the FAO Reference Soil Group (IUSS Working Group WRB 2015).

Soil analyses were undertaken to characterise the soils. Soil available water capacity (AWC) was determined in 0.1-m increments to 1 m depth (El-Naggar et al. 2020). Soil samples for physical analysis were collected in August 2016 using undisturbed cores sampled at 20-cm intervals to 1 m depth. Sampling details for the site are described in El-Naggar et al. (2021), with general physical methods in Drewry et al. (2021c). Analyses undertaken at the Manaaki Whenua – Landcare Research Soil Physics Laboratory were: particle size distribution, bulk density, porosity, and air permeability (Gradwell and Birrell 1979). Soil samples for chemical analysis were collected in August 2016. Analyses were: Olsen P, total carbon (C) and N, pH, cation exchange capacity, and exchangeable K, Mg, and Ca, with analyses undertaken at the Massey University School of Agriculture and Environment laboratory.

Crop and irrigation management

The site was cropped with a series of monoculture crops for harvest, including peas (Pisum sativum), beans (Phaseolus vulgaris), spring wheat (Triticum aestivum), and oats (Avena sativa), followed by Barkant turnips (Brassica rapa), plantain (Plantago lanceolata), chicory (Cichorium intybus), ryegrass (Lolium perenne), red clover (Trifolium pratense), and white clover (Trifolium repens) in pasture mixes grazed by animals, over the period 2016–2021 (Table 1). The cropped site is typically grazed with sheep when not in crop or is fallow, termed ‘mixed cropping’.

The 86-m irrigator had one span containing 31 sprinklers, each with a spray radius of 5 m and a flow rate of 26.3 m³/h. The pivot provided a mechanism to apply known and fixed amounts of irrigation. Soil water balance (SWB) scheduling used the FAO56-ET method to calculate daily soil water deficits to determine crop water requirements (El-Naggar et al. 2020).

Experimental design and drainage flux meters

An electromagnetic (EM) survey of the site was conducted to identify the two management zones. Zone maps were derived using geostatistical methods, as reported by El-Naggar et al. (2021). Briefly, electromagnetic and gamma radiometric data were interpolated into regular 5-m resolution grids using ordinary kriging. These were then used to derive zone maps using unsupervised clustering; further details are described elsewhere (Hedley et al. 2016; El-Naggar et al. 2021). The VRI pattern was developed for this site based on this information.

Twelve DFMs were constructed for each of two management zones using an existing design (Gee et al. 2002, 2009; Norris et al. 2023). The wick length was 0.06 m and the control tube (convergence ring) was approximately 0.015 m in height. Each DFM was installed 3–4 m apart with its top at least 0.65 m below the soil surface. The trial site was designed with three plots applied to each of the two soil zones, with four replicate DFMs grouped in each plot (Fig. 1). The DFMs were located on the 76 m arc from the centre of the irrigator.

Fig. 1.

Experimental layout of the site with soil types and management zones.

The 24 DFMs were installed (28 April to 26 May 2016) to collect drainage water below 0.65 m soil depth. For installation, a hole was prepared 0.20 m in diameter, and 1.8 m deep, with a hand auger. Stones were encountered at approximately 1.2–1.5 m depth in Zone 2. At installation, the 0.20-m section of pipe above the wick was repacked with soil taken from the same depths. The wick in the funnel section was wet up and then diatomaceous earth added and packed tightly. Finally, 99% pure silica sand was added, and water carefully poured over this until it measured approximately 0.01 m in depth above the sand, removing air bubbles. Pipes with collection tubes were installed beneath approximately 0.65 m to facilitate drainage water collection, as the DFMs remain in the soil. Further information is available in Karakkattu et al. (2020).

Drainage sampling measurements and water nutrient concentrations

Drainage volumes were measured and collected regularly using a customised pumping system. Drainage was measured on 74 occasions between 26 May 2016 and 23 June 2021. Water quality samples were collected from drainage on 40 sampling dates. Water quality samples were collected every 1–3 weeks during the first year, then approximately 2-monthly or every 80–100 mm of rainfall until approximately late 2019. From then, water quality samples were collected seasonally, due to resourcing constraints. Water samples were analysed for ammonia-N, NO_x-N (nitrate-N + nitrite-N) and reactive phosphorus by the Manaaki Whenua – Landcare Research Environmental Chemistry Laboratory, with methods described in Manaaki Whenua – Landcare Research (2021). All samples were refrigerated at 4°C during storage, and filtered at the laboratory with Whatman No. 42 filter papers prior to analysis. The ammonia-N (NH₄-N) was determined colorimetrically using the indophenol reaction with sodium salicylate and hypochlorite; NO₃-N by Cd reduction and NEDD colorimetry; and NO_x-N is reported (nitrate-N + nitrite-N). For reactive P (i.e. filterable reactive P), the orthophosphate reacted with ammonium molybdate and antimony potassium tartrate under acidic conditions, then a molybdenum blue complex was formed after ascorbic acid reduction. All analyses were performed with a QuikChem 8500 flow injection analyser (Lachat Instruments, Milwaukee, WI, USA).

Soil water balance modelling

For the SWB, a daily time-step model was developed. This model used the FAO56-ET Penman-Monteith method (Allen et al. 1998) and weather data derived from the Palmerston North CliFlo climate station (http://cliflo-niwa.niwa.co.nz/) (NIWA and AgResearch climate station 21963), which was located only 50 m from the site, to calculate the daily ETo (reference evapotranspiration), ETc (crop evapotranspiration) and soil water deficit (SWD) values. The daily SWDs were used to determine crop water requirements for irrigation scheduling and crop management (El-Naggar et al. 2020, 2021). Separate SWBs were constructed per zone using farm and nutrient management (Table 1, Supplementary Fig. S1 in Supplementary material), crop, irrigation, and the AWC data collected per zone.

Statistical analysis

The measured drainage from Zone 1 was found to be much greater than the SWB-modelled drainage. There were fewer drainage and nutrient samples from Zone 2 after 2017, than for Zone 1, so zone-specific scaling was required. The approach was to estimate the parameters of a linear, quadratic, or cubic zone-specific model to estimate the SWB estimates from the DFM measurements. Overall, the cubic model was best and was adopted. The models are zone-specific and the cubic version of the correction is:

where y_corr is the corrected DFM cumulative drainage to match the SWB estimate, for zone i. The uncertainty is ε, assumed to be a Gaussian with a zero mean. Using the zone-specific cubic model to correct the field-measured DFM cumulative flow, the cumulative drainage for all the sample dates was evaluated.

Concentration values below the limit of detection were replaced by half the detection limit. Evaluation indicated the temporal variation is likely to be non-linear over time. A generalised additive mixed model (GAM) with a random effect for the replicate, and smoothing for the date to account for between- and within-year trends, was used. The GAM model was fitted with each nutrient treated as a Gamma with log-link. The long-term temporal trend was accounted for using a (tensor spline) smooth function of the number of years since the start of 2016, and the intra-annual temporal trend accounted for using a (tensor spline) smooth function of the day of the year. Both splines are zone-specific. A random effect was included for the group within each zone, to account for possible random variation at a site. The model used for estimation of the load during each hydrological year included terms for cumulative rainfall and irrigation. The general GAM model form was:

where μ_i,j = E(Y_i,j) and Y_i,j = EF(μ_i,j, ϕ). Y_i,j is a response variable (one of NO_x-N, ammonia-N, or P), EF(μ_i,j, ϕ) is a Gamma distribution with mean μ_i,j and scale parameter ϕ, for zone i. The f_j(Years since 2016) and f_k(Day of year) are zone-specific smooth functions of the number of years since 1 January 2016 and the day or the year, respectively (Wood 2017), which is cyclic.

The terms for cumulative rainfall and irrigation include an interaction with the date. This form is required since rainfall and irrigation have a large proportion of zeros (i.e. no rainfall and/or no irrigation). It is not clear whether adding these terms would improve the model for load, since omitting the terms would be partially offset by the smooth terms f_1,1 and f_2,1.

The fitted GAM used the mgcv::gam() procedure (Wood 2017) using maximum likelihood, so models could be compared. Models using an additive random effect (group, zone, or DFM) were evaluated using the R-squared, the proportion of deviance explained, and the Akaike Information Criterion (AIC) (Table 2). The model with the random effect associated with DFMs was improved over zone and group (also repeated for the other response variables), so was chosen as the final model, explaining 53.5%, 73%, and 69.7% of the total deviance for NO_x-N, ammonia-N, and reactive P, respectively.

The nutrient loads were calculated from cumulative drainage and the nutrient concentrations from the models developed. For annual statistics of drainage and load, the period 1 July to 30 June was used, as it is common practice in the Southern Hemisphere to use this hydrological year. All statistical analyses were carried out using the statistical package ‘R’ ver. 4.2.0 (R Core Team 2022).

Site soil characteristics

The median sand content of Zone 1 ranged from 55 to 95% and increased with depth, while the median sand content of Zone 2 was 38–60% and increased with depth (Table 3). Median silt (Zone 1: 3–31%; zone 2: 29–46%) and clay contents in both zones decreased with depth (Table 3). Zone 1 median cation exchange capacity values were less than in Zone 2 (Table 3). Zone 1 Olsen P values were greater than in Zone 2; bulk density, and total C and N varied between the zones and depths (Table 3). The air permeability, porosity, pH, exchangeable K, Mg, Ca, and sulfate values per zone and depth are presented in the supplementary material (Table S2). Air permeability was greater in Zone 1.

Drainage

The relationship between SWB modelled and field-measured cumulative drainage over time for the zones using the matching procedure in the methods is shown in Fig. 2. For the modelled drainage on an annual basis, the drainage in Zone 1 was less variable between years than in Zone 2. The cumulative rainfall and irrigation are presented in Fig. 3.

Fig. 2.

Relationship between SWB cumulative drainage and predicted cumulative drainage using the cubic matching procedure, over time, for the management zones.

Fig. 3.

Cumulative rainfall, irrigation, measured and modelled drainage using the matching procedure (mm) for the two management zones. The drainage values are calibrated using the cubic polynomial model described in the methods, so values can be negative.

Nutrient concentrations and loads

There is some evidence in Zone 1 drainage concentrations of an intra-annual trend, and this trend is consistent between the DFMs (Figs 4 and S2). In Zone 2, the DFMs had fewer samples beyond 2019, so given the shorter time (2016–2019), the intra-annual trend is not clear for Zone 2 (Fig. S2 in Supplementary material).

Fig. 4.

Drainage nutrient concentrations over time (points), with each curve corresponding to the GAM model estimate for a specific DFM, using a random effect for each DFM. The shaded region around each line indicates plus-or-minus one standard error for each prediction. Arrows correspond to dates when fertiliser were added in each zone. Note the vertical scale is logarithmic.

The nutrient concentrations over time using the GAM model are presented in Fig. 4. Ammonia-N data indicate an annual effect within Zone 1, but this is not clear for Zone 2. The difference in the predictions between DFMs in Zone 1 are quite small, but likely to be greater in Zone 2. The data for NO_x-N suggests an annual effect within Zone 1, but not for Zone 2. For reactive P, the results suggest an annual effect within Zone 1, with large differences between DFMs (Fig. 4). For Zone 2, the reactive P differences between DFMs is not as large and the annual effect is less pronounced, in part likely due to sparse sampling from 2019. Overall, the adopted model gives a reasonable account of nutrient concentrations over time and between DFMs.

The timing of fertiliser application is shown by arrows shown in Fig. 4. Although variations in nutrient can be observed in Fig. 4, the relationship between the dates of fertiliser application and nutrient concentration variation is not clear. This is due, at least in part, to the seasonality in rainfall, irrigation, and the fertiliser application. A formal analysis of the relationship between these quantities would require considerably more time series data than are available.

The nutrient loads are presented in Fig. 5, showing the variability of individual points for each DFM, while the box plot corresponds to the value pooled over all DFMs (in each year and zone). Although variable, Zone 1 had mostly greater NO_x-N loads than Zone 2 (Table 4). The dominant form was nitrate rather than ammonia-N (Table 4). The load of NO_x-N during 2016 for both zones was greater than subsequent years (Fig. 5, Table 4). Zone 2 had greater variability in ammonia-N and NO_x-N loads per DFM than Zone 1. Zone 1 had much greater variability in reactive P loads per DFM than Zone 2. There were no obvious patterns in annual nutrient loads with annual rainfall and irrigation. However, the ammonia-N annual loads in Zone 2 showed a general decreasing pattern with the decreasing annual rainfall during the period.

Fig. 5.

Nutrient loads (kg N/ha/year or kg P/ha/year) of ammonia-N, NO_x-N and reactive phosphorus (P) for the two management zones, per year (1 July to 30 June). The points are for each DFM. The box plot is generated for the pooled value over all DFMs. For boxplots, the median is shown, and the boxes represent the inter-quartile range, (25th and 75th percentile), and whiskers represent 1.5 times the inter-quartile range. Note the scale is logarithmic.

As described in the methods, the model used for estimation of loads during each hydrological year included terms for cumulative rainfall and irrigation. A likelihood ratio test of models, with and without these terms, shows that both have a highly significant beneficial effect for reactive P (P < 0.01), but the benefits for NO_x-N and ammonia-N are less convincing (P = 0.14 and P = 0.017, respectively). The benefit of including these terms is adding high temporal-frequency information, but when the load was estimated over a full hydrological year, that high-frequency information was largely lost. Thus, estimates of load over a full year differ by little, if the effects of rainfall and irrigation were included or excluded. Table S1 in the Supplementary material shows the absolute difference in the estimated median load for each hydrological year for all responses. There is no obvious pattern, except that perhaps differences are larger in the sandy soil.

In summary, the median load of NO_x-N across all 5 years was 119 and 120 kg N/ha/year for Zones 1 and 2, respectively (Table 5). The mean load of NO_x-N across all years was 133 and 121 kg N/ha/year for Zones 1 and 2, respectively. The mean load of reactive P across all 5 years was 0.17 and 0.14 kg P/ha/year for Zone 1 and Zone 2, respectively (Table 5). Table 5 presents the coefficient of variation (COV; ratio of the standard deviation to the mean). Since the mean was greater than the median, the results in Table 5 show that the distributions of each nutrient are right-skewed. The high COV values for ammonia-N and reactive P in zone 2 (i.e. high standard deviation relative to the mean) in Table 5 reflect the high variation in the nutrient load between hydrological years.