Early season N₂O emissions under variable water management in rice systems: source-partitioning emissions using isotope ratios along a depth profile

Soil moisture strongly affects the balance between nitrification, denitrification and N₂O reduction and therefore the nitrogen (N) efficiency and N losses in agricultural systems. In rice systems, there is a need to improve alternative water management practices, which are designed to save water and reduce methane emissions but may increase N₂O and decrease nitrogen use efficiency. In a field experiment with three water management treatments, we measured N₂O isotope ratios of emitted and pore air N₂O (δ¹⁵N, δ¹⁸O and site preference, SP) over the course of 6 weeks in the early rice growing season. Isotope ratio measurements were coupled with simultaneous measurements of pore water ${\mathrm{NO}}_{\mathrm{3}}^{-}$, ${\mathrm{NH}}_{\mathrm{4}}^{+}$, dissolved organic carbon (DOC), water-filled pore space (WFPS) and soil redox potential (Eh) at three soil depths. We then used the relationship between SP × δ¹⁸O-N₂O and SP × δ¹⁵N-N₂O in simple two end-member mixing models to evaluate the contribution of nitrification, denitrification and fungal denitrification to total N₂O emissions and to estimate N₂O reduction rates. N₂O emissions were higher in a dry-seeded + alternate wetting and drying (DS-AWD) treatment relative to water-seeded + alternate wetting and drying (WS-AWD) and water-seeded + conventional flooding (WS-FLD) treatments. In the DS-AWD treatment the highest emissions were associated with a high contribution from denitrification and a decrease in N₂O reduction, while in the WS treatments, the highest emissions occurred when contributions from denitrification/nitrifier denitrification and nitrification/fungal denitrification were more equal. Modeled denitrification rates appeared to be tightly linked to nitrification and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ availability in all treatments; thus, water management affected the rate of denitrification and N₂O reduction by controlling the substrate availability for each process (${\mathrm{NO}}_{\mathrm{3}}^{-}$ and N₂O), likely through changes in mineralization and nitrification rates. Our model estimates of mean N₂O reduction rates match well those observed in ¹⁵N fertilizer labeling studies in rice systems and show promise for the use of dual isotope ratio mixing models to estimate N₂ losses.

1 Introduction

Atmospheric nitrous oxide (N₂O) concentrations continue to rise, and with a global warming potential 298 times that of CO₂, N₂O is a significant contributor to global warming (IPCC, 2007; Ravishankara et al., 2009). Agriculture is estimated to be responsible for roughly 60 % of anthropogenic N₂O emissions (Smith et al., 2008). Considering this, the quantification of field-scale N₂O emissions has been the focus of many studies in the last decades and much progress has been made in identifying agricultural management practices, soil and climate variables that influence emissions (Mosier et al., 1998; Verhoeven et al., 2017; Venterea et al., 2012). However, it remains difficult to quantitatively determine the microbial source processes of emitted N₂O in the field, and knowledge gaps remain in our understanding of how N₂O production and reduction processes change with both time and depth. More specific knowledge of process dynamics is therefore needed to inform and improve biogeochemical models.

Studying N cycling in rice systems offers a unique opportunity to study processes of N₂O production and reduction. Firstly, there is a strong need to develop alternative water management practices with a shortened paddy flooding period in order to save water and mitigate methane (CH₄) emissions. However, such systems can cause an increase in N₂O emission that may partially offset the decrease in CH₄ emission (Devkota et al., 2013; Miniotti et al., 2016; Xu et al., 2015). Hence, water management practices should be improved based on a better understanding of the spatiotemporal origin of N₂O emissions and inorganic N precursors, nitrate and ammonium. Secondly, the complex hydrology and variable soil moisture conditions between soil layers and within the time course of a growing season may induce a patchwork of conditions favorable for nitrification versus denitrification versus N₂O reduction. For example, it is not clear if low N₂O emissions under more moist conditions are the result of lower N₂O production due to substrate limitation (i.e., low nitrification rates and hence low ${\mathrm{NO}}_{\mathrm{3}}^{-}$) or rather increased N₂O reduction. To date, few studies have looked at N₂O processes at depth, and it is not known how moisture and nutrient stratification affect the balance between N₂O production and consumption processes and ultimately surface emissions. Analysis of soil N₂O concentrations along a profile should help answer this. Thirdly, rice cropping systems typically suffer from a lower nitrogen use efficiency (NUE) than other major cereal crops, often attributed to high gaseous NH₃ and N₂ losses (Cassman et al., 1998; Dedatta et al., 1991; Aulakh et al., 2001; Dong et al., 2012). In improving the NUE, a better estimate of N₂O reduction to N₂ is needed to design strategies that reduce N₂ losses without increasing N₂O emission.

N₂O is predominately produced (1) as a byproduct during nitrification, where ${\mathrm{NH}}_{\mathrm{4}}^{+}$ is oxidized to ${\mathrm{NO}}_{\mathrm{3}}^{-}$ via hydroxylamine (NH₂OH) (this step of nitrification is sometimes referred to as hydroxylamine oxidation; Schreiber et al., 2012; Hu et al., 2015), or (2) as an intermediate in the denitrification pathway during which ${\mathrm{NO}}_{\mathrm{3}}^{-}$ is reduced to N₂ (Firestone et al., 1989) or (3) during nitrifier denitrification by specific ammonia-oxidizing bacteria that oxidize ${\mathrm{NH}}_{\mathrm{4}}^{+}$ to NH₂OH and then to ${\mathrm{NO}}_{\mathrm{2}}^{-}$, with a small fraction of ${\mathrm{NO}}_{\mathrm{2}}^{-}$ then being reduced to NO and N₂O (Kool et al., 2010, 2011; Wrage et al., 2001). N₂O may also be produced from additional biotic and abiotic processes, such as fungal denitrification, coupled nitrification–denitrification, dissimilatory nitrate reduction to ammonium, chemodenitrification or hydroxylamine decomposition (Butterbach-Bahl et al., 2013; Heil et al., 2015; Zhu-Barker et al., 2015). Due to the prevalence of anaerobic conditions and the use of ${\mathrm{NH}}_{\mathrm{4}}^{+}$-based fertilizers, fungal denitrification and coupled nitrification–denitrification are likely to increase in flooded rice systems. N₂O is consumed during the final step of denitrification, where N₂O is reduced to N₂ by the N₂O reductase pathway. This can occur sequentially within denitrifying organisms, or N₂O produced elsewhere from other processes or incomplete denitrification can be later reduced by denitrifiers. The final and dominant product of denitrification is N₂. While N₂ emissions are not of concern for global warming, the quantification of gross denitrification rates is of environmental concern because the loss of N via this process may represent a loss of N from the system and indicate reduced fertilizer N efficiency. Gross denitrification rates are difficult to measure in situ without the use of isotope tracers due to the high atmospheric background of N₂; thus, denitrification and N₂ emissions remain relatively unconstrained aspects of N budgets.

The measurement of N₂O isotope ratios at natural abundance is a tool to differentiate between in situ N₂O source processes and N₂O reduction (Toyoda et al., 2011; Ostrom and Ostrom, 2011; Wolf et al., 2015; Baggs, 2008), i.e., N₂O source partitioning. The evolution of analytical techniques now allows us to measure not only the bulk δ¹⁵N-N₂O, but also the intermolecular distribution of the δ¹⁵N within N₂O, called site preference (SP), and the δ¹⁵N of N₂O precursors, nitrate (${\mathrm{NO}}_{\mathrm{3}}^{-}$) and ammonium (${\mathrm{NH}}_{\mathrm{4}}^{+}$). The δ¹⁸O of N₂O and its precursors may also be used to constrain processes (Lewicka-Szczebak et al., 2016, 2017; Kool et al., 2009). Analytical methods of interpretation remain, however, only semiquantitative due to uncertainty and overlap in isotope effects (ε, η or Δ) for individual processes or cumulative processes and/or multiple N and O sources for which the determination of δ¹⁵N and δ¹⁸O remains expensive and time consuming. Theoretically, the O in N₂O derives from O₂ during nitrification and from ${\mathrm{NO}}_{\mathrm{3}}^{-}$ during denitrification or a combination during nitrifier denitrification (Kool et al., 2007, 2010; Snider et al., 2012, 2013; Lewicka-Szczebak et al., 2016). However, in the case of nitrifier denitrification and denitrification, intermediates in the reduction pathway (${\mathrm{NO}}_{\mathrm{2}}^{-}$ and NO) can extensively exchange O atoms with H₂O (Kool et al., 2007). Such exchange lowers the measured δ¹⁸O-N₂O values because of the influence of relatively depleted δ¹⁸O from H₂O, potentially leading to an underestimation of denitrification and N₂O reduction (Snider et al., 2013; Lewicka-Szczebak et al., 2016). Indeed, it has been shown that the ε¹⁸O for denitrification should be calculated relative to H₂O not ${\mathrm{NO}}_{\mathrm{3}}^{-}$, as almost 100 % O exchange occurs (Lewicka-Szczebak et al., 2014, 2016). The use of δ¹⁵N values is theoretically more straightforward, and there is also a much richer body of literature on ε¹⁵N for various processes, which was recently compiled and reviewed by Denk et al. (2017). The authors report a mean isotope effect for ¹⁵N during ${\mathrm{NH}}_{\mathrm{4}}^{+}$ oxidation to N₂O of $-\mathrm{56.6}\pm \mathrm{7.3}$ ‰ and of $-\mathrm{42.9}\pm \mathrm{6.3}$ ‰ for ${\mathrm{NO}}_{\mathrm{3}}^{-}$ reduction to N₂O. Additionally, accurate measurement of the δ¹⁵N of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ at sufficient temporal resolution remains time consuming. In comparison, the SP is thought to be independent of the initial substrate δ¹⁵N values and shows distinct values for two clusters of N₂O production, namely 32.8±4.0 ‰ for nitrification/fungal denitrification/abiotic hydroxylamine oxidation and $-\mathrm{1.6}\pm \mathrm{3.8}$ ‰ for denitrification/nitrifier denitrification (Decock and Six, 2013a; Denk et al., 2017). Abiotic N₂O production from NO has also been reported with an SP of 16 ‰ (Stanton et al., 2018).

The reduction of N₂O to N₂ enriches the pool of remaining N₂O that is measured in δ¹⁵N and δ¹⁸O and thus changes the δ¹⁵N-N₂O, δ¹⁸O-N₂O and SP (Decock and Six, 2013a; Zou et al., 2014). If the δ value of N₂O_initial (prior to reduction) can be reasonably estimated from graphical and mixing model approaches, then the subsequent enrichment of N₂O can be used to estimate N₂O reduction rates and thereby total denitrification rates. This is important because N₂O reduction is a crucial but exceptionally poorly constrained process within the N cycle (Lewicka-Szczebak et al., 2017). Fractionation during N₂O reduction may follow the dynamics of open or closed systems (Fry, 2007; Mariotti et al., 1981).

Our goal was to collect a high-resolution in situ N₂O isotope ratio dataset that could be used to (a) determine the stratification of N₂O production and reduction processes in relation to water management, (b) semiquantitatively assess N₂O and N₂ loss rates among rice water management treatments, and (c) push forward current natural abundance N₂O isotope source-partitioning methods and interpretation at the field scale. We compared three rice water management practices: direct dry seeding followed by alternate wetting and drying (DS-AWD), wet seeding followed by alternate wetting and drying (WS-AWD) and wet seeding followed by conventional flooding (WS-FLD). Isotope data were determined at three depths, simultaneously with soil environmental and nutrient data and soil N₂O and dissolved N₂O concentrations. We hypothesized that N₂O emissions would be highest in the AWD treatments due to greater contributions from nitrification and less N₂O reduction, following the order DS-AWD > WS-AWD > WS-FLD. We also hypothesized that N₂ emissions are controlled by the availability of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ coming from nitrification and high soil moisture. We considered that ${\mathrm{NO}}_{\mathrm{3}}^{-}$ would be higher under WS-AWD but soil moisture would be higher under WS-FLD; therefore, we predicted N₂ emissions to follow in the order WS-AWD > WS-FLD > DS-AWD. Lastly, we hypothesized that longer periods of lowered soil moisture in the DS-AWD and WS-AWD treatments would result in greater production of N₂O at depth and this higher production would increase surface emissions.

2 Materials and methods

2.1 Field experiment

A field experiment consisting of three water management regimes was conducted at the Italian Rice Research Center (Ente Nazionale Risi), Pavia, Italy (45^∘14^′48${}^{\prime \prime }$ N, 8^∘41^′52${}^{\prime \prime }$ E). Experimental work focused only on the early growing season, lasting from the 13 May until 30 June 2016. It is in this period that the highest N₂O losses and N cycling dynamics had been previously observed and the largest differences among water management practices occurred. The experiment was conducted in the fifth year of alternative water management in an existing experimental platform. During the first 3 years the paddies were maintained as dry-seeding + flooding, wet-seeding + flooding and intermittent irrigation as described in Miniotti et al. (2016), Peyron et al. (2016) and Said-Pullicino et al. (2016). In the fourth year, the intermittent irrigation treatment was changed to wet seeding + alternate wet–dry (Verhoeven et al., 2018). In the current study dry-seeding + flooding treatment was shifted to dry-seeding + alternate wet–dry; the other treatments remained as in the fourth year. Irrigation and water management details are provided below. The soil at the site has been classified as coarse silty, mixed, mesic Fluvaquentic Epiaquept (USDA-NRCS, 2010). The mean soil texture in the upper 30 cm of the experimental plots was 26 % sand, 62 % silt and 11 % clay with a mean bulk density of 1.29 g cm⁻³. At the end of the 2015 growing season, mean total organic C and total N were 1.07 and 0.11 % and pH 5.9 (1:2.5 H₂O) and 5.2 (1:2.5 0.01 M CaCl₂), respectively. Annual and growing season mean temperatures in 2016 were 10 and 23 ^∘C, respectively (Fig. S1 in the Supplement). Annual and growing season cumulative precipitation was 618 and 258 mm, respectively. Data for both values were retrieved from a regional weather station operated by the Agenzia Regionale per la Protezione dell'Ambiente-Lombardia, located approximately 200 m from the field site (ARPA).

