Table 2 shows the data used in this study for model application with temporal and spatial coverage. The meteorological data, daily precipitation and temperature, were obtained from the UK Met office (Met Office, 2012). Records of continuous daily water discharge at the Hore and the Hafren were obtained from the National River Flow Archive (NRFA). Discharge records after 2009 were obtained from the Plynlimon experimental site in the UK Center for Ecology and Hydrology (CEH). Weekly water quality data from the Plynlimon experimental site were also provided by CEH (Neal et al., 2013, Norris et al., 2017).

The INCA-C model was employed to reproduce and simulate carbon dynamics in the Plynlimon catchment. INCA-C is a dynamic, semi-distributed catchment scale process-based model of DOC concentrations and fluxes in surface waters (Futter et al., 2007). The INCA model was initially developed to simulate nitrogen (Whitehead et al., 1998a,b) and phosphorus (Wade et al., 2002a) at a catchment scale. Several sub-models were later added, including carbon (Futter et al., 2007), sediment transport and soil erosion (Lazar et al., 2010), pathogens (Whitehead et al., 2016), an organic contaminant model (Lu et al., 2016), and microplastics (Nizzetto et al., 2016).

The INCA-C model simulates daily DOC concentrations, flux and water flow in streams. Fluxes and concentrations of DOC are simulated by solving mass balance equations for terrestrial processes while simultaneously solving flow equations. INCA-C requires a daily time series of soil moisture deficit (SMD), hydrological effective rainfall (HER), air temperature and precipitation as inputs. The latter two input datasets were estimated using a semi-distributed rainfall-runoff model, PERSiST (Precipitation, Evapotranspiration and Runoff Simulator for Solute Transport model, Futter et al., 2014). PERSiST simulates water fluxes in a catchment scale, with a series of user-defined buckets indicating different hydrological processes such as snowmelt, surface runoff generation, soil water storage, groundwater and streamflow. It is particularly designed to provide input time series for the INCA models. In addition, spatial data describing catchment land cover and soil properties, the location of point sources, abstraction and effluent discharge are required to run INCA-C. INCA-C then provides daily time series of flow, DOC and DIC (dissolved inorganic carbon) concentrations at each reach boundary, as well as profiles along with descriptive statistics of these variables at selected sites. INCA-C is calibrated using the time series of observed and modeled daily discharge and point estimates of in-stream DOC concentrations to minimize the differences between observed and modeled values.

In this study, climate scenarios generated by the weather@home system (Massey et al., 2015, Guillod et al., 2018) were applied to explore the potential impacts of climate change on hydrology, the carbon in catchments. The modeling system- weather@home (hereafter, w@h) consists of a global climate model (GCM) with a nested regional climate model (RCM) driven by sea surface temperatures as well as other forcings. It enables the production of a large number of ensembles of weather events with the help of volunteers distributed computing across the world. In this model application, weather@home2 (hereafter w@h2, which includes an improved land surface scheme) was used to generate 30-year-long weather time series of rainfall and temperature projections for three different periods. The data set is specifically designed to support a risk-based approach for the study of extreme events such as drought and heavy precipitation, both of which require the spatial consistency of hydro-meteorological data sets (Guillod et al., 2018). The scenarios in each future time slice all follow the Representative Concentration Pathway 8.5 (RCP8.5) and sample the range of sea surface temperatures and sea ice changes from CMIP5 (Coupled Model Intercomparison Project Phase 5) models. A number of daily and (or) monthly output variables are available from the RCM including temperature, precipitation, surface air humidity, mean sea level pressure, and potential evaporation. To assess the worst possible future condition of water quality, the RCP8.5 was selected as it is the most severe scenario among the current available future scenarios. Of particular relevance for water quality modeling is the changes in intense rainfall, which can affect leaching and drought, which influences chemical retention times and river flows.

From a 25 × 25 km grid over Europe, 100-time series are available for each of time slices which represent baseline (BL, 1975–2004), near future (NF, 2020–2049), and far future (FF, 2070–2099) climatic scenarios. Outputs from grid cells include Plynlimon catchment in Wales were averages to generate input time series for water quality modeling application (Figure 3).

The precipitation data was bias-corrected in order to reproduce the patterns and distributions of observed data better, and a linear approach was applied (Guillod et al., 2018). The bias correction method was conducted by using spatially-variable monthly bias-correction factors, obtained by the ratio of the modeled and observed monthly average precipitation. However, very little bias was observed from temperature, so temperature bias was not explicitly corrected. More details about weather@home2 data and bias-correction method can be found in Guillod et al. (2017) and Guillod et al. (2018). Figure 4 shows the comparison of bias-corrected precipitation (purple) and raw precipitation from w@h2 (gray) with observed precipitation. A dry bias during summer months on precipitation was well-corrected.

