# A hillslope-scale aquifer-model to determine past agricultural legacy and future nitrate concentrations in rivers

Luca Guillaumot, Jean Marçais, Camille Vautier, Aurélie Guillou, Virginie Vergnaud, Camille Bouchez, Rémi Dupas, Patrick Durand, Jean-Raynald de Dreuzy, Luc Aquilina

► To cite this version:

Luca Guillaumot, Jean Marçais, Camille Vautier, Aurélie Guillou, Virginie Vergnaud, et al. A hillslope-scale aquifer-model to determine past agricultural legacy and future nitrate concentrations in rivers. Science of the Total Environment, 2021, 800, pp.149216. 10.1016/j.scitotenv.2021.149216. insu-03320536

HAL Id: insu-03320536

https://insu.hal.science/insu-03320536

Submitted on 16 Aug 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
Journal Pre-proof

A hillslope-scale aquifer-model to determine past agricultural legacy and future nitrate concentrations in rivers

Luca Guillaumot, Jean Marçais, Camille Vautier, Aurélie Guillou, Virginie Vergnaud, Camille Bouchez, Rémi Dupas, Patrick Durand, Jean-Raynald de Dreuzy, Luc Aquilina

<figure><img src="image_2.png" /></figure>

PII: S0048-9697(21)04289-3

DOI: https://doi.org/10.1016/j.scitotenv.2021.149216

Reference: STOTEN 149216

To appear in: *Science of the Total Environment*

Received date: 7 May 2021

Revised date: 16 July 2021

Accepted date: 19 July 2021

Please cite this article as: L. Guillaumot, J. Marçais, C. Vautier, et al., A hillslope-scale aquifer-model to determine past agricultural legacy and future nitrate concentrations in rivers, *Science of the Total Environment* (2018), https://doi.org/10.1016/j.scitotenv.2021.149216

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
A hillslope-scale aquifer-model to determine past
agricultural legacy and future nitrate concentrations
in rivers

Luca Guillaumot¹,²*, Jean Marçais³, Camille Vautier¹, Aurélie Guillou¹,⁴, Virginie Vergnaud¹,
Camille Bouchez¹, Rémi Dupas⁵, Patrick Durand⁵, Jean-Raynald de Dreuzy¹,⁶, Luc Aquilina¹

¹ Univ Rennes, CNRS, Géosciences Rennes, UMR 6118 35000 Rennes, France

² Water Security Research Group, Biodiversity and Natural Resources Program, International Institute for Applied Systems Analysis (IIASA), Luxenburg, Austria

³ INRAE, UR Riverly, F-69625 Villerbanne, France

⁴ Université Savoie Mont Blanc, Polytech-Annecy Chambéry, Le Bourget du Lac 73370, France

⁵ INRAE, Agrocampus Ouest, UMR 1069 SAS, 35000 Rennes, France

