The effect of rainfall variability on Nitrogen dynamics in a small agricultural catchment

Throughout history, extreme storms and droughts have had serious impacts on society and ecosystems globally. Rainfall variability particularly has been identified as a primary manifestation of climate change. However, so far little has been done to explore the effect of rainfall variability on water quality. This study aims to investigate the effect of rainfall variability on nitrogen (N) dynamics and its potential influence on water quality. The transport of water and nitrate was simulated for a small agricultural catchment in central Germany using the fully coupled surface-subsurface model HydroGeoSphere. Rainfall time series with specific climatic characteristics were generated using a stochastic rainfall generator. N transformation and transport were compared for four scenarios (wet, normal, dry and extremely dry conditions) in order to identify the impact of inter-annual rainfall variability. The results suggest that higher annual precipitation enhances N transformation and transport, whereas lower annual precipitation is conducive to the N retention capacity. The retention capacity declines severely when vegetation suffers from drought stress, suggesting that vegetation plays a vital role in the response of N dynamics to extreme droughts. The linear regressions between selected parameters of the rainfall generator and N loads/fluxes were analyzed to elucidate the impact of intra-annual rainfall variability. The results indicate that wet/dry conditions and different dry–wet patterns caused by storm duration distributions and inter-storm period distributions can significantly affect N loads and in-stream nitrate concentration, respectively. In the warm season, droughts prompt the accumulation of soil organic nitrogen (SON), but drying-wetting cycles can enhance the extensive transformation of SON. In-stream nitrate concentration dramatically elevates during the rewetting period after a drought. Wet/dry conditions and patterns of precipitation intensity within a year determined by seasonal average rainfall intensity and the probability of drizzle events prominently alter mineralization and plant uptake. Elevated mean rainfall intensity has a small effect on stream water quality. Therefore, mineralization and plant uptake are critical processes governing N dynamics and the impact of nitrate on water quality under varying rainfall conditions. Overall, the study clarifies the effect of rainfall variability on N dynamics in a small agricultural catchment, which provides theoretical support for formulating fertilization strategies and protecting aquatic ecosystems under climate change.

1 Introduction

The hydrological processes are susceptible to meteorological conditions on various spatial and temporal scales (Ionita et al., 2017; Laaha et al., 2017; Zhang et al., 2021). In the past decades, extreme climate events intensified by human-induced climate change have frequently occurred globally (Pall et al., 2011; Min et al., 2011; Williams et al., 2015; Hari et al., 2020), most of which caused water scarcity and poor water quality at regional scales (Zwolsman and van Bokhoven, 2007; Delpla et al., 2009; Whitehead et al., 2009; Stahl et al., 2016; Ballard et al., 2019; Bauwe et al., 2020; Geris et al., 2022). As predominant extreme climate events worldwide, heavy rainstorms and severe droughts share a common characteristic of rainfall variability (Trenberth, 2011; Pendergrass et al., 2017; Hanel et al., 2018). In the context of global warming scenarios, anthropogenic amplification of rainfall variability has been identified (Zhang et al., 2024). Thus, the effect of rainfall variability on water resources has attracted much attention around the world.

During heavy rainstorms, extraordinary rainfall amounts and intensities cause large amounts of rainwater to infiltrate into soils and trigger flash floods in a short time. Increased groundwater flow and enhanced surface runoff stimulate the movement of solutes retained in the soil, which thereby can lead to water quality degradation (Geris et al., 2022). Different from heavy rainstorms, severe droughts driven by precipitation deficits occur during several months and potentially years (Otkin et al., 2018), from which it takes 1–2 years for hydrological components to recover (Hanel et al., 2018). Wilusz et al. (2017) decomposed the relationship between rainfall variability and transit times, a reflection of water velocities that control solute transport, and illustrated that climate change may seasonally alter the ages of water in streams and thereby influence water quality in the future.

Nitrate (NO${}_{\mathrm{3}}^{-}$-N) is a major solute threatening the quality of drinking water and destroying the structure and functions of aquatic ecosystems (Vitousek et al., 2009; Alvarez-Cobelas et al., 2008; Dupas et al., 2017). The nitrate turnover processes established at the catchment scale are expected to change due to climate change (Whitehead et al., 2009; Hesse and Krysanova, 2016; Mosley, 2015), especially due to increasing drought events (Zwolsman and van Bokhoven, 2007; Ballard et al., 2019; Zhou et al., 2022; Winter et al., 2023). Several extreme rainfall events and droughts occurred in Central Europe during the last two decades, which has attracted widespread attention (Ulbrich et al., 2003; Fink et al., 2004; Orth et al., 2016; Thieken et al., 2016; Hanel et al., 2018; Hari et al., 2020; Voit and Heistermann, 2024). Zhou et al. (2022) detected higher soil N surplus (total N input with the crop/plant uptake subtracted) and a decrease in terrestrial N export in agricultural areas located in central Germany during the drought years (2015–2018). The same phenomenon reported in the Nitrate Report 2020 of the Netherlands (RIVM – National Institue for Public Health and the Environment, 2021) indicates that more N was retained in the soil during the drought period compared to the pre-drought period. Notably, the 2018–2019 (consecutive) drought triggered unprecedented tree mortality across multiple species in Central European forests, accompanied by unexpectedly persistent drought legacy effects (Schuldt et al., 2020), from which Winter et al. (2023) drew the conclusion that severe multi-year droughts can reduce the nitrogen (N) retention capacity of catchments. These appear to be opposite conclusions, which can be attributed to different investigation timescales. The former study compared N export between drought years and the pre-drought period. The latter considered the subsequent rewetting period, when most nitrogen accumulated during the drought left the catchment. Leitner et al. (2020) also found that in the year after a summer drought, NO${}_{\mathrm{3}}^{-}$ leaching via soil water seepage was significantly elevated compared to the long-term mean in a temperate mixed forest on karst, which was investigated in wetland-influenced catchments as well (Watmough et al., 2004). These studies demonstrate that rainfall variability profoundly affects N dynamics at both inter-annual and intra-annual timescales. Therefore, it is imperative to shed light on the impact of rainfall variability on water quality in terms of N dynamics.

