Impact of Climate Change on the Hydrology of the Forested Watershed That Drains to Lake Erken in Sweden: An Analysis Using SWAT+ and CMIP6 Scenarios

Abstract

Precipitation and temperature around the world are expected to be altered by climate change. This will cause regional alterations to the hydrological cycle. For proper water management, anticipating these changes is necessary. In this study, the basin of Lake Erken (Sweden) was simulated with the recently released software SWAT+ to study such alterations in a short (2026–2050), medium (2051–2075) and long (2076–2100) period, under two different climate change scenarios (SSP2-45 and SSP5-85). Seven global climate models from the latest projections of future climates that are available (CIMP 6) were compared and ensembled. A bias-correction of the models’ data was performed with five different methods to select the most appropriate one. Results showed that the temperature is expected to increase in the future from 2 to 4 °C, and precipitation from 6% to 20%, depending on the scenario. As a result, water discharge would also increase by about 18% in the best-case scenario and by 50% in the worst-case scenario, and the surface runoff would increase between 5% and 30%. The floods and torrential precipitations would also increase in the basin. This trend could lead to soil impoverishment and reduced water availability in the basin, which could damage the watershed’s forests. In addition, rising temperatures would result in a 65% reduction in the snow water equivalent at best and 92% at worst.

1. Introduction

Climate change is a worldwide known phenomenon that is widely accepted by the scientific community [1]. This phenomenon is expected to alter precipitation and temperatures around the globe in various ways, making climates drier in some places and wetter in others. One of the main impacts will be on the hydrological cycle and water bodies on land, especially in streamflow regimes [2,3]. Rivers, lakes, aquifers, and their basins are strongly affected by changes in evapotranspiration, runoff processes, floods, and erosion [4,5,6,7]. Hence, changes in the temperature and precipitation patterns will alter these processes. Changes that affect snow processes—such as the snow accumulation pattern or glacier melt runoff—would also have a high impact [8,9]. All these variations will incur long-term negative consequences for ecosystems and human activities, especially in drainage basins. Serious detrimental alterations are currently occurring through droughts and floods [10]. The ecosystems that are most vulnerable to such alterations include those formed by long-living water-dependent vegetation, such as forests. Overall, water resource management could be seriously compromised by climate change. Water-dependent ecosystems would be particularly affected. The only way to prepare for the likely future scenario is to foresee the effects of climate change on the hydrological cycle and make decisions accordingly. In recent years, multiple works have addressed these problems in different parts of the world [11,12,13,14,15,16,17].

Hydrologic circumstances vary by location. Therefore, the effects of climate change on local hydrological and ecosystem processes are expected to differ even under the same climate predictions. For example, in northern Europe, precipitation is expected to increase. In flat, humid areas with vegetation, this increase in precipitation may be linked to flooding and bogging down of the land. However, in a dry, sloping area, the increase in precipitation would lead to a very serious increase in soil erosion. In both cases, palliative measures to adapt to the increase in temperature must be different. When planning effective regional water resources, it is important to understand the regional conditions and to develop region-specific assessments of climate change [18].

In the case of Sweden, as in the rest of northern Europe, the average amount of precipitation is expected to increase in the future [19,20]. However, this does not necessarily imply an increase in water availability. On the contrary, the increase in seasonal precipitation is likely to be related to extreme events in terms of the daily rainfall [19,20]. This scenario could paradoxically lead to a reduction in water availability and to soil impoverishment [17]. These phenomena could seriously alter the forest ecosystems that depend on water resources in specific areas and the soil quality.

In the south of Sweden, as in any land areas towards the poles, snow is an important component of the climate and the hydrologic cycle [21]. All basins are partially covered by snow during the cold season. It should be noted that in snow-dominated basins, snow processes are of vital importance. Snowfall, snow accumulation, and snowmelt have a strong effect on the hydrological cycle [22]. Changes in the estimation of these water balance components can alter a hydrological simulation significantly [23]. Therefore, researchers should consider these components when modelling snow-covered basins [24].

The snow water equivalent (SWE) and snowmelt are the most essential snow statistics for hydrologists [25]. The SWE is described as a water column that forms as a result of the melting of unit cross-section snow samples with a height equal to the depth of the snowpack at the measurement point [26]. The SWE is the major parameter determining the magnitude of the snowmelt runoff volume. It is used as a variable in snow runoff analysis to estimate the distribution and quantity of snow [27].

The aim of this study is to model the Lake Erken basin and use the model to predict how climate change is likely to affect the hydrological cycle in this area and their ecosystems. To achieve this aim, the Soil and Water Assessment Tool+ (SWAT+) [28] software was used to create the model. Although its predecessor, SWAT [29] has been widely used around the world, and to date, few published works have used SWAT+. Only three studies have used SWAT+ to analyse the impact of climate change on a watershed. Chawanda et al. [30] used mass balance calibration and reservoir representations to evaluate, at the regional scale, the climate change effects in Southern Africa. Senent-Aparicio et al. [31] analysed the impact of climate change on environmental flows in the northwest of Spain using a new post-processing tool for SWAT+ models. Kiprotich et al. [32] studied the surface runoff response to climate change and land-use change in Nairobi, Kenya. None of these studies used global climate models from the Coupled Model Intercomparison Project 6 [33] but instead used its predecessor, the Coupled Model Intercomparison Project 5.

Even though SWAT was not developed for flood modelling, there have been SWAT flood studies documented [34,35,36,37,38]. Tan et al. offer a review of some of these studies [39], highlighting the problems that SWAT presents when performing calibrations and validations based on extreme events and the need for improvements in SWAT to capture extreme events. Despite this, most extreme performance evaluation studies reviewed found satisfactory results, with a focus on peak flow comparisons. No studies using SWAT+ to simulate extreme events have yet been published.

The most common data used to calibrate and validate SWAT models is the streamflow in the outpoint of a basin. Given the importance of the snow and snow parameters in the Lake Erken basin, in this study, the SWE data have also been used in the calibration and validation of the model, in a complementary way. Finally, the temperature and precipitation data obtained from various climate change scenarios will be used to predict the evolution of the Lake Erken watershed in the coming decades.

