Inﬂuence of Climate Change and Land-Use Alteration on Water Resources in Multan, Pakistan

: This study presents an evaluation of climate and land-use changes induced impacts on water resources of Multan City, Pakistan. Statistical Down Scaling Model (SDSM) and Geographical Information System (GIS) are used for climate change scenario and spatial analyses. Hydrologic Engineering Center’s Hydraulic Modeling System (HEC-HMS) model is used for rainfall-runoff simulation. The investigated results show signiﬁcant changes in climatological parameters, i.e., an increase in temperature and decrease in precipitation over the last 40 years, and a signiﬁcant urban expansion is also observed from 2000 to 2020. The increase in temperature and urbanization has reduced the inﬁltration rate into the soil and increased the runoff ﬂows. The HEC-HMS results indicate that surface runoff gradually increased over the last two decades. Consequently, the depth of the water table in the shallow aquifer has declined by about 0.3 m/year. Projected climate indices stipulate that groundwater depletion will occur in the future. Arsenic levels have exceeded the permissible limit owing to unplanned urban expansion and open dumping of industrial efﬂuents. The results can help an efﬁcient water resources management in Multan.


Introduction
In the context of climate change, the rapid urbanization, population growth, and anthropogenic activities strongly influence the local and global environments. Climate change has been considered the most alarming worldwide issue by climatologists [1]. According to the Intergovernmental Panel on Climate Change (IPCC) fourth assessment report, the global temperature has increased by 0.6-0.8 • C over the last century, which is a record increment over the last 1000 years [2,3]. Recent research stipulates that regular climate changes, intensification of anthropic activities, and rapid urbanization largely affect the water budget of urban areas [4][5][6][7][8]. Global warming has also changed the timing, amount, and distribution of water, causing adverse effects on the water supply system through flooding [9][10][11][12][13][14][15][16] and droughts risks [17]. Similarly, it has been demonstrated that land-use changes such as pavement, urbanization, logging, deforestation, and soil erosion also contribute to influence hydrologic processes such as rate of infiltration, interception, and evapotranspiration [18][19][20]. Due to an increase in impervious surface, the infiltration rate and concentration-time decrease, causing urbanization to yield more surface runoff and peak inflows to river systems [21][22][23].
Multan is among fast-growing cities in Punjab, Pakistan. Over the last 54 years (1960-2013), an average increase of 1.9 • C in the mean lowest temperature has been observed in Punjab province, with a projected increase of 1.4 • C in the following 30 years [24].
The study area is located in the urban hub of the district of Multan (Multan city), which is the pronounced part of Bari Doab, Punjab plains ( Figure 1) [34]. Geographically, it is located in the center of Pakistan, with an elevation of 114 to 135 m above sea level. The area comprises alluvial plains derived by the Indus, Chenab, Sutlej, Jhelum, and Ravi Rivers. The Sutlej and Indus Rivers define the southeastern and western boundaries of the Punjab plains, respectively. The slope of the Punjab plains is to the southwest and varies from 0.366 m per kilometer in the northern part to less than 0.182 m per kilometer in the southern part; the average slopes are 0.274 m per kilometer [34][35][36], and demonstrate that mixed and varied soil patterns exist in the area due to sudden and frequent changes in the streamflow velocities. Most of the soils in Punjab plains have been characterized as moderately to highly permeable, while few areas were found to have a low permeability [37,38].

Hydrological Framework
The study area lies on the eastern bank of River Chenab, which has become the only source of groundwater recharge for the area. This River typically consists of braided channels within the city of Multan. The elevation in the south parts is low and starts rising toward the middle and northeastern parts of the study region, where the elevation ranges from 126 to 135 m. This means that the natural surface slope is rising toward the center and northeastern parts. These features were further analyzed using Arc-Hydro tools to reveal the natural surface hydrologic setting of the study area as shown in Figure 2. In this figure, natural drainage points, drainage lines, fil, adjoining catchment, and catchment polygons are principally highlighted. Appl

Hydrological Framework
The study area lies on the eastern bank of River Chenab, which has become the only source of groundwater recharge for the area. This River typically consists of braided channels within the city of Multan. The elevation in the south parts is low and starts rising toward the middle and northeastern parts of the study region, where the elevation ranges from 126 to 135 m. This means that the natural surface slope is rising toward the center and northeastern parts. These features were further analyzed using Arc-Hydro tools to reveal the natural surface hydrologic setting of the study area as shown in Figure 2. In this figure, natural drainage points, drainage lines, fil, adjoining catchment, and catchment polygons are principally highlighted.

