Reconstruction of Hydrometeorological Data Using Dendrochronology and Machine Learning Approaches to Bias-Correct Climate Models in Northern Tien Shan, Kyrgyzstan

Abstract

Tree-ring-width chronologies for 33 samples of Picea abies (L.) Karst. were developed, and a relationship between tree growth and hydrometeorological features was established and analyzed. Precipitation, temperature, and discharge records were extrapolated to understand past climate trends to evaluate the accuracy of global climate models (GCMs). Using Machine Learning (ML) approaches, hydrometeorological records were reconstructed/extrapolated back to 1886. An increase in the mean annual temperature (Tmean_a) increased the mean annual discharge (Dmean_a) via glacier melting; however, no temporal trends in annual precipitation were detected. For these reconstructed climate data, root-mean-square error (RMSE), Taylor diagrams, and Kling–Gupta efficiency (KGE) were used to evaluate and assess the robustness of GCMs. The CORDEX REMO models indicated the best performance for simulating precipitation and temperature over northern Tien Shan; these models replicated historical Tmen_a and P_a quite well (KGE = 0.24 and KGE = 0.24, respectively). Moreover, the multi-model ensembles with selected GCMs and bias correction can significantly increase the performance of climate models, especially for mountains region where small-scale orographic effects abound.

1. Introduction

Climate change presents a range of challenges for agriculture-, forest-, and water-management practices [1]. It affects the distribution of water resources in Kyrgyzstan and in greater Central Asia (CA), changes precipitation patterns in different seasons and regions, and affects water resources. As key climatic factors, precipitation and temperature play crucial roles in the water cycle [2]. Moreover, the frequent occurrence of extreme weather and climate events such as heatwaves, droughts, heavy rainfalls, and biotic/abiotic catastrophes during recent years are evidence of climate change [2,3,4,5,6]. Therefore, current and future water-related developments are major national concerns and potential issues in international relations. These issues raise the importance of evaluating the magnitude of change in precipitation and temperature [7].

Kyrgyzstan is the third-most-vulnerable country in Eastern Europe and Central Asia to climate change, mainly owing to its climate-sensitive agricultural systems and a lack of an adaptive capacity [8]. Kyrgyzstan is threatened with glacier melting and a lack of freshwater balance, which are accelerated by global warming [9]. Due to climate changes, frequent natural hazards, strained water resources, and fragile ecosystems in Kyrgyzstan, it is imperative to project regional climate change based on emission scenarios for application to vulnerability, impacts, and adaption assessments. Some efforts have been devoted to climate change projections in CA using both global (GCMs) and regional (RCMs) climate models [10]. GCM simulations largely agree on regionwide precipitation increases in CA by the end of this century [11,12]. RCM simulations consistently show notable warming in the cold season in northern CA [12,13]. However, GCMs typically have too coarse a resolution to directly infer the climate of mountain areas at local scales, which are greatly influenced by complex terrain. Therefore, it is common to downscale over regions of interest using statistical techniques or nested RCMs [13]. This is especially true for mountainous Kyrgyzstan, where 94% of the terrain is at least 1000 m a.s.l. [14].

Multi-model ensembles (MMEs) are usually used to minimize the uncertainty of GCMs and generate outputs consistent with local conditions [15,16]. This is because GCMs are not free of biases and either underestimate or overestimate climate parameters, although the uncertainty is smaller in CMIP6 GCMs than other CMIPs [17]. Therefore, choosing a suitable set of GCMs is important to improve water-resource management and climate-change-impact assessments, especially in developing nations with limited human and computational resources [18,19]. Selected GCMs should replicate the observed climate in terms of spatial and temporal variability [20,21,22]. Generally, a preferable GCM ensemble is selected based on past performance and the envelope method [18,23].

Long historical records are required to select a suitable set of GCMs and do bias-correction. However, instrumental records in Central Asia are typically short and discontinuous. One possible solution is to develop proxy records of hydrometeorological data using tree rings that provide reliable annual records of long-term natural variability in discharge, temperature, and precipitation, extending well beyond the instrumental record [24,25]. Due to these advantages, tree rings are considered good proxies for the study of climate and hydrological changes, providing multi-century archives in many regions of the world [26,27,28,29,30,31]. The long-term discharge and hydroclimate information derived from these dendrochronological records can be used in water-resource planning and are important to supplement instrumental observations [32]. Previous dendrochronology studies in dry regions of Central Asia have focused mainly on reconstruction of climatic data, including temperature [33,34], precipitation [35,36], and mass balance of glaciers [37], mainly focusing on climate change. To date, only a few discharge reconstructions have been developed for the Tien Shan Mountains in China [38,39], Kazakhstan [40], and Kyrgyzstan [32,41,42]. Different statistical approaches have been applied in these hydrometeorological reconstructions using tree-ring data. Several studies employed Machine Learning (ML) methods to solve different earth and climate related problems [9,14,43,44,45,46]. Our research is the first to use ML methods in Kyrgyzstan to reconstruct hydrometeorological data by comparing different statistical-regression and ML models.