Table 1 Dates of management activities during the experimental period in the three water management treatments (WS-FLD: water-seeding + conventional flooding; WS-AWD: water-seeding + alternate wetting and drying; DS-AWD: direct dry seeding + alternate wetting and drying).

Water management in the two WS treatments was identical during the first 3 weeks of the growing season (Table 1). Following regional practices for water seeding, paddies were flooded for 6 days at the time of seeding but then drained for ∼2 weeks to promote germination. During this period of “drainage”, paddies were not dry but maintained near saturation by flush irrigation as necessary (31 May and 6 June). Flush irrigation is a practice in which the water inlet channels are opened for a few hours and then the outlet channels are opened a few hours later resulting in temporary soil saturation or even 1–2 cm ponding for 2–4 h. On 10 June, approximately 3 weeks after seeding, treatment differentiation between the WS-FLD and WS-AWD began. At this time the WS-FLD was flooded, while the WS-AWD was only flush irrigated. On 16 June the WS-FLD was allowed to drain slowly in order to facilitate fertilizer application on 21 June. Following fertilizer application, the WS-FLD treatment was re-flooded and both AWD treatments were flush irrigated on 22 June. In the DS-AWD treatment, no flooding or irrigation water was applied prior to 22 June. Soil moisture depended on rainfall, which was 75 mm during the 4 weeks following seeding.

In all treatments, crop residues were incorporated in the spring, before the cropping season. All paddies were harrowed and leveled approximately 1 month prior to seeding in mid-April 2016. All treatments were pre-fertilized with phosphorus and potassium on 13 May (14 and 28 kg ha⁻¹, respectively). A total of 160 kg N ha⁻¹ as urea was applied to all treatments, with one preplant application on 16 May and two in-season applications on 21 June and 14 July (Table 1). Following best management practices for the three water management practices, a smaller preplant urea application was applied in the DS-AWD treatment, followed by a larger application in this treatment at the second and third fertilization. In the DS-AWD treatment, urea was applied at 40, 70 and 50 kg N ha⁻¹, while these rates were 60, 60 and 40 kg N ha⁻¹ for the WS treatments at fertilization 1, 2 and 3, respectively. The WS-FLD and WS-AWD treatments were seeded on 20 May. All treatments were harvested on 15 September.

Each treatment consisted of two paddies, 20 × 80 m, with two plots in each paddy, n=4 (Fig. S2). The experimental design was identical to that of Verhoeven et al. (2018), with the addition of the DS-AWD treatment and some adjustment to plot placement in order to accommodate data logging devices and field equipment. Each paddy was approximately 2 m apart and hydrologically separated by a levee of 50 cm above the soil surface, flanked by an irrigation canal on either side. Sampling for N₂O surface fluxes, pore water parameters (${\mathrm{NO}}_{\mathrm{3}}^{-}$, ${\mathrm{NH}}_{\mathrm{4}}^{+}$, dissolved organic carbon (DOC), dissolved N₂O) and pore air N₂O occurred on 15–17 dates, from the 20 May to the 30 June 2016 (Table S1 in the Supplement). Sampling dates were on average 3 days apart with a greater frequency before and after N application on the 21 June. Subsamples of pore water from 10 to 12 dates were analyzed for δ¹⁵N-${\mathrm{NO}}_{\mathrm{3}}^{-}$, δ¹⁸O-${\mathrm{NO}}_{\mathrm{3}}^{-}$ and δ¹⁵N-${\mathrm{NH}}_{\mathrm{4}}^{+}$.

2.2 Soil environment: temperature, redox potential and moisture

Soil moisture was measured using PR2 capacitance probes (Delta T Devices, UK) at 5, 15, 25, 45 and 85 cm. Water-filled pore space (WFPS) was calculated using bulk density measurements at 5, 12.5 and 25 cm collected at the beginning of the season using a Giddings manual soil auger. Soil temperature was measured in only one plot per paddy (n=2) at three depths (5, 12.5 and 25 cm). Measurements were made manually at the time of surface flux gas measurements. Soil redox potential (Eh) was measured continuously in each plot using sturdy tip probes outfitted with 5 Pt electrodes that were permanently connected to a 48-channel Hypnos-III data logger (MVH Consult, the Netherlands) with two Ag∕AgCl-reference probes. Soil Eh was measured every hour at six depths: 5, 12.5, 20, 30, 50 and 80 cm. We took the average of the 20 and 30 cm readings to derive a 25 cm reading in order to correlate it to other measurements.

2.3 N₂O measurements: surface emissions, pore air and dissolved gas

All N₂O concentration measurements were measured by gas chromatography (GC) on a Scion 456-GC (Bruker, Germany) equipped with an electron capture detector (ECD). A standard curve was derived from 10 replicates of at least five concentrations to determine the standard deviation for a given concentration. For example, the error of the GC was determined to be ±0.012 at 0.3 ppm and ±0.024 ppm at 1.0 ppm. N₂O surface emissions (N₂O_emitted) were measured by the non-steady-state closed chamber technique (Hutchinson and Mosier, 1981). The chamber design and deployment was identical to that of Verhoeven et al. (2018). Gas samples were taken at 0, 10, 20 and 30 min in each chamber and injected into pre-evacuated Exetainers (Labco, UK). At time 0 and 30 min an additional ∼170 mL of sample was taken and injected into gas crimp neck vials sealed with butyl injection stoppers (IVA Analysentechnik, Germany) to be used for isotope analysis. When the accumulation of gas over the course of measurement was less than the standard deviation associated with the highest concentration of the four measurements, the flux was determined to be below detection. Fluxes above the detection limit were calculated by linear or nonlinear regression following the method outlined by Verhoeven and Six (2014). Soil N₂O (N₂O_soil) was sampled using passive diffusion probes installed at 5, 12.5 and 25 cm. The probe design and sampling strategy have been previously described in Verhoeven et al. (2018). In brief, the samples were collected in He flushed and pre-evacuated 100 mL glass crimp neck vials (actual volume 110 mL, IVA Analysentechnik, Germany) and after sampling topped with high-purity He gas to prevent leakage into under-pressurized vials. The final N₂O concentration was determined by gas chromatography, as described above, on a subsample, while the remainder of the sample was retained for isotope analysis. The final N₂O concentration was calculated by accounting for sample dilution based on the pressure after evacuation, after sampling and after topping with He gas. Samples for dissolved N₂O (N₂O_dissolved) were collected by injecting a 5 mL subsample of pore water, collected as described in Sect. 2.4, into N₂ flushed and filled Exetainers that also contained 50 µL of 50 % ZnCl to stop microbial activity. Samples were stored at 4 ^∘C until the end of the experimental campaign and transported back to the lab for analysis; therefore, there was adequate time for the equilibration between the headspace and aqueous phases. The molar concentration of N₂O was calculated by applying the solubility constant of N₂O at the time of analysis (i.e., lab temperature) to Henry's law (Haynes, W. M. and Lide, D. R., 2011/2012; Weiss and Price, 1980; Wilhelm et al., 1977), taking into account the vial volume and headspace.

2.4 Pore water measurements

Two MacroRhizon pore water samplers (Rhizosphere Research Products, the Netherlands) were installed at each depth (5, 12.5 and 25 cm) in every plot. Pore water was then collected in two polypropylene 60 mL syringes at each depth and later pooled together at sample processing. The syringes were attached to the MacroRhizon sample tubes with two-way Luer lock valves and propped open using a wedge, which served to create a low vacuum; the syringes were left to collect water for 2–4 h. Samples were stored at 4 ^∘C and processed within 36 h. During pore water processing ∼15 mL of solution was allocated for analysis of ${\mathrm{NO}}_{\mathrm{3}}^{-}$, ${\mathrm{NH}}_{\mathrm{4}}^{+}$, δ¹⁵N, and δ¹⁸O-${\mathrm{NO}}_{\mathrm{3}}^{-}$; ∼15 mL for δ¹⁵N-${\mathrm{NH}}_{\mathrm{4}}^{+}$; 5 mL for dissolved N₂O; 3–5 mL for dissolved Fe²⁺ and Mn²⁺; and 5 mL for dissolved organic carbon and total dissolved nitrogen (DOC/TDN) analysis. All samples, aside from those for dissolved N₂O, were frozen at −5 ^∘C until analysis. ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{NH}}_{\mathrm{4}}^{+}$ were determined by spectrophotometry following the procedure of Doane and Horwáth (2003). DOC and TDN were determined by first acidifying the water sample to pH < 2 by the addition of concentrated HCl and then analysis on a multi N∕C 2100S : total organic carbon and total nitrogen (TOC/TN) Analyzer (Analytik Jena, Germany).

2.5 Determination of δ¹⁵N, δ¹⁸O and isotope ratios in N₂O_emitted and N₂O_soil

Surface and pore air gas samples were taken in 100 mL glass crimp neck vials (actual volume 110 mL; IVA Analysentechnik, Germany) as described in Sect. 2.3. Pore air gas samples were preconditioned with 1 mL of 1M NaOH solution prior to analysis due to very high CO₂ concentrations in many samples (> 5000 ppm). The intramolecular site-specific isotopic composition of the N₂O molecule was measured using a gas preparation unit (Trace Gas, Elementar, UK) coupled to an isotope ratio mass spectrometer (IRMS; IsoPrime100, Elementar, UK). The gas preparation unit was modified with an additional chemical trap ($\mathrm{1}/{\mathrm{2}}^{\prime \prime }$ diameter stainless steel), located immediately downstream from the autosampler. This pre-trap was filled with NaOH, Mg(ClO₄)₂, and activated carbon in the direction of flow and is designed to further scrub CO₂, H₂O, CO and volatile organic compounds (VOCs) which otherwise would cause mass interference during measurement. Before final injection into the IRMS, the purified gas sample is directed through a Nafion drier and subsequently separated in a gas chromatograph column (5Å molecular sieve).

The IRMS consists of five Faraday cups with m∕z of 30, 31, 44, 45 and 46, measuring δ¹⁵N and δ¹⁸O of N₂O and δ¹⁵N from the NO⁺ fragments dissociated from N₂O during ionization in the source. The ¹⁵N∕¹⁴N ratio of the NO molecule is used to calculate the α (central) position of the initial N₂O, thus allowing measurement of the site-specific isotopic composition of N₂O (SP). Site preference is defined as δ¹⁵N^SP= δ¹⁵N${}^{\mathit{\alpha }}-{{\mathit{\delta }}^{\mathrm{15}}\mathrm{N}}^{\mathit{\beta }}$ with α denoting the ¹⁵N∕¹⁴N ratio of the central N atom and β the ¹⁵N∕¹⁴N ratio of the terminal N atom of the linear NNO molecule. δ¹⁵N^β is indirectly obtained from the rearrangement of

$${{\mathit{\delta }}^{\mathrm{15}}\mathrm{N}}^{\mathrm{bulk}}=({{\mathit{\delta }}^{\mathrm{15}}\mathrm{N}}^{\mathit{\alpha }}+{\mathit{\delta }}^{\mathrm{15}}{N}^{\mathit{\beta }})/\mathrm{2},$$

which represents the average ¹⁵N content of the N₂O molecule.

For IRMS calibration three sets of two working standards (∼3 ppm N₂O mixed in synthetic air) with different isotopic composition (δ¹⁵N${}^{\mathit{\alpha }}=\mathrm{0.954}\pm \mathrm{0.123}$ ‰ and 34.446±0.179 ‰; δ¹⁵N${}^{\mathit{\beta }}=\mathrm{2.574}\pm \mathrm{0.086}$ ‰ and 35.98±0.221 ‰; δ¹⁸O = 39.741±0.051 ‰ and 38.527±0.107 ‰) were used. These standards have been analyzed at the Swiss Federal Laboratories for Materials Science and Technology (EMPA) using a trace gas extractor coupled to a quantum cascade laser absorption spectrometer versus standards with assigned δ values by the Tokyo Institute of Technology (Mohn et al., 2014). These working standards were run in triplicate, evenly spaced throughout a run. Sample peak ratios are initially reported against a N₂O reference gas peak (100 % N₂O, Carbagas, Switzerland) and are subsequently corrected for drift and span using the working standards. Further correction procedures, such as ¹⁷O mass overlap and scrambling, as reported elsewhere, were not applied as the data were inherently corrected by regression between true and measured values of the triplicate working standards. Long-term measurement quality was ensured using a control standard at low N₂O concentration (∼0.4 ppm) treated as a sample. Instrument linearity and stability were frequently checked by the injection of 10 reference gas pulses of either varying or identical height, respectively, with accepted levels of < 0.03 ‰ nA⁻¹. Since instrument linearity could only be achieved for either N₂O or NO, the instrument had been tuned for the former and δ¹⁵N^α subsequently corrected using sample peak height assuming a nonlinearity of 0.1 ‰ nA⁻¹. Such linearity complications have been previously reported using Elementar (Ostrom et al., 2007) and ThermoFinnigan IRMS (Röckmann et al., 2003). Tropospheric air was regularly measured (n=42) and used as a confirmation of correction procedures, yielding consistent and reliable results: δ¹⁵N${}^{\mathrm{SP}}=\mathrm{18.77}\pm \mathrm{1.08}$ ‰; δ¹⁵N${}^{\mathrm{bulk}}=\mathrm{5.96}\pm \mathrm{0.35}$ ‰; δ¹⁵N${}^{\mathit{\alpha }}=\mathrm{15.34}\pm \mathrm{0.70}$ ‰; δ¹⁵N${}^{\mathit{\beta }}=-\mathrm{3.43}\pm \mathrm{0.60}$ ‰; δ¹⁸O = 43.67±0.41 ‰. All ¹⁵N∕¹⁴N sample ratios are reported relative to the international isotope ratio scale AIR-N2 while ¹⁸O∕¹⁶O are reported versus Vienna Standard Mean Ocean Water (V-SMOW). Relative differences are given using the delta notation (δ) in units of ‰:

$$\begin{array}{}\text{(1)}& {\mathit{\delta }}^{Z}X\left[\mathrm{‰}\right]={\displaystyle \frac{R_{\mathrm{sample}}}{R_{\mathrm{reference}}}}-\mathrm{1},\end{array}$$

where R is referring to the molar ratio of ¹⁵N∕¹⁴N or ¹⁸O∕¹⁶O and ^ZX to the abundance of the heavy stable isotope Z of element X.

2.6 Determination of δ¹⁵N-${\mathrm{NO}}_{\mathrm{3}}^{-}$, δ¹⁸O-${\mathrm{NO}}_{\mathrm{3}}^{-}$, and δ¹⁵N-${\mathrm{NH}}_{\mathrm{4}}^{+}$

Pore water ${\mathrm{NO}}_{\mathrm{3}}^{-}$ samples were analyzed for δ¹⁵N and δ¹⁸O at the University of California, Davis, Stable Isotope Facility (SIF, UC-Davis;http://stableisotopefacility.ucdavis.edu/, last access: 15 January 2019), using the denitrifier method developed by Sigman et al. (2001), Casciotti et al. (2002), and McIlvin and Casciotti (2011). δ¹⁵N-${\mathrm{NH}}_{\mathrm{4}}^{+}$ in pore water was determined by micro-diffusion onto acidified disks followed by persulfate digestion (Lachouani et al., 2010; Stephan and Kavanagh, 2009) and lastly by the denitrifier method. For δ¹⁵N-${\mathrm{NH}}_{\mathrm{4}}^{+}$, all steps and analyses were done in-house, including the denitrifier method. Our limit of quantification for δ¹⁵N-${\mathrm{NH}}_{\mathrm{4}}^{+}$ was 0.75 mg L⁻¹ or ∼42 µM ${\mathrm{NH}}_{\mathrm{4}}^{+}$; below this the diffusion gradient was too low for reliable diffusion. Briefly, samples were run in sets of 40 with 24 samples and a combination of 16 standards and blanks. Each run contained at least two δ¹⁵N-${\mathrm{NH}}_{\mathrm{4}}^{+}$ isotope standards (IAEA N2 = 20.3 ‰; IAEA N1 = 0.4 ‰; USGS 25 = −30.4 ‰) at two or three concentrations in duplicate or triplicate in addition to two blanks and two working standards. ${\mathrm{NH}}_{\mathrm{4}}^{+}$ isotope standards were diffused, digested and run through the denitrifier method in parallel with samples and therefore an overall correction and concentration offset was derived and applied for each batch. The denitrifier method was executed using the updated protocol described by McIlvin and Casciotti (2011) using Pseudomonas aureofaciens (ATCC 13985). An IAEA KNO${}_{\mathrm{3}}^{-}$ standard (δ¹⁵N = 4.7 ‰) was included at the denitrifier method step to ensure accurate conversion of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ to N₂O. A propagated error across all steps of δ¹⁵N-${\mathrm{NH}}_{\mathrm{4}}^{+}$ quantification was calculated from the working standards included in each batch (n=18). We excluded three values that were well outside the expected range; our overall precision was 1.9 ‰. The largest sources of error were incomplete diffusion or persulfate digestion. For δ¹⁵N-${\mathrm{NO}}_{\mathrm{3}}^{-}$ and δ¹⁸O-${\mathrm{NO}}_{\mathrm{3}}^{-}$ analyzed at SIF, UC-Davis, the limit of quantification was 0.125 mg L⁻¹ ${\mathrm{NO}}_{\mathrm{3}}^{-}$ or 2.0 µM ${\mathrm{NO}}_{\mathrm{3}}^{-}$, with a precision of 0.4 ‰ and 0.5 ‰ for δ¹⁵N and δ¹⁸O, respectively.

Using N₂O_poreair and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{NH}}_{\mathrm{4}}^{+}$ in pore water we calculated the Δ¹⁵N of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ reduction to N₂O and of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ oxidation to N₂O using Eqs. (2) and (3), respectively.

$$\begin{array}{}\text{(2)}& {\displaystyle }& {\displaystyle }{\mathrm{\Delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}-{\mathrm{NO}}_{\mathrm{3}}={{\mathit{\delta }}^{\mathrm{15}}\mathrm{N}}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}{{\mathit{\delta }}^{\mathrm{15}}\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}\text{(3)}& {\displaystyle }& {\displaystyle }{\mathrm{\Delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}-{\mathrm{NH}}_{\mathrm{4}}}={\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NH}}_{\mathrm{4}}}\end{array}$$

The calculation of Δ¹⁵N_x can be compared to the net isotope effects for nitrification and denitrification-derived N₂O, as found in the literature. In reality the processes in Eqs. (1) and (2) entail a series of sequential reactions each of which has a unique isotope effect (ε_k,1, ε_k,2, ε_k,3, …). It is not possible to measure the isotope values of many of the intermediaries in these reactions series, particularly in in situ field settings; therefore, we report the Δ¹⁵N_x. For the calculation of Δ¹⁵N_x we assume open-system dynamics because all measurements were in situ where substrates, products and intermediaries could be replenished by other processes.

2.7 Determination of N₂O source contribution and N₂O reduction

Two end-member mixing models using SP and δ¹⁸O signatures: closed and open systems

We used two mixing models where N₂O reduction was modeled under “open” and “closed” system dynamics following the theory outlined originally by Fry (2007) and Mariotti et al. (1981), respectively. The two modeling methods are henceforth referred to as open and closed. In reality, the heterogeneity in microbial microhabitat within the soil most likely results in a mixture of closed versus open-system dynamics. Therefore, final data interpretations were made for the average findings across open- versus closed-system dynamics. A schematic of our closed-system approach is given in Fig. 1. For both open and closed methods, two possible scenarios were considered as described by Lewicka-Szczebak et al. (2017); scenario 1 (sc1), where N₂O is produced and reduced by denitrifiers before mixing with N₂O derived from nitrification, or scenario 2 (sc2) where N₂O is produced from both processes, mixed, and then reduced. In both models, N₂O is originally produced from two possible end-members: denitrification/nitrifier denitrification (denoted by subscript den) and nitrification/fungal denitrification (denoted by subscript nit). Our intention was to keep the derivation of end-member values consistent between this study and Lewicka-Szczebak et al. (2017). Our SP end-member values (SP_den and SP_nit) and N₂O reduction fractionation factors (ε¹⁸O_red and εSP_red) were taken directly from Lewicka-Szczebak et al. (2017) (Table 2). For δ¹⁸O-N₂O_(x) end-member values, we could not directly use the values reported in Lewicka-Szczebak et al. (2017) because these were reported relative to δ¹⁸O-H₂O (as δ¹⁸O-N₂O(N₂O/H₂O)) and we did not measure the isotope signature of water in our study. Therefore, δ¹⁸O-N₂O_nit was recalculated using the original mean values (δ¹⁸O-N₂O as opposed to δ¹⁸O-(N₂O/H₂O) of the six studies referenced by Lewicka-Szczebak et al. (2017); this yielded a mean of 36.5 ‰ (Heil et al., 2014; Sutka et al., 2006, 2008; Frame and Casciotti, 2010; Rohe et al., 2014; Maeda et al., 2015). For δ¹⁸O-N₂O_den we adjusted the value used in Lewicka-Szczebak et al. (2017) by an estimate of δ¹⁸O-H₂O of water for our site rather than recalculating it from the four studies originally referenced by Lewicka-Szczebak et al. (2017) (Lewicka-Szczebak et al., 2014, 2016; Frame and Casciotti, 2010; Sutka et al., 2006). We used a δ¹⁸O-H₂O value of −8.3 ‰, as reported by Rapti-Caputo and Martinelli (2009) for an uncontained aquifer of the Po River delta. We chose to do this because some of the mean values used in calculations by Lewicka-Szczebak et al. (2017) were themselves calculated from data originally reported.

Table 2 End-member values used for modeling of the fraction of residual N₂O not reduced (gross rN₂O) and the fraction of N₂O + N₂ attributed to denitrification (gross frac_DEN ) for both open and closed N₂O reduction fractionation dynamics.

^* Lewicka-Szczebak et al. (2017) originally report ${\mathit{\delta }}_{\mathrm{0}}^{\mathrm{18}}\mathrm{O}-{\mathrm{N}}_{\mathrm{2}}\mathrm{O}({\mathrm{N}}_{\mathrm{2}}\mathrm{O}/{\mathrm{H}}_{\mathrm{2}}\mathrm{O})$. Thus, to calculate a pure ${\mathit{\delta }}_{\mathrm{0}}^{\mathrm{18}}\mathrm{O}-{\mathrm{N}}_{\mathrm{2}}\mathrm{O}$, we added the δ¹⁸O-H₂O value used in our study (−8.3 ‰).

Figure 1 Mapping approach scheme used in the closed-system modeling. Adapted from Lewicka-Szczebak et al. (2017).

Closed-system fractionation for N₂O reduction was modeled following the method described in Lewicka-Szczebak et al. (2017) (Fig. 1). A detailed protocol for these calculations can also be found on ResearchGate (https://doi.org/10.13140/RG.2.2.17478.52804, last access: 15 January 2019). In brief, sample SP and δ¹⁸O-N₂O values are used to derive sample-specific intercepts that pass through the sample and reduction line (sc1) or the sample and the mixing line (sc2). A fixed slope for the reduction line can be calculated from εSP_red ∕ ε¹⁸O_red (i.e., in our case, $-\mathrm{5}/-\mathrm{15}$). In sc1, the intercept of the mixing and reduction line represents N₂O that has been produced from denitrification/nitrifier denitrification and partially reduced but not yet mixed with N₂O produced from nitrification/fungal denitrification. In sc2, the intercept of these lines represents N₂O that has been produced by the two end-member pools, mixed but not yet reduced. The Y axis (i.e., SP) value of these respective intercepts can be used in a generalized Rayleigh equation (Eq. 4) to calculate the extent of N₂O reduction, represented by the fraction of residual N₂O not reduced.

$$\begin{array}{}\text{(4)}& {\displaystyle }{\mathrm{SP}}_{\mathrm{resid}.\mathrm{N}\mathrm{2}\mathrm{O}}\approx {\mathrm{SP}}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}-\mathrm{unreduced}}+\mathit{\epsilon }{\mathrm{SP}}_{\mathrm{red}}\cdot \ln \left(r{\mathrm{N}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{net}}\right)\end{array}$$

In sc1 the rN₂O is determined with respect to N₂O from denitrification/nitrifier denitrification only; therefore, to calculate the residual fraction of total production (i.e., N₂+N₂O), we calculate gross rN₂O:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathrm{gross}r{\mathrm{N}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{sc}\mathrm{1}}={\displaystyle \frac{\mathrm{1}}{{\mathrm{fracDenit}}_{\mathrm{net}}/r{\mathrm{N}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{net}}+\mathrm{1}-{\mathrm{fracDenit}}_{\mathrm{net}}}}\\ \text{(5)}& {\displaystyle }& {\displaystyle }\left(\mathrm{sc}\mathrm{1},\mathrm{in}\mathrm{sc}\mathrm{2}r{\mathrm{N}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{net}}=r{\mathrm{N}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{gross}}.\right)\end{array}$$

To calculate the fraction of denitrification of the total initially produced N₂O (emitted as N₂O and N₂), we calculate the gross denitrification fraction:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\left(\mathrm{sc}\mathrm{1}\right)\mathrm{gross}{\mathrm{frac}}_{\mathrm{DEN}\mathrm{sc}\mathrm{1}-\mathrm{closed}}=\\ \text{(6)}& {\displaystyle }& {\displaystyle \frac{{\mathrm{fracDenit}}_{\mathrm{net}}/r{\mathrm{N}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{net}}}{{\mathrm{fracDenit}}_{\mathrm{net}}/r{\mathrm{N}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{net}}+\mathrm{1}-{\mathrm{fracDenit}}_{\mathrm{net}}}}.\end{array}$$

To calculate the fraction of denitrification/nitrifier denitrification to the net N₂O produced, we use Eq. (7). For simplicity and comparison with open-system calculations, we call this DenContribution.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\left(\mathrm{sc}\mathrm{1}\right)\mathrm{net}{\mathrm{frac}}_{\mathrm{DENsc}\mathrm{1}-\mathrm{closed}}={\displaystyle \frac{{\mathrm{SP}}_{\mathrm{sample}}-{\mathrm{SP}}_{\mathrm{nit}}}{{\mathrm{SP}}_{\mathrm{resid}.\mathrm{N}\mathrm{2}\mathrm{O}}-{\mathrm{SP}}_{\mathrm{nit}}}}\\ \text{(7)}& {\displaystyle }& {\displaystyle }={\mathrm{DenContribution}}_{\mathrm{closed}-\mathrm{sc}\mathrm{1}}\end{array}$$

In this case, SP_resid.N2O is the signature of residual bacterial N₂O after partial reduction but before mixing. This was determined from the graphical method (Lewicka-Szczebak et al., 2017). In sc2 both net and gross fractions of denitrification are equal and can be expressed as

$$\begin{array}{}\text{(8)}& {\displaystyle }\left(\mathrm{sc}\mathrm{2}\right){\mathrm{DenContribution}}_{\mathrm{closed}-\mathrm{sc}\mathrm{2}}={\displaystyle \frac{{\mathrm{SP}}_{\mathrm{N}\mathrm{2}\mathrm{O}-\mathrm{unreduced}}-{\mathrm{SP}}_{\mathrm{nit}}}{{\mathrm{SP}}_{\mathrm{den}}-{\mathrm{SP}}_{\mathrm{nit}}}}.\end{array}$$