⁶ Univ Rennes, CNRS, OSUR (Observatoire des sciences de l'univers de Rennes), UMS 3343, 35000 Rennes, France

<footnote>*corresponding author: guillaumot@iiasa.ac.at</footnote>
# 1 Introduction

Intensive agriculture which has developed since the 1960's has caused eutrophication in aquatic environments (Steffen et al., 2015; Withers et al., 2014). In coastal areas, nitrogen excess led to dramatic green algae proliferations having ecosystem, sanitary and economic repercussions (Galloway et al., 2008; Kronvang et al., 2005; Ménesguen and Salommon, 1988). These problems have raised public awareness and led to regulations to reduce the nutrient load in water bodies (Boers, 1996; European Commission, 1991). As a result, since the 1990's, the agricultural inputs of nitrogen have decreased in many European regions (Abbott et al., 2018; Aquilina et al., 2012; Kronvang et al., 2008; Poisvert et al., 2017). Yet, the impacts of mitigation strategies are still difficult to evaluate and even more to predict (Withers et al., 2014). Especially, the fate of the missing nitrogen (input minus river export), which is either stored or removed, is uncertain and blurred by other uncertainties such as data uncertainties and the imbrications of both spatial and time scales (Breemers et al., 2002).

Nitrate concentrations in rivers do not only depend on the anthropic nitrogen inputs, but also on natural processes that occur within catchments (Chen et al., 2018; Małoszewski and Zuber, 1982) such as the relative importance of overland flows, shallower groundwater flows and deeper groundwater flows or such as autotrophic and heterotrophic denitrification in soils, rivers and aquifers. The rate and dynamics of the excess of nitrogen delivery to the water courses depends on agricultural management and soil processes. Nitrogen excess can be stored in soils and aquifers yielding a temporary retention of nitrogen, the so-called nitrogen legacy (Dupas et al., 2020; Ehrhardt et al., 2019; Hrachowitz et al., 2015; Van Meter et al., 2017, 2016). Indeed, aquifers exert control over the long-term nitrogen concentration in rivers (Aquilina et al., 2012; Hamilton, 2012; Meals et al., 2010) because groundwater can transit several decades before
discharging into the river (Ayraud et al., 2008; Gleeson et al., 2016; Kolbe et al., 2016; Marçais et al., 2018). At a global scale, an important part of river discharge comes from groundwater. In regions with a crystalline geology, such as French Brittany, most river water originates from groundwater systems, with transit times ranging from days to decades. Beside the long-term effect of nitrate transfer with groundwater, the nitrate concentration in rivers is also impacted by a fraction of young water. Indeed, rivers are globally fed by a substantial proportion of water less than three months old (Benettin et al., 2017; Jasechko et al., 2016), which mixes with older groundwater. Young waters either runs off to the river without infiltrating into the aquifer or emerges from saturated shallow horizons (Marçais et al., 2017). This fraction of young water drives a short-term, typically seasonal, variability in nitrate concentrations in rivers (Martin et al., 2004; Molénat et al., 2002; Van Der Velde et al., 2010a). Finally, nitrate concentrations in rivers are also controlled by the microbial denitrification (Knowles, 1982), which has the potential to reduce the total load of nitrate. The denitrification occurs under anoxic conditions after oxygen consumption by aerobic microorganisms, nitrate being the second highest energy level support following oxygen. Such conditions can be found in several compartments of the catchment where oxygen is poorly available such as riparian wetlands, which remain waterlogged during a large part of the year (heterotrophic denitrification with organic matter) and in the deeper part of the saturated zone with longer and deeper flow paths (autotrophic denitrification with pyrite) (Aquilina et al., 2018; Green et al., 2016; Kolbe et al., 2019; Korom, 1992; Molénat et al., 2002; Roques et al., 2018a; Tarits et al., 2006; Van Der Velde et al., 2010b; Van Meter and Basu, 2015). Therefore, evaluating the relative effects of these functional properties of the catchment is crucial to understand the dynamics of nitrate and to predict the impacts of mitigation strategies on river quality.
The role of young water fraction, aquifer storage and denitrification remains uncertain because hydrologists can never fully quantify the groundwater contribution to surface water associated to each point of the transit time distribution (Hrachowitz et al., 2016). Moreover, the competition between transport and denitrification, commonly expressed as the Damköhler number, is known to give a functional, integrative view of the nitrate fate in catchments (Green et al., 2010; Ocampo et al., 2006a; Oldham et al., 2013; Takuya et al., 1993; Zarnetske et al., 2012). Nevertheless, the transport and denitrification properties fundamentally arise from the geological, biogeochemical and hydrological properties of the catchment (Pinay et al., 2015). We argue that dealing with this issues requires (1) data constraining the partitioning of groundwater transit times, (2) to define the most important catchment properties and (3) time scales that they control in the aim to (4) reproduce catchment nitrate dynamics and be able to predict concentration trends.

Many studies have been performed at catchment scales to understand and predict nitrate variations in rivers and groundwater systems. They developed and combined different modelling approaches, inferring different representations of transit times within the catchment (Chen et al., 2018; Hrachowitz et al., 2010). Simple conceptual (or lumped) models can be used (Berghuijs and Kirchner, 2017; Foret et al., 2015; Kirchner et al., 2000; Marçais et al., 2015). They generally represented catchments by one or several reservoirs whose water contents and associated concentrations are governed respectively by linear (or Dupuits equation) and distribution functions of transit times (going from perfect mixing assumption to more complex gamma function coupled with a degradation law). Some studies focused on the relationship between temporal variations and spatial distributions of nitrate (Martin et al., 2006; Ocampo et al., 2006a; Pinault and Pauwels, 2001) mainly at hillslope scale. Other authors implemented 2D
or 3D numerical groundwater models and computed transit time distributions and/or nitrate
concentrations from the resulting flow structure (Gburek and Folmar, 1999; Kaandorp et al.,
2018; Van Der Velde et al., 2010a; Wriedt and Rode, 2006). On the other hand, physically-based
and spatially distributed models such as TNT2-STICS (Beaujouan et al., 2002) and INCA (Wade
et al., 2002), also called mechanistic models, take into account most of the processes impacting
nitrate fluxes such as water transfer in soils, crop development and the associated soil/plant
nitrogen transformation. The accuracy of these models' predictions relies on a precise knowledge
on agricultural practices (in time and space), on soil properties and on the control of epistemic
uncertainties resulting from the simplified conceptualization of the system and on parameter
identification issues. Their hydrogeological part is generally described by simpler conceptual
models.

In this paper, we investigate the controlling factors of the nitrate trends and variability in
three rivers of the highly eutrophic French region of Brittany (Abbott et al., 2018; Aquilina et al.,
2012; Beaujouan et al., 2001; Poisvert et al., 2017). This area has been investigated for almost
two decades on nitrate contamination issues (Ayraud et al., 2008; Conan et al., 2003; Fovet et al.,
2015; Martin et al., 2006; Moléat et al., 2002; Roques et al., 2018b; Ruiz et al., 2002) and for
water resources in crystalline aquifers (Goderniaux et al., 2013; Jimenez-Martinez et al., 2013;
Kolbe et al., 2016; Le Borgne et al., 2006; Leray et al., 2014; Wyns et al., 2004). We focus on
the roles of sub-surface aquifers, which are strongly connected to the rivers. The catchment is
represented by a single equivalent hillslope model governed by groundwater flow equation and
water table interception with the surface. The vertically-resolved approach combined with a first-
order degradation law allows to represent a classic water transit time stratification in the aquifer
and the resulting nitrate stratification. The model is informed by long-term river data such as
streamflows and nitrate concentrations, complemented with punctual CFC-derived ages in rivers, springs and wells. This parsimonious modelling approach enables to test a wide range of values of hydrogeological and biogeochemical parameters (around 20 000 simulations for each catchment), without any assumptions on the aquifer structure. Based on the validated models, we identify the main aquifer properties controlling the nitrate dynamics in terms of hydrogeological behavior, reactivity and geological structure. We further predict future nitrate trajectories from 2015 to 2050 following two input scenarios.

## 2 Study area setting

The study focuses on three agricultural catchments (Douron, Ris and Kerharo) located near the coast of Brittany, Western France. Brittany is the first livestock region in France due to massive industrialization of its agriculture, which began in the 1960's (Deschamps et al., 2016). This led to high inputs of organic and mineral nitrogen (Dupas et al., 2018; Poisvert et al., 2017). The excess of nitrogen causes dramatic green algae proliferations in several bays during the summer season (Ménesguen and Salommon, 1988), implying damages for the coastal ecosystems and for the tourism industry (Gambino, 2014). Since the 1990's, efforts have been made to reduce the agricultural nitrogen inputs, leading to a decrease in the nitrogen concentrations in rivers (Abbott et al., 2018; Aquilina et al., 2012; Poisvert et al., 2017). However, nitrogen concentrations in rivers are still elevated and many rivers have a nitrate concentration exceeding 25 mg/L (Abbott et al., 2018).
<figure><img src="image_3.png" /></figure>

**Figure 1:** The Ris and the Kerharo catchments, located on the westernmost point of France, discharge into the bay of Douarnenez, in the Atlantic Ocean. The Douron catchment, located on the north coast of Brittany, discharges into the bay of Loquairec, in the Manche.

The three agricultural catchments studied here (Figure 1) are discharging into bays subjected to green algae proliferation, making them representative of the regional eutrophication issue. In the three catchments, the major agricultural activity is dairy production. Land use includes maize, winter wheat and rapeseed crops in rotation with ley, as well as pastures. Catchment areas range between 30 and 38 km² (Table 1). The climate of the three catchments is temperate and oceanic, with precipitations relatively distributed over the year (on average 1100 mm/yr). Geologically, the studied catchments are underlain mainly by granite and shales. In Brittany, subsurface is divided into the soil layer (a few decimeters), the weathered zone (a few
meters to a few tens of meters), the fractured zone (a few tens of meters to a few hundreds of meters) and the fresh basement (Maréchal et al., 2004; Molénat et al., 2013; Mougin et al., 2008; Wyns et al., 2004). The Ris catchment, as delimited by the gauging station, (northern blue point on **Figure 1**), is mainly composed of granite. The geological map reveals regional fractures aligned to a long a NW-SE axis and some shale bodies in the downstream section (provided by the BRGM, French Geological Survey). The Kerharo catchment is composed of NW-SE fractures and metamorphous shales (micaschist) with pyrite found during borholes drilling at the bedrock interface (few meters depth) (Faillat et al., 1999). The Denton catchment is composed of different granitic bodies and some metamorphous shales in the downstream section.

This subsurface structure strongly connects surface and sub-surface water flows, with many interactions between aquifers, soils and rivers (Mougin et al., 2006; Molénat et al., 2013; Ruiz et al., 2002). Estimated recharge to aquifers is on average close to 400 mm yr⁻¹ (Habets et al., 2008; Le Moigne, 2012; Quintana-Seguí et al., 2008) making the weathered horizon the main capacitive layer for water (Ayraud et al., 2008; Wyns et al., 2004). For the three catchments, the hydrogeological conceptual model is composed of a weathered capacitive layer (< 2 m depth) overlying a fractured draining layer. Surface runoff mainly occurs as saturation excess overland flow in valley bottoms where the aquifer intersects the land surface during the wet season (Ogden and Watts, 2000). Excess infiltration overland flow is generally low in Brittany as precipitations are equally distributed throughout the year. Streamwater therefore mainly consists in surface runoff, lateral transfer within the soil, and groundwater circulation in the weathered parts of the aquifer. The proportion of river water originating from the aquifer varies throughout the year, but is globally high. Previous studies in Brittany found a groundwater contribution to rivers of 55% on average (Mougin et al., 2006). By comparing estimated mean groundwater
recharge from land surface model applied over France (SURFEX platform) (Habets et al., 2008; Masson et al., 2013; Quintana-Seguí et al., 2008) with observed mean streamflow, we found a groundwater contribution to stream of 80% on average for the three studied catchments.

**Table 1:** Main characteristics of the Ris, Kerharo and Douron catchments. Geology is given by maps of the French Geological Survey (BRGM). Climate data (average from 1998-2018) come from the SURFEX platform. Discharge data and nitrate concentrations (average from 1998 to 2010 and from 1998 to 2015 respectively) are given by the water basin agencies. The trend in nitrate concentration in the river is obtained by a linear regression.

<table><thead><tr><th>Property</th><th>Ris</th><th>Kerharo</th><th>Douron</th></tr></thead><tbody><tr><td>Area [km²]</td><td>31</td><td>38</td><td>30</td></tr><tr><td>Slope [%]</td><td>5</td><td>6</td><td>6</td></tr><tr><td>Characteristic length [km]</td><td>1</td><td>1.1</td><td>1.1</td></tr><tr><td>Dominant lithology</td><td>Granite<br>(+fault zone)</td><td>Micaschist</td><td>Granite</td></tr><tr><td>Precipitation [mm/yr]</td><td>1132</td><td>1072</td><td>1145</td></tr><tr><td>PET [mm/yr]</td><td>670</td><td>655</td><td>610</td></tr><tr><td>Discharge [mm/yr]</td><td>569</td><td>473</td><td>530</td></tr><tr><td>Mean nitrate concentration<br>in river [mg/L]</td><td>34</td><td>32</td><td>31</td></tr><tr><td>Trend in nitrate<br>concentration [mg/L/yr]</td><td>-0.46</td><td>-0.92</td><td>-0.55</td></tr></tbody></table>
# 3 Material and methods

## 3.1 Field data

### 3.1.1 Long term monitoring data: discharge and nitrate streamwater time series

To inform the main hydro-biogeochemical processes occurring at the catchment scale, long time series of stream discharge and nitrate concentrations observed in the three catchment rivers were gathered. First, monthly river discharge from 1998 to 2010 are presented (Figure 2). These streamflow data were initially collected by measuring the daily river height (http://hydro.eaufrance.fr/). Average observed streamflow values are given on Table 1 for the 1998-2010 period. Annual discharge standard deviation, representative of the interannual streamflow variability, is on average equal to 0.17 m/yr for the three catchments to compare with the average discharge of 0.52 m/yr. Seasonal variations are quite regular through years and influenced mainly by the low evapotranspiration demand between November and March. Streamflow are approximately 4 times higher during the wet season compared to the dry season indicating a pronounced seasonality and a fast answer to recharge events (Figure 2). The strongest seasonal variability is found for the Kerharo catchment with streamflow ranging each year between 4 and 100 cm/month roughly. In spite of the distance between Ris and Douron catchments (80 km), they show very similar seasonal behaviors. The smaller response time for the Kerharo catchment, illustrated by recessions slope on Figure 2, is attributed to its geomorphological and geological features (Table 1) as this catchment is very close from the Ris catchment, with quite similar climate and soil occupation.
<figure><img src="image_4.png" /></figure>

**Figure 2:** Observed specific streamflow normalized by catchments area (in logarithmic scale) for the three studied catchments along the period of calibration of the model. Monthly data are derived from daily time step streamflow. Grey filled areas correspond to estimated groundwater recharge rates for the Kerharo catchment (section 2.4.1).

Secondly, monthly nitrate concentrations in rivers from 1998 to 2015 are presented (Figure 3). These data were initially collected at bi-monthly time step by local basin agencies. Average nitrate concentration in rivers is around 32 mg/L for the three catchments over the 1998-2015 period (Table 1). For all catchments, concentrations are decreasing during this period (on average -0.6 mg/L/yr), with highest decrease for the Kerharo river (-0.92 mg/L/yr). In addition, nitrate concentrations were occasionally measured between 1976 and 1984 in the Ris river (black dots on Figure 3). These data come from NAÏADES French data base (http://www.naiades.eaufrance.fr/). This data was added to the 1998-2015 nitrate monitoring datasets as they provide an interesting opportunity to assess the long-term behavior of the catchment. As highlighted by the interpolated dashed curve on Figure 3, based on several trends
observed in Brittany (Aquilina et al., 2012; Dupas et al., 2020, 2018; Kolbe, 2017), measured
data did not cover the nitrate peak concentrations occurring between 1990 and 2000.

<figure><img src="image_5.png" /></figure>

**Figure 3:** Observed nitrate concentrations in the river at the discharge stations of the three studied catchments. Initial data are at weekly time step, unless on the Ris river for the period 1976-1984 where punctual data were available (blue acts). Shaded areas corresponds to estimated recharge periods. Dashed curve extrapolates periods without data, with a maximum of concentration occurring between 1990 and 2000 in Brittany.

**3.1.2 Sampling campaigns: CFCs groundwater age tracers**

CFCs and nitrate concentrations were measured in boreholes, springs and rivers in the three studied catchments. Two sampling campaigns were performed, at the end of the wet season (March 2019) and at the end of the dry season (October 2019). Sampling locations are shown in **Figure 1**. For each catchment, we performed three samplings in the river and 7 to 13 samplings in springs and boreholes. The number and the location of springs and boreholes samplings were
constrained by accessibility and landlord authorization. Details about the sampling techniques
used to measure CFCs are presented in the Supplementary Material.

CFCs concentrations were then converted into apparent ages through the use of a Lumped
Parameter Model (Marçais et al., 2015) to be used for calibrating the solute transport parameters
of the hillslope model (section 3.2.2). To do so, we assumed an a priori parametric age
distribution that we convoluted with the input atmospheric CFCs concentration. Exponential
distribution of transit times is compatible with 1D Boussinesq aquifer model and relevant while
seasonal fluctuations do not impact significantly the mean transit time and while the aquifer
saturation is limited (low young water contribution). Thus, exponential distribution is appropriate
to estimate mean water transit time from CFCs sampling in boreholes, spring and rivers
considering all these points are located in convergence zones (Marçais et al., 2015), downstream
of the hillslopes or subject to farmer pumping. Note that CFCs in deep wells (>30m depth)
consistently show very low concentrations revealing long transit times which can be modelled by
gamma distribution. Indeed, when groundwater is old (typically >60 yr considering the
exponential model), the accuracy of the CFC tracers decreases. Thus, these older values will be
considered as a lower boundary and will not be used in the model calibration. For each
catchment, the observed mean transit time, later compared to the simulated mean transit time
(τ<sub>TT</sub>), is an average between March and October CFC-derived ages of all springs and shallow
boreholes. To prevent misinterpretation due to potential CFCs contamination, we only retain
samplings where at least two CFCs among CFC-11, CFC-12, CFC-13 lead to similar apparent
ages (Jurgens et al., 2012).
## 3.2 Modelling approach

The model simulates flow and nitrate transport in three aquifers over the 1955-2010 period. For the three studied catchments, hillslope models were calibrated over the 12-year period from 1998 to 2010 for streamflow and over the 10-year period from 2000 to 2010 for nitrate concentrations in river (section 3.2.3). For Ris and Kerharo catchments, the model is extrapolated until 2015 (see Supplementary Material). Here, we present (1) recharge and nitrate inputs coming from a hydrological and nitrogen soil/plant model (2) the hillslope aquifer model (3) the calibration method and (4) the relevant model properties.

### 3.2.1 Recharge and nitrogen inputs

To estimate recharge and nitrate inputs over the 1955-2015 period, we merged two modeled datasets providing either nitrate and/or groundwater recharge inputs over different time scales. First, we used the TNT2-STICS model providing monthly recharge rates and nitrate concentrations in this recharge water over the 1997-2010 period (Beaujouan et al., 2002). This period is constrained by nitrate measures in river and the monitoring of agricultural practices. TNT2-STICS is a spatially distributed hydrological and nitrogen model running at daily time step and taking into account several soil layers for water transfers. Crop development and the associated soil/plant nitrogen transformation are physically represented along with heterotrophic denitrification in soils and in often saturated areas like humid zones. TNT2-STICS is alimented by climate data (rainfall and potential evapotranspiration), nitrogen inputs by agriculture and crop specific parameters. To model nitrate legacy of past practices, the hillslope model is run on a longer period starting some more than 40 years before, in 1955. We combined TNT2-STICS nitrate inputs with nitrate annual surplus estimates covering the 1955-2015 period (Poisvert et
al., 2017) to retrieve monthly estimates of nitrate inputs over the 1955-1997 period and over the 2010-2015 period, except for the Douron catchment where lack of input data prevent to extend monthly nitrate inputs over the 2010-2015 period (see Supplementary Material). Thus, groundwater recharge from TNT2 constitutes the only water input to the hillslope model. Note that TNT2 already takes into account soil-groundwater exchanges using a simplified representation of the groundwater compartment. We used here the net groundwater recharge (recharge minus “excess groundwater”) corresponding to more than 23 % of the river streamflow for the studied catchments (Beaujouan et al., 2002). Excess infiltration overland flow is considered negligible (section **Error! Reference source not found.**). Finally, the hillslope-scale aquifer model restitutes this recharge to the river as below, aquifer seepage and saturation excess overland flow (section 3.2.2).

**Figure 4** presents the simulated nitrate inputs to the aquifer for the three catchments from 1955 to 2015. They are given in mass per year by multiplying the monthly groundwater recharge rates by the concentrations of nitrate in the recharge. They show a strong increase from 1955 to 1990, followed by a decrease between 2000-2005 after the nitrate directive (1991). Since 2005, the trend is masked by an inter-annual variability stemming from several factors including climate and a lower decrease in agricultural inputs. The three catchments behave slightly differently as illustrated by the variability between the three curves on **Figure 4**. Differences come from land use and agricultural practices.
<figure><img src="image_6.png" /></figure>

**Figure 4:** Annual nitrate mass inputs in the three studied catchments. Values are obtained from TNT2-STICS model and are later imposed to the groundwater model. Note that inputs are given to the model at a monthly time step but are presented here at annual scale and smoothed using a moving average of three years for better visualization.

### 3.2.2 The 2D hillslope model

A two-dimensional hillslope model is defined for each catchment (Figure 5). It assumes that, at the catchment scale, the hydrological system behavior is controlled by two hillslope-scale processes, which are the groundwater flow stratification processes and the saturation-excess runoff processes coming from the water table interception with the surface (Brutsaert, 1994; Fan et al., 2019; Marçais et al., 2017; Matonse and Kroll, 2009; Troch et al., 2003). The structure of the hillslope (length and slope) is derived from averaged geomorphologic properties (Table 1). The bottom of the hillslope is considered horizontal. Its characteristic thickness (E) is defined at the most downstream point of the hillslope (Figure 5). Its hydraulic conductivity (K) and porosity (θ) are assumed uniform. The bottom, left and right boundaries are no-flow conditions (Figure 5).
The right boundary (upstream) corresponds to the topographic divide of the hillslope while the left boundary (downstream) corresponds to the river flowing perpendicularly to the modeled section. When the water table reaches the surface on the upper limit (blue line on Figure 5), groundwater seeps and directly feeds the river. With higher recharge, the seepage front moves upstream along the horizontal axis. Streamflow in the river results from the discharge of the aquifer to the river and upstream from the seepage area. Given the limited extension of the hillslope and the weekly time step of the model, seepage is assumed to be transferred without any delay to the river. Similarly, streamflows are compared with immediate river routing scheme at the outlet of the catchment given the small catchment size (~30 km²). Flow and nitrate transport models are implemented within the well-known MODFLOW and MT3D software suite through the Python Flopy interface (Bakker et al., 2016; Bedekar et al., 2016; Harbaugh, 2005; McDonald and Harbaugh, 1984; Zheng and Wang, 1999). They are presented in details in Supplementary Material. The modelled mean transit times are computed by dividing the groundwater volume by the mean recharge. They are compared to the mean transit time derived from the CFCs (section 3.1.2). Groundwater nitrate inputs, as well as recharge rates, are provided by the TNT2-STICS model. Thus, we only consider here nitrate transport in the aquifer, without taking into account the transport through the soil done by TNT2-STICS. Nitrate is denitrified in the aquifer under favorable anoxic redox conditions with accessible electron donors like pyrite (Green et al., 2016; Kolbe et al., 2019; Korom, 1992; Molénat et al., 2002; Tarits et al., 2006; Van Der Velde et al., 2010a; Van Meter and Basu, 2015). Denitrification inside the aquifer is modeled by an effective first-order reaction (equation (1)):

$$r(\tau) = 1 - \exp\left(-\frac{\tau}{\tau_{NO3}}\right) \quad (1)$$
where $\tau_{NO3}$ is the characteristic denitrification time and $\tau$ is the transit time. Most of the heterotrophic denitrification in often saturated areas such as humid zones is already implicitly taken into account in the TNT2 soil model.

<figure><img src="image_7.png" /></figure>

**Figure 5:** Sketch of the 2D hillslope model. The morphology of the hillslope is defined by its length (L) (here 1000 m), its characteristic thickness (E) (here 80 m) and the slope of its surface (here 0.05). The model is parameterized by its porosity ($\theta$) (here 0.05), its characteristic thickness (E) (here 80 m) and its hydraulic conductivity ($K$) (here 4.10-6 m/s). Transit times are represented in colors on the flow lines (mean transit time is 9.6 yr). The stratification of times is classical for homogeneous hillslope models. Streamflow in the river results from the discharge of the aquifer to the river (upper left corner of the model) and upstream from it as seepage. Given the limited extension of the hillslope, seepage is assumed to be transferred without any delay to the river.

The hillslope model has four unknown parameters, which are its hydraulic conductivity ($K$) [m/s], its thickness (E) [m], its porosity ($\theta$) [without unit] and its characteristic denitrification time ($\tau_{NO3}$) [yr]. From these parameters, we also use a closely related indicator, the mean
saturated thickness ($E_{wet}$) (obtained after simulation). The transmissivity ($KE_{wet}$) partitions the streamflow between groundwater flows and saturation excess overland flows through the mean saturated thickness and through the mean extent of the seepage zone. The time fluctuations of the overland flows and saturated thickness are conditioned by the porosity ($\theta$). The smaller is the porosity, the larger the fluctuations of the seepage zone among seasons. The porosity also intervenes in the mean equivalent water height ($\theta E_{wet}$) which, divided by the recharge, gives the mean transit time ($\tau_{TT}$) (Cornaton and Perrochet, 2007; Danckwerts, 1953; Haitjema, 1995) (equation (2)):

$$ \tau_{TT} = \frac{\theta E_{wet}}{R} \qquad (2) $$

where $R$ is the recharge expressed in [L/T]. Thus, $\theta E_{wet}$ represents the groundwater volume normalized by the hillslope surface.

The characteristic denitrification time ($\tau_{NO_3}$) influences nitrate concentrations but also will adjust itself to the aquifer volume immediately through the calibration. It therefore controls the relative importance of nitrate dilution and removal. All four parameters ($K, E, \theta, \tau_{NO_3}$) are thus interdependent when calibrated from the observed flows, CFCs-derived ages and nitrate concentrations. In-stream denitrification is tightly related to residence time in rivers and diffusion processes within the hyporheic zone (Boulton et al., 1998; Gabriel et al., 2006; Zarnetske et al., 2011). In-stream heterotrophic denitrification is considered negligible within this study because of the small residence time in the short rivers of low Strahler orders (Lefebvre et al., 2007; Montreuil et al., 2010) and the high river oxygenation (Vautier et al., 2020). Low in-stream denitrification might only appear during dry season and would correspond to a limited period of the year when the streamflow is very low (<10% of the annual flow). Nitrate can also be consumed by macrophyte and phytoplankton but this usually occurs in streams with higher
Strahler order where their development is favored (Durand et al., 2011). As a whole, such denitrification processes should remain low and represent a limited extent of nitrate fluxes without noticeable modifications of the results presented here.

Finally, this model is simple and parsimonious in order to be calibrated with a limited amount of data. It is a mechanistic model proposing temporally and spatially dynamic representations of saturation, subsurface flows and surface interception essential to reproduce the conditions of relatively shallow aquifers (10-100 meters) and temperate climates as it prevails in Brittany. Indeed, the wet winter season in Brittany favors high recharge rates and the development of seepage in a landscape of small slopes and limited aquifer capacities (Goderniaux et al., 2013; Kolbe et al., 2016; Merot et al., 2014). The three third-order catchments studied here are in fact made up of small hillslopes of average length around 1 km as determined by the transition from hillslope to river in a classical area-slope relation (Lague et al., 2000) (Table 1). They are relatively flat with mean surface slopes of 5 to 6%.

### 3.2.3 Calibration method

The objective function of the optimization problem was obtained as a combination of classic calibration functions for the three types of observables: streamflow time series, nitrate time series and CFC-derived mean transit times. We used the Nash-log criterion to assess the ability of the model to reproduce monthly streamflow time series (1998-2010, Figure 2). Nash-log equilibrates the relative influence of the low and high flows and is particularly adapted to focus on groundwater flows (Gupta et al., 2009). The monthly nitrate concentrations in rivers (2000-2010, Figure 3) are compared to the simulated data using the root mean square error (Gupta et al., 2009) criterion normalized by the standard deviation of observations (nRMSE). To focus the calibration on the long-term trend, both observed and simulated monthly-river
concentrations are smoothed out by keeping only the best linear function according to a least-squares fit on the data. Punctual nitrate concentrations measured between 1976 and 1984 in the Ris river (blue points on Figure 3) are used a posteriori to assess the relevance of the selected models. Finally, the difference between modelled and observed mean transit time ($\tau_{TT}$ and $\tau_{TT,obs}$) is normalized by a reference age $\tau_{TT,ref}$ of 15 yr. For each catchment, the mean observed transit time ($\tau_{TT,obs}$) corresponds to the average between March and October of CFC-derived ages of all springs and shallow boreholes (12 values).

Finally, the three calibration targets are combined on the basis of a comparable range of variations between 0 and 1, 1 being a perfect match. To this end, the normalized differences $x$ was transformed in $e^{-x^2/2}$ (Vrugt and Sadegh, 2013). The three criteria are combined by the following minimum function to guarantee a minimum adequation for each of the target (equation (3)):

$$J = \min \left( Nash_{to} ; e^{-\frac{nRMSE^2}{2}} ; e^{-\frac{1}{2}\left(\frac{\tau_{TT,obs}-\tau_{TT}}{\tau_{TT,ref}}\right)^2} \right) \quad (3)$$

The threshold $J > 0.7$, considered here as a good fit for the models, was chosen to select the optimal models. For such a value, the error variance of the simulated streamflow, in logarithmic scale, is limited to 30% of the variance of the observed time series. For such a value, root mean square error of the simulated nitrate concentration is smaller than 85% of the standard deviation of the observed signal (thus, smaller than 2.4, 3.8 and 2.5 mg/L for Ris, Kerharo and Douron respectively). For such a value, models with a mean transit time 12 yr higher or lower than the observed one are excluded.

Because of the strong interdependence of the four parameters ($K$, $E$, $\theta$, $\tau_{NO3}$), we chose to calibrate them simultaneously with a systematic sampling of the parameter space, meaning that
parameters values are sampled regularly inside a range of plausible values. The explored
parameters range over broad intervals extending the typical values reported in comparable
geological settings of Brittany and beyond (Table 2). Hydraulic conductivities were sampled in
the range 2 10⁻⁷–10⁻⁴ m/s wider than the range of 10⁻⁶–2 10⁻⁵ m/s derived from previous regional
studies in shallow aquifers of Brittany (Clément et al., 2003; Grimaldi et al., 2009; Kolbe et al.,
2016; Le Borgne et al., 2006; Legchenko et al., 2004; Martin et al., 2006; Roques et al., 2014).
Porosities as high as 50% have been found in shallow weathered zones (Kovacs, 1981; Wright
and Burgess, 1992) while granites and schists have lower porosities (0.1–1%) (Earle, 2015;
Hiscock, 2009; Singhal and Gupta, 2010). Averaged over the whole weathered and fissured
horizons, we considered a 1–20% range for the porosity values. A range of 50 to 300 m for the
aquifer thickness was derived from a regional synthesis for the weathered and fissured zones
built based on well logs (Mougin et al., 2008, 2006). Thus, transmissivity values from 1 10⁻⁵ to
3 10⁻² m²/s were tested. A 1–150 yr range for the characteristic denitrification time covers all
possible values as values lower than 1 yr correspond to almost instantaneous denitrification. A
case without denitrification was also included. Each parameter interval was regularly discretized
(in a logarithmic scale for hydraulic conductivity and denitrification time) with the number of
values given by the right column of Table 2. The number of values was adapted to the width of
the interval. Resulting from the parameter combinations, 19 200 simulations were run for each
catchment. For each simulation the criterion *J* was computed (equation (3)).
<table><caption>**Table 2:** List of the calibrated parameters (left column) with the range of values explored (middle column) and the number of values explored (right column).</caption><thead><tr><th>Parameter</th><th>Range</th><th>Number of values</th></tr></thead><tbody><tr><td>Hydraulic conductivity K [m/s]</td><td>2 10<sup>-7</sup>-1 10<sup>-4</sup></td><td>30</td></tr><tr><td>Porosity θ []</td><td>0.01-0.2</td><td>16</td></tr><tr><td>Thickness E [m]</td><td>50-300</td><td>5</td></tr><tr><td>Denitrification time τ<sub>NO3</sub> [yr]</td><td>1-100 + no denitrification (∞)</td><td>8</td></tr></tbody></table>

### 3.2.4 Model properties

Several properties were computed on the calibrated models to characterize their hydrological and geochemical behaviour. The hydrological behavior was characterized using the characteristic hydraulic response time, the percentage of young water in the river and the relative extent of the seepage zone. The geochemical behavior was characterized using the mean transit time, the Damköhler number and the characteristic denitrification depth. We underline that we use the Damköhler number as a simple indicator to compare chemical and physical processes and not as a full interpretation framework of denitrification processes like what has been proposed in riparian and hyporheic zones (Gu et al., 2007; Harvey et al., 2013; Ocampo et al., 2006a).

The hydraulic response time $τ_H$ is the characteristic delay between a recharge event and the increase of flow in the river. It is defined by **equation (4)** (Gelhar, 1974; Molénat et al., 1999; Townley, 1995):

$$τ_H = \frac{θL^2}{KE_{wet}} \quad (4)$$
where $E_{wet}$ is the mean saturated thickness. The young water percentage is here defined as the contribution to the river of waters infiltrated during the same year. This percentage is dominantly controlled by surface processes as well as shallow groundwater transfers with short transit times if hydraulic conductivity is high. It is obtained by a particle tracking approach (Supplementary Material). The relative extent of seepage can be calculated with the surface area where overflows occur (Figure 5). The overall existence of seepage in Brittany during the winter season is confirmed by field observations (Franks et al., 1998; Merot et al., 2003). Seepage occurs mostly in winter because of the proximity of the water table to the surface and because of the recharge period concentrated on a restricted time range from November to March. The percentage of seepage areas should be smaller than 31% for the three catchments as estimated from geomorphologic and climate data (MEDDE and CUS Sol, 2014). While not a direct target of the calibration process, the existence and relative importance of seepage will be used to confirm the consistency of the calibrated model.

The geochemical behavior of the model is characterized by the mean transit time and the Damköhler number. The mean transit time $\tau_{TT}$ is the mean travel time of an element from its inlet in the catchment to its outlet in the river defined by equation (3). It is straightforwardly given by the ratio of the aquifer volume to the overall recharge. The Damköhler number $Da$ is defined as the ratio of the characteristic denitrification time and transit time (Ocampo et al., 2006b) (equation (5)):

$$Da = \frac{\tau_{TT}}{\tau_{NO3}} \qquad (5)$$

When the $Da$ is smaller than one, the process is reaction limited, while when the $Da$ is larger than one, the process is reaction limited. The denitrification time $\tau_{NO3}$ can be related to a characteristic depth of denitrification in the aquifer using the classical stratification law
illustrated on **Figure 5** according to which the transit time at a position *x* and at a depth *z* is equal to (equation (6)):

$$ \tau = \tau_{TT} \times \ln \left( \frac{E(x)}{(E(x) - z)} \right) \quad (6) $$

where *E(x)* is the saturated thickness at *x* (Vogel, 1967). Assuming that *E(x)* can be approximated by the mean saturated thickness *E*<sub>wet</sub>, we can combine equation (6), (5) and (1) to derive the characteristic denitrification depth *z*<sub>NO3</sub> normalized l v the mean saturated thickness such as (equation (7)):

$$ \frac{z_{NO3}}{E_{wet}} = 1 - (1 - r)^{\frac{1}{D_c}} \quad (7) $$

where *r* is the progress of the reaction ranging between 0 and 1 respectively for a zero denitrification and a complete denitrification. Inally, we will consider in this study the normalized characteristic denitrification depth at *r* = 25%. While it does not mean that denitrification occurs specifically at this depth, it traduces the characteristic denitrification time in terms of characteristic depth, which could ideally be related to lithological changes (Kolbe et al., 2019).

# 4 Results

## 4.1 Catchment hydrological functioning deduced from field data

The preliminary analysis of the field data reveals the existence of rapid transfers from precipitation to river flows. Streamflow increases quickly after the dry season showing a fast response typical of the rapid saturation of lowland areas close to the river (Figure 2). Streamflow decreases slowly after the main recharge period occurring around the winter season, a typical
response of groundwater flow systems. Interpreted as an exponential decrease of streamflow
with time, characteristic response times range between some weeks to a few months. The
exponential decrease typically comes from the drainage of a reservoir where streamflow and its
derivative are linearly proportional (Brutsaert and Nieber, 1977).

Chemical samplings in wells show two well marked components in groundwater (Figure 6). The first component made of “shallow subsurface waters” sampled mostly in shallow wells (<20m) is characterized by high nitrate concentrations (~25-50 mg/L) and transit times of 20 to 25 years as determined by the concentrations of CFCs. It already results from the contribution of different groundwater flowpaths (Supplementary Material). The second component made of “deep water” sampled in deeper wells (from 30 to > 100m) has almost no nitrate and no CFCs revealing transit times higher than 40 years. This value is obtained considering a gamma transit time distribution (using a shape parameter of 0.5) as the exponential distribution gives values between 60 years and hundreds of years (section 3.1.2).

In the rivers and springs, nitrate concentrations (10-30 mg/L) are intermediate between
these two components (Figure 6). Transit times estimated from CFCs sampled in rivers establish
mostly around 20 years and around 20 to 25 years for the springs (Figure 6), similar to shallow
wells. The springs and river concentrations and transit times result from the mixing of surface
water i.e. water from saturated areas and from the top of the watershed during high water levels,
and of the two groundwater components identified above: shallow subsurface (shallow wells)
and deep subsurface discharges (deep wells). Regarding the original data sampled in rivers, the
deduced transit time should be taken as a lower bound because CFC concentrations measured in
the rivers are lower than the current atmospheric concentration, thus exchanges with the
atmosphere can only increase river CFC-concentrations and biased estimations towards younger
transit times. Finally, the fast surface component defined as “young water”, including saturation
excess overland flow and < 1yr subsurface waters, was probably weakly sampled during our
campaign because this component occurs in the most superficial part in downstream areas and
during punctual precipitation events.
<figure><img src="image_8.png" /><figcaption>**Figure 6:** Characteristic groundwater transit times derived from CFC concentrations with an exponential transit time distribution (or a gamma distribution for older and deeper boreholes) and nitrate concentrations from two sampling campaigns performed in March 2019 and October 2019. The two graphs on the bottom represent the average over the three catchments displayed above.</figcaption></figure>
Estimated mean transit times (~20 yr) from rivers, springs and shallow wells are in line with the long-term nitrate data analysis performed by Dupas et al. (2018) who found a delay of 10 yr between the downward trend of inputs and outputs across Brittany. More recently, Dupas et al. (2020) obtained median transit times ranging from 4 to 16 yr, while Martin et al. (2006) suggested a mean groundwater residence time of 14 yr at hillslope scale.

The three end-members (surface young water, shallow subsurface and deep subsurface components) have different effects on the nitrate concentration of the receiving water bodies (springs and rivers). The two groundwater components delay the response of the river to the applied nitrate signals. In the three rivers considered, the decrease of nitrate input in the 1990-2000 period (Figure 4) is shifted in the river to the 2000-2010 period (Figure 3). As a result, river concentrations still decrease in the 2000-2010 period while nitrate input concentrations have stabilized. At smaller seasonal and inter-annual scales (1-5 years), nitrate concentrations strongly vary as a result of the large variations of the surface flow component, of the nitrate inputs between seasons and between successive years, and soil processes variations. Simple correlations are however limited because of the superposition of flow processes operating on different time scales as well as because of their convolution with complex input flow and nitrate forcing conditions.

This preliminary analysis shows the importance of the groundwater transit time stratification and of the aquifer saturation conditions on river nitrate concentrations. Transit time stratification fundamentally comes from the decrease of flow with depth (Bresciani et al., 2014; Vogel, 1967). It is often enhanced by the reduction of hydraulic conductivity with depth but does not require it. Aquifer saturation essentially controls the relative contribution of surface flows through the presence and extent of the seepage area. These two essential processes are the basis
of the hillslope model presented in section 3.2.4 and developed to further quantify the control of
flow organizations on downstream nitrate concentrations.

**4.2 Properties of the calibrated models**

The calibration method described section 3.2.3 has been successively applied to the three catchments. The selection criterion *J* > 0.7 has been met in the three cases for 6 to 36 models of the 19 200 tested models. The objective criterion *J* is on average equal to 0.73, 0.83 and 0.75 for the selected models of Ris, Kerharo and Douron catchment respectively. It goes up to 0.78, 0.91 and 0.78 respectively. The hillslope modeling approach appears to be relevant to capture the main features of the available observations on all three catchments. Even though it is parsimonious, the hillslope model has essential capacities for simulating the previously discussed dynamics and their effects on the transport of inert and reactive transport. The range of values of the calibrated parameters given in *Table 3* along with their properties defined in section 3.2.4 shows the key features of the equivalent water height, hydraulic conductivity and denitrification capacity.
<table><caption>**Table 3:** Range of calibrated values for the four parameters (K, E, θ, τ<sub>NO3</sub>) with the associated range of model properties defined by section 2.2.4 (grey background).</caption><thead><tr><th>Parameters</th><th>Ris</th><th>Kerharo</th><th>Douron</th></tr></thead><tbody><tr><td>K [m<sup>2</sup>/s]</td><td>8 10<sup>-6</sup>–8 10<sup>-5</sup></td><td>1.8 10<sup>-5</sup>–1 10<sup>-4</sup></td><td>6 10<sup>-6</sup>–2.2 10<sup>-5</sup></td></tr><tr><td>θ []</td><td>0.025–0.1</td><td>0.03–0.12</td><td>0.025–0.05</td></tr><tr><td>E [m]</td><td>100–300</td><td>100–100</td><td>100–300</td></tr><tr><td>τ<sub>NO3</sub> [yr]</td><td>100–∞</td><td>50–100</td><td>50–100</td></tr><tr><td>τ<sub>H</sub> [yr]</td><td>0.24–0.65</td><td>0.11–0.47</td><td>0.32–0.78</td></tr><tr><td>Young Water [%]</td><td>7–14</td><td>6–12</td><td>12–15</td></tr><tr><td>Proportion of seepage zone [%]</td><td>2–7</td><td>1–4</td><td>3–8</td></tr><tr><td>θE<sub>wet</sub> [m]</td><td>1.4–12.3</td><td>9.2–14</td><td>5.3–7.9</td></tr><tr><td>τ<sub>TT</sub> [yr]</td><td>15–25</td><td>21–32</td><td>10–15</td></tr><tr><td>Da []</td><td>0.12–0.17<br>0 if no denit.</td><td>0.32–0.55</td><td>0.15–0.2</td></tr><tr><td>25%-denitrification depth<br>[equivalent depth required to<br>achieve denitrification of 0.5%],<br>expressed as a % of total<br>equivalent water heat at</td><td>82–91%<br>(7–9m of eq.<br>water height)</td><td>41–59%<br>(5–6m of eq.<br>water height)</td><td>76–85%<br>(5–6m of eq.<br>water height)</td></tr></tbody></table>

The porosity θ and hillslope thickness E are in the middle and upper ranges of the explored values (Table 2). They display significant uncertainties with factors of variations from 2 to 4 between the highest and lowest interval bounds. Both intervene in the global groundwater volume scaling with their product. The equivalent water height θE<sub>wet</sub> and the mean transit time (equation (2)) are much better defined with uncertainties limited to a factor of variations of at most 1.7 because parameters compensate each other. The mean transit time is efficiently
constrained by the groundwater age tracing of the CFCs and by the observed long-term trend of nitrate concentration in the river (delay between the applied and discharged nitrate, mean concentration and long-term trend between 2000 and 2010).

The calibrated values of the hydraulic conductivity *K* (*Table 3*) are in the upper range of the explored interval (*Table 2*). Their uncertainty is higher than that of *θ* and *E* with a factor of variation between 3.5 and 10. It might however be interpreted as a limited uncertainty given that hydraulic conductivities can vary over several orders of magnitudes (nearly three orders of magnitude in the explored values). Hydraulic conductivity should be small enough to allow the occurrence of seepage and high enough to reproduce the dynamics of streamflow recession. Indeed, hillslopes should not be too conductive for the seepage area to develop in the winter season. Even if the seepage area is restricted to 1-3% on average, the percentage of young water in the river can go up to 6-15% as recharge is positively correlated with the occurrence of seepage (*Table 3*). To the opposite, hillslopes should be conductive enough to capture the relatively quick recessions of characteristic times τ<sub>H</sub> between 0.11 yr and 0.78 yr. The relevance of the hillslope model is confirmed by its possibility to fulfill both constraints with a restricted uncertainty.

The denitrification time $\tau_{NO_3}$ (*Table 3*) is in the upper range of the explored interval (*Table 2*). Values mostly between 50 and 100 years show that denitrification is limited but not negligible. Only one of the qualified solutions for the Ris catchment does not require any denitrification. The occurrence of denitrification is confirmed by the values of the Damköhler number smaller than 1: 0.1-0.2 for Ris and Douron and 0.3-0.5 for Kerharo. The Damköhler number is well constrained showing that denitrification is separated from dilution by the simultaneous analysis of CFCs, which can be interpreted as reference conservative tracers. The
occurrence of denitrification is eventually confirmed by the values of the produced dinitrogen
gas ($N_2$) measured in the samples (Supplementary Material). Indeed, the increase of $N_2$ appears
inversely correlated to the nitrate concentration highlighting the transformation of $NO_3^-$ into $N_2$.
Moreover, high sulfate concentrations are observed in deep wells associated with low nitrate
concentrations (Supplementary Material) suggesting pyrite oxidation and nitrate reduction
(Green et al., 2016; Korom, 1992). Associated with the natural stratification of times in the
aquifer, the small Damköhler numbers indicate that denitrification occurs mostly in the lower
part of the aquifer. The characteristic 25%-denitrification depth ranges between 41% and 91% of
the equivalent water height. Thus, only around 25% of the denitrification is achieved in the first
5 to 9 meters of the equivalent water height.

While globally similar, the three catchments show slight differences. The Kerharo catchment has both higher hydraulic conductivity and porosity than the two other Ris and Douron catchments. It is consistently reduced by smaller hydraulic characteristic times and higher transit times. In the same time Damköhler number indicates a more efficient denitrification. From Figure 3 and 4, one can see that both river concentration and inputs for Kerharo reduced faster from 2000 to 2010 making difficult to interpret the different properties obtained. The higher transit time in Kerharo should lead to a smoother behavior but the decrease from 2000 to 2010 is accentuated by the more efficient denitrification.

As an intermediary summary, the combination of river flows, river nitrate concentrations and CFC concentrations is relevant to model the dominant factors controlling the overall flow, conservative transport and first-order nitrate reactive transport. Catchments have well defined hydraulic parameters. Denitrification is well separated from dilution thanks to the simultaneous analysis of nitrates and CFC which provides the groundwater residence time, heterotrophic
denitrification being further quantified through sulfate and nitrogen excess measurement (Supplementary Material). The Damköhler values show the limited but non-negligible amount of denitrification mostly in the lower part of the aquifer, even though dilution in the full groundwater volume plays a more important role on the resulting nitrate concentrations. In the next section, we investigate the consequences on the extrapolation capacities of the proposed modeling approach.

## 4.3 Retrospective and prospective evolutions of nitrate concentrations in rivers

We compare the observed and simulated nitrate concentrations in rivers obtained for the successfully calibrated models between 2000 and 2010 (inserts on Figure 7). The main trends of the concentration are well reproduced for the three catchments. The inter-annual variations are also well reproduced for the Ris and Douron catchments but they are comparatively too smooth for the Kerharo catchment. The downward trend since 2000 for the Ris and Kerharo catchments and the stability for the Douron catchment reflect the trends observed in the input nitrate concentrations of the groundwater recharge (Figure 4).

The calibrated models remain very close to each other in the calibration period (2000-2010) (Figure 7). They slightly differ when considered over the full simulation period 1955-2015 according to the denitrification time. We recall that the calibrated denitrification and transit times are correlated. Larger denitrification times are compensated by larger aquifer volumes and so by larger transit times. In other words, the lower denitrification is compensated by a higher dilution buffering the input concentrations. It results in systematic lower peaks and slower decreases of the nitrate concentrations for the high denitrification times obtained for the purple and green lines for the Ris catchment, and for the dark blue lines for the Douron and Kerharo catchments
(inserts on Figure 7). More in details, the calibrated models of the Ris catchment are within the range of the 1975-1985 observed concentrations (black dots on Figure 7) thanks to consistent evaluations of the mean transit times with the CFCs age data (see Supplementary Material).

Observed concentrations between 2010 and 2015 for the Ris and Kerharo catchments are also overall well predicted showing the relevance of the approach in a 5-year extrapolation exercise. Note that the additional 2010-2015 data for the Ris catchment tend to exclude the models without enough denitrification. Uncertainty and equifinality on the denitrification parameter and on the mean transit time would be reduced by extending the calibration period as already stated in other studies such as Kirchner (2016a, 2016b).
<figure><img src="image_9.png" /><figcaption>**Figure 7:** Modeled nitrate concentrations (colored lines) and observed nitrate concentrations (black lines and dots) compared to the nitrate concentration in the recharge (dashed lines). Blue lines refer to the retrospective concentrations. Red and green lines stand for the prospective “business as usual” (red) and “nitrate input sudden stop” (green) scenarios. The grey areas illustrate the predictions obtained when the</figcaption></figure>
threshold on the Nash-log criteria is lowered from 0.7 to 0.5. The insert for each catchment represents these calibrated models (colored lines) and the field observations (black line) focused on the calibration period. Data are presented at annual scale.

We further assess the long-term response of each catchment to two well-differentiated nitrate input scenarios from 2010-2015 to 2050 using the calibrated models. Scenarios begin in 2010 for the Douron catchment and in 2015 for the Ris and Kerharo catchments. The first scenario consists in constant nitrate inputs equal to the mean inputs of the 2010-2015 period for the Ris and Kerharo catchments, 2005-2010 for the Douron (dashed lines of Figure 7). It is called the “business as usual” scenario. The second scenario starts like the first scenario until 2025 when the input nitrate concentrations drop to 0 mg/L simulating a sudden ban of nitrate (green dashed lines of Figure 7). It is called the “nitrate input sudden stop” scenario. Both scenarios do not strictly represent actual socio-economic or target scenarios but constitute synthetic experiments to characterize the catchment response times. The “business as usual” scenario is used to assess the current trajectory. The “nitrate input sudden stop” scenario assesses the minimum time required to reach a given river concentration after a sudden drop of inputs. During the simulation period from 2010-2015 to 2050, the typical seasonal signal of the inputs is reproduced from the mean pattern over the existing data chronicle (monthly groundwater recharge rates and associated nitrate concentrations from 1997 to 2010). Consequently, inter-annual variability is removed from the input scenarios while seasonal fluctuations are still reproduced.

The predicted responses of the three sets of calibrated models to the two scenarios are presented by the red and green lines on Figure 7. The limited variability among simulations (at most 15% the mean river concentration of the “business as usual scenario”) shows that the
approach is relevant to predict the overall evolution of the nitrate concentrations. Moreover, lowering the acceptability threshold of the objective criterion *J* from 0.7 to 0.5 leads to similar results as illustrated by grey shaded areas (~100 models for each catchment) on **Figure 7**. This confirms the reliability of the approach. Variability is mainly due to the characteristic denitrification and transit times as detailed previously.

In the “business as usual” scenario, for the Douron river which inputs are relatively stable during the 1990-2010 period, an almost instantaneous quasi-stationary behavior is observed while Ris and Kerharo concentrations still gently decrease until 2050. The gap between the nitrate inputs (dashed lines) and outputs (solid lines) in the Douron case comes from the denitrification. In the two other cases, the evolving gap also comes from the progressive reduction of the legacy aquifer storage.

The persistence of nitrate in the aquifer is further illustrated by its reaction to the “nitrate input sudden stop” experiment as groundwater recharge is cleaned from any nitrate. The nitrate concentrations sharply drop in the river in the first year (6-10% for the Ris, 5-9% for the Kerharo, 11-14% for the Douron) as the direct contribution of nitrate from the surface excess overland flow ceases. The drop is logically higher for the Douron because of its larger seepage zone. The ensuing decrease is smoother and shows the persistence of nitrate in the aquifer. Ten years after the nitrate input stop, the removal rate is 32-51% for the Ris river, 29-47% for the Kerharo river and 56-71% for the Douron river. The faster removal for the Douron comes from its smaller transit time (**Table 3**).

The retrospective analysis shows the quality of the calibration. The prospective analysis highlights the effects of the different processes on the nitrate river concentrations. The interception with the surface controls the short-term release of the nitrate. The aquifer volume
determines the mean transit time of the nitrate and, in turn, the characteristic renewal time, while
the Damköhler number indicates the denitrification potential.

# 5 Discussion

Streamwater nitrate response is the result of past nitrate inputs convolved with catchments intrinsic geo(morpho)logical and biogeochemical properties. As developed in the three following sections, the successfully calibrated models presented in this paper allow us to further:

1. decompose the overall nitrate budget to characterize nitrate storage, removal and discharge;

2. track the nitrate location in the aquifer thanks to the spatially resolved modelling approach;

3. infer the role played by the relative organization of rock properties on denitrification potential, a basis for upscaling to other sites the characterization of nitrate legacy.

## 5.1 Catchment-scale nitrate storage and legacy

We determined the long-term catchment scale nitrate storage by an integrated budget of the nitrate entering, leaving and being degraded within the aquifer based on the calibrated models. The nitrate fluxes and budget are extracted from the model over the full 1955-2010 period for the three studied catchments. Results are presented as the overall quantity of nitrate discharged to the river, stored in the aquifer and removed by denitrification. All three terms are normalized by the integrated inlet quantity of nitrate over the same period (Table 4). We discuss successively the three terms: nitrate still stored in the aquifer, discharged to the river, and degraded in the aquifer.

**Table 4:** Nitrate mass denitrified in the aquifer, discharged to the river and stored in the aquifer integrated over the 1955-2010 period and normalized by the integrated inlet quantity of nitrate derived
from the successfully calibrated models on the three studied catchments. The mass balance is presented for three specific models: the value in bold characters refers to the model which denitrification rate is the median (Med) of the distribution obtained from the calibrated models. The other two sub-columns indicate the results obtained for models having the minimum (Min) and maximum (Max) denitrification rates. Note that, while the median value of denitrification (Med) is necessarily between the minimum and maximum, it does not have to be the case for river discharge and aquifer storage as min and max do not refer to these parameters.

<table><thead><tr><th rowspan="2">Nitrate budget over 1955-2010</th><th colspan="3">Ris</th><th colspan="3">Kerharo</th><th colspan="3">Douron</th></tr><tr><th>Min</th><th>Med</th><th>Max</th><th>Min</th><th>Med</th><th>Max</th><th>Min</th><th>Med</th><th>Max</th></tr></thead><tbody><tr><th>Denitrification %</th><td>0</td><td><strong>9</strong></td><td>11</td><td>13</td><td><strong>21</strong></td><td>22</td><td>10</td><td><strong>14</strong></td><td>14</td></tr><tr><th>River discharge %</th><td>61</td><td><strong>60</strong></td><td>64</td><td>49</td><td><strong>54</strong></td><td>51</td><td>65</td><td><strong>70</strong></td><td>69</td></tr><tr><th>Aquifer storage %</th><td>39</td><td><strong>31</strong></td><td>25</td><td>38</td><td><strong>25</strong></td><td>27</td><td>25</td><td><strong>16</strong></td><td>17</td></tr></tbody></table>

Median case

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

Mean water transit time : 21 yr

<figure><img src="image_12.png" /></figure>

Mean water transit time : 23 yr

<figure><img src="image_13.png" /></figure>

Mean water transit time : 10 yr

The first term is the nitrate still stored in the aquifer, obtained by integrating the concentration of all model's cells in 2010. It is determined more by the river discharge than by the denitrification as the river discharge is 3 to 6 times larger than the denitrification. The proportion of stored nitrate is larger for the Kerharo (27-38%), intermediate for the Ris (25-39%) and smaller for the Douron (17-25%) consistently with the quantity of water stored ($\theta E_{wet}$) increasing by a factor of around two from the Douron to the Kerharo (Table 3).
The second term is the nitrate discharged to the river computed from simulated streamflow and river concentration. It comes both from nitrate leaving the aquifer over the period and from nitrate that were inlet on saturated surfaces, did not enter the aquifer and flew directly to the river. The nitrate discharge rate was largest for the Douron (65-70%) and smallest for the Kerharo (49-54%). It is quite high, representing around two thirds of the total input to the aquifer. Uncertainties on discharged masses from 1955 to 2010 are limited thanks to the constrain provided by observed streamflows and river concentrations from 2000 to 2010. Note that nitrate discharged through the river after less than one year of transit time contributed only slightly to the output (10-21%) (Table 4) like water flows (5-15%) (Table 3). At regional to global scales, young water proportion might generally be larger as suggested by the analysis of Jasechko et al. (2016), who found a median value of 21% and a 10th - 90th percentile range of 4-53% on 254 watersheds around the world regarding the average contribution of river waters younger than three months. Comparing missing nitrogen (input minus river export) in the aquifer relatively to the quantity discharged in the river, we found a ratio of 0.4-1 while a value of 0.3-2.4 was estimated by Dupas et al. (2020) and 0.7 by Aquilina et al. (2012), both at Brittany scale. Typical retention rates including soil and aquifers ranges from 80 to 90% of total nitrogen input (Aquilina et al., 2012; Enhardt et al., 2019; Lassaletta et al., 2012).

The third term comes from the rates of denitrification and is deduced from the nitrate budget. They are comparable across the three different catchments. They amount to non-negligible 9-21% median values even though the denitrification rate is relatively low as shown by the small Damköhler number and high characteristic denitrification times of Table 3 although the estimated denitrification in the aquifer could be slightly overestimated as it integrates potential but limited heterotrophic in-stream denitrification. For the Douron and Kerharo,
denitrification rates reach 10-14% and 13-22% respectively. The Ris presents smaller denitrification rates, which might reflect less capacity for this granitic catchment to remove past accumulated nitrate. While the higher denitrified mass for the Kerharo comes with a longer mean transit time compared to Ris and Douron, the optimal (regarding criterion J) denitrified mass inside one catchment appears to be not correlated to its optimal mean transit time. Indeed, despite quite similar transit times for Ris and Kerharo (21-23 yr), denitrified mass is respectively 9% and 21% while the smaller transit time for the Douron (10 yr) coincides with a 14% denitrified mass. The comparison between the three catchments highlights that the denitrification potential is not only governed by the mean transit time but also by other factors. Denitrification cannot be interpreted as a uniform process within the aquifer volume like what might be done within a simplifying Damköhler framework, according to which any increase of the transit time would be traduced by additional denitrification. As developed below, denitrification additionally depends on the depth of the reduced zone. Above, denitrification would not occur whatever the transit time, while, below, denitrification would occur almost whatever the transit time. In such conditions, the Damköhler number characterizes the stratification of the reduced zone rather than the reaction time per se. While relevant in other hydrological compartments like the hyporheic and riparian zones, it might not be relevant to aquifers as previously hinted (Kolbe et al., 2019; Pinay et al., 2015). We further discuss this conclusion in the next sections by quantifying the impact of the nitrate concentrations stratification and by defining an aquifer denitrification potential.

## 5.2 Stratification of the nitrate storage and denitrification potential

Beyond the overall nitrate budget, the calibrated models can further be used to assess the spatial distribution of the nitrate concentrations and especially their stratification within the
aquifer. **Figure 8a** shows the evolution of the nitrate concentration with depth of the water height at the middle of the hillslope ($x=L/2$) for 3 000 models applied to the Ris catchment and selected randomly from the parameter space (for computational reason). The successfully calibrated models (materialized in orange and yellow) display very close nitrate stratification. The nitrate stratification pattern appears to be well constrained by the calibration data while it was not in the calibration target *J* (on **Figure 8a**, models with high *J* value display similar stratification). This result is obtained freely from the vertically resolved approach. We underline that stratification appears naturally through the increase of the transit time with depth and agrees with the typical nitrate profiles found in similar hydrogeological reservoirs (Faillat et al., 1999; Molenat et al., 2008; Molénat et al., 2002). It does not require any decrease of hydraulic conductivity with depth but a decrease of hydraulic conductivity would enhance the stratification. Lithological stratification is often assumed in other models of shallow Britain aquifers based on field evidences like in the soil and groundwater hillslope-scale two-linear reservoirs model of Fovet et al. (2015) or in the spatialized model of stratified hydraulic conductivity at catchment scale (Kolbe et al., 2016) and hillslope scale (Martin et al., 2006).
<figure><img src="image_14.png" /></figure>

**Figure 8:** Aquifer stratification of nitrate concentrations for the Pis catchment. Nitrate concentrations at the middle of the catchment (x=L/2) (horizontal axis) are represented against the equivalent water thickness taken as the product of the porosity by the area under the water table (vertical axis) (a) for 3000 models in 2005 with different model to data equilibriums represented by the calibration criterion J and (b) for one of the best models as a function of time. Note the reduction of the scale range 0-15 m for (a) to 0-8 m for (b). The mean groundwater transit time profile is added in red on (b). The best model represented on (b) is obtained for an aquifer thickness of 150 m, a porosity of 5%, an hydraulic conductivity of 1.8×10-5 m/s and a denitrification time of 100 yr. The dashed lines on (b) represent the simulated nitrate concentration in 1995 and 2010 when no denitrification is considered.

The nitrate vertical distribution is far from uniform, it is stratified. It globally echoes the input nitrate concentration shown in **Figure 4** in a deformed way as the groundwater transit time increases non-linearly with depth (**Figure 8b**, red line) and because denitrification increases with depth (**Figure 8b**, blue to green lines). The concentration peak is reached in 2005 at 5-6 meters corresponding to the peak of input concentrations of 1990 consistently obtained for a transit time of around 15 years. **Figure 8b** shows the vertical progression of the nitrate concentrations from 1970 to 2010 and the peak appearance in 1995 around two meters and deepening with time. The
decrease of the peak due to the denitrification is limited to around 20% in 15 years. The denitrification is much more effective deeper when the transit time sharply increases (Figure 8b, red line). Nitrate becomes fully degraded at the bottom of the aquifer model where the no-flow boundary condition slows down the transported elements long enough for the nitrate to be degraded. The limited amount of denitrification in most of the aquifer comes from the high characteristic depth at which 25% of the input nitrate are degraded intervening at 82-91% of the full equivalent water height of the aquifer (Table 3). It is also the case for the Douron and Kerharo with somewhat smaller but still high depths. In the case studied, most of the denitrification only occurs close to the basis of the modeled aquifer. We further discuss this important point in the next section.

## 5.3 Emergence of denitrification as a lithological interface process

Knowing where denitrification processes occur within aquifers and the drivers of this localization is crucial to predict denitrification efficiency. In the literature, two main features of aquifers have been previously pointed out for their role on denitrification efficiency: the volume -through residence time- and the existence of a geochemical interface (Kolbe et al., 2019; Liao et al., 2012; Tesoriero et al., 2005; Van Der Velde et al., 2010b; Van Meter and Basu, 2015). Based on the present results, these hypotheses are further questioned and we propose a new hypothesis based on a hydrodynamic interface to explain the denitrification processes within aquifers.

First, it is commonly accepted that the greater the volume, the greater the residence time and the greater the denitrification. However, here no correlations were detected between mean transit time and total denitrified mass when comparing the calibrated models. Indeed, the denitrified percentage on the Kerharo catchment is twice higher than on the Ris while they have
close mean transit times and the denitrified percentage on the Douron is higher than on the Ris while the mean transit time on the Douron is twice smaller (Table 4). In some cases, the percentage of denitrification could thus increase with a decreasing aquifer volume. We hypothesize that it is due to the flow decrease with depth, in other words, a shallower base of the aquifer would favor a higher overall denitrification rate. This assumption brings new insights into the control on denitrification of groundwater flow structure but should be further tested on different sites and at larger scales.

Second, recent studies highlighted the vertical stratification of denitrification reactivity because of geochemical interfaces in aquifers (Green et al., 2010; Kolbe et al., 2019). Autotrophic denitrification in crystalline basements is indeed closely linked to geology and weathering processes and thus highly related to geochemical interfaces. It occurs when oxygen levels are low enough for nitrates to become competitive electron acceptors and when reduced mineral are available as electron colors. Autotrophic denitrification has been commonly associated to sulfur minerals such as pyrite (Böhlke et al., 2002; Pinay et al., 2015) or iron-rich reduced minerals such as biotite (Aquilina et al., 2018). In this framework, denitrification is often considered null in the weathered layer because of prevalent oxic conditions and of weathered reduced minerals. The continuous water flows from the surface bring oxygen and nitrate progressively passivating down mineral surfaces through oxide precipitation, thus slowing down and limiting the denitrification reaction. This interpretation is also supported by experimental studies that have shown that a medium in which no apparent denitrification occurs can provide support for denitrification by reactivating fresh surfaces during rock grinding or during pumping at the field-scale (Roques et al., 2018b; Tarits et al., 2006).
This general framework of denitrification stratification is consistent with the results found in the present study, although we further suggest that a hydrodynamic interface might be a key driver in this stratification. Based on three study sites and successfully calibrated models of nitrate concentrations in both shallow and deep groundwater and river trends, we confirm that nitrate removal in shallow flow paths is very limited, while denitrification occurs deeper. The removal of nitrate is thus slow and does not prevent the build-up of a nitrogen legacy as previously suggested by Kolbe et al. (2019) on another aquifer in Brittany and complementary sites in California. The nitrate stratification is well simulated thanks to the simulated stratification of groundwater flows, showing that the circulation pattern exerts a large control on denitrification. The key point is that the calibrated depth of denitrification (or 25%-denitrification depth, see **Table 3**) is close to the base of the aquifer represented by the lowermost model limit. Similar results would be provided if permeability contrasts or denitrification time-lag were accounted for. Therefore, for the first time, we show here that denitrification processes are occurring deep in the aquifer, driven by the location of the lower boundary of the aquifer and thus by a hydrodynamic interface, which is also supported by the decorrelation between denitrification and mean transit time.

The role of a hydrodynamic interface on denitrification is not inconsistent with previous interpretations of a geochemical interface role and more data would be required to infer roles of both hydrodynamic and geochemical interfaces. This study together with previous studies nevertheless support the idea of co-evolution (Harman and Troch, 2014; Troch et al., 2013; Yoshida and Troch, 2016) where geochemical and hydrogeological properties are related. Denitrification is tightly related to the geochemical reaction along the flow-path but also to the reactive flow history of the aquifer. Geochemical reactions also alter rock properties and
particularly porosity and permeability, as such hydrodynamic and geochemical interfaces are
likely to coincide. Under natural flows, the basement of the aquifer would provide conductive
pathways but not too much developed to keep active electron donor sites. Therefore, lithology
interfaces, with their hydrodynamic and geochemical aspects, likely play a central role in
subsurface denitrification. Rather than an "all or nothing" interpretation, effectively slow and
limited but nonetheless active denitrification processes deep in aquifers should be considered in
nitrogen cycles and depollution studies.

**6 Conclusion: Implications for managers**

This study brings insights about our preliwve capacity regarding water streamflow and
water quality in rivers. A simple physically based groundwater model was used and informed by
streamflow and nitrate concentrations measured in rivers. Additionally, independent
measurements of anthropic atmospheric tracers were made on boreholes, wells and rivers to
derive transit times. The results on three different catchments also bring knowledge about typical
catchment response times and the control of groundwater on river flows and nitrate
concentrations. This study indicates that a significant part (5-15%) of the groundwater recharge
reach rapidly the river (< 1 yr) through the development of seepage areas and rapid flow in the
variously saturated zone. We found that mean transit time in groundwater are around decades
(10-32 yr). Regarding catchments capacity to remove nitrogen excess, an important result is that
denitrification appears limited to the lower part of the aquifer where flows are more limited and
transit times are longer. Eventually, the overall denitrification in the aquifer only reaches 10-20
%. Consequently, an important part of nitrate is still stored as nitrate legacy in groundwater
systems (Table 4). While the parsimonious model developed in this study fulfills its objective, we note that running the simulation until 2020 should reduce uncertainty. A deeper understanding of processes, focusing less on groundwater, should require to include a dynamic coupling between TNT2 and the hillslope-scale aquifer model, and a potential in-stream denitrification. However, it will strongly increase the model complexity and should require to use seasonal nitrate variations in river and additional observed data.

From a management perspective, this study contributes to answer the following key question: “How much time does slow groundwater transfer delay the effects of improved agricultural practices?”. The nitrate concentration in the river results from a combination of fast and slow transfers. Slow transfers (part of subsurface and deep subsurface waters) result from the aquifer storage capacity assessed above. Without any further input of nitrate, the river concentration would decrease by two thirds in 10 years for the Douron and by one third for the Kerharo, showing an important resilience of these catchments, capable of rapidly diminishing the legacy stored nitrate if no further inputs are constantly being added. The existence of rapid transfers (~1 yr) represents a major clue to farmers and managers. It indicates that the nitrate concentration in the river results from the past practices but also from the current practices (and potentially efforts of farmers). A sudden drop of nitrate inputs would lead to 5 to 25% and 15 to 50% decrease in the nitrate concentrations in rivers within 1 to 5 yr for the Kerharo and the Douron, respectively.

This study provides a precise quantitative timing of the effects of potential environmental management frameworks. It shows that the future of pollution, although partly related to past farming practices, is also the responsibility of present practices and environmental measures.
Acknowledgments

This study was funded by the “Région Bretagne” authorities, the "DREAL" Department (Ministery of Environment) and the Water Agency "Agence de l'Eau Loire-Bretagne" through the MORAQUI (MODélisation de la Réactivité des AQUIfères dans les bassins algues vertes) project. We thank the Public offices “EPA Baie de Douarnenez” and the “Communauté de Communes de Morlaix", and the Water Regional Technical Center " CRESEB" for data collection and field work as well as the “Condate-Eau” (CSU- Univ. Rennes 1) analytical platform for the chemical and dissolved gases analysis.

Bibliography

Abbott, B.W., Moatar, F., Gauthier, O., Fo, et, O., Antoine, V., Ragueneau, O., 2018. Trends and seasonality of river nutrients in agricultural catchments: 18 years of weekly citizen science in France. Sci. Total Environ. 624, 845-858. https://doi.org/10.1016/j.scitotenv.2017.12.176

Aquilina, L., Roques, C., Bonson, A., Vergnaud-Ayraud, V., Labasque, T., Pauwels, H., Pételet-
Giraud, E., Pettenat, M., Dufresne, A., Bethencourt, L., Bour, O., 2018. Autotrophic
denitrification supported by biotite dissolution in crystalline aquifers (1): New insights from
short-term batch experiments. Sci. Total Environ. 619-620, 842-853.
https://doi.org/10.1016/j.scitotenv.2017.11.079

Aquilina, L., Vergnaud-Ayraud, V., Labasque, T., Bour, O., Molénat, J., Ruiz, L., de Montety, V., De Ridder, J., Roques, C., Longuevergne, L., 2012. Nitrate dynamics in agricultural catchments deduced from groundwater dating and long-term nitrate monitoring in surface-
and groundwaters. Sci. Total Environ. 435-436, 167-178.

https://doi.org/10.1016/j.scitotenv.2012.06.028

Ayraud, V., Aquilina, L., Labasque, T., Pauwels, H., Molenat, J., Pierson-Wickmann, A.C., Durand, V., Bour, O., Tarits, C., Le Corre, P., Fourre, E., Merot, P., Davy, P., 2008.

Compartmentalization of physical and chemical properties in hard-rock aquifers deduced from chemical and groundwater age analyses. Appl. Geochemistry 23, 2686-2707.

https://doi.org/10.1016/j.apgeochem.2008.06.001

Bakker, M., Post, V., Langevin, C.D., Hughes, J.D., White, J T Starn, J.J., Fienen, M.N., 2016.
Scripting MODFLOW Model Development Using Python and FloPy. Groundwater 54,

733-739. https://doi.org/10.1111/gwat.12413

Beaujouan, V., Durand, P., Ruiz, L., 2001. Modeling the effect of the spatial distribution of agricultural practices on nitrogen fluxes in rural catchments. Ecol. Modell.

https://doi.org/10.1016/S0304-3800(00)00435-X

Beaujouan, V., Durand, P., Ruiz, L., Prousseau, P., Cotteret, G., 2002. A hydrological model dedicated to topography-based simulation of nitrogen transfer and transformation: Rationale and application to the geomorphology-denitrification relationship. Hydrol. Process. 16,

493-507. https://doi.org/10.1002/hyp.327

Bedekar, V., Morway, E., Langevin, C., Tonkin, M., 2016. MT3D-USGS version 1: A U.S.
Geological Survey release of MT3DMS updated with new and expanded transport
capabilities for use with MODFLOW. Groundw. Resour. Progr. 84.

https://doi.org/10.3133/tm6A53

Benettin, P., Bailey, S.W., Rinaldo, A., Likens, G.E., McGuire, K.J., Botter, G., 2017. Young runoff fractions control streamwater age and solute concentration dynamics. Hydrol.
Process. https://doi.org/10.1002/hyp.11243

Berghuijs, W.R., Kirchner, J.W., 2017. The relationship between contrasting ages of groundwater and streamflow. Geophys. Res. Lett. 44, 8925-8935.
https://doi.org/10.1002/2017GL074962

Boers, P.C.M., 1996. Nutrient emissions from agriculture in the Netherlands, causes and remedies. Water Sci. Technol. 33, 183-189. https://doi.org/10.1016/0273-1223(96)00229-6

Böhlke, J.K., Wanty, R., Tuttle, M., Delin, G., Landon, M., 2002. Denitrification in the recharge area and discharge area of a transient agricultural nutrient plume in a glacial outwash sand aquifer, Minnesota. Water Resour. Res. 38, 10-1-10-26.
https://doi.org/10.1029/2001wr000663

Boulton, A.J., Findlay, S., Marmonier, P., Stanley, E.H., Maurice Valett, H., 1998. The functional significance of the hyporheic zone in streams and rivers. Annu. Rev. Ecol. Syst. 29. https://doi.org/10.1146/annurev.ecolsys.29.1.59

Breemen, N. Van, Boyer, E.W., Paustian, K., Seitzinger, S., 2002. The Nitrogen Cycle at Regional to Global Scales. Nitrogen Cycle Reg. to Glob. Scales.
https://doi.org/10.1007/978-94-017-3405-9

Bresciani, E., Davy, P., de Dreuzy, J.-R., 2014. Is the Dupuit assumption suitable for predicting the groundwater seepage area in hillslopes? Water Resour. Res. 50, 2394-2406.
https://doi.org/10.1002/2013WR014284

Brutsaert, W., 1994. The unit response of groundwater outflow from a hillslope. Water Resour. Res. 30, 2759-2763.

Brutsaert, W., Nieber, J.L., 1977. Regionalized drought flow hydrographs from a mature
glaciated plateau. Water Resour. Res. 13, 637-643.

https://doi.org/10.1029/WR013i003p00637

Chen, D., Shen, H., Hu, M., Wang, J., Zhang, Y., Dahlgren, R.A., 2018. Legacy Nutrient Dynamics at the Watershed Scale: Principles, Modeling, and Implications, 1st ed, Advances in Agronomy. Elsevier Inc. https://doi.org/10.1016/bs.agron.2018.01.005

Clément, J.C., Holmes, R.M., Peterson, B.J., Pinay, G., 2003. Isotopic investigation of denitrification in a riparian ecosystem in western France. J. Appl. Ecol. 40, 1035-1048.
https://doi.org/10.1111/j.1365-2664.2003.00854.x

Conan, C., Bouraoui, F., Turpin, N., de Marsily, G., Bigoglio, G., 2003. Modeling Flow and Nitrate Fate at Catchment Scale in Brittany (France). J. Environ. Qual. 32, 2026-2032.
https://doi.org/10.2134/jeq2003.2026

Cornaton, F., Perrochet, P., 2007. Groundwater age, life expectancy and transit time distributions in advective – dispersive systems : 1. Generalized reservoir theory 29, 1267-1291. https://doi.org/10.1016/j.adwatres.2005.10.009

Danckwerts, P. V., 1953. Continuous flow systems. Distribution of residence times. Chem. Eng. Sci. 2, 1-13. https://doi.org/10.1016/0009-2509(53)80001-1

Deschamps, L., Frétière, K., Maillocchon, A., Molina, V., 2016. La Bretagne: première région française pour la production et la transformation de viande, INSEE.

Dupas, R., Ehrhardt, S., Musolff, A., Fovet, O., Durand, P., 2020. Long-term nitrogen retention and transit time distribution in agricultural catchments in western France. Environ. Res. Lett. 15. https://doi.org/10.1088/1748-9326/abe47

Dupas, R., Minaudo, C., Gruau, G., Ruiz, L., Gascuel-Odoux, C., 2018. Multidecadal Trajectory of Riverine Nitrogen and Phosphorus Dynamics in Rural Catchments. Water Resour. Res.
https://doi.org/10.1029/2018WR022905

Durand, P., Breuer, L., Johnes, P.J., Billen, G., Butturini, A., Pinay, G., van Grinsven, H., Garnier, J., Rivett, M., Reay, D.S., Curtis, C., Siemens, J., Maberly, S., Kaste, Ø., Humborg, C., Loeb, R., de Klein, J., Hejzlar, J., Skoulikidis, N., Kortelainen, P., Lepistö, A., Wright, R., 2011. Nitrogen processes in aquatic ecosystems. Eur. Nitrogen Assess. 126–146. https://doi.org/10.1017/cbo9780511976988.010

Earle, S., 2015. Physical geology, BC Open Textbooks. https://doi.org/10.1016/B978-0-444-42758-8.50008-8

Ehrhardt, S., Kumar, R., Fleckenstein, J.H., Attinger, S., Musolff, A., 2019. Trajectories of nitrate input and output in three nested catchments along a land use gradient. Hydrol. Earth Syst. Sci. 23, 3503–3524. https://doi.org/10.5194/ess-23-3503-2019

European Commission, 1991. Fighting water pollution from agricultural nitrates. 91/676/EEC.

Faillat, J.-P., Somlette, L., Sicard, T., 1999. Controles redox et hydrodynamiques des nitrates et zonation hydrochimique verticale des nappes de fissures du socle. Possibilité de généralisation. Bull. la Société géologique Fr. 170.

Fan, Y., Clark, M., Lawrence, D.M., Swenson, S., Band, L.E., Brantley, S.L., Brooks, P.D., Dietrich, W.E., Flores, A., Grant, G., Kirchner, J.W., Mackay, D.S., McDonnell, J.J., Milly, P.C.D., Sullivan, P.L., Tague, C., Ajami, H., Chaney, N., Hartmann, A., Hazenberg, P., McNamara, J., Pelletier, J., Perket, J., Rouholahnejad-Freund, E., Wagener, T., Zeng, X., Beighley, E., Buzan, J., Huang, M., Livneh, B., Mohanty, B.P., Nijssen, B., Safeeq, M., Shen, C., van Verseveld, W., Volk, J., Yamazaki, D., 2019. Hillslope Hydrology in Global Change Research and Earth System Modeling. Water Resour. Res. 1737–1772. https://doi.org/10.1029/2018WR023903
Fovet, O., Ruiz, L., Faucheux, M., Molénat, J., Sekhar, M., Vertès, F., Aquilina, L., Gascuel-Odoux, C., Durand, P., 2015. Using long time series of agricultural-derived nitrates for estimating catchment transit times. J. Hydrol. 522, 603–617.
https://doi.org/10.1016/j.jhydrol.2015.01.030

Franks, S.W., Gineste, P., Beven, K.J., Merot, P., 1998. On constraining the predictions of a distributed model: The incorporation of fuzzy estimates of saturated areas into the calibration process. Water Resour. Res. 34, 787–797. https://doi.org/10.1029/97WR03041

Gabriel, D., Roschewitz, I., Tscharntke, T., Thies, C., 2006. Biodiversity at different spatial scales: Plant communities in organic and conventional agriculture. Ecol. Appl. 16. https://doi.org/10.1890/1051-0761(2006)016[2011-ECADSS]2.0.CO;2

Galloway, J.N., Townsend, A.R., Erisman, J.W., Bekunda, M., Cai, Z., Freney, J.R., Martinelli, L.A., Seitzinger, S.P., Sutton, M.A., 2008. Transformation of the nitrogen cycle: Recent trends, questions, and potential solutions. Science (80-.). https://doi.org/10.1126/science.1146574

Gambino, M., 2014. Les mutations des systèmes productifs français : le modèle breton, à revisiter. Fr. les Mutations des systèmes Prod. 371–382.

Gburek, W.J., Folmar, G.J., 1999. Patterns of contaminant transport in a layered fractured aquifer. J. Contam. Hydrol. 37, 87–109. https://doi.org/10.1016/S0169-7722(98)00158-2

Gelhar, L.W., 1974. Stochastic analysis of phreatic aquifers. Water Resour. Res. 10, 539–545.
https://doi.org/10.1029/WR010i003p00539

Gleeson, T., Befus, K.M., Jasechko, S., Luijendijk, E., Cardenas, M.B., 2016. The global volume and distribution of modern groundwater. Nat. Geosci. 9, 161–164.
https://doi.org/10.1038/ngeo2590
Goderniaux, P., Davy, P., Bresciani, E., De Dreuzy, J.R., Le Borgne, T., 2013. Partitioning a regional groundwater flow system into shallow local and deep regional flow compartments. Water Resour. Res. 49, 2274–2286. https://doi.org/10.1002/wrcr.20186

Green, C.T., Böhlke, J.K., Bekins, B.A., Phillips, S.P., 2010. Mixing effects on apparent reaction rates and isotope fractionation during denitrification in a heterogeneous aquifer. Water Resour. Res. 46. https://doi.org/10.1029/2009WR008903

Green, C.T., Jurgens, B.C., Zhang, Y., Starn, J.J., Singleton, M.J., Esser, B.K., 2016. Regional oxygen reduction and denitrification rates in groundwater from a multi-model residence time distributions, San Joaquin Valley, USA. J. Hydrol. 543, 155–166. https://doi.org/10.1016/j.jhydrol.2016.05.018

Grimaldi, C., Thomas, Z., Fossey, M., Fauvel, Y., Meot, P., 2009. High chloride concentrations in the soil and groundwater under an oak hedge in the West of France: an indicator of evapotranspiration and water movement. Hydrol. Process. 23, 1865–1873. https://doi.org/10.1002/hyp.7316

Gu, C., Hornberger, G.M., Mills, A.L., Herman, J.S., Flewelling, S.A., 2007. Nitrate reduction in streambed sediments: Effects of flow and biogeochemical kinetics. Water Resour. Res. 43, 1–10. https://doi.org/10.1029/2007WR006027

Gupta, H. V., Kling, H., Yilmaz, K.K., Martinez, G.F., 2009. Decomposition of the mean squared error and NSE performance criteria: Implications for improving hydrological modelling. J. Hydrol. 377, 80–91. https://doi.org/10.1016/j.jhydrol.2009.08.003

Habets, F., Boone, A., Champeaux, J.L., Etchevers, P., Franchistéguy, L., Leblois, E., Ledoux, E., Le Moigne, P., Martin, E., Morel, S., Noilhan, J., Seguí, P.Q., Rousset-Regimbeau, F., Viennot, P., 2008. The SAFRAN-ISBA-MODCOU hydrometeorological model applied
over France. J. Geophys. Res. Atmos. 113, 1–18. https://doi.org/10.1029/2007JD008548

Haitjema, H.M., 1995. On the residence time distribution in idealized groundwatersheds. J. Hydrol. 172, 127–146. https://doi.org/10.1016/0022-1694(95)02732-5

Hamilton, S.K., 2012. Biogeochemical time lags may delay responses of streams to ecological restoration. Freshw. Biol. 57, 43–57. https://doi.org/10.1111/j.1365-2427.2011.02685.x

Harbaugh, A.W., 2005. MODFLOW 2005, The U. S. Geological Survey Modular GroundWater Model: the GroundWater Flow Process. U.S. Geol. Surv. Tech. Methods 253.

Harman, C., Troch, P.A., 2014. What makes Darwinian phylogeny “darwinian”? Asking a different kind of question about landscapes. Hydrol. Earth Syst. Sci. 18, 417–433. https://doi.org/10.5194/hess-18-417-2014

Harvey, J.W., Böhlke, J.K., Voytek, M.A., Scott, O., Tobias, C.R., 2013. Hyporheic zone denitrification: Controls on effective reaction depth and contribution to whole-stream mass balance. Water Resour. Res. 49, 6298–6316. https://doi.org/10.1002/wrcr.20492

Hiscock, K., 2009. Hydrogeology. Principles and Practice., John Wiley & Sons. https://doi.org/10.1007/s1265-016-5360-8

Hrachowitz, M., Benettin, P., Breukelen, B.M., Fovet, O., Howden, N.J.K., Ruiz, L., van der Velde, Y., Wade, A.J., 2016. Transit times—the link between hydrology and water quality at the catchment scale. Wiley Interdiscip. Rev. Water 3, 629–657. https://doi.org/10.1002/wat2.1155

Hrachowitz, M., Fovet, O., Ruiz, L., Savenije, H.H.G., 2015. Transit time distributions, legacy contamination and variability in biogeochemical $1/\alpha$ scaling: how are hydrological response dynamics linked to water quality at the catchment scale? Hydrol. Process. 29, 5241–5256. https://doi.org/10.1002/hyp.10546
Jasechko, S., Kirchner, J.W., Welker, J.M., McDonnell, J.J., 2016. Substantial proportion of global streamflow less than three months old. Nat. Geosci. 9, 126–129.
https://doi.org/10.1038/ngeo2636

Jimenez-Martinez, J., Longuevergne, L., Le Borgne, T., Davy, P., Russian, A., Bour, O., 2013. Temporal and spatial scaling of hydraulic response to recharge in fractured aquifers: Insights from a frequency domain analysis. Water Resour. Res. 49, 3007–3023.
https://doi.org/10.1002/wrcr.20260

Jurgens, B.C., Bohlke, J.K., Eberts, S.M., 2012. An Excalibur workbook for interpreting groundwater age distributions from environmental tracer data. U.S. Geol. Surv. Tech. Methods Rep. 4-F3 60.

Kaandorp, V.P., de Louw, P.G.B., van der Velden, V., Broers, H.P., 2018. Transient Groundwater Travel Time Distributions and Age-Ranked Storage-Discharge Relationships of Three Lowland Catchments. Water Resour. Res. 54, 4519–4536.
https://doi.org/10.1029/2017WR00261

Kirchner, J.W., 2016a. Aggregation in environmental systems-Part 1: Seasonal tracer cycles quantify young water tracers, but not mean transit times, in spatially heterogeneous catchments. Hydrol. Earth Syst. Sci. 20, 279–297. https://doi.org/10.5194/hess-20-279-2016

Kirchner, J.W., 2016b. Aggregation in environmental systems-Part 2: Catchment mean transit times and young water fractions under hydrologic nonstationarity. Hydrol. Earth Syst. Sci. 20, 299–328. https://doi.org/10.5194/hess-20-299-2016

Kirchner, J.W., Feng, X., Neal, C., 2000. Fractal stream chemistry and its implications for contaminant transport in catchments. Nature 403, 524–527.
https://doi.org/10.1038/35000537
Knowles, R., 1982. Denitrification. Microbiol. Rev. 46, 866-871. https://doi.org/10.1016/B978-008045405-4.00264-0

Kolbe, T., 2017. Temporal and spatial structures of denitrification in crystalline aquifers.

Kolbe, T., de Dreuzy, J.-R., Abbott, B.W., Aquilina, L., Babey, T., Green, C.T., Fleckenstein, J.H., Labasque, T., Laverman, A.M., Marçais, J., Peiffer, S., Thomas, Z., Pinay, G., 2019. Stratification of reactivity determines nitrate removal in groundwater. Proc. Natl. Acad. Sci. 116, 2494-2499. https://doi.org/10.1073/pnas.1816892116

Kolbe, T., Marçais, J., Thomas, Z., Abbott, B.W., de Dreuzy, J.R., Rousseau-Gueutin, P., Aquilina, L., Labasque, T., Pinay, G., 2016. Coupling 3D groundwater modeling with CFC-based age dating to classify local groundwater circulation in an unconfined crystalline aquifer. J. Hydrol. 543, 31-46. https://doi.org/10.1016/j.jhydrol.2016.05.020

Korom, S.F., 1992. Natural denitrification in the saturated zone: A review. Water Resour. Res. https://doi.org/10.1029/92WR00252

Kovacs, G., 1981. Seepage hydraulics. Elsevier. https://doi.org/10.1016/0022-1694(84)90254-3

Kronvang, B., Andersen, H.E., Tørgesen, C., Dalgaard, T., Larsen, S.E., Bøgestrand, J., Blicher-Mathiasen, G., 2008. Effects of policy measures implemented in Denmark on nitrogen pollution of the aquatic environment. Environ. Sci. Policy. https://doi.org/10.1016/j.envsci.2007.10.007

Kronvang, B., Jeppesen, E., Conley, D.J., Søndergaard, M., Larsen, S.E., Ovesen, N.B., Carstensen, J., 2005. Nutrient pressures and ecological responses to nutrient loading reductions in Danish streams, lakes and coastal waters. J. Hydrol. 304, 274-288. https://doi.org/10.1016/j.jhydrol.2004.07.035

Lague, D., Davy, P., Crave, A., 2000. Estimating Uplift Rate and Erodibility from the Area-
Slope Relationship: Examples from Britanny ( France ) and Numerical Modelling. Phys. Chem. Earth Part A Solid Earth Geod. 25, 543–548. https://doi.org/10.1016/S1464-1895(00)00083-1

Lassaletta, L., Romero, E., Billen, G., Garnier, J., García-Gómez, H., Rovira, J. V., 2012. Spatialized N budgets in a large agricultural Mediterranean watershed: High loading and low transfer. Biogeosciences 9, 57–70. https://doi.org/10.5194/bg-9-57-2012

Le Borgne, T., Bour, O., Paillet, F.L., Caudal, J.P., 2006. Assessment of preferential flow path connectivity and hydraulic properties at single-borehole and cross-borehole scales in a fractured aquifer. J. Hydrol. 328, 347–359. https://doi.org/10.1016/j.jhydrol.2005.12.029

Le Moigne, P., 2012. SURFEX scientific documentation. Geosci. Model Dev. 237.

Lefebvre, S., Clément, J.C., Pinay, G., Thenail, C., Durand, P., Marmonier, P., 2007. 15N-nitrate signature in low-order streams: Effects of land cover and agricultural practices. Ecol. Appl. 17. https://doi.org/10.1890/06-1496.1

Legchenko, A., Baltassat, J.M., Bobacay, A., Martin, C., Robain, H., Vouillamoz, J.M., 2004. Magnetic resonance sounding applied to aquifer characterization. Ground Water 42, 363–373. https://doi.org/10.1111/j.1745-6584.2004.tb02684.x

Leray, S., de Dreuzy, J.R., Aquilina, L., Vergnaud-Ayraud, V., Labasque, T., Bour, O., Le Borgne, T., 2014. Temporal evolution of age data under transient pumping conditions. J. Hydrol. 511, 555–566. https://doi.org/10.1016/j.jhydrol.2014.01.064

Liao, L., Green, C.T., Bekins, B.A., Böhlke, J.K., 2012. Factors controlling nitrate fluxes in groundwater in agricultural areas. Water Resour. Res. 48. https://doi.org/10.1029/2011WR011008

Małoszewski, P., Zuber, A., 1982. Determining the turnover time of groundwater systems with
the aid of environmental tracers. 1. Models and their applicability. J. Hydrol. 57, 207-231.
https://doi.org/10.1016/0022-1694(82)90147-0

Marçais, J., de Dreuzy, J.R., Erhel, J., 2017. Dynamic coupling of subsurface and seepage flows solved within a regularized partition formulation. Adv. Water Resour. 109, 94-105.
https://doi.org/10.1016/j.advwatres.2017.09.008

Marçais, J., de Dreuzy, J.R., Ginn, T.R., Rousseau-Gueutin, P., Leray, S., 2015. Inferring transit time distributions from atmospheric tracer data: Assessment of the predictive capacities of Lumped Parameter Models on a 3D crystalline aquifer model. J. Hydrol. 525, 619-631.
https://doi.org/10.1016/j.jhydrol.2015.03.055

Marçais, J., Gauvain, A., Labasque, T., Abbott, B.W., Finay, G., Aquilina, L., Chabaux, F., Viville, D., De Dreuzy, J.-R., 2018. Dating groundwater with dissolved silica and CFC concentrations in crystalline aquifers. Sci. Total Environ. 636, 260-272.

Maréchal, J.C., Dewandel, B., Subrahmanyam, K., 2004. Use of hydraulic tests at different scales to characterize fracture network properties in the weathered-fractured layer of a hard rock aquifer. Water Resour. Res. 40, 1-17. https://doi.org/10.1029/2004WR003137

Martin, C., Aquilina, L., Gascuel-Odoux, C., Molénat, J., Faucheux, M., Ruiz, L., 2004. Seasonal and interannual variations of nitrate and chloride in stream waters related to spatial and temporal patterns of groundwater concentrations in agricultural catchments. Hydrol. Process. 18, 1237-1254. https://doi.org/10.1002/hyp.1395

Martin, C., Molénat, J., Gascuel-Odoux, C., Vouillamoz, J.M., Robain, H., Ruiz, L., Faucheux, M., Aquilina, L., 2006. Modelling the effect of physical and chemical characteristics of shallow aquifers on water and nitrate transport in small agricultural catchments. J. Hydrol. 326, 25-42. https://doi.org/10.1016/j.jhydrol.2005.10.040
Masson, V., Le Moigne, P., Martin, E., Faroux, S., Alias, A., Alkama, R., Belamari, S., Barbu, A., Boone, A., Bouyssel, F., Brousseau, P., Brun, E., Calvet, J.C., Carrer, D., Decharme, B., Delire, C., Donier, S., Essaouini, K., Gibelin, A.L., Giordani, H., Habets, F., Jidane, M., Kerdraon, G., Kourzeneva, E., Lafaysse, M., Lafont, S., Lebeaupin Brossier, C., Lemonsu, A., Mahfouf, J.F., Marguinaud, P., Mokhtari, M., Morin, S., Pigeon, G., Salgado, R., Seity, Y., Taillefer, F., Tanguy, G., Tulet, P., Vincendon, B., Vionnet, V., Voldoire, A., 2013. The SURFEXv7.2 land and ocean surface platform for coupled o offline simulation of earth surface variables and fluxes. Geosci. Model Dev. 6, 929–943. https://doi.org/10.5194/gmd-6-929-2013

Matonse, A.H., Kroll, C., 2009. Simulating low streamflows with hillslope storage models. Water Resour. Res. 45, 1–13. https://doi.org/10.1029/2007WR006529

McDonald, M.G., Harbaugh, A.W., 1974. A modular three dimensional finite difference groundwater flow model: U.S. Geological Survey Open File Report 83 875, 528 p.

Meals, D.W., Dressing, S.A., Davlopert, T.E., 2010. Lag Time in Water Quality Response to Best Management Practices: A Review. J. Environ. Qual. 39, 85–96. https://doi.org/10.2139/qc.2009.0108

MEDDE, GIS Sol, 2014. Enveloppes des milieux potentiellement humides de la France métropolitaine. Notice d’accompagnement. Programme de modélisation des milieux potentiellement humides de France, Ministère de l’Ecologie, du Développement Durable et de l’Energie, Groupement d’I.

Ménesguen, A., Salommon, J., 1988. Eutrophication modelling as a tool for fighting against Ulva coastal mass blooms. Pap. Present. Comput. Model. Ocean Eng. Venice.

Merot, P., Corgne, S., Delahaye, D., Desnos, P., Dubreuil, V., Gascuel, C., Giteau, J., 2014.
Assessment, impact and perception of climate change in the western part of France: The CLIMASTER project 23, 96–107. https://doi.org/10.1684/agr.2014.0694

Merot, P., Squividant, H., Aurousseau, P., Hefting, M., Burt, T., Maitre, V., Kruk, M., Butturini, A., Thenail, C., Viaud, V., 2003. Testing a climato-topographic index for predicting wetlands distribution along an European climate gradient. Ecol. Modell. 163, 51–71. https://doi.org/10.1016/S0304-3800(02)00387-3

Molénat, J., Davy, P., Gascuel-Odoux, C., Durand, P., 1999. Study of three subsurface hydrologic systems based on spectral and cross-spectral analysis of time series. J. Hydrol. 222, 152–164. https://doi.org/10.1016/S0022-1694(99)00107-9

Molénat, J., Durand, P., Gascuel-Odoux, C., Davy, P., Gruau, G., 2002. Mechanisms of nitrate transfer from soil to stream in agricultural watershed of French Brittany. Water. Air. Soil Pollut. 133, 161–183.

Molénat, J., Gascuel-Odoux, C., Aquilino, L., Ruiz, L., 2013. Use of gaseous tracers (CFCs and SF6) and transit-time distribution spectrum to validate a shallow groundwater transport model. J. Hydrol. 480, 1–9. https://doi.org/10.1016/j.jhydrol.2012.11.043

Molenat, J., Gascuel-Odoux, C., Ruiz, L., Gruau, G., 2008. Role of water table dynamics on stream nitrate export and concentration in agricultural headwater catchment (France). J. Hydrol. 348, 363–378. https://doi.org/10.1016/j.jhydrol.2007.10.005

Montreuil, O., Merot, P., Marmonier, P., 2010. Estimation of nitrate removal by riparian wetlands and streams in agricultural catchments: Effect of discharge and stream order. Freshw. Biol. 55. https://doi.org/10.1111/j.1365-2427.2010.02439.x

Mougin, B., Allier, D., Blanchin, R., Carn, A., Courtois, N., Gateau, C., Putot, E., 2008. SILURES Bretagne (Système d' Information pour la Localisation et l' Utilisation des
Ressources en Eaux Souterraines).

Mougin, B., Carn, A., Jegou, J.-P., Quéméner, G., 2006. SILURES Bretagne (Système d' Information pour la Localisation et l' Utilisation des Ressources en Eaux Souterraines).

Ocampo, C.J., Oldham, C.E., Sivapalan, M., 2006a. Nitrate attenuation in agricultural catchments: Shifting balances between transport and reaction. Water Resour. Res. 42, 1–16.
https://doi.org/10.1029/2004WR003773

Ocampo, C.J., Sivapalan, M., Oldham, C.E., 2006b. Field exploration of coupled hydrological and biogeochemical catchment responses and a unifying conceptual model. Adv. Water Resour. 29, 161–180. https://doi.org/10.1016/j.advwatres.2005.02.014

Ogden, F.L., Watts, B.A., 2000. Saturated area formation in nonconvergent hillslope topography with shallow soils: A numerical investigation. Water Resour. Res. 36, 1795–1804.
https://doi.org/10.1029/2000WR900001

Oldham, C.E., Farrow, D.E., Peiffer, S., 2013. A generalized Damköhler number for classifying material processing in hydrological systems. Hydrol. Earth Syst. Sci. https://doi.org/10.5194/hes-17-i133-2013

Pinault, J.-L., Pauwels, H., 2001. Inverse modeling of the hydrological and the hydrochemical behavior of hydro systems: Application to nitrate transport and denitrification. Water Resour. 37, 2179–2190.

Pinay, G., Peiffer, S., De Dreuzy, J.R., Krause, S., Hannah, D.M., Fleckenstein, J.H., Sebilo, M., Bishop, K., Hubert-Moy, L., 2015. Upscaling Nitrogen Removal Capacity from Local Hotspots to Low Stream Orders’ Drainage Basins. Ecosystems 18, 1101–1120.
https://doi.org/10.1007/s10021-015-9878-5

Poisvert, C., Curie, F., Moatar, F., 2017. Annual agricultural N surplus in France over a 70-year
period. Nutr. Cycl. Agroecosystems 107, 63–78. https://doi.org/10.1007/s10705-016-9814-x

Quintana-Seguí, P., Le Moigne, P., Durand, Y., Martin, E., Habets, F., Baillon, M., Canellas, C., Franchisteguy, L., Morel, S., 2008. Analysis of near-surface atmospheric variables: Validation of the SAFRAN analysis over France. J. Appl. Meteorol. Climatol. 47, 92–107. https://doi.org/10.1175/2007JAMC1636.1

Roques, C., Aquilina, L., Boisson, A., Vergnaud-Ayraud, V., Labasque, T., Longuevergne, L., Laurencelle, M., Dufresne, A., de Dreuzy, J.R., Pauwels, H., Bour, O., 2018a. Autotrophic denitrification supported by biotite dissolution in crystalline aquifers: (2) transient mixing and denitrification dynamic during long-term pumping. Sci. Total Environ. 619–620, 491–503. https://doi.org/10.1016/j.scitotenv.2017.11.104

Roques, C., Aquilina, L., Boisson, A., Vergnaud-Ayraud, V., Labasque, T., Longuevergne, L., Laurencelle, M., Dufresne, A., de Dreuzy, J.R., Pauwels, H., Bour, O., 2018b. Autotrophic denitrification supported by biotite dissolution in crystalline aquifers: (2) transient mixing and denitrification dynamic during long-term pumping. Sci. Total Environ. 619–620, 491–503. https://doi.org/10.1016/j.scitotenv.2017.11.104

Roques, C., Aquilina, L., Bour, O., Maréchal, J.C., Dewandel, B., Pauwels, H., Labasque, T., Vergnaud-Ayraud, V., Hochreutener, R., 2014. Groundwater sources and geochemical processes in a crystalline fault aquifer. J. Hydrol. 519, 3110–3128. https://doi.org/10.1016/j.jhydrol.2014.10.052

Ruiz, L., Abiven, S., Martin, C., Durand, P., Beaujouan, V., Molénat, J., 2002. Effect on nitrate concentration in stream water of agricultural practices in small catchments in Brittany: II. Temporal variations and mixing processes. Hydrol. Earth Syst. Sci. 6, 507–513. https://doi.org/10.5194/hess-6-507-2002
Singhal, B.B.S., Gupta, R.P., 2010. Applied hydrogeology of fractured rocks: Second edition, Springer Science & Business Media. https://doi.org/10.1007/978-90-481-8799-7

Steffen, W., Richardson, K., Rockström, J., Cornell, S.E., Fetzer, I., Bennett, E.M., Biggs, R., Carpenter, S.R., De Vries, W., De Wit, C.A., Folke, C., Gerten, D., Heinke, J., Mace, G.M., Persson, L.M., Ramanathan, V., Reyers, B., Sörlin, S., 2015. Planetary boundaries: Guiding human development on a changing planet. Science (80-). 347. https://doi.org/10.1126/science.1259855

Takuya, K., Yoshishige, H., Nobuharu, K., Tadakatsu, Y., Yasuo, C., 1993. Dispersion effect on the apparent nitrogen isotope fractionation factor associated with denitrification in soil; Evaluation by a mathematical model. Soil Biol Biochem. https://doi.org/10.1016/0038-0717(93)90134-W

Tarits, C., Aquilina, L., Ayraud, V., Pauwels, H., Davy, P., Touchard, F., Bour, O., 2006. Oxidoreduction sequence related to flux variations of groundwater from a fractured basement aquifer (Ploemeur area, France). Appl. Geochemistry 21, 29-47. https://doi.org/10.1016/j.apgeocnem.2005.09.004

Tesoriero, A.J., Spruill, T. R., Mew, H.E., Farrell, K.M., Harden, S.L., 2005. Nitrogen transport and transformations in a coastal plain watershed: Influence of geomorphology on flow paths and residence times. Water Resour. Res. 41, 1-15. https://doi.org/10.1029/2003WR002953

Townley, L.R., 1995. The response of aquifers to periodic forcing. Adv. Water Resour. 18, 125-146. https://doi.org/10.1016/0309-1708(95)00008-7

Troch, P. a, Paniconi, C., Emiel van Loon, E., 2003. Hillslope-storage Boussinesq model for subsurface flow and variable source areas along complex hillslopes: 1. Formulation and characteristic response. Water Resour. Res. 39, 1316.
https://doi.org/10.1029/2002wr001728

Troch, P.A., Carrillo, G., Sivapalan, M., Wagener, T., Sawicz, K., 2013. Climate-vegetation-soil interactions and long-term hydrologic partitioning: Signatures of catchment co-evolution. Hydrol. Earth Syst. Sci. 17, 2209–2217. https://doi.org/10.5194/hess-17-2209-2013

Van Der Velde, Y., De Rooij, G.H., Rozemeijer, J.C., Van Geer, F.C., Broers, H.P., 2010a. Nitrate response of a lowland catchment: On the relation between stream concentration and travel time distribution dynamics. Water Resour. Res. 46, 1–17. https://doi.org/10.1029/2010WR009105

Van Der Velde, Y., De Rooij, G.H., Rozemeijer, J.C., Van Geer, F.C., Broers, H.P., 2010b. Nitrate response of a lowland catchment: On the relation between stream concentration and travel time distribution dynamics. Water Resour. Res. 46, 1–17. https://doi.org/10.1029/2010WR009105

Van Meter, K.J., Basu, N.B., 2015. Catchment legacies and time lags: A parsimonious watershed model to predict the effects of legacy storage on nitrogen export. PLoS One 10, 1–22. https://doi.org/10.1371/journal.pone.0125971

Van Meter, K.J., Basu, N.B., Van Cappellen, P., 2017. Two centuries of nitrogen dynamics: Legacy sources and sinks in the Mississippi and Susquehanna River Basins. Global Biogeochem. Cycles 31, 2–23. https://doi.org/10.1002/2016GB005498

Van Meter, K.J., Basu, N.B., Veenstra, J.J., Burras, C.L., 2016. The nitrogen legacy: Emerging evidence of nitrogen accumulation in anthropogenic landscapes. Environ. Res. Lett. 11. https://doi.org/10.1088/1748-9326/11/3/035014

Vautier, C., Abhervé, R., Labasque, T., Laverman, A.M., Guillou, A., Chatton, E., Dupont, P., Aquilina, L., de Dreuzy, J.R., 2020. Mapping gas exchanges in headwater streams with
membrane inlet mass spectrometry. J. Hydrol. 581, 124398.

https://doi.org/10.1016/j.jhydrol.2019.124398

Vogel, J.C., 1967. Investigation of groundwater flow with radiocarbon, in: Isotopes in Hydrology.

Vrugt, J.A., Sadegh, M., 2013. Toward diagnostic model calibration and evaluation: Approximate Bayesian computation. Water Resour. Res. 49, 4335–4345.

https://doi.org/10.1002/wcr.20354

Wade, A.J., Durand, P., Beaujouan, V., Wessel, W.W., Raat J., Whitehead, P.G., Butterfield, D., Rankinen, K., Lepisto, A., 2002. A nitrogen model for European catchments: INCA, new model structure and equations. Hydro Earth Syst. Sci. 6, 559–582.

https://doi.org/10.5194/hess-6-559-2002

Withers, P.J.A., Neal, C., Jarvie, H.P., Dundy, D.G., 2014. Agriculture and eutrophication: Where do we go from here? Sustain. 6, 5853–5875. https://doi.org/10.3390/su6095853

Wriedt, G., Rode, M., 2006. Modeling nitrate transport and turnover in a lowland catchment system. J. Hydrol. 328, 157–176. https://doi.org/10.1016/j.jhydrol.2005.12.017

Wright, E.P., Burgess, W.G., 1992. The hydrogeology of crystalline basement aquifers in Africa, Geological Society, London, Special Publications.

Wyns, R., Baltassat, J.M., Lachassagne, P., Legchenko, A., Vairon, J., Mathieu, F., 2004. Application of proton magnetic resonance soundings to groundwater reserve mapping in weathered basement rocks (Brittany, France). Bull. la Soc. Geol. Fr. 175, 21–34.

https://doi.org/10.2113/175.1.21

Yoshida, T., Troch, P.A., 2016. Coevolution of volcanic catchments in Japan. Hydrol. Earth Syst. Sci. 20, 1133–1150. https://doi.org/10.5194/hess-20-1133-2016
Zarnetske, J.P., Haggerty, R., Wondzell, S.M., Baker, M.A., 2011. Dynamics of nitrate production and removal as a function of residence time in the hyporheic zone. J. Geophys. Res. Biogeosciences 116, 1–12. https://doi.org/10.1029/2010JG001356

Zarnetske, J.P., Haggerty, R., Wondzell, S.M., Bokil, V.A., González-Pinzón, R., 2012. Coupled transport and reaction kinetics control the nitrate source-sink function of hyporheic zones. Water Resour. Res. https://doi.org/10.1029/2012WR011894

Zheng, C., Wang, P., 1999. MT3DMS: A modular three dimensional multispecies transport model for simulation of advection, dispersion, and chemical reactions of contaminants in groundwater systems. Technical report, Waterways Experiment Station, US Army Corps of Engineers. A Modul. three-dimensional multi-species ... 239.
Credit Author Statement

**Luca Guillaumot** : Data Collection, Conceptualization, Methodology, Software, Data Curation, Writing - Original Draft, Writing - Review & Editing, Visualization

**Jean Marçais** : Methodology, Writing - Review & Editing, Visualization

**Camille Vautier** : Data Collection, Writing - Original Draft

**Aurélie Guillou** : Data Curation, Data Collection

**Virginie Vergnaud** : Data Curation

**Camille Bouchez** : Writing - Review & Editing

**Rémi Dupas** : Writing - Review & Editing

**Patrick Durand** : Methodology, Software, Writing - Review & Editing

**Jean-Raynald de Dreuzy** : Conceptualization, Methodology, Writing - Review & Editing

**Luc Aquilina** : Conceptualization, Methodology, Writing - Review & Editing, Project administration, Funding acquisition
**Declaration of interests**

[x] The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

[ ] The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:

<figure><img src="image_15.png" /></figure>
Abstract

The long-term fate of agricultural nitrate depends on rapid subsurface transfer, denitrification and storage in aquifers. Quantifying these processes remains an issue due to time varying subsurface contribution, unknown aquifer storage and heterogeneous denitrification potential.

Here, we develop a parsimonious modelling approach that uses long-term discharge and river nitrate concentration time-series combined with groundwater age data determined from chlorofluorocarbons in springs and boreholes. To leverage their informational content, we use a Boussinesq-type equivalent hillslope model to capture the dynamics of aquifer flows and evolving surface and subsurface contribution to river nitrate transport was modelled with a depth-resolved high-order finite-difference method and denitrification by a first-order law. We applied the method to three heavily nitrate loaded catchments of a crystalline temperate region of France (Brittany). We found that mean water transit time ranged 10-32 years and Damköhler ratio (transit time / denitrification time) ranged 0.12-0.55, leading to limited denitrification in the aquifer (10-20%). The long-term trajectory of nitrate concentration in rivers appears determined by flows stratification in the aquifer. The results suggest that autotrophic denitrification is controlled by the accessibility of reduced minerals which occurs at the base of the aquifer where flows decrease. One interpretation is that denitrification might be an interfacial process in zones that are weathered enough to transmit flows and not too weathered to have remaining accessible reduced minerals. Consequently, denitrification would not be controlled by the total aquifer volume and related mean transit time but by the proximity of the active weathered interface with the water table. This should be confirmed by complementary studies to which the developed methodology might be further deployed.
<figure><img src="image_16.png" /></figure>

Graphical abstract
Highlights

* A parsimonious model constrained by discharge, nitrate and CFC concentrations.

* Aquifer volume and young water contribution govern nitrate trends in rivers.

* Groundwater transit time stratification explains nitrate trends in rivers.

* Weak denitrification (~10-20%) occurs only in the deepest part of the aquifer.