Introduction

The exponential growth of the human population concomitant with intensive development of fertilized agricultural and industrial activity since the 1950s caused sharp increases in nitrogen loadings to rivers and surficial aquifers worldwide1,2. Dissolved reactive nitrogenous species (e.g., NO3−, NH4+) are prevailing pollutants in many rivers and aquifers, stemming primarily from agricultural activities, municipal waste sources, and combustion derived nitrogen (N) deposition3. The impact of reactive nitrogenous species on water and ecosystems (e.g., eutrophication, hypertrophication) and human health (e.g., methemoglobinemia, cancer, thyroid disease)4,5,6 occurs in series, known as N-cascading1, and is of existential concern.

The global nitrogen cycle is a subject of considerable debate and concern that anthropogenic loadings are driving it beyond Earth’s natural resilience boundaries7. Reactive nitrogen cascades through aquatic ecosystems differentially1, since some systems accumulate N, whereas others transform it through diverse biogeochemical N-cycling processes like nitrification, denitrification, N2-fixation, dissimilatory nitrate reduction (DNRA), ammonification, and biological assimilation at rates dependent on the environmental conditions. These complexities make it difficult to unravel what point-based NO3− concentrations in rivers or their connected aquifers embody at any point in time, apart from regulatory pollutant exceedances, and quantitative knowledge of the roles of N-cycling processes in the terrestrial aquatic environments remains deficient.

Stable isotopes of NO3− (δ15N, δ18O) are used, particularly in the last decades, to help identify the sources of high N-pollution in aquatic systems8,9,10,11,12. Contemporary preparative and isotope techniques enable fast low-cost isotopic analysis of 15N/14N and 18O/16O ratios to ppb concentrations to incorporate comparative isotopic information from pristine aquatic environments[13, and references therein]. Distinctive δ15N and δ18O biplot clusters were initially proposed to “assign” organic and inorganic NO3− origins3,14,15, however, NO3− sources rarely exhibit unique combinations of nitrogen (N) and oxygen (O) isotopic values as shown by ambiguous or widely overlapping clusters (e.g., soil, manure/sewage, NH4+ fertilizer), which make straightforward assignments problematic. Accordingly, others16,17,18 recommend including additional isotopes (e.g., δ11B), other chemicals (e.g., pharmaceuticals, food additives), and/or biological markers (e.g., chlorophyll-a, fecal coliforms) to better differentiate sources of nitrogen pollution in freshwaters, particularly as related to municipal sewage, animal waste, industrial and atmospheric sources.

One common misperception in the interpretation of NO3− isotopes in aquatic systems is that its isotopic composition is a conservative tracer of N source(s), despite isotope fractionations during numerous biogeochemical transformations that significantly alter the δ15N and δ18O values of nitrate19. Instead, the δ15N and δ18O of NO3− in aquatic systems are more frequently a mixture of time and seasonally dependent N-source(s) and cumulative isotopic fractionations that occur during transport and biogeochemical transformation of dissolved N-species20,21,22. It is well-known that biological nitrate uptake (assimilation) and denitrification processes preferentially process the light isotopes, resulting in heavy isotope enrichment of residual nitrate14. Despite uncertainties in the assessment of sources of nitrate pollution using δ15N and δ18O, these stable isotopes retain fundamental source and process information to help decipher the complexity of N-biogeochemistry in aquatic environments21,23,24.

Common patterns of N-cycling processes using global spatial data sets of NO3− and its stable isotopes are established for soils and plant matter25,26,27,28, but larger-scale spatiotemporal studies of nitrate in rivers and adjacent aquifers are lacking20,29,30,31. Here we synthesized global river and adjacent shallow aquifer nitrate isotope (15N and 18O of NO3−) data sets from the scientific literature as well as new data (n = ~5200) to (a) provide a first-order global assessment of spatiotemporal patterns of NO3− and its isotopes in rivers and aquifers and (b) evaluate whether NO3− transformations are impacted by key environmental factors, such as temperature, climate, and seasonality. The investigation of the origin of nitrate pollution in rivers and groundwater at the local scale usually requires a deeper consideration of additional information such as N flux information from homogenized land-use types32; however, the detailed local evaluation was beyond the scope of this work. This synthesis provides a global foundational perspective to evaluate the potential and limitations of nitrate stable isotopes to track nitrogen pollution sources in aquatic systems, especially rivers, and to promote a data-based framework for further improving our understanding of the transformation mechanisms of nitrogen for sustainable management and remediation of N-contaminated waters.

Results

Overall assessment of global river and groundwater nitrate data sets

