Modeling Nitrate Export From a Mesoscale Catchment Using StorAge Selection Functions

StorAge Selection (SAS) functions describe how catchments selectively remove water of different ages in storage via discharge, thus controlling the transit time distribution (TTD) and solute composition of discharge. SAS‐based models have been emerging as promising tools for quantifying catchment‐scale solute export, providing a coherent framework for describing both velocity‐driven and celerity‐driven transport. Due to their application in headwaters only, the spatial heterogeneity of catchment physiographic characteristics, land use management practices, and large‐scale validation have not been adequately addressed with SAS‐based models. Here, we integrated SAS functions into the grid‐based mHM‐Nitrate model (mesoscale Hydrological Model) at both grid scale (distributed model) and catchment scale (lumped model). The proposed model provides a spatially distributed representation of nitrogen dynamics within the soil zone and a unified approach for representing both velocity‐driven and celerity‐driven subsurface transport below the soil zone. The model was tested in a heterogeneous mesoscale catchment. Simulated results show a strong spatial heterogeneity in nitrogen dynamics within the soil zone, highlighting the necessity of a spatially explicit approach for describing near‐surface nitrogen processing. The lumped model could well capture instream nitrate concentration dynamics and the concentration–discharge relationship at the catchment outlet. In addition, the model could provide insights into the relations between subsurface storage, mixing scheme, solute export, and the TTDs of discharge. The distributed model shows results that are comparable to the lumped model. Overall, the results reveal the potential for large‐scale applications of SAS‐based transport models, contributing to the understanding of water quality‐related issues in agricultural landscapes.

