Disentangling the Influence of Landscape Characteristics, Hydroclimatic Variability and Land Management on Surface Water NO3‐N Dynamics: Spatially Distributed Modeling Over 30 yr in a Lowland Mixed Land Use Catchment

Nitrate (NO3‐N) mobilization is generally controlled by available sources, hydrological connectivity, and biogeochemical transformations along the dominant flow paths. However, their spatial heterogeneity and complex interactions often impede integrated understanding of NO3‐N dynamics at the catchment scale. To fully integrate spatiotemporal information for NO3‐N simulations, a grid‐based model, mHM‐Nitrate, was applied to a 68 km2 lowland, mixed land use catchment (Demnitzer Millcreek, DMC) near Berlin. The model successfully captured the spatiotemporal distribution of flow and NO3‐N between 2001 and 2019, but was less successful in 1992–2000 due to land management changes. Re‐optimization of relative parameters was subsequently conducted for this period to understand management effects. The simulated results revealed landscape characteristics and hydroclimatic variability as the main controlling factors on respective spatial and temporal patterns. The combined effects of vegetation cover and fertilizer inputs dictated the spatial distribution of water and NO3‐N fluxes, while wetness condition determined the temporal NO3‐N dynamics by regulating hydrological connectivity and NO3‐N mobilization. Denitrification was also closely coupled with hydroclimatic conditions, which accounted for ∼20% of NO3‐N inputs. In contrast, restoration of riparian wetlands had a modest impact on NO3‐N export (∼10% reduction during 2001–2019), suggesting further interventions (e.g., reducing fertilizer application or increased wetland areas) are needed. Our modeling application demonstrated that mHM‐Nitrate could provide robust spatially distributed simulations of hydrological and NO3‐N fluxes over a long‐term period and could successfully differentiate the key controlling factors. This underlines the model's value in contributing to an evidence base to guide future management practices under climate change.

Spatial distributed water quality modeling is one way to formalize such knowledge, as impacts of various controlling factors can be adequately incorporated via regional parameterization (Orban et al., 2010;Rozemeijer et al., 2016). A number of models have been developed to simulate both catchment water balance and NO 3 -N dynamics; ranging from statistical regression models to detailed physically based models (reviewed in Billen et al., 2011). According to a review of 257 distributed modeling applications (Wellen et al., 2015), ∼70% of these studies were based on semi-distributed models, particularly the well-known SWAT (Neitsch et al., 2011), INCA (Wade et al., 2002), HSPF (Bicknell et al., 1997), and HYPE models (Lindström et al., 2010), which have proven robust tools in terms of water and nutrient fluxes estimation, spanning across a range of catchment-scales (Bekiaris et al., 2005;Pisinaras et al., 2010) to the continental-scale (Abbaspour et al., 2015). However, the strategy of such semi-distributed models, disaggregating a catchment into several homogenous units (sub-catchments or Hydrological Response Units) based on topography (or other landscape characteristics), may lead to information loss regarding class locations and interactions between neighboring classes (Rathjens & Oppelt, 2015;Yang et al., 2018). Alternatively, fully distributed models, which divide the catchment into hydraulically connected elements such as grids or triangular elements, have also drawn interest in the modeling community (Wellen et al., 2015), given their capability to utilize raster information such as DTMs, as well as soil and land-use maps (Ouyang et al., 2010). Increasingly, grid-based water quality models have been developed and applied, generally with detailed physics-based description of hydrological fluxes and nitrogen transformations (e.g., Birkinshaw & Ewen, 2000a;Birkinshaw & Ewen, 2000b;Wriedt & Rode, 2006). For example, Birkinshaw and Ewen (2000a) successfully reproduced the NO 3 -N transport by coupling the NO 3 -N transformation model NITS with the physically based distributed hydrological model SHETRAN. However, application of such physically based models have generally been limited to data-rich catchments at smaller spatial scales (<20 km 2 ), because the subsurface representation based on finite difference cells and complex physical equations results in a computationally demanding calibration, and the need for extensive input data (e.g., prior knowledge of subsurface characteristics) that are rarely available at larger catchment scales (Rathjens et al., 2015;Schoumans et al., 2009). Ultimately, conceptual models, which balance the spatial representativeness and model complexity by simplifying descriptions of dominant physical processes, have many advantages for mesoscale catchment studies (e.g., AGNPS/ AnnAGNPS [Young et al., 1989], TNT2 [Beaujouan et al., 2001], DNMT [Liu et al., 2005]). It has been shown that such models were able to capture similar output simulations to fully physically based models, with much lower computational demands and data needs (Jackson et al., 2008;Quinn, 2004).
Recently, the fully distributed process-based mHM-Nitrate model was developed for mesoscale catchments (Yang et al., 2018). The new approach integrates hydrological and water quality concepts taken from two widely used models: the mesoscale Hydrologic Model (mHM; Samaniego et al., 2010) and the Hydrological Predictions for the Environment model (HYPE; Lindström et al., 2010). A grid-based structure with reservoir-based descriptions of water storage-flux relationships provides comprehensive insights into spatial patterns of the dominant processes governing both water and nitrogen fluxes, while retaining an intermediate model complexity. Moreover, the model uses a multiscale parameter regionalization technique that facilitates full assimilation of spatial information by parameterizing at finer scales and upscaling to coarser simulation scales (Samaniego et al., 2010).
In this study, we applied the mHM-Nitrate in a data-rich, mixed land-use, lowland catchment (68 km 2 ) in Germany over a 30 yr period to investigate the long-term dynamics of hydrological and NO 3 -N fluxes, and key controlling factors. Three main research questions were addressed: 1. Can mHM-Nitrate reproduce the spatial distribution of hydrological and NO 3 -N fluxes? 2. Did the spatially distributed results capture the temporal dynamics of NO 3 -N fluxes over a 30 yr period, including extreme wet and dry periods? 3. What were the key controlling factors of hydrological partitioning and biogeochemical transformations, and how did they regulate the spatiotemporal patterns of NO 3 -N?