2. Methodology

2.1. Conceptual Model

First, a SWAT+ model was developed, then calibrated and validated with observed streamflow and SWE satellite data. In addition, seven global climate models (GCMs) were weighted using real historical climate data as the reference. The future data were used under two climate change scenarios to simulate the water cycle of Lake Erken during a short-term period, a medium-term period and a long-term period. Finally, various parameters related to floods were calculated with Indicators of Hydrological Alteration in Rivers (IAHRIS) [40]. The aim of this step was to compare the historical floods in Lake Erken with the projected floods that would be altered because of the climate change effects. The methodology of this article is represented in Figure 1.

2.2. Description of the Study Area

Lake Erken is located in the eastern part of Sweden (59°50′37″ N, 18°35′38″ E) at an altitude of 10 m above sea level (Figure 2). The surface area of the lake is approximately 24 km². It is a shallow lake; the mean depth is around 9 m, and the deepest point is 21 m underwater. The lake has a residence time of 7.4 years. It can be described as a moderately eutrophic lake with an intermediate level of productivity. Its surface is usually ice-covered during winter, whereas during the summer, the water is stratified [41,42]. Both of these phenomena have been largely studied in Lake Erken [43,44,45].

The Erken basin is oriented toward the east and does not present significant slopes. It is a relatively small drainage basin (141 km²) that is covered by forest, without any significant anthropic activity [46]. In the basin, there are exclusive areas of both deciduous forest and evergreen forest; in most of the basin, both types of vegetation share the space, forming a mixed forest that dominates the watershed (Figure 2c). This forest plays a vital role in the water cycle. Forests regulate processes such as evapotranspiration, runoff, and water retention in forested basins [47,48,49,50,51]. The main large water inflow point to the lake is Kristineholm (Figure 2b). This input is where the water discharge for calibration and validation is measured.

In the south of Sweden, the climate is humid and continental, with a warm summer (17.3 °C average in July and August) and a weak winter—despite the country’s high latitude. The length of daylight varies from 18 h in June to 6 h in December. According to the historical data used in this study for 1990–2014, in the Erken basin the average temperature in July and August was 17.3 °C and in February it was −4 °C. Precipitation does not change much during the year but is slightly higher during autumn; the annual average is 519 mm. Precipitation in the form of snowfall occurs between December and March, and the basin is covered by snow for 75 to 100 days a year.

2.3. SWAT+ Model

The Soil and Water Assessment Tool (SWAT) [29] software provides a hydrological model that is used worldwide. It has been successfully used in a range of scenarios with different climatic conditions, land management practices, and temporal and spatial scales. During the last 20 years, SWAT has been implemented periodically to meet the diverse requirements of the scientific community around the world. However, the current framework has reached the limit of its potential development. The most recent version, SWAT+, improves the runoff routing capabilities while preserving the model’s computational efficiency and ease of use [28].

Water dynamics are represented by fluctuations in the hydrologic response units (HRUs) in both the SWAT and SWAT+ models. Each HRU is a unique combination of land-use, slope, soil, and management activities, which are connected by a geographic information system (GIS) interface. Using this GIS interface, in both models, the modelled basin is divided into various sub-basins, which are further sub-divided into HRUs [52]. The concept of water balance, represented by Equation (1), is applied in the model as the watershed’s primary driver of all hydrology.

$$S{W}_t=S{W}_o+\sum \left(V_i-Q_i-E_i-P_i-Q{R}_i\right)\times ∆t$$

(1)

where SW_t and SW_o represent the final and initial soil water content (mm/day); V_i represents the precipitation (mm/day); Q_i represents the surface runoff (mm/day); E_i represents the evapotranspiration (mm/day); P_i represents the percolation (mm/day); QR_i represents the return flow (mm/day); and ∆t represents the time interval (day). The i term refers to the index.

Despite using the same equations, in SWAT+ the elements of the watershed—such as aquifers, land-use units, HRUs, ponds, and reservoirs—are defined as spatial objects. This feature enhances the flexibility of the configuration and discretisation of the basin compared with the earlier SWAT [28].

2.4. Model Setup

The use of a hydrological model with SWAT+ requires specific geographical information about the area of interest. The data required include a digital elevation map (DEM), a land cover map, and a soil map. The model also requires meteorological input data. For this study, the DEM was obtained from the Shuttle Radar Topography Mission (SRTM). The SRTM uses a single-pass space-borne interferometric SAR system, which operates in both the C-band (5.6-cm) and the X-band (3-cm) frequencies to collect data about the earth’s surface elevation [53]. The land cover map was obtained from Glob Cover 2015, which provides a 300-m resolution [54]. The soil data were gathered from the Food and Agriculture Organisation of the United Nations’ Harmonised World Soil Date. This collection contains information for 16,000 map units with two soil layers, namely 0–30 and 30–100 cm deep [55].

For the meteorological input data, daily precipitation and temperature data were obtained from the Erken Laboratory meteorological station (59°51′30.7080″ N, 18°24′17.5536″ E), situated on the Malma islet (Figure 2b). This station provides automated measurements daily temperature and precipitation data and other information. The Hargreaves method was used to determine potential evapotranspiration [56]; this method requires only the precipitation and temperature data. According to Oudin et al. [57], hydrological models that employ parsimonious temperature-based methods perform similarly to models that use more data-demanding methods.

The calibration and validation of the simulated streamflow required discharge data. The discharge data available has been measured daily at Kristineholm, the largest input of Lake Erken, since mid-2006. In this study, the SWE was also calibrated. The SWE data from the Copernicus Global Land Service were used. It offers a 0.05° spatial resolution, calculated by integrating passive microwave radiometer brightness temperature readings from the Special Sensor Microwave Imager/Sounder with synoptic weather station network snow-depth data [58].

2.5. Calibration and Validation of the SWAT Model