Deposition scenarios from the Emissions Monitoring and Evaluation Programme (EMEP) were applied for this study to include the impacts of historical, present, and future acidic deposition. Historical sulfate deposition time series from 1970 to 2010 on a 50 × 50 km grid resolution was taken from the output of the EMEP/MSC-W model, as well as the future projection on the depositions until 2030 was also provided under the Gothenburg Protocol Currently Legislated Emissions (CLE) scenario (Posch et al., 2007). The chemical transport model developed at the Meteorological Synthesizing Center–West (MSC-W) is called the EMEP/MSC-W model. After 2030, the deposition is assumed to remain constant until 2090. More information on EMEP data can be found in Schöpp et al. (2003).

Reduction in acidic depositions is believed to have helped acidification recovery of European rivers and lakes (Skjelkvåle et al., 2003). As explained in the previous section, organic matter solubility in the INCA-C model is controlled by the concentrations of strong acid anions in soil solution (Futter et al., 2009, 2011). As soil solution data are not readily available, sulfate deposition data are assumed to be a reasonable surrogate for soil solution strong acid anions concentrations. Deposition for the forest land cover type was applied as an input to the modeling study. A smoothed sulfate deposition time series was created using estimates of average monthly wet sulfate deposition. The linearly interpolated monthly sulfate loadings were multiplied to daily precipitation depth to produce monthly average deposition values. These monthly deposition estimates were then assumed to be representative of the soil solution concentration throughout the soil profile. Estimated deposition showed that the peak was in the 1980s and rapidly declined since then. For the calibration, a single deposition series was provided. For the climate projections, 300 deposition time series were produced according to each input precipitation for each climate realizations for each of w@h2 time slices (baseline, near future and far future).

The Monte Carlo simulation produced an ensemble of 10,000 model runs and those results were assessed based on observed values of flow and water quality (DOC) at reach two for both the Hore and Hafren catchments (Lower Hore and Lower Hafren). The best model parameter set was selected for the Hore and Hafren model and used for the rest of the study. A combination of model evaluation measures was used to determine the best performance model including (a) visual assessment of time series, (b) model performance statistics and (c) total DOC loads estimation, and selected model performance was described.

Figure 5 shows the calibration and validation results of flow and DOC concentrations at Lower Hafren and Lower Hore. It is seen that the models are successfully representing a temporal pattern of DOC concentration and flow. The model simulations did not always capture the magnitude of peaks in DOC concentration. However, the models were able to capture the timing of peaks and seasonal patterns of DOC concentrations.

Table 3 shows the performance indices of the INCA model in Plynlimon catchments (1995–2013). The performances were measured by Nash-Sutcliffe Efficiency (NSE), log NSE for streamflow and percent bias (PBIAS) for DOC concentration. Perfect matching was represented by the value of one for NSE. A negative value of NSE means that the average observed value is a better predictor than the simulated values of model simulation. NSE is one of the most popular indices for evaluating models, but it is sensitive to high or low extreme values. Therefore, log-transformed discharge (log NSE) was evaluated to consider the variability balanced at low values (low flows). The NSE values of daily simulated streamflow in both the Hore and Hafren catchments show acceptable ranges for a whole period of model evaluation (1995–2013), 0.68 and 0.66, respectively. The log NSE values were improved from the values of NSE in this study, which is 0.73 and 0.68, respectively, for the Hore and Hafren catchments. The simulated daily flow provided an acceptable reproduction of the observed flow, including a low flow period, particularly given the flow in Plynlimon is generally small.

Concerning the DOC simulation, PBIAS is <10% bias for both catchments, and INCA-C results can be considered satisfactory. PBIAS measures the average tendency of the modeled values to be larger or smaller than observed counterparts. The optimum value is 0, and negative values indicate a model overestimation bias while positive values mean an underestimation of modeled values. As suggested by Moriasi et al. (2007), satisfactory simulation can be indicated by PBIAS from ±25% up to ±70% for water quality variables. The sulfate deposition data application was helpful to the model calibrations. Model calibrations without the sulfate deposition data were shown −42% and −16% of PBIAS values for the Hore and Hafren, respectively.

The monthly in-stream DOC flux (kg/day) was calculated to demonstrate whether the model reproduced similar seasonal loads. Figure 6 showed the simulated and observed monthly DOC loads, which yielded similar results and well-reproduced seasonal patterns. The precipitation-driven annual DOC fluxes showed a great agreement between observed and simulated DOC fluxes (Figure 7). The degree of fit in terms of fluxes load with an R² of 0.71 (for Hore catchments) and 0.77 (for Hafren) are acceptable.

