Are maps of nitrate reduction in groundwater altered by climate and land use changes?
- 1Department of Hydrology, Geological Survey of Denmark and Greenland, GEUS, Copenhagen 1350, Denmark
- 2Department of Agroecology, Aarhus University, Tjele 8830, Denmark
- 3Department of Bioscience – Lake Ecology, Aarhus University, Silkeborg 8600, Denmark
Correspondence: Ida Karlsson Seidenfaden (firstname.lastname@example.org)
Nitrate reduction maps have been used routinely in northern Europe for calculating the efficiency of remediation measures and the impact of climate change on nitrate leaching. These maps are, therefore, valuable tools for policy analysis and mitigation targeting. Nitrate reduction maps are normally based on output from complex hydrological models and, once generated, are largely assumed constant in time. However, the distribution, magnitude, and efficiency of nitrate reduction cannot necessarily be considered stationary during changing climate and land use as flow paths, nitrate release timing, and their interaction may shift. This study investigates the potential improvement of using transient nitrate reduction maps, compared to a constant nitrate reduction map that is assumed during land use and climate change, both for nitrate loads and the spatial variation in reduction. For this purpose, a crop and soil model (DAISY) was set up to provide nitrate input to a distributed hydrological model (MIKE SHE) for an agricultural catchment in Funen, Denmark. Nitrate reduction maps based on an observed dataset of land use and climate were generated and compared to nitrate reduction maps generated for all combinations of four potential land use change scenarios and four future climate model projections. Nitrate reduction maps were found to be more sensitive to changes in climate, leading to a reduction map change of up to 10 %, while land use changes effects were minor. The study, however, also showed that the reduction maps are products of a range of complex interactions between water fluxes, nitrate use, and timing. What is also important to note is that the choices made for future scenarios, model setup, and assumptions may affect the resulting span in the reduction capability. To account for this uncertainty, multiple approaches, assumptions, and models could be applied for the same area. However, as these models are very time consuming, this is not always a feasible approach in practice. An uncertainty of the order of 10 % on the reduction map may have major impacts on practical water management. It is, therefore, important to acknowledge if such errors are deemed acceptable in relation to the purpose and context of specific water management situations.
Nitrate loads from agricultural areas are recognized to cause harmful impacts on groundwater and surface water resources, including eutrophication in aquatic ecosystems (Diaz and Rosenberg, 2008). This is also the case in the Baltic Sea drainage basin (Reusch et al., 2018), including Denmark, where nitrate load from agriculture constitutes one of the major water resources management challenges. Nitrate is removed by a set of natural near-surface removal processes, including plant uptake and soil retention. Furthermore, the natural removal of nitrate in the groundwater and the surface water must also be considered when assessing the impacts of nitrate leaching from agricultural areas on aquatic ecosystems. This removal takes place via natural biogeochemical reduction processes often referred to as denitrification. It can be expressed as a percentage removal, and depending on the actual hydrobiogeochemical conditions, the denitrification may mainly occur in groundwater or in surface water systems such as lakes or wetlands (Quick et al., 2019; Huno et al., 2018).
In the groundwater zone, nitrate reduction takes place when nitrate containing water migrates from aerobic to anaerobic conditions and inherent reduced compounds are available (Postma et al., 1991; Hansen et al., 2014a). For Quaternary sediments, these reduced compounds are mainly organic carbon and pyrite and ferrous ion from clay minerals (Ernstsen and Mørup, 1992; Postma et al., 1991). This transition zone between aerobic and anaerobic conditions is denoted as the redox interface. The amount of nitrate reduction occurring in the groundwater will then depend on the flow paths and the depth to the redox interface. In areas with Quaternary sediments characterized by groundwater-dominated flow patterns and a relatively shallow redox interface, the nitrate reduction in groundwater can be the dominant removal process (Hansen et al., 2009). For example, Højberg et al. (2015) estimated that, on average, 63 % of the nitrate leaching in Denmark is removed by nitrate reduction in groundwater.
Heterogeneities in geology and drainage systems are responsible for substantial local spatial variations in nitrate reduction. However, the spatial variation in nitrate reduction in the groundwater system has, so far, only been investigated in a handful of studies (e.g. Wriedt and Rode, 2006; Kunkel et al., 2008; Højberg et al., 2015; Merz et al., 2009; Knoll et al., 2020; Tesoriero et al., 2015). Different approaches have been used in these studies from nitrate groundwater modelling (Højberg et al., 2015; Merz et al., 2009; Wriedt and Rode, 2006), data-driven machine learning (Knoll et al., 2020), or statistical modelling (Tesoriero et al., 2015). An approach for utilizing and illustrating the results and the spatially varying nitrate removal fractions (percentages) is through a nitrate reduction map (Hansen et al., 2014a). A nitrate reduction map is typically produced by using a complex hydrological model, including simulations of root zone nitrate leaching and groundwater and surface flow and transport, and has been applied in several catchments in Denmark and in catchments surrounding the Baltic Sea region (Højberg et al., 2017; Wulff et al., 2014; Hansen et al., 2014b). Hansen et al. (2014a) produced nitrate reduction maps for a 101 km2 catchment in Denmark, showing that the nitrate reduction may vary from 20 % to 70 % between neighbouring agricultural fields located only a couple of hundred metres apart. Similarly, Højberg et al. (2015) and Andersen et al. (2016) estimated very large variations in the nitrate reduction between different regions in Denmark and in the Baltic Sea drainage basin, respectively.
The efficiencies of remediation measures at different locations on nitrate loadings can easily be calculated with a nitrate reduction map, and the measures can be spatially targeted to locations, where the natural removal is relatively small and the mitigation effect is relatively large (Hansen et al., 2017; Refsgaard et al., 2019). Similarly, a nitrate reduction map can be used to transform climate change and other land use change impacts on nitrate leaching from agricultural areas to a catchment response (Olesen et al., 2019). Using nitrate reduction maps based on a single model run is clearly a much faster method than running multiple complex hydrological simulation models for large ensembles of scenarios and is, therefore, a practical tool for policy analysis (Højberg et al., 2015; Andersen et al., 2016). A severe problem in this respect is, however, that the nitrate reduction maps may not be constant in time, as the reduction taking place at a given location depends on the resulting flow pathways (Hansen et al., 2014b). It is, therefore, very relevant to investigate the potential error arising when nitrate reduction maps are assumed to be constant in time. No studies have been reported on that issue, even though the link between climate change, land use change, and nitrate reduction has been established in previous studies (e.g. Mas-Pla and Menció, 2019; Ortmeyer et al., 2021; Olesen et al., 2019; Sjøeng et al., 2009; Fleck et al., 2017). Ortmeyer et al. (2021) used a water balance model combined with a lumped-parameter nitrate mass model for an area in Germany, finding that nitrate concentrations in the groundwater increased towards the end of the century by up to 89 % as a result of changes in temperature, evapotranspiration, and precipitation. Mas-Pla and Menció (2019) found that climate change, in turn, affects groundwater recharge and, thus, the dilution of nitrate in the subsurface in a study in Catalonia. Meanwhile, Paradis et al. (2016) found that new agricultural practices under changing climate conditions led to substantial nitrate increases on an island in eastern Canada.
The objectives of the present study are to assess (i) how nitrate reduction maps showing spatially varying nitrate removal fractions in the groundwater zone are affected by changes in climate and land use and (ii) the errors in nitrate loading made by assuming nitrate reduction maps to be constant. The analyses are performed using a complex hydrological simulation model for a Danish catchment to calculate nitrate reduction maps for the present conditions and for scenarios of climate and land use change. The reduction in this catchment has previously been shown to be dominated by saturated zone reduction processes (Hansen et al., 2009).
The study site is located in the central part of Denmark on the island of Funen. It consists of the 486 km2 upstream part of the Odense River basin, where the Kratholm discharge station marks the outlet (Fig. 1). The catchment is drained by a 200 km river network, with the outlet located at the Odense Fjord to the northeast. Land use in the area is predominantly agricultural (68 %), mainly with pig farms followed by dairy and plant production farms (Fig. 1). Forest constitutes only 5 %, urban areas are 8 %, 1 % is water bodies, and the remaining area is either fallow or grasslands (Nielsen et al., 2000). The soil map and parameters consist of 10 soil types that are created by Børgesen et al. (2013) and Greve et al. (2007), based on national databases. The soil type is dominated by clayey soils (71 %), with smaller areas of sand (see Karlsson et al., 2015, 2016, for more information). As a result, the agricultural area is heavily drained. The geology is mainly a result of previous glaciations like till deposits. Aquifers are generally confined, and the phreatic groundwater tables are shallow. The discharge station at Kratholm (ST45.21) has one of the best nutrient time series in Denmark, starting in the 1980s, with near-daily sampling from 1989 (Windolf et al., 2016). The station, therefore, provides a long and near-complete dataset for nutrient modelling, as well as an extensive water discharge time series (Trolle et al., 2019). In 2005–2009, the average discharge amounts to 4.6 m3 s−1, and the transport in the stream (load) is approximately 14 kg NO3-N ha−1 yr−1, as calculated from measurements of mean concentration and mean water discharge. A decreasing trend in nitrate loads has been observed previously, during 2000–2013 by Windolf et al. (2016), possibly due to the implementation of mitigation measures in the catchment. Also present are three other discharge stations in the catchment (from downstream, ST45.01, ST45.28, and ST45.20).
Measurements of the redox depth are available from 226 boreholes in the area. The redox depths were mainly interpreted based on sediment colour, as described by, for example, Ernstsen and Mørup (1992), and a few were interpreted by measurements of the reduced compounds. They show a shallow redox interface, where 50 % of the measurements have a redox depth less than 4.5 m below terrain, while 90 % of the depths are located in the upper 12.9 m. An old redox depth map with 1 km resolution, based on measurements and geological interpretation, shows depths between 1 to 5 m at many locations and between 5 to 15 m in other locations (Ernstsen et al., 2010). A recent redox depth map of Denmark was created in 2019, where measurements and system variables were used in a machine learning environment to create a detailed redox depth map at 100 m resolution. This newer map also indicates that the redox depth in the study area is predominantly shallow, with 1–10 m depth, and very few sites of 10–15 m depth (Koch et al., 2019a, b).
The DAISY–MIKE SHE modelling system is used to describe nitrate transport in the catchment. DAISY is used to quantify the leaching of nitrate from the root zone while MIKE SHE is used to simulate the transport and degradation of nitrate in the saturated zone. Both models are forced by the same daily inputs of precipitation, temperature, and potential evapotranspiration. Both MIKE SHE and DAISY have been used extensively in Denmark, and MIKE SHE forms the basis of the national nitrate and groundwater model (Højberg et al., 2010, 2013, 2015; Troldborg et al., 2010; Bruun et al., 2003; Hoang et al., 2010). MIKE SHE is a fully coupled integrated groundwater–surface water model, and this integration is important for assessing the feedback between the unsaturated and saturated zone, especially under a changing climate. However, MIKE SHE does not simulate crop development and nitrate leaching from the root zone, and therefore, information from an agrological model, like DAISY, is necessary.
Calibration is done in the following three phases: (1) the MIKE SHE model is automatically calibrated using discharge and hydraulic heads (Sect. 3.1), and the spatial distribution of the redox interface is then assessed for the calibrated model (see the subsection titled “Estimation of the redox interface”); (2) DAISY adopts calibrated soil parameter from MIKE SHE, and additional flow and crop parameters are manually calibrated to the catchment-scale water balance and crop yield (Sect. 3.2); (3) the combined DAISY–MIKE SHE nitrate model is manually calibrated using the nitrate arrival percentage (Sect. 3.3) by adjusting the depth of the redox interface layer. The modelling results are then used to create the nitrate reduction maps (Sect. 3.4) for a range of land use and climate scenarios (Sect. 3.5).
3.1 MIKE SHE setup and calibration (phase 1)
The calibration and setup of the MIKE SHE model has been reported in detail by Karlsson et al. (2016) but will be described here for clarity. MIKE SHE is a fully distributed hydrological model that is built as a modular system with configurable complexity of the different flow compartments (Abbott et al., 1986; Graham and Butts, 2005). To obtain consistency with the flow calculations in DAISY, the present MIKE SHE is also based on a 1D finite difference description of the unsaturated zone, using the Richards' equation, and parameterized by van Genuchten (1980) formulations of the retention curve and the unsaturated hydraulic conductivity curve. A 3D finite difference description of the saturated zone, based on Darcy's equation, was selected, including a linear reservoir formulation of the tile drainage flow. The study is heavily tile drained, and it can be assumed that drainage will always be present when it is needed in the agricultural areas. However, the actual site-specific location of tile drains is unknown, and therefore, the drains are specified across the entire catchment at a depth of −0.5 m. Drain flow is, however, only activated when groundwater levels rise above drain level. Apart from representing tile drainage, the drainage system also represents small ditches and streams that are too small to incorporate into the river system, following the approach of Troldborg et al. (2010). River flow is described using the MIKE 1 module, where a relatively simple routing method (Muskingum) is used. Additionally, physical formulations for evapotranspiration (Kristensen and Jensen, 1975) and 2D overland flow (DHI, 2019) are selected. The MIKE SHE model was calibrated using the built-in autocalibration scheme of AutoCal (Madsen, 2000), which uses the global search function termed the population simplex evolution method. Calibration was carried out against data from four discharge stations and 455 groundwater wells with hydraulic head measurements from the period 2004–2007 and validated in the periods 2000–2003 and 2008–2009. The multi-objective function consists of the water balance error (mean daily error) for the four discharge stations and RMSE for the hydraulic head measurements and three of the discharge stations (one discharge station, ST45.28, is omitted here due to human regulation on the flow system). After a sensitivity analysis on 28 free parameters with 43 tied parameters, a total of five parameters were chosen for calibration. Of these, we used one soil parameter in the unsaturated zone, one drainage parameter, and three saturated-zone parameters.
Estimation of the redox interface
After the calibration of MIKE SHE, the estimation of the redox interface was done using a method developed by Hansen et al. (2014a). This five-step method is used to determine the redox interface for each cell in MIKE SHE. The method is based on the assumption that the present location of the redox interface is a result of the cumulative oxygen percolation through the soil column since the last ice age in the Holocene 11 700 ka.
The redox interface is assumed to have been at ground level at the end of the glaciation and to have migrated downwards by an unknown number of millimetres per yearly recharge. Following the procedure of Hansen et al. (2014a), the first step is, therefore, to find the average yearly recharge by running a model simulation without anthropogenic influences (abstraction and tile drainage). In the second step, the different redox capacities in soils are accounted for; here the capacity of sandy soils are multiplied with a factor of 3 compared to values specified for clayey soil types, thus applying the classification from Børgesen et al. (2013) and Greve et al. (2007). The third step generates the redox interface expressed through Eq. (3) as follows:
where Redoxdepthi is the redox depth (metres) calculated at each grid (i), and fluxi is the groundwater recharge (metres per year; hereafter m yr−1) which is multiplied by the migration constant f (year; Hansen et al., 2014a). The upper part of the unsaturated zone is assumed to have no redox capacity due to very fast air diffusion, which is accounted for using a minimum redox depth, min.redoxdepth (metres). To account for unrealistically high values of the redox depth, a maximum redox depth is also estimated based on the principles of Hansen et al. (2014a). Hereby, the spatially distributed redox interface layer is quantified and incorporated into MIKE SHE (step four). The final step (five) is the calibration of the depth of the redox interface layer, also described as the calibration of the nitrate model, by adjusting the location of the layer. Step five will be described in Sect. 3.3.
3.2 DAISY setup and calibration (phase 2)
DAISY provides a one-dimensional finite difference description of the soil–water–crop–atmosphere processes (Abrahamsen and Hansen, 2000; Hansen et al., 1991), where flow is described using Richards' equation and, to a minor extent, by macro pore flow for loamy soils (B. Hansen et al., 2012). Nitrate transport is driven by the convection–dispersion algorithm in the soil matrix. The one-dimensional DAISY model can be used to represent an entire catchment by combining multiple columns of DAISY computations and evaluating the summed water balance. Calibration on the catchment scale is, however, not straightforward and must be done manually and in an incremental manner due to the large number of 1D model columns.
The DAISY water balance module used in this study is based on a previous calibration of the catchment (Børgesen et al., 2013), where root zone leaching and groundwater abstraction are compared with river discharge (Refsgaard et al., 2011; Børgesen et al., 2013) and, subsequently, extended on, as described in Karlsson et al. (2016) and in the following section. The model setup for the Odense catchment contains roughly 12 000 1D DAISY columns, and the model is set up so that each column represents unique combinations of soil type, climate, crop rotation, and groundwater depth. The DAISY model uses the same climate input and soil parameter setup as MIKE SHE, and the sensitive and calibrated unsaturated soil parameter from MIKE SHE was, therefore, transferred to DAISY. The water balance performance of DAISY was evaluated in the same calibration (2004–2007) and validation periods (2000–2003 and 2008–2009) as MIKE SHE.
In this study, the DAISY model was also used to simulate nitrate leaching for each soil column that represents a unique combination of soil type, climate, crop rotation, and groundwater depth. Crops are fertilized with mineral and organic nitrogen, depending on the farm type and soil type. The crop-recommended nitrogen rate based on soil type and crop sequence from the Danish Agricultural Agency (Plantedirektoratet, 2005), for the years 2004–2007, was used to set up the fertilization scheme. Nitrate leaching input is simulated on a daily basis and based on the leaching from the permutated crop rotations simulated for the dominating soil type within a 200 m × 200 m square grid (Karlsson et al., 2016). Because of the close feedback mechanism between nitrogen yields and nitrate leaching, the simulated mean nitrogen yields were recalibrated to observed annual mean nitrogen yields on Funen (Statistikbanken, 2015) for the dominating soil type for the period 2004–2007. The calibration is conducted by adjusting the crop parameters, following the methodology of Styczen et al. (2004). Nitrogen concentrations of yields were extracted from table values of mean nitrogen contents for different crops (Møller et al., 2005). For crop rotations including clover grass and peas, nitrogen biological fixation is calculated using Høgh-Jensen et al. (2004), and nitrogen atmospheric deposition is included as input to the soil using the standard DAISY settings (given in S. Hansen et al., 2012) for dry and wet deposition.
In the simulations under climate change, the effect of CO2 concentration change in the atmosphere has an impact on the light-saturated photosynthesis rate, which is a crop parameter in the DAISY crop model code. In order to deal with this feedback mechanism, the procedure to change the parameter was adopted from Børgesen and Olesen (2011).
3.3 Calibration of the nitrate model (phase 3)
Following the approach from Hansen et al. (2014b), the nitrate model is then constructed by combining the two models, i.e. DAISY and MIKE SHE. Daily values of nitrate flux from DAISY serve as input to MIKE SHE, where nitrate transport is simulated by converting nitrate input to particles using the particle-tracking module. Each time the accumulated input of nitrate reaches 0.5 kg N within the model cell (200 m × 200 m), a particle is released from the water table and is allowed to follow the groundwater flow. If the particle penetrates the redox interface, the nitrate is assumed to be removed completely and instantaneously by denitrification (Postma et al., 1991; Hansen et al., 2014a). Remaining particles will emerge in discharge zones typically located in stream valleys and leave the catchment at the river outlet (Kratholm station). The nitrate arrival percentage (NAP) is found to be the cumulative amount of nitrate leaving the catchment divided by the amount released at the water table.
The nitrate model is then calibrated by adjusting the depth to the redox interface through the calibration of f and min.redoxdepth to obtain the observed NAP (step five). As the calibration of these two parameters may result in non-uniqueness, all possible combinations (realizations) of the two parameters resulting in observed NAP are identified. For all realizations, the cumulative distribution of the redox depth is found at the location where observations of redox depth are available from boreholes and the cumulative distribution of the entire catchment. These two graphs are subsequently compared with the cumulative distribution of the actual measured redox depth in boreholes. The realization with the best representation of the fractional distribution of the observed redox depth for both on-site and, especially, catchment-scale values is chosen for the final redox depth parameters.
The reason for comparing calculated redox depths to cumulative distributions for actual measurement locations and the entire catchment distribution is due to several issues. First, measured redox depths are very local point measurements, and large variations in space (within a few metres) are often reported (e.g. Hansen et al., 2008; Ernstsen, 1996). Thus, a measurement may not be representative for the area or model scale, as numerous measurements together are more likely to represent to the correct fractional distributions in the catchment. Furthermore, the calculated redox depth may be applicable on the catchment scale, i.e. the scale on which it is also calibrated, but less trustworthy on the location scale.
3.4 Estimation of reduction map and map correction
The reduction map quantifies the nitrate reduction potential for each model grid (Hansen et al., 2014a). The number of representative particles released at each cell that is subsequently reduced is divided by the total number of particles released. The reduction map is, therefore, based on results from the nitrate model that is run with the calibrated redox interface. In the reduction map, a grid cell with the value of 100 % indicates that all particles released in the cell are subsequently reduced, while a value of 50 % indicates that half of all particles (nitrate) reach the surface waters unreduced. Therefore, it provides valuable information on what the environmental impact is, in terms of what the nitrate loading of farming practices is, for specific parts of the landscape.
Unfortunately, during the release of particles in the MIKE SHE model, numerical issues occasionally cause some particles to remain stuck in the unsaturated zone. MIKE SHE is a commercial modelling tool, and therefore, there is no possibility to access the modelling code in order to correct this numerical error or in any other way account for this model limitation. Therefore, it was necessary to introduce a correction scheme. The actual fate of these stuck particles (reduced/non-reduced) is unknown. At an early stage, the assumption was made that the captured particles, if they had moved correctly through the system, would be subject to a fate similar to the non-captured particles, i.e. that the relationship between reduced/non-reduced particles was the same. If this assumption is valid, then the calculation of the reduction potential in each grid cell is the same with/without the stuck particles. Unfortunately, this assumption may not always be valid. Furthermore, the arrival percentage estimated by the two methods is not the same, as not all particles are released in the complex particle arrival count, and the only way to calibrate the nitrate model is to use these data. For the two methods to be comparable, it is, therefore, necessary to exclude the particles that are stuck in the unsaturated zone. The correction factor is, therefore, introduced to eliminate the particles that are stuck and prevents them from causing changes in the reduction map. The correction uses a simple linear equation, where a correction factor is manually fitted so that the arrival percentage (originating from the reduction map multiplied by the nitrate input) matches the particle arrival percentage. These corrections are done individually for all reduction maps, and the correction causes a change in the reduction in the range of −7 % to 9 %, with a mean of 2 %.
3.5 Climate and land use scenarios
One emission scenario, the Intergovernmental Panel on Climate Change (IPCC) AR4 Special Report on Emission Scenarios (SRES) A1B (Nakicenovic et al., 2000), was chosen as the basis for this study. Since the study was conducted, newer generations of emission scenarios have been developed by the IPCC (van Vuuren et al., 2011) and are known as the Representative Concentration Pathways (RCPs). The A1B scenario is generally comparable to the RCP6.0 emission scenario (medium scenario).
In this study, realizations from four climate model combinations, i.e. the global climate model–regional climate model (GCM-RCM) couplings, were selected from the ENSEMBLES project (Hewitt and Griggs, 2004), where results from the period 2080–2099 were extracted and used as input to the hydrological model. The reference evapotranspiration is calculated using the Food and Agriculture Organization (FAO) Penman–Monteith formula, adapted by Allen et al. (1998), and based on the climate model outputs for minimum and maximum temperature, incoming long- and short-wave solar radiation, relative humidity, and wind speed. Following the recommendations in Allen et al. (1998) and Seaby et al. (2013), the variables needed for the Penman–Monteith formula, e.g. net radiation (calculated from the net incoming short- and long-wave radiation), water vapour pressure, height-adjusted wind speed, and atmospheric pressure, were calculated from these outputs. Precipitation data from both this period and the reference period, 1990–2009, were bias corrected (downscaled) using the distribution-based scaling (DBS) method, which is a direct method that preserves the dynamics and non-stationary nature of the raw climate model results (Seaby et al., 2013), while reference evapotranspiration was downscaled using a bias removal method (Seaby et al., 2013). The four selected realizations represent wet (+19 % increase in precipitation; ECHAM-HIRHAM5), dry (−11 % decrease in precipitation; ARPEGE–RM5.1), and warm (+3.4 ∘C temperature increase; HadCM3-HadRM3) realizations and a model representing a median projection, with +10 % in precipitation and +2.1 ∘C in temperature (ECHAM5–RCA3). The change factors can be seen in Table 1. Both climate models and bias corrections are described in more detail in Karlsson et al. (2016).
The four climate model realizations were combined with four land use scenarios and the baseline (present) land use scenario. The land use scenarios were created during workshops with researchers, farming industries, environmental protection agencies, and government representatives. During the workshops, participants identified possible paths of development for the land use in Denmark, considering the balance of agricultural value on one hand and priorities in the society on the other hand (e.g. environmental concerns or recreational use). From the workshop, four scenarios that describe agricultural management in the period 2080–2099 were created, as follows:
LU1 – agriculture for nature, where the agricultural area is reduced to 40 % of the land area through afforestation and increasing grass areas, and fertilization rates are generally reduced (−40 %).
LU2 – extensive agriculture, with a small 3 % point reduction in agricultural area, resulting in 64 % farmland; however, some of the intensive farm types (with high fertilization rates) are converted to less intensive farm types with less fertilization (total change of −60 %).
LU3 – high-tech agriculture, also with a small decrease in agricultural area of 3 %, but with crop production that is assumed to increase, resulting in an insignificant change in the needed fertilizer inputs (0 %).
LU4 – market-driven agriculture, where forest and some extensive farm types are converted into intensive farming, resulting in an agricultural area of 70 %. At the same time, fertilization rates are increased to reach maximum production levels (+20 %).
More information on the land use scenarios can be found in Olesen et al. (2014) and Karlsson et al. (2016).
All 20 combinations of future climate projections (4) and land use (5) were specified as input to the hydrological model. The model was run for both a future (2080–2099) and reference period (1990–2009), resulting in 40 scenarios. Additionally, the model was run with observed climate for the period 1990–2009 (five scenarios using observed land use and the four land use scenarios). Hence, a total of 45 model simulations were analysed (Table 2). In this paper, the following terminology is used:
Observational period – results from a hydrological model and the root zone model DAISY run forced with observational data in the period 1990–2009. This period covers both the calibration period (2004–2007) and validation periods (2000–2003 and 2008–2009), all of which are driven by observational data.
Reference period – this period is used to specifically describe the climate-model-driven hydrological results from the period 1990–2009. Future climate model runs are always compared with results from this period for the relevant climate model to ensure that climate model biases do not dominate the results.
Future period – the future scenario refers to climate-model-forced runs for the period 2080–2099.
Baseline – the term baseline refers to results from the specific model run combination (scenario 1; Table 3), where the current land use scenario (land use 0 – LU0) is combined with the observational climate data.
4.1 Model evaluation
The first phase of the calibration scheme was the MIKE SHE model calibration. A detailed presentation and evaluation of the water quantity performance of MIKE SHE have previously been described in Karlsson et al. (2014, 2016). MIKE SHE was reported to have good performance in the calibration and validation periods (Fig. S1 in the Supplement), with a water balance error of 3 % during calibration and 6 %–10 % in the validation periods, respectively, and an RMSE value of 1.5 m3 s−1 in the calibration period and 1.2–1.5 m3 s−1 during the validation period for the main station.
The second phase of the calibration scheme was the calibration of DAISY. Due to the complicated setup with a manual simultaneous calibration of more than 12 000 1D DAISY models, the water balance performance of the DAISY model was reported to be somewhat poorer, with water balance errors of 16 % in the calibration period and 3 %–23 % in the validation periods. The water quantity performance of DAISY is also described in more detail in Karlsson et al. (2016). The results from the subsequent manual calibration of the nitrate leaching in DAISY, using the N yields for the region, can be seen in Table 3. For all the crops in the region, the DAISY model is able to reproduce the observed harvested N within a margin of 0–5 kg of harvested nitrogen per hectare. The DAISY model is, therefore, able to represent the observed values of nitrate yields to a satisfactory level on the catchment scale. The catchment average root zone leaching from DAISY was 40 kg NO3-N ha−1 yr−1. As there are not direct measurements of leaching, this value cannot be directly verified. But Hansen et al. (2018) have reported values of the same order of magnitude (35.1 kg NO3-N ha−1 yr−1) for a different time period (2000–2010) using the NLES leaching model (Kristensen et al., 2003, 2008), while Højberg et al. (2015, 2017) reported values of 37.5 kg NO3-N ha−1 yr−1.
The third and final phase of the calibration is the calibration of the nitrate model or the calibration of the redox interface depth using the observed nitrate arrival percentage (NAP). The observed NAP of the catchment was estimated to be 35 %–39 %. The span is a result of the choice of time period; however, a NAP of 37 % was selected as the final calibration target. Several different combinations of the migration constant (f) and the minimum redox depth (min.redoxdepth) resulted in a NAP of 37 %. The cumulative distribution of the redox depths for all combinations are compared to the cumulative distribution of the observed redox depth (Fig. S2), and the parameter combination yielding the best representation of the fractional redox depth distribution is identified. Based on this analysis, the best combination with the correct NAP was found to be f=0.01 yr−1 and min.redoxdepth = 3 m. Figure 2 (left) shows the resulting depth to the redox interface for this combination. Generally, the depth to the redox interface is shallow, with a mean depth of 4.0 m (standard deviation of 5.6 m). This corresponds well with the previously reported redox depth (Koch et al., 2019b). With the calibrated redox interface depth, the resulting nitrate load in the river is 15 kg NO3-N ha−1 yr−1. This corresponds well with the observed loading of approximately 14 kg NO3-N ha−1 yr−1 measured at Kratholm station.
4.2 Baseline, reference, and future nitrate reduction map
The baseline reduction map was generated based on the nitrate input from the DAISY model for the present period and the current land use (Fig. 2; right). The mean nitrate reduction fraction in the catchment is 65 %, with a standard deviation of 38 %. The spatial distribution on the map shows high reduction potential in the uplands at the border of the model, likely where infiltrating water has a longer travel path to the river, and in areas where the redox interface is shallow. Lower reduction potential is seen in the areas near the stream network and in lowland areas with deep redox interface.
Utilizing the different nitrate input from the land use/climate model scenarios, 45 nitrate models (observational, reference, and future; Table 2) were run with the calibrated redox interface. The depth to the redox interface is assumed to be constant in time and is not updated for each scenario. This assumption is acceptable because of the slow migration of the interface, where the migration constant predicts 1 m of downward movement every 100 years. Following the procedure described in Sect. 3.4, this generates 45 nitrate reduction potential maps. The statistics for all 45 reduction potential maps are shown in Table 4, and a selection of these reduction maps can be seen in Fig. S3.
All scenarios show similar patterns of reduction zones with high and low removal fractions when compared to the baseline scenario map (Fig. 2). However, although the general pattern of the reduction maps is comparable, there are spatial differences between the maps. The redox depth remains constant, and all nitrate crossing this interface is assumed to be reduced, regardless of the amount. Therefore, the reason for these changes in the nitrate reduction potential map results from changes in the flow path and/or changes in the nitrate input from DAISY, reflected in the differences between drain/interflow versus groundwater flow to streams and the timing of the nitrate release from the root zone.
4.3 Impact of land use change on reduction maps
To investigate the impact of land use change on the reduction maps, only land use is changed while climate remains constant, and it is shown as the difference between land use change scenarios and the baseline scenario (Fig. 3; top row). The water balance of the models, and hence the groundwater level and water flow paths, is affected mainly by a possible change in evapotranspiration that is introduced with new crop rotation systems and vegetation types. Drains are still present in the entire catchment, regardless of land use, and are only active when they are required (submerged below the groundwater table). At the same time, the land use changes result in a nitrate input distribution and timing that may differ from the current land use (LU0). The changes are minor and give both higher and lower reduction potential, depending on changes in land use and vegetation for the individual grids (Fig. 4; top row).
The changes in the average catchment water flow components as a result of land use change is shown in Fig. 4 (left). This shows the change in each component, drain flow, overland flow, and base flow when changing from a scenario run with LU0 to the scenarios with LU1–4. These changes are minor for the overall water balance on the catchment scale for all land use change scenarios (Fig. 5; left).
The drain flow component is a primary conductor for non-reduced nitrate because it represents fast and shallow flows above the redox interface. It is, therefore, relevant to look at the spatial distribution of the changes in drain flow for the scenarios. While the bias corrections ensure that the climate models reproduce the overall mean and variances of the observed climate, they do not necessarily ensure consistency in the temporal structure of precipitation. Hence, the overall net precipitation (precipitation minus actual evapotranspiration) may change slightly across the climate models, and we have, therefore, plotted the change in drain flow fraction (Fig. 5), which is defined as the drain flow divided by net precipitation, instead of the drain flow component itself. As for the changes in the reduction map (Fig. 3), the differences between the reference and land use scenario drain flow fractions are very small and sporadic. The same is found for changes in recharge and groundwater head presented in Figs. S4 and S5 (top row). Averaged across the catchment, the change in land use gives rise to a maximum of 1.2 % change in the reduction potential (Fig. 3). For all 45 runs, the general statistics also show that changing the land use (horizontally; Table 4) does not change the mean reduction potential by more than a maximum of ±2 %.
4.4 Impact of land use and climate change on reduction maps
Figure 3 (bottom row) shows the results for the future climate compared to the corresponding reference period. For the future simulations, not only the water flow paths are affected by the change in amount and distribution of net precipitation but the vegetation uptake is also considerably different due to increasing potential evapotranspiration. The impact of these changes is evaluated in Fig. 4 (right), where the changes in the flow components are shown from the reference to the future period. The climate model projections have important impacts on the distribution of the water balance components, and substantial differences are found among the models. Even if the land use remains constant, the timing and the amount of nitrate leakage are subject to significant variations as a result of climate change. For the reduction maps (Fig. 3), these changes result in larger differences when compared to the maps produced for the reference period. Averaged across the catchment, changes of up to 7.3 % points are recorded for the four scenarios, and for all 45 scenarios, the change is up to ±8 %.
From a water balance perspective, the four future climate model projections vary greatly but may be grouped in two overall categories. The first category includes the wet models, ECHAM5-HIRHAM5 and ECHAM5-RCA3 (Table 1). Both show a small decrease in annual reference evapotranspiration. At the same time, precipitation is projected to increase in all seasons, resulting in an annual increase of 15 %–30 %. As a result, the net precipitation (Fig. 4) and groundwater recharge (Fig. S4) increase considerably. Since there is a limit to how fast water can be pushed through the deeper groundwater systems, the fast flow components (overland and drain flow) both increase substantially.
Even though the drain flow produced with the wet climate model projections generally increases, the drain flow fraction, defined as the drain flow divided by net precipitation (Fig. 5), shows two overall signals. In the eastern uphill area, the change in the drain flow fraction is positive, while in the central and western areas, in the river valleys, the change in the drain fraction is negative. This implies that less net precipitation, relatively speaking, is channelled through the drainage system than in the reference period in the river valleys (brown areas), while relatively more net precipitation is captured by drains in the eastern uphill locations. This highlights how the change in the distribution of the flow components also changes the spatial pattern across the catchment. The underlying explanation for the pattern recognized in the wet models is primarily found to be a moderate increase in upwelling water to the drainage system in the lowlands in the future compared to the substantial increase in drainage in the uphill areas.
The increase in groundwater recharge (Fig. S4) results in increasing groundwater levels (Fig. S5). However, the change varies across the catchment, and since flow is controlled by the gradient in hydraulic head, the non-homogeneous changes in heads will result in changes in the flow direction. The general tendency is that groundwater levels increase most in the upstream parts of the catchment, while remaining the same in the valleys near the stream. Hence, the gradients will become steeper, which results in an increase in the ratio of horizontal to vertical flow. Thus, this promotes near-surface flow paths in the subsurface.
As the fast flow components, like drain and overland flow, mainly carry non-reduced water, and the subsurface reduction depends on water being transported below the redox interface for nitrate reduction to take place, the increase in both fast flow components and the shallower flow paths for the infiltrated water is expected to give rise to a lower reduction potential, as is indeed observed in Fig. 3 and Table 4.
The second climate model category (dry models) includes ARPEGE-RM5.1 and HadCM3-HadRM3 (Table 1), which both show a significant increase in annual reference evapotranspiration of 10 %–20 % (Karlsson et al., 2016). With respect to precipitation, the two models show a slight decrease or no change in annual values, whereas lower precipitation is found during summer and autumn for both models. Therefore, the net precipitation decreases for both models, leading to a decrease in the fast flow components (overland and drain flow) for the future projections (Fig. 4). The reduction in net precipitation also results in a reversal of the distribution of the drain flow fraction (Fig. 5) and a decreasing groundwater recharge (Fig. S4).
These changes are reflected in the spatial distribution of the change in groundwater level (Fig. S5). As found for the wet models, the general tendency is that groundwater levels change more in the upstream parts of the catchment compared to the valleys near the stream. Therefore, the gradients will become less steep, resulting in a decrease in the ratio of horizontal to vertical flow. This leads to a higher degree of slower and deeper groundwater flow paths and, potentially, to more nitrate crossing the redox interface.
The decrease in the fast flow components and the deeper flow path supports the fact that the two models project a higher reduction potential in the future. This is indeed found for one of the models, HadCM3-HadRM3, which generally has the largest increase in reduction potential of the four models, from the average reduction of 63 % in the reference to 67 % in the future simulation (LU0; Fig. 3 and Table 4).
However, the other model, ARPEGE-RM5.1, shows a completely opposite signal with a lower reduction potential (Fig. 3). Initially, this anomaly did not seem to be explainable by the changing flow paths as the changes (Figs. 4, S4, and S5) are very similar to the other dry model. However, as mentioned previously, changes in the timing and quantities of the nitrate root zone leakage from DAISY are also a determining factor for the reduction map, and these dynamics are not always apparent on average maps like the ones shown in Figs. 3, 4, S4, and S5. A closer inspection of the monthly changes in the flow components and nitrate leakage provides a possible explanation for this phenomenon (calculations not shown).
The ARPEGE-RM5.1 model shows the absolute largest increase in nitrate leakage in January. Even though the model is generally dry, the January precipitation is projected to increase in the future. The dynamics in the model show a shift in flow components for this month towards a smaller amount of drain flow and recharge, while the overland flow increases, probably as a result of larger, prolonged rainfall events. This shift, combined with the very large increase in nitrate leakage, leads to more non-reduced water in the system for the future period than for the control period. As this combination of larger nitrate leakage and flow components shifts is not seen in the HadCM3-HadRM3 model, it is most likely a result of the lower reduction potential for ARPEGE-RM5.1 in the future, in spite of the overall drying signal from the model.
4.5 Effect on the nitrate flux using different reduction maps
To evaluate the impacts of using a fixed reduction map versus different maps calculated for each scenario explicitly considering differences caused by changes in land use and climate change, the total catchment nitrate load (nitrate arrival at Kratholm) was calculated using two different approaches for all scenarios.
A fixed reduction map (the baseline reduction map) is used and combined with the projected nitrate leaching from the root zone.
Targeted reduction maps, i.e. different reduction maps for each scenario (as calculated above), are used and combined with the projected nitrate leaching from the root zone.
For both approaches, the nitrate arrival is calculated for the observational, reference, and future periods. The resulting nitrate arrivals for approach 2 are then, for each case, compared to the corresponding scenario in approach 1, thereby illustrating the effect of using a targeted reduction map compared to a baseline reduction map (Fig. 6).
The grey bars in Fig. 6 show the relative change (compared to baseline), when applying targeted reduction maps for the four land use scenarios but using observed climate. Note, therefore, that the grey bars are the same for all four plots. The effect of the targeted reduction maps versus the fixed reduction map manifests in a limited spread in the estimated nitrate arrivals by only 1 %–3 %. This implies that, for a case of changing land use, the fixed reduction map is a reasonable approximation of the reduction potential in the catchment.
The green bars in Fig. 6 represent the same effect as for the grey bars (targeted versus fixed reduction maps); however, here the climate is the reference climate for each of the four climate models. For all models, the effect is larger here than when using the observed climate. This phenomenon is mainly due to the inherent bias of the reference climate model simulations. This becomes clear when looking at the first bars of each plot (red bars) which denote the change when only the observed climate is replaced by the reference climate but the current land use setup (LU0) is maintained. Even though the climate model output is bias corrected such that the general statistics in the reference period resemble those of the observations, there may still be differences in the temporal structure of the climate model outputs, which may impact the hydrological simulations. This is an important issue to bear in mind when analysing the future signal presented by the blue bars. However, the effect of land use can still be approximated by comparing the results from LU0 to LU1–4 for the different climate models. Again, the magnitude of the effect is between 1 %–3 % change in nitrate arrival; however, the models do not agree on which land use causes the largest change (being either LU1 or LU4).
For the blue bars, the differences from the fixed reduction map to the targeted reduction maps become a combination of bias correction limitations, land use change, and climate change. Even though it is not possible to completely separate the signal of these three components, a cautious estimation can again be achieved by subtracting the blue bars with the red bar result from LU0 in the reference period so that the removal of the signal from the climate model bias is attempted. For a dry model like ARPEGE-RM5.1, the effect (land use and climate) on the nitrate arrival is in the range of −5 % to −9 %, while for HadCM3–HadRM3, it is between −4 % and 4 %. The largest effects are found for the wet models, with changes ranging from 20 % to 27 % for ECHAM–HIRHAM, followed by ECHAM–RCA3 with 9 %–12 %. This shows that the consequences of using a fixed reduction map may be considerable, in particular with large changes in climatic conditions. Across all climate models, the average absolute effect on nitrate arrival is 10 % when using a fixed reduction map compared to a targeted reduction map.
5.1 Nitrate reduction maps are not constant in time
Our analysis clearly demonstrates that nitrate reduction maps are a result of complicated interactions between climate, vegetation, geology, and farm management, leading to a diversity of potential nitrate inputs, distributions, timing, flow paths, and reduction capability. This implies that nitrate reduction maps calculated for the present climatic conditions and flow patterns will differ from nitrate reduction maps under future climate and land use conditions. The main factor causing this is the differences in the precipitation/evapotranspiration regime that results in differences in how large of a fraction of the water percolating from the root zone reaches the stream via shallow flow routes above the redox interface, such as overland flow and runoff via drainpipes, and how large of a fraction of water takes a route through deeper groundwater zones and crosses the redox interface. Therefore, nitrate reduction maps also differ between a wet year and a dry year in the present climate.
Compared to climate effects, the impacts of land use change are minor because land use change does not affect the flow regime to the same extent that variability and change in climate do. This novel finding has not been recognized in previous studies of nitrate reduction maps, such as Hansen et al. (2014b), Højberg et al. (2015), Andersen et al. (2016), and Refsgaard et al. (2019). Nevertheless, it is important to note that the distribution of drains in the catchment was not changed during the land use change scenarios. Hence, uniform drain distribution and parameterization are assumed where the drainage component covers both natural (ditches and small canals) and agricultural drains. A different approach could have been adopted in the land use change scenarios by changing the drainage efficiency (by adjusting drainage parameters) in re- or deforested areas. Unfortunately, little information is available to guide this fine-tuning. However, the influence of land use change on the water balance and, therefore, the nitrate reduction is expected to increase if this effect is accounted for.
To encapsulate the range of uncertainty and influences from climate and land use scenarios in this setup, all scenario combinations were used. However, one could speculate that not all combinations of land use and climate change scenarios may be equally likely or plausible in the future, as decisions on land use application made by local farmers or through national regulations are made concurrently to adapt to or mitigate changes in climatic conditions.
5.2 Water management implications
The advantages of assuming a fixed reduction map are that any projected nitrate input, regardless of climate and land use scenario, can be multiplied with the reduction map and, thus, provide a projected nitrate outflow estimate with little effort and time spent. However, as shown in the present study, the assumption that the reduction map is constant in a changing environment is problematic. The good question then is how large errors are made by assuming a fixed reduction map and how this should be dealt with in water management practices.
The analysis indicates that assuming fixed reduction maps leads to small errors when dealing with land use change impacts but may lead to substantial errors (mean of 10 % on catchment nitrate load) when climate change projections are included. Land use change impacts may, however, be underestimated as a result of the uniform drainage setup. A 10 % error on the reduction map may potentially have major impacts on practical water management. Considering, for instance, the baseline scenario in Table 4, where the average nitrate reductions vary between 55 % and 67 % reduction, this implies that the net impact of a 100 kg N reduction in leaching from the root zone will vary between 45 and 33 kg (i.e. 30 %). Such changes are larger than the effects of sophisticated mitigation measures (Hansen et al., 2017), and whether such errors are acceptable depends on the purpose and context in specific water management situations. Thus, using fixed reduction maps may well be justifiable for initial screening purposes, while targeted reduction maps, explicitly calculated for specific scenarios, may be required for the design of remediation measures that will have significant socioeconomic impacts for stakeholders. The uncertainty of using a fixed reduction map for future scenarios should, of course, be seen in the context of the inherent uncertainties of the nitrate reduction maps (Hansen et al., 2014b).
5.3 Uncertainties and limitations
When using a hydrological model for simulating the impacts of changes in catchment conditions compared to those existing in the calibration period, split sample validation tests are not sufficient to document a model's capability to simulate hydrological changes. Experience shows that models that are used for making predictions beyond conditions for which they are calibrated (such as land use or climate change in the present study), often suffer from model structural uncertainties (Refsgaard et al., 2012), equifinality, or parameters that may not be transferrable in time (Thirel et al., 2015). In such situations, the more comprehensive and data-demanding differential split sample tests are recommended (Refsgaard et al., 2014; Klemeš, 1986). Due to lack of data, such tests were beyond the scope of the present study. Instead, the model structural uncertainty for the present case was assessed using a multimodel approach with two additional hydrological models (Karlsson et al., 2016), suggesting that the signal coming from climate change was dominating over the model structural uncertainty as far as hydrological change is concerned. Therefore, we argue that the inevitable uncertainties arising from model use beyond calibration conditions are most likely not so large that they affect our conclusions.
The present study was carried out for a groundwater-dominated catchment characterized by till deposits, confined aquifers, and relatively shallow redox interfaces and phreatic groundwater tables. Furthermore, this catchment has a relatively uniform soil type distribution that is dominated by clayey soils. We consider the conclusions to be applicable to catchments with similar hydrogeological conditions, while they cannot be used in groundwater-dominated catchments characterized by alluvial plains without fast flow components such as overland flow and drainpipe/ditch flow and aerobic groundwater systems without a redox interface. Similarly, the conclusions cannot be transferred to surface-water-dominated catchments, where the nitrate reduction takes place in streams and lakes, although we suspect that nonlinearities may cause similar effects here.
In this study, the location of the redox interface was assumed to be the same for both present and future scenarios. While this assumption may be reasonable in relation to the slow natural migration of the redox interface caused by the percolation of oxygen, the redox interface migration may be escalated by nitrate application on the land surface, as reported by, e.g., Böhlke et al. (2002) and Wriedt and Rode (2006). This issue cannot be addressed within the framework of this study, but the effect on the redox interface may be substantial, especially for land use scenarios with high nitrate application.
The indication that errors can be up to 10 % is based on only a single case study with one catchment, one model, and a limited number of land use and climate change scenarios. While similar results may be found when applying the same approach for catchments governed by the same dominant flow processes and land use types like the one investigated in this study, the error must be expected to be site and context specific, and therefore, this causes projection uncertainties that should be addressed along with other known sources of uncertainty such as climate model projections, land use projections, and parameter uncertainties, including the effect of equifinality, geological uncertainty, and hydrological model structural uncertainty (Hansen et al., 2014b; Karlsson et al., 2016).
During the calibration and validation of the model, a decrease in model performance was registered in the water balance for the validation period. This could be caused by non-optimal parameter estimates, and there is always a risk of equifinality during the model calibration. In this case, the risk was minimized by using an extensive dataset of both discharge, hydraulic head, redox depth, and nitrate flux during the calibration of the different model steps. However, a multimodel set, as used in this study, may still be prone to the risk of equifinality. Parameter estimations are done here in a stepwise fashion for each of the models, and the catchment-scale calibration of DAISY, along with the particle tracking approach, limits the evaluation of the performance of the nitrate component to mean catchment values. The dynamic of the nitrate system is, thus, impossible to verify. To account for this, a full solute transport solution would be necessary but was unfortunately not possible in the framework of this study. However, it would be a relevant next step in investigating uncertainties and improving the model verification.
During simulations it was found that some particles were not correctly released due to model error, and we, therefore, were forced to make the assumption that the particles would be distributed in a similar manner as non-trapped particles. However, the validity of this assumption is associated with considerably uncertainty. The correction led to mean changes and reduction of 2 %, but the resulting impact on the true reduction map is unknown, as we do not know the actual travel path of the trapped particles.
Nitrate reduction maps are valuable tools used for the calculation of remediation and climate change effects on nitrate leaching and are generally considered constant in time, even though the timing of nitrate leaching and flow paths may change. In this study, we investigate the potential consequences for estimating the nitrate climate and land use change impact projections when assuming a fixed reduction map. For an agriculture-dominated catchment in Denmark, the DAISY model was used to provide nitrate leaching input, while the hydrological model MIKE SHE was used to simulate the flow regime and nitrate flow path through particle tracking. In total, four land use scenarios and four climate change projections were evaluated. The main findings of the study are as follows:
Changing climate conditions lead to reduction map changes of around 10 %, while effects from land use changes were minor. However, land use effects may be underestimated due to drainage formulations in non-agricultural areas.
The magnitude of the changes in the reduction map found here may, however, be influenced by both the model setup (e.g. drainage) and the model errors (e.g. particle flow paths) and assumptions (e.g. fixed redox interface). Furthermore, the span of the chosen land use and climate change scenarios analysed and the flow regime in the study catchment may also influence the results.
The error will, therefore, be specific for the study site and context, and it should, consequently, be tackled along with other sources of uncertainty, like geological, parameter, and model structure uncertainties that are not evaluated in this study.
The hydrological model is based on the commercial software of MIKE SHE, and therefore, the model code is not available.
Data from the study are available upon request and through the public server at https://doi.org/10.22008/FK2/YEMDIS (Seidenfaden, 2022). An exception to this is the observed climate data input, which is currently the property of the Danish Meteorological Institute (DMI) but will be made publicly available through https://www.dmi.dk/frie-data/ (DMI, 2022) before end of 2023.
The supplement related to this article is available online at: https://doi.org/10.5194/hess-26-955-2022-supplement.
IKS performed model simulations, wrote part of the paper, and produced the figures. TSO and JCR contributed to the formulation of the conceptualization, methodology, and simulation assistance, as well as the writing of the paper. DT, CDB, and JEO contributed to the methodology and the writing of the paper. CDB and JEO, additionally, provided the input data and the results for the paper.
The contact author has declared that neither they nor their co-authors have any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
The present study has been carried out as part of the project of Centre for Regional Change in the Earth System (CRES; http://cres-centre.dk/frontpage, last access: 7 February 2022).
This research has been supported by the Strategiske Forskningsråd (Danish Strategic Research Council; grant no. DSF-EnMi 09-066868).
This paper was edited by Rohini Kumar and reviewed by Pia Ebeling, Fanny Sarrazin, and one anonymous referee.
Abbott, M. B., Bathurst, J. C., Cunge, J. A., O'Connell, P. E., and Rasmussen, J.: An introduction to the European hydrological system – systeme hydrologique europeen, “She”, 2: Structure of a physically-based, distributed modeling system, J. Hydrol., 87, 61–77, 1986.
Abrahamsen, P. and Hansen, S.: Daisy: an open soil-crop-atmosphere system model, Environ. Model. Softw., 15, 313–330, 2000.
Allen, R. G., Pereira, L. S., Raes, D., and Smith, M.: FAO Penman–Monteith equation, in: Crop evapotranspiration – Guidelines for computing crop water requirements, edited by: FAO, FAO Irrigation and drainage paper 56, Food and Agriculture Organization of the United Nations, Rome, Italy, 15–18, ISBN 92-5-104219-5, 1998.
Andersen, H. E., Blicher-Mathiesen, G., Thodsen, H., Andersen, P. M., Larsen, S. E., Stålnacke, P., Humborg, C., Mörth, C.-M., and Smedberg, E.: Identifying Hot Spots of Agricultural Nitrogen Loss Within the Baltic Sea Drainage Basin, Water Air Soil Pollut., 227, 38, https://doi.org/10.1007/s11270-015-2733-7, 2016.
Böhlke, J. K., Wanty, R., Tuttle, M., Delin, G., and Landon, M.: Denitrification in the recharge area and discharge area of a transient agricultural nitrate plume in a glacial outwash sand aquifer, Minnesota, Water Resour. Res., 38, 10-11–10-26, https://doi.org/10.1029/2001WR000663, 2002.
Børgesen, C. D. and Olesen, J. E.: A probabilistic assessment of climate change impacts on yield and nitrogen leaching from winter wheat in Denmark, Nat. Hazards Earth Syst. Sci., 11, 2541–2553, https://doi.org/10.5194/nhess-11-2541-2011, 2011.
Børgesen, C. D., Jensen, P. N., Blicher-Mathiesen, G., Schelde, K., Grant, R., Vinther, F. P., Thomsen, I. K., Hansen, E. M., Kristensen, I. T., Sørensen, P., and Poulsen, H. D.: Udviklingen i kvælstofudvaskning of næringsstofoverskud fra dansk landbrug for perioden 2007–2011. Evaluering af implementerede virkemidler til reduktion af kvælstofudvaskning samt en fremskrivning af planlagte virkemidlers effekt frem til 2015, DCA Report 031, DCA – Nationalt Center for Fødevarer og Jordbrug, Tjele, Denmark, URL: https://pure.au.dk/ws/files/68362856/dcarapporten31.pdf (last access: 7 February 2022), 2013.
Bruun, S., Christensen, B. T., Hansen, E. M., Magid, J., and Jensen, L. S.: Calibration and validation of the soil organic matter dynamics of the Daisy model with data from the Askov long-term experiments, Soil Biol. Biochem., 35, 67–76, https://doi.org/10.1016/S0038-0717(02)00237-7, 2003.
DHI: DHI MIKE SHE [Computer software], User Guide, DHI – Water & Environment, Hørsholm, Denmark, https://www.mikepoweredbydhi.com/products/mike-she (last access: 7 February 2022), 2019.
Diaz, R. J. and Rosenberg, R.: Spreading dead zones and consequences for marine ecosystems, Science, 321, 926–929, https://doi.org/10.1126/science.1156401, 2008.
DMI: Climate data, DMI Open data, https://www.dmi.dk/frie-data/, last access: 7 February 2022.
Ernstsen, V.: Reduction of Nitrate By Fe2+ in Clay Minerals, Clays Clay Miner., 44, 599–608, https://doi.org/10.1346/CCMN.1996.0440503, 1996.
Ernstsen, V. and Mørup, S.: Nitrate reduction in clayey till by Fe(II) in clay minerals, Hyperfine Interact., 70, 1001–1004, https://doi.org/10.1007/BF02397497, 1992.
Ernstsen, V., Højberg, A. L., Jakobsen, P. R., von Platen, F., Tougaard, L., Hansen, J. R., Blicher-Mathiasen, G., Bøgestrand, J., and Børgesen, C. D.: Beregning af nitratreduktionsfaktorer for zonen mellem rodzonen og frem til vandløbet, Data og metode for 1. generationskortet, GEUS report 2006/93, Danmarks og Grønlands Geologiske Undersøgelse, Copenhagen, Denmark, 2006.
Ernstsen, V., Jakobsen, P. R., and von Platen, F.: Et første landsdækkende redoxkort, Vand Jord, 4, 159–160, http://vand-og-jord.dk/et-foerste-landsdaekkende-redoxkort/ (last acces: 8 February 2022), 2010.
Fleck, S., Ahrends, B., Sutmöller, J., Albert, M., Evers, J., and Meesenburg, H.: Is Biomass Accumulation in Forests an Option to Prevent Climate Change Induced Increases in Nitrate Concentrations in the North German Lowland?, Forests, 8, 1–28, https://doi.org/10.3390/f8060219, 2017.
Graham, D. N. and Butts, M. B.: Flexible, integrated watershed modelling with MIKE SHE, in: Watershed Models, edited by: Singh, V. P. and Frevert, D. K., CRC Press, 245–272, https://doi.org/10.1201/9781420037432.ch10, 2005.
Greve, M. H., Greve, M. B., Bøcher, P. K., Balstrøm, T., Madsen, H. B., and Krogh, L.: Generating a Danish raster-based topsoil property map combining choropleth maps and point information, Geografisk Tidsskrift – Danish Journal of Geography, 107, 1–12, https://doi.org/10.1080/00167223.2007.10649565, 2007.
Hansen, A. L., Christensen, B. S. B., Ernstsen, V., He, X., and Refsgaard, J. C.: A concept for estimating depth of the redox interface for catchment-scale nitrate modelling in a till area in Denmark, Hydrogeol. J., 22, 1639–1655, https://doi.org/10.1007/s10040-014-1152-y, 2014a.
Hansen, A. L., Gunderman, D., He, X., and Refsgaard, J. C.: Uncertainty assessment of spatially distributed nitrate reduction potential in groundwater using multiple geological realizations, J. Hydrol., 519, 225–237, https://doi.org/10.1016/j.jhydrol.2014.07.013, 2014b.
Hansen, A. L., Refsgaard, J. C., Olesen, J. E., and Børgesen, C. D.: Potential benefits of a spatially targeted regulation based on detailed N-reduction maps to decrease N-load from agriculture in a small groundwater dominated catchment, Sci. Total Environ., 595, 325–336, https://doi.org/10.1016/j.scitotenv.2017.03.114, 2017.
Hansen, A. L., Donnelly, C., Refsgaard, J. C., and Karlsson, I. B.: Simulation of nitrate reduction in groundwater – An upscaling approach from small catchments to the Baltic Sea basin, Adv. Water Resour., 111, 58–69, https://doi.org/10.1016/j.advwatres.2017.10.024, 2018.
Hansen, B., Dalgaard, T., Thorling, L., Sørensen, B., and Erlandsen, M.: Regional analysis of groundwater nitrate concentrations and trends in Denmark in regard to agricultural influence, Biogeosciences, 9, 3277–3286, https://doi.org/10.5194/bg-9-3277-2012, 2012.
Hansen, J. R., Ernstsen, V., Refsgaard, J. C., and Hansen, S.: Field scale heterogeneity of redox conditions in till-upscaling to a catchment nitrate model, Hydrogeol. J., 16, 1251–1266, https://doi.org/10.1007/s10040-008-0330-1, 2008.
Hansen, J. R., Refsgaard, J. C., Ernstsen, V., Hansen, S., Styczen, M., and Poulsen, R. N.: An integrated and physically based nitrogen cycle catchment model, Hydrol. Res., 40, 347–363, 2009.
Hansen, S., Jensen, H. E., Nielsen, N. E., and Svendsen, H.: Simulation of nitrogen dynamics and biomass production in winter wheat using the Danish simulation model DAISY, Fertiliz. Res., 27, 245–259, https://doi.org/10.1007/BF01051131, 1991.
Hansen, S., Abrahamsen, P. T., Petersen, C., and Styczen, M.: Daisy: Model Use, Calibration, and Validation, T. ASABE, 55, 1317, https://doi.org/10.13031/2013.42244, 2012.
Hewitt, C. D., and Griggs, D. J.: Ensembles-based predictions of climate changes and their impacts, Eos Trans. Am. Geophys. Union, 85, 1–566, 2004.
Hoang, L., v. Griensven, A., v. d. Keur, P., Troldborg, L., Nilsson, B., and Mynett, A.: Comparison of the SWAT model versus DAISY-MIKE SHE model for simulating the flow and nitrogen processes, in: The International SWAT conference, Seoul, Korea, 306–319, https://swat.tamu.edu/media/33774/swat2010-proceedings.pdf (last access: 7 February 2022), 2010.
Høgh-Jensen, H., Loges, R., Jørgensen, F. V., Vinther, F. P., and Jensen, E. S.: An empirical model for quantification of symbiotic nitrogen fixation in grass-clover mixtures, Agricult. Syst., 82, 181–194, https://doi.org/10.1016/j.agsy.2003.12.003, 2004.
Højberg, A. L., Troldborg, L., Nyegaard, P., Ondracek, M., Stisen, S., and Christensen, B. S. B.: DK-model2009 – Sammenfatning af opdateringen 2005–2009, Geological Survey of Denmark and Greenland, 1–39, https://vandmodel.dk/media/8078/dk-model2009_sammenfatning.pdf (last access: 7 February 2022), 2010.
Højberg, A. L., Troldborg, L., Stisen, S., Christensen, B. B. S., and Henriksen, H. J.: Stakeholder driven update and improvement of a national water resources model, Environ. Model. Softw., 40, 202–213, https://doi.org/10.1016/j.envsoft.2012.09.010, 2013.
Højberg, A. L., Windolf, J., Børgesen, C. D., Troldborg, L., Tornbjerg, H., Blicher-Mathiasen, G., Kronvang, B., Thodsen, H., and Ernstsen, V.: National kvælstofmodel – Oplandsmodel til belastning og virkemidler. Metode rapport, revideret udgave september 2015 (National nitrogenmodel – Catchment model for load and measures. Method report, revised version September 2015), https://www.geus.dk/media/7744/national-kvaelstofmodel-oplandsmodel-til-belastning-og- (last access: 7 February 2022), 2015.
Højberg, A. L., Hansen, A. L., Wachniew, P., Żurek, A. J., Virtanen, S., Arustiene, J., Strömqvist, J., Rankinen, K., and Refsgaard, J. C.: Review and assessment of nitrate reduction in groundwater in the Baltic Sea Basin, J. Hydrol.: Reg. Stud., 12, 50–68, https://doi.org/10.1016/j.ejrh.2017.04.001, 2017.
Huno, S. K. M., Rene, E. R., van Hullebusch, E. D., and Annachhatre, A. P.: Nitrate removal from groundwater: a review of natural and engineered processes, J. Water Supply, 67, 885–902, https://doi.org/10.2166/aqua.2018.194, 2018.
Karlsson, I. B., Sonnenborg, T. O., Refsgaard, J. C., Trolle, D., Børgesen, C. D., Olesen, J. E., Jeppesen, E., and Jensen, K. H.: Significance of hydrological model choice and land use changes when doing climate change impact assessment, in: EGU General Assembly 2014, Vienna, Austria, 2014.
Karlsson, I. B., Sonnenborg, T. O., Seaby, L. P., Jensen, K. H., and Refsgaard, J. C.: The climate change impact of a high-end CO2-emission scenario on hydrology, Clim. Res., 64, 39–54, https://doi.org/10.3354/cr01265, 2015.
Karlsson, I. B., Sonnenborg, T. O., Refsgaard, J. C., Trolle, D., Børgesen, C. D., Olesen, J. E., Jeppesen, E., and Jensen, K. H.: Combined effects of climate models, hydrological model structures and land use scenarios on hydrological impacts of climate change, J. Hydrol., 535, 301–317, https://doi.org/10.1016/j.jhydrol.2016.01.069, 2016.
Klemeš, V.: Operational testing of hydrological simulation models, Hydrolog. Sci. J., 31, 13–24, https://doi.org/10.1080/02626668609491024, 1986.
Knoll, L., Breuer, L., and Bach, M.: Nation-wide estimation of groundwater redox conditions and nitrate concentrations through machine learning, Environ. Res. Lett., 15, 064004, https://doi.org/10.1088/1748-9326/ab7d5c, 2020.
Koch, J., Ernstsen, V., and Højberg, A. L.: Dybden til redoxgrænsen, 100 m grid, GEUS, Copenhagen [data set] https://data.geus.dk/MetaVis/Klik.jsp?id=2209&lang=da (last access: 7 February 2022), 2019a.
Koch, J., Stisen, S., Refsgaard, J. C., Ernstsen, V., Jakobsen, P. R., and Højberg, A. L.: Modeling Depth of the Redox Interface at High Resolution at National Scale Using Random Forest and Residual Gaussian Simulation, Water Resour. Res., 55, 1451–1469, https://doi.org/10.1029/2018wr023939, 2019b.
Kristensen, K., Jørgensen, U., and Grant, R.: Genberegning af modellen N-LES. Baggrundsnotat til Vandmiljøplan II – slutevaluering (Recalculation of the model N-LES. Background note to Aquatic Environmental Plan II – final evaluation), https://www.dmu.dk/1_viden/2_Publikationer/3_Ovrige/rapporter/VMPII/Genberegning_af_modellen_NLES.pdf (last access: 7 February 2022), 2003.
Kristensen, K., Waagepetersen, J., Børgesen, C. D., Vinther, F. P., Grant, R., and Blicher-Mathiesen, G.: Reestimation and further development in the model N-LES – N-LES3 to N-LES4, Series: DJF Plant Science No. 13, Det Jordbrugsvidenskabelige Fakultet, Aarhus Universitet, Aarhus, 1–30, https://pure.au.dk/ws/files/2423945/DJF_rapport_Markbrug_139 (last access: 7 February 2022), 2008.
Kristensen, K. J. and Jensen, S. E.: A Model for Estimating Actual Evapotranspiration from Potential Evapotranspiration, Hydrol. Res., 6, 170–188, 1975.
Kunkel, R., Eisele, M., Schäfer, W., Tetzlaff, B., and Wendland, F.: Planning and implementation of nitrogen reduction measures in catchment areas based on a determination and ranking of target areas, Desalination, 226, 1–12, https://doi.org/10.1016/j.desal.2007.01.231, 2008.
Madsen, H.: Automatic calibration of a conceptual rainfall-runoff model using multiple objectives, J. Hydrol., 235, 276–288, https://doi.org/10.1016/S0022-1694(00)00279-1, 2000.
Mas-Pla, J. and Menció, A.: Groundwater nitrate pollution and climate change: learnings from a water balance-based analysis of several aquifers in a western Mediterranean region (Catalonia), Environ. Sci. Pollut. Res. Int., 26, 2184–2202, https://doi.org/10.1007/s11356-018-1859-8, 2019.
Merz, C., Steidl, J., and Dannowski, R.: Parameterization and regionalization of redox based denitrification for GIS-embedded nitrate transport modeling in Pleistocene aquifer systems, Environ. Geol., 58, 1587–1599, https://doi.org/10.1007/s00254-008-1665-6, 2009.
Møller, J., Thøgersen, R., Helleshøj, M. E., Weisbjerg, M. R., Søegaard, K., and Hvelplund, T.: Fodermiddeltabel – Sammensætning og foderværdi af fodermidler til kvæg, volume 112, Report no. 69, Dansk Landbrugsrådgivning, Landscentret, 64 pp., ISSN 0906-1398, https://docplayer.dk/104824139-Rapport-112-fodermiddel-tabel-sammensaetning- (last access: 7 February 2022), 2005.
Nakicenovic, N., Alcamo, J., Davis, G., de Vries, B., Fenhann, J., Gaffin, S., Gregory, K., Grübler, A., Jung, T. Y., Kram, T., Rovere, E. L. L., Michaelis, L., Mori, S., Morita, T., Pepper, W., Pitcher, H., Price, L., Riahi, K., Roehrl, A., Rogner, H., Sankovski, A., Schlesinger, M., Shukla, P., Smith, S., Swart, R., van Rooijen, S., Victor, N., and Dadi, Z.: Emission Scenarios, in: Emission Scenarios, edited by: Nakicenovic, N. and Swart, R., Cambridge University Press, Cambridge, UK, https://www.ipcc.ch/site/assets/uploads/2018/03/emissions_scenarios-1.pdf (last access: 11 February 2022), 2000.
Nielsen, K., Stjernholm, M., Olsen, B. Ø., Müller-Wohlfeil, D.-I., Madsen, I.-L., Kjeldgaard, A., Groom, G., Hansen, H. S., Rolev, A. M., Hermansen, B., Skov-Petersen, H., Johannsen, V. K., Hvidberg, M., Jensen, J. E., Bacher, V., and Larsen, H.: Areal Informations Systemet – AIS, National Environmental Research Institute, Denmark [data set], https://www.dmu.dk/1_Viden/2_Miljoe-tilstand/3_samfund/AIS/2_Rapport/ais_rapport.pdf (last access: 7 Feburay 2022), 2000.
Olesen, J. E., Jeppesen, E., Porter, J. R., Børgesen, C. D., Trolle, D., Refsgaard, J. C., Sonnenborg, T., and Karlsson, I. B.: Scenarier for fremtidens arealanvendelse i Danmark, Vand og Jord, 21, 126–129, 2014.
Olesen, J. E., Børgesen, C. D., Hashemi, F., Jabloun, M., Bar-Michalczyk, D., Wachniew, P., Zurek, A. J., Bartosova, A., Bosshard, T., Hansen, A. L., and Refsgaard, J. C.: Nitrate leaching losses from two Baltic Sea catchments under scenarios of changes in land use, land management and climate, Ambio, 48, 1252–1263, https://doi.org/10.1007/s13280-019-01254-2, 2019.
Ortmeyer, F., Mas-Pla, J., Wohnlich, S., and Banning, A.: Forecasting nitrate evolution in an alluvial aquifer under distinct environmental and climate change scenarios (Lower Rhine Embayment, Germany), Sci. Total Environ., 768, 144463, https://doi.org/10.1016/j.scitotenv.2020.144463, 2021.
Paradis, D., Vigneault, H., Lefebvre, R., Savard, M. M., Ballard, J. M., and Qian, B.: Groundwater nitrate concentration evolution under climate change and agricultural adaptation scenarios: Prince Edward Island, Canada, Earth Syst. Dynam., 7, 183–202, https://doi.org/10.5194/esd-7-183-2016, 2016.
Plantedirektoratet: Vejledning om gødskning- og Harmoniregler, Danish Plant Directorate, https://lbst.dk/landbrug/goedning/vejledning-om-goedsknings-og-harmoniregler/ (last access: 7 February 2022), 2005.
Postma, D., Boesen, C., Kristiansen, H., and Larsen, F.: Nitrate Reduction in an Unconfined Sandy Aquifer: Water Chemistry, Reduction Processes, and Geochemical Modeling, Water Resour. Res., 27, 2027–2045, https://doi.org/10.1029/91wr00989, 1991.
Quick, A. M., Reeder, W. J., Farrell, T. B., Tonina, D., Feris, K. P., and Benner, S. G.: Nitrous oxide from streams and rivers: A review of primary biogeochemical pathways and environmental variables, Earth-Sci. Rev., 191, 224–262, https://doi.org/10.1016/j.earscirev.2019.02.021, 2019.
Refsgaard, J. C., Stisen, S., Højberg, A. L., Olsen, M., Henriksen, H. J., Børgesen, C. D., Vejen, f., Kern-Hansen, C., and Blicher-Mathiesen: Vandbalance i Danmark – Vejledning i opgørelse af vandbalance ud fra hydrologiske data for perioden 1990–2010, GEUS report 2011/77, Copenhagen, Denmark, https://doi.org/10.22008/gpub/28804, 2011.
Refsgaard, J. C., Christensen, S., Sonnenborg, T. O., Seifert, D., Højberg, A. L., and Troldborg, L.: Review of strategies for handling geological uncertainty in groundwater flow and transport modeling, Adv. Water Resour., 36, 36–50, https://doi.org/10.1016/j.advwatres.2011.04.006, 2012.
Refsgaard, J. C., Madsen, H., Andréassian, V., Arnbjerg-Nielsen, K., Davidson, T. A., Drews, M., Hamilton, D. P., Jeppesen, E., Kjellström, E., Olesen, J. E., Sonnenborg, T. O., Trolle, D., Willems, P., and Christensen, J. H.: A framework for testing the ability of models to project climate change and its impacts, Climatic Change, 122, 271–282, https://doi.org/10.1007/s10584-013-0990-2, 2014.
Refsgaard, J. C., Hansen, A. L., Højberg, A. L., Olesen, J. E., Hashemi, F., Wachniew, P., Wörman, A., Bartosova, A., Stelljes, N., and Chubarenko, B.: Spatially differentiated regulation: Can it save the Baltic Sea from excessive N-loads?, Ambio, 48, 1278–1289, https://doi.org/10.1007/s13280-019-01195-w, 2019.
Reusch, T. B. H., Dierking, J., Andersson, H. C., Bonsdorff, E., Carstensen, J., Casini, M., Czajkowski, M., Hasler, B., Hinsby, K., Hyytiäinen, K., Johannesson, K., Jomaa, S., Jormalainen, V., Kuosa, H., Kurland, S., Laikre, L., MacKenzie, B. R., Margonski, P., Melzner, F., Oesterwind, D., Ojaveer, H., Refsgaard, J. C., Sandström, A., Schwarz, G., Tonderski, K., Winder, M., and Zandersen, M.: The Baltic Sea as a time machine for the future coastal ocean, Sci. Adv., 4, eaar8195, https://doi.org/10.1126/sciadv.aar8195, 2018.
Seaby, L. P., Refsgaard, J. C., Sonnenborg, T. O., Stisen, S., Christensen, J. H., and Jensen, K. H.: Assessment of robustness and significance of climate change signals for an ensemble of distribution-based scaled climate projections, J. Hydrol., 486, 479–493, 2013.
Seidenfaden, I. K.: Redoxmaps from CRES2016, GEUS, GEUS dataverse [data set], https://doi.org/10.22008/FK2/YEMDIS, 2022.
Sjøeng, A. M. S., Kaste, Ø., and Wright, R. F.: Modelling future NO3 leaching from an upland headwater catchment in SW Norway using the MAGIC model: II. Simulation of future nitrate leaching given scenarios of climate change and nitrogen deposition, Hydrol. Res., 40, 217–233, https://doi.org/10.2166/nh.2009.068, 2009.
Statistikbanken: Statistical regional registrated annual mean yields, Statistikbanken [data set], https://www.statistikbanken.dk/jord3 (last access: 7 February 2022), 2015.
Styczen, M., Hansen, S., Jensen, L. S., Svendsen, H., Abrahamsen, P., Børgesen, C. D., Thirup, C., and Østergaard, H. S.: Standardopstillinger til Daisy-modellen. Vejledning og baggrund, Institut for Vand og Miljø, DHI, 1–62, https://daisy.ku.dk/about-daisy/dansk-standart/daisystaabirapport_240406.pdf (last acess: 7 February 2022), 2004.
Tesoriero, A. J., Terziotti, S., and Abrams, D. B.: Predicting Redox Conditions in Groundwater at a Regional Scale, Environ. Sci. Technol., 49, 9657–9664, https://doi.org/10.1021/acs.est.5b01869, 2015.
Thirel, G., Andréassian, V., and Perrin, C.: On the need to test hydrological models under changing conditions, Hydrolog. Sci. J., 60, 1165–1173, https://doi.org/10.1080/02626667.2015.1050027, 2015.
Troldborg, L., Højberg, A. L., Nyegaard, P., Stisen, S., Christensen, B. S. B., and Ondracek, M.: DK-model2009 – Modelopstilling og kalibrering for Fyn, https://vandmodel.dk/media/8082/dk-model2009_fyn.pdf (last access: 7 February 2022), 2010.
Trolle, D., Nielsen, A., Andersen, H. E., Thodsen, H., Olesen, J. E., Børgesen, C. D., Refsgaard, J. C., Sonnenborg, T. O., Karlsson, I. B., Christensen, J. P., Markager, S., and Jeppesen, E.: Effects of changes in land use and climate on aquatic ecosystems: Coupling of models and decomposition of uncertainties, Sci. Total Environ., 657, 627–633, https://doi.org/10.1016/j.scitotenv.2018.12.055, 2019.
van Genuchten, M. T.: A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils, Soil Sci. Soc. Am. J., 44, 892–898, https://doi.org/10.2136/sssaj1980.03615995004400050002x, 1980.
van Vuuren, D. P., Edmonds, J., Kainuma, M., Riahi, K., Thomson, A., Hibbard, K., Hurtt, G. C., Kram, T., Krey, V., Lamarque, J.-F., Masui, T., Meinshausen, M., Nakicenovic, N., Smith, S. J., and Rose, S. K.: The representative concentration pathways: an overview, Climatic Change, 109, 5–31, https://doi.org/10.1007/s10584-011-0148-z, 2011.
Windolf, J., Tornbjerg, H., Hoffmann, C. C., Poulsen, J. R., Blicher-Mathiesen, G., and Kronvang, B.: Successful reduction of diffuse nitrogen emissions at catchment scale: example from the pilot River Odense, Denmark, Water Sci. Technol., 73, 2583–2589, https://doi.org/10.2166/wst.2016.067, 2016.
Wriedt, G. and Rode, M.: Modelling nitrate transport and turnover in a lowland catchment system, J. Hydrol., 328, 157–176, https://doi.org/10.1016/j.jhydrol.2005.12.017, 2006.
Wulff, F., Humborg, C., Andersen, H. E., Blicher-Mathiesen, G., Czajkowski, M., Elofsson, K., Fonnesbech-Wulff, A., Hasler, B., Hong, B., Jansons, V., Mörth, C.-M., Smart, J. C. R., Smedberg, E., Stålnacke, P., Swaney, D. P., Thodsen, H., Was, A., and Żylicz, T.: Reduction of Baltic Sea Nutrient Inputs and Allocation of Abatement Costs Within the Baltic Sea Catchment, Ambio, 43, 11–25, https://doi.org/10.1007/s13280-013-0484-5, 2014.