To locate the most influential parameters for the streamflow, we performed a sensitivity analysis. Then, the sensitive parameters were adjusted with a daily automatic calibration for the 2007–2015 period. Both analyses were performed in Toolbox [59], a free software designed to perform SWAT+ model sensitivity analysis, calibrations, and more. Afterwards, we used the 2016–2020 daily data to validate the SWAT+ model.

For the sensitivity analysis, Toolbox uses the Sobol method [60]. Within an ensemble, it divides the overall output variance into the variation produced by each parameter. For automatic calibration, Toolbox uses a dynamically dimensioned search (DDS) [60]. The Nash–Sutcliffe efficiency coefficient was applied as the target function in this investigation. The goodness of fit was determined using Moriasi’s four recommended statistics. In addition, for both calibration and validation, the performance was evaluated using Moriasi’s four recommended statistics [61]. These are the coefficient of determination (R²) (Equation (2)), the Nas–Sutcliffe efficiency (NSE) (Equation (3)), the standard deviation of measured data (RSR) (Equation (4)), and the per cent bias (PBIAS) (Equation (5)). The criteria proposed by Moriasi [61] were applied. Both calibration and validation have been performed on a daily scale.

$${\mathrm{R}}^2={\left[\frac{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}\left({\mathrm{O}}_{\mathrm{i}}-\overline{\mathrm{O}}\right) \left({\mathrm{S}}_{\mathrm{i}}-\overline{\mathrm{S}}\right)}{\sqrt{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{O}}_{\mathrm{i}}-\overline{\mathrm{O}}\right)}^2} \sqrt{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{S}}_{\mathrm{i}}-\overline{\mathrm{S}}\right)}^2}}\right]}^2$$

(2)

$$\mathrm{NSE}=1-\frac{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{O}}_{\mathrm{i}}-{\mathrm{S}}_{\mathrm{i}}\right)}^2}{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{O}}_{\mathrm{i}}-\overline{\mathrm{O}}\right)}^2}$$

(3)

$$\mathrm{RSR}=\frac{\sqrt{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{O}}_{\mathrm{i}}-{\mathrm{S}}_{\mathrm{i}}\right)}^2}}{\sqrt{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{O}}_{\mathrm{i}}-\overline{\mathrm{S}}\right)}^2}}$$

(4)

$$\mathrm{PBIAS} =\frac{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}\left({\mathrm{O}}_{\mathrm{i}}-{\mathrm{S}}_{\mathrm{i}}\right)}{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{O}}_{\mathrm{i}}}\times 100$$

(5)

In Equations (2)–(5), the observed and simulated waterflow data are ${\mathrm{O}}_{\mathrm{i}}$ and ${\mathrm{S}}_{\mathrm{i}}$; the average observed and simulated water flow values are $\overline{\mathrm{O}}$ and $\overline{\mathrm{S}};$ and n refers to the total dataset.

To enhance the model, we performed a manual calibration of the SWE with the snow-related parameters. Since Toolbox does not allow for calibrating the SWE, we performed the calibration manually in SWAT+ Editor, which is a tool for editing, running, and saving changes to a SWAT model. Then, we validated each test in a spreadsheet by applying Equations (2)–(5) to the SWE.

2.6. Global Climate Models and Climate Change Scenarios