Study Catchment
The study site is the 68 km 2 mixed-land-use Demnitzer Millcreek catchment (DMC), located ∼55 km SE of Berlin, Germany. The catchment experiences a typical mid-continental climate, with relatively low precipitation (∼569 mm/year) and slightly higher potential evapotranspiration (650-700 mm). Seasonal rainfall is similar in total amount but different in distribution, with more intense and convective events in summer, while more frequent, but low-intensity frontal events occur in winter. It is a flat headwater catchment (average slope <2%; Figure 1a). Streamflow is generally groundwater-dominated and strongly seasonal (winter maxima), with low runoff coefficients . Flows are sensitive to changes in climate or water management, and can even be discontinuous with channels drying for months during dry summers .
The geology is dominated by ground moraine tills in the north and glacio-fluvial deposits in mid-southern sections ( Figure 1b). Brown-earth soils with ∼70% sand characterize much of the catchment (Figure 1c), leading to high vertical hydraulic conductivity . Arable land and pasture are the main land uses, covering >60% of the area and dominating the upstream catchment, while forestry increases downstream, accounting for 36% of the area (Figure 1d). Several wetlands are distributed sporadically in the mid and south-sections of the catchment, with a major one located in the central part of catchment and traversed by the main stream (termed as "central wetland" in the following text). There are also several small dispersed urban settlements in DMC, but their impact on discharge and nutrients is limited due to the low population (∼5,000 residents) and adequate wastewater treatment facility (Wu et al., 2021).
Land management practices in the catchment have gradually changed since 1990. Following the re-unification of Germany, synthetic fertilizer applications were generally reduced and farming became less intensive in many of the German catchments (Ehrhardt et al., 2019), which was likely to have also occurred in DMC though exact data are unavailable. The channel network in northern part of catchment (including arable land and central wetland) was historically deepened and straightened for artificial drainage. These drainage practices ceased in 1990, which initially led to the gradual re-naturalization of channel characteristics (Bösel, 2018). In 2000, the central wetland was restored proactively by excavating backwaters connected to the main stream, and installing in-channel bunds to reduce channel depth and slow flow velocities. More modest, but natural wetland expansion has occurred since 2004 when beavers recolonized the catchment, with numerous dams particularly in the central wetland, increasing the area of inundation. Such changes have increased water residence time and nutrient retention in restored wetland areas Wu et al., 2021).

The mHM-Nitrate Model
The mHM-Nitrate model was recently developed by integrating a nitrogen sub-model into the multi-scale conceptual mHM platform (Kumar et al., 2013;Samaniego et al., 2010). As a fully distributed model, grid-based hydrological and nitrogen fluxes are simulated in terrestrial and in-stream phases at each timestep. Figure S1 shows the model structure and key components. For detailed model descriptions please refer to Yang et al. (2018) as only a brief summary follows.

Hydrological Sub-Model
The hydrological component of the mHM model has a flexible reservoir-based conceptual structure ( Figure S1), which has been successfully applied to multiple river basins at different scales (Samaniego et al., 2010;Thober et al., 2019). The main hydrological processes include canopy interception (Dickinson, 1984), rain-snow partitioning and snow pack dynamics (Hundecha & Bárdossy, 2004), evapotranspiration, infiltration, deep percolation, and runoff generation (Bergstrom, 1995). Runoff processes are classified as fast interflow, slow interflow, and baseflow in the hydrological sub-model, which represent overland/near-surface stormflow, shallow subsurface flow, and deep groundwater flow (GW flow), respectively. These processes are simulated in each grid.
Generated discharge is then routed in-stream to its neighboring downstream cell using the Muskingum flood routing algorithm (Gill, 1978).

Nitrogen Sub-Model
The nitrogen sub-model describes the nitrogen mass balance, including its transport and transformation under various physical fluxes and biogeochemical interactions. Fully integrated with the hydrological sub-model, nitrogen transport aligns with soil water dynamics, including infiltration into the soil, percolation to deeper soil layers and groundwater, and ultimately export to the stream. Dissolved nitrogen is categorized into two classes, the organic (DON) and inorganic nitrogen (DIN, equivalent to NO 3 -N in our study), and both forms are assumed to be fully mixed in each conceptual reservoir.
The description of nitrogen transformations is mainly adopted and modified from HYPE, a well-tested model in various water quality studies (Lindström et al., 2010). Four different nitrogen pools were established in each soil layer of terrestrial grids, known as active and inactive soil organic nitrogen (SON A , SON i ), and dissolved organic and inorganic nitrogen (DON, DIN). Within these pools, nitrogen transformations are updated at each timestep, including denitrification (from DIN to gaseous nitrogen leaving the system), mineralization (from DON/SON A to DIN), dissolution (exchanges between SON A and DON) and degradation (from SON I to SON A ). In-stream denitrification (DIN reduction in the channel network) and assimilation (from DIN to DON) are also simulated in the grids where channels occur. The initial states of all these biogeochemical processes and the input sources (wet deposition, fertilizer/manure applications, and plant residues) should be pre-defined based on the prior knowledge of the catchment.

