# Strong hydroclimatic controls on vulnerability to subsurface nitrate contamination across Europe

Rohini Kumar, Falk Heße, P. Suresh C. Rao, Andreas Musolff, James W. Jawitz, Fanny J. Sarrazin, Luis Samaniego, Jan H Fleckenstein, Oldrich Rakovec, Stephan Thober, et al.

► To cite this version:

Rohini Kumar, Falk Heße, P. Suresh C. Rao, Andreas Musolff, James W. Jawitz, et al. Strong hydroclimatic controls on vulnerability to subsurface nitrate contamination across Europe. Nature Communications, 2020, 11 (1), pp.6302. 10.1038/s41467-020-19955-8. hal-04760679

HAL Id: hal-04760679

https://hal.science/hal-04760679v1

Submitted on 31 Oct 2024

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>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

ARTICLE

https://doi.org/10.1038/s41467-020-19955-8

OPEN

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

Check for updates

# Strong hydroclimatic controls on vulnerability to subsurface nitrate contamination across Europe

R. Kumar<sup>1,8</sup>, F. Heße<sup>1</sup>, P. S. C. Rao<sup>2,3</sup>, A. Musolff<sup>1</sup>, J. W. Jawitz<sup>4</sup>, F. Sarrazin<sup>1</sup>, L. Samaniego<sup>1</sup>, J. H. Fleckenstein<sup>1,5</sup>, O. Rakovec<sup>1,6</sup>, S. Thober<sup>1</sup> & S. Attinger<sup>1,7</sup>

Subsurface contamination due to excessive nutrient surpluses is a persistent and widespread problem in agricultural areas across Europe. The vulnerability of a particular location to pollution from reactive solutes, such as nitrate, is determined by the interplay between hydrologic transport and biogeochemical transformations. Current studies on the controls of subsurface vulnerability do not consider the transient behaviour of transport dynamics in the root zone. Here, using state-of-the-art hydrologic simulations driven by observed hydroclimatic forcing, we demonstrate the strong spatiotemporal heterogeneity of hydrologic transport dynamics and reveal that these dynamics are primarily controlled by the hydroclimatic gradient of the aridity index across Europe. Contrasting the space-time dynamics of transport times with reactive timescales of denitrification in soil indicate that ~75% of the cultivated areas across Europe are potentially vulnerable to nitrate leaching for at least one-third of the year. We find that neglecting the transient nature of transport and reaction timescale results in a great underestimation of the extent of vulnerable regions by almost 50%. Therefore, future vulnerability and risk assessment studies must account for the transient behaviour of transport and biogeochemical transformation processes.

<footnote n="1"><sup>1</sup>UFZ-Helmholtz Centre for Environmental Research, Leipzig, Germany. <sup>2</sup>Lyles School of Civil Engineering, Purdue University, West Lafayette, IN, USA.</footnote>

<footnote n="3"><sup>3</sup>Agronomy Department, Purdue University, West Lafayette, IN, USA. <sup>4</sup>Soil and Water Sciences Department, University of Florida, Gainesville, FL, USA.</footnote>

<footnote n="5"><sup>5</sup>Bayreuth Center of Ecology and Environmental Research, University of Bayreuth, Leipzig, Germany. <sup>6</sup>Faculty of Environmental Sciences, Czech University of Life Sciences, Prague, Czech Republic. <sup>7</sup>University of Potsdam, Potsdam, Germany. <sup>8</sup>email: rohini.kumar@ufz.de</footnote>
large-scale vulnerability assessment of subsurface nitrate contamination. Our study therefore focuses on Europe-wide vulnerability assessment, i.e., the potential for (excess) nitrate leaching from the root zone to deeper in the subsurface (i.e., vadose zone below rooting depth). We demonstrate the oversimplified (static) nature of previous vulnerability assessment approaches by highlighting the relevance of the transient nature of transport dynamics, and we discuss its ramifications for future assessment and subsequent policy decisions. Our continental-scale analysis demonstrates strong spatiotemporal heterogeneity of hydrologic transport dynamics pronounced throughout the European landscapes, and we show that the (static) vulnerability assessment approach that does not account for such transient features greatly underestimates the extent of vulnerable areas prone to subsurface contamination by excess nitrate leaching.

## Results and discussion

**Space-time variability of hydrologic transport times.** Our continental-scale hydrologic simulations show large space-time heterogeneity in the inferred TTDs, which illustrates the complex, non-linear, and transient nature of transport dynamics in the root zone (Fig. 1; see also Supplementary Video). The large spread among the simulated daily TTDs in three exemplary locations (Fig. 1a–c) illustrates the pronounced (space-time) heterogeneity of hydroclimatic factors (e.g., precipitation, soil-water storage, infiltration and evapotranspiration fluxes), and characterizes the different transport dynamics inferred across Europe. These locations represent the humid, sub-humid (transitional), and semi-arid climate regimes, aridity indices ($\phi$ = ratio between mean potential evapotranspiration and mean precipitation) of 0.25 (UK), 1.15 (France), and 3.25 (Spain), respectively. A consistent shift towards longer travel times is noticed when moving from humid to semi-arid regions. The daily travel-time interquartile range (TT<sub>IQR</sub>) increases jointly with the median travel-time (TT<sub>50</sub>) throughout Europe (Fig. 1d–f and see Supplementary Fig. 1). The humid location (UK; Fig. 1a–d) exhibits marked seasonality with higher TT<sub>50</sub> in summer, which is likely due to a regular seasonal pattern of hydrologic states/fluxes, specifically soil moisture and evapotranspiration. In contrast, the daily dynamics of the TT<sub>50</sub> at the semi-arid location (Spain; Fig. 1c, f) are more erratic and episodic in nature. This location is characterized by infrequent rainfall that, when combined with high evapotranspiration losses, leads to highly variable soil-water storage. The location in France marks a transition zone between humid and semi-arid locations, with a seasonal pattern during wet years and an erratic pattern during dry years (Fig. 1b, e). The distinct behavior of the TT dynamics simulated among different locations broadly agrees with past theoretical understandings, even though these previous efforts used synthetic datasets<sup>29,37</sup>.

**Synthesis of hydrologic transport times across Europe.** Figure 2a, b summarizes the spatial patterns of the temporal mean $\mu$(TT<sub>50</sub>) and standard deviation $\sigma$(TT<sub>50</sub>) of the daily soil-water travel times across Europe. The ranges for both the mean and the standard deviation span factors of 6 to 7 across the continent, as 99% of the values lie between 100 and 700 days for $\mu$(TT<sub>50</sub>) and 52 and 320 days for $\sigma$(TT<sub>50</sub>). Fifty percent of the study domain has $\mu$(TT<sub>50</sub>) and $\sigma$(TT<sub>50</sub>) values that exceed 365 and 120 days, respectively. Regions with the shortest soil travel times ($\mu$(TT<sub>50</sub>) ≤ 180 days) are located in areas of high and frequent rainfall (e.g., alpine, northern UK and northern Spain— Pyrenees). The longest travel times ($\mu$(TT<sub>50</sub>) ≥ 540 days) are found in dry areas with less frequent rainfall in southern Spain and the eastern European regions adjacent to the Black Sea. The spatial patterns (Fig. 2a, b) suggest a strikingly high spatial similarity

Despite >15 years of water quality protection implementation under the EU Water Framework Directive (EU-WFD¹), the most recent EU-WFD report² concludes that the majority of European water bodies do not meet the European Union's minimum target, with threats coming from a wide range of pollutants. Among these, excess nitrate from agricultural areas was highlighted as a major concern³⁻⁸. Consequently, the European Nitrate Directive⁹—itself an integral part of the EU-WFD—designates nitrate vulnerable zones (NVZs) as areas at risk from agricultural nitrate pollution and requires prompt actions to improve nitrate management. A number of indices have been developed to delineate these zones¹⁰⁻¹³. While these indices differ in their conceptual and implementation modes, they are often based on a weighted combination of temporally invariant environmental parameters (e.g., terrain slope, land cover and subsurface properties, mean precipitation). A framework for the delineation of such NVZs based on an integrated understanding of the complex and dynamic interplay between hydrologic transport and biogeochemical turnover is still missing.