As noted, all GCMs suffer from some bias, which can be problematic for the assessment of climate scenarios, especially when dealing with extreme weather events [47,48,49]. Thus, we applied a bias-correction technique [13] to correct the climatology of the GCMs. Here, we used 33 new tree-ring-width chronologies from the Kashka-Suu River basin in the northern Tien Shan Mountains. We searched for the best relationship among tree-ring chronologies, annual precipitation, annual average air temperature, and annual average discharge, to develop a reliable 127-year mean annual-discharge reconstruction for the Kashka-Suu River, annual precipitation, and annual average temperature. Finally, we compared the hydrometeorological record with bias-corrected (BC) GCMs.

As such, the objective of this study is to show the importance of GCM spatial resolution and the new parameterization of the physical processes in a region of complex orography, using historical simulations from the Climate Model Intercomparison Project phase 5 (CMIP5), CMIP6, and Coordinated Regional Climate Downscaling Experiment (CORDEX). Moreover, we assessed the importance of GCM bias correction for this mountainous region.

2. Materials and Methods

2.1. Materials

2.1.1. Study Site

The Aksay cone (42°33′ N; 74°29′ E) is the largest debris-flow fan in Ala-Archa National Park (Figure 1) and in the Kashka-Suu River basin [42]. It is located on the northern slope of the Kyrgyz range, 35 km southwest of Bishkek. The catchment area is ca. 28.3 km², and elevation ranges from 4895 to 2250 m a.s.l. In the upper catchment (above 3600 m a.s.l.), two valley glaciers exist: (i) the Uchitel glacier with an accumulation area on the western slopes of the Semenov Tianshanskiy (4895 m a.s.l.) and Korona (4691 m a.s.l.) summits; and (ii) the Aksay glacier, with an accumulation area located on southwestern slopes of the Korona and the western slope of the Dvurogaya (4814 m a.s.l.) summits [42]. The terminus of Uchitel glacier is located at 3640 m a.s.l. and that of Aksay glacier is at an elevation of 3330 m a.s.l. [42].

Long-term hydrometeorological records exist at Baytik station (42°42.8′ N; 74°32.6′ E; 1300 m a.s.l.), located at roughly 10 km downstream of Aksay valley [42]. Annual records (covering the period of 1915–2013) report annual mean precipitation of 517 mm, ranging between 18.1 mm in January to 87.1 mm in May. Mean annual temperature is 6.5 °C, with mean maximum values in July (18.1 °C) and lowest monthly mean values in January (−4.8 °C) [42].

2.1.2. Data

In this research, we used annual precipitation (P_a) and annual mean air temperature (Tmean_a) from the meteorological station at Baytik, covering the period of 1915–2013; mean annual discharge (P_a), monthly average air temperature (Tmean_m), monthly maximum air temperature (Tmax_m), and monthly minimum air temperature (Tmin_m) from meteorological station Ala-Archa, covering the period of 1979–2013; and monthly mean discharge (Dmean_m) from the hydrological station at Baytik, for the period of 1928–2013.

We developed 33 new tree-ring-width chronologies for the Kashka-Suu River basin in the northern Tien Shan Mountains from Picea abies (L.) Karst trees, covering the period of 1886–2013. All observed data were collected by the Institute of Water Problems and Hydropower of the National Academy of Science of the Kyrgyz Republic (IWPH of the NAS KR) in 2014, funded by the Swiss National Science Foundation (SNF) in the framework of the DEFenCC (project no. 152301; Future Debris Flows and lake outburst floods in Tien Shan: possible impacts of projected Climate Change) project. Eight hundred trees damaged by debris flows were sampled together with 33 undamaged trees located between Baytik and Ala-Archa meteorological stations (Figure 1). In our assessment, we used only the 33 undamaged trees with sample depth shown in Figure 2, where the mean tree age was 65 years.

Eleven CMIP5 and CMIP6 GCMs were used in the comparison between the CMIPs phases over the northern Tien Shan Mountains, as recommended in [50].

Table 1. Description of CMIP5, CMIP6 GCMs, and CORDEX RCMs used in our study and their spatial resolution.