to rainfall inputs could be several orders of magnitude faster than the response time of instream solute concentration to solute inputs (Hrachowitz et al., 2016). For these catchments, transport models explicitly based on the transit time are required to capture velocity-driven transport phenomena.
TTDs could be incorporated into conceptual catchment-scale transport models using different approaches. For example, the functional form of TTDs could be predefined and assumed to be time-invariant (Ilampooranan et al., 2019;Rinaldo et al., 2006;Van Meter et al., 2017. However, experimental data and numerical studies have indicated that TTDs (e.g., for discharge) are time-variant for many hydrological systems (Botter et al., 2010(Botter et al., , 2011Heibüchel et al., 2020;Kaandorp et al., 2018;Kim et al., 2016;Rodriguez et al., 2018;Van der Velde et al., 2012;J. Yang, Heidbüchel, et al., 2018). To capture the time-invariance and ease the problem of parameterization of TTDs, Van der Velde et al. (2010) used TTDs obtained from forward physically based groundwater models with particle tracking. In large catchments, however, a rigorous mathematical framework for the representation and parameterization of time-variant TTDs is needed. This is because the application of physically based groundwater models in large catchments is not always possible due to a lack of data and/or computational capacity.
Recent studies have shifted the focus from TTDs to StorAge Selection (SAS) functions (Harman, 2015(Harman, , 2019Rinaldo et al., 2015;Van der Velde et al., 2012). SAS functions describe how water parcels of different ages in storage are mixed to produce outflows, thus, affecting the dynamic of TTDs. SAS functions can be considered as a generalization of TTDs (Harman, 2019). SAS functions are more stable in time and easier for parameterization than TTDs ( Van der Velde et al., 2012). SAS functions could be combined with storage-discharge functions to provide a coherent framework for describing both celerity-driven water flow dynamics and velocity-driven solute transport mechanisms (Harman, 2019;Hrachowitz et al., 2016). For this purpose, there have been several SAS-based models developed for modeling solute (or isotope) transport from plot to catchment scales (e.g., Benettin et al., 2013Benettin et al., , 2017Benettin, Kirchner, et al., 2015;Bertuzzo et al., 2013;Harman, 2015;Lutz et al., 2017;Queloz et al., 2015;Wilusz et al., 2017). An in-depth discussion of SAS-based models was provided by Hrachowitz et al. (2016). In general, these studies have proven the effectiveness of the chosen models in capturing catchment-scale flow and transport phenomena for small catchments (with an area less than 10 km 2 ). However, application and validation of SAS-based models have not been done at larger spatial scales (e.g., drainage areas of about 100 km 2 ).
At larger scales, the catchment's landscape, meteorological conditions, and land use management practices are often heterogeneous. As a result, the catchment responses, especially the nutrient processes within the root zone, could be highly heterogeneous (e.g., X. Yang et al., 2019). While the SAS concept implicitly represents the heterogeneity in flow pathways, the spatial heterogeneity of biogeochemical processes, catchment characteristics, and meteorological conditions have not yet been adequately addressed within the SAS concept. In other words, effects of spatial heterogeneity and a thorough testing of the concept for larger scales have not yet been addressed. The main focus of this research is to fill those gaps.
The mesoscale Hydrologic Model (mHM)-Nitrate model (X. Yang, Jomaa, et al., 2018) is a grid-based water quality (nitrate) model with the hydrological and water quality concepts taken from two widely used models, the mHM (Kumar et al., 2013;Samaniego et al., 2010) and the HYdrological Predictions for the Environment model (Lindström et al., 2010). The mHM-Nitrate model could provide valuable insights into the spatial variability of water and nitrate dynamics in the root zone (X. Yang et al., 2019;X. Yang, Jomaa, et al., 2018). The model is able to account for the spatial heterogeneity in land used management practices, soil type, and meteorological forcing explicitly.
In the mHM-Nitrate model, the groundwater storage consists of active and inactive storages with an assumption of complete mixing between these storages. From our point of view, this assumption does not properly represent velocity-driven transport. Under the complete mixing assumption, a part of the solute input to groundwater would be instantaneously transported to the stream. In other words, there is no time lag between input and output signals. In catchments with velocity-driven transport, however, the time lag between input and output signals could be up to decades (Ehrhardt et al., 2019;Meals et al., 2010). In addition, the subsurface could be far from a completely mixed storage compartment (e.g., J. Yang, Heidbüchel, et al., 2018). Even if the assumption of complete mixing between the active and inactive storages is valid, the applicability of this assumption with nonconservative solutes is still unknown, particularly with application in a single compartment system with high temporal resolution. Few studies (e.g., Hrachowitz et al., 2013Hrachowitz et al., , 2015 have successfully employed either the well-mixed or the partial mixing assumption between the active and the inactive storages to capture velocity-driven transport of conservative solutes (Hrachowitz et al., 2016). Despite nitrate being a reactive solute, leached nitrate out of the soil zone in the mHM-Nitrate model is treated as a nonreactive solute. However, nitrate leaching out of the soil zone will be subject to additional removal processes along its flow path to the stream, for example, via denitrification in the shallow and deep aquifers (e.g., Fukada et al., 2003;Hiscock et al., 1991;Knoll et al., 2020;Kolbe et al., 2019;Rivett et al., 2008;Smith et al., 2004). Nevertheless, the existing mHM-Nitrate model (X. Yang, Jomaa, et al., 2018) provided a promising tool for further development using the SAS concept.
The objectives of this study are to (1) replace the description of the nitrate submodel for the subsurface (below the soil zone) in the mHM-Nitrate model with a time-variant SAS-based model, and by that (2) present a first test of a SAS-based transport model at a mesoscale catchment. The proposed model, hereinafter referred to as the mHM-SAS model, provides a unified approach for modeling both celerity-driven and velocity-driven transport at the catchment scale based on SAS functions. The model accounts for nitrate losses along its flow path from the bottom of the soil zone to the catchment outlet. In this study, we provide not only a detailed implementation of the SAS-based concept at catchment scale (lumped approach) but also an insight into the potential application of the SAS-based concept at a spatially more resolved grid scale (distributed approach).

The mHM-Nitrate Model
The mHM-Nitrate model is a grid-based hydrological and water quality (nitrate) model (Kumar et al., 2013;Samaniego et al., 2010;X. Yang, Jomaa, et al., 2018). Each grid cell consists of a series of leaky storage reservoirs, representing water storage in the soil zone, unsaturated zone, and saturated (groundwater) zone (Figure 1(a)). The soil zone has a depth of around 2 m, representing the root zone (the terms "soil zone" and "root zone" are used interchangeably in this study). The soil zone consists of three soil layers and the saturated zone is divided into active and inactive groundwater storages. The mHM model is parameterized using the Multiscale Parameter Regionalization approach to account for the subgrid variability of catchment properties and to avoid overparameterization (Samaniego et al., 2010).
The mHM-Nitrate model allows a spatially explicit representation of agricultural management practices (e.g., crop rotation, fertilizer application). Within the soil zone, the model tracks the fate of nitrogen in different pools: dissolved inorganic nitrogen (DIN), dissolved organic nitrogen (DON), active organic nitrogen (SON A ), and inactive organic nitrogen (SON I ) (Figure 1(c)). The model assumes that nitrate-nitrogen (N-NO 3 ) is equivalent to DIN. Nitrogen can be transformed between different nitrogen pools via mineralization, dissolution, and degradation within the soil zone. Only DIN and DON are transported with water to the unsaturated, saturated zone, and eventually to the stream. The mHM-Nitrate model does not consider the transformation from DON to DIN and denitrification occurring below the soil zone. In the saturated zone (groundwater), the active and inactive groundwater storages are assumed to be well mixed. The inactive groundwater storage, whose storage volume is set to be land use dependent, is assumed to be much larger than the active storage. It should be noted that the inactive groundwater storage did not exist in the original mHM model and was introduced by X. Yang, Jomaa, et al. (2018) when implementing nitrate transport in the model. Parameters that characterize the transformation of nitrogen between different nitrogen pools (Figure 1(c)) and the denitrification rate in the soil are land use-dependent parameters. They are modified in space and time according to the environmental conditions (soil moisture and soil temperature). For a more detailed description of the mHM-Nitrate model, the reader is referred to X. Yang, Jomaa, et al. (2018).

The Proposed mHM-SAS Model
The mHM-SAS model uses (1) the mHM concept for simulating hydrological processes, (2) the mHM-Nitrate concept for describing the nitrogen dynamics within the soil zone, and (3) the transit time formulation of transport based on SAS functions for representing nitrate transport and removal below the soil zone (Figures 1(b) and 1(c)). In contrast to mHM-Nitrate, mHM-SAS considers the unsaturated and saturated zones over the entire catchment as a single hydrological unit, in the following referred to as the SAS compartment. Hydrological fluxes into and out of the SAS compartment were simulated by the hydrologic routines of the original mHM model (Samaniego et al., 2010). Similar to the mHM and mHM-Nitrate model, we assume that evapotranspiration only occurs within the soil zone, in other words, there is no evaporation from the SAS compartment. Hydrological flux out of the SAS compartment is the summation of groundwater flows to the stream (baseflow) and other shallow subsurface flows (interflow). The spatially distributed hydrologic and nitrate fluxes from the soil zone to the SAS compartment are spatially lumped over the entire catchment. In the SAS compartment, we only track the fate of nitrate in the DIN pool (representing mainly N-NO 3 ).
The SAS compartment is a hydrological system with inflow   J t [L 3 T -1 ] and discharge   In the SAS concept, the system is conceptualized as storage of different water parcels with different ages (and solute concentration (Benettin & Bertuzzo, 2018;Botter et al., 2011;Harman, 2015;Van der Velde et al., 2012). Similarly, discharge is characterized by the TTD   . The volume of water with age younger than T, the so-called age-ranked storage   The transit time formulation of transport based on the SAS concept can be described as follows. The inflow to the SAS compartment and its associated transport component (nitrate) are considered to have an age of zero at the time of entry. Therefore, the transit time in this study refers to the elapsed time since water enters the SAS compartment to the time it leaves the SAS compartment via discharge. The residence time refers to the elapsed time since water enters the SAS compartment to the time of calculation that it is still inside the SAS compartment. Changes in the stored water volume with the age younger than T are induced by inflow , and aging. This can be described by the water age balance equation (e.g., Benettin & Bertuzzo, 2018;Botter et al., 2011;Harman, 2015;Van der Velde et al., 2012): where 0 T S [L 3 ] is the initial age-ranked storage.
The key element of the SAS concept is providing a functional relationship between the age distribution in storage and discharge. Several forms of SAS functions have been proposed and discussed, for example, SAS functions could be in the form of (1) absolute age T (Botter et al., 2011), (2) age-ranked storage   , T S T t (Harman, 2015), or (3) normalized age-ranked storage   , S P T t (Van der Velde et al., 2012). In this study, we used the SAS function as a function of normalized age-ranked storage as it is easy to be parameterized. In other words, the relation between the age distribution of storage and discharge is expressed as with S P varies from 0 to 1.  -3 ] is the solute concentration associated with input J at time t T  , subsurface denitrification is assumed to follow an exponential decay function with k [T -1 ] is the first-order denitrification rate constant, and   In this study, point sources such as discharge from wastewater treatment plants (WWTPs) and direct runoff from sealed areas are added directly to the catchment outlet. The flow-weighted mean concentration was used to calculate solute concentration at the catchment outlet.

Parameterization of the SAS Function
The mHM-SAS allows users to select either the power law (Benettin & Bertuzzo, 2018) or the beta distribution function (Benettin & Bertuzzo, 2018;Van der Velde et al., 2012, 2015J. Yang, Heidbüchel, et al., 2018) as an approximation of the SAS functions: where  and a, b are parameters of the power law (plaw) and beta (beta) distribution functions, respectively (with α, a, b ∈ (0, +∞)) and Γ is the gamma function. Different selection schemes of discharge from storage can be represented by varying parameters of the power law or beta distribution functions within defined ranges (Table 1). Van der Velde et al. (2015) suggested using the one-parameter beta function   , ,1 s beta P a to represent the young-water selection preference and   ,1, s beta P b to represent the old-water selection preference) instead of the two-parameter beta function. However, J. Yang, Heidbüchel, et al. (2018) demonstrated for a small agricultural catchment that both parameters of the beta function might vary in a wide range. Hence, we opted for the two-parameter representation of the beta function.
The SAS function (selection preference scheme for discharge) could vary in time, representing the dynamics of subsurface mixing under different conditions (J. Yang, Heidbüchel, et al., 2018). The time-variance of SAS functions could be related to subsurface storage volume, for example, using a linear functional relationship ( Van der Velde et al., 2015). This "one-to-one" relation, however, might not be sufficient to characterize the dynamics of the selection preference scheme for discharge due to hysteretic behavior of the system (e.g., different selection preference schemes corresponding to the same storage; Benettin, Bailey, et al., 2015;Benettin, Kirchne, et al., 2015). J. Yang, Heidbüchel, et al. (2018) found that the selection preference scheme for discharge depends not only on the current storage but also on the antecedent inputs (e.g., inflow to the SAS compartment during the previous time steps). Furthermore, the authors found that the selection preference scheme for discharge could be lumped together according to a seasonal hydrological situation such as wetting and drying phases during a year. This makes sense for catchments with significant seasonal variation in storage and meteorological forcing conditions as in J. Yang, Heidbüchel, et al. (2018).
In this study, we introduce a new approach for determining the transition between different selection preference schemes for discharge. In this approach, we assume that the young water fraction of streamflow increases with increasing catchment wetness as new fast shallow flow paths are activated, creating a different NGUYEN ET AL. Young-water selection preference Both young-water and old-water selection preference -0 < a, b < 1

Table 1 Selection Preference Schemes and the Corresponding Parameter Ranges of the Power Law and Beta Functions
selection scheme (e.g., Dupas et al., 2017;J. Yang, Heidbüchel, et al., 2018). The catchment wetness is reflected in both antecedent inflow and outflow. Therefore, we propose using the following ratio for determining changes in the selection preference scheme: if 1 t r  : preference for young water (Table 1) if 1 t r  : preference for old (and young) water (Table 1) where i t [T -1 ] is the current time and n [T] is the time window. The ratio t r [-] is a time-variant factor due to the temporal variations of inflow and outflow. The ratio t r explicitly considers the antecedent inflow and implicitly considers the changes in storage. For example, 1 t r  ( 1 t r  ) indicates that the storage is filling (emptying). In this study, the period with 1 t r  is referred to as the wet period while the period with 1 t r  is referred to as the dry period. An example for the relation between t r and storage is shown in Section 3.4. The advantages of relating the selection preference scheme to the ratio t r are (1) information about the minimum and maximum storage is not required, (2) the initial storage does not affect the selection preference scheme, and (3) no prior knowledge about the seasonal changes of storage is needed. It should be noted that the proposed two selection preference schemes, preference for young and preference for old (and young), were shown to be (1) sufficient for describing subsurface mixing ( Van der Velde et al., 2015) and (2) the dominant selection schemes in a subcatchment of the studied catchment (J. Yang, Heidbüchel, et al., 2018).

Study Area and Data
The study catchment is that of the upper Selke River (gauge Silberhütte), which is part of the Bode catchment, a terrestrial environmental observatory within the TERENO network of observatories in Germany (Wollschläger et al., 2017;X. Yang, Jomaa, et al., 2018). The study site covers an area of about 100 km 2 with elevation ranging from 335 m to 595 m above mean sea level (Figure 2(a)). Forest and agricultural land (pasture and arable land) are the dominant land uses/land covers in the area, accounting for 61% and 36% of the total area, respectively ( Figure 2(b)). The main crops planted in the area are winter wheat, triticale, winter barley, rye, rapeseed, and corn (Jiang et al., 2014;X. Yang, Jomaa, et al., 2018). Spodic Cambisols from hard argillaceous and silty slates account for about 70% of the study area while Dystric Cambisols from acid igneous and metamorphic rocks account for 26% of the study area ( Figure 2(c)). The geology of the study area is predominantly characterized by Mississippian wacke/shale, covering 99% of the area (X. Yang, Jomaa, et al., 2018). The aquifers of the study area are relatively shallow (Dupas et al., 2017;J. Yang, Heidbüchel, et al., 2018).
The study area has an average annual precipitation of 765 mm. The average monthly temperature in the area ranges from -3.1°C in December to 16.7°C in July. The area has a strong seasonal runoff regime with high flows during the cold season (November-April; average discharge average Q = 1.7 m 3 /s) and low flows during the warm season (May-October; average Q = 0.5 m 3 /s). About 77% of the total runoff is generated during the cold season. Diffuse nitrogen (N) from fertilizers applied to agricultural fields (with an average application rate of about 130-190 kg N ha -1 year -1 ) is the main source of instream N (Jiang et al., 2014;Kistner, 2007). Contribution from WWTPs to instream N is negligible during high flow periods. During low flow periods, however, N from the WWTPs can account for up to 20% of the total N in the stream.
Input data were obtained from different sources. Daily weather data (precipitation, temperature) and potential evapotranspiration were obtained from the Deutscher Wetterdienst. Digital elevation model, land use map of 30 m resolution, and soil map at a scale of 1:1,000,000 were provided by the Federal Institute for Geosciences and Natural Resources, Germany (BGR). Agricultural practices (fertilizer/manure application, crop rotation) were obtained by field survey/interview (X. Yang, Jomaa, et al., 2018). Daily discharge and weekly instream nitrate concentration were taken from the State Agency for Flood Protection and Water Management of Saxony-Anhalt.

Selection of the SAS Function and Initial Conditions
In this study, the beta functions were used to represent SAS functions because of their flexibility to represent more mixing schemes than the power law function (Table 1). In addition, beta functions have been found to be good approximations of model-derived SAS functions within a subcatchment of the study area (J. Yang, Heidbüchel, et al., 2018). Beta functions with time-variant parameters (a and b) were used to represent the temporal dynamics of the SAS functions. Two SAS functions were defined according to the wetness condition indicated by t r : wet SAS for the wet period ( 1 t r  ) and dry SAS for the dry period ( 1 t r  ). t r was calculated with n = 90 days, which was based on an iterative approach using a range of values to select the one that most suitably represents the seasonal patterns of the selection functions. However, it could be treated as a model parameter and determined via model calibration.
An initial nitrate concentration C 0 of 1.5 mg/L in subsurface water was selected based on the average instream nitrate concentration. The initial subsurface storage S 0 indicates not only the subsurface storage at the beginning of the simulation but also the subsurface storage capacity in general. A reliable estimation of the subsurface storage, which actively participates in the transport process requires extensive data (e.g., Hale et al., 2016). As these data were not available for the study area, we consider the initial storage as a model calibration parameter. The initial age-ranked storage ( 0 T S ) is assumed to linearly increase from 0 to S 0 over the age range [0, 10] years and all water in storage is assumed to have the same initial concentration C 0 .

Model Calibration and Uncertainty Analysis
The Elementary Effect Test (EET, Campolongo et al., 2007;Morris, 1991;Pianosi et al., 2016) has been proven as an effective tool for parameter sensitivity analysis for the study area (X. Yang   x is the average of the absolute elementary effect in r trajectories (Campolongo et al., 2007), which is calculated as follows: where j x is the vector of parameter values in the jth trajectory,  (Campolongo et al., 2007;Morris, 1991): In the EET, higher values of i   indicate higher sensitivities of the respective parameter while higher values of i  indicate stronger interactions of that respective parameter with other parameters. In this study, the Sensitivity Analysis For Everybody (Pianosi et al., 2015) toolbox was used to perform the EET for 54 global parameters, including the initial subsurface storage (S 0 ). Parameter sensitivity analyses were carried out separately for discharge and instream nitrate concentration at the catchment outlet.
For parameter optimization, we performed 20,000 simulations with parameters generated from Latin Hypercube Sampling (LHS). LHS has been demonstrated as an efficient global sampling procedure for optimization problems with a large number of parameters (Abbaspour et al., 2004). The best simulation was selected based on the following multicriteria objective function (OF ): where NSE and ln NSE are the Nash-Sutcliffe Efficiency (Nash & Sutcliffe, 1970) and its logarithmic transformation, obs x and sim x are the observed and simulated values of discharge Q or instream nitrate concentration C, and obs x and ln obs x are the mean and the logarithmic transformation of the observed variables, respectively. The NSE and ln NSE were used to ensure accurate modeling of both high and low values of discharge and nitrate concentration. In addition to the NSE and ln NSE, the percentage bias (PBIAS, Equation 18) was also used to evaluate the best simulation: The model prediction uncertainty is defined as a function of parameter uncertainty. Parameter uncertainty is characterized by the 95% prediction uncertainty (95PPU) estimated from behavioral simulations obtained from 20,000 Latin Hypercube simulations. We classified simulations with an objective function value great-er than 0.65 as behavioral. These behavioral simulations are well above the satisfactory level as suggested by others (e.g., Moriasi et al., 2007). The 95PPU was calculated based on the 2.5% and 97.5% levels of the cumulative distribution of the output variable at every simulated time step (e.g., Abbaspour et al., 2004;Beven & Binley, 1992). The goodness of the 95PPU is evaluated by the p-factor (the percentage of observed data bracketed by the 95PPU) and r-factor (the average thickness of the 95PPU band divided by the standard deviation of the observed data) (Abbaspour et al., 2004). The closer the p-factor to 1 and r-factor to 0, the better the 95PPU. In this study, the model was run at a daily time step with a 2-year warm-up (1993)(1994), 10-year calibration (1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004), and 10-year validation period (2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014). The spatial resolution of each grid cell is 1 km 2 .

Parameter Sensitivity Analysis and Calibrated Parameter Values
Parameter sensitivity analyses were carried out separately for discharge and instream nitrate concentration at the catchment outlet. Results show that discharge generation is most sensitive to soil parameters (soil10, soil11, soil13), followed by interflow parameters (intfl5, intfl4, intfl2), percolation (percol), evapotranspiration (pet1), and the snow parameter (tsnow) (Figure 3(a)). Instream nitrate concentration is sensitive to both hydrological and nitrate parameters (Figure 3(b)). With nitrate parameters, the denitrification rate constants in the soil zone (denis na , denis a ) and below the soil zone (k) are the most sensitive parameters. It is seen that the initial subsurface storage (S 0 ) is also listed among the most sensitive parameters, indicating the potential impact of subsurface storage capacity on catchment-scale nitrate export. Most of the parameters of the SAS functions (a wet , b wet , a dry ) are identified as sensitive parameters. This shows that the selection preference scheme for discharge is highly relevant to the solute export dynamics. Regarding the interaction NGUYEN ET AL.
10.1029/2020WR028490 10 of 25 between parameters, the results show that more sensitive (higher   ) parameters tend to have higher interaction (higher  ) with other parameters. In this study, the 15 most sensitive parameters for discharge and instream nitrate concentration are selected for optimization ( Figure 3 and Table 2). In addition, the parameter b dry is also selected due to its high sensitivity ranking among nitrate parameters. The sign (") indicates that the information in this cell is identical to the cell immediately above it. PTF is the pedotransfer function.

Discharge and Instream Nitrate Concentration
Visual assessment and model performance indices (NSE, ln NSE, PBIAS) show that the proposed model was satisfactorily calibrated and validated for discharge at the catchment outlet (Figure 4(a)). Nonetheless, some high discharge events were underestimated and overestimated (Figure 4(a)). The underestimation and overestimation of individual high discharge events could be attributed to the uncertainty in the rainfall data, for example, underestimation and overestimation of rainfall in some regions. From the flow duration curve (Figure 4(b)), it is seen that low flows were well represented by the model. In addition, the 95PPU band covers 96% of the observed values (p-factor = 0.96) with a narrow band (r-factor = 0.59).  off events with high nitrate concentration were underestimated while runoff events with low concentration were overestimated (Figures 4(c) and 4(d)). The C-Q relationship shows that high concentrations are associated with high flows and low concentrations are associated with low flows. Therefore, the underestimation of high concentrations during high flow conditions could be attributed to (1) the unaccounted direct transport of nitrate from the agricultural field to stream via direct surface runoff and/or (2) the activation of preferential subsurface flow paths that are only activated during extreme events. The overestimation of low concentrations, however, only occurs during some years, especially during the validation period (Figure 4(d)). This could be due to the overestimation of N from WWTPs in some years that was set constant in time in this model (X. Yang, Jomaa, et al., 2018). This is the reason for a lower ln NSE found during the validation period compared with the calibration period. Overall, the model performance statistics for the nitrate concentrations are within acceptable ranges and the 95PPU covers 60% of the observed values (p-factor = 0.60). The 95PPU band for instream nitrate concentration (r-factor = 0.99) is higher than that for discharge (r-factor = 0.59) because the water quality simulation is subjected to additional uncertainty from the hydrological simulation.

Spatial Nitrogen Dynamics and Nitrogen Balance
The mHM-SAS model could provide detailed insights into the spatial nitrate dynamics within the soil zone ( Figure 5). In general, the catchment experiences a strong spatial variability in diffuse nitrate external input (mainly from fertilizer and manure application, Figure 5(b)), mineralization ( Figure 5(c)), plant-nutrient uptake ( Figure 5(d)), denitrification in the soil zone ( Figure 5(e)), and nitrate leaching ( Figure 5(f)). This is expected as the mHM-SAS and mHM-Nitrate (X. Yang et al., 2019;X. Yang, Jomaa, et al., 2018) models use the same concept for describing soil nitrogen processes. As indicated from the input data, diffuse nitrate inputs from agricultural lands are the most dominant diffuse N sources, a consistent pattern all over Europe (Leip et al., 2011). The simulated spatial patterns of N-fluxes due to mineralization, plant uptake, denitrification, and nitrate leaching mainly follow the spatial patterns of diffuse nitrate inputs (with a correlation coefficient > 0.95), higher rates in agricultural areas and lower rates in forest areas. Denitrification rate in agriculturally dominated soil (agricultural fraction > 0.5) is generally higher than in forest dominated soil (forest fraction > 0.5), on average 2.7 times higher. In agriculturally dominated areas, it is seen that a NGUYEN ET AL.

10.1029/2020WR028490
13 of 25 significant amount of nitrate is leached out of the soil zone despite a high rate of nitrate removal by plant uptake and denitrification. This is the major reason for the higher nitrate concentration that is observed in the groundwater zone below agricultural areas compared to forest areas (e.g., Knoll et al., 2020).
For a long-term nitrate balance within the soil zone, the model suggests that most of the nitrate input (59.1 kg N ha -1 year -1 ) and mineralization (20.6 kg N ha -1 year -1 ) to the DIN soil pool was removed via plant uptake (45.9 kg N ha -1 year -1 ), followed by soil denitrification (23.2 kg N ha -1 year -1 ), and finally leaching to the deeper subsurface (11.2 kg N ha -1 year -1 ). In agriculturally dominated areas, denitrification in the soil zone is the largest nitrogen loss pathway, which is common for European agricultural soils (Velthof et al., 2009) but also observed elsewhere . Modeling results indicate that there is almost no long-term accumulation of nitrate in the soil zone. The simulated rates of mineralization (20.5-20.7 kg N ha -1 year -1 ) and denitrification (18.8-31.1 kg N ha -1 year -1 ) in this study are within the measured range reported by Heumann et al. (2011) (mineralization rate: 14-187 kg N ha -1 year -1 ) from different soil types in Lower Saxony, Germany and by Hofstra and Bouwman (2005) (denitrification rate: 8-51 kg N ha -1 year -1 ) from 336 agricultural soils located worldwide. The simulated yearly average N surplus (nitrate input + mineralization -plant uptake) from the optimal parameter set is 33.8 kg N ha -1 year -1 . This is comparable with the calculated N surplus (33 kg N hav -1 year -1 ) from other studies for the area (Häußermann et al., 2019;Winter et al., 2021). In terms of storage, the mobile N storage (DIN) in the soil and groundwater is around 22.3 kg N ha -1 and 4.6 kg N ha -1 , respectively. The total organic N (DON, SON A , and SON I ) in the soil is around 188.8 kg N ha -1 , accounting for 89% of total nitrogen in the soil. This fraction of inorganic nitrogen in the soil is comparable with the value reported by Stevenson (1994) (over 90% of N in most of the soil is in the form of organic N). It should be noted that the total inorganic N storage in the soil zone is highly dependent on the initial storage of the SON I pool. In contrast, instream N is not sensitive to the initial storage of the SON I pool because the degradation and/or mineralization rates are very slow.
A substantial part (about 32%) of the N surplus was leached out of the soil zone to the SAS compartment ( Figure 5(f)). Within the SAS compartment, nitrate is further denitrified along its transport pathways and is removed via discharge. The long-term nitrate balance (from 1995 to 2014) from the optimal parameter set shows that about 37% of the leached nitrate ( Figure 5(f)) was removed via denitrification and 62% (54% during the wet periods and 8% during the dry periods) was exported to the stream. In the study area, different magnitudes of denitrification potential in groundwater have been reported based on measured groundwater chemistry data, ranging from high to low nitrate reduction potential (Hannappel et al., 2018). We thus conclude that the modeled denitrification rate below the soil zone is acceptable.

Subsurface Mixing and Transport
The results show that the selection preference for discharge has a consistent seasonal pattern (Figures 6(a) and 6(b)). In general, it is seen that the system preferentially selects (1) young water in storage for discharge during wet periods with high subsurface storage and (2) both young and old water in storage for discharge during dry periods with low subsurface storage. The dominance of young water in discharge during wet periods is mainly attributed to the activation of fast shallow flow paths under high wetness conditions (Dupas et al., 2017;J. Yang, Heidbüchel, et al., 2018). The infiltrated rainfall takes a relatively short time to travel via these flow paths, providing streamflow with dominant young water. This results in a much smaller median transit time compared to the median residence time during the wet periods (Figure 6(c)). The preference for young water during selected dry periods is due to the fact that occasional rain events with high intensity lead the activation of fast shallow flow paths. When there is no rainfall or rainfall with low intensity, stream discharge is mainly composed of older water due to the deactivation of the fast shallow flow paths and a dominance of slow deep flow paths. As a result, the median transit time (TT 50 ) for discharge in the dry periods is considerably longer than that in the wet periods (Figure 6(c)).
The temporal activation and deactivation of different flow paths affect the age composition of discharge, the young and old water fraction in discharge, and ultimately the dynamics of nitrate in discharge. This is because shorter transit times indicate less time for denitrification and a dominance of young water in discharge indicates a pronounced effect of recent agricultural activities on instream water quality. It is seen that instream nitrate concentration in the wet periods is higher than that in the dry periods (Figure 6(d)).
Higher nitrate concentrations in the wet periods are due to higher fractions of young water (with higher nitrate concentrations). Lower nitrate concentrations in the dry periods are due to a mixture of old water (with low nitrate concentration due to a long time of denitrification) and young water (with low nitrate concentration) (Figures 6(a) and 6(d)). Subsurface mixing and denitrification also result in a lower temporal variability of instream nitrate concentration compared to that of leached nitrate (Figure 6(d)).

Transport and Reaction Time Scales
To explore the interplay between transport and reaction rate on nitrate export, we use the Damköhler number (D a , Ocampo et al., 2006).  of reaction over transport while D a values <1 indicate the dominance of transport over reaction. The simulated average D a number for the wet periods is 1.62 compared to the average D a number of 8.03 for the dry periods. This shows that nitrate transport during the dry periods is characterized by a much more pronounced dominance of reaction over transport.
In the study area, the simulated mean transit time (MTT) over the simulation period (1995-2014) is 2.34 years. This is comparable with the MTT estimated based on stable isotope data for the Meisdorf gauging station (2.19 years) located further downstream of the study area (Lutz et al., 2017). In their study, the MTT was calibrated using stable isotope data and young water fraction under the assumption of a gamma-shaped TTD. In a recent study for the study area, Winter et al. (2021) assumed that TTDs follow a lognormal distribution. Parameters of the lognormal were determined from a comparison of the long-term changes in annual N input and flow normalized nitrate concentrations observed in the output. Their results indicate that the MTT is about 2.12 years, which is comparable to our finding. However, the calculation of the MTT in this study is subjected to a certain degree of uncertainty as described below. Nevertheless, a similar MTT obtained from our study compared to other data-driven approaches validate our findings and thus illustrate the potential of a robust application of the proposed model.
In this study, we found that the variables MTT and D a which take into account the oldest water are highly sensitive to the actual age of the oldest water. Information on the age of the oldest water cannot be determined from the observed instream nitrate time series or the model which is calibrated against that (e.g., Stewart et al., 2010). This is because nitrate in old waters was effectively denitrified to the level that is below the lower detection limit. The aforementioned MTT (2.34 years) was calculated with an assumption that the maximum age in storage is 10 years, older water is merged to the "old" water pool (with the age of 10 years and an average volume of 46% the total subsurface storage). In other words, it means that old water (water with age ≥10 years or D a ≥ 22.2) is assumed to be well mixed. Under the assumption that the oldest water in storage is not restricted to a certain age (the oldest water becomes older as the simulation time increases), the MTT of discharge shifts to 4.03 years. In terms of instream nitrate concentration and TT 50 , the two aforementioned assumptions give almost identical results (Section 3.6). Similar results (instream nitrate concentration and TT 50 ) are obtained if the maximum age in storage is limited to 1 year. For solute export, the results indicate that when reaction strongly dominates transport, mixing within the old water storage with very low nitrate concentration compared to young water does not affect solute concentration in the outflow. Figure 7 shows the relation between the TTD of discharge and nitrate load in discharge on the typical wet and dry days. It is seen that the age of discharge on the dry day is much older than on the wet day. In both dry and wet days, the majority of nitrate in discharge is from the young water fraction of discharge. The cumulative nitrate load does not follow the cumulative TTD because of denitrification. On the dry (wet) day, about 75% (85%) of nitrate in discharge is younger than the median transit time. On the dry day, a very small fraction (<1%) of nitrate in discharge is transported by water that is older than a year despite a high overall fraction of old water (about 40% of discharge is older than a year). This result further confirms that in the study area, a detailed representation of mixing inside the old water pool (>1 year) and the age distribution in this pool are not necessary for representing instream nitrate dynamics when the reaction rate is high.

Time Lags in Catchment Response
To understand the time lag between nitrogen input and catchment solute export, a hypothetical scenario was set up. In this scenario, all nitrogen inputs to the soil (fertilizer, manure, atmospheric deposition, and NGUYEN ET AL.   (December 15, 2002) and dry day (August 9, 2003). The x axis is represented on a log-scale for better visualization. The data were derived from the optimal parameter set. It should be noted that the oldest water in discharge during the dry day could be much older than 10 years, however, it was merged to the water pool with the age of 10 years. residues) are stopped after a certain time (Figure 8). The time lag between input nitrogen and instream nitrate concentration signals can be due to biogeochemical (soil) and hydrological (groundwater) time lags (Ilampooranan et al., 2019). In this study, the biogeochemical time lag corresponds to the biogeochemical reaction time scale in the soil zone while the hydrological time lag corresponds to the travel time of nitrate in the subsurface. The time lag between input nitrogen and leached nitrate concentration signals reflects the biogeochemical time lag while time lag between leached nitrate and in instream nitrate concentration signals reflects the hydrological time lag. Figure 8 shows that there is an increase of instream nitrate concentration following a decrease and a complete cessation of all input N. However, the delay between leached nitrate and instream nitrate concentration signals is not clear. This indicates that the biogeochemical time lag in the study is more pronounced than the hydrological time lag. This is because of a short transit time (J. Yang, Heidbüchel, et al., 2018), a dominance of young water fractions in discharge, and a high reaction rate as mentioned in the previous sections. In general, the time lag between N input and output in this catchment is very fast (within a year). A recent data-driven analysis by Winter et al. (2021) for the same area also shows an immediate decrease of instream nitrate concentration and loading following a drastic decrease of N input after 1990. In the lower Selke catchment, longer time lags (up to 12 years) were estimated due to drier conditions and deeper aquifers (Winter et al., 2021). This indicates a high variability of N-related time lags in the region. Dupas et al. (2020) also showed that N lag times in various catchments in western France vary over a wide range (from 2 to 14 years), depending on hydrogeological settings of the area. Evidence from these other studies suggest that N lag times in this study are within a reasonable range.
The biogeochemical time lag could also be explained as follows. About one-half years after the complete cessation of all input N in the modeling scenario, most of the active organic N (SON A and DON) was converted to in DIN, which was subsequently removed by denitrification, plant uptake, and leaching. The only N source that remains was SON I , which was slowly converted to DIN due to low degradation and/or mineralization rates. SON I has been recognized as a long-term source for mineralization and nitrate leaching in other areas (Van Meter et al., 2016). Due to a high denitrification rate in groundwater, the DIN leached from the soil to groundwater was mostly denitrified, resulting in very low instream nitrate concentrations (Figure 8).

Model Comparison
The proposed mHM-SAS model uses the time-variant SAS functions to describe subsurface mixing dynamics and time-variant TTDs of discharge. Within this approach, both celerity-driven and velocity-driven transport mechanisms are taken into account. Transport of reactive solutes (e.g., nitrate) and different types of subsurface mixing behaviors could be incorporated into the model. Compared with the approach that uses the hydrologically inactive and active groundwater storage with either a complete or a partial mixing assumption (e.g., Hrachowitz et al., 2013Hrachowitz et al., , 2015X. Yang, Jomaa, et al., 2018), our approach provides a more general description for subsurface mixing and does not require the separation/estimation of the active and inactive storages. In addition, applications of these studies (Hrachowitz et al., 2013X. Yang, Jomaa, et al., 2018) are limited to either the transport of conservative solute or the transport of reactive solute without considering biogeochemical processes that could occur along subsurface flow paths. Similarly, the approach with single subsurface storage and a complete mixing assumption has been successfully applied for modeling transport of conservative solute (Remondi et al., 2018), however, the application of this approach with reactive solute is unknown. In this study, we have demonstrated that the SAS concept could be applied for describing transport of reactive solute while considering biogeochemical processes that occur along subsurface flow paths.
Compared with the mHM-Nitrate, the mHM-SAS model has more parameters to represent additional processes (denitrification and preferential mixing behavior in the subsurface) instead of simply increasing the model degrees of freedom. Results from the two models for the upper Selke catchment shows that they have a comparable model performance for instream nitrate concentration at the catchment outlet. From an optimization viewpoint, if two models represent the same number of processes, the model with a higher number of parameters is expected to have higher (or at least equal) model performance (e.g., for instream nitrate concentration). However, this is not the case with the mHM-SAS and mHM-Nitrate model. Therefore, the two models (mHM-SAS and mHM-Nitrate) should not be compared using only the simulated instream nitrate concentration. The mHM-SAS model could be considered as a generalization of the mHM-Nitrate model with a more realistic representation of subsurface mixing and nitrate transport. By explicitly representing the travel times and denitrification in the subsurface, the mHM-SAS could represent both hydrologic and biogeochemical legacy behaviors (e.g., Kumar et al., 2020;Van Meter et al., 2016, 2017. The mHM-SAS model is therefore of particular interest for understanding time lags between agricultural nutrient inputs and catchment responses.

Model Capabilities and Limitations
The simulated spatial nitrate patterns have highlighted the necessity of a spatially explicit representation of nitrate dynamics within the soil zone. This could help to identify critical source areas and to advise better management practices. In the catchment-scale application of the SAS approach, the spatial patterns in nitrate leaching from the soil zone are not explicitly considered in the transport process. This SAS approach transfers the transport problem into the time domain and only considers the dynamic distribution of transit times due to the heterogeneity of subsurface transport pathways. In other words, this approach provides insights into the time origin of discharge and the solutes in discharge instead of their spatial origin.
The mHM-SAS model could provide insights into the functioning of the catchment (subsurface mixing) and the internal dynamics of discharge (TTD) and solute in discharge unlike traditional conceptual water quality models (Hrachowitz et al., 2016). The tested catchment is characterized by a small and reactive catchment storage that leads to a fast reaction time of instream nitrate concentration to changes in the input. In catchments with larger groundwater storages and transit times, the long-term effects of biogeochemical and hydrological legacies can play out very differently (Ehrhardt et al., 2019, Van Meter et al., 2018. Here, our modeling approach could serve as an investigation tool for quantifying the long-term memory effects of historical agricultural practices on the present surface water quality status. Understanding the temporal dynamics of subsurface mixing and TTD also allows us to identify when instream water quality is more vulnerable to input contaminants and to develop better management practices. Despite the aforementioned model capabilities, the model is still a simplified representation of the real system and further developments are suggested. The current version of the model does not consider travel times in the soil (root) zone, which could be an important source for hydrologic legacy (Kumar et al., 2020). In this study, temporal changes of the selection preference scheme are not continuous, in other words, the relation between the factor t r and parameters of the beta function is quite abrupt. The selection preference scheme could gradually shift from young to old water affinity when the catchment changes from wet to dry (J. Yang, Heidbüchel, et al., 2018). In addition, instream nitrate removal is not explicitly considered. Instead, it is lumped with subsurface nitrate denitrification. However, the travel time in the stream network is of different magnitudes compared to the travel time in the subsurface, therefore, separation of these processes is required for the areas where instream nitrate removal is significant. In our study area, instream nitrate removal is negligible (X. Yang, Jomaa, et al., 2018). The current lumped version of the mHM-SAS model does not consider the spatial variability of leached nitrate out of the soil zone. To preserve the spatial information of leached nitrated from the root zone in the transport process or to answer the question about the spatial origin of discharge at the catchment outlet, a spatially more resolved, grid-based application of the SAS concept is required. This also applies when the model is transferred to larger basins with a distinct spatial heterogeneity in subsurface properties that does not allow for an effective lumped parameterization.

Toward a Spatially Distributed SAS-Based Model
In this section, we evaluate the potential applicability of a spatially distributed SAS-based model. In this approach, SAS functions are applied at the grid cell level. All grid cells are assumed to have the same SAS parameters (including the initial conditions), which are the optimal values obtained from the lumped approach (Table 2). This assumption reflects a case with homogeneous hydrogeological settings where outflow from each grid cell is directly discharged to the stream (no subsurface flow between grid cells). Changes in the mixing scheme in each model grid cell are defined by the antecedent inflows and outflows as described in Section 2.3. Results show that the simulated instream nitrate concentration and the median transit time of discharge at the catchment outlet from the two approaches are almost identical (Figure 9). It should be noted that for other catchments, the results from the two approaches could be different even though the SAS functions at the grid cell level (distributed approach) and the basin level (lumped approach) are the same (e.g., Kirchner et al., 2001). In this study, the results could indicate that (1) the spatial information about nitrate fluxes from the root zone and (2) subsurface connectivities between grid cells in the study area have only minor effects on catchment nitrate export and the catchment-scale median transit time of discharge. Satisfactory results from the distributed approach also show that the assumption about homogeneous SAS parameters could be a valid assumption for the study area.
In the spatially distributed approach, the model can provide spatial information about, for example, the age of storage (residence time, RT) and discharge (transit time, TT) ( Figure 10). This information has significant implications for the understanding of flow and transport of contaminants. It is seen that even though the spatial patterns of residence times, which are characterized by the median of the median RT, are far from homogeneous (Figure 10(a)). In this example, the spatial patterns of the residence time are mainly controlled by the spatial pattern of recharge, the median RT 50 is inversely correlated with the recharge rate with a correlation of -0.94. The recharge rate is further controlled by precipitation, land cover, topography, and soil  properties. In this study, it is seen that shorter residence times are observed in highland areas while longer residence times are observed in lowland areas. Shorter (or longer) residence times indicate faster (or slower) responses of stream water quality to changes in agricultural practices. However, this information could be less relevant for stream water quality in the study area because of the dominance of reaction over transport. The spatially distributed approach also allows us to explore the spatial patterns of the transit time of discharge ( Figure 10(b)). It is seen that even though the mixing scheme is spatially homogeneous, the transit time of discharge is highly heterogeneous. In general, the spatial pattern of the transit time of discharge ( Figure 10(b)) follows the spatial pattern of the residence time (Figure 10(a)). Shorter transit times indicate higher vulnerabilities of stream water quality to input contaminants. The evolution of the transit times along the river network is shown in Figure 10(c). Changes in the transit time of discharge along the river network are expected because the main river receives discharges from tributaries with different TTDs along its course. The temporal variation of the RT and TT (lower panel, Figure 10) indicates that the TT of discharge has a higher temporal variation than the RT. This is due to the seasonal changes in the mixing scheme.
A major disadvantage of most distributed conceptual hydrological models is that lateral subsurface flow and transport between model grid cells is usually neglected, for example, Variable Infiltration Capacity model (Liang et al., 1994), mHM-Nitrate (X. Yang, Jomaa, et al., 2018), GROWA-DENUZ-WEKU (Kunkel et al., 2017). Thus, the maximum flow path length is limited to the cell size. This could be true if the grid resolution is low (e.g., cell sizes are as large as the subsurface catchment size), subsurface flow occurs only within a given grid cell. However, if the grid resolution is high (small cell sizes), water and solutes from the upstream grid cell can be transported to downstream grid cells and mixed with water and solute in these grid cells before entering the river. The response of instream solute concentration to the input signal from a cell located at a distance could be slower than the response to the input signal from a cell located nearby. In other words, there would be legacy effects due to the longer transit times of nutrients from regions, which are not directly connected to the stream network (Figure 11). In a fully spatially distributed approach, which implicitly accounts for lateral subsurface flows between grid cells, the aforementioned flow and NGUYEN ET AL. transport mechanisms could in principle be represented. For example, transport of water and solute from a grid cell located far away from the river could be conceptualized with a selection preference for older water compared to the selection preference for younger water for the cell located near the river (Figure 11). Mixing occurring along longer flow paths could be conceptualized as mixing in the river, where the flow contributions from all distant and close grid cells are eventually combined (e.g., Kirchner et al., 2001). This example ( Figure 11) indicates the potential of application of the fully distributed SAS-based model for representing lag times of N inputs and outputs due to hydrologic legacy. For a reactive solute such as nitrate, the distributed approach would also allow to vary denitrification rates between grid cells.
Despite these potential advantages of a fully distributed approach, several challenges would have to be overcome in its implementation. For example, the functional relationships between grid cell characteristics (e.g., meteorological forcing, hydrogeological properties, and location of the grid cell) and parameters of the SAS functions need to be addressed. In addition, the fully distributed model will significantly increase the number of model parameters (e.g., parameters of the SAS function could be changed in space and time), which could lead to the problem of overparameterization. The distributed approach will also require more computational and storage capacity compared to the lumped approach. Furthermore, additional field data would be required to constrain or verify the spatially resolved output from the model to ensure model robustness. A potential alternative approach could be using a regional-specific SAS function based on the region characteristics (e.g., hydrogeological conditions) or alternatively a few different SAS functions for a set of geologically meaningful subsurface units instead of a fully distributed approach. In addition, the advancement of physically based groundwater models as tools to evaluate processes more mechanistically as well as an increasing amount of field data from experimental catchments could help to alleviate some of these verification problems.

Conclusions and Outlook
SAS functions have emerged as a novel tool for modeling solute transport at the catchment scale. However, a thorough representation of the spatial heterogeneity of catchment characteristics (e.g., land use, soil, NGUYEN ET AL. 10.1029/2020WR028490 21 of 25 Figure 11. A hypothetical example for the representation of transport of conservative solute from different grid cells to a river with the SAS-based approach. In this example, both grid cells are assumed to have the same initial storage (S 0 = 500 mm), initial concentration (C 0 = 0 mg/L), impulse input signal (C = 1 mg/L) at time t = 0, and constant input and output fluxes (Q in = Q out = constant = 1 mm/day). SAS, StorAge Selection. and topography) in such models and a systematic testing of SAS-function-based models at larger scales (e.g., mesoscale-catchments) have not been done to date. In this study, we took a step in this direction and integrated a SAS-function-based nitrate transport model into a fully distributed soil nitrate model (mHM-Nitrate) at both the catchments as well as grid cell scales, resulting in the mHM-SAS model. Seasonal variations in the age selection schemes of the catchments as represented by shifting SAS functions were implemented in the model based on antecedent inflows and outflows to the subsurface compartment of the model (i.e., entire catchment or on the grid cell level). For the first time, to the best of our knowledge, the SAS concept has been evaluated in a mesoscale catchment (100 km 2 ) with heterogeneous catchment characteristics (land cover, land use management practices, and soil types). Key results show that • Denitrification below the soil zone could be significant and should be accounted for (e.g., the upper Selke in this study). • The lumped SAS-based approach could well represent streamflow and solute export dynamics of a mesoscale heterogeneous catchment with realistic reaction rates and transit times. • Both lumped and distributed SAS-based approaches yield comparable results in terms of instream nitrate dynamics and median transit times of discharge at the catchment outlet. • Temporal activation and deactivation of different flow paths control the transit time of discharge and solute export dynamics of the catchment. • Knowledge about the age of the oldest water in storage or discharge is not required for characterizing solute export dynamics from a highly reactive system. • Temporal changes in the SAS functions could be related to the antecedent inflow and outflow ratio, which does not explicitly require prior knowledge about subsurface storage (e.g., minimum, maximum, and seasonal changes). • Heterogeneity in the recharge rates controls the spatial patterns of transit times.
This study has demonstrated the general applicability of SAS-function-based solute transport models to mesoscale catchments. However, the application of the SAS concept at this scale is still in an early stage. Testing of the SAS concept in other catchments with different settings is needed. The mHM-SAS model can be considered as the first prototype for a parsimonious SAS-function-based solute transport model for larger catchments. However, the proposed general integration framework could easily be applied to other distributed water quality models. X. Yang, and M. Rode for providing input data for setting up the mHM-SAS model. R. Kumar acknowledges the Initiative and Networking Fund of the Helmholtz Association through the project Advanced Earth System Modelling Capacity (www.esm-project.net). Source codes of the mHM-SAS model and relevant data for reproducing the work are available online at https://git.ufz. de/nguyenta/mhm-sas and https://git. ufz.de/yangx/mHM-Nitrate. We thank two anonymous reviewers for their constructive comments that helped to considerably improve the quality of the manuscript.