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

Brittany is a 27 000 km² 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⁻²) 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 (Agreste Bretagne 2019). Industrial agriculture began to develop there after 1945 and the European Union's Common Agricultural Policy (1962), this development was more rapid and intense there than in other French regions (Gascuel-Odoux et al 2010). The agricultural N surplus was only 33 kg ha⁻¹ yr⁻¹ in 1955, peaked at 104 kg ha⁻¹ yr⁻¹ in 1989 and decreased to 44 kg ha⁻¹ yr⁻¹ in 2015 (Poisvert et al 2017). Organic and mineral N fertilizers represented 69% of N inputs in 1955, 85% in 1989 and 80% in 2015. The rest came from biological fixation by legumes and atmospheric deposition.

Zoom In Zoom Out Reset image size

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.eau-loire-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². Because the study catchments ranged from 104–10 135 km², 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 et al 2015)) and no long-term point-source data were available.

The 16 catchments selected ranged from 104–10 135 km², had specific cumulative discharge of 195–689 mm yr⁻¹ 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).

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 long-term 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:

$$\begin{equation*}Retention = 1 - \frac{{\mathop \sum \nolimits_{i = 1976}^{2015} Flux\,i}}{{\mathop \sum \nolimits_{i = 1976}^{2015} Surplus\,i}}\end{equation*}$$

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 long-term N surplus—N load data in a 270 km² 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)

Cumulative N surplus inputs from 1976–2015 ranged from 1873–4424 kg N ha⁻¹ among the 16 catchments, with a mean of 2960 kg N ha⁻¹( i.e. 74 kg ha⁻¹ yr⁻¹ on average for the 40-year period) (table 2). During the same period, cumulative N flux ranged from 398–1820 kg N ha⁻¹, with a mean of 911 kg N ha⁻¹( i.e. 23 kg ha⁻¹ yr⁻¹). 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⁻¹, with a mean of 2050 kg N ha⁻¹( i.e. 51 kg ha⁻¹ yr⁻¹ 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 (figure 2). The N surplus increased from 18–44 kg ha⁻¹ yr⁻¹( mean of 34 kg ha⁻¹ yr⁻¹) in 1976 to a maximum of 74–176 kg ha⁻¹ yr⁻¹( mean of 113 kg ha⁻¹ yr⁻¹) in 1988. It then decreased to 35–74 kg ha⁻¹ yr⁻¹( mean of 49 kg ha⁻¹ yr⁻¹) in 2015. The flow-normalized N flux increased from 2–23 kg ha⁻¹ yr⁻¹( mean of 11 kg ha⁻¹ yr⁻¹) in 1976 to a maximum of 12–57 kg ha⁻¹ yr⁻¹( mean of 30 kg ha⁻¹ yr⁻¹) in 1996. It then decreased to 11–41 kg ha⁻¹ yr⁻¹( mean of 21 kg ha⁻¹ yr⁻¹) in 2015.

Zoom In Zoom Out Reset image size

The mean ± standard deviation of total N content measured in the topsoil was 2.20 ± 0.87 g kg⁻¹ from 2000–2004 and 2.29 ± 0.77 g kg⁻¹ from 2010–2014, indicating a statistically significant increase of 0.09 g kg⁻¹ in a 10-year period (Student's t-test, p < 0.05). Assuming a typical soil bulk density of 1.2 t m⁻³( Ellili et al 2019), this increase in total N content from 7920 ± 3132 to 8244 ± 2772 kg N ha⁻¹ represented a mean accumulation of 32.4 kg ha⁻¹ yr⁻¹ in the first 30 cm of agricultural topsoil.

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).

Zoom In Zoom Out Reset image size

The response time (10% transit time) ranged from 1.2–10.8 years, with a median 5%–95% credibility interval of 9.4 years (table 2). The median transit time ranged from 3.8–15.7 years, with a median 5%–95% credibility interval of 7.4 years. The estimated equilibrium time (90% transit time) ranged from 12.7–25.2 years, with a median 5%–95% credibility interval of 20.2 years. The mode (peak) of the TTD ranged from 1.6–13.8 years, with a median 5%–95% credibility interval of 9.8 years. Despite of a variety of TTD shapes (figure 4), response times, median transit times, mode and equilibrium time were strongly correlated (r = 0.83–0.99, p < 0.05, figure SI 3).

Zoom In Zoom Out Reset image size

The long-term (1976–2015) 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).

Zoom In Zoom Out Reset image size

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, Boyer et al 2002, Howarth et al 2006). 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⁻¹ yr⁻¹, and the mean of all agricultural exports (harvest) was 76 kg N ha⁻¹ yr⁻¹, resulting in an estimated N surplus of 68 kg N ha⁻¹ yr⁻¹( 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 (Howarth et al 2006), 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⁻¹ yr⁻¹) represented 64% of the estimated annual retention observed over a 40-year period (51 kg ha⁻¹ yr⁻¹). 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² Kervidy-Naizin research catchment in Brittany, where denitrification was estimated to represent only 17% of retention (Durand et al 2015) 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).

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 (Van Meter et al 2017, 2018, 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).