To fill the gap, the present study explored the impact of rainfall variability on N dynamics and its potential influence on water quality across inter-annual and intra-annual timescales. To characterize rainfall variability, a stochastic rainfall generator (Robinson and Sivapalan, 1997) was employed to generate rainfall time series with different climatic characteristics. The research was conducted in a small agricultural catchment located in central Germany, where a hydrological model was previously established utilizing the fully coupled surface-subsurface numerical simulator HydroGeoSphere (Yang et al., 2018). The framework of N dynamics was modified from the ELEMeNT approach (Exploration of Long-tErM Nutrient Trajectories; Van Meter et al., 2017). The research is divided into two main components. First, three representative years (with high, normal, and low annual precipitation amounts, respectively) were chosen from the past two decades in central Germany as the target scenarios. A fourth scenario with low annual precipitation amounts coupled with reduced plant uptake represents a case where vegetation is partially destroyed by extreme droughts. The statistical analyses of N loads and fluxes and a comparison across different scenarios were conducted to reveal the influence of inter-annual rainfall variability. Second, rainfall time series generated by adjusting the values of specific rainfall-generator parameters one at a time (with other parameters held constant) replace the observed rainfall data for the simulation period and are used to drive the flow and nitrogen transport models. The responses of the N loads and fluxes to the parameters (e.g., the amplitudes of the seasonal variations in the storm duration and inter-storm period) were thoroughly analyzed to clarify the effect of intra-annual rainfall variability. The study will provide theoretical support for formulating fertilization strategies and protecting aquatic ecosystems in the context of climate change.

2 Data collection

2.1 Study site

The studied catchment, Schäfertal, is located on the lower reaches of the eastern Harz Mountains, central Germany, with an area of 1.44 km². Since 1968, this first-order catchment has been subject to broad hydrogeological investigations, analyses, and research (Altermann and Mautschke, 1970; Altermann et al., 1977; Borchardt, 1981; Altermann and Steininger, 1994; Wollschläger et al., 2016). The valley bottom contains riparian zones with pasture and a small channel (Fig. 1a). The hillslopes on both sides of the channel have an average slope of 11°, mostly used as farmland. The farmland undergoes intensive agriculture, primarily winter wheat growing (Yang et al., 2018). Two small portions near the western edge are mostly forest. The land use types generally remain stable unless economic and ecological goals change, e.g., shifting between planting and pasture (Wang et al., 2023). A meteorological station, which is 200 m from the catchment outlet, provides records of precipitation, air and soil temperatures, radiation, and wind speed. The catchment outlet in the stream, where a gauge station was built, is considered the sole exit that allows water and solutes to leave the catchment and enter the downstream catchments. This is because a subsurface wall (∼55 m long and ∼7 m deep) was erected across the valley to block subsurface flow. The gauge station provides discharge data at 10 min intervals, aggregated to daily means in this study. Nitrate concentration data were measured by sampling near the gauge station at 14 d to monthly intervals (Dupas et al., 2017), covering the period 2001–2010.

The aquifer is thin, with the thickness varying from ∼ 5 m near the valley bottom to ∼2 m at the top of the hillslopes (I-I cross section, Fig. 1a). Thirteen wells, each ∼2 m deep, were constructed and fitted with automated data loggers to record groundwater levels. The groundwater levels exhibit pronounced seasonal variations, rising to the land surface during winter and receding to depths of ∼3 m below the ground during summer. The groundwater storage is low (∼500 mm, stored water volume divided by catchment area). Most of the groundwater converges towards the channel vicinity, with the upper sections of the hillslopes typically in an unsaturated state (Yang et al., 2018). Luvisols and Gleyic Cambisols are the aquifer materials of the hillslope. The valley bottom is dominated by Gleysols and Fluvisols (Fig. 1b) (Anis and Rode, 2015). Generally, the aquifer is comprised of two horizontal layers: a top layer of approximately 0.5 m thickness, with higher permeability, higher porosity, and the developed root zone from crops; and a base layer with less permeability due to the high loam content (Yang et al., 2018). The bedrock underlying the aquifer is comprised of greywacke and shale (Graeff et al., 2009). Owing to the aquifer material and unknown hydraulic properties, the bedrock is often regarded as impermeable.

Figure 1 (a) The catchment “Schäfertal” indicated by the red line (background image from © Google Maps), with a cross-sectional view for the flow and saturation (Yang et al., 2018). (b) The distribution of soil type in the catchment. (c) The measured daily precipitation (P), discharge (Q), and the simulated actual evapotranspiration (ET) (Yang et al., 2018).

2.2 Measured data

The studied catchment exhibits a temperate and humid climate with pronounced seasonality. The humid climate in wet regions is quantified by an aridity index of 1.0. The ET is the main driver of the hydrologic seasonality, as the precipitation is uniformly distributed across the year. According to the meteorological data records from 1997 to 2010, the mean annual precipitation (P) and discharge (Q) are 610 and 160 mm, respectively (Fig. 1c). Based on the 14-year water balance ($P=\text{ET}+Q$), the actual mean annual evapotranspiration (ET) is 450 mm. The in-stream nitrate concentration (C_Q) was measured at fortnightly to monthly intervals (Dupas et al., 2017), covering the period 2001–2010. The N surplus, which is the annual amount of nitrogen remaining in the soil after the consumption by plant uptake from the external N input, was estimated as 48.8 kg N ha⁻¹ yr⁻¹ during 1997–2010 for this catchment (Bach and Frede, 1998; Bach et al., 2011). The C_Q and N surplus are used to calibrate the N transport model (Wang et al., 2023). Adequate data from numerous investigations and previous research support the exploration of N dynamics in the agricultural catchment.

3 Methodology

The hydrological model of Schäfertal catchment was established using HydroGeoSphere (Therrien et al., 2010) in the previous study (Yang et al., 2018). The framework of N dynamics (Yang et al., 2021), modified from the ELEMeNT modeling approach (Exploration of Long-tErM Nutrient Trajectories, Van Meter et al., 2017), was applied to track the fate of N in soil and groundwater in the present study. Based on the hydrological and transport model, the effect of topographic slope on the export of nitrate (Yang et al., 2022) and the spatial-temporal variation of nitrogen retention (Wang et al., 2023) were explored. In the present study, the two models were employed to investigate the effect of rainfall variability on N dynamics, for which rainfall time series with climatic characteristics, substituting for the precipitation data during the simulation period, were generated by a stochastic rainfall generator (Robinson and Sivapalan, 1997; Wilusz et al., 2017). These models are described below.

3.1 Flow modeling

HydroGeoSphere is a 3D control volume finite element simulator that can model fully coupled surface-subsurface hydrological processes by the dual nodes approach. It can not only describe crucial hydrological processes such as dynamical evapotranspiration, snow, sublimation, snowmelt, and freeze and melt of pore water, but also simulate 2D overland flow by Manning's equation and the diffusive-wave approximation of the St. Venant equations, 3D variably saturated underground flow by Richards' equation and Darcy's law, flow in porous media, as well as the transport of reactive solutes. In addition, HydroGeoSphere allows the simulation of 1D surface flow in a channel network and water exchange flux between the channel domain and the subsurface domain (Yang et al., 2015). HydroGeoSphere has been frequently used to model catchment hydrological processes and solute transport in many previous studies (e.g., Therrien et al., 2010; Yang et al., 2018). Please refer to Therrien et al. (2010) for the governing equations and technical details.

