Simulated future changes of extreme nutrient loads in a mesoscale agricultural watershed in Bavaria / Simulierte zukünftige Änderungen der Extremwerte für Nährstofffrachten in einem mesoskaligen landwirtschaftlichen Einzugsgebiet in Bayern

Figure 1 shows the location of the upper part of the Altmühl watershed in Germany, which encompasses a total basin area of 980 km² and includes 130 km of river length from the source of the Altmühl to the gauge in Treuchtlingen (10^o 54'48.91”E, 48^o 57'11.31”N). The Treuchtlingen gauge is the outlet in this study. The elevation ranges from 660 m to 406 m at the outlet. The soils in the upper part of the study area are mainly loamy clay and loamy gravelly sand, with pockets of gravelly sand. Below the Altmühl Lake are mainly clayey silt soils along the floodplain. Clay loam soils are predominant near the outlet, with a few areas of Karst in the most south-western tip of the basin. The watershed consists primarily of agricultural (60%), forested (30%), and urban areas (5%). Detailed field data available from the Bavarian State Office for Agriculture, show the agricultural land in 2008 to be mainly winter and summer cereals (26,065 ha), pasture (24,620 ha), maize (7942 ha), rape (2690 ha), alfalfa and clover (1369 ha) with pockets of other vegetable and fruit crops (Table 1). The Altmühl watershed receives approximately 700 mm of precipitation annually. Evapotranspiration is 475 mm and runoff in the basin is 175 mm (StLfW, 1996), the remaining 50 mm is most likely lost to deep aquifer recharge due to the Karst areas in the southern watershed.

The hydrological model Soil and Water Assessment Tool (SWAT; Arnold et al., 1998) was applied to the watershed to examine streamflow, as well as NO₃^--N and TP loads for both a historic and a future period. SWAT was developed by the United States Department of Agriculture to reflect the impacts of agricultural management practices on streamflow, agricultural nutrient, and sediment yields in surface waters. SWAT is a semi-distributed, process-based hydrological model run on a daily time step (Gassman et al., 2007; Arnold et al., 2012). SWAT uses the concept of hydrological response units (HRUs) to model the hydrological processes in the watershed. HRUs are areas of similar land use, soil, and slope in each sub-basin.

ArcSWAT version 510 was run on an ArcGIS 9.3.1 (ESRI 2009, California, USA) platform. The setup for the water-shed was based on a 50 m Digital Elevation Model.

In SWAT, the simulated processes include infiltration, surface runoff, evaporation, plant water uptake, lateral flow, and percolation to shallow and deep aquifers. Flow, sediment, and nutrients are summed across the HRUs in a sub-basin, and the flows and pollutant loads are then routed through channels, ponds, and reservoirs to the watershed outlet. The volume of surface runoff is estimated using the modified SCS curve number (CN) method (USDA, 1972), which takes into account precipitation, abstractions, and soil storage; the lower the CN, the more permeable the surface. The CN is adjusted at each time step, depending on how much soil moisture is available. The potential evapotranspiration was estimated using the Penman–Monteith method.

The crop growth was modeled with the EPIC sub-model (Williams et al., 1984) that bases the phenological development of the plant on the accumulated plant heat units (PHUs), which are a function of the minimum and maximum air temperatures.

The SWAT model has three major forms of nitrogen that it models in mineral soils: (1) the organic pool associated with humus, (2) a plant-available pool in the soil solution, and (3) an insoluble inorganic component (Neitsch et al., 2011). There are five main pools associated with these nitrogen forms: two inorganic pools (NH₄⁺ and NO₃-) and three organic pools (fresh plant residue, stable humic substances, and active humic substances). For details of the transformation processes, see Neitsch et al. (2011).

In SWAT, the P cycle contains three major sources of P in mineral soils: the organic pool associated with humus, a plant-available pool in the soil solution, and an insoluble inorganic component.

SWAT models six different pools of P in the soil; three pools are associated with the inorganic forms of P (solution, active, and stable) and the other three with the organic P forms (fresh, stable, and active). The fresh organic pool is associated with the crop residues and microbial biomass. The active and stable organic P pools are both associated with the soil humus, which is partitioned into these two pools to allow for the P to transform from humic substances to mineralized substances. The P mineralization algorithms are net mineralization algorithms, which take into account immobilization.

