Downscaling Climatic Variables at a River Basin Scale: Statistical Validation and Ensemble Projection under Climate Change Scenarios

Abstract

Climatic statistical downscaling in arid and topographically complex river basins remains relatively lacking. To address this gap, climatic variables derived from a global climate model (GCM) ensemble were downscaled from a grid resolution of 2.5° × 2.5° down to the station level. For this purpose, a combination of multiple linear and logistic regressions was developed, calibrated and validated with regard to their predictions of monthly precipitation and daily temperature in the Jordan River Basin. Seasonal standardized predictors were selected using a backward stepwise regression. The validated models were used to examine future scenarios based on GCM simulations under two Representative Concentration Pathways (RCP4.5 and RCP8.5) for the period 2006–2050. The results showed a cumulative near-surface air temperature increase of 1.54 °C and 2.11 °C and a cumulative precipitation decrease of 100 mm and 135 mm under the RCP4.5 and RCP8.5, respectively, by 2050. This pattern will inevitably add stress to water resources, increasing management challenges in the semi-arid to arid regions of the basin. Moreover, the current application highlights the potential of adopting regression-based models to downscale GCM predictions and inform future water resources management in poorly monitored arid regions at the river basin scale.

1. Introduction

The temporal increase in surface and air temperatures coupled with alterations in precipitation patterns due to climate change has affected the water balance and reduced the available water resources in various regions [1,2]. These hydro-climatic shifts are having profound multiscale impacts on agriculture and food security, thus highlighting the need to better understand the impacts of global climate change on regional and local water resources [3,4,5]. Proper water resources management is imperative and requires knowledge of precipitation, temperature, air humidity, wind speed and other variables and how they vary under a changing climate. Thus, future projections of climatic variability are needed to improve water resources management through informed mitigation measures and adaptation strategies.

Global climate models (GCMs) are now the main resource for obtaining these future projections over a variety of regional and temporal scales due to their ability to capture the many physical processes underpinning climate systems [6]. Nevertheless, GCM outputs are still prone to biases and their ability to capture subgrid-scale characteristics is restricted. Additionally, their ability to provide physical atmospheric dynamics and hydrological processes at a regional scale is limited [7]. While general circulation models (GCMs) are considered the primary tool for projecting changes in global climatic conditions, they provide a coarse spatial resolution that often falls short of representing changes observed in the local climate [8,9,10]. Accordingly, when a GCM output is adopted for the assessment of climate change impacts, large biases are often generated in the simulated hydrological processes at the local (~typical watershed) level [11,12]. This necessitates the downscaling of these outputs to higher resolutions to ensure a reliable assessment of the local hydrological impacts of climate change.

Two approaches are commonly used in the downscaling process namely, dynamical and statistical [13]. In dynamical downscaling (DD), regional climate models (RCMs) or mesoscale weather models are utilized to downscale the GCM output, allowing the inclusion of small-scale details and enhancing the reliability of the results [14,15]. In this context, an RCM is broadly similar to an atmospheric GCM but with higher resolution. RCMs integrate the complex topography, the land-sea contrast, finer surface heterogeneities, and more detailed descriptions of physical processes with the large-scale information supplied by the GCM output to generate climate information at higher resolution [16]. Since the RCM is nested within a GCM and the boundary conditions needed for the RCM to downscale large-scale conditions are provided by the GCM, the overall quality of the dynamically downscaled RCM output depends on the accuracy of the large-scale GCM and its biases [17]. The main drawbacks of dynamical downscaling are the high complexity and computational cost [18]. Lately, some studies have combined the traditional downscaling approach with machine learning so that a collection of statistical models can emulate the downscaling [19].

Statistical downscaling (SD) develops empirical relationships between larger-scale, and local-scale observed weather variables based on statistical methods that link between large surface predictors and local surface predictands [20,21,22,23,24]. In SD, either large-scale GCM outputs or RCM outputs are used as predictors to obtain local variables or predictands. The development of statistical relationships between local and large scales may implicitly include topography, vegetation and hydrological processes [25]. SD methods are computationally inexpensive and usually require less time and effort in comparison to DD. However, SD relies on the critical assumption of stationarity, which assumes that the relationship between the coarse GCM-simulated climate and the local climate will remain valid in the future [1,26]. Another limitation of the SD method is its dependency on long-term time series of observed data to be able to derive statistical relationships [27]. In addition, SD methods might poorly assess the variance and might be inaccurate in reproducing extreme events [28]. Furthermore, and similar to dynamical downscaling, SD results are dependent on the GCM output and their biases. Lately, deep learning techniques (in particular convolutional neural networks, CNNs) have recently emerged as a promising approach for statistical downscaling due to their ability to learn spatial features from huge spatiotemporal datasets.

Nevertheless, because each SD approach relies on certain assumptions and approximations, the findings are frequently compromised with bias and restrictions [29,30]. Certain presumptions raise questions about the validity of downscaled forecasts and might restrict the range of applications for which downscaling techniques are appropriate [31]. The assumptions that underlie the conclusions for various approaches need to be assessed because there is not a single statistical downscaling technique that works best for all applications and geographical areas, even while certain methods are better for particular uses. As a result, end users may choose the best approach for each application depending on the advantages and disadvantages of the approach, the information demands (such as the required geographical and temporal resolutions), and the resources that are available (data, knowledge, computer resources, etc.).