A major challenge for such a framework is to capture the hydrologic transport capacity or the intrinsic vulnerability to subsurface contamination by diffuse pollutants¹¹,¹⁴. Subsurface transport is particularly elusive and uncertain due to its complex flow patterns. To account for this uncertainty, research has focused on the statistical characterization of transport dynamics through travel-time distributions (TTDs), which capture the journey of water and dissolved solutes through a given subsurface compartment¹⁵⁻²⁰. Much work has been based on steady-state TTDs, however, more recently, studies have started to acknowledge the transient nature of TTDs²¹⁻²⁷. Typically, such studies have focused on empirical observations at the catchment scale or at a limited number of densely gauged small-scale catchments. While transport dynamics have recently been investigated at larger-scales²⁸, there are, however, no studies that systematically examine the transient nature of travel times, identify the main driving forces, and connect them to the reactive behavior of (diffuse) pollutants at regional to continental scales. This information would be relevant for management and decision making.

To address this gap, we provide a Europe-wide assessment of hydrologic transport behavior as an integrated measure of the intrinsic vulnerability to subsurface contamination by diffuse pollutants (e.g., nitrate). We then use the case of widespread nitrate contamination across arable lands in Europe to show the unrecognized importance of the transient nature of hydrologic transport in previous vulnerability and risk assessments. Our analysis is based on state-of-the-art continental-scale hydrologic simulations driven by meteorological observations over the period of 1950–2015 combined with the recent theoretical developments for characterizing the transient nature of hydrologic transport dynamics²¹,²²,²⁶,²⁷ at high spatial and temporal resolutions (0.25° and daily timescale; see “Methods”). We focus on the root zone because it is the interface between the land surface and deeper subsurface. This zone is the most dynamic and active part of the subsurface and acts as both a hydrologic and a biogeochemical filter, determining the delivery and turnover of surface inputs and the partitioning of flow paths to near and deeper subsurface waters²⁹⁻³¹. The rooting depth varies across European landscapes depending on, among other geophysical attributes, vegetation types and groundwater table depth³². However, our focus in this vulnerability assessment is limited to arable landscapes, we therefore account for the dynamics of the first meter of soil that mostly coincides with the rooting zone for arable lands³³. We use the dimensionless Damköhler number³⁴⁻³⁶ to link the hydrologic and biogeochemical timescales (see “Methods”) and provide an objective measure for the
<figure><img src="image_4.png" /><figcaption>**Fig. 1** Transient features of the hydrologic transport dynamics. Illustration of the daily travel-time distributions (a-c), and the corresponding temporal evolution of the daily median $TT_{50}$ and interquartile ranges $TT_{IQR}$ (d-f) for three distinct locations in the UK, France, and Spain representing the sub-humid, humid or transitional, and semi-arid hydroclimatic settings across Europe, respectively.</figcaption></figure>

between $\mu(TT_{50})$ and $\sigma(TT_{50})$. This was confirmed through a linear regression analysis (Fig. 2d), resulting in $R^2 = 0.72$ and 0.97 for raw and binned data ($p$-value < 0.00001). Consequently, European regions that have, on average, longer soil travel times are also more variable or episodic in time, and vice versa. This result also means that the coefficient of variation (CV) of the daily $TT_{50}$, which is the regression slope between $\mu(TT_{50})$ and $\sigma(TT_{50})$, is remarkably consistent across Europe (~0.4). An analogous consistency in CV values has been also reported in a previous study³⁷ for catchment-scale travel times across different hydroclimatic settings with synthetic datasets. The temporal variability in median travel times is strongly controlled by variability in climate, and we find that the value of $TT_{50}$ CV = 0.4 is in good agreement with the temporal variability of daily potential evapotranspiration (CV = 0.34, see Supplementary Fig. 2).

The rather well-organized spatial patterns of $\mu(TT_{50})$ and $\sigma(TT_{50})$ follow the hydroclimatic gradient observed across Europe (Fig. 2c), and here, the latter pattern is represented through the aridity index ($\phi$) that primarily controls the partitioning of incoming rainfall and energy into outgoing water fluxes (i.e., evapotranspiration vs. runoff). Approximately 70–73% of the variance in the Europe-wide estimates of $\mu(TT_{50})$ and $\sigma(TT_{50})$ can be explained solely by the spatial heterogeneity of $\phi$ (Fig. 2e, f). The spatial correlation structure for these travel-time characteristics and the aridity index are nearly identical (Supplementary Fig. 3). The close relation of $\mu(TT_{50})$ and $\sigma(TT_{50})$ to $\phi$ emphasizes that the high soil moisture and frequent rainfall conditions in

humid regions lead to a high connectivity and fast displacement of water in the soil column²⁹. In contrast, the (semi-)arid regions with infrequent rainfall and high soil-water deficits due to high evapotranspiration losses generally exhibit long travel times.

Other site-specific landscape attributes related to terrain, soil and vegetation characteristics show a weaker correspondence to the spatial heterogeneity of $\mu(TT_{50})$ and $\sigma(TT_{50})$ compared to $\phi$ and, therefore, constitute secondary controls (see Supplementary Fig. 4). Notably, the combined effect of individual secondary factors appears to be stable across the range of $\phi$ values, as is visible in the scatter of points in Fig. 2e, f; as well as by the overlying standard deviation estimates of the respective binned data (for every $\phi$ interval of 0.15). For the (average) binned estimates, we find an almost perfect linear correlation of $\phi$ with $\mu(TT_{50})$ and $\sigma(TT_{50})$ ($R^2 = 0.96–0.98$; $p$-value < 0.00001). These results are also found for the extremes of the daily TT distributions. For example, $TT_{10}$ and $TT_{90}$, which are indicative of young and old water fractions, respectively, also show a high spatial correlation with $\phi$ (see Supplementary Figs. 5 and 6). This result underpins the dominant role of the hydroclimatic factor ($\phi$) in shaping the dominant features of hydrologic transport timescales inferred across Europe.

**Vulnerability to nitrate leaching across the cultivated areas of Europe.** Thus far, our analysis has focused on the hydrologic transport dynamics that represent the intrinsic vulnerability¹⁴ of
<figure><img src="image_5.png" /><figcaption>b</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>c</figcaption></figure>

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

<figure><img src="image_8.png" /><figcaption>e</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>f</figcaption></figure>

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

**Fig. 2** Synthesis of the hydrologic transport time dynamics across Europe. The transient feature of the daily median travel times (TT<sub>50</sub>; blue lines in Fig. 1d, e, f) is summarized as their temporal mean, μ(TT<sub>50</sub>) and standard deviation, σ(TT<sub>50</sub>) (a, b) for the period 1985–2015. The strong spatial correspondence between μ(TT<sub>50</sub>) and σ(TT<sub>50</sub>) is evident through the point-wise correlation analysis (d). The prevailing hydroclimatic feature in the form of aridity index, φ in (c), the ratio between the mean potential evapotranspiration (Ē<sub>p</sub>) and mean precipitation (P̄), is identified as the dominant factor controlling the spatial heterogeneity of the transient transport characteristics as μ(TT<sub>50</sub>) and σ(TT<sub>50</sub>) simulated across Europe (e, f). On each of the scatter plots, along with point-wise cloud data, there is also the corresponding bin estimates given as the mean and one standard deviation values of grouped data for every φ interval of 0.15. Binning is performed with the aim of seeking generalized relationships (i.e., after reducing the noise in scatter due to outliers) and specifically in (e, f) to depict the role of secondary (landscape-related) factors through the representation of the binned standard deviation estimates that are relatively stable and present across the whole range of φ values.