The inorganic P in solution is the available form of P that plants can take up. This pool is in rapid equilibrium with the active pool. The primary movement of soluble P in the soil is by diffusion due to a concentration gradient in the soil. Organic and inorganic P may be transported through the attachment to soil particles. The details of the processes are provided in Neitsch et al. (2011).

A suite of RCM data generated by the Ouranos Consortium on Regional Climatology and Adaptation to Climate Change was available for this research (Velázquez et al., 2013). Each simulation from the RCM was driven by a coupled GCM for the time periods 1970–2000 and 2041–2070. In total, seven coherent sets of climate variables of temperature, precipitation, relative humidity, solar radiation, and wind speed were available to drive the SWAT hydrological model (Figure 2). The regional climate models were: Rossby Centre Regional Atmospheric Climate Model (RCA), Regional Atmospheric Climate Model (RAC), Canadian Regional Climate Model (CRC). The driving Global climate models were: Bergen Climate Model (BCM), ECHAM (ECM) version 5 (members 1, 2 and 3), and Hadley Centre Coupled Model (HCM), and Canadian General Circulation Model (CGC).

The global climate models had the A1B greenhouse gas emission SRES projection (Nakicenovic et al., 2000), except for the CGC, which used the A2 scenario. Both of these SRES represent the higher greenhouse gas contribution scenarios. The A1B scenario has CO₂ emissions peaking around 2050, at 16 GtC and fall to about 13 GtC by 2100. In the A2 scenario, the global CO₂ emissions reach 29 GtC by 2100. The 1990 global CO₂ emissions were 6 GtC.

Temperature for each member of the ensemble was bias-corrected using a monthly correction factor based on the difference between the ensemble mean of the 30-year mean monthly minimum and maximum air temperature and the 30-year monthly means of the daily-observed minimum and maximum air temperatures. A bias-correction method for precipitation was applied to all climate simulations using the Local Intensity Scaling (Schmidli et al., 2006) at a sub-daily time step. Finally, the RCM outputs were scaled to a finer resolution of a 1 km grid with the scaling tool SCALMET (Marke, 2008), which preserves energy and mass at the scale of the RCM grid. For more detailed explanations of the climate simulations post-processing, see Muerth et al. (2013).

A semi-automated inverse modeling procedure called Sequential Uncertainty Fitting algorithm (SUFI-2) version 4.3.2 (Abbaspour, 2011) was used for calibrating the SWAT-simulated outputs to the available time series data of streamflow, NO₃^--N, and TP loads. SUFI-2 is a stochastic procedure, drawing independent parameter sets using Latin Hypercube Sampling (LHS). For a detailed description of SUFI-2 procedures, see Abbaspour et al. (2004). The SWAT model was calibrated sequentially, at the yearly, monthly, and daily time steps, for streamflow, NO₃^--N, and TP (Arnold et al., 2012). SWAT was first calibrated (1964–1974) at the outlet gauge (Treuchtlingen) for surface flow (validated from 1975–1984). Subsequently, NO₃^--N, and TP were calibrated (1982–1983) at the monthly time step at the Thann gauge (and validated in 1984). The SWAT model was set up to meet the satisfactory objective criteria listed in Moriasi et al. (2007), whereby the Nash-Sutcliffe efficiency (NS) was used as the main objective function. The details of the calibration are presented in Mehdi et al. (2015).

The available observed station data and simulated climate data were used from 1970–2000 (“reference period”), and 2041–2070 (“future period”), respectively. However, a 5-year warm up period was removed from all SWAT simulations; therefore, outputs were examined and presented from 1975–2000 and 2046–2070, respectively.

All climate simulations were bias-corrected, which allowed a direct comparison to be made between the SWAT simulations using the observed climate and using the future climate simulations. However, to remove the errors brought about by the SWAT model, as well as by the climate models, comparisons were made from SWAT outcomes using the climate simulations from the reference period and SWAT outcomes using the future climate simulations.

Once the SWAT model was set up, it was run with input from the observed station climate data from 1970–2000. Next, all else remaining the same, SWAT was run in turn with each of the climate simulations from the reference period, by replacing the observed climate. The SWAT model was subsequently run with each of the climate simulations for the future period. The output variables examined were always streamflow, NO₃^--N, and TP.

Independent t-tests were carried out (p<0.05) between the output variables (by season) simulated from the reference and the future period.