The hydrological model of Schäfertal catchment is briefly recapitulated in the following. More details are provided in Yang et al. (2018). In the model, the subsurface domain between the surface and bedrock was depicted using 3D prisms with side length ranging from 30 to 50 m. The surface domain was filled by the uppermost 2D triangles of the generated 3D prismatic mesh, while the channel was delineated by 1D segments with specified widths and depths. According to different aquifer materials and the difference in their permeability in the vertical direction, the subsurface domain was separated into ten property zones with zonal hydraulic conductivity and porosity values. The surface domain and channel were parameterized with a Manning roughness coefficient representing the land use. Spatially uniform and temporally variable precipitation (Fig. 1c) was applied to the entire surface domain. ET is computed as a combination of plant transpiration from the root zone and evaporation down to the evaporation depth (Therrien et al., 2010). A critical depth boundary condition was assigned to the outlet in the channel domain to simulate the discharge (Q) of the catchment. In order to eliminate the influence of initial conditions, a preliminary model run was performed from 1997 to 2007. The simulated results at the end of the period were taken as initial conditions for the actual simulations.

Key parameters that significantly influence the hydrological processes were selected for the calibration (Anis and Rode, 2015; Graeff et al., 2009) using the software package PEST (Doherty and Hunt, 2010). PEST uses the Marquardt method to minimize a target function by varying the values of a given set of parameters until the optimization criterion is reached. The groundwater levels measured in groundwater wells and the discharge measured at the gauging station were used as target variables. Considering high CPU time (∼1 d to run the model for the period from 1997 to 2007), the calibration period spans 1 year, from October 2002 to October 2003. Subsequently, the calibrated model was verified by reproducing time-variable groundwater levels for the wells over the entire simulation period (Yang et al., 2018).

3.2 Nitrogen transport in soil and groundwater

The transformation and transport of nitrogen in the subsurface are tracked by the framework of N dynamics modified from ELEMeNT. ELEMeNT is a comprehensive model maintaining a landscape memory, that is, considering the effect of not only current conditions but also past land use and nutrient dynamics on current fluxes (Van Meter et al., 2017). External nitrogen input goes through transformation and transport in the soil, and subsequently filters into groundwater, and gets exported to the surface water body. Throughout the nitrogen cycle, various forms of nitrogen undergo complex biogeochemical processes.

The framework includes a N source zone forming in shallow soil and a groundwater zone (Fig. 2). There are two assumptions in the N source zone: (1) the total N load in the N source zone is comprised of an organic N (SON) pool and an inorganic N (SIN) pool; (2) the external N input contributes only to the SON. The external N input represents atmospheric deposition, biological fixation, animal manure from the pasture area, and N fertilizer from the farmland. The SON is further distinguished as active SON with faster reaction kinetics and protected SON with slower reaction kinetics. Both SON(a) and SON(p) are transformed into SIN by mineralization. SIN is further consumed by plant uptake and denitrification, and leaches from soil to groundwater as dissolved inorganic N (DIN), representing mainly nitrate in the studied catchment (Yang et al., 2018; Nguyen et al., 2021). DIN can further undergo denitrification until being exported to the stream. The framework is able to capture the main processes of nitrogen transformation and transport in soil and groundwater (Yang et al., 2018).

Figure 2 The framework simulating the transformation and transport of nitrogen in soil N source zone and groundwater zone, modified based on (Yang et al., 2021).

The governing equations to calculate N fluxes follow the ones of the framework in Yang et al. (2021). The land-use dependent protection coefficient h (Van Meter et al., 2017) determines the amount of external N input that contributes to the protected SON (SON(p)), and the residual contributes to the active SON (SON(a)). Mineralization and denitrification are described as first-order processes. Based on the results of Wang et al. (2023), the effect of wetness is considered in mineralization, using:

$$\begin{array}{}\text{(1)}& {\displaystyle }{\text{MINE}}_{\mathrm{a}}=k_{\mathrm{a}}\cdot f\left(\text{temp}\right)\cdot f\left(\text{wetness}\right)\cdot {\text{SON}}_{\mathrm{a}}\text{(2)}& {\displaystyle }{\text{MINE}}_{\mathrm{p}}=k_{\mathrm{p}}\cdot f\left(\text{temp}\right)\cdot f\left(\text{wetness}\right)\cdot {\text{SON}}_{\mathrm{p}}\text{(3)}& {\displaystyle }{\text{DENI}}_{\mathrm{s}}={\mathit{\lambda }}_{\mathrm{s}}\cdot \text{SIN}\text{(4)}& {\displaystyle }{\text{DENI}}_{\mathrm{g}}=\mathit{\lambda }\cdot \text{DIN}\end{array}$$

where MINE_a, MINE_p (kg ha⁻¹ d⁻¹) are the mineralization rates for active SON and protected SON, respectively. DENI_s and DENI_g (kg ha⁻¹ d⁻¹) are the denitrification rates for SIN and DIN, respectively. k_a, k_p, λ_s, and λ (d⁻¹) are coefficients for the four first-order processes. f(temp) and f(wetness) are factors representing the constraints by soil temperature and wetness (Lindström et al., 2010; Wang et al., 2023), respectively. Plant uptake rate UPT follows the equation used in the HYPE model (Lindström et al., 2010):

$$\begin{array}{}\text{(5)}& {\displaystyle }\text{UPT}=min({\text{UPT}}_{\mathrm{P}},\mathrm{0.8}\cdot \text{SIN})\text{(6)}& {\displaystyle }\begin{array}{rl}{\text{UPT}}_{\mathrm{P}}& =p\mathrm{1}/p\mathrm{3}\cdot \left({\displaystyle \frac{p\mathrm{1}-p\mathrm{2}}{p\mathrm{2}}}\right)\\ & \cdot e^{-(\text{DNO}-p\mathrm{4})/p\mathrm{3}}/\left(\mathrm{1}+\left({\displaystyle \frac{p\mathrm{1}-p\mathrm{2}}{p\mathrm{2}}}\right)\right.\\ & {\left.\cdot e^{-\left(\text{DNO}-p\mathrm{4}\right)/p\mathrm{3}}\right)}^{\mathrm{2}}\end{array}\end{array}$$

where UPT and UPT_P (kg ha⁻¹ d⁻¹) are the actual and potential uptake rates, respectively. The logistic plant growth function is considered in the equation of potential uptake rates. DNO is the day number within a year. p1, p2, p3 are three parameters depending on the type of plants, whose units are (kg ha⁻¹), (kg ha⁻¹), and (day), respectively. p4 is the day number of the sowing date. The leaching rate adopts a first-order process considering soil saturation and groundwater velocity, using:

$$\begin{array}{}\text{(7)}& {\displaystyle }\text{LEA}=k_{\mathrm{l}}\cdot f_{\mathrm{w}}\cdot f_q\cdot \text{SIN}\text{(8)}& {\displaystyle }{f}_{\mathrm{w}}={\displaystyle \frac{S-S_{\mathrm{r}}}{\mathrm{1}-S_{\mathrm{r}}}}\text{(9)}& {\displaystyle }{f}_q=\text{MIN}\left({\displaystyle \frac{q}{q_{\mathrm{ref}}}},\mathrm{1}\right)\end{array}$$

where LEA is the leaching flux of SIN from the N source zone to the groundwater, k_l is a leaching coefficient (d⁻¹), f_w and f_q are two factors representing the constraints of soil saturation and groundwater velocity to the leaching process, respectively. S is the soil saturation, and S_r is the residual saturation. q (m d⁻¹) is the groundwater Darcy flow rate, q_ref (m d⁻¹) is a reference Darcy flow rate. This formulation of LEA is modified from the one used in Van Meter et al. (2017).

In these hydrogeochemical processes, a portion of N is retained in the catchment as the biogeochemical legacy in soil or the hydrological legacy in groundwater or permanently leaves catchments by denitrification, which does not degrade the water quality of the catchment during a certain period. Therefore, the N retention is used to quantify a catchment's capacity to prevent nitrogen from entering surface water bodies during a certain period (Wang et al., 2023), which is the fraction of the N retained in the catchment and consumed via denitrification to the total external N input, calculated using (Ehrhardt et al., 2021):

$$\begin{array}{}\text{(10)}& \mathrm{Retention}=\mathrm{1}-{\displaystyle \frac{N_{\mathrm{out}}}{N_{\mathrm{in}}}}=\mathrm{1}-{\displaystyle \frac{{\int }_{{\mathrm{t}}_{\mathrm{1}}}^{t_{\mathrm{2}}}{N}_{\mathrm{outlet}}}{{\int }_{t_{\mathrm{1}}}^{t_{\mathrm{2}}}{N}_{\mathrm{input}}}}\end{array}$$

where the N_outlet is the N mass leaving the catchment through the outlet during the time period (t₁−t₂).

Due to the lack of spatiotemporal variation information of the external N input, its value was fixed at 180 kg ha⁻¹ yr⁻¹ according to Nguyen et al. (2021), where the N balance was simulated in the upper Selke catchment covering the Schäfertal catchment. The protection coefficient h was selected as 0.3 according to the values reported in Van Meter et al. (2017). p4 as the sowing date of plant growth activities was set to 63 d (early March) (Yang et al., 2018). The DIN was transported in the coupled surface water-groundwater system, with longitudinal and transverse dispersion coefficients of 8 and 0.8 m, respectively. Other parameters relative to N dynamics were calibrated by PEST in empirical ranges. The measured C_Q at the gauge station and the N surplus of 48.8 kg N ha⁻¹ yr⁻¹ were used as the target variables. The calibration period spans from March 2001 to August 2003, during which successive C_Q data were obtained. The CPU time was ∼2 h for a single iteration.

The N source zone serves as a boundary condition at the top of the aquifer for simulating DIN transport in the groundwater. The bedrock is treated as impermeable for water and nitrate. The catchment outlet is the only boundary allowing the exit of nitrate from the catchment. In order to minimize the influence of the initial conditions, a preliminary transport simulation (together with the flow simulation) was performed with zero loads in the SON and SIN pools and zero DIN concentration in the catchment, such that a quasi-steady state for the SON and SIN pools can be reached at the end of the preliminary simulation. The resulting data (N loads and concentrations) were used as initial conditions for the actual simulations.

3.3 Stochastic rainfall generator

In order to investigate the effect of rainfall variability on N dynamics, rainfall time series with different climatic characteristics were generated by a stochastic rainfall generator. The stochastic rainfall generator is a stochastic model of rainfall time series originally created to investigate the timescales of flood frequency response to rainfall, incorporating any combinations of storm, within-storm as well as between-storm, and seasonal variabilities of rainfall intensity (Robinson and Sivapalan, 1997). It can output rainfall series representing different rainfall patterns under climate change. The equations can be found in Sect. S1 of the Supplement.

We adopted the source code in Python v2.7 of the rainfall generator from Wilusz et al. (2017). The stochastic model derives monthly average storm durations and inter-storm periods from seasonal parameters: averages (γ_s, δ_s), amplitudes (α_γ, α_δ), and phases (τ_γ, τ_δ). $a_{\mathrm{1}}^{\mathrm{1}}\sim a_{\mathrm{1}}^{\mathrm{4}}$, $b_{\mathrm{1}}^{\mathrm{1}}\sim b_{\mathrm{1}}^{\mathrm{4}}$, a₂ and b₂ are the parameters characterizing the dependence of average rainfall intensity on storm duration for the first season (January–March), the second season (April–June), the third season (July–September), and the fourth season (October–December). In addition, there is an isolated parameter P_drizzle in the source code, the probability of drizzle events. In the model, storms with precipitation smaller than 4 mm are drizzle events, while storms with precipitation greater than 4 mm can be identified as synoptic frontal events.

Rainfall time series representing different climate conditions were generated following these steps below:

1. To represent contrasting hydroclimatic conditions, three years were selected from the 1997–2022 record (mean 607.9 mm; range 408.2–916.3 mm): the wet year (2007, P=916.3 mm), the normal year (2008, P=588.7 mm), and the dry year (2018, P=444.1 mm);

2. For each of the representative years, a set of parameters was determined so that the generated rainfall time series can fit the actual rainfall data of this year best. This inverse process was conducted using PEST;

3. The three sets of best-fit parameters (Table 1) were used again to generate 100 stochastic realizations (rainfall time series), for the wet, normal and dry year, respectively. These realizations may deviate from the actual rainfall data in terms of daily rainfall values, still being statistically identical with the rainfall pattern of the representative years.

Table 1 Three parameter sets of the rainfall generator for the wet, normal and dry year, respectively.

3.4 Simulation scenarios

3.4.1 Inter-annual rainfall variability

