Proper Interpretation of Dissolved Nitrous Oxide Isotopes, Production Pathways, and Emissions Requires a Modelling Approach

Stable isotopes ( 15N and 18O) of the greenhouse gas N2O provide information about the sources and processes leading to N2O production and emission from aquatic ecosystems to the atmosphere. In turn, this describes the fate of nitrogen in the aquatic environment since N2O is an obligate intermediate of denitrification and can be a by-product of nitrification. However, due to exchange with the atmosphere, the values at typical concentrations in aquatic ecosystems differ significantly from both the source of N2O and the N2O emitted to the atmosphere. A dynamic model, SIDNO, was developed to explore the relationship between the isotopic ratios of N2O, N2O source, and the emitted N2O. If the N2O production rate or isotopic ratios vary, then the N2O concentration and isotopic ratios may vary or be constant, not necessarily concomitantly, depending on the synchronicity of production rate and source isotopic ratios. Thus prima facie interpretation of patterns in dissolved N2O concentrations and isotopic ratios is difficult. The dynamic model may be used to correctly interpret diel field data and allows for the estimation of the gas exchange coefficient, N2O production rate, and the production-weighted values of the N2O source in aquatic ecosystems. Combining field data with these modelling efforts allows this critical piece of nitrogen cycling and N2O flux to the atmosphere to be assessed.


Introduction
Nitrous oxide (N 2 O) is a powerful greenhouse gas, 298 times more potent than CO 2 over a 100-year time line [1]. Atmospheric N 2 O concentrations have been increasing at a rate of 0.25%/year over the last 150 years [2]. Consequently, the global N 2 O budget has been the subject of intensive research efforts over the past few decades. N 2 O is produced through multiple microbial pathways: hydroxylamine oxidation during nitrification and as an obligate intermediate during denitrification and nitrifier-denitrification. Because these pathways of N 2 O production have different stable isotopic enrichment factors, isotopic analysis of N 2 O can potentially distinguish N 2 O produced through different pathways or from different sources [3]. Identifying N 2 O sources will provide insights on the fate of N at the ecosystem-scale (e.g., [4][5][6]). The isotopic ratios of N 2 O produced in soil environments (e.g., [7][8][9][10][11]), and in aquatic environments (e.g., [12][13][14][15][16][17][18]) have been measured to some extent. Although N 2 O production in rivers and estuaries is a significant portion of the global N 2 O budget (approximately 1. 5 TgN/year, [19]), few studies report isotopic data for rivers [5,20,21].
In ice-free aquatic ecosystems, the d 15 N and d 18 O of dissolved N 2 O is affected by gas exchange with the atmosphere. As a result, the isotopic ratios of dissolved N 2 O are not equal to those of the N 2 O produced within the aquatic ecosystem and continue to change as atmospheric exchange (both ingassing and outgassing) occurs. In addition, isotopic fractionation during influx and efflux causes the isotopic ratios of N 2 O flux emitted to the atmosphere to be different than that of the dissolved N 2 O [22]. Thus, the simple method of calculating the instantaneous isotopic ratios of the N 2 O flux by taking measured dissolved isotopic ratios, adding an equilibrium isotope fractionation, and applying them to measured flux rates is inappropriate. Adjustments of measured isotopic ratios are necessary to understand the isotopic ratios of both produced and emitted N 2 O.
In this paper, we present a dynamic model of the stable isotopic composition of both the dissolved and emitted N 2 O in aqueous systems. We apply this model to two different measured diel patterns of the isotopic ratios of N 2 O in an aquatic ecosystem. We use the model to elucidate the relationship between the isotopic ratios of source, dissolved, and emitted N 2 O, to allow for improved interpretation of dissolved N 2 O isotope data. Ultimately, a process-based understanding on N cycling with aquatic ecosystems may be developed based on interpretation of N cycling processes.  15 N isotopomers can be measured (e.g., [23][24][25]), the gas exchange fractionation factors are not affected by the intramolecular distribution of 15 N [22]. Many laboratories cannot measure the intramolecular distribution of 15 N and analysis of the bulk 15 N: 14 N ratio of N 2 O is more common [26]. Here, we confine our analysis to bulk 15 N: 14 N ratios and use 15 N 2 O to represent the average abundance of the two 15 N isotopomers. The same approach could easily be extended to consider each isotopologue separately.