the system to subsurface contamination. Here, we complement the transport timescales with the biogeochemical turnover timescales of nitrate in soil to determine the extent of vulnerable regions to nitrate leaching across the cultivated areas of Europe. We focus on two competing processes of excess nitrate removal (after consideration of plant uptake) from the soil by contrasting the timescales of nitrate leaching (hydrologic transport) and denitrification (biogeochemical turnover). Denitrification rates are poorly constrained due to a lack of reliable observations, especially at large scales. To acknowledge this uncertainty, we consider different characteristic denitrification timescales³⁸, 〈RT₅₀〉, defined here as the 50% removal of the initial substrate, to allow comparability with the median transport times (TT₅₀). We consider a range of 〈RT₅₀〉 values between 0.5 and 5 years in our analysis³⁹⁻⁴³ (see “Methods” for more details). 〈RT₅₀〉 represents the effective timescale encapsulating relevant environmental factors, such as soil moisture, temperature, and organic carbon content, that affect the site-specific reaction rates⁴²,⁴⁴. In the following text, we analyze two cases and conduct a nitrate vulnerability assessment corresponding to the static (time-averaged) and transient behaviors of the hydrologic transport (TT) and denitrification (RT) timescales.

We connect the transport and denitrification timescales through the dimensionless Damköhler number (here defined as $D_a = \frac{TT_{50}}{RT_{50}}$) that enables us to assess the interplay between these

two competing processes across the geographical domain³⁴⁻³⁶,³⁸. When $D_a < 1$, transport (leaching) dominates over reaction (denitrification), and vice versa. Our static vulnerability assessment case, based on the range of 〈RT₅₀〉 (0.5–5 years) and the averaged transport times μ(TT₅₀), results in $D_a$ values ranging between 0.05 and 4.0 across the majority of European cultivated areas (Fig. 3a). To interpret this $D_a$ number, we rely on a series of prior studies. First, a prior study³⁴ presents a variety of field observational datasets that demonstrate an empirical (non-linear) relationship between $D_a$ and nitrate removal. Subsequent studies³⁶,⁴⁵ showed how this empirical relationship can be described by a parametric model based on the exponential function. Using the latter, the above value of $D_a = 0.05$ would imply more than 90% of the nitrate leaching from soil, while a value of $D_a = 4.0$ would correspond to less than 10% of the nitrate leaching from soil. Furthermore our results suggest that the majority of the cultivated areas in Europe would be highly vulnerable to nitrate leaching ($D_a \le 1$) under the higher 〈RT₅₀〉 value (≥ 2 years), while very few locations would be classified as vulnerable in the case of the lowest 〈RT₅₀〉 value of 0.5 year due to the dominance of the denitrification timescale.

We further elaborate on the moderate case of an 〈RT₅₀〉 = 1 year, which is a highly plausible value, as inferred by the values of denitrification timescales in previous large-scale³⁹,⁴²,⁴⁶ and catchment-scale studies⁴⁰,⁴³. Using this 〈RT₅₀〉 value and the
<figure><img src="image_11.png" /></figure>

Static case

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

Transient case

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

Cropland [circa 2000, %]

D$_{a}$ [-]

<5 5-30 30-50 50-70 >70

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

Frequency of D$_{a}$(t) ≤ 1

**Fig. 3** Subsurface nitrate vulnerability assessment across Europe under static and transient considerations of hydrologic transport and denitrification timescales. **a** Box-plots summarizing the spatial distribution of the Damköhler number ($D_a = \frac{TT_{50}}{RT_{50}}$) for a range of effective denitrification timescales $\langle RT_{50} \rangle$ and the averaged transport times $\mu(\text{TT}_{50})$ estimated across the cultivated areas of Europe$^{72}$. Boxplot is displayed with the horizontal bar at the median, the box indicates the first and third quartiles and the whiskers indicate $\pm$1.5 times the interquartile range. Grid cells with at least 5% of the cropland areas are considered cultivated areas in the analysis. **b** Europe-wide $D_a$ estimates for the static vulnerability assessment case corresponding to the effective reaction timescale, $\langle RT_{50} \rangle$ and the average transport timescales, $\mu(\text{TT}_{50})$. Areas with $D_a \le 1$ are vulnerable to subsurface nitrate leaching and -42% of the total cultivated area falls within this category under the static vulnerability assessment case. **c** Summary of the transient daily $D_a(t)$ estimates for every cultivated grid cell, arranged according to their corresponding aridity index values ($\phi$). Summary statistics are presented as quantiles of the daily $D_a$ values (Q0.1,1,5, ...,95,99,99.9) for the actual estimates (in background) and the LOESS (locally weighted smoothing)-derived smooth statistics (in a transparent foreground for better clarity). **d** Frequency of the daily $D_a(t) \le 1$ estimated based on the time-varying $TT_{50}(t)$ and $RT_{50}(t)$. Considering the transient nature of $D_a$, nearly 75% of cultivated areas of Europe become vulnerable to nitrate leaching ($D_a(t) \le 1$) for at least one-third of the year (i.e., frequency estimates $\ge$ 0.33 or 4/12).

average transport time $\mu(\text{TT}_{50})$, we find ~42% of the cultivated areas across Europe to be vulnerable to subsurface nitrate leaching, with $D_a \le 1$ (Fig. 3b). Figure 3b shows that the cultivated areas in the humid and transitional climate zones ($\phi \le 1.5$), specifically in the western part of Europe, i.e., France, Germany,

Italy, UK and Ireland, are dominated by hydrologic transport ($D_a \le 1$) and are thus vulnerable to nitrate leaching. In contrast, the agricultural areas in the Iberian Peninsula and eastern European countries present higher dominance of denitrification over the hydrologic transport processes. The vulnerable areas
delineated by our approach based on the Europe-wide $D_a$ map are remarkably consistent with the (nitrate) leaching risk potential map published by the European Commission⁴⁷ and this map was established following a static indices-based approach (see Supplementary Fig. 7 for more details).

We now contrast the above static characterization of nitrate vulnerability with the results of the transient case based on the time-varying daily $TT_{50}(t)$ and $RT_{50}(t)$. Here, we account for the spatiotemporal variability of the environmental factors, $f_E(t)$, that affect the daily dynamics of the denitrification timescale, $RT_{50}(t)$. We include the effect of varying soil moisture, air temperature, and organic carbon content following established parameterization schemes³⁹⁻⁴¹,⁴³,⁴⁸ (see “Methods” for more details). This allowed us to construct the temporal variability of $RT_{50}(t)$ for each grid cell by explicitly accounting for the space/time variability of the environmental factor ($f_E$) (see “Methods”). The daily mean and variability of the environmental factor ($f_E$) follows the hydroclimatic gradient of the aridity index $\phi$ observed across Europe (see Supplementary Fig. 8). We observe lower mean and variability of the daily $f_E$ in (semi-)arid areas, resulting from the relatively drier soil moisture conditions that persist for long periods, than in humid areas that have on average higher $f_E$ because of wetter and more strongly seasonal soil moisture dynamics (Supplementary Fig. S8).

