Abstract
Pressure on water resources has reached unprecedented levels during the last decades because of climate change, industrialization, and population growth. As a result, vulnerability to inappropriate water availability and/or quality is increasing worldwide. In this paper, a Soil & Water Assessment Tool (SWAT) model of the Carp River watershed located in the city of Ottawa, Ontario was calibrated and validated. The model was then used to evaluate the individual and coupled impacts of urbanization and climate change on water quantity (discharge) and quality (nitrogen and phosphorus loads). While most of the watershed is currently rural, the headwaters will undergo rapid urbanization in the future, and there are concerns about possible negative impacts on water quantity and quality. Seven scenarios were developed to represent various watershed configurations in terms of land use and climate regime. Future climate time series were obtained by statistically downscaling the outputs of nine regional climate models, run under representative concentration pathways (RCP)4.5 and RCP8.5. The impacts were evaluated at the main outlet and at the outlet of an upstream sub-watershed that would be most affected by urbanization. Results show that climate change and urbanization's impacts vary greatly depending on the spatial scale and geographic location. Globally, the annual average discharge will increase between 6.75 and 9.34% by 2050, while changes in annual average nitrogen and phosphorus loads will vary between −1.20 and 24.84%, and 19.15 and 23.81%, respectively. Local impacts in sub-watersheds undergoing rapid urbanization would be often much larger than watershed-scale impacts.
HIGHLIGHTS
A global-local approach was used to assess the impacts.
Local impacts can be an order of magnitude higher than watershed-scale impacts.
Impacts at the local and watershed scale might have opposite directions of change.
Climate change will likely drive the flow and nitrogen load in the Carp River watershed.
Urbanization will control the phosphorus load.
INTRODUCTION
During the last decades, pressure on natural resources has been steadily increasing, mainly driven by industrialization and population growth. The United Nations estimated that the world population would increase by 2 billion in the next 30 years, from 7.7 billion in 2020 to 9.7 billion in 2050 (United Nations 2019). This demographic change will be associated with higher standards of living, increased demand for land and water, and increased pollution. This will inevitably affect water supply, possibly disrupting domestic and industrial water uses, as well as environmental sustainability. The consequence will be lower access to essential water-related services, decreased economic prosperity, and increased conflicts around water bodies. For instance, Munia et al. (2020) reported that the number of people exposed to water stress would increase by 50–100%, depending on the emission scenario and population growth rate.
The Mississippi Valley Conservation Authority (MVCA) is the responsible authority for the sustainable management of the Carp River watershed located in Ontario, Canada (Figure 1). The watershed is currently mainly in a rural state; however, the headwaters are expected to undergo rapid urbanization in the future. The MVCA anticipates significant urban developments in the upcoming decades. It is expected that both urbanization and climate change will induce substantial changes in the watershed's hydrology, which in turn will affect its ecosystem health through changes in water quantity and quality. The results of this study could support authorities in planning and managing water resources and providing information for sound decision-making.
The combined effects of climate change and land-use change on water quantity and/or quality have been examined by many authors (Tu 2009; Dunn et al. 2012; El-Khoury et al. 2015; Karlsson et al. 2016; Kundu et al. 2017; Čerkasova et al. 2018; Choukri et al. 2020; Guo et al. 2020; Wang et al. 2020). These studies differ in the location of the study area, the hydrological model, the type of assessed water quantity and quality variables, and the way climate change and land-use scenarios were generated. These studies involved numerical experiments where calibrated hydrological models incorporated climatic inputs representing one or several global warming scenarios and using current and/or projected land-use maps. Although these studies confirmed the impact of urban development and climate change on streamflow, the evaluation was only limited to the effects at the main outlet (global level). In this study, the impacts at the local scale were also considered by estimating the change at the outlet of the sub-watershed that was most affected by urbanization.
Since 2005, several studies related to the preparation of the Carp River model were conducted: City of Ottawa (2016), Greenland International Consulting Ltd (2009, 2010, 2011, 2014), Stantec Consulting Ltd (2006), and CH2M Hill (2005, 2006). All these works focused on flood level analysis and restoration plans, but not climate change or urbanization. A limited number of studies addressing the combined effects of climate change and urbanization in Canada can be found in the literature. Of particular interest is El-Khoury et al. (2015), who assessed the impact of climate change and deforestation on the South Nation watershed in opposition to urbanization.
Hydrological impacts are generally assessed by simulating hydrological processes and water exchange within an area. The impact of natural variability and anthropogenic factors on a hydrologic system can be assessed by changing the model's climatic inputs or watershed configuration (Yu 2015). The application of rainfall–runoff models varies from straightforward water balance models to more complex watershed models, such as the Annualized Agricultural Non-Point Source model (Young et al. 1989) and the Soil Water & Assessment Tool (SWAT) (Arnold et al. 1998), which can simulate loads and concentrations of various water quality variables. The selection of the appropriate model depends on several factors such as scale of the problem, computational cost, and robustness of the model. There is no single best model, as there are many plausible solutions (Ogden 2021). Therefore, a typical model selection is more guided by familiarity with the model than appropriateness (Addor & Melsen 2019). The SWAT model was selected in this research, as it has been widely used to assess the hydrology in small- and large-scale catchments: Huo & Li (2013), Khoiab & Thomb (2015), and Jajarmizadeh et al. (2017) in Asia; Arnold et al. (2012), El-Khoury et al. (2015), and Havrylenko et al. (2016) in America; Abbaspour et al. (2015) in Europe; and Mengistu et al. (2019), Choukri et al. (2020), and Muluneh (2020) in Africa. The SWAT is also used to investigate the relative impact of climate, land-use & land-cover (LULC) changes, and management practices on the available water resources (Arnold et al. 1998; Neitsch et al. 2009; Tirupathi & Shashidhar 2020). In addition, its structure and applicability to address different problems make it more flexible. The SWAT also allows for the simulation of water quantity and quality processes within the same model (e.g., flow, N, and P) without the contribution of other programs. For instance, the precedent model of the municipality developed by the municipality combines several hydrologic platforms and a separate hydraulic model.
Several approaches have been reported in the literature for generating climate change projections: Wang et al. (2020) utilized the Delta-change method, which consisted of adding an arbitrarily chosen (but realistic) variation to historical observations to represent global warming. Others used the outputs of global and regional climate models (RCMs) without any processing (Tu 2009; Čerkasova et al. 2018), or in combination with a statistical downscaling approach such as the Delta-change method (e.g., Dunn et al. 2012) and bias correction (El-Khoury et al. 2015; Karlsson et al. 2016; Kundu et al. 2017; Choukri et al. 2020; Guo et al. 2020). The Delta-change approach and the bias correction method have their advantages and weaknesses. They are straightforward in the generation of the climate output. However, data generated from raw global data outputs with the Delta method cannot capture changes in climate variability, and bias correction does not remove distortion because the approaches may be too simplistic to characterize climate trends (El-Khoury et al. 2015). In addition, the choice of the downscaling method is very determinant, as it significantly affects the outputs. For instance, Miralha et al. (2020) have demonstrated that the choice of bias correction method affects the direction of change and the magnitude of nutrient loads (N and P) and hydrological processes. Chen et al. (2012) have found that for the same global climate model (GCM), the simulated runoffs vary significantly when using rainfall provided by different statistical downscaling techniques as the input to the hydrological models. Therefore, reasonable care must be taken in the selection of downscaling method to represent future hydrological scenarios that are more realistic.
This paper aims to evaluate the individual and coupled impacts of urbanization and climate change on water quantity (discharge) and quality (nitrogen and phosphorus loads) in the carp watershed. The outputs of nine RCMs were downscaled using quantile mapping. The results are examined at two scales: the watershed scale (i.e., at the outlet of the watershed) and the local scale (i.e., at the outlet of the sub-watershed which undergoes the most drastic urbanization).
MATERIALS AND METHODS
Study area
The Carp River is an important natural ecosystem located west of the City of Ottawa, ON, Canada. It is 42 km long and drains an area of approximately 265 km². It has its headwaters in the Glen Cairn area of Kanata and flows north into the Ottawa River at Fitzroy Harbour. The Carp River is one of the main sources of the City of Ottawa's drinking water. The annual average precipitation of the historical period of simulation (1990–2018) is 933 mm. The daily minimum, mean, and maximum of the observed streamflow at outlet 8 (Figure 2) during the same period are 0.032, 2.95, and 23.7 m3/s, respectively. The highest elevation has an altitude of 203 m above sea level, while the outlet lies at 59 m. The average elevation measures 117 m above sea level. Based on the latest provincial land-use update in 2017 (Ministry of Natural Resources and Forestry 2019), the watershed is classified into 10 different land-use types with a dominance of agricultural land, covering more than half of the total area (52%). Urbanized area represents 10% (Figure 1). Mixed forests consisted of two types (FRST and FOMI) and only differ on biomass die-fraction. Tables 1 and 2 show the percentage of coverage of each land-use/-cover and soil type within the watershed.
Land use/cover type . | SWAT code . | Area (%) . |
---|---|---|
1. Agricultural land | AGRC | 52.34 |
2. Wooded wetland | WEWO | 18.02 |
3. Evergreen needle leaf forest | FOEN | 9.93 |
4. Urban | URBN | 9.90 |
5. Forest – deciduous | FRSD | 4.16 |
6. Forest-mixeda | FRST | 3.25 |
7. Barren or sparsely vegetated | BSVG | 0.82 |
8. Water | WATR | 0.69 |
9. Mixed foresta | FOMI | 0.55 |
10. Cropland/woodland mosaic | CRWO | 0.34 |
Land use/cover type . | SWAT code . | Area (%) . |
---|---|---|
1. Agricultural land | AGRC | 52.34 |
2. Wooded wetland | WEWO | 18.02 |
3. Evergreen needle leaf forest | FOEN | 9.93 |
4. Urban | URBN | 9.90 |
5. Forest – deciduous | FRSD | 4.16 |
6. Forest-mixeda | FRST | 3.25 |
7. Barren or sparsely vegetated | BSVG | 0.82 |
8. Water | WATR | 0.69 |
9. Mixed foresta | FOMI | 0.55 |
10. Cropland/woodland mosaic | CRWO | 0.34 |
aTwo types of mixed forests defined in the SWAT with two different codes (FRST and FOMI).
Soil name . | SWAT code . | Area (%) . |
---|---|---|
Loam*1 | Be1-2a-4649 | 51.29 |
Loam*2 | Be1-2a-4648 | 40.30 |
Sandy Loam | Po2-1-2b-4971 | 6.33 |
Clay | Gm3-3a-3070 | 2.08 |
Soil name . | SWAT code . | Area (%) . |
---|---|---|
Loam*1 | Be1-2a-4649 | 51.29 |
Loam*2 | Be1-2a-4648 | 40.30 |
Sandy Loam | Po2-1-2b-4971 | 6.33 |
Clay | Gm3-3a-3070 | 2.08 |
Note: The two types of loam differ in their available water capacity: 0.08 mm/mm (*1) and 0.175 mm/mm (*2).
Model setup
The SWAT (Arnold et al. 1998) is a physically-based, semi-distributed, and continuous model used to evaluate and predict the impacts of climate change and land-use management on water quantity and quality. The SWAT operates on a daily time step and can simulate continuously over a long period. The procedure used by the SWAT consists of dividing the watershed into several sub-watersheds, which are further subdivided into hydrologic response units (HRUs). An HRU is a unique combination of topographical, land-use management, and soil characteristics (Arnold et al. 2012).
S is expressed in terms of a CN, a dimensionless watershed parameter ranging from 0 to 100. Groundwater is represented by two reservoirs (a shallow and a deep aquifer), which interact with surface water through infiltration, percolation, evaporation, lateral flow, and baseflow. The SWAT also simulates the complete nutrient cycle (transformation and movement) for nitrogen and phosphorus as well as the degradation of any pesticides applied in an HRU. Nutrient levels are computed on a mass basis, while it allows the input data as concentrations.
The SWAT model was set up with ArcSWAT, a GIS preprocessor of the SWAT model that combines the spatial and climatic inputs and converts them into a SWAT model. ArcSWAT requires a few subjective user inputs, including the location of sub-watershed outlets and the HRUs definition method (i.e., how soil, land use, and slope are classified and used to define HRUs). Options for HRU definition are:
- (1)
multiple HRUs per sub-watershed representing all combinations of the slope, land-use, and soil classes;
- (2)
a single HRU per sub-watershed corresponding to the dominant slope class, dominant land-use, and dominant soil classes; and
- (3)
a single HRU corresponding to the most extensive combinations of the slope, land-use, and soil classes.
In this paper, the first option was selected. Given the relatively flat topography of the area, a single slope class was considered.
Different land-use changes and climate change scenarios were developed by downscaling the outputs of nine RCMs; the projected climate time series were used to force the SWAT model, and changes in water quantity and water quality parameters were calculated.
Input data
The model setup requires three types of input data (described in Table 3). Observed water quantity and quality time series of the parameters of interest are required for calibration and validation. In this paper, time series of streamflow, total nitrogen, and total phosphorus were used. The initial land-use classes were reviewed and reclassified based on land-use types defined in the SWAT database. Observed daily precipitations and temperature were obtained from the Environment and Climate Change Canada (Environment and Natural Resources Canada 2020) at two locations in the area of Ottawa Airport. Other meteorological parameters (humidity, wind speed, and solar radiations) were retrieved from the Watch-Forcing Data ERA-Interim (Weedom et al. 2014). Streamflow, Total N, and Total P load data were available for different periods: measured streamflow at outlet 8 covers the 1990–2015 period, while records of total N and total P for outlets 8 and 30 were available for 2000–2018 and 2000–2007 (Figure 4). These observed water quantity and quality data were used for model calibration and validation. Table 4 presents the characteristics of the record stations used.
Data . | Description . | Temporal scale . | Spatial scale . | Period . | Source . |
---|---|---|---|---|---|
Topography | Elevation – LiDAR | N/A | 1 m | 2015 | City of Ottawa |
Land use and cover | Categories of land occupation (wood, urban, and agricultural) | N/A | 30 m | 2017 | Ministry of Natural Resources and Forestry of Canada (2019) |
Soil | Soil types and physical properties | N/A | 10 km | 2007 | Food and Agricultural Organization (2019) |
Weather | Precipitation and temperature | Daily | Point | 1990–2018 | Environment and Climate Change Canada (ECCC) |
Solar radiation, relative humidity, and wind speed | Daily | 0.5° | 1990–2018 | WFDEI | |
Water quantity | Streamflow | Daily (one measurement per month) | Point | 1990–2015 | HYDAT – Environment Canada |
Water quality | N, P | Monthly | 2000–2018 2000–2007 | Ministry of Environment (Ontario) |
Data . | Description . | Temporal scale . | Spatial scale . | Period . | Source . |
---|---|---|---|---|---|
Topography | Elevation – LiDAR | N/A | 1 m | 2015 | City of Ottawa |
Land use and cover | Categories of land occupation (wood, urban, and agricultural) | N/A | 30 m | 2017 | Ministry of Natural Resources and Forestry of Canada (2019) |
Soil | Soil types and physical properties | N/A | 10 km | 2007 | Food and Agricultural Organization (2019) |
Weather | Precipitation and temperature | Daily | Point | 1990–2018 | Environment and Climate Change Canada (ECCC) |
Solar radiation, relative humidity, and wind speed | Daily | 0.5° | 1990–2018 | WFDEI | |
Water quantity | Streamflow | Daily (one measurement per month) | Point | 1990–2015 | HYDAT – Environment Canada |
Water quality | N, P | Monthly | 2000–2018 2000–2007 | Ministry of Environment (Ontario) |
N/A, not applicable.
Model performance measures
NS ranges from –∞ to 1. The optimal value is 1, corresponding to the situation where the plot of observed data fits the simulation (Khoiab & Thomb 2015) perfectly, while values less than 0 indicate that the observed data mean is a more accurate predictor than the simulated output (Jajarmizadeh et al. 2017). The target interval for PBIAS is zero to ±25%, with zero corresponding to the optimum. Positive values indicate the model is underestimating the observations, while negative values indicate an overestimation.
The statistical measures are compared to the criteria illustrated in Table 5 (Moriasi et al. 2007) to classify the model's performance.
. | Station name (Climate ID) . | Longitude . | Latitude . | Distance to the watershed boundary (km) . |
---|---|---|---|---|
1 | Ottawa CDA (6105976) | −75.72 | 45.38 | 16.5 |
2 | Ottawa Macdonald-Cartier Int'l (6106000) | −75.67 | 45.32 | 22 |
A | – | −75.25 | 45.75 | 30 |
B | – | −75.75 | 45.75 | 42 |
C | – | −75.25 | 45.25 | 11 |
D | – | −75.75 | 45.25 | 11 |
. | Station name (Climate ID) . | Longitude . | Latitude . | Distance to the watershed boundary (km) . |
---|---|---|---|---|
1 | Ottawa CDA (6105976) | −75.72 | 45.38 | 16.5 |
2 | Ottawa Macdonald-Cartier Int'l (6106000) | −75.67 | 45.32 | 22 |
A | – | −75.25 | 45.75 | 30 |
B | – | −75.75 | 45.75 | 42 |
C | – | −75.25 | 45.25 | 11 |
D | – | −75.75 | 45.25 | 11 |
. | Satisfactory . | Good . | Very good . |
---|---|---|---|
NS | 0.5–0.7 | 0.7–0.8 | 0.8–1.0 |
PBIAS (%) | 15–25 | 10–15 | Less than 10% |
. | Satisfactory . | Good . | Very good . |
---|---|---|---|
NS | 0.5–0.7 | 0.7–0.8 | 0.8–1.0 |
PBIAS (%) | 15–25 | 10–15 | Less than 10% |
Sensitivity analysis, calibration, and validation
Calibration consists of modifying model parameters to achieve a good fit between the simulated output and the observed data. Validation checks the model's performance using the calibrated parameters and an independent data set (not used in the calibration process). The simulation period is split into calibration and validation periods. Table 6 presents how the period of data was used for the calibration and the validation.
Variable . | Period . | ||
---|---|---|---|
Calibration . | Validation . | ||
Discharge | Outlet 8 | 1990–2005 | 2007–2015 |
N | Outlet 8 | 2000–2012 | 2013–2018 |
Outlet 30 | 2000–2003 | 2004–2006 | |
P | Outlet 8 | 2000–2012 | 2013–2018 |
Outlet 30 | 2000–2003 | 2004–2006 |
Variable . | Period . | ||
---|---|---|---|
Calibration . | Validation . | ||
Discharge | Outlet 8 | 1990–2005 | 2007–2015 |
N | Outlet 8 | 2000–2012 | 2013–2018 |
Outlet 30 | 2000–2003 | 2004–2006 | |
P | Outlet 8 | 2000–2012 | 2013–2018 |
Outlet 30 | 2000–2003 | 2004–2006 |
The mean of the variations in the objective function estimates the sensitivity. The t-stat is the regression coefficient of a parameter divided by its standard error. The p-value for each parameter tests the null hypothesis that the regression coefficient is equal to zero. The higher the absolute value of t-stat and the smaller the value of p-value, the more sensitive the parameter is (Abbaspour 2007). A small p-value, typically lower than 0.05, means the sensitivity of the parameter is statistically significant (Verhagen et al. 2004), while a value of 0.05 indicates that there is a 95% probability that a parameter change will affect the dependent variable (Abbaspour et al. 2009).
Both the sensitivity analysis and the calibration were conducted using the SUFI2 (Sequential Uncertainty Fitting) algorithm implemented in the SWAT-CUP program (Abbaspour 2007). The parameters included in the sensitivity analysis are shown in Table 7.
Parameters . | Description . | Variation range . |
---|---|---|
r__CN2.mgt | Initial SCS CNII value | −0.2 to 0.5 |
v__ESCO.hru | Soil evaporation compensation factor | 0.01–1 |
v__GWQMN.gw | Threshold water depth in the shallow aquifer for flow | 0–2,500 |
v__ALPHA_BF.gw | Baseflow alpha factor | 0–1 |
a__REVAPMN.gw | Threshold depth of water in the shallow aquifer for ‘revap’ to occur (mm) | 0–500 |
v__SOL_AWC.sol | Available water capacity | 0–0.8 |
a__OV_N.hru | Manning's n value for overland flow | 0.01–20 |
v__GW_REVAP.gw | Groundwater ‘revap’ coefficient | 0.02–0.15 |
a__SLSUBBSN.hru | Average slope length | 10–130 |
v__SOL_K.sol | Saturated hydraulic conductivity | 0–1 |
a__RCN.bsn | Concentration of N in rainfall | 0–15 |
a__ERORGN.hru | Organic N enrichment ratio | 0–5 |
a__SHALLST_N.gw | Concentration of nitrate in groundwater contribution to streamflow from sub-basin (mg N/l) | 0–1,000 |
a__N_UPDIS.bsn | N uptake distribution parameter | 0–100 |
a__NPERCO.bsn | N percolation coefficient | 0–1 |
a__ANION_EXCL.sol | Fraction of porosity (void space) from which anions are excluded | 0.01–1 |
a__P_UPDIS.bsn | P uptake distribution parameter | 0–100 |
a__PPERCO.bsn | P percolation coefficient | 10–17.5 |
a__ERORGP.hru | Organic P enrichment ratio | 0–5 |
a__GWSOLP.gw | Concentration of soluble P in groundwater contribution to streamflow from sub-basin (mg P/l) | 0–1,000 |
a__PSP.bsn | P sorption coefficient | 0.01–0.7 |
a__PHOSKD.bsn | P soil partitioning coefficient | 100–200 |
Parameters . | Description . | Variation range . |
---|---|---|
r__CN2.mgt | Initial SCS CNII value | −0.2 to 0.5 |
v__ESCO.hru | Soil evaporation compensation factor | 0.01–1 |
v__GWQMN.gw | Threshold water depth in the shallow aquifer for flow | 0–2,500 |
v__ALPHA_BF.gw | Baseflow alpha factor | 0–1 |
a__REVAPMN.gw | Threshold depth of water in the shallow aquifer for ‘revap’ to occur (mm) | 0–500 |
v__SOL_AWC.sol | Available water capacity | 0–0.8 |
a__OV_N.hru | Manning's n value for overland flow | 0.01–20 |
v__GW_REVAP.gw | Groundwater ‘revap’ coefficient | 0.02–0.15 |
a__SLSUBBSN.hru | Average slope length | 10–130 |
v__SOL_K.sol | Saturated hydraulic conductivity | 0–1 |
a__RCN.bsn | Concentration of N in rainfall | 0–15 |
a__ERORGN.hru | Organic N enrichment ratio | 0–5 |
a__SHALLST_N.gw | Concentration of nitrate in groundwater contribution to streamflow from sub-basin (mg N/l) | 0–1,000 |
a__N_UPDIS.bsn | N uptake distribution parameter | 0–100 |
a__NPERCO.bsn | N percolation coefficient | 0–1 |
a__ANION_EXCL.sol | Fraction of porosity (void space) from which anions are excluded | 0.01–1 |
a__P_UPDIS.bsn | P uptake distribution parameter | 0–100 |
a__PPERCO.bsn | P percolation coefficient | 10–17.5 |
a__ERORGP.hru | Organic P enrichment ratio | 0–5 |
a__GWSOLP.gw | Concentration of soluble P in groundwater contribution to streamflow from sub-basin (mg P/l) | 0–1,000 |
a__PSP.bsn | P sorption coefficient | 0.01–0.7 |
a__PHOSKD.bsn | P soil partitioning coefficient | 100–200 |
Replace or v_: the existing parameter value is to be replaced by the given value.
Absolute or a_: the existing parameter value is added to a given value.
Relative or r_: (1+a given value) multiply the existing parameter value.
Land-use projections
The projected land-use map was elaborated by the City of Ottawa in collaboration with the MVCA based on projections of the next 30 years (by 2050). The land-use scenario consisted of a change in the land-use area ratios and the type of urban area: from residential-medium density to commercial in the future. This conversion allows predicting impacts in severe conditions. The two land-use types differ in several characteristics, including the fraction of total impervious area, fraction directly connected to impervious, mass fraction of total N and total P in suspended solid load from the impervious area. No change in management operations in agricultural lands was completed between current and future projected conditions.
Table 8 presents the percentages of change between the current land uses and projected land uses within the watershed (map shown in Figure 3).
Land-use type . | Change (%) . |
---|---|
1. Agricultural land | (−2.63) |
2. Wooded wetland | (−0.54) |
3. Evergreen needle leaf forest | (−0.13) |
4. Urban | (+3.64) |
5. Forest – deciduous | (−0.17) |
6. Forest-mixed (FRST) | (−0.17) |
7. Barren or sparsely vegetated | (+0) |
8. Water | (+0) |
9. Mixed forest (FOMI) | (−0.02) |
10. Cropland/woodland mosaic | (+0) |
Land-use type . | Change (%) . |
---|---|
1. Agricultural land | (−2.63) |
2. Wooded wetland | (−0.54) |
3. Evergreen needle leaf forest | (−0.13) |
4. Urban | (+3.64) |
5. Forest – deciduous | (−0.17) |
6. Forest-mixed (FRST) | (−0.17) |
7. Barren or sparsely vegetated | (+0) |
8. Water | (+0) |
9. Mixed forest (FOMI) | (−0.02) |
10. Cropland/woodland mosaic | (+0) |
bold: increase; italic: decrease
Climate change scenarios
Developing a climate change scenario involves considerable inherent uncertainties (Shrestha et al. 2020). As these scenarios will be used in planning and decision-making processes, it is crucial to estimate a range of plausible future climate conditions in studies of climate change impact. Therefore, using a single RCM or carbon emission scenario is not sufficient for assessment purposes, since it cannot provide an uncertainty range. Consequently, a combination of three RCMs with two emission scenarios were used for this study. Scenarios are defined based on IPCC AR5 representative concentration pathway (RCP)4.5 and RCP8.5 (Collins et al. 2013). The RCP4.5 represents medium emissions, while RCP8.5 provides future climate data with high emissions assuming high population and relatively slow income growth with modest rates of technological change and energy intensity improvements resulting in the long term to high-energy demand and greenhouse gas (GHG) emissions without any climate change policies (Wayne 2013). Various RCMs were used to generate the projected future climate data at smaller scale driven by GCMs from several institutes as listed in Table 9.
No. . | Institute . | Institute acronym . | Driving GCMs . | RCMs . |
---|---|---|---|---|
1 | Canadian Centre for Climate Modeling and Analysis | CCCma | CanESM2 | CANRCM4 |
2 | Institut Catholique des Hautes Études Commerciales | ICHEC | EC-EARTH | RCA4, HIRHAM5 |
3 | Met Office Hadley Centre | MOHC | HadGEM2-ES | WFR, RegCM4 |
4 | Max Planck Institute for Meteorology | MPI-M | MPI-ESM-LR | WFR, RegCM4 |
5 | National Oceanic and Atmospheric Administration – Geophysical Fluid Dynamics Laboratory | NOAA-GFDL | GFDL-ESM2M | WFR, RegCM4 |
No. . | Institute . | Institute acronym . | Driving GCMs . | RCMs . |
---|---|---|---|---|
1 | Canadian Centre for Climate Modeling and Analysis | CCCma | CanESM2 | CANRCM4 |
2 | Institut Catholique des Hautes Études Commerciales | ICHEC | EC-EARTH | RCA4, HIRHAM5 |
3 | Met Office Hadley Centre | MOHC | HadGEM2-ES | WFR, RegCM4 |
4 | Max Planck Institute for Meteorology | MPI-M | MPI-ESM-LR | WFR, RegCM4 |
5 | National Oceanic and Atmospheric Administration – Geophysical Fluid Dynamics Laboratory | NOAA-GFDL | GFDL-ESM2M | WFR, RegCM4 |
No. . | Scenarios . | Climate data frame . | Land-use map . | Description . |
---|---|---|---|---|
1. | S0o | Historical observations (1990–2018) | 2017 | Current model (Baseline) |
2. | S0m | Historical model (1990–2018) | 2017 | Baseline scenario for climate change and combined effects scenarios |
3. | S1 | Historical observations (1990–2018) | 2050 | Land-use change |
4. | S0M45 | RCP4.5 (2021–2050) | 2017 | RCP4.5 climate change |
5. | S0M85 | RCP8.5 (2021–2050) | 2017 | RCP8.5 climate change |
6. | S1M45 | RCP4.5 (2021–2050) | 2050 | Combined land-use change and climate change under RCP4.5 |
7. | S1M85 | RCP8.5 (2021–2050) | 2050 | Combined land-use change and climate change under RCP8.5 |
No. . | Scenarios . | Climate data frame . | Land-use map . | Description . |
---|---|---|---|---|
1. | S0o | Historical observations (1990–2018) | 2017 | Current model (Baseline) |
2. | S0m | Historical model (1990–2018) | 2017 | Baseline scenario for climate change and combined effects scenarios |
3. | S1 | Historical observations (1990–2018) | 2050 | Land-use change |
4. | S0M45 | RCP4.5 (2021–2050) | 2017 | RCP4.5 climate change |
5. | S0M85 | RCP8.5 (2021–2050) | 2017 | RCP8.5 climate change |
6. | S1M45 | RCP4.5 (2021–2050) | 2050 | Combined land-use change and climate change under RCP4.5 |
7. | S1M85 | RCP8.5 (2021–2050) | 2050 | Combined land-use change and climate change under RCP8.5 |
The 1985–2100 precipitation, temperature (maximum and minimum), solar radiation, relative humidity, and wind speed data sets were simulated using the RCMs under RCP4.5 and RCP8.5 scenarios. Data were downloaded from North America – CORDEX (Coordinated Downscaling Experiment: Mearns et al. 2017). The time series for the five variables were extracted for the period of 2021–2050 and downscaled using the quantile mapping, also known as quantile–quantile or quantile matching.
Seven scenarios
Based on the climate change and land-use change projections, the seven following experiments, presented in Table 10, were developed: the first scenario (S0o) is the reference (baseline scenario) with observed climate, while S0 m is generated with simulated climate data of RCMs. S1 represents the impact of urbanization on the current model, while S0M45 and S0M85 represent climate change impacts under RCP4.5 and RCP8.5, respectively. The effects of combined urbanization and climate change are assessed with S1M45 and S1M85 under RCP4.5 and RCP8.5, respectively. S1 outputs are compared to S0o, while the scenario involving climate change will be compared to S0 m. This is to ensure the comparison between the two scenarios is consistent with the climate data frame (observations or simulations).
The local/watershed-scale approach consisted of evaluating the impact of each scenario at the main outlet (watershed scale) and the outlet of sub-watershed 46 (local scale). This upstream sub-watershed was selected because it is the most affected by urbanization.
Monthly averages are calculated by averaging the values for a specific month during the simulation period excluding the warmup period, as the model does not include it in the output time series. Annual maximum boxplots are created based on the maximum discharge observed at each year of simulation.
RESULTS AND DISCUSSION
Parameter sensitivity
A total of 55 sub-watersheds were delineated (Figure 4).
Parameter selection and sensitivity analysis
Recommendations of similar studies in the literature (e.g., Abbaspour 2007; Arnold et al. 2012; Abbaspour et al. 2015; El-Khoury et al. 2015; Khoiab & Thomb 2015) led to the selection of 22 parameters to use in the calibration. First, a sensitivity analysis was performed to assess the influence of these parameters on the objective function, which in this case is the sum of the NS coefficient for streamflow, total nitrogen, and total phosphorus in their respective calibration periods. The sensitivity analysis provides information about the most important process drivers in the study region, depending on local characteristics (Abbaspour et al. 2018). The results are shown in Table 11, where it can be observed that the most sensitive parameter is N_UPDIS with a p-value and t-stat of 0.00 and 35.17, respectively. On the other hand, GWQMN.gw is the parameter with the least influence on the performance of the model. In this case, the uptake distribution of N seemed to affect the model performance the most.
Parameter name . | t-stat . | p-value . | Rank . |
---|---|---|---|
A__N_UPDIS.bsn | 35.17 | 0.00 | 1 (most sensitive) |
A__ERORGP.hru | −23.95 | 0.00 | 2 |
A__ERORGN.hru | 19.98 | 0.00 | 3 |
R__CN2.mgt | −9.30 | 0.00 | 4 |
A__OV_N.hru | 6.37 | 0.00 | 5 |
V__ESCO.hru | −2.04 | 0.04 | 6 |
A__RCN.bsn | −1.37 | 0.17 | 7 |
V__GW_REVAP.gw | −1.29 | 0.20 | 8 |
A__REVAPMN.gw | 0.93 | 0.35 | 9 |
A__P_UPDIS.bsn | 0.85 | 0.40 | 10 |
V__SOL_K(..).sol | −0.78 | 0.43 | 11 |
A__PPERCO.bsn | −0.77 | 0.44 | 12 |
A__SLSUBBSN.hru | 0.65 | 0.52 | 13 |
R__GWSOLP.gw | 0.64 | 0.52 | 14 |
V__SOL_AWC.sol | 0.48 | 0.63 | 15 |
A__PHOSKD.bsn | 0.47 | 0.64 | 16 |
A__PSP.bsn | 0.47 | 0.64 | 17 |
A__ANION_EXCL.sol | −0.43 | 0.67 | 18 |
R__SHALLST_N.gw | −0.30 | 0.76 | 19 |
V__ALPHA_BF.gw | −0.26 | 0.79 | 20 |
A__NPERCO.bsn | 0.15 | 0.88 | 21 |
A__GWQMN.gw | −0.03 | 0.98 | 22 (least sensitive) |
Parameter name . | t-stat . | p-value . | Rank . |
---|---|---|---|
A__N_UPDIS.bsn | 35.17 | 0.00 | 1 (most sensitive) |
A__ERORGP.hru | −23.95 | 0.00 | 2 |
A__ERORGN.hru | 19.98 | 0.00 | 3 |
R__CN2.mgt | −9.30 | 0.00 | 4 |
A__OV_N.hru | 6.37 | 0.00 | 5 |
V__ESCO.hru | −2.04 | 0.04 | 6 |
A__RCN.bsn | −1.37 | 0.17 | 7 |
V__GW_REVAP.gw | −1.29 | 0.20 | 8 |
A__REVAPMN.gw | 0.93 | 0.35 | 9 |
A__P_UPDIS.bsn | 0.85 | 0.40 | 10 |
V__SOL_K(..).sol | −0.78 | 0.43 | 11 |
A__PPERCO.bsn | −0.77 | 0.44 | 12 |
A__SLSUBBSN.hru | 0.65 | 0.52 | 13 |
R__GWSOLP.gw | 0.64 | 0.52 | 14 |
V__SOL_AWC.sol | 0.48 | 0.63 | 15 |
A__PHOSKD.bsn | 0.47 | 0.64 | 16 |
A__PSP.bsn | 0.47 | 0.64 | 17 |
A__ANION_EXCL.sol | −0.43 | 0.67 | 18 |
R__SHALLST_N.gw | −0.30 | 0.76 | 19 |
V__ALPHA_BF.gw | −0.26 | 0.79 | 20 |
A__NPERCO.bsn | 0.15 | 0.88 | 21 |
A__GWQMN.gw | −0.03 | 0.98 | 22 (least sensitive) |
Model performance
The calibration was performed using the 22 parameters identified in the previous section. The optimal range for each parameter as well as the fitted values are shown in Table 12. The outputs of the calibrated model are shown in Figure 5.
Parameters . | Optimal range . | Fitted values . |
---|---|---|
r__CN2.mgt | [0.11713, 0.12033] | 0.11873 |
a__GWQMN.gw | [−1347.63, −1136.19] | −1241.9 |
v__GW_REVAP.gw | [0.097516, 0.100094] | 0.098805 |
v__ESCO.hru | [0.997944, 1.010913] | 1.0044 |
v__ALPHA_BF.gw | [0.482128, 0.54378] | 0.51295 |
v__SOL_AWC().sol | [0.974373, 1.001545] | 0.98796 |
a__REVAPMN.gw | [187.1984, 212.1717] | 199.69 |
v__SOL_K().sol | [0.577949, 0.600627] | 0.58929 |
a__OV_N.hru | [7.371342, 8.932384] | 8.1519 |
a__SLSUBBSN.hru | [135.1207, 137.2007] | 136.16 |
a__RCN.bsn | [14.5449, 16.35953] | 15.452 |
a__N_UPDIS.bsn | [−20.5479, −16.0446] | −18.296 |
r__SHALLST_N.gw | [−135.764, −86.191] | −110.98 |
a__ERORGN.hru | [−0.13804, 0.148926] | 0.005445 |
a__NPERCO.bsn | 1.144399, 1.161047] | 1.1527 |
a__ANION_EXCL.sol | [0.763784, 0.786918] | 0.77535 |
a__PSP.bsn | [0.570577, 0.594221] | 0.5824 |
a__PHOSKD.bsn | [116.2567, 120.4304] | 118.34 |
a__P_UPDIS.bsn | [89.41322, 99.38607] | 94.400 |
a__PPERCO.bsn | [12.34229, 12.44181] | 12.392 |
a__ERORGP.hru | [0.011739, 0.019001] | 0.01537 |
r__GWSOLP.gw | [297.5318, 319.4124 | 308.47 |
Parameters . | Optimal range . | Fitted values . |
---|---|---|
r__CN2.mgt | [0.11713, 0.12033] | 0.11873 |
a__GWQMN.gw | [−1347.63, −1136.19] | −1241.9 |
v__GW_REVAP.gw | [0.097516, 0.100094] | 0.098805 |
v__ESCO.hru | [0.997944, 1.010913] | 1.0044 |
v__ALPHA_BF.gw | [0.482128, 0.54378] | 0.51295 |
v__SOL_AWC().sol | [0.974373, 1.001545] | 0.98796 |
a__REVAPMN.gw | [187.1984, 212.1717] | 199.69 |
v__SOL_K().sol | [0.577949, 0.600627] | 0.58929 |
a__OV_N.hru | [7.371342, 8.932384] | 8.1519 |
a__SLSUBBSN.hru | [135.1207, 137.2007] | 136.16 |
a__RCN.bsn | [14.5449, 16.35953] | 15.452 |
a__N_UPDIS.bsn | [−20.5479, −16.0446] | −18.296 |
r__SHALLST_N.gw | [−135.764, −86.191] | −110.98 |
a__ERORGN.hru | [−0.13804, 0.148926] | 0.005445 |
a__NPERCO.bsn | 1.144399, 1.161047] | 1.1527 |
a__ANION_EXCL.sol | [0.763784, 0.786918] | 0.77535 |
a__PSP.bsn | [0.570577, 0.594221] | 0.5824 |
a__PHOSKD.bsn | [116.2567, 120.4304] | 118.34 |
a__P_UPDIS.bsn | [89.41322, 99.38607] | 94.400 |
a__PPERCO.bsn | [12.34229, 12.44181] | 12.392 |
a__ERORGP.hru | [0.011739, 0.019001] | 0.01537 |
r__GWSOLP.gw | [297.5318, 319.4124 | 308.47 |
Scenarios . | Impact . | Annual maximum flow . | Annual mean flow . | Annual nitrogen load . | Annual phosphorus load . | ||||
---|---|---|---|---|---|---|---|---|---|
Value (m3/s) . | Change (%) . | Value (m3/s) . | Change (%) . | Value (×105 Kg) . | Change (%) . | Value (×103 Kg) . | Change (%) . | ||
S0o | Watershed scale | 10.75 | – | 35.40 | – | 8.71 | – | 6.62 | – |
Local | 0.57 | – | 1.97 | – | 0.19 | – | 0.62 | – | |
S0 m | Watershed scale | 12.63 | – | 3.15 | – | 8.96 | – | 6.65 | – |
Local | 0.68 | – | 0.18 | – | 0.19 | – | 0.62 | – | |
S1 | Watershed scale | 10.83 | 0.73 | 3.09 | 1.57 | 8.55 | − 1.88 | 8.38 | 26.49 |
Local | 0.60 | 4.78 | 0.18 | 9.45 | 0.14 | − 29.00 | 1.08 | 73.56 | |
S0M45 | Watershed scale | 13.11 | 3.76 | 3.32 | 5.49 | 11.61 | 29.62 | 6.72 | 1.07 |
Local | 0.67 | − 1.05 | 0.18 | 2.61 | 0.24 | 28.24 | 0.62 | − 0.55 | |
S0M85 | Watershed scale | 12.75 | 0.90 | 3.38 | 7.52 | 9.14 | 2.03 | 6.35 | − 4.49 |
Local | 0.66 | − 3.54 | 0.18 | 1.13 | 0.19 | − 0.96 | 0.55 | − 11.19 | |
S1M45 | Watershed scale | 13.13 | 3.97 | 3.36 | 6.75 | 11.19 | 24.84 | 8.24 | 23.81 |
Local | 0.69 | 1.61 | 0.19 | 11.10 | 0.16 | − 12.92 | 1.03 | 66.15 | |
S1M85 | Watershed scale | 12.79 | 1.29 | 3.44 | 9.34 | 8.85 | − 1.20 | 7.93 | 19.15 |
Local | 0.68 | − 0.36 | 0.20 | 13.07 | 0.14 | − 24.80 | 0.98 | 57.97 |
Scenarios . | Impact . | Annual maximum flow . | Annual mean flow . | Annual nitrogen load . | Annual phosphorus load . | ||||
---|---|---|---|---|---|---|---|---|---|
Value (m3/s) . | Change (%) . | Value (m3/s) . | Change (%) . | Value (×105 Kg) . | Change (%) . | Value (×103 Kg) . | Change (%) . | ||
S0o | Watershed scale | 10.75 | – | 35.40 | – | 8.71 | – | 6.62 | – |
Local | 0.57 | – | 1.97 | – | 0.19 | – | 0.62 | – | |
S0 m | Watershed scale | 12.63 | – | 3.15 | – | 8.96 | – | 6.65 | – |
Local | 0.68 | – | 0.18 | – | 0.19 | – | 0.62 | – | |
S1 | Watershed scale | 10.83 | 0.73 | 3.09 | 1.57 | 8.55 | − 1.88 | 8.38 | 26.49 |
Local | 0.60 | 4.78 | 0.18 | 9.45 | 0.14 | − 29.00 | 1.08 | 73.56 | |
S0M45 | Watershed scale | 13.11 | 3.76 | 3.32 | 5.49 | 11.61 | 29.62 | 6.72 | 1.07 |
Local | 0.67 | − 1.05 | 0.18 | 2.61 | 0.24 | 28.24 | 0.62 | − 0.55 | |
S0M85 | Watershed scale | 12.75 | 0.90 | 3.38 | 7.52 | 9.14 | 2.03 | 6.35 | − 4.49 |
Local | 0.66 | − 3.54 | 0.18 | 1.13 | 0.19 | − 0.96 | 0.55 | − 11.19 | |
S1M45 | Watershed scale | 13.13 | 3.97 | 3.36 | 6.75 | 11.19 | 24.84 | 8.24 | 23.81 |
Local | 0.69 | 1.61 | 0.19 | 11.10 | 0.16 | − 12.92 | 1.03 | 66.15 | |
S1M85 | Watershed scale | 12.79 | 1.29 | 3.44 | 9.34 | 8.85 | − 1.20 | 7.93 | 19.15 |
Local | 0.68 | − 0.36 | 0.20 | 13.07 | 0.14 | − 24.80 | 0.98 | 57.97 |
Overall, the model shows a good performance for the discharge at outlet 8 in calibration: NS is in the satisfactory range (0.5<NS=0.67<0.7) according to Moriasi et al. (2007), while the PBIAS is excellent (close to 0). In validation, a relatively satisfactory performance was reached with an NS of 0.48 and a PBIAS of +16.9% (Figure 5(a)). However, the graphical comparison shows the model tends to underestimate the peaks and overestimate low flows. This is a common problem in rainfall–runoff modeling, which is linked to the quality of climatic inputs and the spatial discretization of the study area (Ghosh & Hellweger 2012; Sharma & Regonda 2021). In addition, the spatial repartition of rain gauges is a factor (Lobligeois et al. 2014) as a sparse number of rain gauges (as in this paper) or very large sub-watershed may lead to some extreme precipitation events being overlooked or extreme occurring in one small area that happens to have a rain gauge assumed to occur over a larger area.
As expected, a lower performance was obtained for water quality parameters, especially for P. According to the NS (0.54 for outlet 8 and 0.55 for outlet 30), the simulation of the N at both outlets showed acceptable results in calibration. Meanwhile, both PBIAS were excellent: 5.2 and −9.9%. Unsurprisingly, discharge is better simulated than nutrients by the model. A few studies (El-Khoury et al. 2015; Zhou et al. 2015; Čerkasova et al. 2018) have also obtained better performance for streamflow calibration. The reason for the low performance in predicting water quality parameters is the limited available data for a large number of parameters, highlighted in Hoang et al. (2019) and Noori et al. (2020). There are also greater uncertainties in nutrient data associated with errors in streamflow measurements, sample collection, storage, and analysis (Harmel & Smith 2007). For instance, in this study, the length of validation period for discharge was 11 years (2007–2015), while only 6 and 3 years of data were used for N and P loads at outlets 8 and 30, respectively.
All the simulations for water quality parameters showed unsatisfactory performance in the validation phase except for P at outlet 30, where PBIAS was 16.5%. The graphical comparison shows the model better simulates the average value of N at outlet 8 (Figure 5(b)), but not the peaks. It is important to note that these results are in the same range (or even better) as some published studies such as El-Khoury et al. (2015) with NS=0.24 and 0.28 in the calibration of nitrite-nitrogen and organic P, respectively; Shrestha et al. (2021) with NS=−0.71 and PBIAS=60.3% for P calibration.
Impacts of the scenarios on discharge, N load, and P load
This section presents the impacts of each scenario in terms of annual maximum discharge and monthly average discharge, N load, and P load at the main outlet (global) and at the outlet of sub-watershed 46 (local). This sub-watershed, located upstream (Figure 4), was selected because it was one of the most affected by the developments, with an important difference in the urbanization ratio compared to the global watershed. This allows an observable difference in the impacts at the two scales when comparing them. It also helps to investigate how different the effects between the global and the local scales would be. Figure 6 illustrates the repartition of the land uses before and after the developments in terms of watershed and local scales.
Impact of urbanization (S0o vs. S1) on discharge, N load, and P load
Figure 7 shows the annual maximum discharge, while Figures 8 and 9 present the average monthly outputs related to discharge, N, and P for the experiments S0o and S1. Globally, urbanization slightly increased the averages of discharge and P load by 1.9 and 23%, respectively, while the average N load decreased by 4.4%. The N load decreases particularly in July–August, where higher precipitation amounts are observed. In terms of local impacts, the variation of the monthly averages has the same direction of change for all the three variables except that local impact is more significant. Annual averages in local increased by 9.6 and 69.5% for discharge and P load, respectively, and decreased by 30.2 for N load (higher than the previous percent changes). The high contrast change in P load compared to the N load is not surprising, since the model performance in simulating P was low (Figure 5(c) and 5(e)).
Peak values for discharge and P are observed during the summer–spring period (April–August), while the same is observed in July–August for N. This is due to the increase of total precipitations (combination of snow melting and rainfall). In terms of water quantity, extreme events are more likely to occur around April. They are expected to vary by +4.78% at the local scale, while a relatively small difference of +0.73% is estimated at the watershed scale (Figure 7).
Impacts of climate change (S0 m vs. S0M45/85) on discharge, N load, and P load
Unlike the land-use change scenario, the change in climate change scenario is calculated and compared to S0 m output instead of S0o, although the S0o results are also presented on graphs only as a reference.
In terms of annual maximum, the discharge has experienced a relatively low change under the two RCPs, as presented in Figure 10. On average, at the watershed scale, the maximum discharge reached 12.63 m3/s in S0 m, and the percentages of change under RCP4.5 and RCP8.5 scenarios are +3.76 and +0.90%, respectively. However, an opposite change is observed at the local scale with −1.06 and −3.54% for S0M45 and S0M85, respectively. Annual average discharge presented a similar trend in S0M45 at the watershed scale and the local scale with an increase of 5.49 and 2.61%, respectively, while in S0M85, the change was 7.52% (watershed scale) and 1.13% (local) increase. In terms of water quality, overall the climate change seems to have a significant impact on N in S0M45 with an increase of 29.62% (+265 tons) globally and 28.24% (+5 tons) locally, while the change is less important in S0M85 with +2.03% (+18 tons) and −0.96% (−0.18 tons) globally and locally, respectively. Monthly results are illustrated in Figure 11 and Figure 12.
Combined impacts of land-use and climate changes (S0 m vs. S1M45/85) on discharge, N load, and P load
Separately assessing the impacts of land-use change and climate change provides a better understanding of the level of impact caused by each of the two scenarios on the water quantity and quality. However, it is also important to determine and quantify the impact when the two effects are coupled. Figure 13 shows that the maximum discharge is not significantly affected by the combined effects of climate and land-use changes. The difference between the means (cross mark inside the boxes) of the scenarios S0 m, S1M45, and S1M85 is relatively small, both in watershed and local scales, although the RCP4.5 has a higher impact compared to the RCP8.5.
In contrast, the annual discharge is more sensitive to the RCP8.5 than RCP4.5. The annual average changed from 3.15 m3/s in the baseline conditions to 3.44 m3/s (+9.34%) under RCP8.5, while RCP4.5 caused an increase of 0.21 m3/s (+6.75%). This is also the same at the local scale, where an increase of 13.07 and 11.10% is observed under RCP8.5 and RCP4.5, respectively. The coupled effect caused a reduction of the annual N load except in S1M45, where a significant increase of 24.84% was observed at the watershed scale. The local impact under RCP8.5 is also likely to be important, as the load decreased by 24.80% compared to scenario S0 m. However, the results indicated the annual load at the main outlet under RCP8.5 was not significantly impacted, and only a 1.20% decrease (−10.8 tons) of the baseline condition load is expected. The quantity of P faced a watershed-scale significant increase and was even higher at the local outlet (in terms of percentage of change). The load increase at the main outlet was estimated at 1.6 tons (23.81%) in S1M45 and 1.3 tons in S1M85 (19.15%), while in the upstream sub-watershed, it was 0.411 tons (66.15%) and 0.360 tons (57.97%). Monthly results are illustrated in Figure 14 and Figure 15.
Summary of global and local impacts
Discussion
Results suggest that water quantity was more sensitive to climate change than urbanization in the Carp River watershed when the factors are taken individually (S1 and S0M45/85). When climate change alone is considered in the watershed, the increase in maximum discharge is expected to be within 0.90–3.76% at the main outlet. Urbanization will only lead to a 0.76% increase in the maximum discharge. Similarly, urbanization is expected to increase the average discharge by 1.57%, while this change is expected to be between 5.49 and 7.42% under climate change. This is not surprising as similar studies on areas located in Canada have come up with the same conclusion: Kaykhosravi et al. (2020) in Montreal and Toronto; El-Khoury et al. (2015) for the South Nation watershed in Ontario. For instance, the impact caused by the land-use and climate changes in the South Nation watershed was, respectively, +1.2 and +11.2%. Although the percentage of the urban area increased in the urbanization scenario, the agricultural lands are expected to remain dominant with 49.7% of the total area against only 13.5% for urban zones in 2050. This land-use repartition minimized the upstream developments’ effect and could explain the small impact observed in the Carp River watershed. A different urbanization rate or a different location of the developments within the same watershed may lead to very different results.
Monthly output analysis indicated extreme events were likely to occur at the end of the winter (around April) and during the summer. This is very common in cold areas. This study has shown that snow contribution to the total precipitations is more than rainfall, as the peaks of April are higher than during the summer, meaning that floods were expected to be more severe in April. According to the model, this characteristic of the watershed is not expected to change by 2050 because the predicted monthly hydrograph in the future conditions (S1M45/85) follows the same variation as the historical period (1990–2018). The variation of the peak flow during April is expected to be within −27.1 to 2.38%. Therefore, no significant increase is expected, but on the contrary, peaks may decrease up to 27.1% on average over the next 30 years. Although the peaks of July (due to rainfall events) are less severe, unlike April peaks, they were predicted to increase on average by 8.60%.
In terms of water quality, the sources of pollution of N and P are generally related to land uses such as urban and agricultural. Based on the impacts observed in the urbanization scenario, N seemed more sensitive to agricultural lands during July–August when agricultural activities occur. As an important amount of fertilizers is used for these activities, it is unsurprising to observe a drastic increase in the monthly graphs for all the scenarios. The sensitivity is also demonstrated by a 1.88% decrease in the annual total load at the main outlet when the ratio of the agricultural land is reduced by 2.63% (see Table 8). Similarly, a more significant decrease was observed at the local scale when the agricultural land ratio dropped from 49.7 to 10.6%, resulting in a reduction of the impact from 1.88 to 29.00%.
In the same perspective, as the increase of urban areas led to an increase of the annual P load (Figure 16(a)), the pollution of P seems to be dominated by the urban areas. This is consistent with what is observed at the local scale where a higher urban area ratio caused a larger increase (from 26.49% globally to 73.56% locally when urban ratio passed from 13.5 to 47.9%, Figure 6).
Considering the most severe conditions in the future (S1M85), a decrease in the annual average of N load (−1.20%) is expected in opposition to an increase of 24.84% in the scenario S1M45. Therefore, a contrasting impact on N is likely to occur under the two RCPs. As a justification, the increase of the temperatures in the future conditions, especially in RCP8.5, influences the evaporation and transpiration (indispensable for the growth of plants) from natural (forests and water bodies) and agricultural lands. Given that plant cells control the openings where water is released to the atmosphere, higher temperatures stimulate them to open, enhancing water release. Consequently, the rate of transpiration, defined as the water loss from living plant surfaces (Manashi 2016), goes up with the temperature (USGS 2020). Plant growth is influenced by the availability of N, a vital element in the photosynthesis process. The higher the plant growth rate, the higher the demand in N. Considering the expected increase in transpiration and N availability, the plant's demand for N is likely to increase in the future. Consequently, this will decrease the annual N load, as observed in the S1M85 scenario with a −1.20% reduction.
The above justification is also valid for P, as the increase of the load at the main outlet is smaller in the hottest scenario (S1M85) than that of the coolest (S1M45) with 19.15 and 23.81%, respectively. The difference is that P is less soluble than N, and consequently, it is less sensitive to temperature. N and P are both partially lost in plant uptake, but in N's case, the losses also include denitrification and volatilization. Given that higher temperatures enhance their decomposition and accumulation (Hong & Tang 2014; Geng et al. 2017), denitrification and volatilization are promoted, and consequently, N loss may be more significant, which is the case in the Carp River watershed. In addition, even though having similar sources of pollution, N and P seemed to have different transport pathways. Due to its solubility, N is more likely to be retained within the watershed by infiltration into the groundwater and/or in the atmosphere by volatilization, while the P will mostly go into the surface water, especially for urbanized areas. This was demonstrated for the Mississippi River watershed in the USA (Hobbie et al. 2017) where only 22% of P inputs was retained vs. 80% of the N inputs. Nieder et al. (2018) also stated that nitrate is the most common chemical contaminant in the world's groundwater aquifers, and surface water is particularly affected by P, confirming the findings in the Carp River watershed.
The contrast in N and P variation (Figure 16) observed in our study area differs from some published studies such as Mehdi & Lehner (2014), El-Khoury et al. (2015), Eum et al. (2016), and Luo et al. (2020) and with the different direction of change. For instance, El-Khoury et al. (2015) found out that P will increase regardless of the type of conditions (LUC, CC, or the two combined), which contradicts what is observed in the Carp River watershed. This divergence indicates the impact of climate change, land-use change, or both must be evaluated case by case. It cannot be generalized even for watersheds located in the same region because every study area has different characteristics in terms of topography, land-use occupation, and climate conditions. Therefore, water quality is considerably affected by the entire watershed characteristics (Liu et al. 2018).
CONCLUSION
This study examined the influence of the land-use change and climate change individually and their coupled impact on the watershed of the Carp River using a hydrological model SWAT. The model was applied using input data from several sources and different resolutions, and overall, the calibrated model showed a satisfactory performance. The impacts were evaluated by developing realistic scenarios based on the 2050 year land-use map projected by the local authority and future climate conditions under RCP4.5 and RCP8.5. Results showed that climate change and urbanization's impacts vary greatly depending on the spatial scale and geographic location. In addition to that, the following points list the main other findings that emerged from this study: (a) climate change is likely to be the most dominant factor affecting the flow (peaks and trend) and N, while urbanization will control the P quantity in the future; (b) global and local impacts may significantly differ and it is crucial to check the impacts at both levels; (c) the local impact can be more significant than the watershed-scale impact; (d) impacts are not additive; and (e) for impacts assessment study, evaluating one effect alone without considering the other is incorrect and leads to a wrong understanding of the variation. In terms of quantity and quality:
Expected global increase of the annual maximum flow from 1.29 to 3.97%.
Expected increase of the annual average flow (global) from 6.75 to 9.34%.
Expected global variation of the annual average N load from −1.20 to 24.84%.
Expected global variation of the annual average P load from 19.15 to 23.81%.
As a recommendation, P management should emphasize reducing watershed inputs and transport from urban zones. In contrast, N management should reduce watershed inputs and transport from cultivated zones when they are active.
ACKNOWLEDGEMENTS
The authors gratefully acknowledge the data and financial support of the MVCA and financial support from MITACs. The authors also expressed their gratitude to the City of Ottawa for providing us with high-resolution topographic data and for sharing with us the future urban development map of the municipality.
AUTHOR'S CONTRIBUTION
B.-S.Z. is involved in conceptualization, data curation, formal analysis, methodology, software, visualization, writing – original draft, and writing – review & editing. O.S. is involved in conceptualization, supervision, resources, visualization, and writing – review & editing. M.S. is involved in conceptualization, supervision, resources, visualization, writing – review & editing, and funding acquisition. N.N. and K.S. are involved in investigation, data curation, and project administration.
FUNDING
This work was supported by the MVCA through the Mitacs Accelerate (No. IT15546) program.
DATA AVAILABILITY STATEMENT
All relevant data are included in the paper or its Supplementary Information.