Quantifying Climate Change Variability for the Better Management of Water Resources: The Case of Kobo Valley, Danakil Basin, Ethiopia

Abstract

Alterations in the hydrological cycle due to climate change are one of the key threats to the future accessibility of natural resources. This study used 12 GCM climate models from CMIP6 to evaluate future climate change scenarios by applying model performance measures and trend analysis in Kobo Valley, Ethiopia. The models were ranked based on their ability to analyze the historical datasets. The result of this study showed that the outputs of the FIO-ESM-2-0 CIMP6 model had a good overall ranking for both precipitation and temperature. After bias correction of the model-based projections with the observed data, the average annual precipitation in the average scenario (SSP2-4.5) decreased by 4.4% and 13% in 2054 and 2084, respectively. Similarly, in the worst-case scenario (SSP5-8.5), by the end of 2054 and 2084, decreases of 4% and 12.8%, respectively, were predicted. The average annual maximum temperature under the SSP2-4.5 scenario increased by 1.5 °C in 2054 and by 2.1 °C in 2084. The average annual maximum temperature under the worst-case (SSP5-8.5) scenario increased by 1.7 °C in 2054 and by 3.2 °C in 2084. In the middle scenario (SSP4.5), the average annual minimum temperature increased by 2.2 °C in 2054 and by 3 °C in 2084. The average annual minimum temperature under the worst-case (SSP5-8.5) scenario increased by 2.6 °C in 2054 and by 4.3 °C in 2084. The seasonal variability in precipitation in the studied valley will decrease in the winter and increase in the summer. A decrease in precipitation combined with an increase in temperature will strengthen the risk of drought events in the future.

1. Introduction

Global climate change encompasses changes in the atmospheric composition as well as changing interactions between the atmosphere and other various geological, chemical, and biological factors within Earth’s system [1,2,3,4]. According to [5], global climate change is driven by natural or human forces that can lead to changes in the likelihood of the strength of extreme weather and climate events. Changes in climate variables have additional influences on Earth’s system, such as changes in land use/land cover, water resource availability, and environmental hazards [6,7,8].

The results of human-induced pressure can directly or indirectly alter the composition of the atmosphere, and consequently, the natural climate variability is observed over comparable referenced time scales [9]. Natural processes such as volcanic eruptions [10] and human activities such as land use and deforestation [11] can have a significant role in changing greenhouse gases. Different natural and anthropogenic factors are responsible for changing precipitation patterns [12]. Greenhouse gas concentrations are used as input statistics for General Circulation Models (GCMs) [13].

According to the Intergovernmental Panel for Climate Change (IPCC) 2013 reports, changes in rainfall and temperature in many parts of Africa are causing fluctuations in freshwater and affecting the quality and quantity of water accessibility. The GCM and analyses of their products provide historical and projected values of climate variables.

Several meteorological stations lack any recorded data or their records are over short periods. Therefore, satellite data have reduced the challenges associated with the shortage of station observations and enhanced the spatial and temporal continuity in meteorological data in data-scarce regions, thereby creating great progress in the study and understanding of climate change. After the 1990s, satellite data entered data assimilation systems, which further improved the accuracy of meteorological data [14]. Currently, many scholars have used satellite or reanalysis products to represent the spatiotemporal variability in climate change across the world. In this study, we used satellite-based data, reanalysis models, and the combination of satellite data with reanalysis products to investigate the spatiotemporal variability in climate change in Kobo Valley.

The global climate change projections were carried out using GCMs, which provide future climate change variables at large spatial scales. The Intergovernmental Panel for Climate Change (IPCC) released five series of reports on climate change. However, the CMIP6 models are more advanced in terms of modeling groups, the number of projection scenarios, and the number of various experiments involved [15,16,17,18,19]. The CMIP6 models have a broader range of complexity as compared with the previous assessment phases due to several improved physical processes and their spatial resolution. Climate change scenarios provide critical results in climate models through their description and the involvement of human-induced as well as natural processes for emissions trajectories [20].

The Intergovernmental Panel for Climate Change (IPCC) released their sixth assessment report (AR6) based on the latest CMIP6 models [21]. A significant improvement in the CMIP6 models compared with the previous CMIP5 models is the inclusion of socioeconomic development factors with GHG emission scenarios (representative concentration pathways (RCPs)) [22]. The CMIP6-GCMs present aspects of improvement over previous generations, such as higher spatial resolution and better parameterization schemes of the physical and biogeochemical processes of the climate system [3].

General Circulation Models (GCMs) are one of the primary tools for understanding future climate projections. Shared Socioeconomic Pathways (SSPs) characterize a more realistic socioeconomic development by considering different social, economic, technological, and political scenarios [23]. Several studies have revealed that CMIP6 models show better simulation performance in terms of hydrological responses and impact assessment [17,24,25,26,27,28,29]. However, the spatial scale of these GCMs is generally hundreds of kilometers in order to represent the Earth’s system, including land, oceans, and the atmosphere. Therefore, capturing local-scale details should be improved in order to generate promising results for users working on regional-scale studies [30].

The seasonal and hydrological variability in climate change has been experienced in different parts of Ethiopia [31,32,33,34]. However, different studies have identified varying scales and magnitudes of their impacts. For example, the precipitation projections in ref. [31] show a slightly (statistically insignificant) increasing trend for the near (long)-term periods in the Upper Blue Nile Basin. Natural resources in Ethiopia are highly vulnerable to climate variability because of their topography and anthropogenic factors such as land degradation, increasing population, and natural resource management practices [35,36].