Wet Year (WY), Normal Year (NY), and Dry Year (DY) (Table 1) were set as the simulation scenarios. In addition, a fourth scenario consistent with the Dry Year but coupled with reduced plant uptake was considered as Extreme Dry Year (EDY). This EDY scenario was used to account for the extreme drought occurring in 2018 that caused vegetation to die back and impacted N dynamics to some extent (Winter et al., 2023). In the EDY scenario, the plant uptake was assumed to decrease down to 36 % of the original value used in the DY scenario, according to the classification and the occurrence period of the drought (Liu et al., 2010). To isolate the causal effects of rainfall variability on N dynamics, the time-invariant fertilization rate was used across all scenarios. For each of the scenarios, the water flow and N transport model were conducted 100 times by substituting the rainfall data during the simulation period with the generated 100 stochastic rainfall time series (Fig. S1, Supplement). Average annual N loads and fluxes in each scenario were calculated to ensure that they were not controlled by a single realization and were statistically meaningful. Additionally, we simulated the fluctuations of C_Q in each scenario with time over three years. Finally, the average annual N loads, fluxes, and C_Q as well as the variation of C_Q were cross-compared and analyzed to elucidate the effect of inter-annual rainfall variability on N dynamics.

3.4.2 Intra-annual rainfall variability

In order to explore the effect of intra-annual rainfall variability on N dynamics, simple linear regressions were performed between stochastic rainfall-generator parameters and each of the simulated annual N loads, fluxes, and C_Q, in which the parameters of NY served as a reference. The seasonal averages of storm duration and inter-storm period (γ_s and δ_s), the amplitudes of the seasonal variations in storm duration and inter-storm period (α_γ and α_δ), the average rainfall intensity of four seasons (E₁–E₄) and squared coefficient of variation of the average rainfall intensity (CV²), as well as the probability of drizzle events (P_drizzle), were selected as the experimental parameters. E₁–E₄ and CV² can be approximately calculated by Eqs. (S5) and (S6), using relative parameters and seasonally averaged storm duration (γ_s). In total, 10 parameters (Table 1) were subjected to the assessment of their effect on N dynamics.

For each of the experimental parameters, the rainfall generator generated 100 rainfall time series randomly, whose annual precipitation amounts are within the historical range, with this experimental parameter varying randomly within the adjustable range but other parameters being fixed to the best-fit values for the NY scenario. The ranges of the experimental parameters used in the linear regression analyses are listed in Table 1. Finally, the response of the annual N loads, fluxes and C_Q to different rainfall parameters can be analyzed.

4 Results

The calibrated N transport model showed good performance in fitting the in-stream nitrate concentration (Fig. 4a), with a Nash-Sutcliffe efficiency (NSE) of 0.79. The modeled N surplus of 51.87 kg ha⁻¹ yr⁻¹ is comparable to the measured value of 48.8 kg ha⁻¹ yr⁻¹. The calibrated best-fit values for the transport parameters are listed in Table S1 (see Supplement).

Figure 3 The simulated 14-year N mass balance in the entire catchment.

Figure 3 illustrates the 14-year N mass balance simulated by the calibrated N transport model in the entire catchment. In the soil source zone, the total N consisted of SON (552 kg ha⁻¹, including SON(a) of 90 kg ha⁻¹ and SON(p) of 462 kg ha⁻¹) and SIN (48 kg ha⁻¹). The load of SON accounts for 92 % of the total N, which corresponds to the research result that the organic N fraction is greater than 90 % (Stevenson, 1995). As for N transformation, the mineralization rate of 173 kg ha⁻¹ yr⁻¹ is within the range (14–187 kg ha⁻¹ yr⁻¹) reported by Heumann et al. (2011) for study sites located in Germany. 71 % of the SIN was absorbed by plants at a rate of 123 kg ha⁻¹ yr⁻¹, which is very close to the value (around 120 kg ha⁻¹ yr⁻¹) suggested in Nguyen et al. (2021) for the same area. The denitrification flux of 38 kg ha⁻¹ yr⁻¹ (4 kg ha⁻¹ yr⁻¹ in the soil source zone, 34 kg ha⁻¹ yr⁻¹ in the groundwater) is within the range (8–51 kg ha⁻¹ yr⁻¹) investigated for 336 agricultural areas around the world by Hofstra and Bouwman (2005). The SIN entered the groundwater zone at a rate of 44 kg ha⁻¹ yr⁻¹, which is within the range (15–60 kg N ha⁻¹ yr⁻¹) reported in Nguyen et al. (2021). Approximately 23 % of the leachate of SIN was exported to the stream, where the rest of the leachate (77 %) was removed by denitrification. The export and denitrification fluxes are comparable to the values of 39 % and 61 % reported in Nguyen et al. (2021). Based on measured groundwater chemistry data from nearby catchments and multiple estimates of denitrification potential that have been proved and reported (Hannappel et al., 2018), we conclude that the simulated transformation and transport of N are acceptable.

Figure 4 The fluctuations of simulated in-stream nitrate concentration (C_Q) under (a) the historical rainfall (1997–2007) and the scenarios of (b₁) WY, (b₂) NY, (b₃) DY, and (b₄) EDY (2008–2010). The grey areas are formed by the realizations of each scenario.

Figure 5 The simulated annual in-stream nitrate concentration (C_Q) for scenarios WY, NY, DY, and EDY. The whisker represents the concentration ranging from the 5th to 95th percentiles, with the triangle indicating the median and the lines marking the maximum, average, and minimum of the data set from top to bottom. The same notations are used below.

4.1 Effect of Inter-annual rainfall variability

The variations of C_Q in these scenarios are shown in Fig. 4b₁–b₄. The range of C_Q generally decreases when the rainfall pattern shifts from WY, via NY, to DY, but the range of EDY is wider than that of WY. Figure 5 shows the range of the simulated annual in-stream nitrate concentration for the four scenarios. WY produced higher C_Q, with the median concentration reaching nearly 5 mg L⁻¹. The concentrations generally decrease when the rainfall pattern transforms from WY, via NY, to DY. However, when vegetation dieback occurred, most concentration values in the EDY scenario were far higher than those of WY. As shown in the results, the median concentration of EDY reached 8 mg L⁻¹. The soil organic nitrogen (SON(a) and SON(p)) loads increase when the rainfall pattern shifts from WY, via NY, to DY/EDY (Fig. 6a and b). The SON loads in DY and EDY were identical because SON transformation is independent of vegetation state. There were no remarkable differences in soil inorganic nitrogen (SIN) load among WY, NY, and DY (Fig. 6c). However, the SIN load in EDY was the highest among the four scenarios.

Figure 6 The simulated loads of (a) active soil organic nitrogen (SON(a)), (b) protected soil organic nitrogen (SON(p)), and (c) soil inorganic nitrogen (SIN) for scenarios WY, NY, DY, and EDY.