Forcing Data Sets
The required forcing meteorological data sets are precipitation, temperature, and potential evapotranspiration.
The daily data sets of the first two were obtained from the three nearest weather stations (Lindenberg, Manschnow, and Müncheberg) operated by the German Weather Service (DWD), all within 19 km distance of the catchment. Potential evapotranspiration was estimated via the Penman-Monteith method.
The external NO 3 -N inputs consist of wet deposition, plant residual matter, fertilizer, and manure application. NO 3 -N in precipitation was defined as 2 mg/L according to the EMEP data for the period 1995(EMEP, 2001, which resulted in an annual wet deposition ranging from 8 to 17 kg/ha·yr in DMC over the past 30 yr. The remaining inputs were based on a prior survey in DMC and typical values from the similar German catchments (Yang et al., 2018). The total NO 3 -N inputs ranged from 80 to 120 kg/ha·yr in agricultural area over the study period ( Figure S2).

Calibration and Validation Data Sets
Discharge data has been collected at Demnitz Mill ( Figure 1a) since 1992 at various timesteps. The time series was constructed for model calibration and validation, by resampling the temporal resolution to weekly before 2001 and daily after 2001. Also, eight groundwater wells were established in DMC, with groundwater levels monitored since 2001 and resampled to daily timestep for subsurface storage validation.
Water quality monitoring was conducted via grab sampling at five sites along the stream (Figure 1a): the entrance and outlet of the central wetland (Peat North and Peat South), the same location as the water level gauge (Demnitz Mill), and catchment outlet (DM 27 and Berkenbrück). The NO 3 -N time series also started since 1992 and is generally weekly in temporal resolution, with several gaps at some sites due to funding limitations.

Model Parametrisation
For the description of the influential parameters and the range of all parameters please refer to Tables 1 and S1. Generally, the parametrization strategy for the hydrological simulation was adopted from mHM, including the parameter ranges for calibration. The parameters related to snow/rain partition, soil moisture distribution, and interflow generation were assigned based on simplified land use type (impervious, pervious, and forest area), while the parameters of GW flow generation were allocated based on geology classification ( Figure 1b). The other parameters remained identical in all grids (Samaniego et al., 2010). The Muskingum routing was realized based on three general parameters corrected by slope, channel length, and pervious area fraction.
As for nitrogen simulations, six parameters were used to reflect the six major NO 3 -N biogeochemical processes described above. These processes are dominated by in situ physical, chemical, and biological characteristics, which generally show high spatial heterogeneity and are sensitive to anthropogenic impacts such as agricultural fertilizer application (Lassaletta et al., 2009;Taylor et al., 2016). Therefore, the corresponding parameters were conceptualized as land-use-dependent and assigned based on the land use types (Figure 1d).

Sensitivity Analysis
Although regional parameterization grouped the spatially distributed parameters to reduce the computational demands, the total number of parameters in the hydrological and nitrogen sub-models still reached 50 and 36, respectively, which is computationally intensive. Therefore, we applied the Morris Method (Morris, 1991), also known as the elementary effects (EE) method, to identify insensitive parameters. The parameter ranking was realized according to the EE extracted from multiple perturbation trajectories, by generating each trajectory with a set of starting parameters, and perturbing each parameter p i by a variation Δ i based on a radial one-at-a-time strategy. The EE of ith parameter can be estimated as: where f denotes the objective function (root-mean-square-error, RMSE).
To achieve a better coverage of the feasible parameter space, the Latin-Hypercube sampling method was selected to generate the starting points and variation Δ. A matrix of EE (n * r) was then available for further calculation of the mean (μ i ) and standard deviation (σ i ): where EE is the EE of ith parameter in jth trajectory.
The sensitivity analysis was conducted using the SAFE tool (Sensitivity Analysis for Everybody, Pianosi et al., 2015). As mHM-Nitrate consists of two sub-models, we ranked the parameters separately for each sub-model first. Then, a joint sensitivity analysis was applied for the full version with both sub-models activated. The top 15 ranked parameters of the joint sensitivity analysis were identified as sensitive (Table 1).

Model Calibration
Given that this is the first distributed NO 3 -N modeling in DMC, all parameters were included in the calibration. Due to the expensive computational demands, the dynamically dimensioned search (DDS) was selected for parameter optimization for its approximation capacity of global optimal solutions within limited iterations (Tolson & Shoemaker, 2007). Also, we combined the Nash-Sutcliffe efficiency (NSE) and its logarithmic form for error estimation, which can effectively offset the over-sensitivity towards high values or outliers when using the conventional NSE approach: The objective function was then estimated as the weighted error of discharge (w q = 0.9) and in-stream NO 3 -N (w n = 0.1) according to Yang et al. (2018): After the initial test, the number of DDS iterations was set as 50,000 for a full convergence of optimized parameter sets. Due to the availability of discharge data, the calibration was only conducted at the Demnitz Mill from 2010 to 2019, which includes wettest and driest conditions in the past 30 yr (Wu et al., 2021). We repeated the calibration independently for five times, and manually checked the optimized parameters to ensure they reflect our local knowledge. Then the best parameter set was chosen for further analysis.
Moreover, to evaluate the representativeness of the optimized parameter set, 100,000 additional parameter sets were generated, with insensitive parameters fixed as the best parameter sets above and sensitive parameters re-sampled via Latin hypercube method. The 95% confidence bands of simulated discharge and in-stream NO 3 -N at Demnitz Mill from the top 5% parameter sets were treated as the parametric uncertainty.

Model Validation
The model was first validated at Demnitz Mill for two periods (2001-2009 and 1992-1999), with performance metrics NSE and RMSE applied to evaluate the goodness-of-fitness of both discharge and in-stream NO 3 -N simulations. Then, the sum of the water stored in all subsurface reservoirs in mHM-Nitrate (i.e., three soil layers, unsaturated and saturated zone) was calculated. Comparison between the simulated subsurface storage and measured groundwater levels from eight wells allowed the subsurface flows to be assessed qualitatively. In addition, a spatial validation for NO 3 -N simulation in stream water was conducted at Peat North, Peat South, DM27, and Berkenbrück, due to the available NO 3 -N time series (>100 records) at these sites.

Re-Optimization for Period 1992-2000
Due to unsuccessful transfer of the optimized parameter set to period 1992-2000 for NO 3 -N simulation (see below), the model was separately re-optimized for pre-2001. All the hydrological parameters were fixed and excluded from re-optimization, considering the good performance of flow simulation in this period, and to avoid overparameterization. Given the major differences between the two periods were in the agricultural area and central wetland (i.e., gradually cessation of drainage maintenance since 1990 and the wetland restoration in 2000), the most sensitive NO 3 -N parameters (rates of aquatic denitrification, soil denitrification, and soil mineralization) in arable land, pasture and wetland (deniw2-4, denis2-4, minlr2-4) were selected for optimization. 50,000 parameter sets were generated using Latin hypercube method to fully explore the parameter space. Two criteria, RMSE and percent bias, were used to evaluate the performance of these parameter sets.

Parameter Sensitivity
The sensitivity analysis was conducted separately within hydrological and nitrogen sub-models, and jointly with integrated model. The top 15 sensitive parameters from joint ranking are listed in Table 1. According to the separate sensitivity analysis (Figure 2a), parameters in the soil moisture module (which determine hydraulic conductivity and soil moisture distribution) dominated the hydrological simulations. The parameters for PET, infiltration, and percolation also impacted the hydrological simulation. In terms of separate sensitivity of NO 3 -N parameters (Figure 2b), denitrification and primary production rates were most influential, while parameters of  Table 1 for parameter explanation). Increases in the X and Y axes indicate a stronger influence on model results and interdependence with the other parameters, respectively. the remaining transformation processes (mineralization, degradation, and dissolution) only showed minor effects. The joint sensitivity ranking indicated the dominant role of hydrological parameters in the overall NO 3 -N simulation (Figure 2c), while only two NO 3 -N parameters (in-stream denitrification and soil denitrification in arable land) were identified as sensitive.