Institution, Country	Models-CMIP6	Resolution	Models-CMIP5	Resolution	CAS-CORDEX	Resolution
Australian Research Council Centre of Excellence for Climate System Science, Australia	ACCESS-CM2	1.87 × 1.25	ACCESS13	1.90 × 1.20	-	-
Beijing Climate Center, Beijing, China	BCC-CSM2-MR	1.12 × 1.12	BCC-CSM1.1-M	2.80 × 2.80	-	-
Institute for Numerical Mathematics, Russia	INM-CM5-0	2.00 × 1.50	INMCM4.0	2.00 × 1.50	-	-
Institute Pierre Simon Laplace(IPSL), France	IPSL-CM6A-lR	2.50 × 1.27	IPSL-CM5A-lr	3.70 × 1.90	-	-
Japan Agency for Marine-Earth Science and Technology (JAMSTEC), Japan	MIROC6	1.40 × 1.40	MIROC5	1.40 × 1.40	-	-
Max Planck Institute for Meteorology, Germany	MPI-ESM1-2-HR	0.94 × 0.94	MPI-ESM-MR	1.90 × 1.90	RegCM4-3.v5	0.44 × 0.44
	MPI-ESM1-2-LR	1.87 × 1.86	MPI-ESM-LR	1.86 × 1.87	REMO2015.v1	0.22 × 0.22
Meteorological Research Institute, Japan	MRI-ESM2-0	1.12 × 1.12	MRI-ESM1	1.10 × 1.10	-	-
Norwegian Climate Centre, Norway	NorESM2-MM	1.00 × 1.00	NorESM1-M	1.89 × 2.50	REMO2015.v1	0.22 × 0.22
Met Office Hadley Centre, UK	HadGEM3-G	1.00 × 1.00	HadGEM2-ES	1.25 × 1.87	REMO2015.v1	0.22 × 0.22
	-	-	HadGEM2-ES	1.25 × 1.87	RegCM4-3.v5.	0.44 × 0.44
National Centre for Meteorological Research, France	CNRM-CM6-1	1.00 × 1.00	CNRM-CM5	1.40 × 1.40	ALARO-0.v1	0.22 × 0.22

lists the GCMs and their basic description. CMIP5 introduced radiation concentration pathways (RCPs), which explored various potential greenhouse gas emissions scenarios. CMIP6 uses socioeconomic shared pathways (SSPs) for climate projection, which considers potential shifts in Earth’s environment and global economic and demographic shifts. Six RCMs from the CAS-CORDEX regional experiment (Table 1) were used in our research. The historical simulation and future projection of GCM monthly precipitation, Tas, Tmax, and Tmin were collected from data portals [51,52,53].

2.2. Methods

Width and density of annual rings of trees are formed, in particular, under the influence of current and past climatic conditions—e.g., temperature, air humidity, precipitation, components of solar radiation, and wind regime [41]. Based on multi-century proxies of annual tree-ring widths, since the early work of Andrew Douglass [54], attempts have been made to develop relationships between the time series of various natural processes and indicators of annual growth of tree rings.

Samples were prepared and measured following standard dendrochronological procedures [55]. Tree rings were first counted and then measured with a precision of 0.01 mm using a digital LINTAB positioning table connected to a Leica stereomicroscope and the TSAP Win software [56]. Tree rings from undisturbed trees were converted into width indices by standardizing raw data to detect characteristic common years [57]. In our work to reconstruct hydrometeorological data and validate the reconstruction, we followed standard dendrochronology procedures [58], including selecting and coring the trees, and then measuring ring width.

To reconstruct hydrometeorological data, various ML approaches were evaluated to describe the relationships between tree-ring chronologies and hydrometeorological parameters (

Table 2. Correlation coefficients for comparisons between the hydrometeorological observations and tree-ring widths. Coefficients in red bold are significant at p-level 0.05.

Station Type Name		Meteorological, Baytik		Meteorological, Ala-Archa						Hydro-Logical Gauging, Baytik	-
	Variable	Pa	Tmeana	Pa	Tmeana	Tmaxa	Tmina	Tmaxma	Tminma	Dmeana	Tree-Ring Width
MeteorologicalBaytik	Pa	1	−0.07	0.83	−0.21	−0.36	0.44	−0.2	0.22	0.04	0.47
	Tmeana	−0.07	1	−0.05	−0.05	−0.01	0.19	0.57	0.48	0.36	−0.26
MeteorologicalAla-Archa	Pa	0.83	−0.05	1	−0.21	−0.44	0.45	−0.3	0.09	−0.09	0.6
	Tmeana	−0.21	−0.05	−0.21	1	0.19	−0.12	0.02	−0.08	0.24	−0.26
	Tmaxa	−0.36	−0.01	−0.44	0.19	1	−0.13	0.22	−0.01	0.16	−0.46
	Tmina	0.44	0.19	0.45	−0.12	−0.13	1	−0.24	0.39	−0.24	0.46
	Tmaxma	−0.2	0.57	−0.3	0.02	0.22	−0.24	1	0.3	0.34	−0.39
	Tminma	0.22	0.48	0.09	−0.08	−0.01	0.39	0.3	1	0.29	−0.02
Hydrological Gauging Baytik	Dmeana	0.04	0.36	−0.09	0.24	0.16	−0.24	0.34	0.29	1	0.14
-	Tree-ring width	0.47	−0.26	0.60	−0.26	−0.46	0.46	−0.39	−0.02	0.14	1

). ML algorithms can analyze incoming information and look for explicit and hidden patterns in these data, so, thus, they represent an extremely powerful modeling approach that facilitates the reproduction of extremely complex dependencies. In the ML approach employed for the parameters, 30% of the data were allocated for the test and 70% for training. Meteorological data from Ala-Archa were not used because the dataset was too short (n = 34 year); we used all these data to train and test the ML approach.