Dynamic Isotope Model for Dissolved N 2 O
A simple three box model (SIDNO, Stable Isotopes of Dissolved Nitrous Oxide) was created using Stella modelling software (version 9.1.4, http://www.iseesystems.com) in order to study the relationships between the isotopic ratios of source, dissolved and emitted N 2 O (model file is available at https://github.com/ jjvenky/SIDNO and by contacting the corresponding author). This model is an adaptation of the isotopic gas exchange portion of the PoRGy model [27], which successfully modelled diel isotopic ratios of O 2 resulting from photosynthesis, respiration, and gas exchange in aquatic ecosystems. One key difference is photosynthetically produced O 2 in PoRGy has a d 18 [28,29] though certain waters can exhibit significant N 2 O reduction to N 2 [30,31]. The masses and magnitude of the flows of 15 [22]. These two a values are related to the equilibrium fractionation factor: a eq~Rgas =R dissolved~aev =a in (0.99925 for d 15 N and 0.99894 for d 18 O) and are independent of temperature over the range of 0 0 C to 44.5 0 C [22].
The d values of tropospheric N 2 O are 6.72 %+ 0.12% for d 15 N and 44.62% + 0.21% for d 18 O [38]. Therefore, at equilibrium, dissolved N 2 O has dissolved d values slightly greater than these at 7.48% and 45.73%, respectively.
In the model, net N 2 O flux between the atmosphere and dissolved phase was calculated using the thin boundary layer approach as: where the N 2 O flux is calculated in mol/m 2 /h, k is the usermodifiable gas exchange coefficient (m/h), is the partial pressure of tropospheric N 2 O (assumed to be 320 ppbv from data provided by the ALE GAGE AGAGE investigators, [39,40] where T is temperature in kelvins. Gas exchange is a two-way process. The net N 2 O flux rate (the difference between the invasion and evasion rates) depends on the dissolved N 2 O concentration. When a solution is at equilibrium with the atmosphere, the invasion and evasion rates will be equal, and the net flux will be zero. As

Test of Model Performance
To test the ability of SIDNO to reproduce observed isotopic data, input parameters (N 2 O production rate, N 2 O d values, and k) were set to replicate a series of experiments designed to derive fractionation factors for N 2 O gas exchange [22]. In these experiments, degassed water was exposed to N 2 O gas of known isotopic ratios in a sealed container to varying degrees of saturation.   [22]. The coefficient of determination for experimental data and SIDNO model outputs were comparable to those of [22], R 2~0 :77 for d 15 [38]. The trajectories on these plots ( Figure 2C

Modelling Scenarios with Steady State Production of N 2 O
The SIDNO model can be used to probe the stable isotope dynamics of N 2 O in a variety of situations that may be encountered in aquatic environments to elucidate the relationship between the N 2 O source (a function of N cycling processes), dissolved (the easily measured component), and emitted (of consequence for greenhouse gas production and global N and In the steady-state production of N 2 O (constant rate and d values), by definition, the d values of N 2 O production must be the same as those of the emitted N 2 O. As a result, the d values of the dissolved N 2 O cannot equal that of the source (or emitted) N 2 O at steady state because the dissolved N 2 O must be offset from the emitted N 2 O by at least the a ev values. As the steady-state production rate was increased, the steady-state N 2 O concentration increased and the dissolved d values approached but did not equal the source (Figure 3). Even at moderate supersaturations (v1000%) the effect of atmospheric N 2 O equilibration on the d values of dissolved N 2 O cannot be ignored.
At steady state, the d values of the emitted N 2 O must be equal to the source; the large difference between source/emitted and dissolved N 2 O underscores the importance of adjusting the measured d values of dissolved N 2 O in order to determine aquatic contributions of N 2 O to the atmosphere or N 2 O sources. This is critical when using dissolved measurements of N 2 O to constrain the global isotopic N 2 O budget, but not been done in most studies, e.g., [16,[42][43][44] but see [45].

Modelling Scenarios with Variable Production of N 2 O
The relationship between the d values of source, dissolved, and emitted N 2 O are much more complicated when N 2 O production is variable rather than when it is constant. N 2 O production may vary with respect to production rate and/or d values; in many aquatic environments, N 2 O production is not likely to be constant. The N 2 O production processes, nitrification and denitrification, are sensitive to redox conditions, which can be highly variable, due to diel changes in dissolved O 2 concentration, flow regime, etc. For example, [34] observed diel changes in the denitrification rate in the Iroquois River and Sugar Creek (Midwestern USA) and found that the denitrification was consistently greater during the day than night. The relative importance of nitrification and denitrification can change in response to the diel oxygen cycle: e.g., [46] observed a change from daytime nitrification to nighttime denitrification in a subtropical eutrophic stream. Coupling of N 2 O and O 2 diel cycles has been observed in agricultural and waste-water treatment plant (WWTP) impacted rivers [36]. Since fractionation factors and substrates are different for nitrification and denitrification, ecosystem-scale fractionation factors may be rate and process dependent, and the d values of N 2 O production in a given ecosystem may not be constant over a diel cycle.
To simulate the diel variability, various scenarios were modelled by adjusting either production rate and/or the associated d values. The variabilities in these input parameters were driven by a sine function with a 24 h period similar to a dissolve O 2 curve. In all scenarios, the chosen range of production rates was based on published N 2 O flux rates (Table 1) and varied from 1 to  (Table 2), which was between the diel variation in N 2 O flux observed by [33] and [46]. Temperature was held constant at 20 0 C. The value of k was varied between 0.1 and 0.3 m/d (Table 2), within the range observed in other river studies ( Table 1). The combination of production rates and k values were chosen to produce N 2 O between 150% and 500% saturation (

Model Scenario #2: Constant Production Rate, Variable Isotopic Composition of Source
In scenario #2, when the N 2 O production rate was held constant and the d values of the source varied with time (from 250% to 210% for d 15 N and from 10% to 30% for d 18 O), the d values of the dissolved N 2 O was also much farther from that of the source than the dissolved N 2 O due to the effects of atmospheric exchange and the emitted N 2 O varies linearly between the two source values. In contrast, the dissolved N 2 O is parallel but offset from the line connecting the two sources ( Figure 5, Table 3

Model Scenario #3: Constant Production Rate, Variable Isotopic Composition of Source
To examine the effects of varying k on the scenario of constant N 2 O production with variable isotopic signature of the source, k was reduced from 0.3 m/h (scenario #2) to 0.1 m/h (scenario #3; Figure 6, Table 3  To simulate a system alternating between two N 2 O production processes, such as differing relative contributions of nitrification and denitrification, with different rates of N 2 O production and d values, the model was run with both production rate and its d values variable with time (scenarios #4, #5, and #6). The production rate and d values were adjusted so that the maximum rate coincided with the lowest source d values in scenarios #4 and #6 and so that maximum rate coincided with the highest source d values in scenario #5.

Model Scenario #4: Variable Production Rate, Variable Isotopic Composition of Source
For scenario #4, the resulting N 2 O concentrations were identical to those in model scenario #1, with the maximum concentration lagging approximately 2.75 h behind the maximum production rate (Figure 7, Table 3). The relationship between the d values of the dissolved and emitted N 2 O was more complex than in other scenarios. The lag time between the maximum source d values and those of dissolved and emitted N 2 O (when the production rate was minimum) was 3.75 h; however, the lag time between the minimum source d values and those of the dissolved and emitted N 2 O (when the production rate was maximum) was only 2.25 h. The difference between the emitted and source N 2 O was 3.1% to 8.0% for d 15 N and 1.3% to 3.4% for d 18 O. The d values of emitted N 2 O were closer to those of the source during periods of high production rates (and thus higher concentrations) than periods of low production rates. However, the flux-weighted average d values of emitted N 2 O were equal to the average production-weighted source d values because the system was at dynamic steady state.

Model Scenario #5: Variable Production Rate, Variable Isotopic Composition of Source
The isotopic counterpoint to scenario #4 is adjusting the timing of maximum N 2 O production to coincide with the highest d values of production (scenario #5). All other parameters were the same as scenario #4 (Table 3). The resulting pattern for the d values of dissolved N 2 O was very different than scenario #4 (Figure 8, Table 3

Grand River
The ability of SIDNO to reproduce measured patterns of N 2 O concentration and d values in a human-impacted river was also assessed. The Grand River is a seventh-order, 300 km long river that drains 6800 km 2 in southern Ontario, Canada, into Lake Erie, see [36,37,51]. There are 30 WWTPs in the catchment and their cumulative impact can be observed via the increase in artificial sweeteners in the river [52].
Samples were collected approximately hourly for 28 h at two sites in the central, urbanized portion of the river: sites 9 and 11 in [51,52]. The upstream site, Bridgeport, is where the river enters the urban section of the river at the city of Waterloo and is immediately above that city's WWTP. Blair is 26.6 km downstream of Bridgeport and below the cities of Waterloo and Kitchener. It is also 5.5 km downstream of the Kitchener WWTP. Average river depth at both sites was 30 cm. Values of k were determined by best-fit modelling of diel O 2 and d 18 O-O 2 values at the sites [36]. N 2 O concentration analyses were performed on a Varian CP-3800 gas chromatograph with an electron capture detector and isotopic ratio analyses were performed on a GV TraceGas pre-concentrator coupled to a GV Isoprime isotope ratio mass spectrometer, see [5] for analytical details.
Data from upstream and downstream of large urban wastewater treatment plants on the Grand River show diel patterns in N 2 O saturation and d values (Figures 10 and 11). At the Bridgeport site, the diel patterns of N 2 O saturation and d 15 N values were opposite of each other, that is, when N 2 O saturation was highest around sunrise the d 15 N values were lowest and when For both Bridgeport and Blair data, the cause of poorer fits for d 18 [45,46], may also manifest itself in d 18 O-N 2 O values because of the strong O isotope fractionation factor during denitrification [55].  Figures 4 with 6). Qualitatively, the d values of emitted N 2 O will be similar to the source if the equilibration time of dissolved N 2 O is small relative to the period of source variability (e.g., 24 h period due to diel changes in N cycling [36,45]). Assuming homogeneous N 2 O release upstream, the equilibration time can be approximated from a decay curve as 3z k , where z is mean depth [49]. If z is small and/or k is high, the equilibration time will be short and the d values of the emitted N 2 O will be close to the source. With decreasing k (or increasing equilibration time), the d values of emitted N 2 O will lag farther behind and will always have a smaller range of d values than the source. At the most extreme case, the variability in the d values of emitted N 2 O will be reduced to nearly    [53]). Nevertheless, these values can be measured and biogeochemical relationships between N species, redox, and N 2 O can be used as supporting information for process separation (e.g., [5,12,36,45]  Ultimately, the dynamic model SIDNO may be used to estimate k, N 2 O production rate and d values of the N 2 O source, an indication of the production pathway and N cycling, in aquatic ecosystems via inverse modelling. If physical properties, such as depth and temperature are known, SIDNO may be used to fit the measured field data (N 2 O concentration and d values) by adjusting the N 2 O source parameters. SIDNO can also be used to explore the dynamics between dissolved, source, and emitted N 2 O to query production scenarios and design field campaigns for studies of N cycling processes.