Model Performance at the Calibration Site
The hydrological simulations at Demnitz Mill performed well in all three periods, with NSE >0.6 and no deterioration in performance for validation periods (Figure 3; Table 2). The NO 3 -N simulations also successfully captured the seasonal and long-term dynamics of in-stream NO 3 -N at Demnitz Mill after 2000. The simulated uncertainty range derived from Latin hypercube samplings bracketed the best optimized simulation from DDS, and most of the discharge and NO 3 -N observations, further giving us confidence in the robustness of the model calibration schemes. However, the optimized parameter set failed to captured in-stream NO 3 -N at Demnitz Mill during 1992-2000 (though flow dynamics were well simulated), with significant underestimation of almost all high peaks and low NSE value (0.50; see Figure S4).

Spatial Validation of Subsurface Storage and In-Stream NO 3 -N at Uncalibrated Sites
The GW level monitored at 8 sites since 2001 provided an opportunity to evaluate the robustness of subsurface hydrological simulations. From a qualitative perspective, the simulated subsurface storage changes matched the long-term dynamics of measured groundwater levels, including the regular inter-annual variability and the significant increase during the wettest period (2011-2012; Figure 4). The spatial distribution of simulated subsurface storage also showed a high affinity with the data from the corresponding GW wells (average R > 0.7, p < 0.01; Table S3). The only exception was GW7, located near a small inundated area, which could not be captured by the model, resulting poorer correlation.
Spatial validation was also applied for in-stream NO 3 -N simulation, using the in-stream NO 3 -N records from 5 monitoring sites across the catchment (Figure 1a). The NSE of the NO 3 -N simulations at these sites ranged from 0.54 to 0.67 (Table 3), indicating that the model captured the spatial variability of   Mill During 2010(Calibration, c), 2001(Validation 1, v1) and 1992-2000 in-stream NO 3 -N reasonably well. The lowest NSE was found at DM27, which could be explained by the limited number of samples there.

Re-Optimization for Period 1992-2000
To identify a more appropriate parameter ranges for NO 3 -N simulation in 1992-2000, 50,000 random parameter sets were evaluated using two criteria RMSE and percent bias. As is shown in Figure 5, the model performance was sensitive to both aquatic and soil denitrification rates, while mineralization rates had nearly no impact. Among the denitrification parameters, those in arable land (deniw2, denis2) showed more dominant impacts, while the rest were less important (deniw3,4 and denis3,4). Although the domination of deniw2 and denis2 masked the rest parameters to some extent, it is clear that decreasing denitrification rates can lead to better model performance for this earlier period. Further, the parameter PDF against RMSE suggested the ideal parameter ranges of denis2 and denis2. So accordingly, the denitrification rates in arable land, pasture, and wetland were manually lowered, while mineralization remained unchanged (see re-optimized values in Table S2). The comparison of parameter sets performance before and after optimization is shown in Figure S4. Significantly improved performance was achieved at Demnitz Mill with the new parameter set, as NSE increased from 0.50 to 0.66 (Table 4). The simulated in-stream NO 3 -N now captured both low concentrations in the growing seasons and high peaks in winter. The NSE at other sites also reached >0.65 (Table 4), which validated the spatial representativeness of the optimized parameter set.
Given the reasonable validation performance of hydrological and NO 3 -N simulations, the outputs generated by the original parameter set for post-2001, along with the re-optimized parameter set for pre-2001 were integrated and used for further analysis.

Spatial Distribution of Hydrological Fluxes and Biogeochemical Transformations
For an overview of spatially distributed model outputs, the main hydrological fluxes, corresponding NO 3 -N concentrations and biogeochemical transformations were averaged over the whole period . As shown in Figure 6, most hydrological fluxes exhibited moderate spatial variability, mainly depending on the vegetation cover. As the dominant hydrological processes in DMC, actual evapotranspiration (AET) were up to 42 mm/yr higher from forestry areas (Figure 6a). Consequently, reduced vertical hydrological transport was simulated compared to the non-forested areas, which is covered by crops, grass, and shrubs. Infiltration and subsequent interflow generation were 29.0% and 31.4% lower respectively in southern part of DMC (mostly underlain by forest) compared to the northern part (mostly underlain by crops; Figure 6b). GW flow exhibited a similar spatial pattern, though its parameterisation was based on catchment geology (Figure 1b). The geological impact was masked by the strong differences of water recharge from soil layers with different vegetation covers, which further led to the vegetation-dependent distribution of GW flow (27.8% lower in non-forestry area).  The patterns of NO 3 -N concentrations and biogeochemical processes were also highly land use-dependent. Figure 6d shows the comparison of NO 3 -N inputs under different land use. In arable land and pasture with fertilizer/ manure applications, the inputs reached up to 120 kg/ha·yr, while those in forest and wetland, which merely consist of plant residues, were limited <25 kg/ha·yr. Highly dependent on the available pools, biogeochemical transformations mainly followed the spatial pattern of NO 3 -N inputs, with ∼30-50 kg/ha·yr of soil mineralization and denitrification in agricultural areas compared to <10 kg/ha·yr in the non-agricultural areas. Aquatic denitrification was different, as it is determined by both discharge and available NO 3 -N, increasing along the main stream and reaching a maximum near Demnitz Mill. Despite the significantly stronger denitrification in northern part of DMC, the discrepancy on NO 3 -N inputs in different land uses was not offset. Consequently, higher NO 3 -N concentrations in interflow and GW flow were simulated in the agricultural area, with respective average values over 30 yr of 10.6 and 11.9 mg/L in agricultural area and 2.5 and 0.5 mg/L in the remaining predominantly forested area.