The relationship between hydrometeorological data and tree-ring width was assessed using Pearson correlation analysis. All statistical procedures were evaluated at p < 0.05 levels of significance. We used linear regression and ML algorithms to estimate the dependent discharge, temperature, and precipitation variables from a set of potential tree-ring predictors. We used Artificial Neural Network (ANN), Hyperparameter Tuning (HPT), Xgboost Regressor (XgbR), Random Forest Regressor (RFR), K-Nearest Neighbors Regressor (KNN), Decision Tree Regressor (DTR), Lasso Regression (LaR), and Linear Regression (LR) [46]. Once the best model was fully evaluated, it was applied to the full period of tree-ring data to generate the reconstruction. Models were evaluated based on statistical methods for assessing errors between actual and reconstructed data, including root-mean-square error (RMSE) and mean absolute error (MAE) [14]. We compared values of RMSE and MAE with the various models (ANN, HPT, XgbR, RFR, KNN, DTR, LR, and LaR).

A Mann-Kendall Trend Test [59] was used to determine whether significant trends were present in various time-series data. Specifically, trends of annual air temperature, annual precipitation, and annual discharge at Baytik meteorological station were assessed. Statistical analysis to compare model outputs with observations is necessary. RMSE, Pearson correlation coefficient (r), and Nash–Sutcliffe model efficiency (NSE) are widely used to evaluate the goodness of fit between simulations and observations [2]; however, their suitability as metrics has been questioned [2]. The Kling–Gupta efficiency (KGE) considers the correlation coefficient, bias, and variability between simulations and observations, and this method has proven to be more efficient than the commonly used NSE, while clearly being time-sensitive [2]. In addition, KGE solves problems caused by interactions between these components, such as underestimation of variability. In a perfect model with no data errors, the KGE value is equal to 1. Therefore, for better comparison between the observations and GCMs, we applied the KGE approach [60]. KGE can be described as

$$\mathrm{KGE}=1-\sqrt{{(\mathrm{r}-1)}^2+{(\mathsf{\alpha }-1)}^2+{(\mathsf{\beta }-1)}^2}$$

(1)

$$\mathsf{\alpha }={\mathsf{\sigma }}_{\mathrm{s}}/{\mathsf{\sigma }}_{\mathrm{o}}$$

(2)

$$\mathsf{\beta }={\mathsf{\mu }}_{\mathrm{s}}/{\mathsf{\mu }}_{\mathrm{o}}$$

(3)

where r is the linear correlation between simulated and observed data, α is the ratio of standard deviation of observed and simulated data, β is the ratio of mean values of simulations and observations, s subscript represents the simulations, and o represents the observations.

In addition, the Taylor diagram was then applied to compare the climate model assessment studies [61]. Taylor diagrams show correlation statistics of multiple climate models, such as the Pearson correlation coefficient, RMSE, and standard deviation, and, thus, can be used to summarize model accuracy and consistency [62].

The bias-correction technique developed by Bruyère et al. [13] is applied in this study to correct the climatology of the GCMs and allow synoptic and climate variability to change. Formulas are as follows:

$$\mathrm{GCM}=\overline{\mathrm{GCM}}+{\mathrm{GCM}}^{’}$$

(4)

$$\mathrm{Obs}=\overline{\mathrm{Obs}}+{\mathrm{Obs}}^{’}$$

(5)

$${\mathrm{GCM}}_{bc}=\overline{\mathrm{Obs}}+{\mathrm{GCM}}^{’}$$

(6)

where overbar terms are the mean annual cycle, primed terms are perturbations from the mean annual cycle, and ${\mathrm{GCM}}_{bc}$ is the revised (bias-corrected) GCM output. Obs represents observations from meteorological stations.

3. Results

3.1. Relationship between Discharge, Precipitation, Temperature, and Tree-Ring Width