The highest and the lowest average mineralization levels were produced in WY and DY, respectively (Fig. 7a). DY and EDY have identical mineralization fluxes. Note that the fluxes of plant uptake, denitrification (in soil and groundwater), and leaching are negative values representing sink terms for the SIN pool (Fig. 7b, c and d). These fluxes decrease when the rainfall pattern transforms from WY, via NY, to DY, due to decreasing soil moisture. Because of the accumulated SIN load (Fig. 6c) that resulted from the lowest level of plant uptake, denitrification and leaching fluxes in EDY with low annual precipitation are still even larger than those of WY. The N export flux to the stream generally follows the same trend as leaching, with relatively higher fluxes occurring in WY and EDY and lower fluxes in NY and DY (Fig. 7e). Based on these results, it can be preliminarily concluded that annual precipitation and vegetation state critically influence N transformation and transport.

Figure 7 The simulated fluxes of (a) mineralization, (b) plant uptake, (c) denitrification (in soil and groundwater), (d) leaching, and (e) N export for scenarios WY, NY, DY, and EDY.

4.2 Effect of Intra-annual rainfall variability

The determination coefficients (R²) of the linear regressions between the rainfall-generator parameters and the C_Q, N loads, and fluxes are illustrated in Fig. 8a. Larger R² values can be observed for the amplitudes of the seasonal variations in the storm duration and inter-storm period (α_γ and α_δ), the average rainfall intensity of the third season (E₃) as well as the probability of drizzle events (P_drizzle). These parameters are the four most important parameters (Fig. 8b) that influence N dynamics, and are subject to further discussion.

Figure 8 (a) The determination coefficients (R²) of the linear regressions between the parameters of the rainfall generator and in-stream nitrate concentration (C_Q), N loads, and fluxes. (b) The importance of the parameters influencing the N dynamics.

4.2.1 Storm duration

The parameter α_γ determines how the monthly average storm duration is distributed over the course of a year (Fig. 9a). Generally, larger α_γ values represent that storms with longer duration are more likely to occur in mid-year when ET is high and the inter-storm duration is longer (yellow zone, Fig. 9e). This pattern not only influences the drying and rewetting cycles in summer, but also supplies more water during this period, thereby leading to an overall wet year (e.g., year 2007, α_γ=1.764). A decrease in α_γ causes storms with longer durations to shift towards the beginning/end of the year, when ET is low and the inter-storm duration is shorter (blue zone, Fig. 9e), so that the year becomes drier (e.g., year 2018, α_γ=0.55).

Results suggest that the storm duration distribution during a year significantly affects N transformation and transport. Lower SON loads (Fig. 9a; $k=-\mathrm{28.04}$, p<0.05) and higher SIN loads (Fig. 9b; k=13.31, p<0.05) occurred when longer-duration storms concentrate in mid-year. This is because wetter catchment conditions promote the mineralization of SON into SIN, especially during warm periods when mineralization is primarily constrained by soil moisture rather than temperature. The wetter conditions associated with larger α_γ values increased the leaching flux of SIN from soil to groundwater (Fig. 9c; $k=-\mathrm{4.77}$, p≥0.05), which elevated the average C_Q (Fig. 9d; k=0.7, p<0.05). When shorter-duration storms dominate mid-year, N transformation and transport were inhibited.

Figure 9 The effect of varying monthly average storm duration distributions over the course of a year on (a) SON load, (b) SIN load, (c) leaching flux (LEA), and (d) in-stream nitrate concentration (C_Q). (e) The monthly average inter-storm period distribution over the course of a year (α_δ is kept constant when α_γ varies). As α_γ values decrease, longer-duration storms shift gradually from the period with high ET and longer inter-storm periods (yellow zone) to the period with low ET and shorter inter-storm periods (blue zone).

4.2.2 Inter-storm period

α_δ regulates the variation in monthly average inter-storm periods throughout the course of a year (Fig. 10a). Larger α_δ values cause longer inter-storm periods to cluster primarily in mid-year, coinciding with longer storm durations (blue zone, Fig. 10e), while shorter inter-storm periods concentrate at the beginning/end of the year, associated with shorter storm durations (yellow zone, Fig. 10e). This forms a “wet–dry–wet” climate pattern that promoted SON accumulation (Fig. 10a; k=3.39, p<0.05) but reduced SIN load (Fig. 10b; $k=-\mathrm{1.34}$, p<0.05), leaching flux (Fig. 10c; k=1.73, p<0.05), and C_Q (Fig. 10d; $k=-\mathrm{0.18}$, p<0.05). Conversely, decreasing α_δ redistributes longer inter-storm periods to the year's boundaries and shorter inter-storm periods to mid-year, thereby forming a “dry–wet–dry” climate pattern over a year. It resulted in a lower SON load and a higher SIN load, as mid-year warmth and moisture promoted the transformation and transport of SON accumulated in the previous season. Consequently, leaching flux and C_Q significantly increased.

Figure 10 The effect of varying monthly average inter-storm period distributions over the course of a year on (a) SON load, (b) SIN load, (c) leaching flux (LEA), and (d) in-stream nitrate concentration (C_Q). (e) The monthly average storm duration distribution over the course of a year (α_γ is kept constant when α_δ varies). As α_δ values decrease, longer-period inter-storms shift gradually from the period with longer storm durations (blue zone) to the period with shorter storm durations (yellow zone).

4.2.3 Average Storm Intensity

Figure 11 shows the responses of N loads, fluxes, and in-stream nitrate concentration (C_Q) to seasonal average rainfall intensity (E₁–E₄). Fundamentally, increased rainfall intensity creates more humid conditions that intensify biogeochemical processes (e.g., mineralization, plant uptake). Based on determination coefficients (R²), linear regression slopes (k), and the significance, E₂ and E₃ generally exert stronger effects on N loads and fluxes than E₁ and E₄. Because mineralization was not constrained by temperature during the warm seasons (Fig. 11b₄ and c₄), elevated $E_{\mathrm{2}}/E_{\mathrm{3}}$ enhanced the transformation of SON(a) and SON(p), leading to significant decreases in their loads (Fig. 11b₁, c₁, b₂ and c₂). Although SIN load theoretically should increase with rainfall intensity, as the only SON sink, it decreased in the second season (Fig. 11b₃) due to substantially increased plant uptake (Fig. 11b₅). Warmer temperatures and enhanced moisture during the growing season promoted vigorous nutrient absorption by vegetation (Fig. 11b₅ and c₅). In addition, soil denitrification increased markedly during the third season as average rainfall intensity increased (Fig. 11c₆), which was due to favorable microbial conditions. In the fourth season, soil denitrification decreased slightly with increasing average rainfall intensity (Fig. 11d₆), which was attributed to low temperatures and decreased SIN loads (Fig. 11d₃). In warm periods, enhanced average rainfall intensity increased actual evapotranspiration (ET) (Fig. S2b₁ and c₁). The linear relationships between seasonal-average rainfall intensity and recharge, Q, N export, and C_Q are weak, as indicated by low R² values and small slopes (k). Although several regressions reach statistical significance (p<0.05), their effects are minor (Figs. 11b₇, S2a₂ and d₄).