Here, SP_{N2O−undreduced} is the signature of N₂O mixed from nitrification/fungal denitrification and denitrification/nitrifier denitrification but before reduction. This was determined from the graphical method (Lewicka-Szczebak et al., 2017).

To predict rN₂O in open systems, we set up a series of mass balance equations using our measured N₂O flux or N₂O_poreair concentrations and measured δ¹⁸O and SP values. We used the same end-member values listed in Table 2 for all equations. As above, we can model the interaction between mixing and reduction assuming sc1 (Eqs. 9–11) or sc2 (Eqs. 9, 12, 13). In Eqs. (9)–(13), we use k_nit, k_den and k_red to represent the gross process rates or concentrations of N₂O attributable to nitrification, denitrification and N₂O reduction, respectively.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathrm{N}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{flux}}\left(\mathrm{or}{\mathrm{N}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{poreair}}\right)=k_{\mathrm{nit}}+k_{\mathrm{den}}-k_{\mathrm{red}}\mathrm{note}:k_{\mathrm{den}}\\ \text{(9)}& {\displaystyle }& {\displaystyle }=\mathrm{total}\mathrm{denitrification}({\mathrm{N}}_{\mathrm{2}}\mathrm{O}+{\mathrm{N}}_{\mathrm{2}})\text{(10)}& {\displaystyle }& {\displaystyle }\left(\mathrm{sc}\mathrm{1}\right)\mathrm{SP}-{\mathrm{N}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{measured}}={\displaystyle \frac{\left\{\begin{array}{c}{\mathrm{SP}}_{\mathrm{nit}}{k}_{\mathrm{nit}}\\ +\left({\mathrm{SP}}_{\mathrm{den}}-\mathit{\epsilon }{\mathrm{SP}}_{\mathrm{red}}\left(\frac{k_{\mathrm{red}}}{k_{\mathrm{den}}}\right)\right)\\ \left(k_{\mathrm{den}}-k_{\mathrm{red}}\right)\end{array}\right\}}{k_{\mathrm{nit}}+k_{\mathrm{den}}-k_{\mathrm{red}}}}{\displaystyle }& {\displaystyle }\left(\mathrm{sc}\mathrm{1}\right){\mathit{\delta }}^{\mathrm{18}}\mathrm{O}-{\mathrm{N}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{measured}}\\ \text{(11)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{\left\{\begin{array}{c}\left({\mathit{\delta }}^{\mathrm{18}}{{\mathrm{ON}}_{\mathrm{2}}\mathrm{O}}_{\mathrm{nit}}\right)k_{\mathrm{nit}}\\ +\left({\mathit{\delta }}^{\mathrm{18}}{{\mathrm{ON}}_{\mathrm{2}}\mathrm{O}}_{\mathrm{den}}-{{\mathit{\epsilon }}^{\mathrm{18}}\mathrm{O}}_{\mathrm{red}}\left(\frac{k_{\mathrm{red}}}{k_{\mathrm{den}}}\right)\right)\\ \left(k_{\mathrm{den}}-k_{\mathrm{red}}\right)\end{array}\right\}}{k_{\mathrm{nit}}+k_{\mathrm{den}}-k_{\mathrm{red}}}}{\displaystyle }& {\displaystyle }\left(\mathrm{sc}\mathrm{2}\right)\mathrm{SP}-{\mathrm{N}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{measured}}={\displaystyle \frac{\left({\mathrm{SP}}_{\mathrm{nit}}{k}_{\mathrm{nit}}+{\mathrm{SP}}_{\mathrm{den}}{k}_{\mathrm{den}}\right)}{k_{\mathrm{nit}}+k_{\mathrm{den}}}}\\ \text{(12)}& {\displaystyle }& {\displaystyle }-\mathit{\epsilon }{\mathrm{SP}}_{\mathrm{red}}\left(\mathrm{1}-{\displaystyle \frac{{\mathrm{N}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{flux}}}{k_{\mathrm{nit}}+k_{\mathrm{den}}}}\right){\displaystyle }& {\displaystyle }\left(\mathrm{sc}\mathrm{2}\right){\mathit{\delta }}^{\mathrm{18}}\mathrm{O}-{\mathrm{N}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{measured}}\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{\left({\mathit{\delta }}^{\mathrm{18}}{{\mathrm{ON}}_{\mathrm{2}}\mathrm{O}}_{\mathrm{nit}}\right)k_{\mathrm{nit}}+\left({\mathit{\delta }}^{\mathrm{18}}{{\mathrm{ON}}_{\mathrm{2}}\mathrm{O}}_{\mathrm{den}}\right)k_{\mathrm{den}}}{k_{\mathrm{nit}}+k_{\mathrm{den}}}}\\ \text{(13)}& {\displaystyle }& {\displaystyle }-{{\mathit{\epsilon }}^{\mathrm{18}}\mathrm{O}}_{\mathrm{red}}\left(\mathrm{1}-{\displaystyle \frac{{\mathrm{N}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{flux}}}{k_{\mathrm{nit}}+k_{\mathrm{den}}}}\right)\end{array}$$

These two sets of equations (Eq. 9, 10, 11) or (Eq. 9, 12, 13), representing each scenario, were applied to measured surface fluxes to produce process rates in g N₂O-N ha⁻¹ d⁻¹ or were applied to N₂O_poreair concentrations to produce concentrations of N₂O in µg N₂O-N L⁻¹. By rearranging these process rates or concentrations, we can calculate gross rN₂O, frac_DEN and the contribution of denitrification to N₂O using Eqs. (14)–(16).

$$\begin{array}{ll}\text{(14)}& {\displaystyle }& {\displaystyle }\mathrm{gross}{\mathrm{frac}}_{\mathrm{DEN}\mathrm{sc}\mathrm{1},\mathrm{sc}\mathrm{2}-\mathrm{open}}={\displaystyle \frac{k_{\mathrm{den}}}{k_{\mathrm{nit}}+k_{\mathrm{den}}}}\text{(15)}& {\displaystyle }& {\displaystyle }\mathrm{gross}r{\mathrm{N}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{sc}\mathrm{1},\mathrm{sc}\mathrm{2}-\mathrm{open}}={\displaystyle \frac{k_{\mathrm{nit}}+k_{\mathrm{den}}-k_{\mathrm{red}}}{k_{\mathrm{nit}}+k_{\mathrm{den}}}}{\displaystyle }& {\displaystyle }{\mathrm{DenContribution}}_{\mathrm{sc}\mathrm{1},\mathrm{sc}\mathrm{1}-\mathrm{open}}={\displaystyle \frac{\left(k_{\mathrm{den}}-k_{\mathrm{red}}\right)}{\left[{\mathrm{N}}_{\mathrm{2}}\mathrm{O}\right]}},\left[{\mathrm{N}}_{\mathrm{2}}\mathrm{O}\right]\\ \text{(16)}& {\displaystyle }& {\displaystyle }={{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}_{\mathrm{flux}}\mathrm{or}{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}_{\mathrm{poreair}}\end{array}$$

Plausible solutions for k_red, k_den and k_red were estimated based on minimizing the sum of squares between the modeled and measured N₂O flux (or concentration), δ¹⁸O and SP values using a generalized reduced gradient (GRG) nonlinear algorithm in the Solver function of Excel. Example calculations for the open-system modeling are given in the Supplement. Solutions with a minimum sum of squares over 500 were considered implausible (8.3 % of solutions) (Table S2). Both models produced some non-plausible solutions, i.e., fractional contributions over 1 or under 0. Only solutions with a gross rN₂O, gross frac_DEN and DenContribution between 0 to 1 and an open-system minimum sum of squares < 500 were retained. In sc1, roughly 75 % of solutions met these criteria. For sc2, less than 10 % of solutions in the open system met this criterion; therefore, we did not proceed to analyze and discuss solutions from sc2 (Table S2 and Fig. S3).

2.8 Statistical analyses

Response variables were analyzed using a linear mixed effects analysis of covariance (ANCOVA) model with treatment, date and depth (if applicable) as fixed effects and plot as a random effect. The longitudinal position in the field (Y position) measured in meters from the central driveway (Fig. S2), was used as a covariate to account for potential heterogeneity in the longitudinal direction. In the case of non-normally distributed data, data were transformed to obtain a normal distribution of residuals. Due to the non-normal distribution of many variables, Spearman correlations were used to analyze the relationship between N₂O_emitted fluxes, isotope ratios, and soil environmental and substrate variables. Post hoc analysis of treatment and depth within a given day was performed using the lsmeans function with a Tukey adjustment for multiple comparisons. For the analysis of modeling results we eliminated the 25 cm depth due to poor data availability. All data analysis was done in R version 3.3.2.

3 Results

3.1 Yield

At the end of the growing season yield was measured in the larger plots in which sampling plots were situated. The DS-AWD treatment had a significantly lower yield, 6.6 t ha⁻¹, relative to 8.9 and 8.2 t ha⁻¹ in the WS-FLD and WS-AWD, respectively (Table 1).

3.2 N₂O fluxes and dissolved and pore air N₂O concentrations

3.2.1 Temporal patterns in N₂O fluxes and concentrations

After the first basal fertilization (16 May) and prior to the second topdressing fertilization (21 June), emissions were significantly higher in the DS-AWD treatment than in WS-AWD and WS-FLD on 8 and 6 of the 11 sampling days, respectively (Fig. 2). During this time four peaks in emissions were observed in the DS-AWD treatment: on 20 May, 1–3, 7–9, and 20 June, averaging 39.5±5.1 g N₂O-N ha⁻¹ d⁻¹. A peak in emissions following the second fertilization (21 June) was observed in all treatments; in the DS-AWD treatment, emissions peaked at 108.2±4.2 g N₂O-N  ha⁻¹ d⁻¹ on 23 June while in the WS-AWD and WS-FLD treatments, emissions peaked 1 day earlier reaching 49.4±17.9 and 77.67±10.6 g N₂O-N ha⁻¹ d⁻¹, respectively. In the WS-AWD treatment, emissions remained slightly elevated following this fertilization until the end of the monitoring campaign, while in the DS-AWD and WS-FLD, emissions declined after 22 or 23 June, respectively.

Figure 2 N₂O surface emissions, log₁₀ of dissolved and pore air N₂O concentrations, ${\mathrm{NH}}_{\mathrm{4}}^{+}$, ${\mathrm{NO}}_{\mathrm{3}}^{-}$, DOC, Eh, and WFPS throughout the field measurement period in the three water management treatments (WS-FLD: water-seeding + conventional flooding; WS-AWD: water-seeding + alternate wetting and drying; DS-AWD: direct dry seeding + alternate wetting and drying). The dashed vertical line indicates the date of fertilization (60 kg urea-N ha⁻¹). Blue shaded areas represent periods of flooding; shaded areas that last only 1 day indicate “flush irrigation”, i.e., flooding for < 6 h. The error bars represent the standard error of the mean. Red and orange dashed vertical lines represent the date of seeding and fertilization in each treatment, respectively.

If we exclude N₂O_dissolved measurements from the DS-AWD treatment following the second fertilization (i.e., after the 22 June, when concentrations reached as high as 594.4±112.6 µg N₂O-N L⁻¹ at 5 cm), concentrations throughout the profile of all treatments remained under 20 µg N₂O-N L⁻¹. Due to the large differences between dates and treatments, we present the concentrations on a log₁₀ scale (Fig. 2) and a non-transformed scale (Fig. S4). Peak concentrations in the WS treatments occurred at 5 cm on the first day of measurement, reaching 17.7±5.1 and 18.5±2.8 µg N₂O-N L⁻¹ in the WS-AWD and WS-FLD, respectively. In comparison, in the DS-AWD treatment, peak concentrations prior to the second fertilization were observed at 25 cm on 3 June, reaching 18.5±8.3 µg N₂O-N L⁻¹.

As with dissolved N₂O, pore air N₂O concentrations were highly variable between treatments and between sampling days and are again presented on a log₁₀ scale (Fig. 2) and non-transformed scale (Fig. S4). In both WS treatments, the highest concentrations were observed on the first day of measurement, 20 May, reaching 2903.3±1103.6 and 1321±998.0 µg N₂O-N L⁻¹ at 5 cm in the WS-FLD and WS-AWD, respectively. Elevated concentrations of N₂O_poreair were also observed in the DS-AWD on the first day of measurement but were 70.1 µg N₂O-N L⁻¹ at 5 cm (roughly 40× lower than in WS-FLD on this date). Maximum concentrations in the DS-AWD treatment were observed 2 days after the second fertilizer application, reaching 1902.2 µg N₂O-N L⁻¹; in contrast, no change was observed in the WS treatments following this fertilizer application. In all treatments the majority of N₂O_poreair concentrations were orders of magnitude lower than these peaks. There was a tendency of lower N₂O_poreair concentrations in the DS-AWD treatment relative to the WS treatments; this pattern was most evident at 5 cm (Fig. 2). However, treatment differences in N₂O_poreair were not significant (p=0.08, Table S3), and there was a significant date × treatment interaction.

Table 3 Spearman correlations of N₂O_emitted with N₂O_emitted isotope ratios, N₂O driving variables and N₂O_poreair isotope ratios measured at 5 cm in the three water management treatments (WS-FLD: water-seeding + conventional flooding; WS-AWD: water-seeding + alternate wetting and drying; DS-AWD: direct dry seeding + alternate wetting and drying). Significance indicated by ${}^{****}$ < 0.0001, ${}^{***}$ < 0.001, ${}^{**}$ < 0.01, ^* < 0.05.

3.2.2 Relation of N₂O fluxes and concentrations with soil environment, substrates and N₂O isotope ratios

We evaluated the correlation of N₂O_emitted with Eh, WFPS, ${\mathrm{NO}}_{\mathrm{3}}^{-}$, ${\mathrm{NH}}_{\mathrm{4}}^{+}$, dissolved and pore air N₂O concentrations, and N₂O isotope ratios at 5 cm (Table 3). Among these variables, N₂O emissions in the WS treatments were negatively correlated with pore water ${\mathrm{NH}}_{\mathrm{4}}^{+}$ and DOC in the WS-AWD treatment. In the DS-AWD treatment, emissions positively correlated with N₂O_poreair, WFPS and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and negatively with N₂O isotope ratios. Examining the isotope ratios of N₂O_emitted, we observed that N₂O_emitted was negatively correlated with δ¹⁸O-N₂O_emitted in all treatments, negatively with δ¹⁵N-N₂O_emitted in the DS-AWD treatment, and negatively with SP-N₂O_emitted in the WS-FLD and DS-AWD. Interestingly, a positive correlation between N₂O_emitted and SP-N₂O_emitted was observed in the WS-AWD treatment. Relative to the DS-AWD, the WS treatments had fewer significant correlations between N₂O isotope ratios, soil environment or pore air N₂O isotope ratios. DOC was positively correlated with δ¹⁵N-N₂O_emitted in the WS-AWD and with δ¹⁸O-N₂O_emitted in the WS-FLD. SP-N₂O_emitted was positively correlated to Eh and negatively to WFPS in the WS-AWD treatment. In comparison, in the DS-AWD treatment, N₂O_emitted isotope ratios were positively correlated to that of N₂O_poreair for all three isotopes. Furthermore, N₂O isotope ratios in the DS-AWD treatment were negatively correlated with N₂O_poreair concentrations, WFPS, ${\mathrm{NO}}_{\mathrm{3}}^{-}$ (δ¹⁵N-N₂O_emitted only) and N₂O_dissolved (δ¹⁸O-N₂O_emitted and SP-N₂O_emitted only). It should be noted that N₂O_dissolved in the DS-AWD treatment was not measurable at the 5 cm depth in 10 of the 16 sampling dates due to low soil moisture and low pore water volumes.

3.3 Spatiotemporal patterns of N₂O isotope ratios

3.3.1 δ¹⁵N-N₂O

A consistent temporal pattern of higher N₂O_poreair concentrations and N₂O_emitted fluxes in association with lower δ¹⁵N was observed in the DS-AWD treatment. In the WS treatments, high N₂O_emitted fluxes on 23 June following the second fertilization were associated with lower δ¹⁵N (Fig. 3); this was not the case for a high flux in the WS-AWD on 17 June. N₂O_poreair at 5 cm in the WS-AWD treatment tended to be higher in concentration and lower in δ¹⁵N relative to other depths; however, in general a consistent relationship between concentration and δ¹⁵N was less evident in the two WS treatments. On average, the δ¹⁵N of N₂O_emitted was lower relative to N₂O_poreair in the DS-AWD treatment. In contrast, in the WS treatments N₂O_emitted was depleted in ¹⁵N relative to N₂O_poreair at all depths only immediately before and after the second fertilization. In these treatments, δ¹⁵N-N₂O_poreair was generally lower at 5 cm relative the other depths but tended to increase and reach similar values as the other depths over the experimental period. As a result, N₂O_emitted was often enriched in ¹⁵N relative to N₂O_poreair at 5 cm in these treatments, particularly in the WS-AWD treatment.

Figure 3 Time course of δ¹⁵N-N₂O, δ¹⁸O-N₂O and SP-N₂O in N₂O_emitted and N₂O_poreair across the three depths and water management treatments (WS-FLD: water-seeding + conventional flooding; WS-AWD: water-seeding + alternate wetting and drying; DS-AWD: direct dry seeding + alternate wetting and drying). The error bars represent the standard error of the mean. Red and orange dashed vertical lines represent the date of seeding and fertilization in each treatment, respectively.

3.3.2 δ¹⁸O-N₂O

As with δ¹⁵N, δ¹⁸O isotope ratios spanned a large range, particularly in the emitted N₂O (Fig. 3). δ¹⁸O-N₂O_poreair in the DS-AWD followed a temporal pattern similar to δ¹⁵N,s and similarly, δ¹⁸O was generally lower in N₂O_emitted relative to N₂O_poreair. The highest δ¹⁸O-N₂O_poreair was seen in the DS-AWD treatment at moderate N₂O_poreair concentrations where δ¹⁸O isotope ratios were higher than other concentrations in the DS-AWD or any concentration in the WS treatments. These samples were also nearly always taken from 12.5 or 25 cm. In all treatments, lower δ¹⁸O values were observed in N₂O_poreair and N₂O_emitted on the first day of sampling: global mean of 35.1±1.1 and 29.6±1.7 ‰ relative to 46.9±0.4 and 43.9±1.7 ‰, respectively. Otherwise, no distinct pattern with depth, time or concentration was observed in the WS treatments.

3.3.3 SP-N₂O

The SP of N₂O_emitted ranged from 4.5±0.4 ‰ to 25.6±8.1 ‰, from 2.9±1.0 ‰ to 37.2 ‰ (un-replicated) and from 5.8±0.6 ‰ to 40.6±12.4 ‰, in the DS-AWD, WS-AWD and WS-FLD treatments, respectively (Fig. 3). In contrast to δ¹⁵N and δ¹⁸O isotope ratios, the SP-N₂O_poreair tended to increase with time but only in the WS treatments. As with δ¹⁵N-N₂O and δ¹⁸O-N₂O, moderate- and lower-concentration N₂O_poreair samples showed higher SP values relative to higher-concentration N₂O_poreair samples. For example, 2 days after the second fertilizer application (23 June), SP values decreased in conjunction with increased N₂O_poreair concentrations in the DS-AWD treatment. On this date, mean SP values at 5 cm demonstrated the largest treatment differences with values of 0.7±4.5 ‰, 27.6±2.1 ‰ and 39.9±2.7 ‰ in the DS-AWD, WS-AWD and WS-FLD treatments, respectively. On this date, the pattern between the treatments was consistent throughout the three depths.

3.3.4 Relationships between N₂O isotope ratios

Considering all depths and emitted data together, δ¹⁸O-N₂O significantly and positively correlated with δ¹⁵N-N₂O and SP across all treatments. The slope of δ¹⁸O-N₂O vs. δ¹⁵N-N₂O was 0.67, 0.28 and 0.52 (Fig. S5) and 0.67, 0.54 and 0.31 for SP vs. δ¹⁸O-N₂O in the DS-AWD, WS-AWD and WS-FLD treatments, respectively (Fig. 4a). There was no correlation between SP and δ¹⁵N-N₂O in the two WS treatments, but a positive correlation for the DS-AWD was found, with a slope of 0.62 (Fig. 4b). Examining these relationships by depth, we saw the strongest relationship and highest slope in the N₂O_emitted and at 25 cm for δ¹⁸O-N₂O vs. δ¹⁵N-N₂O (Fig. S5). While the SP vs. δ¹⁸O-N₂O showed no correlation among the surface fluxes in the WS treatments, the two isotope ratios were positively correlated in N₂O_poreair at all depths and treatments (Fig. S6). A contrasting relationship between SP and δ¹⁵N-N₂O was observed for the WS-FLD treatment in the N₂O_emitted and N₂O_poreair where the two isotope ratios were negatively correlated in N₂O_emitted and positively in N₂O_poreair (Fig. S7).

Figure 4 Graphical two end-member mixing plot after Lewicka-Szczebak et al. (2017) where sample values are plotted in SP × δ¹⁸O−N₂O space A and two-end mixing plot after Toyoda et al. (2011) where sample values are plotted in SP × δ¹⁵N-N₂O space B. In panel (a) the black dots indicate the mean literature end-member values used in our modeling scenarios and the boxes represent a range of values derived from the literature attributed to each process; see Sect. 2.7 and Table 2. To calculate the range of N₂O potentially produced by nitrification or denitrification in (b), we used the mean isotope effects, ${\mathit{\epsilon }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}/{\mathrm{NO}}_{\mathrm{3}}}$ and ${\mathit{\epsilon }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}/{\mathrm{NH}}_{\mathrm{4}}}$, reported in Denk et al. (2017) to represent denitrification and nitrification derived N₂O, respectively, and then added the minimum and maximum δ¹⁵N-${\mathrm{NO}}_{\mathrm{3}}^{-}$ and δ¹⁵N-${\mathrm{NH}}_{\mathrm{4}}^{+}$ values observed in each treatment (Table S1.4). The linear relationship between each isotope pair is indicated in italics for all points together and for N₂O_poreair only. The three water management treatments were as follows: WS-FLD – water-seeding + conventional flooding; WS-AWD – water-seeding + alternate wetting and drying; DS-AWD – direct dry seeding + alternate wetting and drying.