However, the climate community often views these models as opaque, especially when it comes to applications related to climate change [32]. Despite these shortcomings, SD remains a valuable and widely used tool. SD is often categorized into three types: weather classification, regression, and weather generators [13] (

Table 1. Statistical downscaling methods.

Method	Description	Drawbacks	Examples
Weather classification	- Arrange days into discrete weather states using clustering techniques [13] - Relationships between large-scale categories and local climate variables can be identified and predicted using regression, resampling or Monte Carlo techniques	Subjectivity in creating classification states	Principal components, neural networks such as radial basis function (RBF), multilayer perceptron (MLP), analog and fuzzy c-mean clustering
Regression	- Most common approach that represents the relations, either linear or non-linear, between predictands and large-scale atmospheric forcing or GCM that are the predictors [33]	Variance underestimation, especially of daily precipitation, because of the non-normality of the process	Examples: multiple linear regression (MLR), positive coefficient regression (PCR), principal component regression (PCR), stepwise regression (SR), and canonical correlation analysis (CCA)
Weather generators	- Produce a series of data that fit to observed data. The models produced are linked to random number generation algorithms to yield stochastic realizations of daily weather series that resemble real weather data - Offer the ability to generate data of any time length with similar properties as those characterizing observations [34]	- Requires long sequences of daily data and is sensitive to missing data [35]. - Replicate the mean and variance but not the actual sequence of events [36] - Miss the temporal autocorrelation of real weather	Examples: K-nearest neighbor (KNN), Markov chains, conditional random fields (CRF) and Gamma distribution are examples of weather generator methods

). Each of these approaches makes use of various downscaled parameters, GCMs, methods to select predictors, and performance measures (Tables S1 and S2, Supplemental Material).

Among the SD methods, MLR has been widely used in climate change assessment studies [37,38,39,40,41,42] and is reportedly the most capable of reproducing various statistical characteristics of the observed data [43,44]. Several open questions, however, remain concerning the application and assessment of SD at a river basin scale with complex topology. To bridge this gap, in this study, we address the following questions: can MLR downscaling significantly improve the projected climate variables at a local scale compared to the coarse GCM fields? How should an MLR model trained on historic reanalysis data be applied to downscale future GCM projections and deal with their biases? The application domain for this study is the Jordan River Basin (JRB), with its complex topography. A question that naturally follows is thus: what is the vulnerability of water resources in the arid JRB basin to climate change? The overarching goal of this study is to decrease the uncertainties in projecting climate variability toward a better definition of strategies that can meet the challenge of water scarcity in a region where climate change is expected to exacerbate existing water shortages [45,46].

2. Materials and Methods

2.1. Study Area and Data

The study area encompasses the JRB, which is a transboundary river shared among five riparians (Figure 1): Israel, Jordan, Lebanon, Palestinian Authority, and Syria, with a total catchment area of 18,300 km². It is classified into three sections: the upper part (in Israel, Lebanon, and Syria) encompassing the Hasbani, Dan, and Banias Rivers that flow into Lake Tiberias, the Yarmouk River that is fed by springs and wadis mostly in Syria, and the lower part including the Zarqa River with much wastewater from Amman [47,48] to eventually discharge into the Dead Sea.

The study area, which is characterized by a complex topography and several microclimates, is located within the climate zone of the eastern Mediterranean and is thus vulnerable to global climate change effects [49]. In general, the Middle East region is mostly comprised of arid and semi-arid lands [50], which is particularly affected by climate change and expected to witness further water shortages due to population growth coupled with a predicted increase in temperature and a decrease in precipitation [51,52,53]. The situation is exacerbated by a historical conflict that is reflected in disagreements over shared water resources in the basin due to scarcity and shortage of water.