Figure 11 The responses of N loads, mineralization, plant uptake, and denitrification (in soil) fluxes, as well as in-stream nitrate concentration (C_Q) to the average rainfall intensity of the four seasons (E₁–E₄). The sign and magnitude of the slopes (k) in these linear relationships denote the direction and intensity of the response of N dynamics to the variations in average rainfall intensity, respectively. Asterisks indicate the significance of the regression slopes (p<0.05); ns denotes non-significant relationships (p≥0.05).

4.2.4 Probability of drizzle event

Figure 12 illustrates the responses of N loads, fluxes, discharge (Q), and in-stream nitrate concentration (C_Q) to the probability of drizzle events (P_drizzle). P_drizzle characterizes the pattern of precipitation intensity over the course of a year. A low P_drizzle value indicates a predominance of moderately intense rainfall evenly distributed throughout the year, resulting in consistently humid conditions (e.g., year 2007, P_drizzle=0.137). Conversely, a high P_drizzle value represents a situation with high-frequency drizzle and low-frequency extreme precipitation, forming drought-flood extremes. With the increase of P_drizzle (climate transition from consistently humid conditions to extreme drought-flood patterns), SON load decreased due to the enhanced transformation (Fig. 12a₁ and b₁). Elevated plant uptake (Fig. 12a₂) resulted in reduced SIN load (Fig. 12c₁). Denitrification and leaching fluxes were reduced (Fig. 12b₂ and c₂), while C_Q declined with increasing P_drizzle (Fig. 12c₃). The linear relationships between P_drizzle and N export and between P_drizzle and C_Q are not significant (Fig. 12a₃ and b₃; p≥0.05).

Figure 12 The responses of N loads, fluxes, discharge (Q) and in-stream nitrate concentration (C_Q) to the probability of drizzle events (P_drizzle). The sign and magnitude of the slopes (k) in these linear relationships denote the direction and intensity of the response of N dynamics to the probability of drizzle events (P_drizzle), respectively. Asterisks indicate the significance of the regression slopes (p<0.05); ns denotes non-significant relationships (p≥0.05).

5 Discussion

5.1 N transformation and transport upgrade in wet years

The comparison of C_Q, N loads, and fluxes across four scenarios (WY, NY, DY, and EDY) reveals the effect of inter-annual rainfall variability on N transformation and transport, as well as water quality. Mineralization, the crucial process of the transformation from SON to SIN, exhibits strong soil moisture dependence. Consequently, the highest average mineralization rates in the WY scenario promoted N transformation through plant uptake, denitrification, and leaching. In contrast, low mineralization in the DY scenario led to SON accumulation, thereby restraining overall N dynamics. Additionally, the transport of SIN relies on groundwater discharge and subsurface flow paths. Increased discharge from precipitation changes flushes more nitrogen into surface water bodies (Mitchell et al., 1996; Creed and Band, 1998) via fast shallow flow paths with short transit times (I-I cross section, Fig. 1a) (Yang et al., 2018). Consequently, N export and C_Q variations align with wetness changes among WY, NY, and DY. These results demonstrate that high annual precipitation enhances N transformation and transport, while low annual precipitation promotes catchment N retention. Our simulation results echo previous findings that N retention is influenced by changes in precipitation (Dumont et al., 2005; Howarth et al., 2006).

5.2 Stream water quality deteriorates during extreme droughts

SIN significantly accumulated in soil when the vegetation suffered from drought stress (RIVM – National Institue for Public Health and the Environment, 2021; Winter et al., 2023). This is because in general plant uptake is the major sink for SIN (Yang et al., 2021; Zhou et al., 2022; Wang et al., 2023). Due to the accumulation of SIN, denitrification and leaching in the EDY scenario were higher than those in the DY scenario. Thus, N export to the stream in the EDY scenario was higher than that in the DY scenario, even though their stream discharges were the same (Fig. 13). This resulted in higher C_Q in the EDY scenario than in the WY scenario.

Figure 13 The simulated discharge (Q) for scenarios WY, NY, DY, and EDY.

During the 2018–2019 drought in Central Europe, the observed peaks of in-stream nitrate concentration at a mesoscale catchment in central Germany were significantly higher than the previous concentration peaks (Fig. 14). Based on a data-driven analysis, Winter et al. (2023) indicated that nitrate loads in the 2018–2019 drought were up to 73 % higher than the long-term average loads. They demonstrated that such increases were attributed to decreased plant uptake and subsequent flushing of accumulated nitrogen during the rewetting period. Our results confirm their findings that reduced plant uptake in extreme droughts causes C_Q to elevate significantly (Fig. 4b₄) and the N retention capacity to decline. These results suggest that vegetation state plays a vital role in the increased risk of N pollution during extreme droughts.

Figure 14 Measured daily averaged nitrate concentrations (C, red line) and discharge (Q, black line) in the nearby Upper Selke catchment in central Germany (Winter et al., 2023). The yellow area marks the 2018–2019 drought occurring across Central Europe.

5.3 The distribution of storm durations controls the N loads

Our study indicates that the wet/dry conditions in mid-year are key to N transformation and transport throughout the whole year. The combination of longer storm durations and extended inter-storm periods in mid-year represents an aspect of extreme climate events where intense rainfall events alternate with droughts (Zhang et al., 2024). Without temperature limitations, alternating longer storm durations and longer inter-storm periods increase the extensive SON transformation, which is consistent with the results of laboratory studies that mineralization increases after successive drying and rewetting cycles in soils (Cabrera, 1993). It is concluded that drying-rewetting events induce prominent changes in N dynamics (Fierer and Schimel, 2002). However, the drier conditions induced by shortened storm durations and longer inter-storm periods in mid-year cause the accumulation of SON, which is consistent with the finding by Zhou et al. (2022) that the drought in the warm season contributes to catchment nitrogen retention. Thus, storm duration distributions within a year can significantly affect N loads in a catchment.

5.4 The distribution of inter-storm periods changes stream water quality

Different dry–wet patterns over the course of a year have an impact on N dynamics, particularly on in-stream nitrate concentration. In the “dry–wet–dry” pattern, accumulated SON from the dry season is transformed and transported during the subsequent wetting period, which leads to higher C_Q. During the 2018–2019 drought, accumulated N loads were flushed during the rewetting period and thereby caused elevated nitrate loads (Fig. 14). In the “wet–dry–wet” pattern, humid conditions with low temperatures at the year's boundaries and warm droughts in mid-year both cause SON loads to accumulate throughout the entire year and C_Q to decrease. Zhou et al. (2022) attributed reduced C_Q to limited terrestrial export loads. In summary, inter-storm period distributions cause noticeable variations in water quality in terms of nitrate concentration.