For Baytik meteorological station, we used the long-term observations; however, for Ala-Archa meteorological station and Baytik hydrological-gauging station, we recalculated some parameters. By using P_a, Tmean_m, Tmax_m, and Tmin_m, we calculated annual maximum air temperature (Tmax_a), annual minimum air temperature (Tmin_a), annual mean of the monthly maximum air temperature (Tmax_ma), and annual mean of the monthly minimum air temperature (Tmin_ma) for Ala-Archa meteorological station. By using Dmean_m, we calculated annual mean discharge (Dmean_a) for Baytik hydrological-gauging station. To analyze the relationship between the hydrometeorological-observation variables and tree-ring widths, Person correlation coefficients were computed (Table 2).

Our analysis revealed significant positive correlations between tree-ring chronologies and annual precipitation (P_a) at Baytik meteorological station (r = 0.47, p < 0.05, n = 98), Ala-Archa meteorological station (r = 0.60, p < 0.05, n = 34), and Baytik hydrological-gauging station (r = 0.14, p < 0.05, n = 85); significant negative correlations were found for Tmin_a and Tmax_ma at Ala-Archa meteorological station.

3.2. Hydrometeorological-Data Reconstruction and Validation

The evaluation metrics indicate that the best results were obtained for RFR, KNN with HPT, and XgbR with HPT (red bold metrics in

Table 3. Evaluation metrics amongst selected hydrometeorological reconstructed variables. Coefficients in red bold are significant at p-level 0.05.

Station, Variable, Unit	Metrics	LR	XgbR with HPT	RFR with HPT	KNN with HPT	LaR	DTR with HPT	ANN
Meteorological, Baytik, Pa, mm	MAE	70.8	60.9	69.4	59.7	70.6	77.8	70.6
	RMSE	89.2	77.5	86.1	75.9	89.3	101.4	89.3
Meteorological, Ala-Archa, Pa, mm	MAE	39.2	31.2	25.6	38.1	42.5	33.3	38.1
	RMSE	49.7	41.3	33.1	50.4	53.8	44.7	50.4
Meteorological, Ala-Archa, Tmaxa, °C	MAE	1.3	0.9	1.1	1.1	1.3	1.2	1.2
	RMSE	1.6	1.2	1.4	1.3	1.7	1.4	1.4
Meteorological, Ala-Archa, Tmina, °C	MAE	1.6	1.2	1.3	1.2	1.6	1.4	1.5
	RMSE	2.0	1.5	1.7	1.6	2.0	1.8	1.9
Meteorological, Ala-Archa, Tmaxma, °C	MAE	0.7	0.6	0.6	0.5	0.7	0.5	0.7
	RMSE	0.8	0.7	0.7	0.6	0.8	0.7	0.8
Meteorological, Ala-Archa, Tminma, °C	MAE	0.6	0.5	0.4	0.7	0.6	0.7	0.6
	RMSE	0.8	0.6	0.5	0.8	0.8	0.8	0.8
Hydrological Gauging, Baytik, Dmeana, m3 s−1	MAE	0.5	0.2	0.4	0.5	0.5	0.5	0.5
	RMSE	0.6	0.1	0.5	0.6	0.6	0.6	0.6

). Noteworthy as well is that XgbR with HPT requires more computer time to calculate, compared to RFR and KNN with HPT.

Based on the strong correlation between tree-ring chronologies and observed hydrometeorological variables (Table 2), the annual mean discharge from 1886 to 1928 was reconstructed for the Baytik hydrological-gauging station (Figure 3b) using KNN with HPT annual precipitation reconstructed from 1886 to 1915 for the Baytik meteorological station using XgbR with HPT (Figure 3a). All meteorological variables from 1886 to 1979 for the Ala-Archa station were reconstructed using KNN with HPT and XgbR with HPT.

Based on the Mann-Kendall Trend Test, there was no significant (p > 0.05) temporal trend in annual precipitation (Figure 3a). However, a positive trend was observed for average annual temperature (p < 0.05) and average annual discharge (Figure 3b). The increase in average annual discharge (p < 0.05) can be partially explained by the average annual temperature increase. Given that there was no significant temporal trend in annual precipitation, the increased annual discharge most likely emanated from melting of the Uchitel and Aksay valley glaciers. Here, we can see that the reconstructed data are slightly smoothed, but fluctuations occur within the error of the reconstruction.

3.3. Performance Evaluation of CMIP5, CMIP6 GCMs, and CORDEX RCMs

The performance of CMIP5, CMIP6 GCMs, and CORDEX RCMs in reproducing P_a and Tmean_a, based on KGE and RMSE metrics, is shown in

Table 4. Statistical summary of the comparisons between the CMIP5, CMIP6 GCMs, and CORDEX RCMs simulations and observations at Baytik and Ala-Archa meteorological stations. KGE: Kling–Gupta efficiency, RMSE: root-mean-square error, P_a: annual precipitation, and Tmean_a: mean annual air temperature; Baytik station metrics: standard deviation (SD) of the P_a –104 mm/year, SD of the Tmean_a −0.7 °C/year; Ala-Archa station metrics: SD of the P_a –99 mm/year, SD of the Tmean_a −0.5 °C/year.

	Models	Meteorological Station Baytik				Meteorological Station Ala-Archa
		Pa		Tmeana		Pa		Tmeana
		KGE	RMSE, mm/year	KGE	RMSE, °C/year	KGE	RMSE, mm/year	KGE	RMSE, °C/year