Temporal Dynamics of Water and NO 3 -N Balances
The annual water and NO 3 -N balances in DMC were calculated based on the average values of all grids (Figure 7). The water balance generally followed precipitation inputs, which showed significant increases in storage and groundwater recharge in wet years when annual precipitation exceeded 600 mm (1993, 1998, 2002, 2007, 2010, and 2017), while subsurface storage remained relatively constant or moderately decreased in the remaining years ( Figure 7b). The strongest replenishment of subsurface storage occurred during 2010 (113.5 mm). But this was followed by a prolonged dry period, and a gradual decrease of subsurface storage to a deficit of 190.8 mm by the end of 2016. AET was the dominant component that contributed to the water balance (>80%), while runoff generation always accounted for <25% (Figures 7a and 7b). The temporal variability of all hydrological processes was also aligned with annual precipitation, as significant positive correlations were found between annual precipitation and AET/interflow/GW flow. However, contributions from three hydrological components varied due to varying inter-annual wetness. for example, AET accounted for only 67%-82% of annual precipitation in wet years, while in dry years such as 2014-16, simulated AET could exceeded precipitation. In contrast, a higher proportion of precipitation contributed to interflow generation in wet years and the following year. GW flow also followed a similar long-term pattern to interflow but with a significant lag, as the peaks of GW flow generation occurred in 1-2 yr after precipitation peaks. Accordingly, the higher proportion of runoff generation did not occur in wet years but one year afterwards, ranging from 21.5% to 34.4% of precipitation inputs.
Biogeochemical transformations also followed inter-annual wetness conditions. The only exception is plant uptake, which remained independent from annual precipitation (31.2-45.3 kg/ha·yr) and accounted for 65.0%-73.8% of the total NO 3 -N inputs. In contrast, denitrification was negatively correlated with annual precipitation. In wet years 2007, 2010, and 2017, denitrification was simulated as 42.5, 39.6 and 39.4 kg/ha·yr and accounted for 33.2%-37.7% of annual inputs, while in dry years 2014 and 2015, NO 3 -N loss by denitrification was only 19.0 and 19.5 kg/ha·yr (<17% of annual inputs). The