3.4 ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{NH}}_{\mathrm{4}}^{+}$ concentrations and isotope ratios

3.4.1 Spatiotemporal trend in ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{NH}}_{\mathrm{4}}^{+}$ concentration and δ¹⁵N and δ¹⁸O isotope ratios

In all treatments, pore water ${\mathrm{NH}}_{\mathrm{4}}^{+}$ concentrations were highest at 5 cm relative to the other depths (Fig. 2). In the DS-AWD treatment concentrations were almost zero prior to the second fertilization, remaining below 0.85 mg ${\mathrm{NH}}_{\mathrm{4}}^{+}$-N L⁻¹ across all depths. Following this fertilization, concentrations increased at all depths, most notably at 5 cm. An opposing pattern was observed in the WS treatments where ${\mathrm{NH}}_{\mathrm{4}}^{+}$ was nearly always significantly higher than in DS-AWD for each corresponding depth leading up to the second fertilization but dropped to near zero following the fertilization. Nitrate concentrations were exclusively less than 1.5 mg NO₃-N L⁻¹ in both WS treatments throughout the experimental period. In sharp contrast, ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentrations in the DS-AWD were at times more than 75 times higher than in WS treatments, peaking on 1 June at 113.6±22.4 mg NO₃-N L⁻¹. Following this spike, concentrations steadily declined and dropped to zero following the second fertilization.

3.4.2 δ¹⁵N-${\mathrm{NO}}_{\mathrm{3}}^{-}$, δ¹⁵N-${\mathrm{NH}}_{\mathrm{4}}^{+}$ and isotope enrichment factors: ${\mathrm{\Delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}/{\mathrm{NO}}_{\mathrm{3}}}$ and ${\mathrm{\Delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}/{\mathrm{NH}}_{\mathrm{4}}}$

Concentrations of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ or ${\mathrm{NH}}_{\mathrm{4}}^{+}$ were often too low for isotope measurements. Hence, we could only obtain sufficient replication for statistical analysis across depths and treatments on 5 days for ${\mathrm{NO}}_{\mathrm{3}}^{-}$ (24 and 27 May; 1, 14 and 23 June) and 2 days for ${\mathrm{NH}}_{\mathrm{4}}^{+}$ (24 May and 23 June) (Fig. S9). Daily mean δ¹⁵N-${\mathrm{NO}}_{\mathrm{3}}^{-}$ ranged from −4.3 ‰ to 28.3 ‰ across all treatments and depths. In the DS-AWD treatment a consistent depth pattern was observed with ¹⁵N enrichment of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ at 25 cm > 12.5 cm = 5 cm. δ¹⁵N-${\mathrm{NO}}_{\mathrm{3}}^{-}$ increased with time at 5 cm, rising from $-\mathrm{4.3}\pm \mathrm{1.5}$ ‰ to 22.0±4.9 ‰. Significant treatment and depth differences were observed on 24 and 27 May and 1 June, but no differences were observed on later dates (14 or 23 June). Following the second fertilizer application, δ¹⁵N-${\mathrm{NO}}_{\mathrm{3}}^{-}$ values in the DS-AWD treatment rose by approximately 10 ‰ at all depths. Daily mean δ¹⁵N-${\mathrm{NH}}_{\mathrm{4}}^{+}$ ranged from −6 ‰ to 15.2 ‰ (Fig. S9). Averaging across the experimental period and depths, mean δ¹⁵N values of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{NH}}_{\mathrm{4}}^{+}$ were similar (8.4 ‰ and 7.0 ‰, respectively; Table S5). There was no evident temporal or depth trend in δ¹⁵N-${\mathrm{NH}}_{\mathrm{4}}^{+}$ in any of the treatments. The only significant difference was lower δ¹⁵N-${\mathrm{NH}}_{\mathrm{4}}^{+}$ in the DS-AWD on 23 June. δ¹⁵N-${\mathrm{NO}}_{\mathrm{3}}^{-}$ values positively correlated to N₂O_poreair concentrations in the DS-AWD and WS-FLD treatments and were negatively correlated to ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentrations and to δ¹⁵N-${\mathrm{NH}}_{\mathrm{4}}^{+}$ in the DS-AWD treatment (Table 4). δ¹⁵N-${\mathrm{NH}}_{\mathrm{4}}^{+}$ was negatively correlated to N₂O_poreair concentrations and ${\mathrm{NH}}_{\mathrm{4}}^{+}$ concentrations and positively to δ¹⁵N-N₂O_poreair in the DS-AWD treatment.

Table 4 Spearman correlations between δ¹⁵N-${\mathrm{NO}}_{\mathrm{3}}^{-}$ and δ¹⁵N-${\mathrm{NH}}_{\mathrm{4}}^{+}$ with N₂O_poreair concentration, δ¹⁵N-N₂O_poreair, ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{NH}}_{\mathrm{4}}^{+}$ concentrations in the three water management treatments (WS-FLD: water-seeding + conventional flooding; WS-AWD: water-seeding + alternate wetting and drying; DS-AWD: direct dry seeding + alternate wetting and drying).