Concentrations of nitrate in rivers and adjacent aquifers deviated from normal distributions and were highly skewed (skewness index > 2), hence median values were retained33. Of all samples, rivers had substantively lower median NO3− concentrations (0.3 ± 0.2 mg L−1, n = 2902) compared to groundwater (5.5 ± 5.1 mg L−1, n = 2291) (Supplementary Table 1). Around 2% of river water samples exceeded the WHO (World Health Organization) threshold of 50 mg L−1 (as NO3−), whereas in groundwater the exceedance was far higher (~34%). The Kruskal–Wallis test indicated that nitrate concentration differences between rivers and groundwater subsets were significant (p-value < 0.05) (Supplementary Table 2).

The distribution patterns for δ15N and δ18O of NO3− in rivers and groundwater deviated from normal distributions (Supplementary Fig. 1). The isotope data were not as skewed as NO3− concentrations (skewness index < 2 for the number of records >30033), hence average values were retained to compare these isotope variables among them. Outliers were mostly extreme values falling within the upper and lower tails of the distributions of the Quantile–Quantile (Q–Q) plots (Supplementary Fig. 1). The δ15N values of nitrate in rivers were slightly lower (+7.1 ± 3.8‰, n = 2902) than aquifers (+7.7 ± 4.5‰, n = 2291, p-value < 0.05). Kendall14 presented an early compilation of nitrate isotope data in aquatic systems and reported groundwater nitrate has considerably higher δ15N values than rivers, however, others found more positive δ15N nitrate values in surface waters than nearby groundwater15. Other case studies affirmed higher mean δ15N values for NO3− in shallow groundwater compared to local rivers in urban waters (Manila, Bangkok, and Jakarta)34 or under different land-use-based catchments12. Overall, the δ18O values of nitrate in rivers and groundwaters worldwide were indistinguishable from each other (+2.3 ± 6.2‰, n = 2790, vs +2.3 ± 5.4‰, n = 2235, p-value > 0.05, respectively) (Supplementary Tables 1 and 2).

Isotope variation with latitude

The global data set included river and groundwater samples spanning from 67o N to 38o S LAT and from 145o E to 123o W LON (Fig. 1). There were more samples between 30° N to 60° N than other latitudinal ranges, reflecting the preponderance of data from North America, Europe, and East Asia. However, Central Asia, Central, and South America, Oceania, and Africa are poorly represented for nitrate isotope data sets due to a lack of water N pollution studies in those areas. The δ15N and δ18O of nitrate were negatively correlated with latitude for rivers (δ15N-NO3− = −0.07 ± 0.01 × Latitude + 10.14 ± 0.46, p-value < 0.05, n = 521, R2 = 0.07; δ18O-NO3− = −0.12 ± 0.02 × Latitude + 9.47 ± 0.79, p-value < 0.05, n = 506, R2 = 0.06) and groundwater (δ15N-NO3− = −0.15 ± 0.03 × Latitude + 14.59 ± 1.16, p-value < 0.05, n = 301, R2 = 0.08; δ18O-NO3− = −0.30 ± 0.03 × Latitude + 16.05 ± 1.43, p-value < 0.05, n = 299, R2 = 0.18) (Supplementary Fig. 2 and Supplementary Table 3). When data from the Northern Hemisphere were considered separately, the correlation between δ15N and δ18O of NO3− and latitude yielded a similar regression slope and intercept.

