==== Front Nat Commun Nat Commun Nature Communications 2041-1723 Nature Publishing Group UK London 19955 10.1038/s41467-020-19955-8 Article Strong hydroclimatic controls on vulnerability to subsurface nitrate contamination across Europe http://orcid.org/0000-0002-4396-2037Kumar R. rohini.kumar@ufz.de 1 Heße F. 1 Rao P. S. C. 23 http://orcid.org/0000-0002-0115-1359Musolff A. 1 http://orcid.org/0000-0002-6745-0765Jawitz J. W. 4 Sarrazin F. 1 http://orcid.org/0000-0002-8449-4428Samaniego L. 1 Fleckenstein J. H. 15 http://orcid.org/0000-0003-2451-3305Rakovec O. 16 Thober S. 1 Attinger S. 17 1 grid.7492.80000 0004 0492 3830UFZ-Helmholtz Centre for Environmental Research, Leipzig, Germany 2 grid.169077.e0000 0004 1937 2197Lyles School of Civil Engineering, Purdue University, West Lafayette, IN USA 3 grid.169077.e0000 0004 1937 2197Agronomy Department, Purdue University, West Lafayette, IN USA 4 grid.15276.370000 0004 1936 8091Soil and Water Sciences Department, University of Florida, Gainesville, FL USA 5 grid.7384.80000 0004 0467 6972Bayreuth Center of Ecology and Environmental Research, University of Bayreuth, Leipzig, Germany 6 grid.15866.3c0000 0001 2238 631XFaculty of Environmental Sciences, Czech University of Life Sciences, Prague, Czech Republic 7 grid.11348.3f0000 0001 0942 1117University of Potsdam, Potsdam, Germany 9 12 2020 9 12 2020 2020 11 630221 4 2020 9 11 2020 © The Author(s) 2020Open 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/.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. Excess fertilizer use causes subsurface contamination. Here, the authors conduct an assessment of water quality vulnerability across Europe, finding that 75% of agricultural regions are susceptible to nitrate contamination for least one-third of the year, two times more than using standard estimation procedure. Subject terms Element cyclesHydrologyissue-copyright-statement© The Author(s) 2020 ==== Body Introduction Despite >15 years of water quality protection implementation under the EU Water Framework Directive (EU-WFD1), the most recent EU-WFD report2 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 concern3–8. Consequently, the European Nitrate Directive9—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 zones10–13. 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 pollutants11,14. 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 compartment15–20. Much work has been based on steady-state TTDs, however, more recently, studies have started to acknowledge the transient nature of TTDs21–27. 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-scales28, 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 dynamics21,22,26,27 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 waters29–31. The rooting depth varies across European landscapes depending on, among other geophysical attributes, vegetation types and groundwater table depth32. 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 lands33. We use the dimensionless Damköhler number34–36 to link the hydrologic and biogeochemical timescales (see “Methods”) and provide an objective measure for the 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 (ϕ = 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 (TTIQR) increases jointly with the median travel-time (TT50) throughout Europe (Fig. 1d–f and see Supplementary Fig. 1). The humid location (UK; Fig. 1a–d) exhibits marked seasonality with higher TT50 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 TT50 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 datasets29,37.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 TT50 and interquartile ranges TTIQR (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. Synthesis of hydrologic transport times across Europe   Figure 2a, b summarizes the spatial patterns of the temporal mean μ(TT50) and standard deviation σ(TT50) 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 μ(TT50) and 52 and 320 days for σ(TT50). Fifty percent of the study domain has μ(TT50) and σ(TT50) values that exceed 365 and 120 days, respectively. Regions with the shortest soil travel times (μ(TT50) ≤  180 days) are located in areas of high and frequent rainfall (e.g., alpine, northern UK and northern Spain— Pyrenees). The longest travel times (μ(TT50) ≥ 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 between μ(TT50) and σ(TT50). This was confirmed through a linear regression analysis (Fig. 2d), resulting in R2 = 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 TT50, which is the regression slope between μ(TT50) and σ(TT50), is remarkably consistent across Europe (~0.4). An analogus consistency in CV values has been also reported in a previous study37 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 TT50 CV = 0.4 is in good agreement with the temporal variability of daily potential evapotranspiration (CV = 0.34, see Supplementary Fig. 2).Fig. 2 Synthesis of the hydrologic transport time dynamics across Europe. The transient feature of the daily median travel times (TT50; blue lines in Fig. 1d, e, f) is summarized as their temporal mean, μ(TT50) and standard deviation, σ(TT50) (a, b) for the period 1985–2015. The strong spatial correspondence between μ(TT50) and σ(TT50) 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 (Ep¯) and mean precipitation (P¯), is identified as the dominant factor controlling the spatial heterogeneity of the transient transport characteristics as μ(TT50) and σ(TT50) 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 rather well-organized spatial patterns of μ(TT50) and σ(TT50) follow the hydroclimatic gradient observed across Europe (Fig. 2c), and here, the latter pattern is represented through the aridity index (ϕ) 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 μ(TT50) and σ(TT50) can be explained solely by the spatial heterogeneity of ϕ (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 μ(TT50) and σ(TT50) to ϕ 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 column29. 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 μ(TT50) and σ(TT50) compared to ϕ 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 ϕ 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 ϕ interval of 0.15). For the (average) binned estimates, we find an almost perfect linear correlation of ϕ with μ(TT50) and σ(TT50) (R2 = 0.96–0.98; p-value <  0.00001). These results are also found for the extremes of the daily TT distributions. For example, TT10 and TT90, which are indicative of young and old water fractions, respectively, also show a high spatial correlation with ϕ (see Supplementary Figs. 5 and 6). This result underpins the dominant role of the hydroclimatic factor (ϕ) 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 vulnerability14 of 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 timescales38, RT50, defined here as the 50% removal of the initial substrate, to allow comparability with the median transport times (TT50). We consider a range of RT50 values between 0.5 and 5 years in our analysis39–43 (see “Methods” for more details). RT50 represents the effective timescale encapsulating relevant environmental factors, such as soil moisture, temperature, and organic carbon content, that affect the site-specific reaction rates42,44. 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 Da=TT50RT50) that enables us to assess the interplay between these two competing processes across the geographical domain34–36,38. When Da < 1, transport (leaching) dominates over reaction (denitrification), and vice versa. Our static vulnerability assessment case, based on the range of RT50 (0.5–5 years) and the averaged transport times μ(TT50), results in Da values ranging between 0.05 and 4.0 across the majority of European cultivated areas (Fig. 3a). To interpret this Da number, we rely on a series of prior studies. First, a prior study34 presents a variety of field observational datasets that demonstrate an empirical (non-linear) relationship between Da and nitrate removal. Subsequent studies36,45 showed how this empirical relationship can be described by a parametric model based on the exponential function. Using the latter, the above value of Da = 0.05 would imply more than 90% of the nitrate leaching from soil, while a value of Da = 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 (Da ≤ 1) under the higher RT50 value (≥ 2 years), while very few locations would be classified as vulnerable in the case of the lowest RT50 value of 0.5 year due to the dominance of the denitrification timescale.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 (Da=TT50RT50) for a range of effective denitrification timescales RT50 and the averaged transport times μ(TT50) estimated across the cultivated areas of Europe72. Boxplot is displayed with the horizontal bar at the median, the box indicates the first and third quartiles and the whiskers indicate  ±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 Da estimates for the static vulnerability assessment case corresponding to the effective reaction timescale, RT50 and the average transport timescales, μ(TT50). Areas with Da ≤ 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 Da(t) estimates for every cultivated grid cell, arranged according to their corresponding aridity index values (ϕ). Summary statistics are presented as quantiles of the daily Da 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 Da(t) ≤ 1 estimated based on the time-varying TT50(t) and RT50(t). Considering the transient nature of Da, nearly 75% of cultivated areas of Europe become vulnerable to nitrate leaching (Da(t) ≤ 1) for at least one-third of the year (i.e., frequency estimates ≥  0.33 or 4/12). We further elaborate on the moderate case of an RT50 = 1 year, which is a highly plausible value, as inferred by the values of denitrification timescales in previous large-scale39,42,46 and catchment-scale studies40,43. Using this RT50 value and the average transport time μ(TT50), we find ~42% of the cultivated areas across Europe to be vulnerable to subsurface nitrate leaching, with Da ≤ 1 (Fig. 3b). Figure 3b shows that the cultivated areas in the humid and transitional climate zones (ϕ ≤ 1.5), specifically in the western part of Europe, i.e., France, Germany, Italy, UK and Ireland, are dominated by hydrologic transport (Da ≤ 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 Da map are remarkably consistent with the (nitrate) leaching risk potential map published by the European Commission47 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 TT50(t) and RT50(t). Here, we account for the spatiotemporal variability of the environmental factors, fE(t), that affect the daily dynamics of the denitrification timescale, RT50(t). We include the effect of varying soil moisture, air temperature, and organic carbon content following established parameterization schemes39–41,43,48 (see “Methods” for more details). This allowed us to construct the temporal variability of RT50(t) for each grid cell by explicitly accounting for the space/time variability of the environmental factor (fE) (see “Methods”). The daily mean and variability of the environmental factor (fE) follows the hydroclimatic gradient of the aridity index ϕ observed across Europe (see Supplementary Fig. 8). We observe lower mean and variability of the daily fE 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 fE because of wetter and more strongly seasonal soil moisture dynamics (Supplementary Fig. S8). Next we analyze the variability of the daily Da(t) that depicts the interplay between the hydrologic transport timescale (TT50) and the reactive timescale (RT50) behaviors for the excess nitrate removal (leaching vs. denitrification) through the soil. Here, in Fig. 3c we summarize the daily Da(t) as the quantile estimates for every cropland cell arranged according to their respective aridity index (ϕ) value. For example, for a given cropland cell located in a climate region of ϕ = 1 and having the 50th percentile (or median; Q50) value of the daily Da = 1.0, this number would indicate that the given cell would be prone to excess nitrate leaching (i.e., Da ≤ 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 Da 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 Da range (Q95-Q5) for cropland cells with ϕ of 1 to 3. Interestingly, while the averages (medians) of the daily Da values are usually high (above 1) in (semi-) arid regions, the Da 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 Da statistics (median and Q95-Q5) and the aridity index (ϕ) across the European croplands (R2 ≥ 0.92). We now analyze how the temporal dynamics of the daily Da(t) affects our vulnerability assessment. The results of the frequency analysis based on the time-varying Da ≤ 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 approach47, 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 conditions49–51. 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 (mHM52,53) 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 locations52–54. 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 models55,56. 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 studies52,53,57, 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 studies53,54,58 (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)59, the multiscale parameterization approach implemented in mHM allows for the explicit treatment of fine-scale, sub-grid variability of landscape features52,53. 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 datasets27 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 model43,60. 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)21,25,26,61–63. Specifically, we adopt the notion by Botter et al.,21,61 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 J (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 follows21: 1 pI(t−tin∣tin)=I(t)θ(tin)S(t)exp−∫tintI(t′)+E(t′)S(t′)dt′ with t − tin (t > tin) represents the time from the moment the water parcel enters the control volume (tin) until now (t). The partition function θ(tin) indicates the portion of the water parcel that enters the control volume at tin and leaves as the exfiltration flux (specifically, I as opposed to E) and is expressed as follows: 2 θ(tin)=∫tin∞I(t)S(t)exp−∫tintI(t′)+E(t′)S(t′)dt′dt. This partition function θ 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 time64,65. The complementary approach in which the ages of the water particles exiting the system are tracked back in time relates to an age concept65, and both forward and backward TTs can be related through the Niemi’s continuity equation66. Recent studies22,23,25,63,67,68 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 (TT50) and interquartile range (TTIQR), as well as the tails of the distribution such as the 10th (TT10) and 90th (TT90) percentile estimates. One of the critical considerations of the above TTD formulation is the choice of the mixing schemes (or StorAge Selection; SAS functions25,61,62,69,70) 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 compartments26; and such an approach provides a meaningful way to simulate non-random sample dynamics as shown in a recent study70. 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 timescale38 corresponding to a given percentage of removal of initial substrate—here taken as 50% (RT50) to allow for comparability with the corresponding hydrologic transport timescale (TT50). Following the first-order kinetics adopted in this study, a RT50 of 1 year, for example, would correspond to a denitrification rate constant of -ln(0.5)/RT50 = 0.69 y−1. Following previous large-scale studies39,42,46, RT represents an effective timescale RT50 that encapsulates the relevant environmental factors, such as soil moisture, temperature, and organic carbon content, that affect the site-specific reaction behaviors48. 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 μ(TT50) and contrast them with the effective RT50 values. In the transient case, we account for the daily dynamics of TT50(t) and RT50(t) for each grid cell. The transient nature of the site-specific RT50(t) is constructed from for the spatiotemporal variability of the environmental factors, fe(t), that affect the daily variability of the denitrification process. Specifically, we account for the dynamic reduction factors caused by varying soil moisture fS(t) and temperature fT(t) conditions. Following established parameterization approaches39–41,43,48, we define the following time-varying, dimensionless power functions fS(t) and fT(t) varying between [0,1] to reflect the status of anoxic and optimal temperature conditions required for the denitrification process48: 3 fS(t)=0S(t)