The Spatiotemporal Variability of Temperature and Precipitation Over the Upper Indus Basin: An Evaluation of 15 Year WRF Simulations

: Investigating the trends in the major climatic variables over the Upper Indus Basin (UIB) region is difficult for many reasons, including highly complex terrain with heterogeneous spatial precipitation patterns and a scarcity of gauge stations. The Weather Research and Forecasting (WRF) model was applied to simulate the spatiotemporal variability of precipitation and temperature over the Indus Basin from 2000 through 2015 with boundary conditions derived from the Climate Forecast System Reanalysis (CFSR) data. The WRF model was configured with three nested domains (d01–d03) with horizontal resolutions increasing inward from 36 km to 12 km to 4 km horizontal resolution, respectively. These simulations were a continuous run with a spin ‐ up year (i.e., 2000) to equilibrate the soil moisture, snow cover, and temperature at the beginning of the simulation. The simulations were then compared with TRMM and station data for the same time period using root mean squared error (RMSE), percentage bias (PBIAS), mean bias error (MBE), and the Pearson correlation coefficient. The results showed that the precipitation and temperature simulations were largely improved from d01 to d03. However, WRF tended to overestimate precipitation and underestimate temperature in all domains. This study presents high ‐ resolution climatological datasets, which could be useful for the study of climate change and hydrological processes in this data ‐ sparse region.


Introduction
The Himalaya-Karakoram-Hindukush (HKHK) region has significant strategic importance due to large amounts of snow and glacial ice [1]. Often referred to as Earth's "third pole" [2] (to indicate its massive ice and snow storage comparable to the North and South poles), it provides freshwater supply for agriculture, energy production, and drinking purposes for 1.5 billion people living in its catchments [3]. The HKHK region is also hydrologically important due to the presence of major rivers such as the Indus, Ganges, and Brahmaputra.
The Indus Basin Irrigation System (IBIS) is one of the largest irrigation systems in the world, which supports not only over 90% of Pakistan's agricultural production [3], but also the majority of the water needs for Pakistan [4]. Over the catchment area up to the Tarbella Reservoir (known as UIB), the contribution from snow and glacier melt runoffs accounts for about 70-80%, and the Karakoram Mountains alone contribute to more than 50% of runoff [3]. Despite the hydrological importance, there is no consensus on the assessment of the stability of the glaciers in the region. Previous studies have shown conflicting results for glaciers located in the Karakoram ranges, concluding that they are stable, retreating, or even advancing (the so-called Karakoram Anomaly [5,6]). The Indus River Basin (IRB) has experienced increasing temperatures and features immense spatial variability in climate [7].
Because the Indus River is largely dependent on cryosphere melt, it is highly vulnerable to climate change [8]. The Global Climate Risk Index (CRI) 2018 [9] has ranked Pakistan in the top ten countries most affected by climate change. Global climate models (GCMs) are the basic tools to simulate and project climate change, but their low spatial resolution makes them unable to resolve mesoscale flow patterns [10]. Several studies have downscaled GCM simulations over the IRB (for example, [11][12][13]); they found that temperature will continue to increase, but precipitation will be highly uncertain over this region. In addition, GCMs do not have the capability to resolve the orographic precipitation in complex mountainous regions such as UIB. This uncertainty in projected changes in climatic variables, especially precipitation, causes substantial uncertainty in runoff estimations.
In addition, the UIB is a data-scarce region where very few rain gauges are installed. These gauges are unevenly distributed and primarily at low-elevation valley locations, raising concerns about their representativeness of higher elevation orographic effects [14]. The lack of sufficient meteorological observations is usually the most challenging for flood forecasting and water resources management over this region. Complex topography, coupled with the challenges of field study in this region, has led to considerable uncertainty in assessing glacial mass balance and even meteorological trends. Moreover, the available global reanalysis datasets are useful to assess the large scale flow patterns over this region, but their coarse resolution cannot effectively characterize the complex orography and other local dynamics [14,15].
Regional climate models such as the Advanced Research Weather Research & Forecasting model (ARW-WRF; hereafter WRF) can be applied to simulate climatic parameters in complex terrain at high resolution. High-resolution atmospheric modeling can fill the gap of observational data scarcity to advance the understanding of atmospheric variability. The WRF model has successfully been applied worldwide [16][17][18][19][20] and to the Tibetan Plateau [21][22][23][24][25]. Gao et al. [25] applied the WRF model to perform 32 years (1979-2011) of simulations over the Tibetan Plateau at 30 km grid cell resolution, but this resolution is still somewhat coarse for capturing orographic precipitation patterns [23]. Maussion et al. [14] also used the WRF model to perform 11 years (October 2000-September 2011) of simulations over the Tibetan Plateau at 30 km and 10 km grid cell resolution in the outer and inner domain, respectively. Norris et al. [15] emphasized the importance of finer grid cell resolution over this region to simulate climatic variables, especially precipitation. However, investigating the longterm high-resolution meteorological variability (at 4 km) over the Karakoram Mountain Range is lacking, which is very important for the glaciological and hydrological communities.
In view of the above, we used the WRF model Version 3.8.1 [26] to dynamically downscale 16 years (2000-2015) over the Indus Basin. Our purpose in this study is to investigate spatiotemporal variability in the climate over the UIB and to assess the ability of WRF to capture observed variations using an analysis of correlation, bias, and trend. Potential issues with robustness stemming from record length are captured in our results via rigorous significance testing on all trends. The WRF model has been configured with three nested domains (d01-d03) with the boundary conditions derived from the National Oceanic and Atmospheric Administration (NOAA) National Center of Environmental Prediction (NCEP) 6 hourly Climate Forecast System Reanalysis (CFSR) data [27], which has a 38 km horizontal resolution. This study provides an analysis of the WRF model's applicability over the complex terrain of the UIB at an annual and seasonal scale.