Data Collection
The observed daily predictors of National Centre for Environmental Prediction (NCEP), and daily predictors of Global Climate Model (GCM) were gathered from Canadian Government (http://climate-scenarios.canada.ca/?page=pred-hadcm3; accessed on

Hydrological Framework
The study area lies on the eastern bank of River Chenab, which has become the only source of groundwater recharge for the area. This River typically consists of braided channels within the city of Multan. The elevation in the south parts is low and starts rising toward the middle and northeastern parts of the study region, where the elevation ranges from 126 to 135 m. This means that the natural surface slope is rising toward the center and northeastern parts. These features were further analyzed using Arc-Hydro tools to reveal the natural surface hydrologic setting of the study area as shown in Figure 2. In this figure, natural drainage points, drainage lines, fil, adjoining catchment, and catchment polygons are principally highlighted.

Data Collection
The observed daily predictors of National Centre for Environmental Prediction (NCEP), and daily predictors of Global Climate Model (GCM) were gathered from Canadian Government (http://climate-scenarios.canada.ca/?page=pred-hadcm3; accessed on

Data Collection
The observed daily predictors of National Centre for Environmental Prediction (NCEP), and daily predictors of Global Climate Model (GCM) were gathered from Canadian Government (http://climate-scenarios.canada.ca/?page=pred-hadcm3; accessed on 20 September 2020). Historical daily observed maximum and minimum temperature, precipitation, and other climate variables such as air humidity and wind speed were collected from Pakistan Metrological Department (PMD). Data for depth to the water table and arsenic concentration level were collected from the Punjab irrigation department. Digital Elevation Model (DEM) and Landsat data were downloaded from the earth explorer, USGS (https://earthexplorer.usgs.gov; accessed on 20 September 2020) site. Figure 3 shows the methodological framework of this research. Firstly, Statistical Downscaling Model (SDSM) is used for climate downscaling. Past temperature and pre-cipitation data are used for SDSM calibration and validation using National Centre for Environmental Prediction (NCEP) and Global Climate Models (GCMs) predictors, respectively. Secondly, ArcGIS packages is implemented for hydrological modeling (watershed delineation), and spatial analyses of land-use land cover (LULC) changes. Thirdly, Hydrologic Engineering Center's-Hydrologic Modeling System (HEC-HMS) model is employed to simulate the runoff by using daily rainfall data from 2000 to 2020. Finally, results from the prior steps are used to evaluate the impact of climate and LULC changes on surface water (i.e., rainfall-runoff), groundwater (i.e., depth to water-table), and water quality (i.e., arsenic concentration level only) of the study area. senic concentration level were collected from the Punjab irrigation department. Digital Elevation Model (DEM) and Landsat data were downloaded from the earth explorer, USGS (https://earthexplorer.usgs.gov; accessed on 20 September 2020) site. Figure 3 shows the methodological framework of this research. Firstly, Statistical Downscaling Model (SDSM) is used for climate downscaling. Past temperature and precipitation data are used for SDSM calibration and validation using National Centre for Environmental Prediction (NCEP) and Global Climate Models (GCMs) predictors, respectively. Secondly, ArcGIS packages is implemented for hydrological modeling (watershed delineation), and spatial analyses of land-use land cover (LULC) changes. Thirdly, Hydrologic Engineering Center's-Hydrologic Modeling System (HEC-HMS) model is employed to simulate the runoff by using daily rainfall data from 2000 to 2020. Finally, results from the prior steps are used to evaluate the impact of climate and LULC changes on surface water (i.e., rainfall-runoff), groundwater (i.e., depth to water-table), and water quality (i.e., arsenic concentration level only) of the study area.

Description of SDSM
The SDSM is a coupled Multiple Linear Regression (MLR) and Stochastic Weather Generator (SWG) model developed by [39][40][41][42][43][44]. In this method, statistical and empirical relations between NCEP's predictors and observed predictands are generated by MLR throughout predictors screening and SDSM calibration processes [45][46][47][48][49]. In this study, the projected climate scenario (Table 1) was determined in the calibration phase through predictor screening. Let us recall that to screen the predictor variables, it is necessary to individually correlate them with the predictands, i.e., precipitation and temperature. Considering this aspect, a conditional sub-model was selected for precipitation to account for wet-day occurrences, whereas an unconditional process was adopted for independent temperature. As a result, predictor variables with a p range of (0 ≤ p ≤ 0.05) were chosen as the best fit correlated predictors for a single predictand. Moreover, using available observed daily data from 1978 to 2000, the model was calibrated and validated for Tmax and Tmin. However, for precipitation downscaling, the NCEP predictors were used for model calibration and GCM predictors for model validation using data of the same period 1990-2001 (details can be found in Section 3.1).

Description of SDSM
The SDSM is a coupled Multiple Linear Regression (MLR) and Stochastic Weather Generator (SWG) model developed by [39][40][41][42][43][44]. In this method, statistical and empirical relations between NCEP's predictors and observed predictands are generated by MLR throughout predictors screening and SDSM calibration processes [45][46][47][48][49]. In this study, the projected climate scenario (Table 1) was determined in the calibration phase through predictor screening. Let us recall that to screen the predictor variables, it is necessary to individually correlate them with the predictands, i.e., precipitation and temperature. Considering this aspect, a conditional sub-model was selected for precipitation to account for wet-day occurrences, whereas an unconditional process was adopted for independent temperature. As a result, predictor variables with a p range of (0 ≤ p ≤ 0.05) were chosen as the best fit correlated predictors for a single predictand. Moreover, using available observed daily data from 1978 to 2000, the model was calibrated and validated for T max and T min . However, for precipitation downscaling, the NCEP predictors were used for model calibration and GCM predictors for model validation using data of the same period 1990-2001 (details can be found in Section 3.1).
Furthermore, the weather generator process was conducted to produce synthetic daily data by using calibrated output file (*.PAR-extension) and observed NCEP selected predictors variables. It allows the verification of calibrated model using observed data and generated artificial time series for the present status of climate [50]. Lastly, the scenario generation was performed following a procedure analogous to the weather generation process. The only difference between weather generation and scenario generation processes is the H3A2a_1961-2099 file that was selected in the GCM directory for scenario generation (instead of the NCEP_1961-2001 file for weather generation). Scenario generation gives the projected climate predicted variables supplied by the global climate models. Addinsoft's XLSTAT Software was used to perform the statistical Mann-Kendall test. In climatologic and hydrologic time series, the Mann-Kendall test is a widely used statistical approach for trend analysis [51]. In this test, the null hypothesis H0 indicates that there is no trend in the time series and that the data are independent and randomly ordered. Conversely, the alternative hypothesis H1 assumes that there is a trend in the time series. In reality, the null hypothesis H0 is rejected in favor of H1 when p-value exceeds the alpha value (0.050) [51]. Kendall's parameters were used to analyze the trend of climatologic and hydrologic data. The positive or negative value of Statistic "S" was used to infer an upward or downward trend, while the statistic Kendall's tau was used to measure the relationship and correlation between two variables.

Supervised Classification of Satellite Imageries and Digital Elevation Model (DEM) Reconditioning
ArcGIS 10.5 was used for the supervised classification of satellite imageries. Specifically, after the image processing, atmospheric corrections, and radiometric calibration steps, the satellite images were used for supervised classification to categorize the study area into different LULC classes [50]. Similarly, Digital Elevation Model (DEM) reconditioning was performed using Arc-Hydro tools for watershed delineation. As a digital visualization of ground surface terrain, DEM contains several hydrological data from which can be extracted relevant surface information such as flow direction or slope, drainage line, and catchment polygons [52].

HEC-HMS Model Setup
HEC-HMS rainfall-runoff simulation model is used to calculate the total outflow of the small urban catchment. The HEC-HMS model simulates the runoff by using rainfall data and other input parameters [53]. The Soil Conservation Service-curve number (SCS-CN) method was used to find the losses for reach's input parameters. Basically, this method develops a relation between lag time (T lag ) and time of concentration (T c ) to calculate the lag time for each time scenario. In particular, it uses the curve number (CN), which is calculated from hydrological soil and land-use land cover data. In reality, an increase in curve number suggests a decrease in infiltration and an increase in runoff. Once the LULC classes are identified, the next step is to determine the percentage of each hydrological soil group (A-D) in the study area. Then, the curve number (in %) for each class is calculated and summed to give the final value of curve number. Finally, based on the length of the reach and characteristics of the subbasins, the Kirpich formula is used to find the time of concentration as follows: where S stands for the slope, L represents the Reach length (in ft), and I a is the Initial abstraction.
In this study, the urban catchment was divided into three subbasins along with their junction and reach for model input parameters, as shown in Figure 4. Rainfall real-time daily data from 2000 to 2020 were used to run HEC-HMS to simulate the total runoff from all subbasins of the urban catchment [53]. The start and end date of the storm event and the time interval were both provided in the manager component of "Control Specifications" in HEC-HMS model. A 1-day time interval was used for each five-year rainfall scenario from 2000 to 2020. For runoff simulation, the entire period from 1 January 2000 to 31 December 2020 is divided into four five-year scenarios. Hydrographs are obtained as simulation output at every unit (junction, reach, and outflow) of the established model; then hydrographs peaks are used to infer the maximum outflow discharge at junctions. data and other input parameters [53]. The Soil Conservation Service-curve number (SCS-CN) method was used to find the losses for reach's input parameters. Basically, this method develops a relation between lag time (Tlag) and time of concentration (Tc) to calculate the lag time for each time scenario. In particular, it uses the curve number (CN), which is calculated from hydrological soil and land-use land cover data. In reality, an increase in curve number suggests a decrease in infiltration and an increase in runoff. Once the LULC classes are identified, the next step is to determine the percentage of each hydrological soil group (A-D) in the study area. Then, the curve number (%) for each class is calculated and summed to give the final value of curve number. Finally, based on the length of the reach and characteristics of the subbasins, the Kirpich formula is used to find the time of concentration as follows:

1900√
(1) where S stands for the slope, L represents the Reach length (in ft), and Ia is the Initial abstraction.
In this study, the urban catchment was divided into three subbasins along with their junction and reach for model input parameters, as shown in Figure 4. Rainfall real-time daily data from 2000 to 2020 were used to run HEC-HMS to simulate the total runoff from all subbasins of the urban catchment [53]. The start and end date of the storm event and the time interval were both provided in the manager component of "Control Specifications" in HEC-HMS model. A 1-day time interval was used for each five-year rainfall scenario from 2000 to 2020. For runoff simulation, the entire period from 1 January 2000 to 31 December 2020 is divided into four five-year scenarios. Hydrographs are obtained as simulation output at every unit (junction, reach, and outflow) of the established model, then hydrographs peaks are used to infer the maximum outflow discharge at junctions.    Figure 5 depicts the simulated results for average monthly temperature. It can be seen that these results are in good agreement with field observations. This substantiates the efficacy of the proposed model for temperature simulation and climate projection in the context of temperature change. However, the simulation and projection of precipitation were found slightly challenging for the proposed model. The reason is that during the conditional setup of SDSM for precipitation, the intensity of local precipitation depends on the wet-dry period, which in turn relies on regional predictors like atmospheric pressure and air humidity [48]. This also explains why previous studies [54] show low-performance simulation outputs for precipitation of different regions. Besides the exception of some peaks, the monthly and seasonal mean and maximum simulated precipitation are consistent with the trends in observed data ( Figure 6). tion were found slightly challenging for the proposed model. The reason is that during the conditional setup of SDSM for precipitation, the intensity of local precipitation depends on the wet-dry period, which in turn relies on regional predictors like atmospheric pressure and air humidity [48]. This also explains why previous studies [54] show lowperformance simulation outputs for precipitation of different regions. Besides the exception of some peaks, the monthly and seasonal mean and maximum simulated precipitation are consistent with the trends in observed data ( Figure 6).  Average precipitation (mm) Figure 6. SDSM calibration and validation using NCEP and GCMs; For mean and maximum precipitation.

Climate Change Scenarios
Mann-Kendall's test results show an increasing trend in the mean monthly temperature and a decreasing trend in the monthly total precipitation. The positive value of S (1975) indicates an upward trend for monthly mean temperature. The positive tau value (0.017) denotes a positive correlation between time and temperature ( Figure 7). On the

Climate Change Scenarios
Mann-Kendall's test results show an increasing trend in the mean monthly temperature and a decreasing trend in the monthly total precipitation. The positive value of S (1975) indicates an upward trend for monthly mean temperature. The positive tau value (0.017) denotes a positive correlation between time and temperature ( Figure 7). On the other hand, the negative value of S (−17169) indicates a downward trend for monthly total precipitation, whereas the negative value of tau (−0.156) suggests a negative correlation between time and precipitation. p-values for monthly temperature and precipitation were found smaller than the alpha value (0.050), rejecting the null hypothesis H0 in favor of hypothesis H1. The acceptance of hypothesis H1 is even substantiated by the trend in temperature and precipitation time-series data. Conversely, Kendall's time series analysis for seasonal temperature rejected the null hypothesis H0 for all seasons ( Table 2). J a n u a r y F e b r u a r y M a r c h A p r i l M a y J u n e J u l y

Climate Change Scenarios
Mann-Kendall's test results show an increasing trend in the mean monthly temperature and a decreasing trend in the monthly total precipitation. The positive value of S (1975) indicates an upward trend for monthly mean temperature. The positive tau value (0.017) denotes a positive correlation between time and temperature (Figure 7). On the other hand, the negative value of S (−17169) indicates a downward trend for monthly total precipitation, whereas the negative value of tau (−0.156) suggests a negative correlation between time and precipitation. p-values for monthly temperature and precipitation were found smaller than the alpha value (0.050), rejecting the null hypothesis H0 in favor of hypothesis H1. The acceptance of hypothesis H1 is even substantiated by the trend in temperature and precipitation time-series data. Conversely, Kendall's time series analysis for seasonal temperature rejected the null hypothesis H0 for all seasons ( Table 2).    The increase in minimum temperature indices was found superior to 2 • C over the past 40 years. This result is consistent with the work of [55] that conducted research on extreme temperature indices in Pakistan and found that the increase in mean minimum temperature in Punjab province was about 1.9 • C over the past 54 years. However, the increment in maximum temperature indices was less significant than that of minimum temperature since a relatively slight variation was observed. After the 2000s, 2010 was the warmest year, with a maximum temperature of 50 • C around the end of May. In the late 1990s, precipitation reduced to its minimum (83 mm/year) over the previous four decades, and droughts occurred in the region. Moreover, a decrease in air humidity and increased wind speed was observed.

Predicted Temperature and Precipitation
The analysis of GCM predictors reveals an increase in rainfall from March to June and a decline from July to February. The spring season (March-May) shows high precipitation for all three-time scenarios as compared with the other three seasons (winter, summer, and autumn). The predicted maximum monthly rainfall was found to vary between 20.42 and 31.50 mm. Although monthly variations differed from one-time scenario to another, the 2080s scenario was found dominant in monthly mean rainfall in March, April, June, July, October, and November. As depicted in Figure 8, the months of April and May 2060s show maximum precipitation. However, the difference in maximum monthly rainfall is not obvious among all scenarios. past 40 years. This result is consistent with the work of [55] that conducted research on extreme temperature indices in Pakistan and found that the increase in mean minimum temperature in Punjab province was about 1.9 °C over the past 54 years. However, the increment in maximum temperature indices was less significant than that of minimum temperature since a relatively slight variation was observed. After the 2000s, 2010 was the warmest year, with a maximum temperature of 50 °C around the end of May. In the late 1990s, precipitation reduced to its minimum (83 mm/year) over the previous four decades, and droughts occurred in the region. Moreover, a decrease in air humidity and increased wind speed was observed.

Predicted Temperature and Precipitation
The analysis of GCM predictors reveals an increase in rainfall from March to June and a decline from July to February. The spring season (March-May) shows high precipitation for all three-time scenarios as compared with the other three seasons (winter, summer, and autumn). The predicted maximum monthly rainfall was found to vary between 20.42 and 31.50 mm. Although monthly variations differed from one-time scenario to another, the 2080s scenario was found dominant in monthly mean rainfall in March, April, June, July, October, and November. As depicted in Figure 8, the months of April and May 2060s show maximum precipitation. However, the difference in maximum monthly rainfall is not obvious among all scenarios.   Figure 9 shows a regular increasing trend in projected seasonal and annual maximum temperatures. The results show that the temperature will increase throughout the three scenarios. A probable increase in simulated highest annual temperature is 3.15 • C, 4.25 • C, and 5.93 • C by the 2040s, 2060s, and 2080s, respectively. However, the maximum increment was observed in the summer season, where the mean maximum temperature reaches 39.98 • C, 40.06 • C, and 40.13 • C by the 2040s, 2060s, and 2080s, respectively. June is the warmest month for all three projected scenarios and shows a significant increase in maximum temperature indices. The highest temperature in this month extended to 53.65 • C, 55.04 • C, and 55.91 • C by the 2040s, 2060s, and 2080s, respectively, and lowest maximum temperature extended to 32.21 • C, 32.51 • C, and 32.84 • C by the 2040s, 2060s, and 2080s, respectively.

Land-Use Land Cover (LULC) Changes
The results of this study expand the work of Manzoor et al. [56] that conducted research on land-use changes in Multan city from 1992 to 2015. They found that the urban land covered which was 1800 ha in 1992, extended by 112.5% and reached 3825.9 ha in 2002. It should be recalled that Multan is the fifth largest city in Pakistan by population. Land-use land cover changes from 2000 to 2020 are shown in Figure 10. It was found that a decade later, i.e., from 2015, a considerable expansion occurred and the urban part increased significantly. The built-up area, that was about 130.03 km 2 in 2000, gradually expanded over the last two decades, resulting in an increase of 25.46 km 2 , to 155.49 km 2 in 2020. On the other hand, the low-vegetation cover decreased from 112.94 km 2 to 105.4 km 2 , and the dense vegetation reduced from 31.8 kkm 2 to 12.21 km 2 . An analogous decreasing trend was also observed in the past, with an agricultural area reducing rate of 0.15% between 1992 to 2002 [56]. Soil characteristics and variability in LULC summarized in Table 3 shows that the major land-use classes of the area are residential land, low and high vegetation, unused or bare soil, and water bodies. The water-covered area increased from 11.36 km 2 (in 2000) to 15.49 km 2 (in 2005) and then decreased to 9.52 km 2 in 2010. It expanded again to 21.08 km 2 (2015) and then finally reduced to 13.03 km 2 (2020) (Figure 10). This happened due to unusual flash floods in the upstream area of the Chenab River basin over the last couple of decades.
Appl. Sci. 2022, 12, 5210 10 of 19 Figure 9 shows a regular increasing trend in projected seasonal and annual maximum temperatures. The results show that the temperature will increase throughout the three scenarios. A probable increase in simulated highest annual temperature is 3.15 °C, 4.25 °C, and 5.93 °C by the 2040s, 2060s, and 2080s, respectively. However, the maximum increment was observed in the summer season, where the mean maximum temperature reaches 39.98 °C, 40.06 °C, and 40.13 °C by the 2040s, 2060s, and 2080s, respectively. June is the warmest month for all three projected scenarios and shows a significant increase in maximum temperature indices. The highest temperature in this month extended to 53.65 °C, 55.04 °C, and 55.91 °C by the 2040s, 2060s, and 2080s, respectively, and lowest maximum temperature extended to 32.21 °C, 32.51 °C, and 32.84 °C by the 2040s, 2060s, and 2080s, respectively.

Land-Use Land Cover (LULC) Changes
The results of this study expand the work of Manzoor et al. [56] that conducted research on land-use changes in Multan city from 1992 to 2015. They found that the urban land covered which was 1800 ha in 1992, extended by 112.5% and reached 3825.9 ha in 2002. It should be recalled that Multan is the fifth largest city in Pakistan by population. Land-use land cover changes from 2000 to 2020 are shown in Figure 10. It was found that a decade later, i.e., from 2015, a considerable expansion occurred and the urban part increased significantly. The built-up area, that was about 130.03 km 2 in 2000, gradually expanded over the last two decades, resulting in an increase of 25.46 km 2 , to 155.49 km 2 in 2020. On the other hand, the low-vegetation cover decreased from 112.94 km 2 to 105.4 km 2 , and the dense vegetation reduced from 31.8 kkm 2 to 12.21 km 2 . An analogous decreasing    Table 4 summarizes the results of runoff generation for different climate/land-use scenarios. For the four scenarios investigated, the amount of runoff gradually increased from 79 to 85 between the 2000s and 2020s. This substantial increase can be correlated with  Table 4 summarizes the results of runoff generation for different climate/land-use scenarios. For the four scenarios investigated, the amount of runoff gradually increased from 79 to 85 between the 2000s and 2020s. This substantial increase can be correlated with the rapid urbanization observed in Multan over the last few decades. The SCS curve number was found to have a direct relation with the rate of infiltration and runoff discharge. The correlation between runoff discharge and curve number is shown in Figure 11. It should, however, be noted that no discharge gauges are installed in the study area to calculate the observed runoff and confirm the HEC-HMS simulation performance. Nevertheless, the simulated results are in good agreement with previously published research, which considered the same LULC, meteorological and soil conditions [44]. Moreover, in the first five-year scenario, the outlet point achieved a maximum peak discharge of 133.3 (m 3 /s) on 25 July 2001. The second five-year scenario reached a maximum discharge of 310.8 m 3 /s on 08 August 2010. The third and fourth five-year scenarios yielded peak discharges of 162.6 m 3 /s and 959.6 m 3 /s on 09 September 2012 and 31 August 2017, respectively. Despite a slight difference in curve number values, similar results were obtained by Rangari et al. [44] for Hyderabad city, where 221.4 mm rainfall transformed into 590.4 m 3 /s. These results indicate that the increase in urbanization reduces the rate of water infiltration and increases runoff generation, which is one of the main causes of groundwater depletion in the study area. 310.8 m 3 /s on 08 August 2010. The third and fourth five-year scenarios yielded peak discharges of 162.6 m 3 /s and 959.6 m 3 /s on 09 September 2012 and 31 August 2017, respectively. Despite a slight difference in curve number values, similar results were obtained by Rangari et al. [44] for Hyderabad city, where 221.4 mm rainfall transformed into 590.4 m 3 /s. These results indicate that the increase in urbanization reduces the rate of water infiltration and increases runoff generation, which is one of the main causes of groundwater depletion in the study area.  Figure 11. Relation between SCS curve number and peak runoff discharge. In summary, using the climatic and Land-use land cover data in HEC-HMS model, it is found that the soil infiltration rate influences overall water balance in the area. Due to the decrease in the rate of infiltration, the groundwater is depleting over time. Extensive urban infrastructure, extreme climate events, and groundwater exploitation also contribute altering the hydrologic cycle. Most of the rainwater is transformed into runoff, while the remaining part evaporates due to arid climate, without recharging the aquifers: causing the depth to water level to drop day by day. Figure 11. Relation between SCS curve number and peak runoff discharge.

Rainfall-Runoff Simulation and Water Balance
In summary, using the climatic and Land-use land cover data in HEC-HMS model, it is found that the soil infiltration rate influences overall water balance in the area. Due to the decrease in the rate of infiltration, the groundwater is depleting over time. Extensive urban infrastructure, extreme climate events, and groundwater exploitation also contribute altering the hydrologic cycle. Most of the rainwater is transformed into runoff, while the remaining part evaporates due to arid climate, without recharging the aquifers: causing the depth to water level to drop day by day.

Coupling Impacts
Both climate and LULC changes play a significant role in water quality and water resources management. The analysis of the combined effect of climate and LULC changes on water resources in the city of Multan indicates three main deleterious impacts. These impacts must be suitably integrated in the water resources management of Multan, particularly in the view of the alarming, predicted climate and LULC changes. Impact 1: reduction of groundwater resource. It was found that the increase in urbanization comes along with the increase of runoff and a decrease in infiltration rate (natural groundwater recharge). The reduction of groundwater recharge in the region has been exacerbated by the decrease of precipitation which was found to have a decreasing trend, particularly between July and February. The direct consequence of the combination of these mechanisms is that the depth of the water table is also declining with time. According to Young et al. [57], the water table in Multan city dropped by 0.3 m/year. Shallow aquifer accesses are found at a depth of 9.1-12.1 m, but the water extraction for drinking purposes is not fit at this depth due to pollutants (See Impact 2). More generally, over the last few decades, the Punjab plains have been subjected to a significant decline in depth to water level due to the high expansion of groundwater pumping [58]. In Punjab, most of the crops are grown using groundwater, where more than 20% of the agricultural area is irrigated by groundwater and 55% of irrigated land receives groundwater and canal water. Punjab's groundwater discharge rate exceeds the recharge rate, resulting in a significant decline in groundwater levels in parts of the Punjab province, specifically in lower Bari Doab [57]. Results show that the depth to water table declined approximately by 3~6 m across the study area's western and eastern parts, respectively, from 2002 to 2014 [57,58]. We speculate that these variations are also influenced by the Monsoon season. In Pakistan, the Monsoon season brings appalling climate conditions, i.e., heavy rainfalls, sudden and flash floods, droughts, sea-level rise, ice retreat, and so forth, influencing the water resources and overall economy. Groundwater table depth was identified between 6 m to 12 m in 2002 and reached almost 9 m to 18 m up to 2014. Variations (rise or fall) in depth to water table for pre-monsoon and post-monsoon are shown in Figures 12 and 13 In the future, modern information technologies, i.e., artificial intelligence [59][60][61][62][63][64], global position system, and remote sensing, should be adopted to investigate the impact on water resources due climate and land-use land variations.     Impact 2: deterioration of water quality. This impact is concomitantly caused by the magnitude of land use and behavior of water cycle. The unplanned and high rate of urbanization in major cities of Pakistan has made it difficult to dispose of solid wastes properly. Multan is among those cities where solid wastes are dumped in the open place or along roadsides without considering water pollution issues. In fact, leachate from waste dumping sites typically possess physicochemical parameters, including color, electrical conductivity (EC), total dissolved solids (TDS), total suspended solids (TSS), biological oxygen demand (BOD), NH 3 -N, PO 4 -P, SO 4 2− , and Fe, beyond the WHO threshold [65]. Consequently, the contamination of surface and groundwater in the study area result from the absence of proper and planned disposal sites [65,66]. Infiltration of water containing complex leachate alters water chemistry in shallows aquifers. Further, the Chenab River, which is the primary source of aquifer recharge in the study area, was found to have a significant level of arsenic (As) contamination above the permissible value (WHO standard) in some areas ( Figure 14). This situation is alarming for the citizens of Multan that are overly concerned and already complained to the government about the arsenic-contaminated water [37,67] identified that arsenic leaching from phosphogypsum of Pak-Arab Fertilizer Factory is deteriorating water quality and affecting agricultural fields, from where it is moved to the food chain with severe health consequences. Water wells at the Nishter hospital Multan, Cantonment board, Metro Plaza, and Basti Khudadad show the highest level of arsenic concentration (i.e., exceeding 91 ppb) ( Figure 14) [58]. Appl. Sci. 2022, 12, 5210 15 of 19 dumping sites typically possess physicochemical parameters, including color, electrical conductivity (EC), total dissolved solids (TDS), total suspended solids (TSS), biological oxygen demand (BOD), NH3-N, PO4-P, SO4 2-, and Fe, beyond the WHO threshold [65]. Consequently, the contamination of surface and groundwater in the study area result from the absence of proper and planned disposal sites [65,66]. Infiltration of water containing complex leachate alters water chemistry in shallows aquifers. Further, the Chenab River, which is the primary source of aquifer recharge in the study area, was found to have a significant level of arsenic (As) contamination above the permissible value (WHO standard) in some areas ( Figure 14). This situation is alarming for the citizens of Multan that are overly concerned and already complained to the government about the arseniccontaminated water [37,67] identified that arsenic leaching from phosphogypsum of Pak-Arab Fertilizer Factory is deteriorating water quality and affecting agricultural fields, from where it is moved to the food chain with severe health consequences. Water wells at the Nishter hospital Multan, Cantonment board, Metro Plaza, and Basti Khudadad show the highest level of arsenic concentration (i.e., exceeding 91 ppb) ( Figure 14) [58]. Impact 3: depletion of surface water. Significant temperature rise has reduced the air and soil moisture and infiltration rate and has increased the rate of evapotranspiration and atmospheric CO2, resulting in water scarcity in the study area. The growing volume of atmospheric CO2 alters global climate systems, intensifying the hydrological regime with negative consequences on the local water system [68,69]. The historical and projected climate change trends, place Multan as one of the world's hottest cities. Climate and landuse land cover (LULC) variations have brought surface water resources deterioration and alteration in the water cycle in many other countries around the globe [70,71]. The discrete evaluation of these changes is difficult in the sense that results may vary from one case to another [72,73]. Climate change vulnerability in Pakistan was recognized by the National Water Policy (2018), the Climate Change Act (2017), and the National Climate Change Policy (2012). Water security has been the main concern under climate change, as well as policy measures ranging from developing awareness and additional water conservation storage [74,75].

Conclusions
The study presented an investigation on the impacts of climate change-land use variation on water resources of Multan city, Punjab, Pakistan, based on an evaluation of past Impact 3: depletion of surface water. Significant temperature rise has reduced the air and soil moisture and infiltration rate and has increased the rate of evapotranspiration and atmospheric CO 2 , resulting in water scarcity in the study area. The growing volume of atmospheric CO 2 alters global climate systems, intensifying the hydrological regime with negative consequences on the local water system [68,69]. The historical and projected climate change trends, place Multan as one of the world's hottest cities. Climate and landuse land cover (LULC) variations have brought surface water resources deterioration and alteration in the water cycle in many other countries around the globe [70,71]. The discrete evaluation of these changes is difficult in the sense that results may vary from one case to another [72,73]. Climate change vulnerability in Pakistan was recognized by the National Water Policy (2018), the Climate Change Act (2017), and the National Climate Change Policy (2012). Water security has been the main concern under climate change, as well as policy measures ranging from developing awareness and additional water conservation storage [74,75].

Conclusions
The study presented an investigation on the impacts of climate change-land use variation on water resources of Multan city, Punjab, Pakistan, based on an evaluation of past climate and land-use changes. Multan, having an arid climate, is pronounced part of lower Bari Doab, which lies to the south of Punjab plains. The following conclusions can be summarized: (1) Results show significant variation in temperature indices, and a regular increment has been observed in minimum temperature indices, which increased significantly over the last 40 years. However, the increment in maximum indices was not regular and varied from year to year. A noticeable decrease in mean annual relative humidity and increase in wind speed were observed over the last four decades; (2) The urban land extended from 130.03 km 2 (45.44%) by 2000 to 155.49 km 2 (55.33%) by 2020. Conversely, the dense vegetation decreased from 31.80 km 2 (11.11%) by 2000 to 12.2 km 2 (4.27%) by 2020, while the light vegetation reduced from 112.94 km 2 (39.47%) to 105.40 km 2 (36.83%) by 2020. An increase in temperature and urbanization (impervious cover) has reduced the rate of infiltration and increased runoff inflows to the river; (3) The evaluated results from HEC-HMS revealed that the transition of rainfall into runoff increases with time. On 31 August 2017, the whole catchment basin's outlet discharge point recorded a maximum surface water runoff of 959.6 m 3 /s, with 306 mm of rainfall. This indicates that infiltration was minimal and that most of the rainfall was converted to surface runoff. Consequently, the groundwater is depleting rapidly due to an imbalance between discharge and recharge of the groundwater; (4) The depth of the water table is declining at the rate of 0.2 m per year. Water quantity and water quality are also deteriorating due to unplanned land-use changes and open dumping of domestic and industrial effluents. At some spots, arsenic concentration was found beyond the WHO threshold, creating health issues for the citizens of Multan. Therefore, climate and land-use changes are both deteriorating water distribution and water quality. Moreover, the effects of urbanization were found more prominent and more critical in the study area.