Table 4 Performances of NO 3 -N Simulation With Original and Modified (*) Parameter Sets
NO 3 -N export via interflow and GW flow generally followed the long-term trends of two runoff components. Interflow export exhibited higher losses in wet years or one-year afterwards, while GW flow export generally showed a 1-2 yr lag. Compared to export in GW flow, losses in interflow were highly variable, and could reach up to 16.0 kg/ha·yr in wet year like 2010 but drop to 1.0 kg/ha·yr in dry year like 2016. Simulated GW flow exports, however, remained at a relatively stable level (4.6-10.9 kg/ha·yr), mainly depending on the hydrological fluxes.
Moderate differences of all NO 3 -N components were found between pre-2001 and post-2001, which could be attributed to the lower denitrification, and higher potential of interflow export prior to wetland restoration. for example, denitrification was simulated as 24.1 and 39.6 kg/ha·yr in year 1998 and 2010, despite annual precipitation being close. Similarly, different interflow export potential can be found between 1993 and 2010, when NO 3 -N export via interflow was similar (16.0 kg/ha·yr in both years), whereas interflow fluxes were very different (107.2

Distinct Long-Term Schemes in Agricultural and Forest Areas
Based on the temporal dynamics (Section 3.3.1) and spatial patterns (Section 3.3.2) of simulated outputs, landscape characteristics and hydroclimatic variability were identified as two main sets of controlling factors on hydrological and nitrogen fluxes. To further reveal their impacts in a more disaggregated way, the model outputs were differentiated based on agricultural/forest areas (see Figure S3 for area selection) and extreme wet/dry periods (2010-11 and 2014-15 respectively). Four distinct conditions are shown in Figure 8.
Consistent with the conclusion above, NO 3 -N fluxes in all components, including soil layers and deeper aquifer, exhibited higher magnitudes under arable land compared to forest. This was evident in both wet and dry periods, indicating the dominant role of NO 3 -N inputs. Similar to the spatial patterns in Section 3.3.1, lower infiltration, percolation, and subsequent interflow/GW flow generation were simulated under forest, regardless of the wetness condition. However, the effects of these hydrological differences on the NO 3 -N distribution are secondary to (partly masked by) the background concentrations in arable land and forest. In other words, landscape characteristics (NO 3 -N inputs and storage), rather than hydrological transport, dictated the spatial pattern of background NO 3 -N values regardless of climatic change.
Hydrological transport, however, had a dominant influence on temporal dynamics -more specifically, during the transition between wet and dry periods in both arable land and forest. In dry periods, the vertical hydrological fluxes (i.e., infiltration through the soil layers and percolation into deeper aquifer) were significantly depressed, leading to a decrease in runoff generation (76.6% and 78.6% under arable land and forest respectively). Additionally, drought conditions profoundly altered the sources of runoff generation. In wet periods, interflow was more dominant, accounting for 60.0%-67.3% of the runoff, while in dry periods, GW flow was dominant and the contribution of interflow dropped to 31.1%-46.0%. The reduced vertical hydrological fluxes correspondingly changed the NO 3 -N distribution in the soil layers and generated runoff, by limiting NO 3 -N transport into deeper soil layers and the underlying aquifer. In Section 3.2.2, we saw that NO 3 -N storage increased in dry years (Figure 7b). Here, Figure 8 shows more clearly that the storage and inputs were restricted to the first two soil layers. In the dry years of 2014-15, simulated NO 3 -N storage in first two soil layers increased by 34.0% and 57.7% in arable land and forest, respectively, and the NO 3 -N concentration in soil layer 1 in arable area reached up to 83.1 mg/L. In contrast, 44.0% lower NO 3 -N concentrations in soil layer 3 were simulated due to the lack of NO 3 -N supply from upper layers. This further led to a 43.4% decrease of NO 3 -N concentration in generated interflow.

Credibility of Modeling Results During 2001-2019
Spatially distributed environmental models have stimulated interest over the past few decades, due to the increasing availability of distributed hydrometeorological/geospatial data and advances in computational resources (Tang et al., 2007;Wellen et al., 2015). Though greater spatial details in models are feasible, model users need to carefully evaluate the results, as the increasing complexity of such models poses challenges in terms of greater potential for overparameterization and uncertainty from model inputs and structure (Tonkin & Doherty, 2005). In our modeling application, good performance was achieved in 2001-2019 with calibration based on hydrological and water quality data from Demnitz Mill. This is similar to many other studies, where calibration is only based on observations at the catchment outlet (Wellen et al., 2015). However, as already noted by many researchers (Fuamba et al., 2019;Moussa et al., 2007;Uhlenbrook & Sieber, 2005), single-gauge-based calibration can be misleading (e.g., higher NSE) compared to calibration to multiple gauges, because equifinality can result in significant but unidentified compensating errors in distributed grids (Tang et al., 2007;van Werkhoven et al., 2008). To assess the spatial representativeness of our modeling, we used in-stream NO 3 -N from 5 distributed sites for validation. Such validation strengthened confidence in the simulation results, as both the high peaks and low values were successfully captured at all sites after 2001 ( Figure S4). In other words, we demonstrate that mHM-Nitrate has the capacity to reproduce reliable spatial distribution of in-stream NO 3 -N based on the observation of a single catchment outlet gauge, and can be a useful tool for sparsely monitored catchments. This is especially meaningful, considering that data is lacking at most sites, and multiple-gauges calibration is usually impossible outside specific research catchments (Moreau et al., 2013).
Despite the apparent success of in-stream NO 3 -N simulations, uncertainty over internal processes remains as a result of data lack for quantitative verification. The main sources of uncertainty come from two aspects: uncertainty over input data and parameter compensation within the model. The uncertainty of inputs (e.g., N-sources) has long been emphasized by the modeling community, given its obvious impact on simulations (Tang et al., 2007;Wriedt & Rode, 2006). Unfortunately, precise records of fertilizer inputs were unavailable for DMC, though the sources and amounts of NO 3 -N inputs in our study were comparable to many modeling applications in similar agricultural landscapes in Europe (e.g., Leip et al., 2011;Yang et al., 2018). The second source of uncertainty, the compensation effects of parameters representing internal processes, is common in distributed modeling applications (Anderton et al., 2002;Arnold et al., 2015;Van der Velde et al., 2010). This is also likely to affect our results given the total number of 86 parameters in mHM-Nitrate. In other words, it is possible that the overestimation of one process compensates the underestimation of other processes, which leads to an overall "good fit" at stream gauges . An alternatively way to constrain the uncertainty is incorporating "soft information" in model assessment Fuamba et al., 2019). In our study, groundwater levels from 8 sites were used to qualitatively evaluate the changes in subsurface storage. The dynamics of simulations and observations matched spatially and temporally (Figure 4), increasing our confidence on subsurface simulations.
To further assess representation of internal processes, comparison with other modeling studies was instructive. In DMC, flow dynamics have been recently simulated with the tracer-aided distributed hydrological model EcH 2 Oiso (see model descriptions in Kuppel et al., 2018). Using water isotopes as tracers and including distributed soil moisture data as calibration targets, the internal hydrological fluxes were simulated with greater model complexity and robustness (Smith et al., 2021). The very similar spatial patterns of AET and discharge from the two models demonstrated the key hydrological fluxes estimated by mHM-Nitrate were reasonable, though we cannot directly compare other internal processes due to the different conceptualization and model structure. For NO 3 -N modeling, our study is the first application in DMC, so comparison is only possible with previous studies in analogous catchments with similar latitude (Table 5). Generally, the NO 3 -N leaching aligned to the agricultural (especially arable) land use, and the leaching rates in DMC were on the same order of magnitude with sites sharing similar proportions of agricultural land cover (e.g., Hesser et al., 2010;Yang et al., 2018). The exception is significantly higher leaching in Rode et al. (2009), which is probably due to dense artificial drainage aiding both mineralization and transport. Similar stimulated NO 3 -N leaching by artificial drainage was also reported in Hu et al. (2007), where agriculture accounted for >90% of the catchment area and NO 3 -N losses reached up to 29 kg/ha·yr. Denitrification rates were unavailable in several studies because they were either included Note. Agri. and Precip. denote agricultural proportion and annual precipitation, respectively.