Study Area
This study had three domains with different scales and grid cell resolutions. The northern parts of the Indus River consist of high mountains, including the HKHK mountain ranges, whereas its southern parts consist of flat plains. The inner domain (d03) covered the whole Indus Basin ( Figure  1), but this study is focused on the Upper Indus Basin (UIB). The climate in the northern parts (UIB) is largely influenced by western disturbances, and this region receives most of the precipitation (i.e., up to 2000 mm) in the winter season [28,29]. However, the precipitation in the Lower Indus Basin (southern parts) is driven by the monsoon precipitation in the summer season, which ranges from 100 to 500 mm [30,31]. In the UIB, the mean monthly temperature varies from 2 to 49 °C, whereas in the Lower Indus Basin, the mean monthly temperature varies from 14 to 44 °C [30]. The elevation of this study region varies from 0 to 8550 m above sea level, which has been derived from a 90 m resolution digital elevation model (DEM) from the Shuttle Radar Topography Mission (SRTM) project. The UIB is a high altitude and complex mountainous area wherein the Karakoram Mountain Ranges and the Hindu Kush and Himalayan Mountains are located. Typically, three major factors may influence the climate over this region: the winter westerlies, the summer monsoon, and a highpressure system formed over high Asia [32]. Most of the precipitation over this region occurs in winter and spring and is influenced by winter westerly flow [33].