We performed statistical downscaling of the National Centers for Environmental Prediction (NCEP) predictors (NCEP data were obtained from https://www.esrl.noaa.gov/ (accessed on 14 October 2018)) to obtain monthly precipitation and daily temperature scenarios at 41 observational stations in the basin and its immediate surroundings. The stations are depicted in Figure 2 and their details are summarized in

Table 2. Stations downscaled.

Station Name	Latitude (°)	Longitude (°)	Elevation (m)	Source	Variable	Data Availability	Country
Amman Hussein College	31.58	35.56	834	JMD	P	January 2000–March 2012	JO
Bal’ama	32.14	36.05	695	JMD	P	January 2000–March 2012	JO
Baqura Met. Station	32.61	35.60	−227	JMD	P	January 1981–April 2009	JO
Damascus International	33.41	36.52	616	NCDC	T	January 1981–December 2017	SYR
Deir Alla Agr. Station	32.12	35.36	−224	JMD	P	January 2000–March 2012	JO
En Nueiyime	32.25	35.55	748	JMD	P	January 1981–April 2009	JO
Ghor Safi	31.03	35.47	−350	NCDC	T	Jul 1983–December 2017	JO
H4 Airbase	32.54	38.20	686	NCDC	T	January 1981–December 2017	JO
Har Kenaan	32.97	35.50	934	NCDC	P, T	January 1981–December 2017	IS
Hosha	32.27	36.04	589	JMD	P	January 1981–April 2009	JO
Husn	32.29	35.53	637	JMD	P	January 1981–April 2009	JO
Irbid School	32.56	35.85	616	JMD	P	January 1981–April 2009	JO
Jaber Mughayyir	32.31	36.13	571	JMD	P	January 1981–April 2009	JO
Jarash	32.17-	35.54	585	JMD	P	January 2000–March 2012	JO
Jerusalem Central	31.77	35.22	815	NCDC	P, T	1981–2014/1981–1999	IS
Jubeiha	32.02	35.58	980	JMD	P	January 2000–March 2012	JO
K. H. Nursery Evap.St(Baq’a)	32.07	35.84	950	JMD	P	January 2000–March 2012	JO
Khanasira	32.24	36.03	810	JMD	P	January 1981–April 2009	JO
Kharja	32.40	35.53	441	JMD	P	January 1981–April 2009	JO
King Hussein	32.36	36.26	683	NCDC	T	January 1983–December 2017	JO
Kitta	32.17	35.51	665	JMD	P	January 2000–March 2012	JO
Kufr Saum	32.41	35.48	423	JMD	P	January 1981–April 2009	JO
Ma An	30.17	35.78	1069	NCDC	T	January 1981–December 2017	JO
Mafraq Airport	32.20	36.14	667	JMD	P	January 1981–April 2009	JO
Midwar	32.17	36.00	760	JMD	P	January 2000–March 2012	JO
Nawasif	32.08	36.16	590	JMD	P	January 2000–March 2012	JO
Prince Feisal Nursery	32.12	35.53	300	JMD	P	January 2000–March 2012	JO
Prince Hasan	32.16	37.15	677	NCDC	T	January 1981–December 2017	JO
Qafqafa	32.20	35.56	930	JMD	P	January 2000–March 2012	JO
Beirut Airport	33.82	35.49	27	NCDC	T	January 1981–December 2017	LB
Ramtha Boys School	32.34	36.01	513	JMD	P	January 1981–April 2009	JO
Rumeimin	32.06	35.48	675	JMD	P	January 2000–March 2012	JO
Ruseifa	32.01	36.02	655	JMD	P	January 2000–March 2012	JO
Sihan	32.08	35.46	495	JMD	P	January 2000–March 2012	JO
Subeihi	32.09	35.42	500	JMD	P	January 2000–March 2012	JO
Sukhna	32.08	36.04	500	JMD	P	January 2000–March 2012	JO
Turra	32.38	36.00	446	JMD	P	January 1981–April 2009	JO
Um El-Jumal Evap .St	32.32	36.37	680	JMD	P	January 2000–March 2012	JO
Um Jauza	32.06	35.44	860	JMD	P	January 1981–March 2012	JO
Um Qeis	32.39	35.41	351	JMD	P	January 1981–April 2009	JO
Wadi Dhuleil Nursery	32.08	36.17	575	JMD	P	January 2000–March 2012	JO

NCDC: National Climatic Data Center, JMD: Jordan Meteorology Department, T: Temperature, P: Precipitation IS: Israel, JO: Jordan, LB: Lebanon, SYR: Syria.

, including the corresponding source of observed data at these stations and the years when data were available. The observed data were obtained from the National Climatic Data Center (NCDC) (https://www.ncdc.noaa.gov/cdo-web/ (accessed on 14 January 2019) and the Jordan Meteorology Department (JMD) and used in the downscaling model, along with the NCEP predictors. The data were also used as a reference for the GCM ensemble bias correction, as detailed later. The NCEP predictors are outputs from atmospheric models that assimilate historic data from surface observation stations, upper-air stations, and satellite-observing platforms [54], yielding results representing those that could be expected to prevail in the actual historical climate record [55].

The GCM output was obtained from the Coupled Model Intercomparison Project phase 5 (CMIP5) archive (https://esgf-node.llnl.gov/ (accessed on 29 April 2019)) [56] for the period 1981–2005 to validate the CMIP data against the historical dataset and for the period 2006–2050 under two Representative Concentration Pathways (RCP): RCP4.5 and RCP8.5 for future projections. RCP4.5 reflects a stabilized scenario, where the total radiative forcing reaches a plateau before 2100 through reducing greenhouse gas emissions [57], whereas RCP8.5 is characterized by increasing greenhouse gas emissions over time [2].

The physical processes of the global climate system in the atmosphere, ocean, cryosphere, and land surface in response to shifting concentrations of greenhouse gases and aerosols are described by the GCMs, which are numerical models. Using three-dimensional global grids, GCMs provide physical and geographical estimations of regional climate and climate change. The CMIP5 makes use of the most recent generation of GCMs that are used in this study to understand both past and future climate changes. These are the models that served as the foundation for the IPCC’s Fifth Assessment Report (AR5) [2].

Several GCMs (

Table 3. GCMs used for downscaling.

Model Name	Institution	Atmospheric Grid Resolution		Scenario	Dates
		Latitude	Longitude
CanESM2	Canadian Centre for Climate Modeling and Analysis	2.7906°	2.8125°	Historical	1981–2005
				RCP4.5	2006–2050
				RCP8.5	2006–2050
GFDL-ESM2M	National Oceanic and Atmospheric Administration (NOAA) Geophysical Fluid Dynamics Laboratory	2.0225°	2.5°	Historical	1981–2005
				RCP4.5	2006–2050
				RCP8.5	2006–2050
HadGEM-CC	Met Office Hadley Centre	1.25°	1.875°	Historical	1981–2005
				RCP4.5	2006–2050
				RCP8.5	2006–2050

) were used to obtain an ensemble to ensure less dependency on one specific GCM [58]. The GCMs were selected based on the availability of data, their spatial resolution, and their individual performance in the region [59,60].

2.2. Re-Gridding and Standardization

As stated above, the study relied on the NCEP predictor datasets to build regression models using past observed data and an ensemble of GCMs under RCP4.5 and RCP8.5 to project precipitation and temperature into the future. Both NCEP and GCM predictors were used in model testing. The individual GCMs and the NCEP data differ in grid resolution and location; thus, all model outputs were re-gridded into the NCEP grid. For this purpose, the inverse distance weighted (IDW) method with a power of 2 was used; it assigns decreasing weights as the distances between locations increase [13].

The IDW is categorized as a deterministic technique and was developed by the National Weather Service in the United States in 1972. This is the result of the computation not having to satisfy any particular statistical assumption, which sets IDW apart from stochastic approaches (like Kriging) [61]. In this study, spatial data is interpolated using the IDW approach, which is based on the idea of distance weighting. By using the known data of sites that are close to the unknown location, it is possible to approximate the unknown spatial temperature and rainfall data. The IDW formulas are given as Equations (1) and (2) [62,63,64,65].

$$Z_p={\displaystyle \sum }_{i=1}^N{w}_i{Z}_i$$

(1)

$$w_i=\frac{d_i^{-\propto }}{{\sum }_{i=1}^N{d}_i^{-\propto }}$$

(2)

where Z_p means the unknown data; Z_i means the data of known stations; N means the number of stations; w_i means the weighting of each station; d_i means the distance between stations; α means the power, and is also a control parameter, generally set at two [66,67].

An ensemble GCM was generated from the re-gridded individual GCMs, and the SD was then applied to the ensemble mean results (since the downscaling model is linear, this is equivalent to downscaling each model and then taking the ensemble mean). Additionally, standardization for NCEP predictors and the GCM ensemble products was implemented to reduce systematic biases in the means and variances.

The systematic biases in the mean and variance of GCM predictor variables with respect to observations are eliminated by commonly used bias correction algorithms [13]. The process entails multiplying by the standard deviation and adding the mean of the corresponding observed or reanalysis data for a predefined baseline period at a timescale of interest after first standardizing the GCM-simulated variables by subtracting the mean and dividing by the standard deviation. Other methods of correcting bias that deal with the direct application of GCM data (particularly temperature and precipitation) include scaling, quantile matching, correction factors, and transfer functions [68,69,70,71,72,73].

The procedure of standardization converts NCEP and GCM predictors to Z-scores through subtraction of the mean and division by the standard deviation of the predictors/for a predefined baseline period for NCEP data. This procedure reduces the bias between the NCEP and the GCM data. However, a limitation of this standardization is that it assumes the bias is restricted to the mean and variance, while bias may also exist in other statistical parameters.

2.3. Selection of Predictors

The selection of significant predictors is a critical factor that could affect the accuracy of estimation. In this study, 12 predictors were selected from the NCEP and GCM data archive (

Table 4. NCEP predictors.

Predictor	Abbreviation
Temperature at 2 m	Temp2m
Pressure	Pressure
U wind component (East/West) at 500 pressure level	UWND.500
U wind component (East/West) at 1000 pressure level	UWND.1000
V wind component (North/South) at 500 pressure level	VWND.500
V wind component (North/South) at 1000 pressure level	VWND.1000
Relative humidity at 500 pressure level	RHUM.500
Relative humidity at 1000 pressure level	RHUM.1000
Specific humidity at 500 pressure level	SHUM.500
Specific humidity at 1000 pressure level	SHUM.1000
Geopotential height at 500 mb pressure level	HGT.500
Geopotential height at 850 mb pressure level	HGT.850

). Backward stepwise regression and correlation analysis were used to obtain potential predictors from the pool of predictors for each station separately and in the case of temperature for each season separately. Backward stepwise regression eliminates predictors that are least contributing to the model skill based on some criterion. The criterion selected for this study is the Akaike information criterion (AIC). A limitation of stepwise regression and correlation analysis is that including all variables may result in choosing predictors that might not have physical meaning; thus, for more accurate analysis, a pool of probable predictors was initially identified such that the selection is based on readily available data from NCEP and the GCM archive and having been used in past downscaling models (Table S2). In addition, the generated MLR models were checked for the potential problem of multi-collinearity between predictors. The statistic used to detect multi-collinearity was the variance inflation factor (VIF) as expressed in Equation (3) [74]:

$$VIF=\frac{1}{1-R_j^2}$$

(3)

where $R_j^2$ is the coefficient of determination when variable (predictor) $X_j$ is regressed on the remaining predictors. A variable is considered to be problematic if its VIF is larger than 10 [74]. For this study, multi-collinear variables with a VIF greater than 10 were removed from the model one at a time, starting with the one with the highest VIF.

2.4. Downscaling

Prior to downscaling, observed data and predictors were divided into a calibration and a validation period. Precipitation data were aggregated into monthly data, and the first step was to classify months as dry or wet, while the second step was to develop a regression model to calculate the precipitation when the month was classified as wet. A 1 mm threshold of precipitation per month was used for defining a dry/wet month [75,76]. For temperature, classification is not needed, and only the regression step is required. The proposed statistical downscaling method, therefore, includes logistic regression followed by multiple linear regression for precipitation and only multiple linear regression for the temperature.

2.4.1. Logistic Regression

A simple logistic regression technique was used to classify precipitation as dry or wet months. Logistic regression is a technique to deal with binary predictands, in this case, a month being dry or wet. It correlates large-scale predictors (from the NCEP or GCM dataset) with rain occurrence. Backward stepwise regression was used to select the most significant of the 12 predictors (Table 4). Equation (4) represents the logistic regression.

$$\ln \left(\frac{P}{1-P}\right)={\beta }_0+{\beta }_1{X}_1+{\beta }_2{X}_2+{\beta }_3{X}_3+\cdots +{\beta }_i{X}_i+\cdots +{\beta }_n{X}_n+\epsilon$$

(4)

where P is the probability of rain occurrence in a given month. It ranges between 0 and 1. ${\beta }_0$ is the intercept, ${\beta }_i$ is the coefficient of the ith independent variable or predictor $X_i$, and $\epsilon$ is the error or residuals of the data. If P is larger than the cut-off value, taken to be 0.5, the month is considered to be a rainy month, and if P is lower than 0.5, then there was no rain occurrence for that month. The performance of the logistic regression is here assessed by McFadden’s R² [77]. Note that precipitation was log-transformed in the MLR model to linearize the relationship between precipitation and the predictors and to eliminate the possibility of getting zeros.

2.4.2. Multiple Linear Regression

In this study, MLR was used to downscale the NCEP predictors to temperature and precipitation following the general form provided in Equation (5). The R software for statistical analysis was used to build the MLR downscaled models [78]. MLR is a least-squares-based method whereby the best-fit line is determined by minimizing the sum of squared errors between the linear model and the observed data. A good MLR explains most of the variance of the dependent variable with a minimum number of independent variables [79]. The performance of the models was evaluated by the coefficient of multiple determination (R²) (Equation (6)) and the correlation coefficient (R) (Equation (7)). The accuracy of the downscaled results was quantified in terms of the root mean squared error (RMSE) of the downscaled values relative to the observed ones (Equation (8)). R² represents how well the regression line approximates the real data, with higher R² representing a better fit and yielding a good predictive model with a low RMSE.

$$Y={\beta }_0+{\beta }_1{X}_1+{\beta }_2{X}_2+{\beta }_3{X}_3+\cdots +{\beta }_i{X}_i+\cdots +{\beta }_n{X}_n+\epsilon$$

(5)

$$R^2=\frac{{\left({\sum }_{i=1}^n\left(P_i-P\right)\left(O_i-O\right)\right)}^2}{{\sum }_{i=1}^n{\left(P_i-\overline{P}\right)}^2.{\sum }_{i=1}^n{\left(O_i-\overline{O}\right)}^2}$$

(6)

$$R=\sqrt{\frac{{\left({\sum }_{i=1}^n\left(P_i-P\right)\left(O_i-O\right)\right)}^2}{{\sum }_{i=1}^n{\left(P_i-\overline{P}\right)}^2.{\sum }_{i=1}^n{\left(O_i-\overline{O}\right)}^2}}$$

(7)

$$RMSE=\sqrt{{\displaystyle \sum }_{i=1}^n\frac{{\left(P_i-O_i\right)}^2}{n}}$$

(8)

where $Y$ is the dependent variable or predictand, ${\beta }_0$ is the intercept, ${\beta }_i$ is the coefficient of the ith independent variable or predictor $X_i$, $\epsilon$ is the error or residuals of the data, $P_i$ is the predicted value, $O_i$ is the observed value, $\overline{P}$ is the mean of the predicted, and $\overline{O}$ is the mean of the observed.

2.5. Bias Correction

Although the model was developed using standardized NCEP reanalysis outputs, standardized GCM ensemble outputs were used to produce the projections for the future under RCPs 4.5 and 8.5. Before projecting into the future, the historical GCM ensemble was used to validate the model. The validation essentially compared historic model predictions when applied using predictors from NCEP versus those from the GCM ensemble. A large bias in the GCM-derived model output was detected, highlighting the need for correction. The bias arises since the SD model was calibrated with the NCEP data rather than the GCM data [58]. However, calibrating with NCEP data has the advantage that the coarse NCEP predictors track the realized historic climate, such that the SD model calibration step only focuses on the downscaling from coarse to fine scale (and not on correcting coarse scale biases). Meanwhile, model validation as well as future projections were also conducted based on historical and future GCM data, which then requires a separate step to correct biases in GCM predictors. Since the GCM time series does not need to produce the same realization of the climate dynamics as NCEP or observations, the correction should focus on the probability distribution of the predictors rather than on a deterministic reproduction of observations. Thus, before projecting into the future, the historical GCM ensemble output was corrected against the observed data using quantile mapping and that correction was adopted for future projections. Quantile mapping uses the quantile-quantile relation to converge the simulated variables’ probability distribution function into the observed one. Some evaluations found that quantile mapping is one of the best-performing methods for correction [80,81]. Quantile mapping uses a statistical transformation to transform the distribution functions of the modeled variables into the observed ones using a mathematical function expressed in Equation (9) [82]:

$$P_0={F_0}^{-1}\left(F_m\left(P_m\right)\right)$$

(9)

where $P_0$ and $P_m$ are the observed and modeled variables, respectively, $F$ is the cumulative distribution function (CDF), $F^{-1}$ is the corresponding quantile function (inverse CDF), $F_0$ is the CDF of the observed data, and $F_m$ is the CDF of the modeled data.

2.6. Scenario Generation

The validated and corrected regression models were then applied to generate future scenarios for the 41 stations utilizing the CMIP5 GCM ensemble data. The study assumes that the relationship between the predictors and temperature or precipitation remains valid under future climate conditions. The generation of future scenarios intends to predict precipitation and temperature between 2006 and 2050. Figure 3 outlines the overall modeling framework used in the downscaling process.

3. Results

3.1. Selection of Predictors

For every station, a different set of selected predictors was used for temperature, precipitation occurrence, and precipitation amount.

Table 5. The frequency of occurrence (as % of times used across all models) of the predictors in the final generated MLR models for temperature.

	Temperature	Precipitation Occurrence	Precipitation Amount
Predictor	Frequency
Temp2m	100%	3%	18%
Pressure	27%	3%	3%
UWND.500	0%	3%	0%
UWND.1000	16%	0%	0%
VWND.500	0%	3%	26%
VWND.1000	0%	3%	68%
RHUM.500	0%	3%	3%
RHUM.1000	11%	6%	56%
SHUM.500	0%	3%	12%
RHUM.1000	11%	6%	56%
SHUM.500	0%	3%	12%
SHUM.1000	0%	0%	0%
HGT.500	29%	97%	79%
HGT.850	0%	94%	41%

Temp2m = Temperature at 2 m, UWND.500 = U wind component (East/West) at 500 mb pressure level, UWND.1000 = U wind component (East/West) at 1000 mb pressure level, VWND.500 = V wind component (North/South) at 500 mb pressure level, VWND.1000 = V wind component (North/South) at 1000 mb pressure level, RHUM.500 = Relative Humidity at 500 mb pressure level, RHUM.1000 = Relative Humidity at 1000 mb pressure level, SHUM.500 = Specific Humidity at 500 mb pressure level, SHUM.1000 = Specific Humidity at 1000 mb pressure level, HGT.500 = Geopotential Height at 500 mb pressure level, HGT.850 = Geopotential Height at 850 mb pressure level.

summarizes the frequency of the chosen predictors for temperature and precipitation. Note that for temperature, we had nine stations. Each had five models (four seasonal models and one annual model). For precipitation, we had 34 stations in total and we developed for each station an individual model for each month.

Note that the 2 m air temperature was selected in all temperature downscaling models, and in some models, it was the only predictor with a correlation reaching up to 0.9, reflecting the strong physical relationship between the observed temperature and the coarse-gridded 2 m air temperature. Interestingly, for precipitation, the coarse grid precipitation was not included in any of the models because it had a very low correlation with the station scale precipitation in the model despite the meaningful physical relationship. This finding is consistent with the results reported by [83]. Meanwhile, the U wind component (East/West) at 500 mb pressure level, the V wind component (North/South) at 500 mb pressure level, the V wind component (North/South) at 1000 mb pressure level, relative humidity at 500 mb pressure level, specific humidity at 500 mb pressure level, specific humidity at 1000 mb pressure level and geopotential height at 850 mb pressure level did not have a good correlation with temperature. Note that for some stations, these predictors had a good correlation with the climatic variable (precipitation or temperature) but were found to cause multi-collinearity and thus had to be removed.

With regards to rain occurrence, the prevailing predictors were geopotential height at 500 mb and 850 mb pressure level. They were used in 97% and 94% of the stations, respectively. For the rain amount, the prevailing predictors were geopotential height at 500 mb in 79% of the stations, V wind component (North/South) at 1000 mb in 68% of the stations, relative humidity at 1000 mb in 56% of the stations, and geopotential height at 850 mb in 41% of the stations. The 2 m air temperature only appeared six times. Similar predictors were observed in other studies, such as [84,85].

3.2. Downscaling

After the selection of predictors and model calibration, downscaling using MLR was conducted, with the data divided into calibration and validation periods. The period of calibration and validation varied between stations due to data availability (Table 2). Predictors were standardized before calibration and validation.

3.2.1. Precipitation

The cumulative change was calculated for every station alone, and then the resulting cumulative change was averaged over all the stations. As for the bias related to elevation, this is indeed a possibility that cannot be tested. The downscaling of precipitation was performed in two steps. First, the logistic regression yielded the months when rain occurred. For most stations, November, December, January, February, March, and April were the wet months and May, June, July, August, September, and October were the dry months, which is reflective of the current weather in the region. Figure 4 shows the differences in McFadden R² for the rain occurrence, with the lowest value at the Jubeiha station (R² = 47%) and the highest value at the Jarash station (R² = 82%). After the division between dry and wet months, the MLR was used to model precipitation amounts in the wet months, with Figure 5 depicting the differences in R² across the stations. Note that the root mean squared error (RMSE) for precipitation ranged between 10.91 mm (at Nawasif) and 55.72 mm (at Kitta), with

Table 6. R² and RMSE for precipitation stations.

Station	R2 (%)	RMSE (mm)	Station	R2 (%)	RMSE (mm)
Har Kenaan	72	40.71	Um Qeis	66	23.06
Ammanhc	63	36.37	Kharja	65	29.04
Balama	55	29.60	Husn	62	15.11
Deir Alla	53	27.47	Nueiyime	55	19.13
Jarash	59	38.11	Ramtha	50	20.98
Jubeiha	54	53.82	Khanasira	48	14.09
Kitta	67	55.72	Mafraq	41	15.78
Midwar	67	33.26	Turra	62	28.11
Nawasif	55	10.91	Hosha	41	15.89
Prince Feisal Nursery	54	28.64	Jaber	45	19.69
Qafqafa	68	41.17	Baqura	55	29.24
Rumeimin	63	33.47	Irbid	62	46.4
Ruseifa	51	13.52	Sukhna	46	15.66
Sihan	53	41.12	Um El Jamal	52	10.65
Subeihi	51	43.30	Um Jauza	52	47.89
Jerusalem	62	54.50	Wadi Dhuleil	44	15.24
Kufr Saum	60	44.30	K H Nursery	49	40.71

presenting the precipitation RMSE for all stations.

3.2.2. Temperature

Daily temperatures were simulated using different MLR models. Five models were developed for each station: one for each season and one representing the whole period using daily data. The models were able to simulate temperature relatively well, with R² ranging between 62 and 95% (

Table 7. Temperature stations and results.

Station	Season	R2 (%)	RMSE (°C)
Beirut Airport	Winter	71	1.09
	Spring	83	1.29
	Summer	62	1.02
	Fall	87	1.09
	One Model	91	1.61
Damascus	Winter	60	1.59
	Spring	85	1.96
	Summer	65	1.38
	Fall	86	1.95
	One Model	91	2.42
H4 Airbase	Winter	74	1.58
	Spring	88	2.08
	Summer	62	1.66
	Fall	90	1.80
	One Model	95	1.93
Ma’an	Winter	77	1.57
	Spring	91	1.72
	Summer	68	1.50
	Fall	89	1.69
	One Model	94	1.80
Prince Hassan	Winter	68	1.64
	Spring	90	1.82
	Summer	71	1.48
	Fall	89	1.72
	One Model	95	1.86
Ghor Safi	Winter	45	1.69
	Spring	81	1.75
	Summer	53	1.30
	Fall	85	1.62
	One Model	91	2.12
King Hussein	Winter	74	1.05
	Spring	89	1.29
	Summer	67	1.02
	Fall	90	1.09
	One Model	95	1.66
Jerusalem	Winter	83	1.60
	Spring	91	1.91
	Summer	76	1.28
	Fall	84	1.73
	One Model	92	1.91
Har Kenaan	Winter	73	1.61
	Spring	87	2.34
	Summer	74	1.28
	Fall	87	2.15
	One Model	91	2.49

). Note that the RMSE for the temperature ranged between 1.02 °C (at Beirut airport and King Hussein) and 2.49 °C (at Har Kenaan), with Table 7 presenting the temperature RSME for all stations.

3.3. Bias Correction

The impact of bias correction on future projections of precipitation and temperature was large, as shown in Figure 6a,b. Evidently, the quantile mapping method was able to reduce the bias between the raw GCM on the one hand and the observed data or NCEP on the other.

3.4. Scenario Generation

The future GCM ensemble was corrected through the bias correction equation generated in the previous step. This corrected GCM ensemble was then used in the multiple linear regression model for projections into the future for RCP4.5 and RCP8.5. Multiple linear regression equations varied for each station and for each season in the case of temperature, with the corresponding equations detailed in Tables S3–S5 of the Supplemental Material. The results, detailed in Table S6 of the Supplemental Material, show that the decrease in precipitation under the RCP4.5 ranged between 0.26 mm/year and 7.17 mm/year, and under the RCP8.5 scenario, the decrease ranged between 1.10 mm/year and 9.53 mm/year, depending on the station, while the increase in temperature ranged between 0.02 °C/year and 0.09 °C/year under the RCP4.5, and between 0.034 °C/year and 0.09 °C/year under RCP8.5. This trend was consistent throughout all the stations of temperature and precipitation for the whole simulation period. The annual changes for all stations for precipitation and temperature, along with the p-Values, are presented in

Table 8. Change in precipitation for RCP 4.5 and RCP 8.5 scenarios.

Station	Change (RCP 4.5) per Year (mm)	p-Value for RCP 4.5 Slope	Change (RCP 8.5) per Year (mm)	p-Value for RCP 8.5 Slope
Har Kenaan	−3.44	<0.05	−7.79	<0.05
Ammanhc	−3.30	<0.05	−6.68	<0.05
Balama	−1.04	0.3	−3.02	<0.05
Deir Alla	−1.56	<0.05	−5.58	<0.05
Jarash	−5.53	<0.05	−5.08	<0.05
Jubeiha	−4.18	<0.05	−6.93	<0.05
Kitta	−7.17	<0.05	−9.53	<0.05
Midwar	−2.35	<0.05	−3.94	<0.05
Nawasif	−1.07	<0.05	−1.71	<0.05
Prince Feisal Nursery	−4.66	<0.05	−4.50	<0.05
Qafqafa	−1.38	0.468	−5.76	<0.05
Rumeimin	−3.66	<0.05	−5.38	<0.05
Ruseifa	−1.49	<0.05	−1.09	<0.05
Sihan	−4.58	<0.05	−4.67	<0.05
Subeihi	−5.18	<0.05	−3.21	<0.05
Jerusalem	−3.44	<0.05	−4.96	<0.05
Kufr Saum	−4.33	<0.05	−5.48	<0.05
Um Qeis	−3.68	<0.05	−5.62	<0.05
Kharja	−3.80	<0.05	−4.40	<0.05
Husn	−4.16	<0.05	−5.27	<0.05
Nueiyime	−3.10	<0.05	−3.55	<0.05
Ramtha	−2.40	<0.05	−2.21	<0.05
Khanasira	−2.03	<0.05	−2.44	<0.05
Mafraq	−1.36	<0.05	−1.55	<0.05
Turra	−2.37	<0.05	−3.18	<0.05
Hosha	−1.19	<0.05	−1.98	<0.05
Jaber	−0.26	0.764	−1.85	0.03
Baqura	−3.43	<0.05	−4.71	<0.05
Irbid	−3.29	<0.05	−4.12	<0.05
Sukhna	−0.92	0.08	−1.83	<0.05
Um El Jamal	−1.34	<0.05	−1.66	<0.05
Um Jauza	−5.60	<0.05	−5.74	<0.05
Wadi Dhuleil	−0.87	0.03	−1.24	<0.05
K H Nursery	−5.26	<0.05	−3.45	<0.05

and

Table 9. Change in temperature for RCP 4.5 and RCP 8.5 scenarios.

Station	Season	Change (RCP 4.5) per Year (°C)	p-Value (RCP 4.5) Slope	Change (RCP 8.5) per Year (°C)	p-Value (RCP 8.5) Slope
Beirut Airport	Winter	0.03	<0.05	0.06	<0.05
	Spring	0.02	<0.05	0.04	<0.05
	Summer	0.04	<0.05	0.04	<0.05
	Fall	0.03	<0.05	0.04	<0.05
	One Model	0.03	<0.05	0.03	<0.05
Damascus	Winter	0.04	<0.05	0.06	<0.05
	Spring	0.03	<0.05	0.04	<0.05
	Summer	0.05	<0.05	0.06	<0.05
	Fall	0.04	<0.05	0.06	<0.05
	One Model	0.04	<0.05	0.05	<0.05
H4 Airbase	Winter	0.05	<0.05	0.08	<0.05
	Spring	0.04	<0.05	0.05	<0.05
	Summer	0.06	<0.05	0.07	<0.05
	Fall	0.04	<0.05	0.06	<0.05
	One Model	0.04	<0.05	0.05	<0.05
MA AN	Winter	0.04	<0.05	0.07	<0.05
	Spring	0.03	<0.05	0.05	<0.05
	Summer	0.08	<0.05	0.09	<0.05
	Fall	0.04	<0.05	0.06	<0.05
	One Model	0.03	<0.05	0.05	<0.05
Prince Hassan	Winter	0.05	<0.05	0.07	<0.05
	Spring	0.03	<0.05	0.06	<0.05
	Summer	0.09	<0.05	0.09	<0.05
	Fall	0.04	<0.05	0.06	<0.05
	One Model	0.04	<0.05	0.05	<0.05
Ghor Safi	Winter	0.02	<0.05	0.05	<0.05
	Spring	0.03	<0.05	0.04	<0.05
	Summer	0.05	<0.05	0.05	<0.05
	Fall	0.04	<0.05	0.05	<0.05
	One Model	0.03	<0.05	0.05	<0.05
King Hussein	Winter	0.04	<0.05	0.07	<0.05
	Spring	0.03	<0.05	0.05	<0.05
	Summer	0.07	<0.05	0.07	<0.05
	Fall	0.04	<0.05	0.06	<0.05
	One Model	0.03	<0.05	0.05	<0.05
Jerusalem	Winter	0.07	<0.05	0.09	<0.05
	Spring	0.04	<0.05	0.06	<0.05
	Summer	0.08	<0.05	0.08	<0.05
	Fall	0.04	<0.05	0.06	<0.05
	One Model	0.03	<0.05	0.05	<0.05
Har Kenaan	Winter	0.06	<0.05	0.08	<0.05
	Spring	0.04	<0.05	0.07	<0.05
	Summer	0.04	<0.05	0.06	<0.05
	Fall	0.04	<0.05	0.06	<0.05
	One Model	0.04	<0.05	0.05	<0.05

, respectively. The temperature and precipitation Parametric Test and the p-Values are detailed in Tables S7 and S8 of the Supplemental Material.

4. Discussion

In general, the results of the logistic model were realistic and, as expected, did not vary much across the region (Figure 4). R² ranged between 41% in the Mafraq Airport station to 72% in Har Kenaan station. Following the separation of the months into dry and wet ones, the MLR was utilized to estimate precipitation during the wet months. Figure 5 shows the variations in R² between the stations. Note that no pattern could be discerned between model skill (in terms of R²) and station elevation or its spatial location. Yet, we observed that stations with more observed data and fewer gaps in their data gave better results, probably due to better calibration. Thus, we are assuming that we will always have rain in the wet months (as predicted by the logistic model) for every station. A comparison of precipitation downscaling results between this study and previous studies [84,86] reveals similar performance.

Regarding the temperature, the highest R² was associated with the models that disregarded seasons, which is an important practical finding since it indicates that the “simplest” model that does not segregate data by season performs the best. The lowest R² was observed for the summer and winter season models, while the models for the fall and spring had higher R² values. A comparison with previous research, namely [84,86], shows that the results are aligned even when the study area and the downscaling method differed; all studies depicted good results for downscaling temperature.

The necessity for bias correction reinforced the need to implement the bias correction step before using the GCM future predictions in an MLR model that is calibrated based on reanalysis data (NCEP in this case). The MLR analysis for future projections showed a change in the climate variables with a decreasing pattern of precipitation and an increasing pattern for temperature under the RCP4.5 and 8.5. A detailed assessment of the future predictions indicates that the predictions pointed to an increase in extreme precipitation and temperature as the skewness and variance of future projections increased compared to the observed data. This trend was consistent throughout all the stations of temperature and precipitation for the whole simulation period. Several recent studies have reported reductions in rainfall and increases in temperature for similar basins. For example, [87] observed decreasing trends of precipitation over Onkaparinga, Australia for the period 2041–2060 (Table 8). Samadi et al. [88] stated that daily temperature will increase and precipitation will decrease in a semi-arid catchment in western Iran, which is consistent with the results presented in Table 9. Hertig and Jacobeit [89] used two statistical downscaling methods, canonical correlation coupled with MLR analysis, to assess the expected Mediterranean precipitation changes for the period 1990–2100 under increased greenhouse gas conditions. They reported mainly negative precipitation changes for the rainy season ranging from October to May, similar to the results obtained here. In addition, a dynamical downscaling [90] conducted on the same catchment area projected that precipitation would decrease and temperatures increase in selected extreme years (2020, 2029, 2040, and 2050) for RCP4.5 and (2017, 2023, 2035, and 2050) for RCP8.5, which parallels the outputs of the conducted SD.

5. Conclusions

In this study, high-resolution statistical downscaling (SD) and GCM simulations with data re-gridding and correction were used to define climate variables under RCP 4.5 and RCP 8.5 scenarios until the year 2050. Precipitation and temperature were downscaled at 41 stations in the Jordan River Basin, with the aim to decrease uncertainties in predicting climate variables toward helping in the development of strategies that can meet the challenge of water scarcity in a region where climate change is expected to exacerbate existing water shortages.

The statistical downscaling approach consisted of adopting MLR models that were developed for each station using coarse historical reanalysis data from NCEP. These MLR models were then used to project future trends using coarse data from an ensemble of GCMs after these data were bias-corrected based on historical records spanning 1981–2005. Aggregated unified annual models performed best for temperature, and aggregated monthly models performed best for precipitation, with correlations reaching 93% and 78%, respectively. While the SD proved effective with temperature downscaling, exhibiting adequate skill, precipitation downscaling will still benefit from further improvements.

The future SD results for the period between 2006 and 2050 showed an increase in temperature and a decrease in precipitation under both the RCP 4.5 and RCP 8.5 scenarios. The annual increase in temperature ranged between 0.02 and 0.09 °C/year under RCP 4.5 and between 0.034 and 0.09 °C/year under RCP 8.5, depending on the station, which represents a cumulative surface temperature increase of 1.54 and 2.11 °C under the RCP 4.5 and RCP 8.5. In contrast, an annual decrease in precipitation is expected, ranging between 0.26 and 7.17 mm/year under the RCP 4.5 and between 1.10 and 9.53 mm/year under RCP 8.5, depending on the station, which represents a cumulative decrease of 100 and 135 mm under the RCP 4.5 and RCP 8.5, respectively, or the equivalent reduction of 10 and 15%, respectively. This pattern will inevitably add stress to water resources, increasing management challenges in semi-arid to arid regions of the basin.

The precipitation and temperature downscaling process used in this study was consistent with previous research [91,92], which demonstrated the superiority of using circulation variables (e.g., geopotential, vorticity, or the wind component) and temperature (e.g., geopotential heights at various levels and specific/relative humidity near the mid-troposphere and specific/relative humidity) in order to establish a satisfactory relationship when downscaling temperature and/or precipitation as opposed to any single predictor. Given that it is easier to apply and has “less preprocessing requirements and computational costs” Tavakol-Davani et al. [93], in the future, the SD is likely to be examined more in the JRB than other downscaling techniques.

At a broader level, SD proved to be efficient in predicting climate change variables and hence can be relied upon by policymakers for climate change analysis and water resources management at a river basin scale. Nevertheless, it is imperative to recognize the need for future work to enhance accuracy, including the comparison with several GCMs and various statistical methods, as well as uncertainties in the downscaling process, the bias correction, and the quality of the observed data.