Fig. 1: Geographical distribution of δ15N and δ18O data for nitrate in river water and groundwater by number (#) of measurements by 0.1° latitude and year of sample collection.
figure 1

The size of each circle corresponds to the concentration of NO3-N in mg L−1. The light and dark gray bars indicate the data were retrieved from literature or the IAEA Coordinated Research Project (CRP), respectively. Most data were from 30° N and 60° N LAT and between 2010 and 2013.

Full size image

Relationship between 15N and 18O of NO3 −

Both river and groundwater nitrate revealed positive correlations between the δ18O and δ15N values (Fig. 2). The δ18O of NO3− was weakly but significantly correlated with δ15N for rivers (δ18O-NO3− = 0.64 ± 0.03 × δ15N-NO3− − 2.51 ± 0.26, p-value < 0.05, n = 2445, R2 = 0.13) and groundwater (δ18O-NO3− = 0.43 ± 0.02 × δ15N-NO3− − 0.99 ± 0.21, p-value < 0.05, n = 2172, R2 = 0.13). However, river waters had a significantly (p-value < 0.05) higher regression slope (0.64 ± 0.03) compared to groundwater (0.43 ± 0.02) (Fig. 2).

Fig. 2: δ18O vs δ15N of NO3− in rivers and groundwater.
figure 2

Symbol size denotes NO3-N concentration in mg L−1. The δ18O values are relative to SMOW and normalized to in situ H2O, as described in Venkiteswaran et al.31. The results of the linear regression models of panels a–d are in Supplementary Table 4.

Full size image

The influence of 18O in water and molecular O2 on the 18O of NO3− produced during nitrification is typically described using a simple isotope mass balance14:

$${\delta^{18}}{\rm{O}}{\hbox{-}}{\rm{NO}^{-}_3}=2/3 \times {\delta^{18}}{\rm{O}}{\hbox{-}}{\rm{H}_2}{\rm{O}} +1/3 \times {\delta^{18}}{\rm{O}}{\hbox{-}}{\rm{O}_2}$$
(1)

Eq. (1). Studies investigating nitrification in aquatic systems found the δ18O of NO3− does not necessarily conform to this model, even when the δ18O of O2 is assumed constant or in equilibrium with air35 (+23.5 − +24.2‰). To reduce the influence of δ18O from local water contributing to nitrate during its transformations, we normalized the δ18O-NO3− to the δ18O-H2O values from the same water sample wherever possible (i.e., “relative to in situ H2O” instead of “relative to SMOW”)31. The normalization to in situ water 18O resulted in an improved but weak correlation between δ18O-NO3−vs H2O and δ15N-NO3− for river waters (δ18O-NO3−vs H2O = −0.21 ± 0.02 × δ15N-NO3− −0.33 ± 0.19, R2 = 0.09, p-value < 0.05, n = 946) and groundwater (δ18O-NO3−vs H2O = −0.10 ± 0.02 × δ15N-NO3− −1.17 ± 0.06, R2 = 0.15, p-value < 0.05, n = 1138). The slopes of the “normalized” regression became negative and decreased from ~0.6 to ~−0.2 for rivers and from ~0.4 to ~−0.1 for aquifers. The river water samples from warmer tropical areas (e.g., Kenya, Ghana, Thailand) drove slopes to more negative δ18O and higher δ15N values. Results of the least-square linear regression analyses of the isotopic variables for river water and groundwater subsets are summarized in Supplementary Table 4.

Variation of nitrate isotopic composition with water temperature

Water temperature (Tw) is often used as a proxy for potential rates of the microbial activity or primary productivity and was used to consider its relationship to the isotopic composition of NO3− in rivers and aquifers36,37. In this assessment, the temperature was the river or groundwater temperature at the time of sampling without further considering overall climatic or seasonal aspects. Detailed results of least-square linear regression analyses of the nitrate isotopic variables for river water and groundwater are depicted in Supplementary Table 5. The δ15N of NO3− in rivers showed a weak but significant positive correlation with Tw (δ15N-NO3− = 0.12 ± 0.04×Tw + 6.55 ± 0.82, R2 = 0.03, p-value < 0.05, n = 253), but for groundwater, there was no correlation (δ15N-NO3− = 0.05 ± 0.05×Tw + 8.45 ± 1.03, R2 = 0.01, p-value > 0.05, n = 189) (Fig. 3). The δ18O of NO3− in rivers had a weak but significant positive correlation with Tw (δ18O-NO3− = 0.10 ± 0.04×Tw + 5.19 ± 0.90, R2 = 0.02 p-value < 0.05, n = 275), whereas in groundwater the correlation was very strong (δ18O-NO3− = 0.48 ± 0.04×Tw − 3.79 ± 0.95, R2 = 0.41, p-value <0.05, n = 170). Regressions against Tw using in situ normalized 18O nitrate values were conducted and resulted in less variable and slightly negative non-zero slopes like Fig. 2.

Fig. 3: Variation of δ18O and δ15N of NO3− in river waters and groundwater with water temperature (Tw) in °C.
figure 3

The linear regression of δ18O-NO3− (‰) vs Tw (°C) was more significant for groundwater (R2 = 0.41) than river waters (R2 < 0.02). The results of the linear regression models for panels a–d are in Supplementary Table 5.

Full size image

Seasonal variation of stable isotopes of NO3 −

To investigate the effect of seasonality on the isotopic composition of nitrate in rivers and groundwater, we used data from temperate climates (C), given that seasonality is strong in this climatic zone compared to the others (see “Methods” section, Fig. 4). The mean δ15N values of nitrate in rivers showed no difference between spring and fall (+6.9 ± 3.6‰, n = 583 to +6.7 ± 2.6‰, n = 606, p-value > 0.05), but a ~1‰ increase in δ15N in winter compared to all other seasons (+7.7 ± 3.8‰, n = 200, p-value < 0.05). The average δ18O values of nitrate in rivers revealed a ~2‰ decrease from spring to autumn (from + 2.2 ± 5.9‰, n = 530 to −0.1 ± 4.5‰, n = 562, p-value < 0.05), and a ~4.5‰ increase from autumn to winter (from −0.1 ± 4.5‰, n = 562 to +4.4 ± 4.5‰, n = 149, p-value < 0.05).

Fig. 4: Variation of δ15N and δ18O of nitrate in rivers and groundwater classified by season.
figure 4

The symbol size expresses the NO3-N in mg L−1. Higher average δ15N and δ18O values in river water nitrate are found in winter. Error bars show average (gray line), median (light to dark gray division), 1st and 3rd quartile, and whiskers extend to 1.5 interquartile range.

Full size image

Shallow aquifers showed no δ15N seasonality (p-value >0.05) in mean NO3− between spring, summer, autumn, and winter (+7.8 ± 4.4‰, +7.7 ± 3.9‰, +8.1 ± 3.8‰, and +8.1 ± 3.5‰, respectively). The average δ18O of NO3− in groundwater was similar in summer, autumn, and winter (+0.6 ± 4.8‰, +0.6 ± 4.6‰, and +0.1 ± 3.3‰, respectively) but was ~2‰ higher in spring (+2.1 ± 6.0‰, n = 268, p-value < 0.05) compared to all other seasons. However, considering that groundwater replenishment usually takes place in the wet period (October–April) in the northern hemisphere, we examined the difference in the isotopic values between the high recharge period and the low recharge period (May–September). There was no significant difference in δ15N and δ18O of NO3− between the two periods, although the δ18O of NO3− was slightly lower in the high recharge period (Supplementary Fig. 3). No clear-cut pattern was observed for the relationship between δ15N and NO3− in rivers or groundwater by season to confirm single predominant biogeochemical processes, including denitrification14 (Supplementary Fig. 4).

Variation of stable isotopes of NO3 with climate

The average δ15N of nitrate in rivers was significantly higher (~3.0‰) in tropical (A) (+10.2 ± 5.1‰, n = 137) vs temperate (C) (+7.5 ± 3.7‰, n = 834) or cold (D) (+7.4 ± 4.0‰, n = 252) climates (Fig. 5). The same pattern was observed for δ18O of NO3− for rivers, having ~2–3‰ higher values in tropical (A) (+6.6 ± 5.9‰, n = 137) vs temperate (C) (+4.4 ± 5.4‰, n = 722) and cold (D) (+3.8 ± 5.6‰, n = 275) climates. Arid climates (B) showed significantly higher δ15N values for nitrate (+8.9 ± 6.0‰, n = 74) compared to temperate (C) or cold climate (D), but similar δ18O values for nitrate (+5.5 ± 8.1‰, n = 74) compared to other climate types. However, the arid climate zone had fewer data compared to the other climate types due to a limited number of river water nitrate studies in those climates.

Fig. 5: Variation of δ15N and δ18O of NO3− in rivers and groundwater by climate.
figure 5

A: Tropical/Equatorial, B: Arid, C: Temperate, D: Cold. Symbol size expresses the NO3-N concentration (mg L−1). Tropical climates had the highest δ15N and δ18O values for nitrate in rivers. Error bars show average (gray line), median (light to dark gray division), 1st and 3rd quartile, and whiskers extend to 1.5 interquartile range.

Full size image

The shallow aquifers had the lowest δ15N values in temperate (C) (+7.8 ± 4.1‰, n = 986) and arid (B) climates (+7.3 ± 4.1‰, n = 522) and the highest in tropical/equatorial (A) (+8.9 ± 5.5‰, n = 44) and cold (D) climates (+9.0 ± 5.2‰, n = 24). However, these differences were indistinguishable (p-value > 0.05), hence no N-isotopic impact solely from climate could be found. Groundwater nitrate δ18O values exhibited the same pattern as rivers, with lower δ18O values in cold (D) climates (−0.1 ± 6.4‰, n = 24), that differed significantly from those in tropical (A) (+4.6 ± 6.0‰, n = 44) climates. The average δ18O of NO3− values in arid (B) (+3.4 ± 4.0‰, n = 505) and temperate (C) (+2.5 ± 4.9‰, n = 505) climates were significantly different from the other two climate types. Cold climates had fewer δ15N and δ18O records for groundwater compared to the other three climate types, which could bias this result. Rivers and groundwater showed no clear-cut relationship between δ15N and NO3− by climate type (Supplementary Fig. 5).

Discussion

The lower median NO3− concentrations in rivers compared to adjacent shallow aquifers can be explained by the decrease of agricultural N inputs at the industrialized countries of the Northern hemisphere due to more severe fertilization regulations, whereas many aquifers do not yet exhibit decreasing nitrate contents. This is due to much longer mean residence times for N in the unsaturated zone and the aquifer38, which is controlled by several factors such as the thickness and hydro-lithological characteristics of the unsaturated zone, the spatiotemporal variations in recharge rate, and the diffusive and dispersion processes through the nitrogen stock in the soil and the unsaturated zones. The lower median NO3− concentrations in rivers can also be explained by a complex combination of riparian and in situ attenuation and uptake by biological productivity, seasonal dilution (e.g., runoff, snowmelt), differential baseflow connectivity, and long-term dispersion processes, which lowers the ambient NO3− concentration in rivers. Given that NO3− sources are difficult to “fingerprint” using only stable isotopes, we suggest that the notable differences in NO3− concentrations between river waters and adjacent aquifers are attributable to unique locational differences in in situ N-cycling processes. As microbes have low biomass in most groundwater systems, NO3− attenuation, and cycling is achieved by closed-system nitrification and denitrification and uncommonly by DNRA or anammox29,39. Shallow aquifers are the recipients of NO3− transport from residual leaching from soils and the unsaturated zone, where NH4+ and other nitrogenous compounds are favorably oxidized but can also be adsorbed or reduced depending on soil retention capacity40,41,42,43. Pronounced denitrification effects have been reported for some unconfined aquifers30,44,45.

Regression analysis of δ18O vs δ15N for river waters and groundwater generated slopes suggesting an “apparent” denitrification trend (~0.5)3,31. However, denitrification is unlikely to be the only process for the observed δ18O:δ15N slope, which is expected to be canonical under closed-system conditions46. The lack of a clear pattern between the isotopes and NO3− concentration indicates that multiple processes that recycle nitrogen occur in river water systems. Granger and Wankel19 proposed that the deviations from the canonical slope of 1 for denitrification can be due to concurrent NO3‐ production catalyzed by nitrification and/or anammox.

A non-zero normalized slope suggests that H2O was not the only factor controlling the oxygen isotopic composition of nitrate (Eq.1). The large decrease in the slopes of the regression after normalization for in situ water confirmed a strong dependence of δ18O of NO3− on the δ18O of water when the NO3− is formed by oxidation of nitrite or ammonium42. However, in situ water normalizations do not consider potential large daily or seasonal variations in the δ18O of in situ gaseous or dissolved O235 involved in the N-cycling processes. Thus, the observed negative deviations of the slopes from “zero” suggest an influence of δ18O-O2 on the isotopic composition of nitrate. The latter may be more profound for rivers than aquifers. Investigations of δ18O-O2 in productive rivers show a preponderance of lower δ18O values (as low as +3.4‰), particularly during the growing season due to photosynthetic processes35,47,48. Tropical rivers typically show high rates of primary production, all-year-round growth, and less seasonal variation in solar irradiance than temperate latitudes49. On a diel basis, δ18O-O2 is lowest during the day due to the production of DO by photosynthesis from substrate water, whereas during the night δ18O-O2 can be higher than the atmospheric equilibrium value due to preferential consumption of light O2 by respiring organisms with inward diffusion of isotopically heavy O2 from the atmosphere50. Thus, the negative in situ normalized slopes of the δ18O regression line for rivers suggest some influence of 18O depleted O2, especially for highly productive rivers where oxygen supersaturation and 18O2 depletion and nitrification coexist. But even if the δ18O-O2 contribution was considered as a fixed atmospheric constant, any kinetic and equilibrium isotope fractionation during oxygen exchange between H2O and NO2− before final oxidation to NO3− may also result in deviations from a 2:1 mixing ratio line (δ18O-H2O:δ18O-O2)42. Therefore, the relationship between δ18O-NO3− and δ18O-H2O can be used to detect isotope fractionation effects related to the incorporation of oxygen atoms during ammonia oxidation to nitrite (NO2−), which readily exchanges oxygen atoms with in situ water during the oxidation of the latter to nitrate51. Conversely, for aquifers the δ18O-O2 values in recharge zone gas often deviate from the air (+23.5‰), but mostly in a positive direction up to >+39‰ as soil O2 is consumed by microbial respiration52,53. The substantial variability in the δ18O of NO3− associated with nitrification is also dependent on soil and water properties and is also a function of residence time54.

The difficulty to elaborate clear relationships between 18O of H2O, O2, and NO3− is also due to the temporal variability of the oxygen variables and underscores the need for further investigations to better understand N-dynamics related to cycling. Unfortunately, oxygen isotope normalization processes that consider all potential reactive oxygen pools (e.g., in situ H218O and 18O2) involved in N-cycling processes have not been fully investigated to our knowledge55. New studies measuring δ18O-NO3− and δ18O-H2O and δ18O-O2 in the same water are therefore required to better explain the patterns and variations of the oxygen isotope composition of nitrate in rivers and groundwater.

The δ15N and δ18O of NO3− in rivers showed a positive correlation with Tw likely reflecting the influence of higher microbial activity and productivity with temperature and thus the biogeochemical N processing rates. Higher δ15N values in rivers are likely driven by enhanced N-cascading. This means that due to increased temperature, metabolic rates accelerate, which progressively enriches 15N in NO3− compared to its source signature14. Average δ15N values of nitrate >7‰ typically suggest that nitrogen has undergone multiple recycling after its initial fixation from N2 (either biological or Haber–Bosch), as these processes bring reactive N in the plant–soil–water continuum with a 15N value of ~0‰56. An observed weak positive correlation of δ18O-NO3− with Tw and a negative correlation with latitude suggests that ambient δ18O-H2O has some impact on nitrate isotope values.

The pattern of decreasing δ15N-NO3− in rivers with latitude agrees with the N isotopic climatic patterns observed for soil organic matter (SOM)25. The pattern of higher δ15N-NO3− values with the water temperature in rivers is observed from cold (D) to tropical (A) climates. Amundson et al.25 observed a weakly positive correlation between soil organic matter (SOM) δ15N and Mean Annual Temperature (MAT) (slope = 0.13, R2 = 0.11, p-value < 0.1, n = 85) for arid and tropical zones. A decrease in SOM δ15N with lower MAT was seen at sites between about 10–40 degrees North or South albeit with a smaller data set (slope = −0.08, R2 = 0.09, p-value < 0.04, n = 49)27. Craine et al.26 observed a positive trend of SOM δ15N with MAT (slope = 0.18, p-value < 0.001, MAT > 9.8 °C) and a negative correlation with Mean Annual Precipitation (MAP) at a global scale that was attributed to soil and organic matter properties26. Martinelli et al.28 reported higher values for δ15N in N-rich soils and foliage from tropical forests compared to temperate forests. Brookshire et al.57 suggests that higher MAT enhances soil N-cycling as seen by 15N enrichment in soil organic matter and the NO3− concentration of tropical rivers.

The seasonality of δ15N and δ18O of nitrate in temperate rivers showed a pattern of isotopic enrichment of nitrate in spring and summer compared to autumn, which could be attributed to higher biological assimilation in surface water systems in the warmer seasons58. The δ15N values during winter compared to other seasons appeared to more closely reflect mixing of multiple organic-pollutant N-sources (e.g., sewage wastes, urban, and livestock wastewaters)9,24 as a result of reduced application and leaching of fertilizers, and restricted microbial and primary productivity during the cold season59. Additionally, in the cold season, the relative fraction of NO3− denitrified vs assimilated is higher, since increased N-loadings in winter or early spring runoff are often associated with high discharge events and snowmelt57. Elevated δ18O of NO3− during winter can also be attributed to δ18O-H2O and δ18O-O2 variation combined with variable soil/water and land-use properties of the different catchments. For example, high δ18O values in nitrate can be attributed to the influence of δ18O-O2 (Eq.1), especially in respiration-influenced ecosystems48, increased atmospheric NO3− deposition (>+30.0‰)14, especially during high discharge events, or artificial NO3− fertilizer (>20.0‰)14 mixed with nitrification-derived NO3− in soils and infiltrating pore waters, especially in agricultural areas54. Full equilibration or O-isotope exchange of NO2− with H2O before being oxidized to NO3− in soils, local respiration, and rapid redox cycling between NO3− and NO2− can also potentially increase δ18O-NO3− values19,54,60.

Shallow phreatic aquifers showed strongly dampened seasonal nitrate and isotopic variability compared to river waters. Given NO3− concentrations are generally buffered in aquifers over longer timeframes (e.g., years), the range of δ15N values in groundwater reflect longer-term climatic and land-use changes or established processes taking place in the soil and the unsaturated zone before NO3− leaches into the groundwater system30,61. The lack of a clear seasonal variation of δ15N and δ18O of NO3− in groundwater is attributed to nitrifying and mixing of nitrate sources introduced at different points or periods in time due to variable residence times (from months to years) in the unsaturated zone38,62. There was no significant difference in δ15N and δ18O of NO3− between the high and low recharge periods, although the δ18O of NO3− was slightly lower in the high recharge period. This could imply a change in recharge source water (e.g., snowmelt vs rains) or the 18O of O2, which combined affect the 18O values of microbially mediated nitrate. Because groundwater temperatures are usually stable and generally reflect MATs37 (and in situ δ18OH2O correlates with air temperature), it was unsurprising to see a correlation of δ18O of NO3− with Tw. This suggests we cannot discount the possibility that 18O variance of groundwater nitrate is controlled by additional processes or isotope fractionations affecting O2 involved in soil and vadose zone nitrification53. However, the relationship between δ18O in nitrate with latitude showed the highest oxygen isotope values tended to be at the warmer tropical/equatorial climates and can be attributed to a systematic change in the δ18O of associated H2O (Eq.1) with decreasing latitude, continentality or altitude63.

In summary, although rivers and shallow aquifers are primary and often well connected receiving environments for land-based anthropogenic N-pollution, shallow aquifers have five times higher median NO3− concentration than rivers, indicating that NO3− is more persistent in aquifers and has long residence times, given lower microbial activity and a limited number of possible N-removal process mechanisms (e.g., denitrification). River waters and groundwater receive nitrate with δ15N and δ18O values spanning the expected natural and anthropogenic ranges (from −10‰ to +26‰ for δ15N and −17‰ to +25‰ for δ18O)3, indicating a globally consistent suite of N-sources. We identified some global drivers of the isotopic variability of nitrate in rivers and groundwater, recognizing that this simplification required aggregating many possible diverse nitrate sources typically observed at the catchment or complicated by seasonal variability. A deeper analysis of latitude and climatic factors is recommended to better distinguish the role of soil properties, precipitation patterns, land-use, and agricultural practices. The δ18O vs δ15N relationships observed suggested combinations of multiple biogeochemical (e.g., nitrification, anammox) and mixing processes taking place in rivers and aquifers in different seasons, which are not easily distinguishable with infrequent or synoptic samplings. River water nitrate had higher δ15N and δ18O values in winter, suggesting a change in N-processing due to lowered temperatures and productivity and perhaps less alteration of the original isotopic composition of multiple organic-related N-sources. Across a wide range of global climate types, river waters and groundwater showed systematically higher δ15N and δ18O values in the arid and tropical/ equatorial climates, mirroring an increase in MAT. This relationship was also seen in the increase in δ15N of nitrate in rivers and groundwater with decreasing latitude and the positive correlation between δ15N and δ18O with Tw.

We have provided a first-order attempt to identify global patterns in the isotopic composition of nitrate related to common environmental factors, such as climate and season impacting N-cascade-related processes in rivers and associated shallow groundwaters. Our findings highlight the importance of water temperature as a driving force of biogeochemical microbial activity and productivity in rivers, however, the impact on the δ15N and δ18O of nitrate differed when examined on a seasonal or on a climate basis. Our findings suggest that higher frequency or seasonal sampling of river water systems is of critical importance to better and more deeply understand annual N-pollution dynamics, and winter sampling appears to better reflect the N-isotopic signals of the anthropogenic sources of nitrate. Further investigations to study N-cycling in river waters at a higher frequency are urgently needed (e.g., seasonal, diel) and systematic data collections of nitrate, particularly from the southern hemisphere climate types, such as arid and polar, for which data are limited or absent, are strongly recommended. Recent technological advances in isotopic assays of nitrate have vastly reduced the analytical and cost barriers to conduct high-frequency δ15N and δ18/17O of nitrate. Other key parameters, such as DO, redox potential, labile organic material, and solar radiation should also be systematically measured or considered to help better understand the controlling factors of the N-transformations in river waters, especially on a seasonal basis. Measuring δ18O of the water and dissolved oxygen from the same samples as δ18O of NO3− will allow a better assessment of N-cycling processes in rivers and groundwater. Understanding the fate of nitrogen in aquifers and the connectivity to rivers requires a fuller local understanding of hydrogeological processes, such as water source, flow paths, and water residence times (e.g., 3H tracers) that control nitrogen transit times and contribution to baseflow (e.g., retardation of N due to the recycling of N in organic matter and subsequent mineralization processes). Our synthesis affirms nitrate isotope techniques are useful to assess the origin of nitrogen pollution in rivers and groundwater. However, we caution that interpretations of nitrate isotopes may be complicated where nitrogen undergoes many biogeochemical transformations, especially in productive systems, that can mask the original isotopic source signal. Examination of N-cascading aspects are critically important to inform the implementation of beneficial land and agricultural management strategies aiming at mitigating increasingly serious nitrogen pollution in Earth’s rivers and aquifer systems. The clear isotopic linkage of the N-cascading in rivers to key climate parameters may have implications for future water quality management under changing climatic conditions.

Methods

Literature data collection

A literature search was conducted using online tools and bibliographic databases (e.g., Web of Science, Google Scholar, Scopus, etc.) for papers containing nitrate and isotope data from rivers and adjacent aquifers. We reported shallow aquifer data because many of the river studies contained data for adjacent aquifers owing to the possible connectivity of groundwater to riverine baseflow, although in most papers any connectivity of adjacent aquifers to the rivers was unconfirmed. Data compilations were restricted to journal articles where both δ15N and δ18O values of nitrate are available, and samples were then classified as River or Groundwater. Scientific papers focusing on precipitation, deep aquifers, soils, tap water, seawater, wastewater effluent, or artificial isotopic tracers were omitted. The list of cited data used is found in Supplementary References.

Data preparation

Established data preparation procedures resulted in a curated data set suitable for statistical analysis64. Data were extracted from tables and supplementary materials, maps, and text (Supplementary Fig. 6). If the original work did not provide location coordinates, approximate latitude and longitude were obtained using Google Earth (https://www.google.com/earth/). The data assessment affirmed that NO3− was the species linked to the reported isotope ratios (15N/14N and 18O/16O) along with verification of units (as NO3− or -N) and concentration (e.g., mg L−1, µmol L−1, μg L−1, etc.). Nitrate concentration data were converted to a common unit (mg L−1 -N) and sampling date format (e.g., 2016-03-01) (Supplementary Fig. 7). The 18O/16O of in situ H2O associated with a nitrate sample, if available, was also included in the data set. The ratios of 15N/14N and 18O/16O of NO3− and 18O/16O of H2O were reported in δ notation in units of per mille (‰), where δ = (Rsample/RAIR or SMOW − 1) and R is the ratio of 15N/14N or 18O/16O in N2 in atmospheric air (AIR) or Standard Mean Ocean Water (SMOW), respectively.

Seasons were identified separately for each Hemisphere (N vs S) using the reported sampling date (month) to ensure consistency. Each sample location was coded using the Köppen climatic classification65. Due to limited data from high latitudes, river and groundwater data were classified into four climatic groups; A: Equatorial/Tropical climate, B: Arid climate, C: Temperate climate, D: Cold climate. Additional categorical and site information associated with the nitrate and isotopic data was obtained from the publications and is included in Supplementary Table 6. Important variables, such as Eh, dissolved O2 (DO), dissolved organic nitrogen (DON), precipitation per catchment, Total Nitrogen (TN), that might contribute to a better understanding of the global N isotope variations of nitrate were either lacking or unavailable. For example, DO values were only available in 7.7% of the river waters data set and thus were excluded from further processing.

Data curation

Excessively high or suspicious values of NO3− and other N-species concentrations were verified by a careful review of the original paper. In some cases, authors were contacted to provide clarification. The δ15N and δ18O values of nitrate were sorted to avoid duplication from publications reusing the same data sets26. Almost all groundwater data (~95%) was from adjacent shallow phreatic aquifers of <20 m depth below ground surface.

Final data set and analytical aspects

The data set comprised 5194 sample points from 70 papers (~86% of the whole data set, Supplementary References) plus new isotopic results. There were 2902 nitrate and associated isotope analyses from rivers and 2291 analyses from nearby shallow groundwater. Most samples were from between 2002 and 2016; fewer data exist from 1968 to 1990s due to analytical barriers particularly for δ18O66.

A total of 374 new river water nitrate and isotope samples from 10 countries were analyzed in the IAEA Isotope Hydrology Laboratory using the Ti(III) reduction method13 to reduce nitrate to N2O gas in septum vials. The N2O in the headspace was measured for 15N and 18O using a CF-IRMS (Isoprime-100 Trace Gas Analyser). The analytical uncertainties were ±0.2‰ and ±0.4‰ for δ15N/δ18O, respectively. For Quality Control, the water samples, that do not give the expected N2O yield on the IRMS based on the determined NO3− concentration, were repeated until they agree within 95% of the expected target before results are accepted. Oxygen isotopes (δ18OH2O) of the same water samples were analyzed by laser spectroscopy at the IAEA using a Los Gatos Research Liquid-Water Isotope Analyzer (TLWIA-912). The analytical uncertainty was ±0.1 for δ18OH2O. Chemical analysis of new samples included nitrate following standard discrete analyser methods67.

Statistical evaluation

To evaluate the data, graphical plots by sample type (river or groundwater) were employed along with frequency distribution diagrams (histogram, Q–Q plot, P–P plot, boxplot) and metrics such as Pearson skewness. To check for normality, Shapiro–Wilk, and Kolmogorov–Smirnov tests were applied. The mean and standard deviation (SD) were used to describe the variation of variables with a skewness index of <2, whereas median and the median absolute deviation (MAD) were used otherwise33. A single-factor analysis of variance (ANOVA) and Kruskal–Wallis non-parametric analysis of variance were performed to determine if differences between the means or the medians of the groups were significant. A conventional probability value (p-value) of 0.05 was used to indicate significance for differences between river and groundwater and seasonal subsets68.

A least-square linear regression model of type I was applied to investigate relationships between δ15N and δ18O of nitrate to other variables. The least-square linear regression type II was tested by assuming that the two isotopic variables (δ15N and δ18O in nitrate) were dependent69. However, the difference in the slopes and intercepts were small in most cases (<10%), which allowed us to use model-I patterns for nitrate isotopes in both rivers and groundwaters. The results of regression models of type II are found in Supplementary Methods (Supplementary Tables 7–8 and Supplementary Figs. 8–15). All linear regression models were evaluated by examining the goodness of fit and significance of the slopes. All statistical tests were done using R v.3.3.270.