WRF Model Configuration
The Advanced Research Weather Research & Forecasting model (ARW-WRF, hereafter WRF) Version 3.8.1 [26] was used to downscale 16 years (2000-2015) of CFSR data dynamically [27], which have approximately a 38 km horizontal grid resolution. These simulations were a continuous run with a spin-up year (i.e., 2000) to equilibrate the soil moisture, snow cover, and temperature at the beginning of the simulation. The motivation for using CFSR data in this study was taken from Bao and Zhang [34], wherein they evaluated several datasets over the Tibetan Plateau. They found CFSR and the Interim European Center for Medium-Range Weather Forecasts (ECMWF) reanalysis datasets (ERA-Interim) to simulate atmospheric changes effectively. They showed that both datasets presented smaller RMS error and mean bias. The main reason for selecting the CFSR dataset, specifically in this study, was its resolution. The CFSR dataset had a higher spatial resolution (0.5° by 0.5°) than ERA-Interim (0.75° by 0.75°).
The WRF model used in this study was configured with three nested domains with gradually increasing horizontal resolution from 36 km through 12 km to convection-permitting 4 km so that the innermost domain did not rely on a cumulus parameterization. The model configuration presented in Figure 1 was specifically chosen to limit the influence of boundary conditions on the results by assuring large margins between the nested domains. In addition, the relaxation zone of five points used in this study was very small compared to the domain sizes. The choice of the innermost domain size stemmed from work by Norris et al. [15], who emphasized that a fine grid cell resolution was required to resolve orographic precipitation in this region. The detailed model strategy and parametrization schemes used in this study are given in Table 1.
Research showed that WRF was highly sensitive to the selection of the land surface model (LSM) and cloud microphysical scheme [15]. Based on limited computational resources, four simulation experiments (Table 1. Section C) were performed for the year 2004 with a combination of three cloud microphysical schemes (Thompson, Morrison, and Goddard) and two land surface models (Noahmulti-parameterization (Noah-MP) and (Community Land Model version 4 (CLM4)) over the UIB. The RMS error and Pearson's r between WRF, stations, and TRMM data were estimated at Skardu station for the year 2004. The results showed that the RMS error between the WRF and station data was lower in all three domains in Simulation 1. In addition, Pearson's r between the WRF and station data in Simulation 1 was slightly higher than the other three simulations. Therefore, the results suggested that out of the tested configuration, Simulation 1 (Thompson and Noah-MP) offered the best performance. This was also consistent with Norris et al. [15], wherein they performed sensitivity analysis over the same region for selected summer and winter days. Morrison and Noah-MP Simulation 3 Goddard and Noah-MP Simulation 4 Thompson and CLM4

Model Validation
We evaluated the WRF temperature and precipitation output with the station meteorological data collected from the Pakistan Meteorological Department (PMD) and Water and Power Development Authority (WAPDA) (details of the stations are given in Table 2). We evaluated the WRF temperature with the available data of ten stations (Table 2) for the time period of 2001 through 2010. Similarly, the stations' data for precipitation were not available for the entire simulation period. Because of the limited availability of station data, we also assessed the WRF precipitation output with the Tropical Rainfall Measuring Mission (TRMM) 3B42 Version 7 gridded precipitation data [42] at a monthly temporal scale. The TRMM 3B42V7 data were available on a 0.25° by 0.25° latitude-longitude grid at 3 hourly temporal resolution. TRMM 3B42V7 data were collected from remote sensing and adjusted based on the monthly gauge data. Despite its coarse resolution and other limitations, TRMM 3B42V7 is considered to be one of the reliable gridded precipitation datasets [4,15]. Krakauer et al. [4] compared different precipitation datasets with the available stations' data over the Indus Basin and found the TRMM dataset to perform best among the remote sensing datasets. Similarly, Ali et al. [43] evaluated the TRMM Multi-satellite Precipitation Analysis (TMPA) precipitation products (3B42V6, 3B42V7, and 3B42RT) with gauge stations over the Hunza Basin in Karakoram Mountainous range. They also found 3B42V7 reasonably better than the other two products. Several studies have also evaluated the WRF precipitation with the TRMM dataset (for example, [14,15]). The WRF precipitation output was evaluated by the root mean squared error (RMSE), percentage bias, and Pearson correlation coefficient (r), whereas WRF temperature output was evaluated by RMSE, mean bias, and Pearson correlation coefficient. The expressions of RMSE, percentage bias (PBIAS), and mean bias for n grid points or n stations are: where and represent model simulations and observed data, respectively.

Trend Analysis
In this study, we performed a trend analysis using the Mann-Kendall (MK) significance test [44,45]. It is a non-parametric test and less affected by the extreme values, which is also widely used to detect trends in hydrologic time series [46,47].

Results
In this section, we evaluate the simulated precipitation in the three domains using the PMD meteorological data and TRMM data. However, the simulated temperature was evaluated only using the PMD meteorological data.

Station-averaged Precipitation Trends
This section describes the extent to which WRF was accurate in reproducing the spatio-temporal variability of precipitation in the three domains. The annual and seasonal mean precipitation trends for WRF, TRMM, and station data for all domains are shown in Figure 2. In addition, the Pearson correlation coefficient was computed between WRF and both observed datasets.
There were twenty-three stations (given at Table 2) operated by PMD and WAPDA above Tarbella Reservoir, and WRF and TRMM data were interpolated at these stations. When WRF, TRMM, and gauge data at these stations were averaged and compared in d01, WRF tended to overpredict the total annual precipitation relative to TRMM (115%) and the gauge data (132%) (left column, Figure 2). In addition, the results showed that the winter season had more bias (PBIAS-206%; RMSE-214 mm) and the summer season had less bias (PBIAS-70%; RMSE-91 mm) to the stations' data in the outer domain. For seasonal precipitation variation from 2001 to 2015, WRF showed a positive correlation with TRMM and station data for all seasons (Table 3). Similarly, WRF annual precipitation was significantly correlated with TRMM and gauge data (p < 0.05). However, trend analyses for WRF, TRMM, and gauge stations located in d01 did not show any statistically significant trends except the summer season, where station data showed a statistically significant trend. Moreover, the annual TRMM showed a statistically significant trend. When WRF and TRMM data were interpolated at the twenty-three gauge locations in d02 and compared, WRF still overpredicted total annual precipitation compared with the TRMM (35%) and gauge data (45%) (middle column, Figure 2). In addition, the results showed that the winter season had more bias (73%) and the summer season had less bias (3%) to the stations' data in the middle domain. WRF was significantly correlated with TRMM and gauge stations (p < 0.05) in all seasons except spring, where WRF was positively, but not significantly correlated with the station data. However, trend analyses for WRF, TRMM, and gauge stations located in d02 did not show any statistically significant trends except the summer season, where station data showed a statistically significant trend. Moreover, the annual TRMM showed a statistically significant trend.
When WRF and TRMM data were interpolated at the twenty-three gauge locations in d03 and compared, WRF still overpredicted total annual precipitation compared with TRMM (15%) and gauge data (23%) (right column, Figure 2). In addition, the results showed that the winter season had more bias (39%) and the summer season has less bias (−10%) to the stations' data in the inner domain (Table 3). Moreover, WRF was significantly correlated with TRMM and gauge stations (p < 0.05) in all seasons except spring, where WRF was positively, but not significantly correlated with the station and TRMM data. However, trend analyses for WRF, TRMM, and gauge stations located in d03 did not show any statistically significant trends except the summer season, where station data showed a statistically significant trend. Moreover, the annual TRMM showed a statistically significant trend. The station-wise comparison of PBIAS (Table 4) also showed that most of the precipitation simulations were largely improved from d01 to d03 at all stations. Besham Qila is located at the lowest elevation, and it had a low PBIAS value (−13%) in the d03; the increase of the resolution did not make any significant difference in the bias value. However, the highest elevated station was Khunjerab, where the precipitation simulations were improved from d01 (bias 350%) to d03 (bias 251%), but this station also had one of the largest biases. Among five high-elevation stations (Khunjerab, Ziarat, Ushkore, Yasin, and Karimabad), Karimabad, Ushkore, and Khujerab stations had the largest bias in d03. Low-elevation stations showed overall smaller biases. Among four low-elevation stations (Besham Qila, Daggar, Kandia, and Phulra), only Kandia showed a high bias (239%), whereas the other three had smaller biases (less than 20%).

Station-Averaged Temperature Trends
This section examines the extent to which WRF was accurate in reproducing the spatio-temporal variability of temperature in the three domains. There were only ten stations' data available for the 10 year time period (i.e., 2001 through 20110). The WRF output was compared with the station data by performing a lapse rate correction to adjust vertically from the WRF grid point elevation to the station elevation, as detailed below (Figure 3). For all three domains, WRF underpredicted the average annual and seasonal temperatures in all seasons. WRF showed a positive correlation with the station data (Table 5) in all seasons. However, most of these trends for temperature were not statistically significant. The important feature was how the simulated temperature and bias changed as the resolution increased. The results showed that RMSE and mean bias were noticeably changed from d01 to d03 in all seasons. This showed that WRF tended to underestimate the 2 m air temperature in d01 (annual mean bias of −7.7°C), d02 (annual mean bias of −3.4°C), and d03 (annual mean bias of −0.24°C), respectively.   Figure 4 shows the change in the mean lapse rate with altitude in all three domains. The average lapse rate for d01, d02, and d03 across the simulated period was estimated to be −6.9167, −6.7863, and −6.8111°C/km, respectively. These lapse rates were used to interpolate the simulated temperature time series vertically to the actual heights of the stations. For each station, the difference between the station height and WRF height was computed, multiplied by the average lapse rate, and added to the simulated temperature values ( Table 6). The station-wise comparison of the mean bias (Table 4) showed that most of the temperature simulations were greatly improved from d01 to d03 at all stations. The results showed that Astore, Chillas, Gilgit, and Skardu stations had the lowest bias in d03. However, Bunji, Chitral, and Gupis had the highest bias in d03.

Discussion and Conclusions
The WRF model was applied to simulate the spatiotemporal variability of precipitation and temperature over the Indus Basin from 2000 through 2015 using boundary conditions derived from the CFSR reanalysis dataset. We found that most of the precipitation trends in WRF, TRMM, and station precipitation at gauge locations were not statistically significant, which was consistent with Ahmad et al. [48], Khattak et al. [49], and Norris et al. [23].
The critical aspect was how the simulated precipitation amount and bias changed as resolution increased. The precipitation simulations were significantly improved from d01 (36 km) to d03 (4 km). The annual PBIAS between WRF and stations' data was decreased from d01 (132%) to d03 (23%). These results seemed broadly consistent with the previous studies [14,15,23,50,51], which also showed a positive influence of increasing resolution on precipitation over complex regions. Moreover, the WRF was significantly correlated with TRMM and station data in most of the seasons (except spring season) in all domains. WRF was able to capture the precipitation trends in most of the seasons, but tended to overestimate precipitation in all domains, which may reflect an underestimation of precipitation by the observed datasets as discussed earlier. Many researchers have revealed that the existing observation dataset network including TRMM 3B42V7 underestimates precipitation amounts by 66% to 300% over the UIB [52][53][54]. The gauges in the UIB in Pakistan suffer from undercatch for winter precipitation, especially snowfall (Norris et al. [15]), which could cause a significant amount of wet bias in this season. However, precipitation occurs in the form of rain in spring and summer seasons over this region, which the rain gauges measure more accurately. This contributed to less discrepancy between WRF precipitation and station observations in the summer season. However, the WRF simulations for temperature were much better than those of precipitation. The mean bias was noticeably reduced with the increasing resolution. The WRF model showed an underestimation of 2 m air temperature in all domains and seasons, which was consistent with prior studies [23] over this region.
The higher resolution simulations show a significantly improved representation of the climatic variables. Therefore, the spatial grid spacing resolution of the inner domain (d03) is high enough to capture the orographic effects in this area. Moreover, the results show that both precipitation and temperature output are largely improved with the increase of the resolution.
Some future research possibilities arise from this discussion. One is that an in-depth sensitivity analysis with different physical parameterization schemes (microphysics, LSM, PBL, and cumulus) may be beneficial in this region. In addition, we used the WRF model to simulate the spatiotemporal variability of precipitation and temperature over the Indus Basin from 2000 through 2015 with boundary conditions derived from the Climate Forecast System Reanalysis (CFSR) data. The CFSR data have been shown to capture the trends overall. However, it would be potentially useful to assess the alternative global reanalysis data products (such as ERA-Interim) to derive the boundary conditions over this region. It is possible that some of those alternative global reanalysis datasets may produce even better results. In this study, we have used the station data and TRMM 3B42V7 dataset to validate the WRF output. However, employing multiple gridded datasets such as GPCC, CRU, and MODIS snow cover product to assess the WRF output would be recommended.
Knowing the limitations of GCMs over complex and rugged terrain, regional climate models were the best available tools to simulate climate change. This study did not intend to analyze the climatic processes that influence the precipitation and temperature, but provided an overview of the ability of WRF to simulate spatiotemporal variability of precipitation and temperature over the UIB. These analyses could provide a better understanding of the limitations and reliability of WRF simulations. In addition, this study presented high-resolution climatological datasets, which could be useful for climate change and other hydrological studies in this region.