Largely reflecting the depth pattern of δ¹⁵N-${\mathrm{NO}}_{\mathrm{3}}^{-}$ in the DS-AWD, the calculated ${\mathrm{\Delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}/{\mathrm{NO}}_{\mathrm{3}}}$ tended to be highest at 5 cm, with a mean of $-\mathrm{7.2}\pm \mathrm{2.7}$ ‰, while mean values at 12.5 and 25 cm were slightly lower ($-\mathrm{9.5}\pm \mathrm{2.0}$ and $-\mathrm{16.0}\pm \mathrm{2.1}$ ‰, respectively; Fig. S9). At 5 cm ${\mathrm{\Delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}/{\mathrm{NO}}_{\mathrm{3}}}$ values in the DS-AWD were significantly higher than in the WS treatments; at 12.5 cm they tended to be higher as well but the difference was not significant. Two days after the second fertilizer application, the ${\mathrm{\Delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}/{\mathrm{NO}}_{\mathrm{3}}}$ in the DS-AWD markedly decreased at all depths to a treatment mean of $-\mathrm{23.6}\pm \mathrm{2.6}$ ‰. In comparison, WS treatment ${\mathrm{\Delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}/{\mathrm{NO}}_{\mathrm{3}}}$ values rose 1 (WS-FLD) or 2 (WS-AWD) days following the fertilization. In the WS-FLD, the increase in ${\mathrm{\Delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}/{\mathrm{NO}}_{\mathrm{3}}}$ values lasted only 1 day; unfortunately, low ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentrations precluded δ¹⁵N-${\mathrm{NO}}_{\mathrm{3}}^{-}$ analysis on many dates making temporal patterns difficult to observe. Mean depths by treatment isotope effects calculated relative to δ¹⁵N-${\mathrm{NH}}_{\mathrm{4}}^{+}$ (${\mathrm{\Delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}/{\mathrm{NH}}_{\mathrm{4}}}$) were $-\mathrm{12.7}\pm \mathrm{3.2}$ ‰, $-\mathrm{24.5}\pm \mathrm{2.6}$ ‰ and $-\mathrm{20.6}\pm \mathrm{2.2}$ ‰ at 5 cm; $-\mathrm{9.9}\pm \mathrm{4.0}$ ‰, $-\mathrm{12.8}\pm \mathrm{2.8}$ ‰ and $-\mathrm{15.9}\pm \mathrm{1.9}$ ‰ at 12.5 cm; $-\mathrm{17.0}\pm \mathrm{5.9}$ ‰, $-\mathrm{6.4}\pm \mathrm{1.7}$ ‰ and $-\mathrm{5.8}\pm \mathrm{2.7}$ ‰ at 25 cm for DS-AWD, WS-AWD and WS-FLD treatments, respectively. Data for ${\mathrm{\Delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}/{\mathrm{NH}}_{\mathrm{4}}}$ was scarce in the DS-AWD treatment due to low ${\mathrm{NH}}_{\mathrm{4}}^{+}$ concentrations; in the WS treatments ${\mathrm{\Delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}/{\mathrm{NH}}_{\mathrm{4}}}$ increased with depth, but these differences were not significant.