Next we analyze the variability of the daily $D_a(t)$ that depicts the interplay between the hydrologic transport timescale ($TT_{50}$) and the reactive timescale ($RT_{50}$) behaviors for the excess nitrate removal (leaching vs. denitrification) through the soil. Here, in Fig. 3c we summarize the daily $D_a(t)$ as the quantile estimates for every cropland cell arranged according to their respective aridity index ($\phi$) value. For example, for a given cropland cell located in a climate region of $\phi = 1$ and having the 50th percentile (or median; Q50) value of the daily $D_a = 1.0$, this number would indicate that the given cell would be prone to excess nitrate leaching (i.e., $D_a \le 1.0$) for nearly half of simulation times (on average six months of a year). The results depicted in Fig. 3c clearly show the increasing range of the daily $D_a$ variability (e.g., between Q95 and Q5) when moving from humid to semi-arid and arid regions. We observe a nearly twofold (100%) increase in the $D_a$ range (Q95-Q5) for cropland cells with $\phi$ of 1 to 3. Interestingly, while the averages (medians) of the daily $D_a$ values are usually high (above 1) in (semi-) arid regions, the $D_a$ values also exhibit high temporal variability and, therefore, frequently fall below a critical value of 1 (i.e., the favorable times for (excess) nitrate leaching). The correlation analysis suggests a strong correspondence between the averaged $D_a$ statistics (median and Q95-Q5) and the aridity index ($\phi$) across the European croplands ($R^2 \ge 0.92$).

We now analyze how the temporal dynamics of the daily $D_a(t)$ affects our vulnerability assessment. The results of the frequency analysis based on the time-varying $D_a \le 1$ suggest that ~75% of the cropland cells across Europe are vulnerable to nitrate leaching for at least one-third of the year (Fig. 3f; see also Supplementary Fig. 9). Our estimate of potentially vulnerable areas is nearly twice the estimate (42%) obtained above under the static consideration of transport and reaction timescales. Importantly, the cultivated areas located in the Iberian Peninsula and eastern European countries, which were not recognized as vulnerable regions under the static assumptions (Fig. 3b), are now regarded as regions that are temporarily vulnerable to subsurface nitrate contamination (Fig. 3c, d). We find that the majority of cropland cells in Europe (>90%) are prone to nitrate leaching for at least two months of the year (see Supplementary Fig. 9). Therefore, our results highlight the limitation of the static vulnerability assessment approach⁴⁷, which leads to an underestimation of nitrate vulnerable regions and has serious implications for nitrate management across Europe.

**Concluding remarks.** In this study, we provide a Europe-wide assessment of transient hydrologic transport behavior in the upper one meter of the subsurface. This approach allows us to quantify the intrinsic vulnerability of the subsurface to contaminant leaching at unprecedented spatial and temporal resolutions at the continental scale. We demonstrate the dominant role of large-scale hydroclimatic factors, as expressed in the aridity index, in determining the spatial heterogeneity of transport characteristics (e.g., temporal mean and variability) and the environmental factors that affect the daily variability of denitrification timescales. We apply an approach based on the dimensionless Damköhler number to characterize the vulnerability of subsurface waters to (excess) nitrate leaching from soil by accounting for the complex and dynamic interplay between the hydrologic transport and biogeochemical (denitrification) reaction timescales. This approach provides a general framework to objectively assess the vulnerability to other agrochemical pollutants in different subsurface compartments (e.g., root zone, deeper vadose zone, and eventually shallow and deep groundwater).

Using this framework as a decision tool for subsurface nitrate contamination assessments, our study closes a pressing gap by making the recent progress in the field of transport dynamics accessible to practitioners, regulators, and decision-makers that aim to safeguard and restore European waters. Our results/ framework can be used to identify hot spots and times vulnerable to (excess) nitrate leaching through the soil, and thereby could be used for assisting (nitrate) management strategies, such as the optimization/regulation of fertilizer applications. Our emphasis on the transient aspects in the process of defining vulnerable zones will become more important due to the projected increases in the frequency and intensity of extreme hydroclimatic events under changing climate conditions⁴⁹⁻⁵¹.

Our results call for improved vulnerability assessment approaches in Europe and other regions of intensive agriculture. Current practices that do not consider transient dynamics could lead to a substantial underestimation of the extent of vulnerable areas and the associated risk. Our study addresses this limitation and thus provides crucial vulnerability criteria, which can then be combined with information on the exposure (i.e., available data on excess nitrate components) for a risk assessment. To this end we provide a showcase example of vulnerability assessment combining information on excess nitrate (see “Methods”). We, therefore, urge modelers and planners to further develop and evaluate tools that use transient dynamics to assess the vulnerability of the subsurface to diffuse pollution (e.g., excessive nutrient surpluses in the root zone). Further avenues of research include the improvement, validation, and quantification of the predictive uncertainty of this kind of vulnerability assessment. Improved vulnerability maps are fundamental for deciding on agricultural subsidies and nitrate management, which are key components of the EU’s common agricultural policy (CAP).

## Methods

**Continental-scale hydrologic simulations.** We use the spatially explicit, process-based, mesoscale hydrologic model (mHM⁵²,⁵³) to perform the continental-scale hydrologic simulations over Europe. The model features the multiscale parameter regionalization (MPR) technique, which explicitly accounts for the spatial heterogeneity of fine-scale terrain, soil, vegetation, and other landscape properties. The mHM-MPR modeling framework provides a unique capability to simultaneously and seamlessly operate the model at multiple scales and locations⁵²⁻⁵⁴. We briefly describe the processes within the root-zone soil compartment, relevant for this study, in the following text. Here, we account for the dynamics of the first meter of soil that coincides with the rooting zone for arable lands. This compartment is modeled as three consecutive layers with the following depths: 0–5, 6–25, and 26–100 cm. In each layer, the incoming water in the form of rainfall plus snowmelt, after accounting for the canopy interception, is partitioned into soil-water storage and exfiltration based on a non-linear (power) function depending on the degree of
saturation of the corresponding layer, following the conceptualizations used in other large-scale models⁵⁵,⁵⁶. The exfiltrated water from the first layer is input to the second layer and from the second to the third layer. The evapotranspiration losses from each layer are modeled as a fraction of potential evapotranspiration and are dependent on water storage deficit-induced stress and the fraction of the vegetation roots in each layer (see Supplementary Fig. 10 for a conceptualization of these root-zone soil moisture accounting processes). Readers interested in a complete model description may refer to previous studies⁵²,⁵⁴,⁵⁷, and www.ufz.de/mhm for the source code and a detailed user manual. A comprehensive overview of underlying datasets, processing steps, as well as model establishment including impact assessment of natural/human intervention activities (e.g., irrigation) is detailed in Supplementary Note 1 and in Supplementary Table 1. Model simulations are produced at a 0.25° spatial resolution and on a daily timescale for the period 1950–2015, using the model parameterizations established in previous studies⁵³,⁵⁴,⁵⁸ (see also the Supplementary Note 1 for more details). We conducted a thorough multivariate evaluation of model performance across the European domain (see Supplementary Note 2 for details). While the choice of spatial resolution used here is constrained by the availability of meteorological forcing datasets (E-OBS; v13.0)⁵⁹, the multiscale parameterization approach implemented in mHM allows for the explicit treatment of fine-scale, sub-grid variability of landscape features⁵²,⁵³. We emphasize that our study focuses on providing a general framework for the characterization of subsurface nitrate vulnerability and applied it to the pan-European landscape. This framework can also be applied/expanded to finer scales (e.g., catchment scale), incorporating more detailed datasets²⁷ and relevant processes to provide crucial insights into local vulnerability and assist policy intervention strategies. Examples are the localized effects of irrigation (see Supplementary Note 1 for details) and artificial (tile) drainage, which can impact soil-water storage and indirectly the root-zone transport dynamics.

Hydrologic components of mHM has been also coupled to a nitrate transport model and previous studies have demonstrated successful applications of this coupled model⁴³,⁶⁰. The formulation of reactive timescales related to denitrification process in soil used in this study follows a similar approach to that of the coupled model (i.e., a first-order denitrification in soil along with a space-time variability of environmental factors; see the Section below on "Nitrate vulnerability assessment" for more details).