CMIP5	ACCESS1-3	0.15	113	0.12	1.7	0.03	113	−0.12	2.3
	BCC-CSM1-1m	0.00	160	0.19	1.4	−0.05	149	−0.98	5.2
	CNRM-CM5	−0.13	160	−0.73	1.6	−0.20	149	−0.62	2.1
	HadGEM2-ES	0.08	116	0.08	0.6	0.07	105	0.02	1.0
	INMCM4.0	−0.12	197	−0.05	4.2	−0.05	182	−0.26	1.0
	IPSL-CM5A-LR	−0.19	158	−0.25	5.5	−0.12	143	−0.37	1.6
	MIROC5	−0.01	136	−0.01	2.5	0.07	117	−1.41	6.3
	MPI-ESM-MR	−0.05	143	0.10	3.0	−0.02	127	0.11	1.4
	MRI_ESM1	−0.08	156	0.17	0.9	0.03	145	−0.69	4.2
	NorESM1-M	−0.11	184	0.17	1.0	−0.01	168	−0.82	4.1
CMIP6	ACCESS-CM2	0.13	112	0.06	1.1	0.03	128	−0.85	4.5
	BCC-CSM2-MR	−0.13	162	0.18	1.1	−0.06	147	−0.85	4.6
	CNRM-CM6	−0.16	306	−0.15	5.0	−0.17	302	−0.04	1.2
	HadGEM3-GC31-MM	−0.24	282	−0.16	5.0	−0.25	274	0.03	1.3
	INM-CM5-0	0.02	127	0.05	1.1	0.05	113	−0.76	4.5
	IPSL-CM6A-LR	0.11	112	−0.10	3.3	0.04	103	−0.20	1.3
	MIROC6	−0.19	177	−0.21	4.5	−0.22	165	−1.85	8.1
	MPI-ESM1-2-LR	−0.17	170	0.06	3.2	−0.17	153	0.08	1.2
	MRI-ESM2-0	0.10	223	−0.05	2.1	0.05	222	−1.13	5.7
	MPI-ESM1-2-HR	−0.11	147	0.19	1.1	−0.16	135	−0.39	3.3
	NorESM2-MM	−0.02	251	0.18	1.3	−0.05	242	−0.90	4.9
CORDEX	ALARO-0.v1_CNRM	−0.15	268	−0.38	4.6	−0.43	326	−1.09	4.8
	RegCm4-3v5_HadCEM2-ES	−0.11	298	−0.01	3.9	−0.22	308	0.01	0.9
	RegCm4-3v5_MPI-ESM-MR	−0.35	315	0.24	2.6	−0.39	321	−0.07	1.7
	REMO2015_HadGEM2-ES	0.06	131	−0.07	4.2	0.07	128	0.11	0.8
	REMO2015_MPI-M-ESM-LR	0.17	137	−0.33	3.9	0.24	131	−0.37	0.8
	REMO2015_NCC-NorESM1-M	−0.08	164	−0.03	3.6	−0.11	160	0.02	1.0

. Metrics are highlighted from best to worst, from dark green to red, respectively (Table 4).

The nodes with the highest correlation coefficients and the smallest values of the average absolute error, which were calculated from the ensemble of simulated fields, were selected for validation, rather than the model nodes closest to the meteorological station. Research on the verification of atmospheric models in mountain regions [63] shows that similar atmospheric circulation features are more important than selecting the closest station within complex terrain.

The results of statistical estimations show higher KGEs for CMIP6 models for both variables, indicating better performance than their older versions in CMIP5, except for a few models. Furthermore, the CORDEX models could replicate historical Tmen_a with KGE = 0.24 and P_a with KGE = 0.24. The upgraded models simulated historical P_a and Tmean_a better than their older versions, except for BCC-CSM2-MR, HadGEM3-GC31-MM, and MIROC6.

3.4. Bias Corrections of GCMs