The climate by the end of the century in Germany is simulated to be warmer, and there will be a noticeable increase in extremes in Europe, such as heat waves, droughts, and precipitation events (Kovats et al., 2014). The climate simulations used in this study showed increases in the mean monthly temperatures ranging from 0.5 to 4^oC with most of the models agreeing on the greatest warming occurring in winter and in summer months (Figure 2a).

The simulated precipitation for most months was higher than in the reference period, especially during December to April, by up to 50%. The months of August and September showed a decrease in precipitation, of up to 80% in August (Figure 2b). Overall, the average annual precipitation in the observed reference period was 705.3 mm, whereas in the future period, the mean of the climate simulations showed this amount to be 752.3 mm. The month with the greatest increase in precipitation for the future was April with an additional 10 mm (average of seven climate simulations).

After the calibration and validation steps, SWAT was able to reproduce the timing of daily dry spells and peak flows. The magnitude of flow events was also modeled satisfactory, but SWAT tended to underestimate the overall flow (PBIAS = 13.8%), especially for the lower streamflow between runoff events. The calibration and validation of SWAT at the daily time step had NS values of 0.57 and 0.68, respectively, whereas the monthly NS values were 0.77 and 0.75, respectively. The NO₃^--N and TP simulations by SWAT corresponded well to the observed timing of the events. However, modeled NO₃^--N was overestimated by SWAT (PBIAS – 11.8%) and TP had overall lower simulated values (PBIAS 33.5%). The NS monthly calibration and validation values for NO₃^--N were 0.77 and 0.72, respectively. For TP, the NS values were 0.47 and 0.52 for the calibration and validation, respectively. Based on the statistics and on the performance of the calibrated model, SWAT was able to satisfactorily simulate stream-flow, NO₃^--N, and TP in the watershed.

The SWAT-simulated streamflow using each of the climate simulations for the future period was compared statistically (p<0.05) to the simulated streamflow using the climate simulations from the reference period (Table 2). The mean streamflow in September was simulated to be significantly lower in the future due to lower future precipitation in August, whereas in spring (April, May, and June), the mean streamflow was significantly higher caused by the greater amounts of precipitation throughout the winter and into the early spring months (Figure 3).

The long-term (1941–2008) average annual flow at the Treuchtlingen gauge was measured to be 5.73 m³/s, with daily flows ranging from a minimum 0.28 m³/s to a maximum 183 m³/s (BLfU, 2012). For the mid-term future, the simulated mean annual streamflow (5.44 m³/s) lies slightly under the historic average flow.

The monthly 10^th percentile flows remain unchanged in the future. However, the simulated mean 90^th percentile flows (statistically representing a 1 in 10-year event) are significantly lower in autumn (Table 3).

Higher mean values for the 90^th percentile flows were simulated during the spring months (Tables 2 and 3).

Comparing the nutrient percentiles in each season between the simulated reference and future climate periods provides an indication of the changes due to possible future climates (Figures 4a and 4b). The simulated nutrient percentiles using the reference and future climates in SWAT were compared with the observed data from the reference time period (mean monthly values from 1975–2000) in the form of boxplots to indicate the deviation by the simulations from the observed data. Ideally, the median of the boxplot from using the reference climate simulations should lie on, or very near, the zero line to indicate no difference between the observed climate and the simulated reference climate outputs. Discrepancies in the simulated observed and simulated reference results are manifold and are due to the imperfect SWAT model setup, the chosen parameter values during the calibration, and due to the observed data. As well, in the climate models, the differences in parameters and in the representation of physical processes and different initial conditions from the members of the same model (i.e., ECM members 1, 2, and 3) will contribute to the spread of the simulated output results.

The SWAT model was able to capture the observed 90^th percentile data well (depicted by the median nutrients in three out of the four seasons lying close to the zero line in Figure 4b), thus indicating greater confidence in the trend toward significantly higher 90^th percentiles of NO_3--N loads in the future in the Altmühl. The future mean monthly NO₃^--N loads are also simulated to be significantly higher in six months of the year (Table 4).

In SWAT, the governing processes simulating NO₃^--N loads transported from the field to the stream are dependent on soil mineralization and decomposition of the fresh organic pools and humus. Net mineralization occurs when the C:N is <20. Immobilization causes a net decrease in NO₃^--N and occurs when the C:N is >30. All N mineralization algorithms in SWAT are net mineralization algorithms that consider immobilization into the equations. Net mineralization is dependent on the soil water availability in relation to the field capacity, and on soil temperature. When soil temperature (>0^oC) and water content increase, more nitrogen is mineralized and added to the soil nitrate pool.