5.5 Average rainfall intensity has only a small effect on stream water quality

Wet conditions caused by high mean rainfall intensities enhance the transformation of SON to SIN (Knapp and Smith, 2001), when mineralization is not constrained by low temperatures (Eqs. 1 and 2). The potential plant uptake is subject to vegetation state (Eq. 6), while the actual plant uptake is further limited by SIN load and soil moisture content (Eq. 5) (Yang et al., 2022; Cramer et al., 2009). Thus, high mean rainfall intensities promote plant absorption of SIN during active growth periods. Since plant uptake is the major sink of SIN, the responses of SIN load to mean rainfall intensity variations are jointly determined by mineralization and plant uptake, resulting in complex SIN load dynamics. SIN load controls denitrification in soil and leaching. However, denitrification is dependent on microbial conditions as well. Leaching is additionally influenced by soil saturation and groundwater velocity (Eqs. 3, 4, and 7). As a result, no significant linear relationship exists between mean rainfall intensity and leaching. It is noteworthy that high-intensity precipitation events with short durations and substantial surface runoff rarely reach the water table, thereby exerting minimal effects on recharge, discharge, and N export (a component of leaching flux) (Fig. S2, Supplement). Therefore, C_Q shows weak and non-significant responses to extreme precipitation.

5.6 Probability of drizzle event alters N dynamics

Extreme dry–wet patterns enhance the transformation of SON more effectively than continuously humid conditions, which echoes laboratory studies conducted by Cabrera (1993) that successive drying and rewetting cycles in soils can intensify mineralization. As the major sink, elevated plant uptake due to increased SIN load in turn reduced SIN load, which resulted in weaker denitrification and leaching. This is because denitrification and leaching are controlled by SIN load (Eqs. 3, 4, and 7). Following the same trend as leaching, C_Q slightly decreased. It follows that the probability of drizzle events prominently alters N dynamics.

On the whole, increased variations in precipitation can alter the transformation and transport of nitrogen (Kane et al., 2008). The results of the present study can help in managing the water quality of agricultural catchments with prominent rainfall variability and protecting aquatic ecosystems during extreme rainfall events and droughts. For instance, a fertilization scheme for inter-storm periods can be formulated to avoid flushing during the rewetting periods after droughts. Also, increased fertilization in the growing season with humid conditions is conducive to vegetation absorption of nutrients while not putting pressure on soil and water bodies.

5.7 Limitations and Outlook

The transport model preserves the main pathway for N-NO${}_{\mathrm{3}}^{-}$ leachate by simplifying the complexities of different N pools and transformations via mineralization, leaching, and denitrification within the soil zone. However, the external N input entering directly into inorganic N, the transformation from protected organic N to active organic N, and the loss of organic N via dissolution have not been included in the model, which may lead to miscalculations of nitrogen loads. As a complex biogeochemical process, denitrification is governed by various factors such as temperature, soil moisture content, and SIN load, rather than being a first-order decay process. The current transport model does not account for spatial heterogeneity in denitrification. These simplifications may introduce uncertainties in the simulated results. In addition, due to the lack of groundwater nitrate concentration data, the nitrate transport model was calibrated using only in-stream nitrate concentration data (Wang et al., 2023). More nitrate concentration data in groundwater and rivers will contribute to establishing a more realistic nitrate transport model.

The plant-soil process was not represented in the model. In the extreme dry scenario, plant death was not actually simulated under high temperature and low soil moisture content, but plants were manually assigned a lower plant uptake potential. The plant growth stages were assumed to be constant. This is because the N uptake process was only mathematically described using the empirical formula (Eq. 6) in the simplified N framework, rather than using a full plant-soil process. However, such simplification was proved to be effective in terms of reproducing reasonable N loads, fluxes and in-stream nitrate concentrations. Therefore, we think it is acceptable to identify the key effect of plant die-off (i.e., reduced N uptake potential) on stream water quality during extreme dry scenarios.

In the future, the alteration of soil properties caused by rainfall variability should be studied. In the model, precipitation and the solute enter the subsurface domain with fixed soil porosity, ignoring the alteration of pores. Microorganisms and soil fauna are highly active, and soil organic and inorganic matter undergo continual biochemical reactions in soil pores, which all change the physical and chemical properties of porous media and thereby impact seepage. These biogeochemical processes are impacted by several factors, such as hydrometeorological conditions (Ondrasek et al., 2019). For example, the activities of microorganisms and animals degrade under drying conditions, which is not beneficial for maintaining and improving soil structures. The severe deficit of soil moisture results in the disintegration of soil aggregates, and even compacted soil forming in the upper layer due to salt accumulation, which ultimately reduces porosity. Additionally, cracked surface soil caused by extreme droughts is detrimental to plant growth and root development. Therefore, altered soil properties induced by complex and extreme climate patterns can substantially influence soil biogeochemical activities. The clarification of these mechanisms will be conducive to exploring the more realistic response of solute transport in groundwater to extreme rainfall events and droughts.

6 Conclusions

Within the context of climate change and the increasing occurrence of extreme rainfall events and droughts, the study pioneered the application of a stochastic rainfall generator (Robinson and Sivapalan, 1997) coupled with the flow and transport model to systematically investigate inter-annual and intra-annual rainfall variability effects on N dynamics in a small agricultural catchment located in central Germany. The principal findings are as follows:

1. Higher annual precipitation enhances N transformation and transport, whereas lower annual precipitation is conducive to the N retention capacity. The retention capacity declines severely when vegetation suffers from drought stress. Vegetation plays a vital role in the response of N dynamics to extreme droughts.

2. Wet/dry conditions determined by storm duration distributions within a year can significantly affect N loads in a catchment. Droughts can prompt the accumulation of SON, but a drying-wetting cycle can enhance extensive SON transformation in the warm season.

3. Different dry–wet patterns formed by inter-storm period distributions during a year can cause noticeable variations in in-stream nitrate concentrations that prominently elevate during the rewetting period after a drought.

4. Elevated mean rainfall intensity contributes to N transformation when mineralization is not constrained by low temperatures, and promotes plant absorption of SIN during the growing season. There is a small effect on stream water quality.

5. The probability of drizzle events, characterizing the pattern of precipitation intensity within a year, prominently alters N dynamics. Extreme dry–wet patterns enhance SON transformation and plant uptake but inhibit denitrification and leaching more effectively than continuously humid conditions.

Overall, the study clarifies the effect of rainfall variability on N dynamics in a small agricultural catchment, which provides theoretical support for formulating fertilization strategies and protecting aquatic ecosystems in the context of climate change.