By using Equations (4)–(6), we conducted bias correction of the CMIP6 GCMs. Then, statistical metrics were calculated (

Table 5. Statistical summary of comparisons between the CMIP6 GCMs with and without bias correction simulations and observations at Baytik meteorological station. KGE: Kling–Gupta efficiency, RMSE: root mean square error, P_a: annual precipitation, and Tmean_a: mean annual air temperature; Baytik station metrics: standard deviation (SD) of the P_a –104 mm/year, SD of the Tmean_a −0.7 °C/year.

	Models	Baytik Meteorological Station
		Pa		Bias Corrected Pa		Tmeana		Bias Corrected Tmeana
		KGE	RMSE, mm/year	KGE	RMSE, mm/year	KGE	RMSE, °C/year	KGE	RMSE, °C/year
CMIP6	ACCESS-CM2	0.13	112	0.19	112	0.06	1.1	0.06	1.0
	BCC-CSM2-MR	−0.13	162	−0.10	125	0.18	1.1	0.19	1.0
	CNRM-CM6	−0.16	306	−0.06	153	−0.15	5.0	0.13	1.1
	HadGEM3-GC31-MM	−0.24	282	−0.16	173	−0.16	5.0	0.12	0.9
	INM-CM5-0	0.02	127	0.04	118	0.05	1.1	0.06	0.9
	IPSL-CM6A-LR	0.11	112	0.11	110	−0.10	3.3	0.01	1.1
	MIROC6	−0.19	177	−0.10	165	−0.21	4.5	0.00	1.1
	MPI-ESM1-2-LR	−0.17	170	−0.10	168	0.06	3.2	0.19	0.9
	MRI-ESM2-0	0.10	223	0.15	146	−0.05	2.1	−0.01	1.0
	MPI-ESM1-2-HR	−0.11	147	−0.05	143	0.19	1.1	0.20	0.9
	NorESM2-MM	−0.02	251	0.06	117	0.18	1.3	0.19	1.0
Mean	MMEs	−0.06	188	0.01	139	0.01	2.6	0.10	1.0

), and a Taylor diagram of the CMIP6 GCMs was produced (Figure 4). Moreover, we calculated the means for multi-model ensembles (MMEs). All models we assessed demonstrated similar performance increases after bias correction.

RMSE for all simulations of P_a and Tmean_a decreased after bias correction (Table 5), especially for CNRM-CM6 (305 mm/year to 153 mm/year and 5 °C/year to 1.1 °C/year). As noted, MMEs are typically used to minimize GCM uncertainty and generate outputs consistent with local conditions. Evaluation metrics (Table 5) for MMEs improved with implementation of BC but still indicate low performance. To decrease uncertainty and errors, we calculated MMEs with only selected GCMs with BC to achieve optimal performance; for Tmean_a: BCC-CSM2-MR, CNRM-CM6, HadGEM3-GC31-MM, MPI-ESM1-2-LR, MPI-ESM1-2-HR, and NorESM2-MM; and for P_a: INM-CM5-0, IPSL-CM6A-LR, MRI-ESM2-0, and NorESM2-MM (Table 5). Statistical metrics for selected GCMs are shown in the Taylor diagram (Figure 4). After selection and application of BC methods, errors decreased for P_a (RMSE declined from 59.9 mm/year to 55.5 mm/year), correlation coefficients increased from −0.1 to +0.3 (Figure 4a), and Tmean_a RMSE remained at the same level of 0.9 °C/year, but correlation coefficients increased from 0.2 to 0.35 (Figure 4b).

We compared the interannual observations at Baytik meteorological station with the simulated climate data (Figure 5). Correlation analysis revealed significant positive correlations between observed and simulated monthly precipitation for all months except April, September, and November, at a 95% confidence level after bias correction and selection of the previously mentioned best models (Figure 5a). The same increase in significant positive correlation was found between the observed and simulated monthly mean air temperature (Figure 5b) for January, March, April, September, and October (95% confidence level) after bias correction and selection MMEs.

4. Discussion

Different GCM and RCM models have various degrees of uncertainty [23]. Therefore, it is important to select GCMs based on their performance in modeling the climate of an area to decrease the uncertainty of climate projections [15]. However, analysis and selection of suitable GCMs require long historical climate records. Given the short and discontinuous instrumental records in Kyrgyzstan, we reconstructed hydrometeorological data using tree-ring records from 33 Picea abies (L.) Karst. trees covering the period from 1886 to 2013. ML approaches showed better accuracy of data reconstruction compared to multiple linear regression model approaches. Of the 33 Picea abies (L.) Karst. trees used to train ML algorithms, only 5 undamaged trees older than 1914 were available for reconstruction of hydrometeorological data, presenting some inferential limitations (Figure 2).

The reconstructed data show increases in mean annual temperature and mean annual discharge (Figure 3b). However, the precipitation trends are not significant (Figure 3a), likely because increasing temperatures are melting glaciers, which contribute to a higher average annual discharge in the Kashka-Suu River.

We used the past-performance-assessment approach to choose a subset of GCMs (MME) to reduce uncertainties related to the individual model. Since the aim of this study was to determine which climate models can best simulate historical precipitation and temperature, our findings provide confidence for future research on climate projections and climate-change-adaptation strategies.