Most studies conducted in Ethiopia have used a single GCM selection criterion and thus may not provide accurate model selection and impact assessment. Moreover, the climate change variability in the Danakil Basin has not been studied with an appropriate climate model projection. Therefore, in this study, we used several model selection criteria to select the best representative models for seasonal and annual hydrological variability. Minimizing model discrepancies is important and was conducted in this study by using multi-model selection criteria of GCMs with the most recent naval approach (state-of-the-art) to evaluate the individual CMIP6 models. The accurate evaluation of an individual model’s performance before selecting the representative GCMs for a particular impact assessment is also important. In this regard, comprehensive studies are required to select improved CMIP6 models for simulating climate change variability over various spatiotemporal scales.

In recent times, extreme events like momentary devastating floods, prolonged droughts, and water scarcity that limit agricultural activities and socioeconomic activities have become frequent phenomena in Kobo Valley. The variability in hydrometeorological processes is worth investigating using the global climate model. This investigation could also evaluate the climate models’ reliabilities, identifying their strengths and restrictions, as well as identifying the best-performing model for a specified location and time. The new set of investigations will also help to evaluate the applications of machine learning potentials for climate change prediction that would suggest a range of strengths and limitations for a specified location.

Therefore, the objective of this study is to quantify the variability and trends in climatic variables in Kobo Valley. In this study, (a) the correlations and biases of four reanalyzed meteorological variables from reanalysis products (ERA5, MERRA-2, JAR-55-mdl-iso, 20CRv3), satellite (CHIRPS), and merged products of the Climate Prediction Center (CPC) are compared based on the daily weather variable data obtained from 12 ground meteorological stations from 1985 to 2014 and (b) the meteorological variables from CMIP6 models are projected to determine future possible trends. Moreover, a comprehensive set of statistical indicators is used to assess how the calculated climate variables from the reanalysis and merged products compare with the corresponding observational data.

2. Material and Methods

2.1. Area Description

Kobo Valley is located in the Danakil Basin, Ethiopia. The valley has an altitude ranging from 1136 m to 3970 m above sea level. It is found on the intermountain plain between 11.84° and 12.34° north and 39.34° and 39.84° east. The area covers 1228 km². Figure 1 shows the location of Kobo Valley.

2.2. Dataset

The observed daily temperature and precipitation data were collected from the Ethiopian Meteorological Agency. However, in the Danakil Basin, there is limited data availability and various gaps in the time series data. The installation and accessibility periods of the gauged stations were divided into three groups as follows: Group 1, precipitation stations with 30 years (1985–2014) of records, Group 2 with 16 years (1998–2014) of records, and Group 3 with 9 years (2006–2014) of records.

In this study, the average monthly, seasonal (summer, winter), and annual data were derived from the daily meteorological datasets. The summer season starts in June and ends in September (JJAS), while winter begins in March and lasts until May (MAM). The satellite data, the Reanalysis model, and the merged results of satellite data with reanalysis model products were analyzed and compared with historical meteorological data.

Table 1. Available meteorological stations in Kobo Valley.

No	Name of the Station	Location (WGS)
		Latitude	Longitude	Precipitation	Temperature
Group 1
1	Alamata	12.42	39.71	**	**
2	Waja	12.30	39.60	**	**
3	Korem	12.51	39.51	**	**
4	Lalibela	12.04	39.04	**	**
5	Sirinka	11.75	39.61	**	**
Group 2
6	Kobo	12.13	39.63	**	**
7	Zobel	12.17	39.75	**	**
8	Kulmesk	11.94	39.20	**
Group 3
9	Dilb	11.96	39.43	**
10	Sanka	11.89	39.48	**	**
11	Hara	12.67	41.42	**	**
12	Robit	12.01	39.62	**

** Selected for both variables.

shows the meteorological stations used for this study and the grouping based on measurement periods.

Several gauging stations have a lack of consistency in their measured data. Various spatial interpolation techniques are available today with varying degrees of complexity to fill data gaps and maintain consistency. The inverse distance method is the most widely used method for filling missing data [37]. In this method, the weights for each sample are inversely proportional to its distance from the point being estimated [38]. We used the inverse distance method to complete the missing precipitation data of stations that have other stations nearby, using the following equation (Equation (1)):

$$\mathrm{Px}={\displaystyle \frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{N}}\frac{1}{{\mathrm{di}}^2}\mathrm{Pi}}{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{N}}\frac{1}{{\mathrm{di}}^2}}}$$

(1)

where Px = estimate of rainfall for the ungauged station in mm, Pi = rainfall values of rain gauges used for estimation in mm, di = distance from each location to the point being estimated in km, and N = No. of surrounding stations.