δ¹⁸O-${\mathrm{NO}}_{\mathrm{3}}^{-}$ was significantly depleted in the DS-AWD treatment relative to both WS treatments (Fig. S9). Prior to the second fertilization, values were remarkably consistent in the DS-AWD at all depths, ranging from 0.1 ‰ to 7.5 ‰. Two days after this fertilizer application, δ¹⁸O-${\mathrm{NO}}_{\mathrm{3}}^{-}$ rose to a mean of 7.6 ‰ across depths. In comparison the ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}-{\mathrm{NO}}_{\mathrm{3}}^{-}$ of both WS treatments was more variable between sampling dates, fluctuating between 12.2 ‰ and 38.8 ‰ and 10.4 ‰ and 32.7 ‰ leading up the second fertilization in the WS-AWD and WS-FLD, respectively. Two days after the second fertilizer application values rose to a mean of 23.7 ‰ and 27.4 ‰ across depths in the WS-AWD and WS-FLD, respectively. We calculated the net isotope effect for δ¹⁸O-N₂O relative to water (${\mathrm{\Delta }}^{\mathrm{18}}{\mathrm{O}}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}/{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$). The ${\mathrm{\Delta }}^{\mathrm{18}}{\mathrm{O}}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}/{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ in all treatments and depths tended to rise over the course of the measurement period, with the most consistent rise observed at 5 cm. Here values rose from a global mean of 43.8±1.0 ‰ on 20 May to 58.5±1.0 ‰ on 30 June. There was a pattern of higher ${\mathrm{\Delta }}^{\mathrm{18}}{\mathrm{O}}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}/{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ in the DS-AWD treatment relative to the two WS treatments. A drop in ${\mathrm{\Delta }}^{\mathrm{18}}{\mathrm{O}}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}/{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ of ∼10 ‰ was observed in all depths on 23 June, 2 days after the second fertilization with urea, in the DS-AWD only.

3.5 SP × δ¹⁸O-N₂O two end-member mixing model to estimate N₂O reduction, source contributions and N₂O reduction

To further quantitatively interpret our isotope ratio data, we employed a graphical two end-member mixing model (Lewicka-Szczebak et al., 2017), based on the relationship between SP and δ¹⁸O-N₂O (Figs. 1 and 4). Data were modeled for open and closed fractionation dynamics under two scenarios. In sc1, reduction of N₂O from the denitrification/nitrifier-denitrification end-member pool occurs prior to mixing with nitrification/fungal denitrification-derived N₂O; in sc2, mixing of N₂O from both end-member pools occurs before reduction. For sc2, our model yielded implausible results for the contribution of denitrification/nitrifier denitrification to N₂O emissions in about 90 % and 20 % of observations under open- and closed-system dynamics, respectively (Table S2). The poorer outcomes from sc2 in the open system indicate that the assumptions underlying this scenario are likely false in open systems or vice versa. In order to have comparable data between open and closed systems, we discuss only results coming from sc1 simulations.

Figure 5 Modeled denitrification/nitrifier-denitrification contribution and gross rN₂O of open (grey bars), closed (blue bars) and mean (purple points and line) systems predicted by a two end-member mixing model using δ¹⁸O-N₂O and SP values. For open- and closed-system dynamics, the shaded bars represent the standard deviation range for each treatment × depth combination. The purple error bars represent the standard deviation around the mean. Red and orange dashed vertical lines represent the date of seeding and fertilization in each treatment, respectively.

Temporal trends in the gross rates of rN₂O (extent of N₂O reduction) predicted by open- and closed-system N₂O fractionation were nearly identical (Fig. 5b). Gross rN₂O was estimated to be higher (i.e., lower N₂O reduction) under closed-system fractionation dynamics. In reality, it can be assumed that neither perfect open or closed systems exist in nature and processes likely reflect a mixture of these dynamics. The use of one or the other case may bias results; therefore, we chose to take the mean of the two systems to estimate N₂O reduction, nitrification/fungal denitrification and denitrification/nitrifier denitrification-derived N₂O emissions (Decock and Six, 2013b; Wu et al., 2016). Due to a disproportionate number of missing values at 25 cm in the two WS treatments, we chose not to include data from this depth in our analysis and discussion. Therefore, further values refer to the mean of open and closed systems and N₂O_emitted or N₂O_poreair at 5 cm and 12.5 cm unless explicitly stated otherwise. Gross rN₂O fractions tended to be higher in N₂O_emitted (treatment means 0.14 to 0.19) relative to the subsurface (treatment means 0.06 to 0.15). While water management treatment had a significant effect on process contributions to N₂O_emitted and N₂O_poreair (Table 5), significant interactions with depth and date were observed. Gross rN₂O fractions in N₂O_poreair were significantly lower in the DS-AWD relative to the WS-FLD on 6 of 15 days, with the WS-AWD falling in between. In the N₂O_emitted, the opposite pattern was mostly observed with gross rN₂O fractions often being higher in the DS-AWD than one or the other WS treatments, significantly so on 4 of 15 days. Aggregated across depths, the contribution of denitrification/nitrifier denitrification to N₂O_poreair were higher in the DS-AWD relative to one or both WS treatments on four dates and lower on three dates (Fig. 5a). The mean contribution of denitrification/nitrifier denitrification to N₂O_emitted ranged from 43 % to 49 % in all treatments (Fig. 6). Denitrification/nitrifier-denitrification contributions to N₂O_emitted were higher in the DS-AWD relative to the WS treatments on 9 and 23 June, and relative to WS-AWD only, they were also higher on 28 June and lower on 21 June.

Table 5 P-value ANCOVA results of modeled residual N₂O not reduced (gross rN₂O), fraction of total N₂+N₂O production coming from denitrification (gross frac_DEN), and the fraction of N₂O attributed to denitrification (DenContribution) derived from N₂O_emitted and N₂O_poreair. The Y position was used a covariate and represents the longitudinal position of each replicate within the field. Numbers in bold indicate significance at p < 0.05.

Figure 6 Estimated contribution of denitrification/nitrifier denitrification and nitrification/fungal denitrification to N₂O surface emissions in the three water management treatments (WS-FLD: water-seeding + conventional flooding; WS-AWD: water-seeding + alternate wetting and drying; DS-AWD: direct dry seeding + alternate wetting and drying). Estimates were derived from the mean of open and closed dynamics in a two end-member mixing model using δ¹⁸O-N₂O and SP values. Red and orange dashed vertical lines represent the date of seeding and fertilization in each treatment, respectively.

4 Discussion

4.1 Patterns of N₂O_emitted, N₂O_poreair and N₂O isotope ratios

In accordance with results from past studies (Miniotti et al., 2016; Peyron et al., 2016; Cai et al., 1997) and in line with our hypothesis, we observed higher N₂O emissions on most days in the DS-AWD relative to the two WS treatments (Fig. 2). A belated divergence in water management between the WS-FLD and WS-AWD (Table 1), in addition to a relatively wet early summer, likely contributed to similar observed soil environmental conditions and N substrates among these two treatments. Therefore, given the similarities in soil conditions, it is not surprising that N₂O fluxes and isotope ratio differences between these two treatments were generally fewer than expected. The lower yield in the DS-AWD treatment likely contributed additional differences in pore water N concentrations because lower N demand in this treatment should have resulted in higher soil N concentrations.

Mean daily δ¹⁵N, δ¹⁸O and SP values of N₂O_emitted and N₂O_poreair per depth and treatment ranged from −27.9 ‰ to 12.3 ‰, 30.9 ‰ to 63.0 ‰ and −14.0 ‰ to 53.2 ‰, respectively (Fig. 3). These values are similar in magnitude to those observed by Yano et al. (2014) in the early growing season of rice, where ranges of −24 ‰ to 6 ‰, 24 ‰ to 66 ‰ and 4 ‰ to 25 ‰ were reported. Our values are also similar in magnitude to those observed in other field studies which have included depth sampling (Koehler et al., 2012; Zou et al., 2014). Relative to these two studies we observed higher δ¹⁵N-N₂O and both higher and lower SP ratios. This was likely due to a higher sampling frequency, which covered more variable soil environments and generally higher soil moisture in our study than in the others. For example, it has been shown that organic matter decomposition and DOC availability in rice systems can decline with the introduction of wet–dry cycles or dry seeding (Said-Pullicino et al., 2016; Yao et al., 2011); thus, it is likely that conditions promoting complete denitrification declined in the AWD treatments. In contrast, saturated conditions favoring complete denitrification certainly prevailed in the WS treatments at times. Working in a denitrifying aquifer, Well et al. (2012) observed very large ranges in δ¹⁵N and SP ratios, varying from −55.4 ‰ to 89.4 ‰ and 1.8 ‰ to 97.9 ‰, respectively.

4.2 Source-partitioning N₂O production

One method to source partition emissions is to calculate net isotope effects and compare these to literature values derived from controlled and pure culture experiments where isotope effects were determined for individual processes. The calculated ${\mathrm{\Delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}/{\mathrm{NO}}_{\mathrm{3}}}$ in the DS-AWD treatment, with depth means of −7.2 ‰ to −16.0 ‰, was consistently much higher (i.e., less strong fractionation) than literature values reported for denitrification of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ (mean: $-\mathrm{42.9}\pm \mathrm{6.3}$ ‰; Denk et al., 2017) (Fig. S9). At 5 cm in the two WS treatments, the mean ${\mathrm{\Delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}/{\mathrm{NO}}_{\mathrm{3}}}$ was lower than in the DS-AWD (−23.2 and −21.5 in the WS-AWD and WS-FLD, respectively) but still nearly 20 ‰ higher than literature values. In a rice system, Yano et al. (2014) observed an ${\mathrm{\Delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}/{\mathrm{NO}}_{\mathrm{3}}}$ of −6.7 ‰, very well within the range of our calculated ${\mathrm{\Delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}/{\mathrm{NO}}_{\mathrm{3}}}$. Similarly, the global mean of our ${\mathrm{\Delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}/{\mathrm{NH}}_{\mathrm{4}}}$ values was −14.8 ‰, thus on average much higher than those reported in the literature for nitrification (−46.9 ‰, Sutka et al., 2006; $-\mathrm{56.6}\pm \mathrm{7.3}$ ‰, Denk et al., 2017). For both isotope effects, similar scenarios may explain our high observed Δ¹⁵N_x (i.e., low fractionation), namely, (i) non-steady-state reactions, for example rapid refreshing of the ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pools or near-complete substrate consumption, or (ii) significant reduction of N₂O serving to increase δ¹⁵N-N₂O values and thereby reducing the net isotope effect.

Considering the moist conditions and high reduction rates, it seems most likely that strong N₂O reduction was the largest contributor to the greater degree of isotopic discrimination observed. To check this, we estimated initial δ¹⁵N-N₂O values before N₂O reduction using our modeled N₂O reduction fraction (rN₂O), measured δ¹⁵N-N₂O values and an ¹⁵N isotope effect during reduction of −6.6 ‰ (Denk et al., 2017) in the Rayleigh equation. We could then estimate amended ${\mathrm{\Delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}/{\mathrm{NO}}_{\mathrm{3}}}$ values if N₂O reduction effects were accounted for from the difference between our initial δ¹⁵N-N₂O estimates and δ¹⁵N-${\mathrm{NO}}_{\mathrm{3}}^{-}$. These calculations yielded a ${\mathrm{\Delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}/{\mathrm{NO}}_{\mathrm{3}}}$ from −25.0 ‰ to −36.5 ‰, −32.6 ‰ to −42.3 ‰ and −29.0 ‰ to −51.1 ‰ in the DS-AWD, WS-AWD and WS-FLD across depths (Table S6). These amended ${\mathrm{\Delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}/{\mathrm{NO}}_{\mathrm{3}}}$ values do decrease and, especially for the WS treatments, come relatively close to literature values for ${\mathrm{\Delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}/{\mathrm{NO}}_{\mathrm{3}}}$ during denitrification. Thus, significant N₂O reduction can likely explain much of the high ${\mathrm{\Delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}/{\mathrm{NO}}_{\mathrm{3}}}$ values observed, particularly in the WS treatments. Yet other factors were also likely at play to some degree. For example, in the DS-AWD, where we observed evidence of significant nitrification, it is quite possible to envision isolated enrichment of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ at anaerobic micro-sites where N₂O is produced, while the bulk soil ${\mathrm{NO}}_{\mathrm{3}}^{-}$ pool remained less enriched. It is also true that we could not always measure δ¹⁵N values of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ or ${\mathrm{NH}}_{\mathrm{4}}^{+}$ because the concentrations were too low; thus, we could not calculate isotope effects. This highlights a persistent dilemma, which is true for all isotope ratios, that we cannot accurately measure isotope ratios at very low concentrations. Hence, until more sensitive methodologies are developed, in situ measurements such as these will always be biased toward higher-concentration scenarios where perhaps the strongest and most interesting effects of substrate enrichment are missed.

The use of any one isotope signature alone is confounded by overlap in the isotope effects between processes, unknown and possibly rapidly changing substrate δ values, and the unknown contribution of N₂O reduction effects. To overcome these drawbacks, graphical interpretations of dual N₂O isotope ratios have been used in field studies to interpret datasets similar to ours (Well et al., 2012; Koehler et al., 2012). For a more quantitative assessment of source partitioning, mixing models using a dual isotope approach can be used (Yano et al., 2014; Toyoda et al., 2011; Koba et al., 2009; Lewicka-Szczebak et al., 2017; Zou et al., 2014). In the subsequent analysis, we employ both approaches using our sample values plotted in SP × δ¹⁸O and SP × δ¹⁵N space (Figs. 4 and S10–S12).

In both SP × δ¹⁸O and SP × δ¹⁵N plots, our sample values mostly fell between the mixing and reduction lines predicted by either isotope relationship (Fig. 4) and somewhat surprisingly showed stronger enrichment, indicative of greater N₂O reduction in the DS-AWD treatment relative to the WS treatments. In the DS-AWD and to a lesser extent in the WS-AWD treatment, high pore air N₂O concentrations were associated with denitrification or nitrifier denitrification, while midrange concentrations were associated with a higher degree of N₂O reduction and the lowest concentrations fell neatly in between. Similarly, in the WS-FLD treatment, denitrification or nitrifier denitrification associated samples almost exclusively coincided with high N₂O_poreair. Most likely the moderate N₂O_poreair concentrations derived from N₂O reduction following high denitrification/nitrifier-denitrification production. This analysis is supported by data showing a trend of enrichment over the course of the measurement period (Fig. S10) and high WFPS values associated with the most enriched N₂O_poreair in the DS-AWD (Fig. S12). All treatments showed an enrichment of SP with time (Fig. S10), but interestingly only in the DS-AWD did δ¹⁸O and δ¹⁵N-N₂O enrich over the course of the experiment. This may reflect an increase over time in δ¹⁵N and δ¹⁸O of ${\mathrm{NO}}_{\mathrm{3}}^{-}$, which was observed in the DS-AWD treatment, albeit not strongly (Fig. S9). More ${\mathrm{NO}}_{\mathrm{3}}^{-}$ was available for denitrification in the DS-AWD treatment; thus, for greater enrichment of this pool to occur, we propose that more ${\mathrm{NO}}_{\mathrm{3}}^{-}$ was trapped in denitrifying micro-sites as the soil dried or O₂ was consumed.

In the WS treatments we observed a minimized trend of N₂O reduction compared to the DS-AWD treatment, more scattered high SP values and more values intermediate to the two end-member pools. These results may partially be explained by greater contributions from abiotic hydroxylamine decomposition (SP ∼34 ‰–35 ‰; Heil et al., 2014) or fungal denitrification (SP ∼35 ‰; Rohe et al., 2014). Zhou et al. (2001) showed that fungal denitrification requires minimal oxygen to proceed; similarly Seo and DeLaune (2010) found that fungal denitrification dominated relative to bacterial denitrification at modest reducing conditions to weakly oxidizing conditions (Eh > 250 mV). Indeed, there is some evidence that high scattered SP values corresponded to more moderate WFPS (70 %–90 %) in the WS-FLD treatment (Fig. S12). Abiotic hydroxylamine decomposition requires nitrification for the production of NH₂OH and iron or manganese (hyrdr)oxides as electron acceptors to proceed (Bremner et al., 1980). Given the moist conditions, nitrification rates were likely low in the WS treatments. Feasible co-occurrence of these species could really only occur directly in the rhizosphere of a flooded rice soil, where O₂ is transported to the immediate root zone by the aerenchyma. Tightly coupled nitrification–denitrification in the rhizosphere of rice plants has been shown before (Arth and Frenzel, 2000) as has coupling of nitrogen–iron transformations (Ratering and Schnell, 2000), but we cannot say the extent to which this may have occurred in our system.

It is necessary to contextualize N₂O isotope data with our measured substrate concentrations and soil environmental data. Based on our observations of low ${\mathrm{NH}}_{\mathrm{4}}^{+}$ concentrations, high ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentrations, an Eh over 400 mV and WFPS often below 60 % (5 cm) or below 85 % (12.5 and 25 cm) in the DS-AWD treatment, we can safely deduce that extensive nitrification of either basal urea fertilizer or of indigenous soil N occurred in this treatment (Fig. 2). Furthermore, the δ¹⁸O-${\mathrm{NO}}_{\mathrm{3}}^{-}$ in the DS-AWD treatment ranged from 0.1 to 14.8 (Fig. 7), thus falling in the range attributed to ${\mathrm{NO}}_{\mathrm{3}}^{-}$ produced from nitrification (Kendall and McDonnell, 1998). Additionally, we observed that both δ¹⁵N-${\mathrm{NO}}_{\mathrm{3}}^{-}$ and δ¹⁵N-${\mathrm{NH}}_{\mathrm{4}}^{+}$ were negatively correlated to substrate concentrations in the DS-AWD treatment, indicative of active consumption of both N substrates (Table 4). In the DS-AWD, there also was a positive correlation between δ¹⁵N-${\mathrm{NO}}_{\mathrm{3}}^{-}$ and N₂O_poreair but a negative correlation between δ¹⁵N-${\mathrm{NH}}_{\mathrm{4}}^{+}$ and N₂O_poreair. The former likely indicates N₂O production via denitrification and subsequent enrichment of the ${\mathrm{NO}}_{\mathrm{3}}^{-}$ pool. The latter is more difficult to interpret, but we attributed this to higher emissions associated with fresh inputs of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ (from urea or mineralization) which should have a δ¹⁵N value around 0 ‰. Together these data show that coupled nitrification–denitrification was responsible for the majority of N₂O emissions. Similar results were also reported by Dong et al. (2012) for an AWD system. The separation of isotope ratios by date, N₂O concentration and WFPS suggests that ${\mathrm{NO}}_{\mathrm{3}}^{-}$ produced early in the growing season was progressively denitrified and reduced over the course of the sampling period. Similarly, N₂O produced early in the growing season may have been progressively reduced.

Figure 7 Relationship of δ¹⁸O-${\mathrm{NO}}_{\mathrm{3}}^{-}$ to δ¹⁵N-${\mathrm{NO}}_{\mathrm{3}}^{-}$ in pore water samples of the three water management treatments (WS-FLD: water-seeding + conventional flooding; WS-AWD: water-seeding + alternate wetting and drying; DS-AWD: direct dry seeding + alternate wetting and drying). After Kendall and McDonnell (2012). The black arrow represents the trajectory of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ reduction effects. The black asterisk signifies the δ¹⁸O value atmospheric O₂ (25.3 ‰) while the dashed black line indicates the range of δ¹⁸O in soil water. δ¹⁸O−H₂O was not directly measured in our study. We assumed a value of −8.3 ‰ taken from an uncontained aquifer in the region by Rapti-Caputo and Martinelli (2009). The symbol colors indicate the concentration of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ in each sample (mg L⁻¹).

4.3 Inferring the extent of N₂O reduction

It has been suggested that the slope of SP/δ¹⁸O, SP/ δ¹⁵N and δ¹⁸O δ∕¹⁵N or their isotope effects can be used to estimate the extent of N₂O reduction (Jinuntuya-Nortman et al., 2008; Well and Flessa, 2009; Lewicka-Szczebak et al., 2017; Ostrom et al., 2007). However, many studies deriving these relationships have taken place under controlled conditions when N₂O supply was often limited. Therefore, fractionation following closed-system dynamics would result in larger fractionation effects on the residual substrate than under open-system dynamics. The positive and significant relationship between all isotopes and across all depths in the DS-AWD treatment suggests an influence of reduction at all depths. In contrast, in the WS treatments we observed no relationship between SP and δ¹⁸O within N₂O_emitted (Fig. S7) and only a weak relationship between SP and δ¹⁵N at 25 cm in the WS-AWD and even a negative relationship between SP and δ¹⁵N in the WS-FLD N₂O_emitted (Fig. S8). The range of observed δ¹⁸O/ δ¹⁵N slopes (0.21 to 0.90; Fig. S5) was substantially lower than those observed in many N₂O reduction studies (1.94 to 2.6; Jinuntuya-Nortman et al., 2008; Ostrom et al., 2007; Well and Flessa, 2009; Lewicka-Szczebak et al., 2017) but closer to the 0.45 slope observed by Yano et al. (2014) in an in situ rice field study. When a significant relationship was observed, overall or N₂O_poreair SP∕δ¹⁵N slopes ranged from 0.49 to 0.83 (Fig. 4b). These slopes are either close to those of other field studies (0.48 to 0.52; Yano et al., 2014; Wolf et al., 2015) or intermediary between field studies and controlled N₂O reduction studies (0.59 to 1.01; Well and Flessa, 2009; Lewicka-Szczebak et al., 2017). From controlled N₂O reduction studies, an SP∕δ¹⁸O slope between 0.2 to 0.4 has been observed (Jinuntuya-Nortman et al., 2008; Well and Flessa, 2009); thus, in this case the N₂O_poreair slopes observed in our study were substantially higher (Figs. 4a and S7). The lower overall SP and δ¹⁸O slope in the WS treatments was due to the inclusion of the N₂O_emitted values, which individually showed no relationship in these treatments.

A deviation in slopes compared to those observed in controlled N₂O reduction studies likely points to a growing influence of open-system dynamics where substrates are continuously refreshed. It has been demonstrated that when mixing processes dominate over reduction processes, the SP∕δ¹⁸O slope rises (Lewicka-Szczebak et al., 2017). It is also plausible that high rates of oxygen exchange during denitrification served to partially mask an increase in δ¹⁸O-N₂O values, resulting in the higher observed SP/δ¹⁸O slopes or lower δ¹⁸O∕δ¹⁵N slopes. To estimate the extent of oxygen exchange with denitrification precursors (NO_x) we plotted ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}-{\mathrm{N}}_{\mathrm{2}}\mathrm{O}/{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}-{\mathrm{NO}}_{\mathrm{3}}^{-}$ by ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}-{\mathrm{H}}_{\mathrm{2}}\mathrm{O}/{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$-${\mathrm{NO}}_{\mathrm{3}}^{-}$ following Snider et al. (2009). The slope of this relationship ranged from 0.7 to 2.1 (data not shown). Thus, we assume oxygen exchange was effectively 100 % across treatments during denitrification. In summary, the observed positive relationships between the isotope pairs is indicative of an influential role of N₂O reduction in the DS-AWD treatment. This is less clear in the WS treatments where relationships were more erratic, suggesting a stronger influence of changing nitrification and denitrification process rates or changing δ¹⁵N of N substrates. It is likely that isotope ratios in the WS treatments were affected by near-complete denitrification to N₂. Well et al. (2012) observed highly variable isotope ratios in a strongly denitrifying aquifer and concluded that N₂O reduction was strongly progressed but variable. However, it should be noted that their system had abundant ${\mathrm{NO}}_{\mathrm{3}}^{-}$ while ours did not. The inconsistent relationships between N₂O_emitted and N₂O_poreair for SP∕δ¹⁵N and SP∕δ¹⁸O in the WS treatments and the stronger enrichment observed in the DS-AWD N₂O_emitted (Fig. 4) demonstrate a disconnection between subsurface N₂O_emitted and N₂O_poreair across treatments. Such results suggest that N₂O reduction may not have had as strong of an influence on the signature of N₂O_emitted as it did on N₂O_poreair, particularly in the WS treatments. A decoupling between subsurface N₂O concentrations and surface emissions and their isotope ratios has been observed in other studies (Van Groenigen et al., 2005; Goldberg et al., 2010a). This phenomenon is most simply explained by emitted N₂O truly coming from a mix of sources and depths, while subsurface N₂O is representative of a much smaller spatial zone and more likely to be dominated by one process. While difficult to practically measure, processes at shallow depths above 5 cm were also likely influential to surface emissions.

4.4 Complementary evidence from a two end-member mixing model approach

To quantitatively estimate the extent of N₂O reduction (gross rN₂O), N₂O production and reduction rates, and the contribution of denitrification to N₂O emissions, we used an open- and closed-system two end-member mixing model based on SP-N₂O and δ¹⁸O-N₂O relationships. As described in Sect. 2.7, we tested our models under two scenarios; in scenario 1 (sc1) N₂O is produced and reduced by denitrifiers before mixing with N₂O derived from nitrification; in scenario 2 (sc2) N₂O is produced from both processes, mixed and then reduced (Fig. 1). While we could estimate gross rN₂O and the fraction of denitrification from both scenarios, sc2 yielded mostly implausible solutions for the contribution of denitrification to N₂O in open systems (Fig. S3 and Table S2). We thus conclude that the assumptions underlying this scenario in open systems were not valid in our system. In a closed system N₂O is progressively consumed and not replenished, resulting in a stronger isotope effect and faster enrichment of the remaining N₂O; thus, a smaller degree of N₂O reduction is needed to achieve an equivalent enrichment as in open systems. Our results for open and closed systems align well with this theory on N₂O fractionation. However, we feel strongly that with in situ measurements in a heterogeneous soil environment, a combination of closed- and open-system dynamics likely exits; therefore, the following interpretation of our data is based on an average of open- and closed-system values. Given the lower moisture and evidence of extensive nitrification occurring in the DS-AWD treatment, we expected a higher contribution of nitrification/fungal denitrification in this treatment, coming from an increase in nitrification. However, this was not the case and denitrification/nitrifier-denitrification contributions tended to be higher in the DS-AWD treatment relative to WS treatments (Figs. 5a, 6). Treatment differences were significant in the surface fluxes; however, there was a significant interaction with sampling day. There was no treatment effect on denitrification contribution in the subsurface (Table 5). The equivalent or higher contributions of nitrification/fungal denitrification in the WS treatments (Fig. 6) are most easily explained by higher fungal denitrification; in their laboratory experiments, Lewicka-Szczebak et al. (2017) also observed relatively high fungal denitrification contributions under very wet conditions. Larger contributions from fungal denitrification would also help explain the less clear reduction trends as fungal denitrifiers are thought to largely produce N₂O as an end-product rather than N₂. It should be noted that due to low surface fluxes or N₂O_poreair, we had fewer data points in the WS treatments. Previous studies have attributed significant amounts of N₂O emissions in paddy systems to nitrification in periods of low soil moisture (Lagomarsino et al., 2016; Verhoeven et al., 2018). Yet such studies were not able to quantitatively source partition emissions. Given our results here, it is possible that N₂O produced either via nitrifier denitrification or coupled nitrification–denitrification has been previously underestimated.

The modeled gross rN₂O fractions indicate high levels of N₂O reduction for all treatments and depths, (rN₂O: 0.06 to 0.19) even in the DS-AWD where soil moisture was frequently below 60 % at 5 cm (Fig. 2). These results are at first surprising, but there is still much we do not know about subsurface N₂O production and consumption. Direct measurements of N₂O reduction at depth are few. Using membrane inlet mass spectrometry, Zhou et al. (2017) detected higher N₂O reduction to N₂ in paddy soil water at 20 cm versus 60 or 80 cm and could relate this to higher DOC concentrations at 20 cm. Other studies suggest high subsurface N₂O reduction based on the inference of declining N₂O concentration accompanied by isotope enrichment moving up a soil profile (Goldberg et al., 2008; Clough et al., 1998; Van Groenigen et al., 2005). We are also methodologically limited by our inability to measure N₂O isotopes at near or complete N₂O reduction because there is too little remaining N₂O to measure. We assume this was more often the case in the WS treatments; therefore, we postulate that the signature of N₂O reduction was stronger in the DS-AWD largely because there was more N₂O left to measure. In their experiments to validate the mixing model we used, Lewicka-Szczebak et al. (2017) found that the model routinely underestimated gross rN₂O rates relative to measured rates in an oxic mineral soil but performed better under anoxic conditions and in an organic soil. Therefore, an underestimation of rN₂O rates, particularly in the DS-AWD treatments, remains possible. However, considering the strong indication of N₂O reduction from other isotope relationships (i.e., SP and δ¹⁵N and δ¹⁵N and δ¹⁸O), we believe that subsurface N₂O reduction rates were simply high in our system, regardless of water management.

Table 6 Spearman correlations between modeled rN₂O-gross, frac_DEN-gross and DenContribution with soil environmental variables and inorganic N substrates and δ¹⁵N-N₂O. Results are for the mean of open- and closed-system dynamics. Subsurface correlations were performed on data aggregated across 5 and 12.5 cm depths. Significance indicated by ${}^{****}$ < 0.0001, ${}^{***}$ < 0.001, ${}^{**}$ < 0.01, ^* < 0.05.

In the subsurface, the contribution of denitrification/nitrifier denitrification to N₂O concentrations was positively correlated to N₂O_poreair concentrations and WFPS in all treatments, indicating an increasing contribution of denitrification/nitrifier denitrification at times of higher N₂O production in conjunction with rising soil moisture (Table 6). In the two AWD treatments, the contribution of denitrification/nitrifier denitrification negatively correlated to δ¹⁵N signature of N₂O_poreair and N₂O_emitted (DS-AWD only). Many studies have demonstrated that high subsurface N₂O production is correlated to depleted δ¹⁵N-N₂O (Goldberg et al., 2008, 2010b; Van Groenigen et al., 2005). These results further support the conclusion that high N₂O_poreair and N₂O_emitted were produced from denitrification/nitrifier denitrification associated with more depleted δ¹⁵N-N₂O. Higher gross rN₂O (less N₂O reduction) was associated with higher N₂O_emitted in all treatments and higher N₂O_poreair (WS-AWD only), demonstrating that higher N₂O resulted not only from increased denitrification/nitrifier denitrification but also from a decrease in N₂O reduction. Interestingly, higher rN₂O in N₂O_emitted of the DS-AWD was also associated with higher WFPS. Such a result can only be explained by a dependency of reduction on N₂O production. Overall, there was a negative relationship between rN₂O and δ¹⁵N-N₂O, yet the relationship was not consistently strong or significant between treatments. A negative relationship supports an isotope enrichment effect with greater N₂O reduction. Considering the above, it appears that maximum N₂O production and emissions occurred during periods of increased contribution from denitrification/nitrifier denitrification, which were accompanied by small declines in N₂O reduction. These relationships were most robust in the DS-AWD treatment. Correlations within the N₂O_emitted dataset were undoubtedly affected by lower data availability, particularly in the WS treatments, and should be taken with caution. Despite the high estimates of N₂O reduction for all treatments, we still observed relevant contributions from nitrification/fungal denitrification on many dates (Fig. 6). Nevertheless, the highest fluxes in the DS-AWD aligned with higher contributions from denitrification/nitrifier denitrification, while the highest fluxes in the WS treatment had nitrification/fungal denitrification contributions of ca. 50 %. In the WS treatments we again postulate that fungal denitrification rates increased because conditions were not ideal for high nitrification. Studies have shown that fungal denitrification and co-denitrification can play a significant role in soil N₂ and N₂O emissions from soil (Long et al., 2013; Laughlin and Stevens, 2002).

From our modeling results we could estimate N₂ production or emissions based on our calculated N₂O reduction rates (Fig. S13). Due to poor data availability and high variability, we could neither confidently estimate N₂ production at 25 cm nor surface N₂ emissions on many dates of the WS treatments, but we have more confidence in the estimates obtained for the DS-AWD treatment. Mean daily N₂ emissions found in our study were 236±53 (n=43), 194±37 (n=41) and 197±35 (n=31) g N ha⁻¹ d⁻¹ in the DS-AWD, WS-AWD and WS-FLD, respectively. To our knowledge only one other study by Yano et al. (2014) has conducted similar calculations to estimate N₂ emissions in rice systems from isotope ratios. The authors also found high rates of N₂O reduction, around 80 % to 85 %, corresponding to an rN₂O of 0.15 to 0.20 and N₂ emissions between 0.1 to 422 µg N m² h⁻¹ (or 0.024 to 101.4 g  ha⁻¹ d⁻¹). Therefore, the estimated extent of N₂O reduction was quite similar to our surface-emitted reduction rates, with somewhat lower N₂ emissions corresponding to somewhat lower N₂O emissions. Using labeled ¹⁵N urea, Lindau et al. (1990) measured N₂ emissions of 254 g ha⁻¹ d⁻¹, while Dong et al. (2012) observed similar rates of 194 g N₂-N ha⁻¹ d⁻¹ for an AWD treatment. Considering that these results only account for N₂ derived from fertilizer, the modeled mean daily N₂ emissions found in our study are plausible. Differences between the treatment means were not significant for N₂O_poreair or N₂O_emitted (p=0.431 and p=0.858) and thus do not indicate a higher potential for N₂ losses in the WS treatments. We must reject our hypothesis that higher ${\mathrm{NO}}_{\mathrm{3}}^{-}$ in the WS-AWD relative to the WS-FLD would drive higher denitrification and N₂ losses because we observed no differences in final modeled N₂ production and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentrations were essentially zero for both WS treatments. Our results show there is promise for estimating N₂ emissions from N₂O isotope ratios using simple models, but the precision of these estimates remains constrained by the limitations discussed below.

All modeling attempts to date rely on isotope signatures and effects determined in laboratory studies, and thus changes in these values in response to environmental or microbial population dynamics in the field remain a large question. As this was an in situ field experiment, conditions were not constant across treatments or throughout the sampling time frame, yet it has been shown that isotope effects, particularly for N₂O reduction, change with shifts in environmental conditions such as increasing water-filled pore space (Jinuntuya-Northman et al., 2008). Therefore, the use of fixed isotope effects in our model is a simplification. Future modeling efforts may be improved by the incorporation of variable isotope effects based on soil moisture or O₂ for example. Careful, controlled experiments across a range of soils with different management histories are necessary to determine if consistent variation in isotope effects in relation to specific environmental parameters can be determined or if such parameters are site specific. The microbial δ¹⁸O signature for denitrification used in our model was calculated relative to δ¹⁸O-H₂O. We therefore assumed complete exchange between N₂O substrates, intermediaries and water during denitrification. We based this on previous work showing that O exchange is high and that the isotope effect between water and N₂O is relatively stable (Lewicka-Szczebak et al., 2016, 2017; Snider et al., 2013; Kool et al., 2007). In reality, results over time and between treatments may have been affected by varying degrees of ¹⁸O exchange between N₂O, intermediaries and water and by variation in δ¹⁸O−H₂O values. We recommend that future studies measure the δ¹⁸O-H₂O to better constrain results. Modeling results would also be more robust if complete δ¹⁵N -N₂O, -${\mathrm{NH}}_{\mathrm{4}}^{+}$ and -${\mathrm{NO}}_{\mathrm{3}}^{-}$ across treatments and times were available, allowing for complementary modeling of SP × ¹⁵N(${\mathrm{N}}_{\mathrm{2}}\mathrm{O}/{\mathrm{NO}}_{\mathrm{3}}^{-}$ or N₂O∕NH₄). As an exploration of model sensitivity to the isotope signature used, we evaluated results across a range of ${\mathit{\delta }}_{\mathrm{0}}^{\mathrm{18}}\mathrm{O}-{\mathrm{N}}_{\mathrm{2}}\mathrm{O}$_nit and ${\mathit{\delta }}_{\mathrm{0}}^{\mathrm{18}}\mathrm{O}-{\mathrm{N}}_{\mathrm{2}}\mathrm{O}$_den values (Fig. S14) before selecting those used in Lewicka-Szczebak et al. (2017). More elaborate modeling efforts could employ iterative simulation techniques where a range of literature values for N₂O signatures and isotope effects are used and drawn to help highlight model sensitivity to specific isotope values and improve its accuracy. Lastly, more work needs to be done to validate results such as those generated here which rely on laboratory derived values, with complementary measurements of microbial community dynamics, such as that by Snider et al. (2015).

5 Conclusions

The relatively dry conditions in the DS-AWD treatment and application of urea fertilizer led to extensive nitrification, subsequent denitrification and denitrification-derived N₂O emissions. Even with evidence of nitrification and relatively aerobic conditions in the DS-AWD treatment, both graphical and two end-member mixing model results indicated significant N₂O reduction in all treatments, graphically most convincingly in the DS-AWD treatment. Treatment differences may also reflect paddy history as this was the fifth year of alternative water management at the site. Yields were also lower in the DS-AWD, which likely lowered N demand and increased soil N concentrations in this treatment. Differences between depths were often more evident in N₂O_poreair, ${\mathrm{NO}}_{\mathrm{3}}^{-}$, ${\mathrm{NH}}_{\mathrm{4}}^{+}$ and DOC concentrations than in N₂O isotope signatures at the various depths, particularly for the WS treatments. In the DS-AWD treatment, isotope signatures of δ¹⁸O-N₂O and SP values demonstrated notably lower values at 5 cm relative to other depths, mostly likely indicating higher N₂O production and less reduction in the upper layer. Overall, the highest N₂O production and emissions were associated with an increasing contribution from denitrification/nitrifier denitrification accompanied by decreases in N₂O reduction in the AWD treatments. Our isotope data suggest that contributions from fungal denitrification to N₂O emissions may have increased in the WS-FLD treatment. The role of fungal denitrification in paddy rice systems should be further investigated with the use of fungal inhibitors. Surface-emitted N₂O reduction rates were similar for all treatments; therefore, our hypothesis of a greater potential for gaseous N₂ losses in the WS-AWD is refuted. Despite the difficulty in obtaining a full dataset for all treatments and the inherent spatiotemporal variability in the original measured fluxes, we came to good agreement with the magnitude of N₂ emissions reported from previous ¹⁵N labeled fertilizer studies. Thus, natural abundance isotope methods do show promise for estimating N₂ emissions and closing N budgets, even without the δ¹⁵N of N substrates. Model results would likely improve with controlled incubations to determine site-specific isotope effects and whether these effects change in a consistent manner with specific environmental conditions. In saturated or partly saturated systems, future studies should probe the disconnection between subsurface and emitted N₂O isotopes by employing methods that allow for larger subsurface spatial integration along vertical and horizontal planes. It appears that to effectively manage N losses in alternative water management paddy systems, the inhibition of nitrification is necessary, particularly very early in the growing season when N availability exceeds crop N demand.