In contrast to earlier climate-model-comparison studies conducted in this region [2,11,12,62,64], our study used meteorological data from Baytik and Ala-Archa stations rather than observation data from the Climatic Research Unit (CRU) time series (TS) v4.01 [65], with a low spatial resolution (0.5° × 0.5°), which did not incorporate data from the meteorological stations used in our study. In Kyrgyzstan, only 7 of the 56 existing meteorological stations participate in that international exchange program.

We assessed the performance of GCMs of CMIP5, CMIP6, and CORDEX to simulate historical annual hydrometeorological data over northern Tien Shan (Kyrgyzstan), for the period from 1886 to 2013, using reconstructed data. The performance was evaluated for both the individual models and their MME means. In addition, the performance of individual models in replicating historical climate was evaluated using RMSE and KGE, considering their capability to simultaneously estimate three statistical indices (correlation, bias, and variability). The results revealed an improvement of most CMIP6 GCMs, compared to their counterparts in CMIP5, for simulating precipitation and temperature over the northern Tien Shan Mountains.

Improvement in model resolution from CMIP5 to CMIP6 improved the model’s performance. This better performance of some of the CMIP6 GCMs may be due to enhanced parameterization. Some models fail to reproduce precipitation and temperature over the northern Tien Shan Mountains. A possible common source of bias is the complex orography in Kyrgyzstan. However, the highest performance in simulating precipitation and temperature was obtained using the CORDEX models, which replicated historical Tmean_a with KGE = 0.24 and P_a with KGE = 0.24. This is due to the high spatial resolution of CORDEX–REMO models (0.22°). Thus, spatial resolution plays a key role for regions with complex orography, both for modeling atmospheric processes and for model validation. Moreover, in our study, the nodes with the highest correlation coefficients and smallest values of average absolute error, which were calculated from the ensemble of simulated fields, were selected for validation rather than those closest to the meteorological-station-model nodes.

As previously mentioned, MMEs are usually used to minimize uncertainty in GCMs and generate outputs consistent with local conditions. Evaluation metrics (Table 5) for MMEs improved with implementation of BC but still had low performance. To assess uncertainty and error, we calculated MMEs with only selected GCMs with BC and showed that the implementation of the suggested methods of BC reduced noises in GCMs.

CMIP6 models were used for bias-correction assessment, due to their improvements from earlier versions in terms of new future-forcing scenarios. Compared to the earlier socio-economic forcing-based RCPs, CMIP6 employs a framework of Shared Socioeconomic Pathways (SSPs). The new pathways include a bundle of future mitigation and adaptation actions on climate change, based on projected economic and social changes—e.g., population and social aspects, resources and economic development, ecosystems, and institutions [20]. However, the performance of the selected GCMs with BC should be improved by doing physical downscaling of the GCMs in such mountainous terrain. For example, RMSE (Table 5) values were reduced from 1.1 °C/year to 0.9 °C/year for selected GCMs but were still barely acceptable with the SD of the observed Tmean_a of 0.7 °C/year. It is known from previous studies that RMSE values that are less than half of the SD of the observed data may be considered acceptable for model evaluation [66,67]. In the future, we plan to conduct bias correction and use future CORDEX models with CMIP6-driven models as soon as they become available and conduct experiments with physical downscaling of the GCMs.

5. Conclusions

The key findings of this research are summarized as follows:

(1)

Instrumental observations in Kyrgyzstan are insufficient for assessing long-term climate and hydrological changes. Due to their annual resolution and sensitivity to climate, tree rings provide reliable proxies that can be used to extend instrumental records, as shown in our findings between climate and discharge changes in drylands of Kyrgyzstan. We also provide qualitative information on long-term hydrologic variability in the region that can inform water managers, stakeholders, and decision-makers.

(2)

ML algorithms that combined RFR, KNN with HPT, and XgbR with HPT performed best, and these were used to reconstruct hydrometeorological data in Kyrgyzstan for the first time.

(3)

Increases in the average annual temperature and mean annual discharge of the Kashka-Suu River were associated with more rapid glacier melting; however, precipitation did not significantly change with time.

(4)

The CORDEX models best simulated precipitation and temperature over northern Tien Shan. These successfully replicated historical Tmean_a (KGE = 0.24) and Pa (KGE = 0.24), due to their high spatial resolution (0.22°), indicating that spatial resolution plays a key role in complex mountain regions both for modeling atmospheric processes and model validation.

(5)

Multi-model ensembles with selected GCMs and bias correction significantly increased performance of climate models.

Since climate models are important tools for assessing future climate change effects, their accuracy must be demonstrated before they can be employed to estimate climate change scenarios. Accurate climate models can then provide the basis for climate change adaptation as well as disaster reduction and prevention in Kyrgyzstan. In the future, similar research in other climatic zones of Kyrgyzstan is planned for more effective implementation of climate change adaptation strategies.