The historical precipitation, maximum, and minimum temperature of the satellite and merged products over the study area (Kobo Valley) were retrieved from the Royal Netherlands Meteorological Institute Climate Explorer website (https://climexp.knmi.nl/) accessed on 28 March 2023. The reanalysis data were retrieved from the Collaborative Reanalysis Technical Environment Intercomparing Project (CREATE-IP) under the designation in the Earth System Grid Federation (ESGF) website (https://esgf-node.llnl.gov/projects/create-ip/) accessed on 28 March 2023.

The satellite products of the Climate Hazards Group InfraRed Precipitation with Station (CHIRPS) data, and the reanalysis products of European Community Medium-range Weather Forecasts v5 (ECMWF-ERA5), NASA Modern Reanalysis Evaluation (NASA-MERRA v2), Japanese Meteorological Agency (JRA-55-mdl-iso), and Twentieth Century Reanalysis V3 (20CRv3) were used. The satellite, reanalysis model, and CPC merged products were used and validated by the observational data for the baseline period (1985–2014). The following table (

Table 2. Data sources of satellite, reanalysis, and merged products.

No	Name	Data Product	Resolution	Source
Precipitation
1	CHIRPS	Satellite	0.25° × 0.25°	NOAA
2	ERA5	Reanalysis	0.25° × 0.25°	ECMWF
3	MERRA-2	Reanalysis	0.5° × 0.625°	NASA
4	JRA-55-mdl-iso	Reanalysis	0.5625° × 0.5625°	JMA
5	CPC	Merged	0.25° × 0.25°	NOAA
Temperature (Tmax, Tmin)
1	ERA5	Reanalysis	0.25° × 0.25°	ECMWF
2	MERRA-2	Reanalysis	0.5° × 0.625°	NASA
3	20CRv3	Reanalysis	1° × 1°	NOAA
	Evapotranspiration
1	ERA5	Reanalysis	0.25° × 0.25°	ECMWF

) shows the data products, resolutions, and their sources for the selected variables.

Estimation of Potential Evapotranspiration (PET)

The other important hydrological variable is evapotranspiration, which is governed by several factors, where the most important are temperature, humidity, and the availability of water [39]. There are different methods to quantify the potential evapotranspiration of a given watershed. The FAO (Food and Agriculture Organization) recognized that the Penman–Monteith method (PM) is a standard and widely acceptable method for potential evapotranspiration (PET) estimation. However, because of detailed input meteorological datasets, the application of the Penman–Monteith method (PM) is often limited in many regions. The model requires several climate variables, including the air temperature (T_max and T_min), relative humidity (RH), solar radiation (R_s), and wind speed at 2 m height (U₂). Most of these climate variables are not available at the meteorological stations in many regions in Ethiopia. In such cases, selecting an alternative method with similar efficiency to that of the standard technique has to be identified [40].

The Hargreaves method is an alternative method for estimating PET if insufficient meteorological data are available and is recommended by the FAO [41]. The Hargreaves method is a simple evapotranspiration model that requires a few input parameters. For this study, the PET was simulated by using the Hargreaves method in Standardized Precipitation–Evapotranspiration Index (SPEI) packages written in R [42]. The Hargreaves equation requires input data on radiation and temperature. The advantage of this method is that no radiation has to be measured since the model works with calculated extra-terrestrial radiation. Equation (2) represents the potential evapotranspiration estimation using the Hargreaves method.

$$\mathrm{PET}=0.0135\ast \mathrm{Rs}\ast \mathrm{Conv}\ast \left(\mathrm{T}+17.8\right)$$

(2)

where PET is potential evapotranspiration (mm day⁻¹), T is the mean temperature of the day (°C), Rs is solar radiation (MJ m⁻² day⁻¹), and Conv is conversion to the ET factor (m² mm MJ⁻¹). This method uses nearest neighbor interpolation to simulate the spatiotemporal distribution of potential evapotranspiration.

2.3. Future Scenarios

The temperature and precipitation data were extracted from the World Climate Research Programme’s (WCRP’s) CMIP6 under a medium forcing scenario (SSP2-4.5) and a strong forcing scenario (SSP5-8.5). For future climate scenarios, the monthly weather data were investigated for 30-year time intervals for the baseline period (1985–2014), the near term (2025–2054), and the long term (2055–2084).

The selection of future projections was based on historical data and was considered in three time steps, which are average monthly, seasonal (JJAS, MAM), and annual. This allowed us to grasp all statistical analyses of the time series climate data for the selected individual model. The primary operations were the analysis of CMIP6-driven weather datasets and the evaluation/calibration with station data. Overall, 12 global climate models were obtained from the WCRP platform (https://esgf-node.llnl.gov/search/cmip6/). The CMIP6 model products were accessed on 28 March 2023. We used several criteria to evaluate the individual model results in terms of the monthly basis, the grid level (gn), the nominal resolution (100 km), the type of source (AOGCM), and the variant level (r1i1p1f1). The selected individual CMIP6 products (acronyms, resolutions, countries, and references) are summarized in

Table 3. Selected CMIP6 climate model products.

No	CMIP6 Model	Resolution	Country
1	ACCESS-ESM1-5	1.875°× 1.25°	Australia
2	BCC-ESM1 (**)	1.1° × 1.1°	China
3	CAMS-CSM1-0	1.1° × 1.1°	China
4	CanESM5	2.8° × 2.8°	Canada
5	CAS-ESM2-0 (**)	1° × 1°	China
6	CESM2-WACCM	1.3° × 0.9°	USA
7	CMCC-ESM2 (**)	0.9° × 1.25°	Italy
8	FIO-ESM-2-0 (**)	0.9° × 1.25°	China
9	MIROC6	1.4° × 1.4°	Japan
10	MRI-ESM2-0 (**)	1.1° × 1.1°	Japan
11	NorESM2-LM	1.9° × 2.5°	Norway
12	TaiESM1	0.9° × 1.25 °	Taiwan

** Selected CMIP6 models for both variables.

.

The process started with extracting the individual outputs of GCMs by location and period. The extraction process was performed with CDO (Climate Data Operator). Since all climate models have coarse resolution, spatial downscaling was performed considering the outputs of the GCM as grid points. The Inverse Distance Weighting (IDW) method was used to extract four surrounding grid points for each GCM output to the observed station.

The individual GCM outputs were evaluated based on their statistical relationships to the observed datasets on time levels of average monthly, seasonal, and annual bases, followed by identifying the best possible probability distributions and their trends. Bias correction methods were applied after the extraction of GCM outputs. Various bias correction methods are available [43]. In this study, the quantile mapping method was applied to minimize errors. The outputs of individual GCMs are biased because of the lack of parameters, inadequate length of recorded data, and low spatial resolution.

2.4. Model Performance

2.4.1. Statistical Analysis

Model performance and trend studies were applied using multiple statistical analyses. We used the correlation coefficient, the root mean square error, the mean absolute error, probability distributions, and trend analysis. The mathematical expressions of these indexes are given below (Equations (3)–(5)) [44]. The correlation coefficient is a value that indicates the strength of the relation between variables, and it ranges from −1 to 1.

$$\mathrm{R} ={\displaystyle \frac{n({\displaystyle \sum }\left(xy\right)-{\displaystyle \sum }\left(x\right){\displaystyle \sum }\left(y\right))}{\sqrt{[n{\displaystyle \sum }x{i}^2-{({\displaystyle \sum }x)}^2][n{\displaystyle \sum }{y}^2-{({\displaystyle \sum }y)}^{2 }]}}}$$

(3)

where R is the correlation coefficient, n is number of given datasets, x is the observed variable, and y is the model variable.

The root mean square error (RMSE) is a frequently used measure of the difference between values predicted by a model and the values observed from the environment that is being modeled.

$$\mathrm{RMSE}=\sqrt{{\displaystyle \frac{1}{\mathrm{n}}}{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{(\mathrm{Xobs},\mathrm{i}-\mathrm{Xmodel},\mathrm{i})}^2}$$

(4)

where n is the number of observations, (Xobs, i) is the observed value, and (Xmodel, i) is the model value.

The mean absolute error (MAE) is the absolute value of the difference between the model value and the observed value. The mean absolute error measures accuracy of continuous variables.

$$\mathrm{AE}={\displaystyle \frac{1}{\mathrm{n}}}{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}|\mathrm{Oi}-\mathrm{Pi}|$$

(5)

where n is the number of observations and |Oi − Pi| is the absolute error. Both the mean absolute error and the root mean square error express average model prediction error in units of the variable of interest. The mean absolute error and the root mean square error can range from 0 to ∞ and are indifferent to the direction of errors. It is noteworthy that the lower the value of RMSE and MAE, the more reliable the prediction of the data.

2.4.2. Trend Test Analysis

The trends in the time series datasets were detected using parametric or non-parametric tests. The Mann–Kendall (MK) test is a widely used method for non-parametric tests because of its independence on the distribution of the data and because it is not affected by missing datasets and their outliers [45]. Non-parametric methods are alternatives that require few assumptions to be made about the data. These methods are most often used to analyze data that do not meet the distributional requirements of parametric methods. To apply the non-parametric MK test, the data time series does not need to be normally distributed, but it must be serially independent and randomly ordered.

2.4.3. Probability Distributions

The climate data were analyzed to identify the best-fit probability distribution for each time level in the study area. Several studies have investigated precipitation analysis and identified the best-fit probability distribution functions such as Normal, Lognormal, Gumbel, Weibull, and Pearson-type distributions [46]. Probability distributions are selected based on Akaike Information Criteria (AIC) and its distribution type (DT). The smaller the AIC value, the better is the model fit [47]. For the observed data and the individual climate models, the best likely probability distributions were selected with the help of the “gamlss” package in R [48]. In this approach, probability distributions of the observed and reanalysis results were analyzed using “fittest” in R packages.

3. Results

3.1. Data Quality Checking

The observed precipitation and temperature data quality were checked by both distribution types and their trends.

Table 4. Distribution type of the observed, satellite, reanalysis, and merged products.

Datasets	Average Monthly		JJAS		MAM		Annual
Precipitation
	AIC	DT	AIC	DT	AIC	DT	AIC	DT
Observed	226.01	NO	372.61	NO	350.56	RG	375.12	NO
CHIRPS	220.49	SEP2	374.79	PE	326.68	RG	369.58	SEP2
ERA5	255.42	SEP2	361.03	SEP3	370.07	NO	403.44	SEP2
MERRA-2	286.94	LO	417.73	RG	409.13	RG	436.03	LO
JAR-55-mdl-iso	277.42	SN2	384.14	RG	376.94	RG	426.53	SN2
Merged (CPC)	229.62	NO	362.63	NO	339.98	RG	378.72	NO
Evapotranspiration
Observed	236.01	NO	368.61	NO	358.56	RG	379.12	NO
ERA5	254.18	NO	312.45	NO	369.56	SEP3	403.28	NO

AIC—Akaike Information Criteria, DT—distribution type, NO—Normal, LO—Logistic distribution, SN2—Skew Normal Type 2, SEP2 and SEP3—Skew Power exponential 1-4, RG—Reverse Gumbel, PE—Power Exponential distributions. The shaded color indicates that the data sources have a similar distribution type to the observed data.

shows the distribution type of the observed, satellite data, reanalysis, and merged products at different time levels.

The trend test results can evaluate the reanalysis products by considering the similarity in their trends. The merged (CPC) products of the precipitation follow the same trend as the observed datasets at all time levels.

Table 5. Trend test results of the observed, satellite, reanalysis, and merged products.

Dataset	Average Monthly		JJAS		MAM		Annual
Precipitation
	TTR	Dir	TTR	Dir	TTR	Dir	TTR	Dir
Observed	NS	-	NS	+	NS	-	NS	-
CHIRPS	NS	+	NS	+	NS	-	NS	+
ERA5	S	D	NS	+	S	D	S	D
MERRA-21	NS	-	NS	+	S	D	NS	-
JAR-55-md l-iso	S	D	S	D	S	D	S	D
Merged (CPC)	NS	-	NS	+	NS	-	NS	-
Maximum Temperature
Observed	S	I	S	I	S	I	S	I
ERA5	NS	+	S	I	S	I	S	I
MERRA-2	S	I	S	I	S	I	S	I
20CRv3	S	I	S	I	S	I	S	I
Minimum Temperature
Observed	S	I	S	I	S	I	S	I
ERA5	S	I	S	I	S	I	S	I
MERRA-2	S	I	S	I	S	I	S	I
20CRv3	S	I	S	I	NS	+	S	I
Evapotranspiration
Observed	S	I	S	S	S	S	S	S
ERA5	S	I	S	S	S	S	S	S

TTR—trend test result, Dir—direction, S—significant, NS—non-significant, I—increasing, D—decreasing, + and – represent the sign.

shows the trend analysis of the observed, satellite, reanalysis, and merged products at different time levels.

The merged (CPC) products of the precipitation follow the same trend as the observed datasets at all time levels. Figure 2 shows the results of the comparison between observed, satellite, and reanalysis products for average monthly precipitation data. The unit of measurement for precipitation is mm.

The historical temperature datasets were compared to the reanalysis products. Figure 3 shows the comparison results between the observed and reanalysis products for the average monthly maximum and minimum temperature datasets. The unit of measurement for temperature is degree centigrade.

The historical potential evapotranspiration (PET) datasets were compared to the reanalysis products. Figure 4 shows the results of the comparison between observed and reanalysis products of the average potential evapotranspiration. The unit of measurement for PET is mm.

Meteorological stations represent the area of the specific locations that are near and around the measured point. The spatial interpolation technique was used to estimate meteorological parameters for other neighboring locations. The Inverse Distance Weight (IDW) method was selected from the spatial interpolation methods available in ArcGIS. Figure 5 shows the spatial distribution of the climate variables over the valley.

3.2. Selection of the Future Climate Model

The selected projected GCMs were evaluated by their distribution type, statistical characteristics, and trends in the observed/historical data of average monthly, seasonal, and annual time intervals.

3.2.1. Identification of Distribution

The observed average monthly, seasonal (summer, winter), and annual precipitation data were identified to have Normal (NO), Normal (NO), Reverse Gamble (RG), and Normal (NO) distributions, respectively. The MRI-ESM2-0, CESM2-WACCM, and FIO-ESM-2-0 models for precipitation and the FIO-ESM-2-0 model for maximum temperature have top ranks among the selected CMIP6 models.

Table 6. Probability distributions of observed and CMIP6 products.

Parameter	Average Monthly		JJAS		MAM		Annual
	AIC	DT	AIC	DT	AIC	DT	AIC	DT
Precipitation
Observed	226.01	NO	372.61	NO	350.56	RG	375.12	NO
ACCESS-ESM1-5	245	LO	368.69	RG	323.79	SEP3	394.64	LO
BCC-ESM1	195.07	RG	336.31	RG	309.1	RG	344.17	RG
CAMS-CSM1-0	236.61	LO	386.66	RG	274.86	SN2	385.70	LO
CanESM5	198.03	SN2	324.01	ST2	316.99	SN2	347.12	SN2
CAS-ESM2-0	300.99	SEP3	405.43	SEP4	425.80	SN2	450.09	SEP3
CESM2-WACCM	228.04	NO	334.50	SEP3	363.33	RG	377.14	NO
CMCC-ESM2	218.37	RG	356.51	SN2	300.02	SN2	367.47	RG
FIO-ESM-2-0	212.07	NO	325.88	LO	328.30	RG	361.17	NO
MIROC6	260.69	LO	386.23	NO	365.01	RG	409.79	LO
MRI-ESM2-0	264.18	NO	382.45	NO	389.56	SEP3	413.28	NO
NorESM2-LM	226.41	SN2	338.30	RG	356.48	SN2	375.50	SN2
TaiESM1	225.76	GU	381.12	NO	316.49	RG	374.85	GU
Maximum Temperature
Observed	67.44	SHASH	55.41	NO	78.19	SHASH	49.74	SEP2
CMCC-ESM2	57.80	GU	98.29	GU	76.85	SEP2	63.42	GU
BCC-ESM1	54.75	SN2	82.04	GU	68.89	NET	32.45	SEP3
MRI-ESM2-0	102.55	SEP3	79.39	NET	129.3	RG	82.77	NO
CAS-ESM2-0	55.71	RG	66.93	RG	83.78	NO	41.10	LO
FIO-ESM-2-0	70.56	GU	74.67	NO	88.46	GU	49.46	NO
Minimum Temperature
Observed	−12.99	SEP2	28.37	LO	46.64	GU	27.92	LO
CMCC-ESM2	73.59	LO	71.91	NO	85.28	SEP2	64.45	LO
BCC-ESM1	51.69	NO	45.17	NO	74.41	NO	42.50	NO
MRI-ESM2-0	40.35	NO	11.94	GU	29.32	LO	−0.36	SN2
CAS-ESM2-0	36.35	NO	25.73	NO	49.65	NO	17.42	NO
FIO-ESM-2-0	56.56	PE	34.39	SEP3	68.57	GU	48.75	GU

DT—distribution type, GU—Gumbel, SHASH—Sinh–Arcsinh, NET—Normal Exponential t, ST—Skew t distributions, NO—Normal, LO—Logistic distribution, SN2—Skew-normal, SEP2 and SEP3—Skew Power exponential 1-4, RG—Reverse Gumbel, PE—Power Exponential distributions. The shaded color indicates that the CMIP6 model results have similar distributions to the observed datasets.

shows the results of observed and historical CMIP6 model distributions.

3.2.2. Trend Test Results

The trend test was also performed at three different levels for the distribution analysis. For the precipitation data series, the CanESM5 model has the same trend as the observed data in all four levels, while the ACCESS-ESM1-5 model has a better fit for three time levels except in JJAS. CS-ESM2-0 has the same trends as the observed data for maximum temperature datasets, while CAS-ESM2-0 and FIO-ESM-2-0 have similar trends as the observed/historical minimum temperature series for all time levels.

The significance was determined based on the p-values at a significance level of 0.05.

Table 7. Trend test of the observed and GCM outputs.

Parameter	Average Monthly		JJAS		MAM		Annual
	TTR	Dir	TTR	Dir	TTR	Dir	TTR	Dir
Precipitation
Observed	NS	-	NS	+	NS	-	NS	-
ACCESS-ESM1-5	NS	-	NS	-	NS	-	NS	-
BCC-ESM1	NS	+	NS	+	NS	+	NS	+
CAMS-CSM1-0	NS	-	NS	-	S	I	NS	-
CanESM5	NS	-	NS	+	NS	-	NS	-
CAS-ESM2-0	NS	+	NS	+	NS	-	NS	+
CESM2-WACCM	NS	+	NS	+	NS	-	NS	+
CMCC-ESM2	NS	-	NS	-	NS	+	NS	-
FIO-ESM-2-0	NS	+	NS	+	NS	+	NS	+
MIROC6	NS	+	NS	+	NS	+	NS	+
MRI-ESM2-0	NS	+	NS	+	NS	-	NS	+
NorESM2-LM	NS	+	NS	+	NS	+	NS	+
TaiESM1	NS	+	NS	+	NS	-	NS	+
Maximum Temperature
Observed	S	I	S	I	S	I	S	I
CMCC-ESM2	NS	+	NS	+	NS	+	NS	+
BCC-ESM1	NS	+	NS	+	NS	+	NS	+
MRI-ESM2-0	NS	-	NS	-	NS	+	NS	+
CAS-ESM2-0	S	I	S	I	S	I	S	I
FIO-ESM-2-0	NS	+	NS	+	NS	+	NS	+
Minimum Temperature
Observed	S	I	S	I	S	I	S	I
CMCC-ESM2	S	I	S	I	NS	+	S	I
BCC-ESM1	NS	+	NS	+	S	I	S	I
MRI-ESM2-0	S	I	S	I	NS	+	S	I
CAS-ESM2-0	S	I	S	I	S	I	S	I
FIO-ESM-2-0	S	I	S	I	S	I	S	I

NB: TTR—trend test result, Dir—direction, S—significant, NS—non-significant, I—increasing, + and – represent the sign.

shows the results of the trend tests of the observed/historical data and the individual CMIP6 models.

3.2.3. Model Performance Measurement

The model performance measures were also calculated at three different levels. For the precipitation data datasets, the TaiESM1 model has a better overall model performance at all time levels. The CMCC-ESM2 products have a better model performance with the observed datasets for maximum temperature, while MRI-ESM2-0 has a good overall performance with the historical data series for minimum temperature at all time levels. The correlation, RMSE, and MAE of the observed/historical data and the individual GCM outputs are given in

Table 8. The model performance results.

Models	Monthly			JJAS			MAM			Annual
	R	RM	MA	R	RM	MA	R	RM	MA	R	RM	MA
Precipitation
ACCESS-ESM1-5	10	4	4	10	3	3	6	3	3	11	4	4
BCC-ESM1	3	8	8	3	7	7	10	8	7	3	8	8
CAMS-CSM1-0	8	1	2	12	5	5	2	7	8	10	1	2
CanESM5	9	10	10	4	10	10	1	5	5	9	10	10
CAS-ESM2-0	2	12	12	2	12	12	4	12	12	2	12	12
CESM2-WACCM	1	7	7	11	8	8	5	1	2	1	7	7
CMCC-ESM2	7	6	6	6	6	6	7	9	9	8	6	6
FIO-ESM-2-0	12	3	3	7	4	4	8	2	1	12	3	3
MIROC6	5	11	11	8	11	11	11	11	11	7	11	11
MRI-ESM2-0	4	5	5	1	2	2	9	10	10	5	5	5
NorESM2-LM	11	9	9	9	9	9	12	6	6	4	9	9
TaiESM1	6	2	1	5	1	1	3	4	4	6	2	1
Maximum Temperature
CMCC-ESM2	2	1	1	1	1	1	3	3	3	3	2	2
BCC-ESM1	3	3	3	2	3	3	2	2	2	1	3	3
MRI-ESM2-0	5	4	4	4	4	4	5	4	4	5	4	4
CAS-ESM2-0	1	5	5	3	5	5	1	5	5	2	5	5
FIO-ESM-2-0	4	2	2	5	2	2	4	1	1	4	1	1
Minimum Temperature
CMCC-ESM2	2	5	5	3	5	5	5	5	5	3	5	5
BCC-ESM1	5	1	1	5	2	2	3	1	1	5	1	1
MRI-ESM2-0	3	2	2	1	1	1	4	2	2	4	2	2
CAS-ESM2-0	1	4	4	2	4	4	1	4	4	1	4	4
FIO-ESM-2-0	4	3	3	4	3	3	2	3	3	2	3	3

R is correlation, RM is the root mean square error, and MA is the absolute mean error. The integer numbers in the table represent the ranks of each model.

. Each model is scored a rank.

The final ranks for the correlation (R), root mean square error (RMSE), and mean absolute error (MAE) of the observed/historical data and the CMIP6 climate models are given in

Table 9. Overall ranks of statistical tests.

CMIP6 Model	R	RMSE	MAE	Sum Rank	Total Rank
Precipitation
ACCESS-ESM1-5	37	14	14	65	5
BCC-ESM1	19	31	30	80	7
CAMS-CSM1-0	32	14	17	63	3
CanESM5	23	35	35	93	9
CAS-ESM2-0	10	48	48	106	11
CESM2-WACCM	18	23	24	65	5
CMCC-ESM2	28	27	27	82	8
FIO-ESM-2-0	39	12	11	62	2
MIROC6	31	44	44	119	12
MRI-ESM2-0	19	22	22	63	3
NorESM2-LM	36	33	33	102	10
TaiESM1	20	9	7	36	1
Maximum Temperature
CMCC-ESM2	9	7	7	23	1
BCC-ESM1	8	11	11	30	3
MRI-ESM2-0	19	16	16	51	5
CAS-ESM2-0	7	20	20	47	4
FIO-ESM-2-0	17	6	6	29	2
Minimum Temperature
CMCC-ESM2	13	20	20	53	5
BCC-ESM1	18	5	5	28	2
MRI-ESM2-0	12	7	7	26	1
CAS-ESM2-0	5	16	16	37	4
FIO-ESM-2-0	12	12	12	36	3

.

The FIO-ESM-2-0 model has the top rank among the selected CMIP6 models for precipitation.

Table 10. Results of the overall model performance measures.

Model	PDF	MPM	Trend	Sum Rank	Total Rank
	Precipitation
ACCESS-ESM1-5	7	5	2	14	6
BCC-ESM1	6	7	2	15	7
CAMS-CSM1-0	7	3	2	12	5
CanESM5	7	9	1	17	8
CAS-ESM2-0	7	11	2	20	12
CESM2-WACCM	1	5	2	8	4
CMCC-ESM2	7	8	2	17	8
FIO-ESM-2-0	1	2	2	5	1
MIROC6	4	12	2	18	10
MRI-ESM2-0	1	3	2	6	2
NorESM2-LM	7	10	2	19	11
TaiESM1	4	1	2	7	3
	Maximum Temperature
CMCC-ESM2	2	1	2	5	1
BCC-ESM1	2	3	2	7	3
MRI-ESM2-0	2	5	2	9	5
CAS-ESM2-0	2	4	1	7	3
FIO-ESM-2-0	1	2	2	5	1
	Minimum Temperature
CMCC-ESM2	1	5	1	7	2
BCC-ESM1	3	2	3	8	4
MRI-ESM2-0	3	1	3	7	2
CAS-ESM2-0	3	4	3	10	5
FIO-ESM-2-0	1	3	1	5	1

shows the summary of probability distribution functions (PDFs), model performance measures (MPMs), and trend (MK-tests) for the selected variables from the CMIP6.

The FIO-ESM-2-0 model has the top rank among the selected CMIP6 models for temperature.

Table 11. Ranks for temperature dataset.

Model	Tmax	Tmin	Sum Rank	Final Rank
CMCC-ESM2	1	2	3	2
BCC-ESM1	3	4	7	3
MRI-ESM2-0	5	2	7	3
CAS-ESM2-0	3	5	8	5
FIO-ESM-2-0	1	1	2	1

shows the sum of ranks for maximum and minimum temperatures from the CMIP6 models.

Based on the selected multi-model selection criteria, the FIO-ESM-2-0 model has the best overall model performance for all climate variables.

3.3. Bias Correction and Future Scenarios

After bias correction of the results of each GCM, the average annual precipitation in the middle scenario (SSP2-4.5) shows a decrease of 4.4% and 13% in 2054 and 2084 respectively. Similarly, in the worst-case scenario (SSP5-8.5), in 2054, there is a decrease in precipitation by 4%, and in 2084, by 12.8%. The average annual maximum temperature under the SSP2-4.5 scenario increases by 1.5 °C in 2054 and by 2.1 °C in 2084. The average annual maximum temperature under the worst case (SSP5-8.5) scenario increases by 1.7 °C in 2054 and by 3.2 °C in 2084. In the middle scenario (SSP4.5), the average annual minimum temperature increases by 2.2 °C in 2054 and by 3 °C in 2084. The average annual minimum temperature under the worst case (SSP5-8.5) scenario increase by 2.6 °C in 2054 and by 4.3 °C in 2084. Figure 6 presents the projected precipitation in the medium and worst-case scenarios derived from FIO-ESM-2-0. In all projection scenarios, the unit of measurement for the monthly precipitation and potential evapotranspiration is mm, while for temperature, it is degree centigrade.

Figure 7 presents the probable monthly maximum temperature in the medium and worst-case scenarios derived from FIO-ESM-2-0.

Figure 8 presents the probable monthly minimum temperature in the medium and worst-case scenarios derived from FIO-ESM-2-0.

Figure 9 presents the probable monthly potential evapotranspiration in the medium and worst-case scenarios derived from FIO-ESM-2-0.

4. Discussion

The results of this study quantify climate change variability for better management of water resources under medium and worst-case scenarios. Reanalysis, satellite, and merged (satellite and reanalysis) climate data products were evaluated using historical observed datasets. The merged products of the Climate Prediction Center (CPC) have a similar distribution as the observed precipitation on seasonal and annual bases. The reanalysis products of maximum temperature from MERRA-2 and 20CRv3 have the same trend test results as the observed datasets at all time levels.

We observed that the MERRA-2 product overestimates the precipitation data and underestimates temperature datasets at all time levels. The ERA5 product has the best fit with the simulated evapotranspiration using the Hargreaves method. The 20CRv3 product is the best-fitted model for maximum temperature.

The timing and magnitude of projected future climate changes are uncertain because of different ambiguities in anthropogenic and natural conditions, and climate sensitivity can change the predicted results. The results of this study showed that the outputs of the FIO-ESM-2-0 CIMP6 model have a good overall ranking for both precipitation and temperature.

This study evaluates the climate change model’s potential for analyzing climate variables and their trends. The merging of satellite data with model products reproduces best the historical observed meteostation data. Climate change scenarios can provide information on how future human-induced factors are expected to alter the composition of the atmosphere and how this may affect global climate conditions. The projected precipitation derived from the FIO-ESM-2-0 CIMP6 model will decrease in the medium and worst-case scenarios, while the projected temperature will increase.

5. Conclusions and Recommendations

The outputs of CMIP6 are preferred over CMIP5 because of several improved climate projection scenarios. Several improvements have been provided by CMIP6 such as better resolution, the number of parameterizations involved, and considering various emission scenarios. The model evaluation criteria are important for appropriate model selection. The GCMs are at a coarse resolution, and spatial downscaling should be performed before analyzing the climate change effects at a watershed level. Appropriate bias correction methods are important to reduce errors in CMIP6 model results. Because of some uncertainty in the variables, the projected climate change trends and their variability may change over time in the study area. Based on the statistical analysis, the FIO-ESM-2-0 model has the best performance for both weather datasets under the selected SSP scenarios.

The average annual precipitation in the middle scenario (SSP2-4.5) shows a decrease of 4.4% and 13% in 2054 and 2084, respectively. The worst-case scenario (SSP5-8.5) shows a decrease in precipitation by 4% in 2054 and by 12.8% in 2084. The average annual maximum temperature under the SSP2-4.5 middle scenario increases by 1.5 °C in 2054 and by 2.1 °C in 2084. Under the worst-case (SSP5-8.5) scenario, it increases by 1.7 °C in 2054 and by 3.2 °C in 2084. In the middle scenario (SSP2-4.5), the average annual minimum temperature increases by 2.2 °C in 2054 and by 3 °C in 2084, while in the worst-case scenario, it increases by 2.6 °C in 2054 and by 4.3 °C in 2084.

The results of the predicted climate variables using the CMIP6 under SSP scenarios also predict an increase in annual average temperature in the future. Moreover, the seasonal amount of precipitation will decrease in the winter but slightly increase in the summer. This will be accompanied by an increase in the annual potential evapotranspiration rate. The conclusions of this research have crucial importance for impact studies of climate change scenarios. These results can support hydrological impact studies based on the CMIP5 projections and would benefit from updating the CMIP6 ensemble to obtain more confident estimations of future hydrological conditions [49].

We obtained satisfactory results for the predicted climate change variability using three model selection criteria. However, some limitations have been observed such as the spatial resolution of the CMIP6 model result, which is coarser. Some studies revealed that CMIP6 models showing the highest warming are unlikely to be representative of the real world, and CMIP6 projections of global surface temperature should not be exclusively relied on for policy-relevant decisions [50]. Several improvements should be made to avoid drawbacks, such as addressing cloud effects, carbon budgets, and net-zero emission objectives. If the model selection criteria were modified, the probability of selecting another CMIP6 model as the optimal one would be high, and the projected results would be different.

A decrease in the estimation errors of CMIP6 projection products will be obtained using daily time steps. Moreover, studies are needed to combine a daily basis future data series to select the appropriate models for future predictions.