Table 5
Reported  and Denitrification (Denitr.) Rates in Modeling Applications (From Catchment to Regional-Scale) in terrestrial NO 3 -N retention or not reported. But the available reports range from 19 to 24 kg/ha·yr, which is consistent with our simulated value (21 kg/ha·yr; Table 5). For mineralization, the reported values span a wide range, as the process is closely coupled with soil properties and land management; for example, 11.5 (Kopáček et al., 2013), 16 (Devito et al., 1999), 65-67 (Ferrant et al., 2011), 115-139 (Devito et al., 1999, and 4-104 (Andersson et al., 2002) kg/ha·yr. Our results are in the same order of magnitude and close to some studies (e.g., Yang et al., 2018), which represent similar characteristic are they are in the same climate and geomorphic unit in NE Germany. This suggestd the mineralization has been reasonably simulated. However, uncertainty needs to be considered as no direct observation is available in DMC for validation.

Limitations in NO 3 -N Simulation Pre-2001
Though the initial calibrated parameter set performed well after 2001, it performed less well in reproducing the in-stream NO 3 -N prior to 2001. The poorer performancepotentially originates from the lack of physics-based description of some biogeochemical processes (e.g., dissimilatory nitrate reduction to ammonia [Rogers et al., 2021] or vegetation uptake ), as they were all highly conceptualized and simplified in mHM-Nitrate. Further, the change in management practices might account for another limitation in parameter transferability, given the gradual cessation of drainage since 1990 and wetland restoration in 2000 in DMC Wu et al., 2021). Many previous studies have highlighted the potential impacts of management on nitrogen transport (Arheimer & Wittgren, 2002;Gärdenäs et al., 2006;Olesen et al., 2019;Taylor et al., 2016).
Ideally, these management practices should be conceptualized and integrated into the model for the pre-2001 simulation; however, conceptualizing and upscaling these practice-related processes are difficult and detailed data would be required (e.g., the location and fluxes of drains), which are generally unavailable at the catchment scale (Rode et al., 2009). Consequently, detailed mechanistic descriptions have often been limited to field-scale modeling (e.g., tile drainage in Gärdenäs et al., 2006), while for catchment-scale modeling, they were simplified as changing the land use, defining a simple new module or re-parameterizing the relative processes (Arheimer & Wittgren, 2002;Olesen et al., 2019;Rode et al., 2009;Taylor et al., 2016). In this study, we used the latter solution: recalibrating the parameters in terms of soil denitrification, aquatic denitrification, and mineralization in arable land, pasture and wetland, to represent the impacts of drainage and wetland restoration. By revealing the PDF of parameter values against model performance ( Figure 5), 50,000 iterations were enough to indicate how these values should be adjusted for period pre-2001, that is, lower soil and aquatic denitrification rates. This matched our local knowledge in DMC Wu et al., 2021), as both reducing drainage and wetland restoration have been shown to increase stream water residence times, therefore increasing the potential for denitrification and NO 3 -N retention (Hansen et al., 2018). It is important to recognize that following the re-unification of Germany, synthetic fertilizer applications were generally reduced in many agricultural areas in the former East Germany (Ehrhardt et al., 2019). Considering the abnormally high NO 3 -N observations before 1996 (Figures 3  and S4), reduction on fertilizer inputs since then was also likely in DMC and may affect our results. However, here we did not change the fertilizer inputs during optimization due to the lack of accurate records, and to prevent introducing more uncertainties. Such correspondence between a modeler's local knowledge and parameter adjustment are necessary, as overparameterization is common for distributed models due to the high-dimensional, nonlinear parametric spaces (Tang et al., 2007).