The many measures used for evaluating model performance showed generally satisfactory simulations, in terms of reproduction of the seasonal patterns of both streamflows and DOC concentrations over a long period, given the uncertainty that characterizes both model results and measurement data values. It suggests the model can be used with confidence in predicting future projections.

Model simulations were forced by a combination of gridded data for sulfate deposition and daily weather inputs and calibrated against point observations of streamflow and water chemistry. The gridded data used here provided complete, gap-free time series, which is a prerequisite for INCA modeling. The EMEP gridded data are widely used when modeling future trajectories of recovery from acidification (e.g., Posch et al., 2007, Futter et al., 2009). Within a region, e.g., the Severn River catchment, long-term temporal variation in sulfate deposition is higher than spatial variation. Thus, it was considered appropriate to use gridded EMEP data as an input to the model. Ledesma and Futter (2017) showed that gridded climate products could give as good or better simulations of streamflow when compared to observational weather data. This finding, combined with the difficulties in obtaining gap-free observational weather data and the need for consistency when applying model parameterizations based on the present-day to future conditions motivated the use of gridded climate data throughout this study. Point observations of streamflow and water quality are also spatially-integrated values as they are the result of processes occurring throughout the catchment.

The criteria for behavior were defined by examining the observed data and identifying a range of values for monthly mean and maximum DOC concentration from 1995 to 2013. Behavior criteria were set to include an increasing trend of DOC concentration and reasonable ranges of monthly mean and maximum DOC concentrations for selected years and months. Summer months (August, July) were chosen and year of 1996, 2002, 2004, and 2008 were all selected to demonstrate increasing trends of DOC concentration. The top 2,000 model runs were retained to proceed with the general sensitivity analysis. Each simulation result consists of the parameter vector itself and the behavioral outcome, i.e., whether the particular parameter vector gave rise to the behavior or not. The criteria used to separate behavioral and non-behavioral parameter sets. Once the parameter set space was partitioned into two groups, the degree of difference between the two groups identified the model parameter influence. If the distribution of non-behavioral values of a model parameter does not separate from the behavioral one, this indicates that the parameter does not have a significant influence on model results. The Kolmogorov-Smirnov two-sample test was used to evaluate model parameter sensitivity. Following refinement of the parameter ranges, the 2,000 simulations for the Hore and the Hafren produce a mixture of behaviors (m) and non-behaviors (n) as expected.

In this study, 15 parameters for each type of land cover were applied to the sensitivity analysis, including parameters related to carbon processing, microclimate, hydraulic, enzymatic latch, and sulfate effects. As three land types (Forest, unimproved grassland, and peatland) were employed in this study, forty-five parameters in total were examined in the whole sensitivity analysis. The parameters and corresponding statistics (d_m,n) that are significantly above the 95% level are ranked in order of importance and listed in Table 4. Eleven out of forty-five parameters are significant at the 95% level or greater in the Hore catchment. Nine for the Hafren catchment are found influential, and the relatively large separations denoted by the K-S test d_m,n explains that these parameters are important for obtaining the correct behavior in a simulation. In this research, nine parameters contribute significantly to simulate an increasing trend of DOC concentrations in Plynlimon. Among the important parameters, the soil temperature response, sulfate effect (b₁) multiplier, and organic layer hydraulic parameters (residence time and retention volume) have a significant impact on both catchments. It is not a surprise that the temperature effect has ranked as the most important parameter to influence model behavior and DOC concentrations. In soil process, sorption and desorption rates of organic carbon are dependent on soil temperature and moisture status, which decide the total mass of DOC and SOC in soil layers and in-stream DOC concentrations. Parameters associated with the effect of sulfate and enzymatic latch revealed to the influential parameters in Plynlimon (Table 4).

According to the results of w@h2 climate scenarios, Plynlimon precipitation is expected to show a slight increase in winter months and a clear decrease in summer months (Figure 3). The trends are gradual from the near future to the far future time period. In the far future, the average increase of December precipitation is about 15% whilst a 43% average decrease is expected in July. The total annual average precipitation of the Plynlimon is expected to decrease by 3% in the near future and 6% in the far future. Projected changes in the seasonal patterns of precipitation are more significant than the total annual average precipitation. The average daily air temperature is projected to increase for all months for both future time periods (Figure 3). The mean temperature changes are from 9.01°C (8.76°C−9.28°C) as the baseline to 10.22°C (9.94°C−10.44°C) for the near future and to 12.28°C (12.02°C−12.48°C) for the far future scenario.