In SWAT, the inputs from fresh organic pools (fertilizer or residue applications) are important sources of nitrogen; however, these pools were not altered during any of the model simulations. The average fertilizer applied to the cropland in the watershed was 67 kg N/ha/yr and 26 kg P/ha/yr, for both the reference and the future period, respectively.

Due to its negative charge, the NO₃- molecule is not sorbed by the soil particles, rather it is soluble and can be carried over land in surface runoff, and it is readily transported by infiltration and lateral flow into the surface layers of the soil, as well as it can move upward in the soil in response to evaporation. April was the month with the highest increase in precipitation, yet May was the month when the difference in mean NO₃^--N load was the largest (0.41 kg/ha) between the reference and the future period (Table 4); this is explained by the lag time for NO₃^--N to be transported by the infiltrated water through the soil and into the river.

Denitrification removes NO₃^--N from the soil under anaerobic conditions. In SWAT, when the water-filled pore space is >60%, denitrification occurs. This process is also temperature-dependent, so higher soil temperatures (>0^oC) coupled with anaerobic conditions enhance denitrification. However, denitrification did not cause any NO₃^--N changes in the future simulations.

Plant growth also removes N from the soil during the growing season as a function of the growth stage. Due to the warmer growing season in the future climate, some crops (e.g., maize) extracted more nutrients from the soil to meet the growth needs, but this did not affect the 10^th or 90^th percentiles. Phosphorus is much less mobile than nitrate. In SWAT, the P forms can be in solution and plant-available or they can be insoluble and bound to soil particles or precipitated out of solution.

Due to the low mobility of P, SWAT accounts for P in surface runoff and leachate only from the top 10 mm into the first soil layer. Organic and mineral P attached to the soil particles can be transported by surface runoff to the stream. This amount is associated with sediment loading.

For the mid-term future, the TP loads simulated during winter were diminished; this was depicted by the decreasing relative changes in the 10^th and 90^th percentiles, and also by the mean TP loads, which were all simulated to be significantly lower from December to February (Table 5 and Figure 5). April and May mean TP loads and the spring 10^th percentile were higher in the future (Figure 5 and Table 5), showing a shift toward fewer low events. The 10^th percentile TP loads were well simulated, and therefore, indicate confidence in this trend.

TP loads in the reach were correlated (R² = 0.74) to streamflow, so that the changes in TP seasonality followed the streamflow pattern closely. The months with the largest increases in flow compared with the reference period were April to June. The lowest flows occurred from August to February (Figure 3). The 90^th percentiles of TP loads decreased significantly in winter and otherwise remained statistically similar to the reference period.

Decomposition and mineralization of fresh organic residues and of the humus add plant-available P to the soil. These processes are controlled by the decay rate constant, which is determined in part by the soil temperature; warmer temperatures increase the decay rate.

For the future period, each month showed that the mean monthly temperature increased compared with the reference period. Warmer temperatures contribute to an accelerated mineralization of nutrients; this happened especially during the winter months when the soil temperature remained above 0^oC more often.

By the mid-term future, at the outlet of the Altmühl water-shed with a size of 99,335 ha, the simulated mean monthly median NO₃^--N loads were 1.0 to 26.8 Mg higher than from 1975–2000; and simulated mean monthly median TP loads were up to 3.1 Mg lower.

The 90^th percentile increase in the future streamflow in spring did not consistently increase in the 90^th percentiles of the nutrient simulations. The lower 10^th percentile in autumn for flow also had no impact on the nutrients 10^th percentiles. Very high nutrient loads are usually due to a combination of high precipitation events with a field management practice, such as fertilization. Very low nutrient loads can be due to low precipitation, combined with rapid plant growth. The dynamics of TP and NO₃^--N transport, are therefore, better described using smaller time steps. Nitrate can be transported after a precipitation event through several pathways (infiltration, lateral flow, surface runoff) and end up within hours or days in the stream, whereas TP is mainly transported by surface runoff and interacts with the top 10 mm of the soil. The transported TP is sorbed to sediments and is transported by erosion or surface runoff processes, which take longer to travel because the particulate P can be deposited and eventually never make it to the stream. TP is more correlated to streamflow and high surface runoff events than NO₃^--N loads are.