Controlling Factors in Spatiotemporal Pattern of NO 3 -N
From the analysis of the simulation results, we identified landscape characteristics, hydroclimatic variations and management practices as three controlling factors for NO 3 -N patterns in DMC. Although these factors have been well documented elsewhere (Inamdar & Mitchell, 2006), their site-dependent characteristics means that transferring knowledge of the relative influence of such factors between different catchments is usually difficult (Musolff et al., 2015). However, by fully incorporating relevant spatial information, distributed modeling provides a framework for differentiating these controlling factors and evaluating them in a quantitative manner, and a fine level of spatial resolution.
From the spatial perspective, vegetation cover was the major factor on catchment hydrology, controlling the depressed fluxes in interflow and GW flow in forest compared to non-forest areas ( Figure 6). This was mainly attributable to the higher ET from forests, also consistent with other modeling applications in DMC (see enhanced transpiration in forest in Smith et al., 2021). However, the heterogeneity of fluxes was modest, and only had a minor impact on NO 3 -N patterns. In contrast, NO 3 -N inputs dictated the spatial distribution of NO 3 -N in all soil layers, aquifers and streams ( Figure 6). For in-stream NO 3 -N draining arable land, the concentrations have remained >10 mg/L over much of the past 30 yr, close to the maximum values suggested by the Drinking Water Directive (Council of the European Union, 1998). Similarly, high concentrations have been reported in many surface water bodies in Europe due to extensive diffuse sources (Bøgestrand et al., 2005;Lassaletta et al., 2009;Velthof et al., 2009). The high N surplus (60 kg/ha·yr in Europe, see Leip et al., 2011) has not only contributed to surface water pollution, but also extensively contaminated GW bodies (European Environment Agency, 2012). As a groundwater-dominated catchment, GW flow in DMC accounted for 40%-70% of runoff generation and its NO 3 -N concentration remained stable (∼12 and 0.5 mg/L in agricultural and forested areas respectively; Figure 8). In other words, NO 3 -N storage in GW reservoir dictated the high background concentrations in arable land due to historical fertilizer inputs. The recovery from the contaminated GW is generally slow due to the long residence times (Tesoriero et al., 2013), which can exceed >25 yrafter intensive management (Pérez-Martín et al., 2016). To note is that, identifying NO 3 -N inputs as the major landscape characteristic for NO 3 -N concentrations in DMC was based on the fact that the soil properties (i.e., sandy and clay proportions) were relatively uniform. In other catchments, subsurface geology and soil properties can be the dominant factors as well (e.g., sandy proportion in Rode et al., 2009).
Hydroclimate only had a minor influence on the overall spatial NO 3 -N pattern due to the near-uniform distribution of precipitation. However, it was the main driver of the temporal dynamics; more specifically, the wetness conditions determined the inter-annual variability of NO 3 -N by regulating the transport capacity. The four schemes summarised in Figure 8, show how NO 3 -N was retained in the soil zone due to the lack of vertical hydrological connectivity in dry periods, and led to a strong decrease in interflow concentrations. The importance of hydrological connectivity in solute transport has been well documented (Bechtold et al., 2003;Hesser et al., 2010;Ocampo et al., 2006). For example, Hesser et al. (2010) found that NO 3 -N dynamics could only be simulated satisfactorily with reasonable simulation of key contributing runoff components. Not only the solute transport, the dominance of hydrological conditions also largely regulated biogeochemical transformations in DMC, demonstrated by the sensitivity ranking of model parameters (Figure 2). Despite soil denitrification and mineralization accounting for 20%-30% of NO 3 -N inputs, their temporal dynamics were mainly controlled by the wetness conditions and the resulting influence on catchment hydrology. On the one hand, the landscape connectivity determined how much NO 3 -N can be transported vertically and laterally, and therefore where N is stored and available for further transformations (Bechtold et al., 2003). On the other hand, all the transformation processes were intrinsically correlated to the soil moisture distributions (Seitzinger et al., 2006). Similar to other catchments (Hesser et al., 2010;Rode et al., 2009), NO 3 -N export in dry years was substantially lower than in wet years (91.5% less in DMC), as an integrated result of depressed hydrological fluxes and biogeochemical processes. However, dry years were still important in the longer term, as nitrogen accumulation in the soil during these periods enhanced NO 3 -N concentrations and fluxes in the subsequent years as wetness increased (e.g., 2002 and 2017; Figure 3). These accumulation-flushing sequence of NO 3 -N observed in DMC (Wu et al., 2021) have been reported for many other catchments (Bechtold et al., 2003;Dupas et al., 2016;Outram et al., 2016). Considering that the risk of NO 3 -N peaks can even be higher due to projected (Papadimitriou et al., 2016) or observed  climate change (more frequent and severe droughts), further understanding of hydrological and biogeochemical re-coupling in re-wetting periods is needed.
Management practices for stream and wetland restoration were the remaining influence on the NO 3 -N dynamics in the study period. Their incorporation in the modeling results is "implicit" here because they were not conceptualized by any model equation but by recalibration to the period before management intervention. By comparing the simulated results between the original and recalibrated parameter sets, we can roughly estimate that the NO 3 -N export would be ∼10% higher if the cessation of tile drainage and wetland restoration were not applied. Considering the restored wetland accounted for 2.5% of the catchment area, such performance was not ideal comparing to the previous observations or simulations on restored wetlands. For examples, a 6% NO 3 -N reduction was achieved by restoring only 0.4% of the catchment area as wetlands in a Swedish catchment (Arheimer & Wittgren, 2002), and a 5% increase of wetland areas was simulated to reduce ∼30% NO 3 -N export in Mississippi River basin (Hansen et al., 2018). Given the high NO 3 -N budgets in arable land, further management to enhance wetland areas is a potential option for reducing surface water pollution. Moreover, reducing agricultural inputs is key to achieving water quality targets. As McLellan et al. (2015) point out, effective reduction of NO 3 -N export is unlikely without reducing fertilizer inputs and removing large areas of land from agricultural production. To aid the development of more sustainable management strategies, coupling catchment nitrogen models with economic models is recommended, as it can evaluate alternative management practices from both NO 3 -N reduction and economic loss perspectives (Peña-Haro et al., 2009). Thus, spatially distributed models like mHM-Nitrate have important potential in landscape analysis to support such management.

Conclusions
Nitrogen transport and transformations at the field scale are quite well understood, but our ability to quantify NO 3 -N fluxes and export at the catchment scale remain limited due to spatially heterogeneous biogeochemical-hydrological processes that are dynamic at a range of temporal scales. To unravel the coupling of NO 3 -N cycling and hydrological processes in a more explicit and quantitative way, a grid-based model, mHM-Nitrate, was applied in a mesoscale catchment (Demnitzer Millcreek, DMC) near Berlin. Calibration against observations at a single gauge were successfully validated for the 2001-2019, and the spatial outputs were further validated against in-stream NO 3 -N concentrations from 5 sites. However, changes in management practices led to degraded performance in pre-2001. A recalibration was conducted for this period based on Latin-Hypercube sampling and our prior site knowledge, which support an increased role for instream de-nitrification following restoration.
Overall, mHM-Nitrate can adequately simulate hydrological and NO 3 -N fluxes with relatively limited data, but it is sensitive to dramatic change of management practices.
Further analyzing the long-term hydrological and NO 3 -N fluxes, we identify landscape characteristics, hydroclimatic variability, and management practices as three dominant controlling factors in DMC. (a) Landscape characteristics, that is, vegetation cover and NO 3 -N inputs, dictated the spatial patterns of flow and NO 3 -N. (b) Hydroclimatic variability, that is, changing inter-annual wetness conditions, determined the hydrological transport capacity, which further drove the NO 3 -N dynamics. The simulated NO 3 -N in soil layers showed that NO 3 -N was retained in soil zone in dry periods due to the depressed vertical transport, therefore leading to lower concentrations in generated runoff. Consequently, inter-annual NO 3 -N dynamics exhibited an accumulation-flushing mechanism in dry-wet years. The sensitivity analysis also revealed that biogeochemical transformations were closely related to hydrological transport. (c) Management practices, that is, the cessation of tile drainage and wetland restoration, reduced NO 3 -N export by 10%. Due to the high potential of NO 3 -N peaks, our study highlighted the need for further research on rewetting periods. It also demonstrated that further management practice will be needed to reduce NO 3 -N losses in DMC, including reducing agricultural inputs, taking land out of production and/or increasing restored wetland areas.

Data Availability Statement
The data used are presented in the tables, figures and Supporting Information. For source codes please refer to https://doi.org/10.5281/zenodo.3891629.