Long-term nitrogen retention and transit time distribution in agricultural catchments in western France

Elevated nitrogen (N) concentrations have detrimental effects on aquatic ecosystems worldwide, calling for effective management practices. However, catchment-scale annual mass-balance estimates often exhibit N deficits and time lags between the trajectory of net N inputs and that of N riverine export. Here, we analyzed 40-year time series of N surplus and nitrate-N loads in 16 mesoscale catchments (104–10 135 km2) of a temperate agricultural region, with the aim to (1) investigate the fate of the ‘missing N’, either still in transit through the soil—vadose zone—groundwater continuum or removed via denitrification, and (2) estimate the transit time distribution of N by convoluting the input signal with a lognormal model. We found that apparent N retention, the ‘missing N’, ranged from 45%–88% of then N net input, and that topsoil N accumulation alone accounted for ca. two-thirds of this retention. The mode of the nitrate-N transit time distribution ranged from 2–14 years and was negatively correlated with the estimated retention. Apparent retention was controlled primarily by average runoff, while the transit time mode was controlled in part by lithology. We conclude that the fate of the soil ‘biogeochemical legacy’, which represents much of the catchment-scale ‘missing N’, is in our hands, since the N accumulated in soils can still be recycled in agroecosystems.


Introduction
Excess reactive nitrogen (N) in the environment causes several detrimental effects, including acidification, climate change and eutrophication, particularly in marine ecosystems (Galloway et al 2003, Pinay et al 2017. Coastal and near-coastal eutrophication contributes to the formation of algal blooms and hypoxia, and these conditions threaten biodiversity, tourism, and fisheries (Turner et al 2008, Andersen et al 2017, Wang et al 2016. Besides alteration to ecosystem structure and function, massive algal blooms, also called 'green tides' , can even result in a direct human health risk with documented cases of mortality due to exposure to toxic H 2 S gas emissions resulting from decomposition of excess algae in shallow tidal areas (Smetacek andZingone 2013, Le Moal et al 2019).
N inputs from agriculture represent most of the reactive N that enters streams and rivers, primarily as nitrate, in economically developed countries (Howarth and Marino 2006, De Vries et al 2011, Li et al 2019. Efforts to limit agricultural N losses to water include reduction of N surplus and maximization of N removal in natural and engineered wetlands (Seitzinger et al 2006). N mass-balance estimates at catchment and regional scales often reveal a N deficit, sometimes called 'missing N' (Van Breemen et al 2002), which commonly exceeds 50% of net inputs (Howarth et al 1996, Alexander et al 2002, Aquilina et al 2012, Lassaletta et al 2012, Ehrhardt et al 2019. Knowing the fate of this missing N has important implications for assessing trends in water quality and managing global and local N cycles. The missing N may have accumulated in the catchment, building a legacy source that causes a time lag in its transfer to surface water, or it may have been removed via denitrification (Chen et al 2018).
The N legacy of a catchment includes a 'biogeochemical legacy' , when N accumulates in soil organic matter and may leach several years after on-field application, and a 'hydrological legacy' , which corresponds to the transit time of water in catchments. Both types of legacies lead to time lags of several years to decades between implementation of measures to curb N losses and decreasing N concentrations in rivers Basu 2017, Vero et al 2017). In a long-term tracer experiment, Sebilo et al (2013) observed that 12%-15% of isotopically labelled N still remained in an agricultural soil 30 years after application, 61%-65% having been taken up by plants and the rest emitted into the environment. Van Meter et al (2016) analyzed long-term soil N content in the Mississippi River Basin and estimated a net accumulation rate of 25-70 kg ha −1 yr −1 . On the other hand, N 'hydrological legacy' is assumed to follow the transit time distribution (TTD) of water, which can be investigated with conservative tracers and both empirical and numerical models (Fovet et al 2015, Hrachowitz et al 2016. Denitrification, the primary process by which N surplus can be permanently removed from a catchment, occurs in soils (Oehler et al 2007), groundwater (Kolbe et al 2019), wetlands (Montreuil et al 2010) and river networks (Pinay et al 2015). Denitrification produces both inert N 2 and the greenhouse gas N 2 O. Multiple studies have observed denitrification at local scales, but upscaling to the catchment scale remains a challenge (Seitzinger et al 2006).
Few studies have simultaneously investigated both N retention and TTD, seemingly because doing so requires having detailed agricultural, soil and hydrological data to quantify budgets accurately and having these data over long periods. Such datasets have been rarely available until recently (Howden et al 2010, Ehrhardt et al 2019. In this article, we analyzed long-term  data from previously estimated N surpluses (Poisvert et al 2017) and newly estimated riverine nitrate export and soil N content variation. The research objectives were to (1) estimate N retention and TTD and (2) relate these estimates to geographic variables (catchment area, annual runoff, lithology).
The study sites include 16 catchments in the Brittany region (western France), which is one of the emblematic areas for N-derived coastal eutrophication in Europe (Smetacek and Zingone 2013). It also has few N point sources (Gascuel-Odoux et al 2010), has experienced large variations in its N surplus in the past few decades (Dupas et al 2018) and has relatively quick catchment response times (Fovet et al 2015), which makes it an ideal study area for an empirical modeling study.

Study area
Brittany is a 27 000 km 2 region located in the Armorican Massif of northwestern France (figure 1). Its geology is dominated by igneous and metamorphic rocks (granite, schist, micaschist), and its topography is relatively flat, with elevation ranging from 0-385 m above sea level (Gascuel-Odoux et al 2010). Its climate is temperate oceanic, with a mean annual temperature of 12 • C and a large rainfall gradient that increases from 700 to 1300 mm from east to west (Frei et al 2020). Its hydrology is characterized by shallow aquifers in weathered and fissured layers of bedrock, a high river density (1 km km −2 ) and riparian wetlands (partly cultivated) that cover ca. 20% of the soil map (Aquilina et al 2012, Abbott et al 2018. Agriculture represents 80% of the land use and is dedicated mostly to animal production. Brittany contains 21%, 58%, 33% and 42% of France's dairy cows, pigs, layer chicken and broiler chicken, respectively, on only 6% of France's agricultural area

Monitoring data and catchment characteristics
The study focused on 16 rivers in Brittany in which nitrate monitoring started in 1976 or earlier and in which nitrate was measured near a gauging station that provided daily discharge data (figure 1). Water quality was monitored on a regular interval (monthly to bimonthly): the number of years from 1976-2015 with at least six sampling dates ranged from 23-40 among the catchments. All nitrate and discharge data were collected by water authorities for regulatory purposes; they are publicly available at http://osur.eauloire-bretagne.fr/ and http://hydro.eaufrance.fr/.
N surplus data were estimated by Poisvert et al (2017) based on agricultural statistics for each 'commune' , a French administrative unit with a mean area of 15 km 2 . Because the study catchments ranged from 104-10 135 km 2 , we ignored uncertainties related to farms with fields in multiple communes or that transferred organic N to neighboring farms. Poisvert et al (2017) calculated the annual N surpluses using a land system budget (Oenema et al 2003, De Vries et al 2011. N inputs consist of mineral and organic fertilizer (farmyard manure and slurry), biological N fixation and atmospheric N deposition. N output are the sum of N exported by each crop type. In this study, N surplus was considered as a net diffuse N source for the catchments. We ignored point-source inputs because they represent a small percentage of total N fluxes in France (i.e. 2% from 2005-2009 (Dupas ) and no long-term point-source data were available.
The 16 catchments selected ranged from 104-10 135 km 2 , had specific cumulative discharge of 195-689 mm yr −1 and a percentage of agricultural land use of 73%-92%. Two of them are dominated (i.e. >66%) by granite bedrock, nine are dominated by schist/mica-schist and five of them have mixed lithology (table 1).

Flux estimation and transit time modeling
We first filled data gaps in the discharge time series using the geomorphology-based SIMFEN model (de Lavenne and Cudennec 2019). Developed for Brittany, SIMFEN simulates discharge in ungauged catchments by transposing net rainfall from neighboring gauged 'donor' catchments and convoluting this net rainfall via a geomorphology-based transfer function. The percentage of discharge gaps ranged 1%-23% (median = 2%) and the Nash-Sutcliffe (NS) efficiency for catchments with >5% missing discharge ranged 0.73-0.84. We then interpolated the monthly concentration time series to a daily time step and then to annual time step using the 'Weighted Regressions on Time, Discharge, and Season' (WRTDS) provided by the EGRET package (Hirsch et al 2010, Hirsch and Decicco 2019) of R software. WRTDS is a load estimation method that use time, discharge and season as explanatory variables, and re-estimates statistical dependencies at each time step using the most relevant dates to represent the concentration dynamics and provide accurate load estimates (Zhang and Hirsch 2019).
Previous studies in Brittany have shown (1) large interannual variations in discharge, with annual discharge in a wet year up to three times that in a dry year (Gascuel-Odoux et al 2010), and (2) that nitrate concentrations increased in wet years, causing interannual variability in nitrate load to exceed that of discharge (Dupas et al 2018, Mellander et al 2018. The present study focused on long-term trends in nitrate loads due to changes in N surplus, not to interannual variations in the hydroclimate. Consequently, we considered the EGRET flow-normalized annual load rather than the actual estimated load. Using a flow-normalized load throughout this study was acceptable because, despite interannual variations, none of the 16 study catchments exhibited a longterm monotonous discharge trend from 1976-2015 (Mann-Kendall test, p > 0.05).
The input-output assessment considered N surplus as an input and flow-normalized nitrate flux as an output. N retention from 1976-2015 was estimated as: The estimated retention included both irreversible removal via denitrification and long-term reversible storage. The TTD model consisted of a dynamic convolution model of the N surplus not retained in the catchments, on an annual time step. The mathematical function used was a lognormal distribution model with a mean (µ) and standard deviation (σ). Kirchner et al (2000) found that distributions with tails longer than those of exponential distributions (e.g. lognormal or gamma) were suitable for modeling solute transport in catchments, and Ehrhardt et al (2019) successfully fitted a lognormal distribution to longterm N surplus-N load data in a 270 km 2 temperate catchment. In essence, we fitted the TTD parameters so that the N surplus convolved with the optimal distribution matches the observed N loads at the outlets. We calibrated the TTD model in a generalized likelihood uncertainty estimation (GLUE) framework (Beven and Freer 2001). We performed 5000 model runs assuming a uniform distribution in the range [0:5] for log µ and [0:2] for log σ. We considered models that exceeded a NS efficiency coefficient of 0.60 to be 'behavioral' , and calculated 5%-95% credibility intervals of outputs of the NS-weighted behavioral models.
Statistical analysis of TTD model outputs consisted of calculating correlations between estimated retention, transit time mode, transit percentiles (10%, 50% (median) and 90%) and selected catchment properties. The 10% transit time was selected as the time needed to detect a response to a change in N surplus, the 50% transit time represents the time needed to transfer half of the N surplus in a given year, and the 90% transit time represents the time needed to reach near-equilibrium. In statistical analyses of transit time percentiles, we considered the NS-weighted median of the behavioral models. The catchment properties selected were catchment area, dominant lithology and long-term runoff.
We assessed the fate of the retained N (i.e. removal via denitrification, or accumulation in soil, the vadose zone and groundwater) in relation to estimates of N accumulation from available soil N test data. Specifically, we compared total N content (modified Kjeldahl method ISO 11 261:1995) in the topsoil (0-30 cm) measured from 30 484 and 36 764 soil tests in Brittany from 2000-2004 and 2010-2014, respectively. Because these soil data are protected by statistical confidentiality, spatial coordinates of the sampling points were not available (Saby et al 2014). The objective of this approximate quantification was to compare the magnitudes of plausible N accumulation and quantified N retention during the study period.
All the statistical analyses and figures were done with the R software v. 3.5.1 (R Core Team 2019)

Long-term nitrogen budget (1976-2015)
Cumulative N surplus inputs from 1976-2015 ranged from 1873-4424 kg N ha −1 among the 16 catchments, with a mean of 2960 kg N ha −1 (i.e. 74 kg ha −1 yr −1 on average for the 40-year period) (table 2). During the same period, cumulative N flux ranged from 398-1820 kg N ha −1 , with a mean of 911 kg N ha −1 (i.e. 23 kg ha −1 yr −1 ). The coefficient of variation of cumulative flux exceeded that of cumulative surplus (42% and 26%, respectively), and the former was similar to that of long-term runoff (42%). Cumulative N flux was not correlated with cumulative N surplus, but it was positively correlated with long-term mean runoff (r = 0.87, p < 0.05). Cumulative N retention ranged from 45%-88% among catchments, with a mean of 67%. This retention represented a range of 864-3505 kg N ha −1 , with a mean of 2050 kg N ha −1 (i.e. 51 kg ha −1 yr −1 on average for the 40-year period).
Both the N surplus and the flow-normalized N flux varied by a factor of three during the study period ( The mean ± standard deviation of total N content measured in the topsoil was 2.20 ± 0.87 g kg −1 from 2000-2004 and 2.29 ± 0.77 g kg −1 from 2010-2014, indicating a statistically significant increase of 0.09 g kg −1 in a 10-year period (Student's t-test, p < 0.05). Assuming a typical soil bulk density of 1.2 t m −3 (Ellili et al 2019), this increase in total N content from 7920 ± 3132 to 8244 ± 2772 kg N ha −1 represented a mean accumulation of 32.4 kg ha −1 yr −1 in the first 30 cm of agricultural topsoil.

Transit time distribution modelling
The GLUE approach yielded 166-1044 behavioral parameter sets among the 16 catchments. Maximum NS efficiencies ranged from 0.71-0.96 among the catchments (figures SI 1 and 2 (available online at https://stacks.iop.org/ERL/15/115011/mmedia)). The modeled N flux behaved similarly to the observed flow-normalized N flux in the long term (figure 3).

Relation to geographic variables and short-term concentration metrics
The long-term  retention rate was negatively correlated with long-term mean runoff (r = −0.88, p < 0.05), meaning that higher retention was observed in catchments with lower mean runoff ( figure 5(a)). The mode of the TTD was not significantly correlated with long-term runoff. Granite and mixed-lithology catchments generally had longer transit times than schist catchments ( figure 5). However, schist catchments had high variability in the TTD mode, and the mixed-lithology catchment Couesnon clustered with its schist-dominated neighbors despite having a substantial percentage (46%) of granite. This suggests that our classification of lithology into three simplified categories explained only some of the variability in the TTD mode. Neither long-term retention nor the TTD mode were significantly correlated with catchment area (figure SI 4).

Fate of the missing nitrogen
The range of long-term retention estimated in 16 Brittany catchments (45%-88%) is similar to several  previous estimates in temperate regions. Billen et al (2011) estimated that N retention from European catchments that discharge into the Baltic, North, Mediterranean and Black Seas, and the Atlantic Ocean, averaged 78% of the net anthropogenic input.
Estimates of retention in the northeastern United States also range from 70%-80% (Howarth et al 1996. Studies at more local scales found N retentions of 92% in a Mediterranean catchment (Lassaletta et al 2012), 88% in a temperate continental catchment (Ehrhardt et al 2019) and 53 ± 24% in 160 French catchments using short-term data from 2005-2009(Dupas et al 2015. We note, however, that these studies use different N input data sources, which explains part of the differences together with differences in catchment properties. This long-term estimate of N retention is more reliable than the previous estimate from a 5-year period in France, because estimates from short-term studies have a higher risk of being influenced by transient states than long-term studies (Dupas et al 2015). This long-term estimate is still subject to several sources of uncertainty, however, including estimation of the N surplus, estimation of N load and no consideration of point sources. By construction, estimating the N surplus is uncertain because it is calculated as a difference: in Brittany from 1976-2015, the mean of all diffuse inputs (fertilizers, atmospheric deposition, biological fixation) was 144 kg N ha −1 yr −1 , and the mean of all agricultural exports (harvest) was 76 kg N ha −1 yr −1 , resulting in an estimated N surplus of 68 kg N ha −1 yr −1 (Poisvert et al 2017). Because of the calculation by difference, uncertainty of 10% in the input term would cause uncertainty of 21% in the estimate of N surplus. This effect should encourage researchers to focus more on relative differences in retention among catchments rather than their absolute values. A previous study estimated that nitrate represented 88% of total N load and point sources represented a mean of 2% of N load on average in 38 catchments in Brittany from 2005-2009(Dupas et al 2015. Because the data on non-nitrate N and point source inputs were not available for our long-term study period, we decided not to correct our N load estimate and not include point-source inputs, assuming that doing so would have as little effect from 1976-2015 as from 2005-2009. Like previous studies in France (Dupas et al 2015) and North America , we found that runoff predicted N retention well (r = −0.88). Howarth et al (2006) hypothesized that high runoff could decrease water residence times in riparian wetlands and low-order streams, thus decreasing denitrification. Dupas et al (2015) also pointed out that, due to the short study period in their own study and that of Howarth et al (2006), this observation could also be due to catchments being in a transient state, with wetter catchments responding faster to decreasing N inputs than drier catchments. The fact that we observed the same relation with runoff over a 40-year period strengthens the hypothesis of Howarth et al (2006) that wet areas are less favorable to permanent N removal via denitrification or long-term accumulation in soils. The lack of correlation between N retention and catchment area results from the secondary role played by riverine processes on N retention compared to hillslope processes in the Brittany region (Casquin et al 2020).
The mean annual topsoil N accumulation observed over a 10-year period in Brittany (32 kg ha −1 yr −1 ) represented 64% of the estimated annual retention observed over a 40-year period (51 kg ha −1 yr −1 ). Biogeochemical legacy could therefore account for much of the apparent retention, as also hypothesized by Worrall et al (2015) and observed by Van Meter et al (2016). Both Van Meter et al (2016) and Worrall et al (2015) nonetheless concluded that N accumulated primarily in subsoils, suggesting that our estimate of legacy soil N in the topsoil may even underestimate the actual biogeochemical legacy. Besides the lack of data on subsoils, several factors make our estimate uncertain: potential inability of a 10-year period to represent long-term N accumulation, possible bias in the location of soils sampled from 2000-2004 and 2010-2014, and lack of data on non-agricultural soil. Gaining access the spatial location of soil samples in the confidential data for research purposes would allow us to correct potential bias in sampling location and test whether areas with higher observed N retention also accumulate more N in the soil.
Despite these uncertainties, the conclusion that N accumulation exceeds N removal agrees with two independent studies in the 5 km 2 Kervidy-Naizin research catchment in Brittany, where denitrification was estimated to represent only 17% of retention  and <10% of N export (Benettin et al2020). This legacy soil N is a potential threat to aquatic ecosystems but also a potential resource for crops that could allow farmers to reduce fertilizer applications. Because much of the N retained is likely to remain in the rooting zone, opportunities exist to recycle it in agrosystems by favoring uptake of the mineralized N by crops. Accumulation of organic N in soils also provides opportunities to sequester organic carbon and contribute to climate change mitigation (Bertrand et al 2019).

Nitrate transit time distribution
Many studies throughout the world have documented a lag time of several years to decades between decrease in N net inputs and improvement in water quality (see Chen et al (2018) for a review), although few have related this 'delay' to TTD percentiles. Our estimated mode and median of nitrate TTD of ca. 10 years is similar to those of previous studies on water TTD using chemical tracers and physical-based models in Brittany (Molenat et al 2002, Martin et al 2004, Ayraud et al 2008, Aquilina et al 2012, Fovet et al 2015. Ehrhardt et al (2019) also found modes of 19 and 12 years in two nested temperate catchments, using a similar method with agricultural N surplus as an input and assuming a lognormal transfer function.
As Van Meter and Basu (2015) noted, '[the existence of] biogeochemical nutrient legacies increases time lags beyond those due to hydrologic legacy alone' . Our estimated N TTD matched those previously estimated for water TTD because we assumed that the missing N was permanently retained or removed from the catchment, although the N that accumulated in soil could also be considered to be in transit (Howden et al 2011, Sebilo et al 2013, Worrall et al 2015. As a result, the actual tail of N TTD in Brittany is probably be thicker than our estimate suggests, which explains our focus on the TTD mode rather than on the equilibrium time. With more detailed information about the variation in N storage in soils and the vadose zone during the study period, and more accurate estimates of mineralization and denitrification rates, it would be possible to decrease the risk of equifinality in models such as ELEMeNT, which explicitly decomposes soil legacy, hydrological legacy and denitrification , Ballard et al 2019.
Comparison of the catchments showed that granite and mixed-lithology catchments generally had longer transit times than schist catchments, although the latter had high variability in transit times. Previous groundwater dating studies have observed the same trend of older water in granite catchments than in schist catchments (Martin et al 2006, Ayraud et al 2008.

Conclusion
This study used long-term time series of N surplus and nitrate-N loads to estimate N retention and TTD in 16 catchments in Brittany, France. Estimated N retention ranged from 45%-88%, and the TTD mode ranged from 2-14 years, which agrees with previous studies in similar contexts. We used a two-step approach, first estimating N retention and then convoluting the remaining, exported N with a lognormal distribution. It provided a long-term mean retention rate and TTD that were useful for capturing the average behavior of catchments during the 40-year study period. We acknowledge the variable nature of both retention and TTD over time and emphasize that our modeling approach should not be used to predict future N loads in Brittany, because the soil N accumulation documented in this study may reach saturation in a not too distant future. Apparent retention was controlled primarily by average runoff, and the TTD mode was controlled in part by lithology. Furthermore, observed variations in soil N storage suggests that the biogeochemical legacy of soils accounted for ca. two-thirds of the catchmentscale N retention. From a methodological perspective, the uncertain fate of this biogeochemical legacy, whether in transit through the catchment, denitrified or taken up by crops, influences the estimated TTD. From a management perspective, this biogeochemical legacy is both a potential threat to aquatic ecosystems and a potential resource that could be recycled in agroecosystems.

Acknowledgments
RD and SE designed the study and carried out the analysis. RD wrote the paper. All authors discussed and interpreted the results and contributed to the text. We would like to thank Dr Robert D Sabo and two anonymous reviewers for their valuable and constructive comments.