To make predictions for the Erken watershed, future climate models that simulate the temperature and precipitation in the coming years are required [62]. Historical data from seven GCMs were downloaded from the Coupled Model Intercomparison Project (CIMP6) website (https://esgf-node.llnl.gov/search/cmip6/, accessed on 26 November 2021) [33]. The downloaded data were compared with the real historical data for precipitation and temperature in the Erken watershed. Based on the analysis of McSweeney et al. [63], seven GCMs whose performance for Europe had been identified as satisfactory were selected: BBC, CanESM5, EC-Earth-Veg, GFDL, INM-CM5, MiroC6, and MRI. These seven GCMs were compared and assessed in this study. We used the methodology proposed by Pulido-Velazquez et al. [55]. The Id indicator (Equation (6)) was calculated by summing the increases in the mean and variance of the control series over the historical series for the period 1985–2014, as follows:

$$\begin{array}{c}I{d}_i={\displaystyle \sum _{n,m=1}^2}I{d}_i\left(V_n{S}_m\right) ;\\ I{d}_i\left(V_n{S}_m\right)=\sum _{j=1}^{12}\frac{{\left(V_n{S}_m\right)}_i^j-{\left(V_n{S}_m\right)}_{Hist}^j}{{\left(V_n{S}_m\right)}_{Hist}^j}\end{array}$$

(6)

where subscript I is a subindex for a particular GCM, V1 is the rainfall variable, V2 is the temperature variable, S1 is the mean monthly value, S2 is the monthly standard deviation, and superscript j is the number of months in a year. This indicator allows a ranking of GCMs to be established in terms of goodness of fit to the observed time series, so that an ensemble of predictions can be proposed, which gives more weight to the best scoring (lower Id index values). Then, the total Id indicator was rescaled by defining the Id* index, so that the sum of the values obtained for all models is equal to 1. Finally, the Ib indices were obtained, which are the complementary values of the Id* indicator, i.e., (1-Id*), rescaled to 1. These Ib indicators were used as weights applied to the series obtained with each model to create an ensemble of predictions.

The historical data of each GCM were bias-corrected using the five statistical transformations available in the “qmap” package for R software. These were the robust empirical quantiles (RQUANT), distribution derived transformations (DIST), empirical quantiles (QUANT), parametric transformations (PTF), and smoothing spline (SSPLIN) [64]. These methods try to align the distribution of simulated data with that of observed climate data and are commonly used in hydrological and climatic studies [65,66,67].

The performance of the various methods was measured using the mean absolute error (MAE) between the observed and corrected data. This procedure is recommended by Gudmundsson et al. [64]. The method that best fitted the observed historical data was selected for the correction of future climate change scenarios.

Lastly, the climate projection was analysed under two different shared socioeconomic pathway (SSP) scenarios [68]. These were the SSP 2-45 and the SSP 5-85 scenarios. In the SSP 2, global emissions are predicted to follow current patterns. This implies significant obstacles for reduction and adaptation, but neither is particularly severe. On the other hand, the SSP 5 illustrates a scenario in which economic development takes precedence over environmental impacts. As a result, the challenges posed by climate change are difficult to meet [69]. The SSP 2-45 and SSP 5-85 are, respectively, the equivalents of the RCP4.5 and RCP8.5 of the CMIP 5, updated with socioeconomic reasons [70]. These two scenarios have been widely used in climate change studies because they allow a comparison between a more positive outlook (SSP 2-45) in which greenhouse gas emissions are intermediate and the effects of climate change are not as severe, and a more extreme outlook (SSP5-85) in which the challenges of climate change are greater.

2.7. IAHRIS Software

Indicators of Hydrologic Alteration (IHA) are commonly used to compare the variations in hydrological regimes between a base scenario and an altered scenario [71]. IAHRIS is a software program created by the Centre for Public Works Studies and Experimentation (CEDEX) [40]. It calculates various hydrological parameters of the flow and uses them to calculate various IHA. IAHRIS is particularly suited for evaluating flood events, since the indicators are divided into ordinary values, maximum extreme values (floods), and minimum extreme values (droughts). In addition, there are recent examples of previous studies that use IAHRIS and SWAT with satisfactory results [72,73]. In this work, IAHRIS was used to compare historical floods in Lake Erken with the projected floods that were altered in the future because of climate change effects. To study extreme peak flow events (floods) IAHRIS offers 8 parameters (Qc, Ql, Qconnec, Q5%, CV (Qc), CV (Q5%), consecutive days of flooding and flood days per month), which are described as follows. To analyse the magnitude and frequency of the floods, we used:

-

Qc: the mean of the maximum annual daily flow.

-

Ql: bed generation flow; this parameter represents the flow that performs most of the work of material relocation and is responsible for the geomorphology of the channel.

-

Qconnec: the maximum flow that ensures the river channel–floodplain connection is represented by the connectivity flow.

-

Q5%: this parameter is the flow corresponding to the average flow curve classified at the 5% exceedance percentile.

To analyse the variability of the floods, we used:

-

CV (Qc): coefficient of variation of the series of maximum annual daily flow rates.

-

CV (Q5%): coefficient of variation of the series of usual floods.

To study the duration of the floods, we used the maximum number of consecutive days with an average daily flow rate greater than the Q5%. Lastly, IAHRIS offers the flood days per month to analyse the seasonality of the floods.

More information about these parameters appears in the IAHRIS methodological reference manual [74].

3. Results and Discussion

3.1. Calibration and Validation of the Model

The sensitivity analysis found the following parameters to be sensible: cn2, alpha, esco, perco, revap_co, epco, awc, flo_min, revap_min, surlag, and k (see

Table 1. Description, range, and adjusted value for all modified parameters after calibration.

Parameter	Description	Range	Adjusted Value
Cn2	The curve number relates to the soil permeability	−20–20%	−18%
Alpha	The baseflow recession constant measures the reaction of groundwater flow to recharge variations	0–1	0.35
Esco	Soil evaporation compensation factor	0–1	0.50
Perco	Percolation coefficient (it depends on the soil moisture)	−0.05–0.05	0.001
Revap_co	Groundwater re-evaporation coefficient	0–1	0.05
Epco	Compensation factor for plant uptake	0.01–1	0.50
Awc	The soil layer’s available water capacity	−20–20%	−0.5%
Flo_min	Minimum aquifer storage to allow return flow	0–5000	230.08
Revap_min	Threshold depth of water in the shallow aquifer needed for re-evaporation or percolation to the deep aquifer to occur	0–5000	1313.18
Surlag	Surface runoff lag coefficient	1–24	1.20
K	Saturated hydraulic conductivity	1–100	99.77
Snomelt_tmp	Snowmelt base temperature (°C)	−3–2.5	0.00
Snomelt_lag	Snowpack temperature lag factor	0–1	0.20

). For the SWE, we manually calibrated two parameters, the snomelt_tmp and snomelt_lag, which were modified. Table 1 shows the description of all modified parameters along with their range and the final value after calibration.

Most of the parameters that Cibin [75] identified as having a major impact on the streamflow simulation were sensible parameters in our model. The parameters related to the slope were not sensible, since this basin presents no significant changes in altitude, nor were the parameters related to the lateral flow. In this area, the most important parameters seem to be related to the groundwater and aquifers (alpha, revap_co, flo_min, and revap_min) and to the soil (cn2, esco, perco, awc, and k). Some soil qualities influence groundwater opposition, which warrants the inclusion of soil parameters as useful. These parameters have commonly been used in various works [76].

Most of the parameters fell within the expected range for a watershed in the Scandinavian peninsula. In addition, the parameters cn2, esco, perco, and snomelt_tmp were the best fitted to the expected values [77]. With the adjusted values for the parameters, the model performed well for the streamflow in the calibration period of 2007–2015 according to Moriasi’s standards [78]. The NSE was greater than 0.6, the PBIAS was close to 0, and the RSE approached 0.7 (

Table 2. Calibration and validation of the streamflow statistical daily values.

	NSE	PBIAS	RSR	R2
Calibration period (2007–2015)	0.7073	0.6192%	0.5410	0.7143
Validation period (2016–2020)	0.7704	−19.3319%	0.4791	0.7945

). To validate the model, we used the adjusted values of the parameters for 2016–2020 and applied the statistical evaluation indices to the daily streamflow. The results were also rather good, except for PBIAS.

During the validation period, the precipitation was low compared to the calibration period, which meant the streamflow was also rather low in those years (Figure 3). During the calibration period (2007–2015), the average annual precipitation is over 600 mm, while during the validation period (2016–2020), the average annual precipitation is below 550 mm. The decrease in PBIAS during the validation period suggests that the simulated flow was exaggerated, although it remained within acceptable limits. Because the calibration period was wetter than the validation period, the model overestimated the flow during dry periods. Despite this, the statistics are still acceptable, and the model is valid for both dry and wet periods. According to the statistical evaluation indices shown in Table 2, the SWAT+ model responded adequately to this alteration.

3.2. Assessment of the GCMs and the bias Correction Method

After determining the indicator Id (Equation (6)), the weights (Ib indices) for the creation of the ensemble of predictions of all models were calculated. The Id index represents the sum of the relative difference between the historical series of the model and the control scenario for the precipitation and temperature data. The Ib indicators were used as weights applied to the series obtained with each model to create an ensemble of predictions. Both indices can be seen in

Table 3. Id indicator for each GCM.

GCMs	Monthly Series
	Precipitation		Temperature		Id	Ib
	Id (Δx)	Id (Δσ)	Id (Δx)	Id (Δσ)
BCC	1.63	2.82	27.70	3.64	35.79	0.13
CanESM5	3.13	2.83	−1.12	3.07	7.91	0.16
EC-Earth-Veg	2.72	2.52	11.56	3.06	19.87	0.14
GFDL	3.27	2.31	27.30	2.92	35.80	0.13
INM-CM5	3.57	4.34	1.21	3.70	12.82	0.15
MiroC6	3.62	2.52	10.71	2.10	18.96	0.14
MRI	4.28	4.23	4.17	3.61	16.29	0.15

Note: Δx is the increase in the average and Δσ the increase in the variance.

.

The five bias-correction methods of the “qmap” R package were applied to correct the historical data; the results are shown in

Table 4. MAE of five bias-correction methods for observed versus corrected daily data (1985–2014).

GCMs	Variables	DIST	PFT	QUANT	RQUANT	SSPLIN
BCC	Precipitation	2.53	2.50	2.56	2.55	2.55
	Min. temp.	4.59	4.47	4.43	4.43	4.38
	Max. temp.	4.30	4.11	4.07	4.07	3.98
Can_ESM5	Precipitation	2.51	2.46	2.57	2.53	2.54
	Min. temp.	4.33	3.71	3.81	3.81	3.73
	Max. temp.	4.29	4.10	4.18	4.18	4.16
EC-Earth-Veg	Precipitation	2.53	2.51	2.58	2.55	2.56
	Min. temp.	4.16	3.83	3.90	3.90	3.90
	Max. temp.	3.86	3.82	3.87	3.87	3.87
GFDL	Precipitation	2.58	2.49	2.56	2.54	2.53
	Min. temp.	4.39	4.20	4.22	4.22	4.23
	Max. temp.	4.01	4.03	4.03	4.04	4.04
INM-CM5	Precipitation	2.55	2.51	2.55	2.53	2.55
	Min. temp.	4.32	3.51	3.64	3.64	3.64
	Max. temp.	4.25	3.84	4.04	4.04	4.05
MiroC6	Precipitation	2.53	2.50	2.56	2.55	2.60
	Min. temp.	3.89	3.68	3.80	3.80	3.80
	Max. temp.	3.91	3.63	3.71	3.71	3.71
MRI	Precipitation	2.51	2.47	2.54	2.52	2.52
	Min. temp.	4.46	4.01	4.09	4.09	4.01
	Max. temp.	4.09	3.91	3.99	3.99	3.99

. The PTF was found to be the best method for bias correction as it presented the lowest MAE for the three climate variables (daily precipitation, minimum temperature, and maximum temperature) in the majority of GCMs. Other studies have also shown that the PTF method provides the optimal correction [65]. Therefore, we applied the PTF method—which uses a parametric transformation of the quantile–quantile relation of observed and modelled values—to all the climate data of the seven GCMs for climate change forecasting.

The corrected data series of each model have been graphically compared with the uncorrected series and the real data in Figure 4. As can be seen in the graphs, both simulated temperature and precipitation improve considerably in most models when correcting the data, achieving results very similar to the observed values. This is especially noticeable in temperature, as there is a very appreciable improvement in all models, as can be seen in the graphs. Regarding precipitation, although in all the GCMs the corrected values are closer to the real values, the CanESM5, INM and MiroC6 models tend to underestimate precipitation during the summer months (June, July and August).

3.3. Climate Change Prevision

After the future temperature and precipitation data were corrected, the established periods were compared with the historical data (

Table 5. Historical and projected annual average precipitation and temperature data for the Erken basin for each GCM.

Model	Scenario	Periods	Temperature (°C)	Precipitation (mm)
BCC	Historic	1990–2014	6.65	553.45
	SSP 2–45	2026–2050	6.66 (+0.01 °C)	492.36 (−11%)
		2051–2075	7.25 (+0.6 °C)	551.00 (0%)
		2076–2100	7.61 (+0.9 °C)	520.79 (−6%)
	SSP 5–85	2026–2050	7.01 (+0.4 °C)	570.37 (+3%)
		2051–2075	8.14 (+1.5 °C)	549.62 (−1%)
		2076–2100	9.18 (+2.5 °C)	614.03 (+11%)
CanESM5	Historic	1990–2014	6.63	517.14
	SSP 2–45	2026–2050	8.79 (+2.2 °C)	591.74 (+14%)
		2051–2075	10.55 (+3.9 °C)	638.93 (+24%)
		2076–2100	10.94 (+4.3 °C)	653.29 (+26%)
	SSP 5–85	2026–2050	9.22 (+2.6 °C)	587.27 (+14%)
		2051–2075	11.25 (+4.6 °C)	698.50 (+35%)
		2076–2100	12.78 (+6.1 °C)	795.41 (+54%)
EC-Earth-Veg	Historic	1990–2014	6.74	574.70
	SSP 2–45	2026–2050	9.29 (+2.5 °C)	571.38 (−1%)
		2051–2075	9.46 (+2.7 °C)	575.15 (0%)
		2076–2100	10.35 (+3.6 °C)	588.11 (+2%)
	SSP 5–85	2026–2050	9.08 (+2.3 °C)	563.15 (−2%)
		2051–2075	10.41 (+3.7 °C)	611.52 (+6%)
		2076–2100	13.04 (6.3 °C)	682.18 (+19%)
GFDL	Historic	1990–2014	6.61	553.09
	SSP 2–45	2026–2050	8.29 (+1.7 °C)	569.19 (+3%)
		2051–2075	8.41 (+1.8 °C)	628.55 (+14%)
		2076–2100	8.84 (+2.2 °C)	597.79 (+21%)
	SSP 5–85	2026–2050	7.89 (+1.3 °C)	592.12 (+7%)
		2051–2075	9.41 (+2.8 °C)	630.93 (+14%)
		2076–2100	10.48 (+3.9 °C)	669.08 (+21%)
INM-CM5	Historic	1990–2014	6.78	568.92
	SSP 2–45	2026–2050	5.68 (−1.1 °C)	520.74 (−8%)
		2051–2075	5.89 (−0.9 °C)	532.45 (−6%)
		2076–2100	7.31 (+0.5 °C)	528.20 (−7%)
	SSP 5–85	2026–2050	6.14 (−0.6 °C)	500.27 (−12%)
		2051–2075	7.40 (+0.6 °C)	551.00 (−3%)
		2076–2100	9.03 (+2.3 °C)	533.60 (−6%)
MiroC6	Historic	1990–2014	6.78	541.29
	SSP 2–45	2026–2050	7.77 (+1 °C)	568.86 (+5%)
		2051–2075	8.47 (+1.7 °C)	564.86 (+4%)
		2076–2100	8.88 (+2.1 °C)	591.92 (+9%)
	SSP 5–85	2026–2050	7.74 (+1 °C)	608.62 (+12%)
		2051–2075	9.17 (+2.4 °C)	649.97 (+20%)
		2076–2100	10.35 (+3.6 °C)	686.85 (+27%)
MRI	Historic	1990–2014	6.64	537.19
	SSP 2–45	2026–2050	8.34 (+1.7 °C)	511.49 (−5%)
		2051–2075	8.80 (+2.2 °C)	527.29 (−2%)
		2076–2100	8.61 (+2 °C)	570.44 (+6%)
	SSP 5–85	2026–2050	8.94 (+2.3 °C)	561.42 (+5%)
		2051–2075	9.49 (+2.9 °C)	566.88 (+6%)
		2076–2100	10.04 (+3.4 °C)	609.89 (+14%)

Note: Values in parentheses show the increase in precipitation or temperature compared to the historical period.

). An increase or decrease in temperature and precipitation values was estimated for the coming years, for the near future (2026–2050), the medium term (2051–2075) and the long term (2076–2100).

There is considerable variability among the climate models, but the general trends are similar (Table 5). The INM-CM5 model is the only one that predicts a general decrease in temperatures in the short term, and a decrease in precipitation in the short and long term. The general trend of the rest of the models in both climate change scenarios is, temperature would experience a significant increase with respect to the historical period. In the intermediate scenario (SSP 2-45), the average temperature would increase to around 1 to 4 °C, and in a high emission scenario (SSP 5-85) the increase could be from 4 to 6 °C, depending on the GCMs. In the near future, the rising of temperature would be milder, but still concerning, around 2 °C in both scenarios. In most of these cases, in the far future period, the increase would be large enough to alter the hydrological cycle. One of the consequences will be increased stratification. The stratification in Lake Erken has been increasing for several decades [44] and other studies predict an intensification of this trend with serious implications [45].

Regarding precipitation, the variability among the GCMs is greater. Similarly, the general trend is that average annual precipitation would increase in a high-emission scenario. The INM-CM5 model shows a decrease in precipitation in both scenarios in the far future. The BCC and EC-Earth-Veg models show a decrease or a very slight increase in precipitation for the near and medium future, but an increase between 11 and 19% in the far future for the worst-case scenario (SSP 5-85). The MRI model shows a slight decrease in the near and medium future in the SSP 2-45 (−2,−5%), but a gradual increase in the SSP 5-85 (up to 14%). The GFDL model shows a gradual increase very similar for both scenarios up to 21%. Lastly, the MiroC6 and CanESM5 models show a gradual increase in precipitation for both scenarios, being higher in the SSP5 scenario. These results are congruent with those of other meteorological studies in Sweden. An increase in temperature and rainfall is expected in the north of Europe as a consequence of climate change [19,20,79,80,81].

The monthly distributions of the average temperature and precipitation were analysed in each period for both scenarios (Figure 5). In this figure, the ensemble precipitation and temperature data for the seven GCMs studied have been plotted. These analyses help in understanding the climate alterations in a general overview. Precipitation was found to increase in every scenario during the cold season (October, November, December, January and February). By contrast, during summer (June, July, and August) in some cases the precipitation was found to be lower, even though the average annual rainfall was higher according to most GCMs than during the historic period. The increase in temperature was similar in both scenarios for the near future but was more pronounced in the high-emission scenario (SSP5-85) for the far future. The annual distribution was unaltered because it increased uniformly during the whole year. However, the increase was slightly larger during summer and winter, especially in December, January and February.

Forecasts of seasonal patterns in precipitation and temperature in Sweden by other studies have been highly similar to our results [79,80,81,82]. Notably, the greatest temperature rise is expected to occur throughout the coldest months. As a result of this rise during December, January, February, and March, the hydrology could be altered, especially the processes related to snowfall and snow accumulation. This is because in all scenarios, the monthly average temperature for those months will rise above 0 °C in the far future (Figure 5). This would prevent snow accumulation in the basin. On the other hand, the increase in temperature in summer would increase the stratification in the lake [45].

The 5% probability exceedance distribution for precipitation was plotted to analyse the extreme meteorological events in the far future (Figure 6), using the ensemble data of the seven GCMS. The graphs show how the torrential precipitation is expected to increase in the period 2076–2100 in both the SSP 2-45 and SSP 5-85 scenarios. The increase is higher in the high emission scenario (SSP 5-85). This could imply that the increase in the annual average precipitation in the future (Table 5) would be related to an increase in the torrential events of precipitations.

3.4. Alterations on the Hydrological Cycle

The expected changes in precipitation and temperature will alter the hydrological cycle. To address these changes, we used the SWAT+ model for Lake Erken to simulate the future hydrological cycle for both scenarios, using the GCMs’ temperature and precipitation data. After simulating every GCMs model under the SSP 2-45 and SSP 5-85 scenario, we ensemble the hydrological components with the weighted average. The results show how the water balance is likely to be altered by climate change (

Table 6. Annual average of the water cycle components in the Erken basin. Ensemble of the GCMs.

Scenario	Periods	Discharge (m3/s)	Water Yield (mm)	Evapotrans. (mm)	SWE (mm)
Historic	1990–2014	0.36	50.27	354.44	4.51
SSP 2–45	2026–2050	0.37 (+3%)	46.68 (−7%)	352.36 (−1%)	2.77 (−38%)
	2051–2075	0.41 (+14%)	52.40 (+4%)	357.34 (+1%)	2.52 (−44%)
	2076–2100	0.43 (+18%)	52.74 (+5%)	358.34 (+1%)	1.58 (−65%)
SSP 5–85	2026–2050	0.43 (+20%)	52.44 (+4%)	363.12 (+2%)	2.32 (−49%)
	2051–2075	0.48 (+33%)	59.46 (+18%)	371.84 (+5%)	1.10 (−76%)
	2076–2100	0.53 (+49%)	66.05 (+31%)	379.36 (+7%)	0.38 (−92%)

Note: Values in parentheses show the increase in the component relative to the historical period.

).

In Table 6, “discharge” represents the water that enters Lake Erken through the different channels of the basin. “Water yield” represents the sum of the surface runoff, lateral soil flow, and tile flow; that is, the amount of water contributing to the flow rate. The SWE was defined in the introduction of this paper.

A gradual increase in the water discharge and water yield is the principal long-term pattern (Table 6). This effect of climate change is related to the increase in precipitation and torrential rainfall expected for the long-term scenario [18,19] (Figure 6). In the intermediate scenario, the increase in discharge could reach around 18% in the far future and 49% in the high-emission scenario. In this scenario, the water yield would also increase to around 31%. Such a sharp increase in discharge and runoff could imply a reduced availability of water resources since it would result in less water retention in the basin as well as soil degradation and erosion [16]. It is important to consider that forests buffer these effects and others caused by extreme torrential rainfall events [50]. Degradation of the watershed forest in future years could lead to even greater discharge, water yield, and soil degradation than those predicted in this study, with even heavier consequences for the ecosystem.

Regarding the SWE, even in the lowest-emission scenario, SWE was reduced by between 38 and 61% from 2026 to 2100. In the worst-case scenario, the reduction could reach 92% in the far future period. This is the expected result related to the temperature rise due to climate change. The increase in temperatures during December, January and February would reduce the snowfall and snow accumulation in the Lake Erken basin. As a consequence, in a high-emission scenario, the SWE would be reduced by nearly 95%. In other words, there would be almost no snow accumulation in the Erken basin during the cold season. The absence of snow accumulation could further favour runoff processes, as water would not be retained in the basin in the form of snow during the cold months.

Finally, a slight general increase in evapotranspiration in the basin is expected, especially in the high-emission scenario, where it could increase by 7% in the far future period. The increase in evapotranspiration is explained by the higher temperatures due to climate change.

3.5. Alterations on the Floods

Eight hydrological parameters related to floods were obtained from IAHRIS for each GCMs. The weighted average of the indicators is in

Table 7. Ensemble of the hydrological parameters related to floods of each GCMS (data obtained from IAHRIS).

Scenario	Periods	Qc (m3/s)	QL (m3/s)	Qconnec (m3/s)	Q5% (m3/s)	CV (Qc) (m3/s)	CV (Q5%) (m3/s)	Duration (days)
Historic	1990–2014	1.69	1.67	2.22	1.21	0.49	0.37	16.08
SSP 2–45	2026–2050	1.73	1.77	2.38	1.18	0.53	0.32	17.15
	2051–2075	1.92	2.04	2.79	1.32	0.57	0.30	22.11
	2076–2100	1.97	2.13	2.94	1.36	0.62	0.33	22.81
SSP 5–85	2026–2050	1.80	1.76	2.31	1.30	0.46	0.27	24.65
	2051–2075	2.18	2.47	3.48	1.47	0.64	0.30	29.82
	2076–2100	2.39	2.95	4.29	1.54	0.75	0.25	34.57

. We compared their values to observe how floods might be altered by climate change in future years. Their magnitude, frequency, variability, and duration were analysed for the different scenarios proposed.

According to the four parameters related to the magnitude of floods (Qc, QL, Qconnec and Q5%), in the Lake Erken basin floods would become progressively higher over time. This increase in magnitude will be greater in the SSP 5-85 scenario. The most affected parameter was Qconnec. This finding implies that the flood plain will flood more frequently, will increase in area, and will receive more water as heavy rainfall events increase in the basin. This pattern could change the ecosystem of the flood plain and its surroundings, altering the succession processes of the riparian forest and affecting the maintenance of the diversity and functionality of macroinvertebrate communities [83,84,85].

The variability of the floods was analysed with two parameters coefficient of variation. On the one hand, the variation of the mean of the maximum annual daily flow (CV (Qc)) would rise. This would imply that the magnitude of the floods would change depending on the year. Some years the magnitude of the floods would be higher, others would be lower, although the tendency would be to increase (Qc). On the other hand, the variation of the series of usual floods would decrease in future periods, although there is no clear trend of progressive decrease. This implies that the flow corresponding to the average flow curve classified at the 5% exceedance percentile will maintain similar values every year in the studied periods. In other words, some years there will be much more abundant floods than others, but the number of floods will be similar in most years.

The duration of the floods also showed extreme increases in the long term for the worst scenario. The increase in the flood duration or flood period could, in turn, raise the death rate of the most vulnerable plant species, as it causes anoxic stress [86]. The increase in floods may be related to the increase in precipitation and the increase in events of torrential precipitation (Figure 6).

Finally, we analysed the seasonality; the results are shown in

Table 8. Ensemble of the number of days in a month having a daily average flow rate greater than the Q5% (exceedance percentile) of each GCM.

Scenario	Periods	Jan	Feb	Mar	Apr	May	Jun	Jul	Aug	Sep	Oct	Nov	Dec
Historic	1990–2014	3.13	1.43	1.75	4.22	1.92	0.39	0.19	0.02	0.18	1.57	3.10	3.82
SSP 2–45	2026–2050	4.29	2.45	2.80	3.07	1.95	0.94	0.03	0.14	0.32	0.85	3.24	4.85
	2051–2075	6.52	3.71	2.99	2.82	1.48	0.28	0.06	0.03	0.18	1.68	4.28	6.50
	2076–2100	7.58	3.73	4.45	3.71	0.98	0.41	0.11	0.21	0.37	1.10	3.82	7.08
SSP 5–85	2026–2050	7.78	5.65	4.95	3.97	2.05	0.64	0.18	0.11	0.25	1.22	4.05	6.92
	2051–2075	9.65	6.06	5.14	2.84	1.20	0.33	0.15	0.24	0.54	1.68	5.15	9.76
	2076–2100	14.29	10.63	7.13	2.80	0.92	0.36	0.05	0.17	0.20	1.29	5.56	12.10

. Changes in the natural seasonal patterns of floods can produce strong distortions in the riparian ecosystem. Loss of synchrony with the life cycles of the affected species is the main reason for such distortions, in addition to many accumulating effects from other environmental variables. Some of these effects could include alteration in biomass production [87]; loss of synchrony with the phenology of multiple plant species, altering dispersal and germination processes [84,87]; and progression of generalist foreign plant species, with a resulting decrease in diversity and an impact on riparian forest productivity [84].

An increase in seasonality is observed in Table 8. During the warm months (April to October), the days with floods are expected to remain similar to the historical period or even to decrease. The greater decrease would be in spring, during the month of April, where there would be 1 or 2 days less of floods in the far future. In the cold months, the days with floods would increase significantly, especially during December, January, and February. The increase in winter floods could be related to increased precipitation and a lack of solid snow in the Erken basin. The snow precipitation that accumulates until the thawing months in these scenarios would contribute directly to the watercourse [88].

4. Conclusions

In this work, we analysed the impact of climate change on the hydrological cycle of the Lake Erken basin under both an intermediate and a high-emission scenario. We also developed a SWAT+ model of the watershed and used precipitation and temperature data obtained from seven GCMs (BCC, CanESM5 GCM, EC-Earth-Veg, GFDL, INM-CM5, MiroC6, MRI). The main conclusions drawn from this study are as follows:

• After calibration and validation, the SWAT+ model of Lake Erken performed notably well regarding the examined statistics. Although SWAT+ is a new tool and has not been as widely used compared to its predecessor (SWAT), our results demonstrated its validity and effectiveness. We conclude that the SWAT+ can be used effectively to create hydrological models in lake basins. Furthermore, observing the concordance of the results obtained with SWAT+ and IAHRIS, it seems to respond acceptably to flood and peak flow events.
• Precipitation and temperature are expected to increase in Sweden in the future because of climate change. The results of this study are in concordance with that forecast. In a high-emission scenario (SSP 5-85), the temperature could increase from 4 to 6 °C in the distant future (2076–2100), especially during winter. Precipitation would increase by around 20%.
• An increase in precipitation linked to torrential events would affect the hydrological cycle of the Lake Erken basin. In the distant future, the discharge can be expected to increase by almost 50% and runoff by 30% in a high-emission scenario. This dynamic could intensify soil erosion and degradation processes, which could damage the watershed’s forests. As the duration of water retention in the watershed decreases, the availability of water resources for vegetation may also decrease.
• The increase in temperature during the winter months would affect the snow-related processes of the basin. The accumulation of snow in the basin would be severely altered in both scenarios. In the intermediate scenario, it would be reduced by 65%, and in the high-emission scenario, by 92%. This would imply the near disappearance of the snowpack that usually covers the basin in the winter months. Such a change would drastically affect all forest ecosystems in the basin.
• In the worst scenario (SSP5-85), the magnitude, frequency and variability of floods are expected to increase. The duration of floods would also increase from 16 days to 23 or 35, depending on the scenario. Regarding the seasonality of floods, the number of days in a month having a daily average flow rate greater than the exceedance percentile would increase in winter, especially in February (from 1 to 10 days), January (from 3 to 14 days) and December (from 4 to 12 days), while in the warmer months they would decrease. These changes may cause damage to the ecosystems of the Lake Erken basin and cause a degradation of its soil due to erosion.

5. Limitations

• The present study has several limitations, some of them due to data availability, others due to the methodology of the study itself or the tools used. In any case, these limitations should be considered: The SWAT program has been widely used, but its most recent version, SWAT+, does not have as many studies to support its use. In addition, many SWAT studies have been conducted simulating peaks or flood flows and although acceptable, SWAT models tend to underestimate these peak flows, as explained in the introduction of this study. There are not enough studies with SWAT+ to know if peak simulation is also underestimated by this version.
• Only one meteorological station has been used in this study, even though a single station is acceptable due to the small size of the basin, according to the World Meteorological Organisation [89].
• The inflow data for Lake Erken at Kristineholm did not start to be collected until mid-2006. In order to calibrate and validate the SWAT+ model, daily observed discharged data of full years are required. As a consequence, our calibration and validation period go from 2007 to 2020 and is considerably shorter than the historical period, for which observed precipitation and temperature data are available.
• CMIP6 includes over 100 global climate models from more than 50 modelling centers. In this study, we have used seven of these models. This allows us to decrease the uncertainty of the results but not to reduce it to zero.