**Derivation of hydrologic transport times.** We follow recent theoretical developments to infer the time-variant nature of the travel-time distributions (TTDs)²¹,²⁵,²⁶,⁶¹⁻⁶³. Specifically, we adopt the notion by Botter et al.,²¹,⁶¹ who provided an elegant expression for deriving the time-variant TTDs for water parcels entering or leaving a control volume based on the temporal evolution of water storages and fluxes under several mixing (age function) schemes. We use the transient formulations of TTDs on a control volume taken as a single-grid cell and soil layer, which is characterized by the daily dynamics of the soil-water storage (S) and incoming flux as effective precipitation f (snowmelt plus rainfall minus canopy interception losses) and outgoing flux as O = (I + E), where I and E represent the exfiltration and evapotranspiration fluxes from a given soil layer, respectively. Under a random sampling scheme of mixing that assigns the same probability to all water particles with different ages in storage to be sampled by outgoing fluxes, the analytical expression for the transient TTD at any time t for water parcels exiting (as the exfiltration flux) can be expressed as follows²¹:

$$P_t(t - t_{in}|t_{in}) = \frac{I(t)}{\theta(t_{in})S(t)} \exp \left( - \int_{t_{in}}^{t} \frac{I(t') + E(t')}{S(t')} dt' \right) \quad (1)$$

with $t - t_{in}$ ($t > t_{in}$) represents the time from the moment the water parcel enters the control volume ($t_{in}$) until now (t). The partition function $\theta(t_{in})$ indicates the portion of the water parcel that enters the control volume at $t_{in}$ and leaves as the exfiltration flux (specifically, I as opposed to E) and is expressed as follows:

$$\theta(t_{in}) = \int_{t_{in}}^{\infty} \frac{I(t)}{S(t)} \exp \left( - \int_{t_{in}}^{t} \frac{I(t') + E(t')}{S(t')} dt' \right) dt. \quad (2)$$

This partition function $\theta$ defines a dimensionless number between [0 and 1]. The above expression is formulated for the TTDs conditioned to the entrance (or injection) time of water particles to the control volume, and therefore, it relates to the concept of life expectancy that tracks the ages of water particles forward in time⁶⁴,⁶⁵. The complementary approach in which the ages of the water particles exiting the system are tracked back in time relates to an age concept⁶⁵, and both forward and backward TTs can be related through the Niemi's continuity equation⁶⁶. Recent studies²²,²³,²⁵,²⁶,⁶³,⁶⁷,⁶⁸ demonstrated the usefulness of the transient TTDs in capturing the overall behavior of hydrological and geochemical responses in experimental and intensively monitored catchments.

We numerically solve the expressions to derive the daily evolution of the Europe-wide TTDs based on the water fluxes and storages of the three soil layers simulated by mHM for the period 1985–2015. The procedure to derive the overall TTD representing the entire root zone (0–1 m) is carried out in two steps. First, we derive the time-varying TTDs for each soil layer separately using the layer-specific storages and outgoing water fluxes at every modeling time-step. The overall TTD of water parcels leaving the entire root zone at a 1 m depth is then estimated by

sequential convolution of the independently estimated probability density functions for the first, second, and third soil layers. Since this procedure is followed for each grid cell and the modeling timestep separately, it leads to a very high computational effort at a continental scale and daily time steps. We then summarize the daily TTDs for each grid cell with statistical measures corresponding to the median (TT₅₀) and interquartile range (TT₁₀R), as well as the tails of the distribution such as the 10th (TT₁₀) and 90th (TT₉₀) percentile estimates.

One of the critical considerations of the above TTD formulation is the choice of the mixing schemes (or Storage Selection; SAS functions²⁵,⁶¹,⁶²,⁶⁹,⁷⁰) for the preference of water parcels with different ages being sampled by outflows. Among several mixing approaches, we consider here a random sampling scheme for deriving TTDs in each soil layer meaning that all water parcels in storage have an equal preference for sampling; consequently, all the following analyses of the TTD characterization are contingent on this selection. Our decision is motivated by the fact that there is no a-priori information available for the mixing schemes (or SAS functions) at a continental scale and that the random sampling scheme is one with the highest entropy. It is important to note that despite the random sampling scheme used here for characterizing the TTDs in individual soil layers, the overall sampling scheme for the entire soil column is far from being random due to difference in soil-water content and evapotranspiration fluxes in different soil compartments²⁶; and such an approach provides a meaningful way to simulate non-random sample dynamics as shown in a recent study⁷⁰.

**Nitrate vulnerability assessment.** Our nitrate vulnerability assessment is based on the leaching of excess nitrate from the root zone after accounting for the plant uptake (and other turnover processes). The mechanisms considered are downward transport and removal by denitrification, and we contrast the corresponding timescales of hydrologic transport (TT) with denitrification reaction (RT). Owing to the lack of reliable observations of RT, especially at large scales, we take a scenario approach and consider a wide range of RT estimates varying between 0.5 and 5 years. We represent the RT as a characteristic reaction timescale³⁸ corresponding to a given percentage of removal of initial substrate—here taken as 50% (RT₅₀) to allow for comparability with the corresponding hydrologic transport timescale (TT₅₀). Following the first-order kinetics adopted in this study, a RT₅₀ of 1 year, for example, would correspond to a denitrification rate constant of -ln(0.5)/RT₅₀ = 0.69 y⁻¹. Following previous large-scale studies³⁹,⁴²,⁴⁶, RT represents an effective timescale 〈RT₅₀〉 that encapsulates the relevant environmental factors, such as soil moisture, temperature, and organic carbon content, that affect the site-specific reaction behaviors⁴⁸.

We consider two cases for the nitrate vulnerability assessment that account for the static and transient behaviors of transport and reaction timescales. In the static case, we use the average estimates of the transport times µ(TT₅₀) and contrast them with the effective 〈RT₅₀〉 values. In the transient case, we account for the daily dynamics of TT₅₀(t) and RT₅₀(t) for each grid cell. The transient nature of the site-specific RT₅₀(t) is constructed from the spatiotemporal variability of the environmental factors, $f_s(t)$, that affect the daily variability of the denitrification process. Specifically, we account for the dynamic reduction factors caused by varying soil moisture $f_s(t)$ and temperature $f_r(t)$ conditions. Following established parameterization approaches³⁹⁻⁴¹,⁴³,⁴⁸, we define the following time-varying, dimensionless power functions $f_s(t)$ and $f_r(t)$ varying between [0,1] to reflect the status of anoxic and optimal temperature conditions required for the denitrification process⁴⁸:

$$f_s(t) = \begin{cases} 0 & S(t) < S_r \\ \left( \frac{S(t) - S_r}{S_r - S_r} \right)^\omega & S_r \le S(t) \le S_m \\ 1 & \text{otherwise,} \end{cases} \quad (3)$$

and,

$$f_r(t) = \beta^{\frac{T(t) - T_r}{T_r}} \quad (4)$$

where $S(t)$ is the soil-water storage at a given time, $t$, relative to the saturation limit ($S_m$), and $S_r$ is the threshold limit below which the denitrification process is completely inhibited, taken here as one-third of the saturation limit⁴¹,⁴³,⁴⁸. The parameter $\omega$ (= 2.5) defines the steepness of the curve⁴¹,⁴³,⁴⁸. $f_r(t)$ represents the effect of increasing ambient temperature $T(t)$ on denitrification. $T_r$ (= 25 °C) represents the reference temperature where $f_r(t) = 1$, and $\beta$ (= 2) is the factor that embeds the dependence of the denitrification rate on ambient temperature³⁹,⁴⁰,⁴⁸. We also consider the effect of spatially varying soil organic carbon content by accounting for the spatial information of potential denitrification rate constants based on the agro-ecosystem Carbon And Nitrogen Dynamics (CANDY) model⁷¹. Across the majority of the cultivated areas in Europe, the organic carbon content varies between 0.5 and 8%, with a median estimate of ~1% (see Supplementary Fig. 11) and the corresponding rates range between 0.02 and 0.16 kg/ha/day for the 1 dm soil layer. We use the relative (spatial) variations in this rate information to derive the spatial multiplier factors ($f_{\text{OC}}$) and we reference these factors to the nominal value of 1.0 for the base carbon content of 0%. Following an established approach⁴¹,⁴³,⁴⁸, we estimate the combined effect of multiple environmental factors for each grid cell and modeling timestep as $f_e(t) = f_s(t)f_r(t)f_{\text{OC}}$. However, this $f_e(t)$
value only describes the environmental condition on a single day t, whereas any solute entering on a given day $t_{in}$ will reside in the root zone for a much longer time period and be, therefore, exposed to a range of environmental conditions. We express this combined effect of the varying environmental conditions over the entire transport time as:

$$ \langle f_E \rangle (t_{in}) = \int_{t_{in}}^{\infty} p_I(t - t_{in} | t_{in}) \left( \frac{1}{t - t_{in}} \int_{t_{in}}^{t} f_E(t') dt' \right) dt. \quad (5) $$

The environmental conditions $f_E(t)$ in the above formulation are, therefore, weighed such that their different impact over the entire transport time is acknowledged. By virtue of the inner integral in the above equation, water, and therefore solutes (excess nitrate), leaving at any given day is only affected by the averaged behavior of the environmental condition up to that day (i.e., from $t_{in}$ to $t$). The combined behavior is then represented as an aggregated behavior of all the different portions of water parcels or solutes leaving the root zone; and this effect is reflected in the outer integral in the above equation. We use as grid-specific environmental factor the normalized, time-varying weighting factor $\langle f'_E \rangle(t)$, to establish the site-specific temporal dynamics of the denitrification timescale $RT_{50}(t)$. The normalization is based on the site-specific averaged $\langle f'_E \rangle(t)$ values that allows for the preservation of the effective $\langle RT_{50} \rangle$ across the study domain.

Our subsurface nitrate vulnerability analysis makes use of the dimensionless Damköhler number ($D_a = \frac{TT_{50}}{RT_{50}}$) to depict the complex interplay between the hydrologic transport (leaching) and biogeochemical turnover (denitrification) timescales$^{34–36,38}$. This dimensionless number ranges between 0 and $\infty$ with $D_a < 1$ indicate the dominance of the hydrologic transport over the reaction timescales, and vice versa. We use this objective measure to delineate the regions across Europe (with $D_a \le 1$) that are vulnerable to nitrate leaching from soil. We synthesize this information for the cultivated area using the cropland map$^{72}$ of the year 2000, which was compiled using extensive resources at national and sub-national level agricultural census statistics as well satellite-based land cover classification datasets$^{72}$. We consider as cultivated grid cells only those that have a cropland fraction of at-least 5%. Based on this threshold, a total number of nearly 8100 cropland cells at 0.25° grid resolution (leading to a total of ~192 million ha of cultivated areas) are considered in our vulnerability assessment analysis.

**Vulnerability assessment contingent on N-surplus.** Our analysis mainly considers the nitrate vulnerability across all cultivated areas without accounting for the spatial heterogeneity of N-surplus (or excess N) estimates, i.e., net nitrogen balance after accounting for fertilizer and atmospheric inputs and plant uptake. Although these estimates are needed to properly characterize the risk of nitrate leaching, they are rarely available at sub-annual timescales. Taking the example of the N-surplus balances for croplands around the year 2000$^{73}$, we note that not all of the cultivated areas across the study domain have high N-surplus, especially the east European regions bordering the Black sea (see Supplementary Fig. 12 for the N-surplus map and the corresponding cropland fractions in Fig. 3b). There is, however, a general tendency of increased N-surplus with the larger fraction of cropland areas (see Supplementary Fig. 12). We find that around 33% of the total cultivated areas have a high N surplus ($\ge 1 \times 10^5$ kg-N per grid cell) and is vulnerable to nitrate leaching under the static vulnerability assessment case (i.e., cropland cells of $D_a \le 1$ based on the averaged $TT_{50}$ and effective $\langle RT_{50} \rangle$ value). For the transient case based on the time-varying $TT_{50}$ and $RT_{50}$, we find around 58% of the total cultivated areas with high N-surplus is temporarily vulnerable ($D_a(t) \le 1$) to nitrate leaching for at least one-third of the year. The analysis conducted here provides a preliminary overview on the nitrate risk assessment, and further improvements in this aspect require more information on the temporal variability of N-surplus.

## Data availability

Data used in this study have been obtained from the following open sources: Terrain elevation GOTONOPO30 and EU-DEM from https://www.usgs.gov/ (https://doi.org/10.5066/f7df6pqs) and http://www.eea.europa.eu/data-and-maps/data/eu-dem; the river database CCM2 v2.1 from https://ccm.jrc.ec.europa.eu; Soils texture maps based in HWSD from https://webarchive.iiasa.ac.at/Research/LUC/External-World-soil-database/HTML/; the land cover product GlobCOVER v2 from http://due.esrin.esa.int/page_globcover.php; the CORINE land cover products (v18.4) from http://land.copernicus.eu/pan-european/corine-land-cover; the hydrogeology map IHME1500 v11 from http://www.bgr.bund.de/ihme1500; the historical forcings E-OBS v13 from https://www.ecad.eu/E-OBS/; and the 2000 cropland and N-excess maps from http://www.earthstat.org. Further details on these and other supporting databases are provided in Supplementary Table I. Finally, the underlying data for drawing the main conclusions of this study are available at http://dx.doi.org/10.5281/zenodo.3531907; and other auxiliary model simulations can be available from the corresponding author upon request.

## Code availability

The underlying model source code along with a test case to validate successful installation of the mesoscale Hydrologic Model (mHM) and the detail user manual are available at https://git.ufz.de/mhm/. The source code for computing the travel-time distributions (TTDs) is adapted after Hesse et al.$^{27}$, which is available at https://git.ufz.de/

mhm/mhm/tree/develop/post-proc/sas/. Further support on the details of processing algorithms can be obtained from the corresponding author upon request.

Received: 21 April 2020; Accepted: 9 November 2020;
Published online: 09 December 2020

## References

1. WFD. Directive 2000/60/ec of the European parliament and of the council of 23 october 2000 establishing a framework for Community action in the field of water policy' or, in short, the EU Water Framework Directive. Off. J. Euro. Commun. L 327, 1–72 (2000).

2. Kampa, E. et al. EEA 2018: European Waters Assessment of Status and Pressures. EEA Report No 7/2018 https://www.eea.europa.eu/publications/state-of-water (2018).

3. Rockström, J. et al. A safe operating space for humanity. Nature 461, 472 (2009).

4. Sutton, M. A. et al. Too much of a good thing. Nature 472, 159–161 (2011).

5. Erisman, J. W. et al. Consequences of human modification of the global nitrogen cycle. Philos. Trans. R. Soc. B-Biol. Sci. 368, 20130116 (2013).

6. Bodirsky, B. L. et al. Reactive nitrogen requirements to feed the world in 2050 and potential to mitigate nitrogen pollution. Nat. Commun. 5, 3858 (2014).

7. Zhang, X. et al. Managing nitrogen for sustainable development. Nature 528, 51–59 (2015).

8. Damania, R., Desbureaux, S., Rodella, A.-S., Russ, J. & Zaveri, E. Quality Unknown: The Invisible Water Crisis (2019). https://openknowledge.worldbank.org/handle/10986/32245.

9. EU. Council of the European Communities, 1991. Directive 91/676/EEC Concerning the Protection of Waters Against Pollution Caused by Nitrates From Agricultural Sources. 12 December 1991 (Council of the European Communities, Brussels 1991).

10. Kumar, P., Bansod, B. K., Debnath, S. K., Thakur, P. K. & Ghanshyam, C. Index-based groundwater vulnerability mapping models using hydrogeological settings: A critical evaluation. Environ. Impact Assess. Rev. 51, 38–49 (2015).

11. Wachniew, P. et al. Toward operational methods for the assessment of intrinsic groundwater vulnerability: A review. Crit. Rev. Environ. Sci. Technol. 46, 827–884 (2016).

12. Arauzo, M. Vulnerability of groundwater resources to nitrate pollution: A simple and effective procedure for delimiting nitrate vulnerable zones. Sci. Total Environ. 575, 799–812 (2017).

13. Machiwal, D., Jha, M. K., Singh, V. P. & Mohan, C. Assessment and mapping of groundwater vulnerability to pollution: Current status and challenges. Earth-Sci. Rev. 185, 901–927 (2018).

14. Ground Water Vulnerability Assessment: Predicting Relative Contamination. Published by By National Research Council, Division on Earth and Life Studies, Commission on Geosciences, Environment and Resources, Committee for Assessing Ground Water Vulnerability (National Academies Press, 1993).

15. McGuire, K. I. et al. The role of topography on catchment-scale water residence time. Water Resour. Res. 41, W05002 (2005).

16. Tetzlaff, D. et al. How does landscape structure influence catchment transit time across different geomorphic provinces? Hydrol. Process. 23, 945–953 (2009).

17. Hrachowitz, M., Soulsby, C., Tetzlaff, D., Dawson, J. J. C. & Malcolm, I. A. Regionalization of transit time estimates in montane catchments by integrating landscape controls., Water Resour. Res. 45, W05421 (2009).

18. Hrachowitz, M., Soulsby, C., Tetzlaff, D. & Speed, M. Catchment transit times and landscape controls—does scale matter? Hydrol. Process. 24, 117–125 (2010).

19. Musolff, A., Fleckenstein, J. H., Rao, P. S. C. & Jawitz, J. W. Emergent archetypic patterns of coupled hydrologic and biogeochemical responses in catchments. Geophys. Res. Lett. 44, 4143–4151 (2017).

20. Sprenger, M. et al. The demographics of water: A review of water ages in the critical zone. Rev. Geophys. 57, 800–834 (2019).

21. Botter, G., Bertuzzo, E. & Rinaldo, A. Transport in the hydrologic response: Travel time distributions, soil moisture dynamics, and the old water paradox. Water Resour. Res. 46, W03514 (2010).

22. van der Velde, Y., de Rooij, G. H., Rozemeijer, J. C., van Geer, F. C. & Broers, H. P. Nitrate response of a lowland catchment: On the relation between stream concentration and travel time distribution dynamics. Water Resour. Res. 46, W11534 (2010).

23. Benettin P., van der Velde Y., Sjoerd E. A., van der Zee T. M., Rinaldo A. & Botter, G. Chloride circulation in a lowland catchment and the formulation of transport by travel time distributions. Water Resour. Res. 49, 4619–4632 (2013).
24. Hrachowitz, M., Savenije, H., Bogaard, T. A., Tetzlaff, D. & Soulsby, C. What can flux tracking teach us about water age distribution patterns and their temporal dynamics? *Hydrol. Earth Syst. Sci.* **17**, 533–564 (2013).

25. Harman, C. J. Time-variable transit time distributions and transport: Theory and application to storage-dependent transport of chloride in a watershed. *Water Resour. Res.* **51**, 1–30 (2015).

26. Benettin, P., Kirchner, J. W., Rinaldo, A. & Botter, G. Modeling chloride transport using travel time distributions at Plynlimon, Wales. *Water Resour. Res.* **51**, 3259–3276 (2015).

27. Heße, F., Zink, M., Kumar, R., Samaniego, L. & Attinger, S. Spatially distributed characterization of soil-moisture dynamics using travel-time distributions. *Hydrol. Earth Syst. Sci.* **21**, 549–570 (2017).

28. Kuffour, B. N. O. et al. Simulating coupled surface–subsurface flows with ParFlow v3.5.0: Capabilities, applications, and ongoing development of an open-source, massively parallel, integrated hydrologic model. *Geoscientific Model Dev.* **13**, 1373–1397 (2020).

29. Harman, C. J. et al. Climate, soil, and vegetation controls on the temporal variability of vadose zone transport. *Water Resour. Res.* **47** (2011).

30. Basu, N. B., Thompson, S. E. & Rao, P. S. C. Hydrologic and biogeochemical functioning of intensively managed catchments: A synthesis of top-down analyses. *Water Resour. Res.* **47**, W00J15 (2011).

31. Ascott, M. et al. Global patterns of nitrate storage in the vadose zone. *Nat. Commun.* **8**, 1416 (2017).

32. Fan, Y., Miguez-Macho, G., Jobbágy, E. G., Jackson, R. B. & Otero-Casal, C. Hydrologic regulation of plant rooting depth. *Proc. Natl Acad. Sci.* **114**, 10572–10577 (2017).

33. Fan, J., McConkey, B., Wang, H. & Janzen, H. Root distribution by depth for temperate agricultural crops. *Field Crops Res.* **189**, 68–74 (2016).

34. Ocampo, C. J., Oldham, C. E. & Sivapalan, M. Nitrate attenuation in agricultural catchments: Shifting balances between transport and reaction. *Water Resour. Res.* **42**, W01408 (2006).

35. Oldham, C. E., Farrow, D. E. & Peiffer, S. A generalized Damkoehler number for classifying material processing in hydrological systems. *Hydrol. Earth Syst. Sci.* **17**, 1133–1148 (2013).

36. Pinay, G. et al. Upscaling nitrogen removal capacity from local hotspots to low stream orders' drainage basins. *Ecosystems* **18**, 1101–1120 (2015).

37. Botter, G. Catchment mixing processes and travel time distributions. *Water Resour. Res.* **48**, W05545 (2012).

38. Pittroff, M., Frei, S. & Gilfedder, B. S. Quantifying nitrate and oxygen reduction rates in the hyporheic zone using 222rn to upscale biogeochemical turnover in rivers. *Water Resour. Res.* **53**, 563–579 (2017).

39. Lin, B.-L., Sakoda, A., Shibasaki, R., Goto, N. & Suzuki, M. Modelling a global biogeochemical nitrogen cycle in terrestrial ecosystems. *Ecol. Model.* **135**, 89–110 (2000).

40. Botter, G., Settin, T., Marani, M. & Rinaldo, A. A stochastic model of nitrate transport and cycling at basin scale. *Water Resour. Res.* **42**, W04415 (2006).

41. Lindstroem, G., Pers, C., Rosberg, J., Strömqvist, J. & Arheimer, B. Development and testing of the HYPE (hydrological predictions for the environment) water quality model for different spatial scales. *Hydrol. Res.* **41**, 295–319 (2010).

42. Van Meter, K. J., Basu, N. B. & Van Cappellen, P. Two centuries of nitrogen dynamics: Legacy sources and sinks in the Mississippi and Susquehanna river basins. *Glob. Biogeochem. Cycles* **31**, 2–23 (2017).

43. Yang, X. et al. A new fully distributed model of nitrate transport and removal at catchment scale. *Water Resour. Res.* **54**, 5856–5877 (2018).

44. Van Meter, K. J., Van Cappellen, P. & Basu, N. B. Legacy nitrogen may prevent achievement of water quality goals in the Gulf of Mexico. *Science* **360**, 427–430 (2018).

45. Green, M. et al. Effective denitrification scales predictably with water residence time across diverse systems. *Nat. Preced.* **4**, https://doi.org/10.1038/ npre.2009.3520.1 (2009).

46. Goll, D. S. et al. Carbon-nitrogen interactions in idealized simulations with JSBACH (version 3.10). *Geoscientific Model Dev.* **10**, 2009–2030 (2017).

47. EC. Directorate-General for Environment (European Commission). *Recommendations for Establishing Action Programmes Under Directive 91/676/ EEC Concerning the Protection of Waters Against Pollution Caused by Nitrates from Agricultural Sources. Part A. Appendix 2, Maps of pedo-climatic zones in Europe.* (2011). https://publications.europa.eu/s/m14g.

48. Heinen, M. Simplified denitrification models: Overview and properties. *Geoderma* **133**, 444–463 (2006).

49. Sinha, E., Michalak, A. & Balaji, V. Eutrophication will increase during the 21st century as a result of precipitation changes. *Science* **357**, 405–408 (2017).

50. Samaniego, L. et al. Anthropogenic warming exacerbates European soil moisture droughts. *Nat. Clim. Change* **8**, 421 (2018).

51. Bowles, T. M. et al. Addressing agricultural nitrogen losses in a changing climate. *Nat. Sustainability* **1**, 399–408 (2018).

52. Samaniego, L., Kumar, R. & Attinger, S. Multiscale parameter regionalization of a grid-based hydrologic model at the mesoscale. *Water Resour. Res.* **46**, W05523 (2010).

53. Kumar, R., Samaniego, L. & Attinger, S. Implications of distributed hydrologic model parameterization on water fluxes at multiple scales and locations. *Water Resour. Res.* **49**, 360–379 (2013).

54. Rakovec, O. et al. Multiscale and multivariate evaluation of water fluxes and states over European river basins. *J. Hydrometeorol.* **17**, 287–307 (2016).

55. Bergström, S. in *Computer Models of Watershed Hydrology* (ed. Singh, V. P.) Computer Models of Watershed Hydrology, 443–476 (Water Resources Publications, 1995).

56. Liang, X., Lettenmaier, D. P., Wood, E. F. & Burges, S. J. A simple hydrologically based model of land surface water and energy fluxes for general circulation models. *J. Geophys. Res. Atmos.* **99**, 14415–14428 (1994).

57. Thober, S. et al. The multiscale routing model mRM v1.0: Simple river routing at resolutions from 1 to 50 km. *Geoscientific Model Dev.* **12**, 2501–2521 (2019).

58. Samaniego, L. et al. Hydrological forecasts and projections for improved decision-making in the water sector in Europe. *Bull. Am. Meteorological Soc.* **100**, 2451–2472 (2019).

59. Haylock, M. R. et al. A European daily high-resolution gridded data set of surface temperature and precipitation for 1950–2006, *J. Geophys. Res.* **113**, D20119 (2008).

60. Yang, X., Jomaa, S. & Rode, M. Sensitivity analysis of fully distributed parameterization reveals insights into heterogeneous catchment responses for water quality modeling. *Water Resour. Res.* **55**, 10935–10953 (2019).

61. Botter, G., Bertuzzo, E. & Rinaldo, A. Catchment residence and travel time distributions: The master equation. *Geophys. Res. Lett.* **38**, L11403 (2011).

62. van der Velde, Y., Torfs, P., Van der Zee, S. & Uijlenhoet, R. Quantifying catchment-scale mixing and its effect on time-varying travel time distributions. *Water Resour. Res.* **48** (2012).

63. Benettin, P. et al. Using SAS functions and high-resolution isotope data to unravel travel time distributions in headwater catchments. *Water Resour. Res.* **53**, 1864–1878 (2017).

64. Cornaton, F. J., Park, Y.-J., Normani, S. D., Sudicky, E. A. & Sykes, J. F. Use of groundwater lifetime expectancy for the performance assessment of a deep geologic waste repository: 1. Theory, illustrations, and implications, *Water Resour. Res.* **44**, W04406 (2008).

65. Benettin P., Rinaldo, A. & Botter, G. Tracking residence times in hydrological systems: Forward and backward formulations. *Hydrol. Process.* **29**, 5203–5213 (2015).

66. Niemi, A. J. Residence time distributions of variable flow processes. *Int. J. Appl. Radiat. Isotopes* **28**, 855–860 (1977).

67. Bertuzzo, E., Thomet, M., Botter, G. & Rinaldo, A. Catchment-scale herbicides transport: Theory and application. *Adv. Water Resour.* **52**, 232–242 (2013).

68. Benettin, P. et al. Linking water age and solute dynamics in streamflow at the Hubbard Brook Experimental Forest, NH, USA. *Water Resour. Res.* **51**, 9256–9272 (2015).

69. Rinaldo, A. et al. Storage selection functions: A coherent framework for quantifying how catchments store and release water and solutes. *Water Resour. Res.* **51**, 4840–4847 (2015).

70. Remondi, F., Kirchner, J. W., Burlando, P. & Fatichi, S. Water flux tracking with a distributed hydrological model to quantify controls on the spatio-temporal variability of transit time distributions. *Water Resour. Res.* **54**, 3081–3099 (2018).

71. Franko, U., Oelschloegel, B. & Schenk, S. Simulation of temperature-, water- and nitrogen dynamics using the model CANDY. *Ecol. Model.* **81**, 213–222 (1995).

72. Ramankutty, N., Evan, A. T., Monfreda, C. & Foley, J. A Farming the planet: 1. Geographic distribution of global agricultural lands in the year 2000. *Global Biogeochem. Cycles.* **22**, GB1003 (2008).

73. West, P. C. et al. Leverage points for improving global food security and the environment. *Science* **345**, 325–328 (2014).

## Acknowledgements

R.K. acknowledge the support from the Initiative and Networking Fund of the Helmholtz Association through the project Advanced Earth System Modeling Capacity (ESM) (www.esm-project.net). P.S.C.R. acknowledges partial financial support from the Lee A. Reith Endowment in the Lyles School of Civil Engineering, Purdue University. Collaborations between the Helmholtz Center for Environmental Research-UFZ, the Purdue University, and the University of Florida are supported by an international Memorandum of Understanding (MOU). Support to F.S. was provided by the Reduced Complexity Models project co-funded by the Helmholtz Association. We also acknowledge the partial funding for collection and processing of datasets under a contract for the Copernicus Climate Change Service (C3S 441 Lot1 NERC; http://edge.climate.copernicus.eu). ECMWF implements this service and the Copernicus Atmosphere Monitoring Service on behalf of the European Commission. We would like to thank people from various organizations and projects for kindly providing us with the data that
were used in this study, which includes among others: USGS, ESA, JRC, NASA, GRDC, BGR, EARTHSTAT, and IIASA. We acknowledge the E-OBS dataset from the EU-FP6 project ENSEMBLES (http://ensembles-eu.metoffice.com) and the data providers in the ECA&D project (http://www.ecad.eu). Data analysis was conducted at the High-Performance Computing (HPC) Cluster EVE, a joint effort of both the Helmholtz Center for Environmental Research-UFZ and the German Center for Integrative Biodiversity Research (iDiv) Halle-Jena-Leipzig.

**Author contributions**

R.K. conceived the idea and designed the study with inputs from F.H., P.S.C.R., A.M., J.W.J., and F.S. R.K. collected required datasets, conducted model simulations, and performed the analysis. R.K. led the writing with F.H. and P.S.C.R. F.H., P.S.C.R., A.M., J.W.J., F.S., L.S., J.H.F., O.R., S.T., and S.A. contributed to the editing and commenting on the manuscript.

**Funding**

Open Access funding enabled and organized by Projekt DEAL.

**Competing interests**

The authors declare no competing interests.

**Additional information**

Supplementary information is available for this paper at https://doi.org/10.1038/s41467-020-19955-8.

**Correspondence** and requests for materials should be addressed to R.K.

**Peer review information** *Nature Communications* thanks Paolo Benettin and other, anonymous reviewers for their contributions to the peer review of this work. Peer review reports are available.

**Reprints and permission information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2020