The calibrated INCA Carbon model was run with w@h2 climate scenarios and acidic deposition scenarios for the baseline, near future and far future for both Hore and Hafren catchments. The flows in Hore and Hafren show a clear seasonality, with low flows in summer months and high flows in winter months (Figure 8). The results of w@h2 driven by the INCA model, future climatic change will enhance this seasonal flow patterns. In particular, summer flows will decrease significantly, while changes in winter flows will be mild. This is probably because the projected lower summer rainfall and increased air temperatures in the Plynlimon area contribute to decrease flows, up to an average of 51% in July under the far future scenario. In winter months, projected flows are expected to increase the range of 1–5% in the near future, and 3–15% in the far future compared to the baseline. However, Plynlimon is a small stream with an average flow with 0.234 m³/s and 0.205 m³/s for Hafren and Hore, respectively (average between 1973 and 2009). It is quite difficult to represent in the model simulation, so the range of each future time slices is overlapping with the range of observed flow period. However, the median lines are showing the degree of changes in each time period (Figures 8A,B).

The INCA model results show a substantial increase in DOC concentrations in the near future and the far future from the baseline for both catchments. The DOC concentrations in the far future are stabilized and slightly reduced in the summer months after the near future. Figure 9 shows the DOC monthly average and maximum concentrations driven by the w@h2 climate scenarios. For the Hore catchment, the simulated baseline values slightly overestimate DOC concentrations from April to July compared to the observed values (Figures 9A,C). The rest of the months are represented in a reasonable range. It is worth noting that the time series of baseline condition started in 1975 and DOC monitoring has started since 1983, which may cause a slight discrepancy between baseline and observed values. The Hafren catchment, observed values lie within the range of the simulated values of baseline (Figures 9B,D).

In the near future for both catchments, there is a significant increase in later summer to autumn (August, September, and October). In the forest-dominant catchment at the lower Hore (Figure 9A), there is a slight reduction in DOC concentration the summer months in far future compared to the near future. However, there is still a big surge from the baseline level. Hafren (Peatland-dominant) catchment shows a different trend (Figure 9B). From September to November, there is an even greater increase in DOC in the far future than the near future while the summer months of DOC concentrations are slightly decreased.

Several droughts and associated low flows in summer months and subsequent re-wetting in early autumn in peatlands could destabilize peatland carbon stocks and enhance in-stream DOC production rates. This process is greatly affected by the combination of temperature and water availability in pore water as well. The w@h2 analyses projected less summer precipitation with warmer temperatures, which seems to cause an enzymatic latch effect in peatland soils. In the model simulations, the soil moisture deficit is a surrogate measure of soil dryness, and it is representative of the difference between the water currently in the soil and the water-holding capacity of the soil. The SMD is a key driver of INCA-C results as in-soil all carbon transformation rates are assumed to be soil moisture dependent. Additionally, in the INCA-C model application presented here, SMD is compared to a drought threshold below which an enzymatic latch effect is triggered. The SMD threshold of 31 mm was applied for this study in Plynlimon, and if SMD values were >31 mm depth, the enzymatic latch effect was triggered. The threshold was achieved from the calibration process, and it is consistent value with a measurement from Plynlimon experiment site (Hudson, 1988). Figure 10 shows the changes in the frequency and duration of soil moisture droughts (when SMD > 31 mm) for the two future scenarios that the enzymatic latch effect was triggered. The w@h2 future scenarios suggest a substantial increase in the duration of each soil moisture droughts (Figure 10B) than its frequency. It can be explained that during the drought, organic carbon is not mobilized. However, subsequent re-wetting after a long and severe drought, DOC production rates are accelerated from the peat. When the peatland soils become saturated and are mobilized, they flush DOC into the stream. Therefore post-drought surge of DOC concentrations become visible in autumn months (Figure 9B).

Maximum DOC concentrations show a substantial increase in the wet season, while a small increase in the summer months is projected at around 20% (Figures 9C,D). The projected sulfate deposition continues to decline until 2030 and remain constant after. The projected changes have inter-annual variability, but no intra-annual variability. Therefore, the sulfate decline could affect the stabilization of DOC concentrations in the far future period. However, the changes in seasonal variability of DOC concentrations are most likely driven by climatic changes and enzymatic latch processes during droughts. While this paper has only presented results for the Plynlimon in Wales, enzymatic latch simulations in peatlands can also be extended to understand the other peatlands catchments in the UK.