Integrated SWAT-MODFLOW Modeling-Based Groundwater Adaptation Policy Guidelines for Lahore, Pakistan under Projected Climate Change, and Human Development Scenarios

Abstract

Urban aquifers are experiencing increasing pressures from climate change, land-use change, and abstraction, consequently, altering groundwater levels and threatening sustainable water availability, consumption, and utilization. Sustainability in such areas requires the adaptation of groundwater resources to these stressors. Consequently, this research made projections about future climate, land use, and abstraction, examines how these drives will affect groundwater levels, and then proposes adaptation strategies to reduce the impact on Lahore’s groundwater resources. The objectives are achieved using an integrated modeling framework involving applications of Soil Water Assessment Tool (SWAT) and MODFLOW models. The results indicated a projected rise in T_min by ~2.03 °C and T_max by ~1.13 °C by 2100 under medium (RCP 4.5) and high-end (RCP 8.5) scenarios, respectively. Future precipitation changes for mid, near and far periods are projected to be −1.0%, 25%, and 24.5% under RCP4.5, and −17.5%, 27.5%, and 29.0% under RCP8.5, respectively. The built-up area in the Lahore division will dominate agricultural land in the future with an expansion from 965 m² to 3716 km² by the year 2100 under R1S1 (R2S2) land-use change scenarios (significant at p = 5%). The future population of the Lahore division will increase from 6.4 M to 24.6 M (28.7 M) by the year 2100 under SSP1 (SSP3) scenarios (significant at p = 5%). Groundwater level in bult-up areas will be projected to decline from 185 m to 125 m by 2100 due to increasing groundwater abstraction and expansion in the impermeable surface under all scenarios. In contrast, agricultural areas show a fluctuating trend with a slight increase in groundwater level due to decreasing abstraction and multiple recharge sources under combined scenarios. The results of this study can be a way forward for groundwater experts and related institutions to understand the potential situation of groundwater resources in the Lahore division and implement adaptation strategies to counteract diminishing groundwater resources.

1. Introduction

Groundwater resources in megacities of the developing world are susceptible to overexploitation and other stresses [1]. The multiple interactive triggers such as climate change, population, agriculture, and industry increase pressure on groundwater resources [2,3,4,5,6]. These pressures alter water table elevation and its temporal variations. Established literature highlights some areas in different regions around the world where human-induced climate change, land-use change, and abstraction are altering the groundwater attributes up to a level at which local systems are unable to cope with their negative impacts [7,8,9,10] and thereby affecting the sustainability of groundwater reserves.

Lahore is a metropolitan city in Pakistan where groundwater is the dominant source for domestic, agricultural, and industrial use. Due to overexploitation, Lahore is observing a fast groundwater level recession in certain areas [4,11,12]. Historical analysis of the climate reveals a rise in temperature in the region [13,14,15,16]. The available literature on water use highlights that more than 75% of users in Lahore have access to the piped water supply system (WSS) through direct connections. The estimated abstraction by users from sources other than the Water and Sanitation Agency (WASA) represents ~30% of the water consumed. Historical population records in Lahore show increasing trends from 6.3 M in 1998 to over 10 M at present. The city’s extensive growth has resulted in significant urban development [17], resulting in a rapid rise in tube wells, abstraction rates, and an increase in water table depth with time [11].

Although the groundwater system in Lahore is part of the groundwater reservoir of the Indus basin, exploitation has formed a depression in the local groundwater table, expanding east and southward. Ref. [4] studied the expansion of depression zones as a function of change in groundwater depth during five consecutive years (2007–2011). The results reveal an alarming situation in some areas where the water table has reached more than 38 m since 2007. The findings depict a gradual expansion in depression to the east and south directions due to an increase in the rate of abstraction and decreasing recharge. A study by [11] discovered the worst situation in a few areas of Lahore city where the water table lowered to about 45 m in depth. Synthesis of the existing literature highlights a gap between investigating studies and policy institutions and the lack of (modeling-based) studies on the integrated impact of multiple stresses (i.e., climate, land use, and abstraction). A handful of studies covering quality aspects of the local groundwater resources are readily available. However, quantity aspects are not fully covered [18,19,20]. Ref. [18] evaluated the impact of past climate change and abstraction on groundwater resources of Lahore with a focus on management issues. Ref. [20] used chemical, isotopes, and numerical techniques. They identified types and sources of recharge for the Lahore aquifer and simulated the impact of abstraction on groundwater levels up to 2019.

There are two basic techniques to groundwater modeling, the first of which is the volume-based approach that makes use of direct groundwater monitoring data from observational wells [21,22,23], while other one used physical-based hydrological models, e.g., 3D groundwater flow modeling [24], semi- and fully distributed hydrological models [25] and remotely assisted simple water balance models [2,25]. Volume based methods driven by observational data are generally more accurate, however, uneven, and sparse distribution of monitoring wells make this approach less useful particularly for Lahore study region. Physical-based hydrological models and remotely driven water balance approaches on other hand can be useful for data sparse regions to monitor the groundwater changes [24,25]. For instance, physical based models include the global scale hydrological model (GHMs: e.g., PCRaster Global Water Balance (PCR-GLOBWB)) [26], and regional scale models (e.g., Soil Water Assessment Tool (SWAT) [27], MODFLOW [28]). SWAT and MODFLOW applications in surrounding areas of study region have well documented in previous studies [2,29,30,31,32,33]. Most of the studies have employed modeling-based methodologies with a focus on groundwater use. However, almost all of them have focused on agricultural areas [2,25,33,34,35,36,37], while ignoring the impacts of Uran areas. A few other observation-based and quantitative studies are also part of the literature that points out the continuously deteriorating situation associated with management issues in Lahore city in context of historical records [4,38,39]. The literature is sufficient for understanding the severity of the problem for current conditions. However, future changes in groundwater levels, particularly in the urban and peri-urban areas of Lahore City, have not yet been investigated in the context of projected climate change and human activities. To achieve sustainability of the groundwater levels in an urban area, we need to develop a sound knowledge of the situation for the future. Consequently, knowledge development could help in proposing a concise groundwater management policy to guide structural and non-structural aspects.

Therefore, this study uses an integrated modeling framework consisted of SWAT and MODFLOW applications to investigate the groundwater changes in context of projected changes in climate change and anthropogenic drivers, e.g., land-use change, and abstraction on the groundwater resources in the Lahore division. We further provided a proposed formulation of a few adaptation options to counteract diminishing groundwater resources.

2. Materials and Methods

2.1. Study Area

The study uses a modeling approach involving SWAT for hydrologic modeling and MODFLOW for hydro-geologic modeling. Considering the different requirements of both models, two segments of the study area were formulated. The geographical distribution of the Ravi River basin (28,000 km²) comprises a 20% upstream area having rugged topography. Mountains in the foothills of the Himalayas form the topography of the upstream basin (4300 m above mean sea level (amsl)) that mostly remains covered with snow, while 80% of the basin area comprising the middle and downstream is plain. The Middle and downstream areas of the basin are agricultural lands [14,40]. River Ravi is one of the major tributaries of the Indus River. The river Ravi covers a total stretch of 720 km through the basin from its origin in the Kailash Mountains and drains into the Chenab River in the South-west area. Peek flows seasonally occur during monsoon months (July to Sept), while the rest of the months observe low flows. The river Ravi basin has a semiarid-tropical climate. Upstream areas receive heavier rainfall than downstream plain areas. The rainfall season mostly shows bimodal patterns, November to January are the low rainfall months, and July to September is heavy rainfall. The annual rainfall in the basin varies from 300 to 1200 mm, and mean temperature changes between 8 to 40 °C.

Lahore division lies between latitude 31°15′ N and 31°42′ N, longitude 74°01′ E and 74°39′ E and comprises of an area around 6800 km², covers four districts such as Lahore, Sheikhupura, Kasur, and Nankana Sahib. Geographically, the Lahore division has a plane topography that changes from 200 to 210 m amsl and a general slope of 1:3000 towards the south and southwest. Based on the latest census Lahore division is home to over 19.4 million population. Groundwater is the dominant source of water supply for domestic (95%) and agriculture (39%) consumption. The Lahore division aquifer majorly comprises unconsolidated alluvial complexes. Silt, sand, and clay in varying proportions are the main constituents of the alluvial complex. These constituents consolidate to form a thick sedimentary complex of more than 400 m thickness. Clay is the dominant constituent of the alluvial complex with small quantities near the Ravi River and gradually increases with distance from the river [41] (Figure 1).

2.2. Climate Data Collection and Evaluation

Scanty observed data is an issue of the study area, therefore reanalysis and interpolated data on three major components; Tmax, Tmin, and rainfall were downloaded from Princeton University (PU)website: (http://hydrology.princeton.edu/data/; access date: January 2022), APHRODITE (APH): (http://www.chikyu.ac.jp; access date: January 2022) and Climate change Prediction Center-National Oceanic and Atmospheric Administration (CPCNOAA) (https://www.ncdc.noaa.gov/cdo-web/; access date: January 2022). Likewise, future climate data on Tmax, Tmin, and rainfall for two representative concentration pathways (RCP4.5&8.5) scenarios and seven regional climate models (RCM) were obtained from the CORDEX South Asia website: (http://www.cordex.org/; access date: January 2022) (

Table 1. Data required for assessment of the impact of future climate, land use, and abstraction on groundwater resources in Lahore.

Data Type	Duration	Resolution	Sources(s)
Meteorological Data
Observed Climate Data Rainfall, Tmax, and Tmin	1982–2015	Daily	Pakistan Meteorological Department (PMD) Chandigarh Meteorological Station (CMS)
Gridded Climate Data Princeton University forcing APHRODITE dataset Rainfall	1982–2015	0.5° × 0.5° (daily)	Princeton University (http://hydrology.princeton.edu/data/; access date: January 2022) Research Institute for Humanity and Nature (http://www.chikyu.ac.jp; access date: January 2022)
NOAA climate dataset Tmax and Tmin	1982–2015	Daily	NOAA’s National Centers for Environmental Information (NCEI) (https://www.ncdc.noaa.gov/cdo-web/; access date: January 2022)
Hydrological Data
River discharge and Canal discharge	2000–2014	Daily	Punjab Irrigation Department (PID)
Regional Climate Models (RCMs)	1982–2100	Daily (0.5° × 0.5°)	CORDEX (http://www.cordex.org/; access date: January 2022)
Spatial Data
Land use data Soil date	2007	1 km	Ref. [40] World Map (https://worldmap.harvard.edu/data/geonode:DSMW_RdY; access date: January 2022)
Digital Elevation Model (DEM)		30 m	(https://earthexplorer.usgs.gov/; access date: January 2022)
Aquifer lithology and hydraulic data			Groundwater division of Water and Power Development Authority (WAPDA)
Population data (Counts, density)	1998–2017		Pakistan Bureau of Statistics (PBS) (http://www.pbs.gov.pk/; access date: January 2022)

). The downloaded data was performance checked. Two GCD and four RCM with outstanding performance were shortlisted for climate change analysis and used as input to the SWAT model. The detailed methodology and performance results are reported in [14].

2.3. Hydrological and Hydrogeological Data Collection and Processing

The daily hydrological data on river flows and canal discharge for fifteen years (2000–2014) were obtained from the Punjab irrigation department (PID) which is passed through quality check before use. The hydro-geological data on groundwater levels measured in Jun and Sept every year was obtained from PID, aquifer lithology from CSIRO, hydraulic conductivity, specific storage, and specific yield was extracted from literature [37].

2.4. Spatial Data Collection and Processing

The land use map of the Indus basin developed for 2007 was provided by [40]. This map is one of the several inputs to the hydrological model (SWAT). Two other land use maps, one developed by the European space agency (ESA), cover the global domain from 1990 to 2015. The other map prepared by [42] (covers the Lahore district from 1999 to 2021 (1999, 2011, 2013, 2021, and 2035). Both maps were obtained and used for land-use change analysis. The ASTER based digital elevation model (DEM) data with spatial resolution of 30 m was obtained from the United States Geological Survey (USGS) website: (https://earthexplorer.usgs.gov/; access date: January 2022). The global soil data were obtained from the Harvard University data archive: (https://worldmap.harvard.edu/data/geonode:DSMW_RdY; access date: January 2022). The population data of four censuses from 1951 to 2017 was obtained from the Pakistan bureau of statistics.

2.5. Future Climate Change Projections

The three most relevant parameters of climate (Tmax, Tmin, and rainfall) were used, because of their significance in the climate system, their most often use in climate change studies, and their relevance to the current work. The temporal changes in Tmax, Tmin, and rainfall in the base period (1982–2014) were calculated using PU and CPC-NOAA datasets. Future projections (2020–2095) were made using four shortlisted RCMs under RCP4.5 and RCP8.5 scenarios. Annual and seasonal scale projections were made where the wet season covers months (Apr to Sept) that typically receive heavy rainfall. The dry season covers months (Oct to Mar) that receive low or no rainfall at all. Their significance was tested using a student t-test, and linear trends were fitted using a linear regression model [43].

2.6. Future Land Use Change Projections

Two land-use maps were used; one developed by the [44] hereafter referred to as map-1 and the second by [45], referred to as map-2. Map-1 classifies the area into eleven classes and map-2 into five classes. Considering the scope of the study, both maps were reclassified into three broad categories: Agriculture, built-up, and water. The land-use changes estimated using map-1 and map-2 were compared with each other. These land-use changes were used to develop linear regression models and future land-use change projections made using map-2 under two transition scenarios: R1S1-Business as usual and R2S2-Conservation. The R2S2 somehow includes the effect of policies by the current federal government in Pakistan, such as the construction of high-rise buildings and tree plantation campaigns. The study area is lacking in research on land-use change projections, with the available literature instead focusing more on historical changes. In this situation, although land use projections are based on simplified assumptions, they are still informative for interested individuals or institutions.

2.7. Projected Groundwater Projections

The groundwater abstraction in urban and agricultural areas of Lahore was estimated independently. For urban areas, the population projections were made beforehand by fitting multiple regression models on population data projected by [46] for some random years (2025, 2050, 2075, and 2100) under SSP1 and SSP3 scenarios developed by [47]. Based on population density, the study area was divided into different zones. Groundwater abstraction was calculated as a multiplication product of per capita water demand (m³·c⁻¹·d⁻¹) and total population. For agricultural areas, the groundwater abstraction was calculated using the water balance approach of SWAT [25]. This study assumes that seasonal abstraction is half that of annual and will remain constant during each season. Therefore, seasonal scale analysis is not part of this study.

2.8. Hydrological Modeling Using SWAT

The Soil and Water Assessment Tool (SWAT) model was used to calculate the recharge and abstraction. This model was developed by the USDA-ARS. It can simulate hydrological and biogeochemical cycles and impacts of stimuli such as climate and land-use changes on these processes at the watershed scale [48,49,50,51]. The SWAT model offers to include a wide range of components such as weather, hydrology, soil characteristics, crop growth characteristics, land management operations, and nutrients load and flows in the runoff. An insightful description of this model is available in [52]. SWAT simulates water balance using Equation (1).

$${\mathrm{SW}}_{\mathrm{t}}={\mathrm{SW}}_{\mathrm{O}}+{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{t}}({\mathrm{R}}_{\mathrm{day}}-{\mathrm{Q}}_{\mathrm{surf}}-{\mathrm{E}}_{\mathrm{a}}-{\mathrm{W}}_{\mathrm{seep}}\times {\mathrm{Q}}_{\mathrm{gw}})$$

(1)

where:

• SW_t final soil water content at time t in (mm H₂O)
• SW_o initial soil water content on day i in (mm H₂O)
• t time in (days)
• R_day precipitation on day i in (mm)
• Q_surf surface runoff on day i in (mm)
• E_a evaporation on day i in (mm)
• W_seep percolation and bypass flow leaving the bottom of soil strata (mm)
• Q_gw return flow on day i (mm)

SWAT calculates the surface runoff using the SCS curve number method. It is a function of land use, soil characteristics, and previous soil moisture conditions. The mechanism of the SCS curve number method has been given in Equation (2).

$${\mathrm{Q}}_{\mathrm{surf}}=\frac{{({\mathrm{R}}_{\mathrm{day}}-{\mathrm{I}}_{\mathrm{s}})}^2}{\left({\mathrm{R}}_{\mathrm{day}}-{\mathrm{I}}_{\mathrm{a}}+\mathrm{S}\right)}$$

(2)

where:

• Q_surf = runoff or excess rainfall (mm).
• R_day = precipitation on a given day (mm).
• I_a = initial abstractions such as surface storage, interception, and infiltration prior to runoff (mm) and
• S = retention parameter (mm).

$$\mathrm{S}=25.4\left(\frac{1000}{\mathrm{CN}}-10\right)$$

(3)

CN = curve number a function of land characteristics, varies between 0 to 100; 0 shows easy conversion of water to direct runoff and 100 shows difficult conversion.

SWAT Model Calibration, Validation, and Performance Evaluation

The model was manually calibrated for eight years (2000–2007) and validated for seven years (2008–2014) at Ravi Syphon and Shahdara gauges. Sensitivity analysis of the SWAT model parameters is required to obtain their effect on simulated components and characterize uncertainties in those components [53]. In total, 32 parameters were considered in the sensitivity analysis. The parameter values changed up to four levels (±10 and ±25%) to check their influence on flow and groundwater components. The three most widely used statistical parameters: the coefficient of determination (R²), the percentage bias (PBIAS), and the Nash-Sutcliffe efficiency (NSE), were considered in this study to evaluate the performance of the SWAT model. An insightful description of these parameters and performance rating for SWAT model evaluation is reported in [54].

2.9. Hydrogeologic Modeling Using MODFLOW

In this study, MODFLOW was used to simulate the impact of climate change, land-use change, and abstraction on groundwater levels. MODFLOW is a three-dimensional finite-difference groundwater flow model. Groundwater flow within the aquifer is simulated in MODFLOW using a block-centered finite-difference approach. The partial differential equation (Equation (4)) describes the groundwater flow in each grid.

$$\frac{\partial }{{\partial }_{\mathrm{x}}}\left({\mathrm{K}}_{\mathrm{xx}}\frac{\partial \mathrm{h}}{\partial \mathrm{x}}\right)+\frac{\partial }{\partial \mathrm{y}}\left({\mathrm{K}}_{\mathrm{yy}}\frac{\partial \mathrm{h}}{\partial \mathrm{y}}\right)+\frac{\partial }{\partial \mathrm{z}}\left({\mathrm{K}}_{\mathrm{zz}}\frac{\partial \mathrm{h}}{\partial \mathrm{z}}\right)-\mathrm{W}={\mathrm{S}}_{\mathrm{s}}\frac{\partial \mathrm{h}}{\partial \mathrm{t}}$$

(4)

In Equation (4), h is the hydraulic head (m) and an independent variable. The K_xx, K_yy, and K_zz represent the hydraulic conductivity (m.day⁻¹) in x, y, and z directions, respectively. The S_S shows specific storage that is a dimensionless quantity, while the letter W represents the source or sink. It is measured in day⁻¹. The W has two signs a positive sign shows the recharge/injection of water into the aquifer, a negative sign means abstraction. In this study, the groundwater flow model developed for Lahore also includes the surrounding agricultural areas. The geographical boundaries of the model area are presented in Figure 1. The model area was divided into uniform grids of 500 m spatial resolution in horizontal (X, Y) directions.

Based on the lithological data available at 149 locations, a total depth of 400 m was modeled with four layers of dominant materials. The model was calibrated and validated for transient conditions using the wet seasons from 2003–2014 with two stress periods each year: dry and wet having a duration of 182 and 183 days, respectively. The stress periods were further subdivided into 10-time steps with a time step multiplier of 1.2 to characterize the temporal variation in piezometric heads. The flow components were simulated using four MODFLOW modules.

MODFLOW Model Calibration and Performance Evaluation

In the study, a built-in Parameter Estimation (PEST) module was used to calibrate and validate the transient model for nine years (2003–2009). The PEST takes control of the MODFLOW and performs iterations as many times as necessary to determine the optimum parameter values. Annual piezometric water levels recorded by WASA were used as a reference for calibration. The aquifer parameters: specific yield, hydraulic conductivity, specific storage, and groundwater abstraction were adjusted for four aquifer layers in calibration to obtain the best-simulated heads. The coefficient of determination (R²) was used for performance evaluation.

2.10. Conceptualization of Scenarios Combination for Future Groundwater Level Projections

The study projects future climate, land use, and abstraction (population) under six combined scenarios (RCP4.5 and RCP8.5; R1S1 and R2S2; SSP1 and SSP3). The RCP and SSP are global-scale scenarios. These scenarios project the future evolution of climate and socio-economic development based on multiple variables such as socio-economic, technological, land use, energy, greenhouse gas emissions, and air pollutants. Unlike other studies [47,55] that combine RCP and SSP scenarios based on CO₂ emission projections, we used similar assumption to combine those scenarios based on the global population projections. Comparison of global population projections made under the SSP scenarios and reported in [56]. Background population projections for the RCP scenarios are reported in [57], revealing close agreements between SSP1 and RCP4.5 and SSP3 and RCP8.5 scenarios. The land-use change scenarios are applicable on the local scale. They were combined simultaneously with the climate and population scenario for the future groundwater level projections (

Table 2. Scenarios combinations for projected groundwater levels.

RCP Scenarios	Land Use Scenarios	SSP Scenarios	Scenario Combinations
RCP4.5	R1S1, R2S2	SSP1	RCP4.5-R1S1-SSP1 RCP4.5-R2S2-SSP1
RCP8.5	R1S1, R2S2	SSP3	RCP8.5-R1S1-SSP3 RCP8.5-R2S2-SSP3

).

3. Results and Discussion

3.1. Climate Change Projections

The mean climate will evolve in the future (2020–2095) as both the annual temperature and rainfall will increase under RCP4.5 and RCP8.5 scenarios, as shown in Figure 2a–c. Compared to the base period (1982–2005), the T_max will increase by 0.38 °C, 0.33 °C, and 0.64 °C in near, mid, and far future, respectively, under the RCP4.5, and 0.21 °C, 0.96 °C, and 2.04 °C in near, middle, and far future, respectively, under RCP8.5. Future Tma_x in the wet season (April to September) will decrease, except in the far future under RCP8.5, the expected absolute change will be −0.36 °C, −0.28 °C, and −0.07 °C in near, mid, and far future, respectively, under the RCP4.5 scenario, and −0.41 °C, −0.29 °C, and +0.45 °C in near, middle, and far future, respectively, under RCP8.5. Future Tma_x in the dry season (April to September) will increase, except near futures under RCP4.5, the expected absolute change will be −1.91 °C, +0.95 °C, and +1.36 °C in near, mid, and far future, respectively, under the RCP4.5 scenario, and +0.84 °C, +2.21 °C, and +3.65 °C in near, middle, and far future, respectively, under RCP8.5 (Figure 2a). The Tmin will increase. The absolute change in future T_min will be 0.66 °C, 0.94 °C, and 1.23 °C in near, mid, and far future, respectively, under the RCP4.5, and 0.85 °C, 2.05 °C, and 3.41 °C in near, middle, and far future, respectively, under RCP8.5. Future T_min in the wet season (April to September) will decrease, except in the far future under RCP8.5, the expected absolute change will be +0.73 °C, +0.70 °C, and +0.86 °C in near, mid, and far future, respectively, under the RCP4.5 scenario, and +0.81 °C, +1.73 °C, and +2.86 °C in near, middle and far future, respectively, under RCP8.5. Future T_min in the dry season (October to March) will increase, except near future under RCP4.5, the expected absolute change will be +0.59 °C, +1.19 °C, and +1.61 °C in near, mid, and far future, respectively, under the RCP4.5 scenario, and +0.89 °C, +2.39 °C, and +3.97 °C in near, middle, and far future, respectively, under RCP8.5 (Figure 2b).

Results show an increase in temperature in the future under the climate change scenarios. A comparison of seasonal scale projections shows that the dry season will be warmer than the wet season except for the near future, during which the T_max will decrease under the RCP4.5 scenario. Both Tmax and Tmin will increase while Tmin will rise more than the Tmax. Previous studies report similar results for the study area and other parts of the world [16,58,59,60,61]. Rising future temperatures will increase agricultural water demand due to an increase in evapotranspiration [62,63,64]. All-purpose domestic and industrial water consumption will also increase [63,65]. Consequently, rising demand will exacerbate pressure on groundwater resources.

Compared to the base period (1982–2005), the mean annual rainfall will increase in the future (2020–2095), except near future, under RCP4.5. Changes in rainfall in the near, mid, and far future will be −1%, 25%, and 24% under the RCP4.5, 17%, 27%, and 29% under RCP8.5. Future rainfall in the wet season will increase under RCP4.5 and RCP8.5 scenarios. The dry season will become drier under both climate change scenarios (Figure 2c). A comparison of the climate change scenarios shows agreement for the future projections, especially in the mid and far future. Rainfall projections are in close agreement with the previous studies [66]. Increasing future annual rainfall depicts more water availability which can offset negative pressure on groundwater resources. Since the increase in rainfall will be in the wet season, which is shorter in length than the dry season, dry seasons, due to longer duration and decreasing, rainfall will further exacerbate the pressure on groundwater resources in Lahore.

3.2. Land Use Change Projections

This section presents an analysis of the land-use change during the base period (2000–2014) and future period (2020–2095) under two scenarios: R1S1-Business as usual and R2S2-conservation in Lahore. The future projections for three broad urban land-use types are shown in Figure 3. Compared to the base period, the dominant land use (agriculture) will decrease from 5392 km² to 2862 km² under R1S1 and 3247 km² under the R2S2 scenario, while the built-up area will expand. The expansion in built-up will be from 1184 km² to 3716 km² and 3329 km² under R1S1 and R2S2, respectively. The change in agriculture and built-up areas will be significant across the entire century (

Table 3. Area under three broad land use and land cover types for the NF, MF and FF, projected under R1S1 and R2S2 scenarios.

Year	Land Use and Land Cover Type (km2)
	R1S1			R2S2
	Agriculture	Built-Up	Water	Agriculture	Built-Up	Water
Base period	5392	1184	240	5392	1184	240
2020	5048 *	1528 *	240	5001 *	1575 *	240
2043	4325 *	2251 *	240	4287 *	2289 *	240
2072	3642 *	2935 *	240	3736 *	2840 *	240
2100	2860 *	3716 *	240	3247 *	3329 *	240

Note. Each future time window was tested for significance relative to the base period. The p value for all future time windows was less than 0.05 significance level. * Shows a significant change (at 5%) in land use relative to the base period.

). The third land-use type, water, will remain the same with minor changes throughout the twenty-first century. An earlier study accomplished in the study area projected expansion in the urban area (built-up) [67]. The future contraction in Agricultural areas and expanding built-up will have implications for groundwater resources [6,33,68]. Besides putting pressure on groundwater, the agricultural (cultivated) area also plays a role in replenishing groundwater resources. Because a portion of the water received from rainfall and irrigation supplies percolates to groundwater, with the future contraction of agricultural areas, both the demand for irrigation and recharge will decrease [33]. The effect of the future expansion of built-up on groundwater resources will be more severe than the decrease in the agricultural area. It is due to a decrease in recharge and an increase in population-triggered water demand. Since the land use projections are driven by simplified assumptions (land use change (expansion or contraction) will only occur on the outskirts of existing landscapes). This assumption might overestimate land use change and its impact on groundwater resources.

3.3. SWAT Model Sensitivity Analysis

Nine parameters are found to be the most influential (Table 2). The CN2.mgt has the most influence and LAT_TIME.hru the least while the other seven parameters fall between the two. In terms of their relationship with simulated parameters, CANMX.hru, TLAPS.sub, ESCO.hru, and CH_K2.rte have a direct relationship with the surface flow components, while CN2.mgt, SOL_AWC.sol, and LAT_TIME.hru have an inverse relationship (Table 2). Only ESCO.hru has a direct relationship with groundwater components. CN2.mgt, TLAPS.sub, CH_K2.rte, and PLAPS.sub have no relationship at all with the groundwater components. Some parameters such as GW_DELAY.gw, SOL_K.sol and PLAPS.sub have no proper relationship with flow components. Besides their relationships, the effects of influential parameters were also observed during sensitivity analysis (Table S1) (supplementary material). For example, CN2.mgt and LAT_TIME.hru reduce instant peaks and smooth the hydrograph. TLAPS.sub and SOL_K.sol affect the occurrence (onset and withdrawal) of hydrograph peaks. Similarly, CANMX.hru and ESCO.hru affect the base flow and discharge while PLAPS.sub influences discharge only and SOL_AWC.sol affects soil moisture storage.

3.4. SWAT Model Calibration and Validation—For Discharge

Figure 4a–c shows a close matchup between observed and simulated hydrographs. Based on the statistics, the performance of the SWAT model varies from good to very good (

Table 4. SWAT model performance for calibration and validation at Ravi syphon and Shahdara gauges.

Performance Statistics	Calibration	Validation
Ravi Syphon gauge (upstream)
Coefficient of determination (R2)	0.77 (very good)	0.72 (good)
Nash-Sutcliffe efficiency (NSE)	0.76 (very good)	0.72 (good)
Percentage bias in volume (PBIAS)	−10.46 (good)	1.68 (very good)
Shahdara gauge (downstream)
Coefficient of determination (R2)	0.75 (very good)	0.81 (very good)
Nash-Sutcliffe efficiency (NSE)	0.74 (good)	0.75 (very good)
Percentage bias in volume (PBIAS)	−0.20 (very good)	0.79 (very good)

). The coefficient of determination (R²) has very good values; 0.77 for calibration and 0.72 for validation at Ravi Syphon, with 0.75 for calibration and 0.81 for the validation period at Shahdara gauge. The NSE also varies from good to very good; 0.76 for calibration and 0.72 for validation at Ravi Syphon with 0.74 and 0.75 at Shahdara gauge. The PBIAS only has a high value (−10.46) at Ravi Syphon for calibration, while it falls in a good range (1.68 and 0.79) at both stations for the validation period. The statistics for calibration and validation show that the model performance at both stations falls in a range from good to very good. The model can be used for further hydrological analysis.

3.5. SWAT Model Calibration and Validation—For Recharge

The SWAT model was calibrated for groundwater recharge as well. The calibration period spans five years (2003–2007) and the validation period seven years (2008–2014). Groundwater recharge, calculated using the water table fluctuation (WTF) method, was used as a reference. The optimized parameters are enlisted in

Table 5. List of sensitive parameters and their calibrated values for the swat model.

Parameter	Description of Parameters	Parameter Range	Calibrated Value	Sensitivity Rank
CN2.mgt	SCS runoff curve number	(35, 98)	47.63	1
CANMX.hru	Maximum canopy storage (mm H2O)	(0, 100)	6.02	2
TLAPS.sub	Temperature lapse rate (°C/km)	(−50, 50)	−3.80	3
ESCO.hru	Soil evaporation compensation factor	(0, 1)	0.35	4
SOL_AWC.sol	Available water capacity of the soil layer (mm H2O/mm soil)	(0, 1)	0.11	5
CH_K2.rte	Effective hydraulic conductivity in main channel alluvium	(−0.01, 500)	9.50	6
GW_DELAY.gw	Threshold depth of water in the shallow aquifer required for return flow to occur	(0, 500)	19.00	7
SOL_K.sol	Saturated hydraulic conductivity (mm/h)	(0, 2000)	42.22	8
PLAPS.sub	Precipitation lapse rate (mm H2O/km)	(−500, 500)	305.10	9
LAT_TIME.hru	Horizontal flow travel time (days)	(0, 180)	8.00	10

. The performance of calibrated SWAT model falls in the range of good to very good, with a minimum difference between annual simulated and reference groundwater recharge (Figure 4d,e). The calibrated model mostly overestimates groundwater recharge during validation. There is a considerable difference between the years 2009 and 2014 (25 mm and 18 mm, respectively). For the remaining years, the model matches well with the reference groundwater recharge. The coefficient of determination (R²) falls in a very good range for calibration (0.84) and validation (0.82), and the NSE falls in the class of satisfactory (0.60) for both calibration and validation periods. The PBIAS shows a very good value (2.0) for calibration and a high value (16) for validation. The SWAT model’s satisfactory performance has also been addressed by previous studies in the study area [2,10,25,30,69]. The success of the calibration of SWAT model depends on the choice of parameters and sensitivity [70]. Both factors largely depend on watershed characteristics and are determined by sensitivity analysis. Discharge in the Ravi River is controlled by a set of parameters that were determined during sensitivity analysis. The influence of selected parameters on discharge is described in Table S3. For the calibrated SWAT model, the CN2, CANMX, ESCO, SOL_AWC, and SOL_K were altered, depicting the dominant control of land surface characteristics and subsurface soil characteristics on river discharge. The calibrated CN2 value 47.63 depicts high infiltration. Small quantities of infiltered water are stored in the soil root zone, as indicated by low SOL_AWC (0.11). The possible reason is that much of the soil water is removed through evaporation, as shown by ESCO (0.35), and recharged to shallow groundwater, as depicted by GW_DELAY (19).

3.6. Groundwater Abstraction Projections

Groundwater abstraction in urban areas is dominated by domestic water use while in agricultural (rural) areas by water demand for agricultural consumption. This study considers only the dominant factors for each sector. Compared to the base period (2000–2014), annual groundwater abstraction will decrease in the 21st century under RCP4.5 and RCP8.5 scenarios (Figure 5). The reduction will be highest in the near future under RCP scenarios. At an annual scale, future abstraction will decrease by −17, −7, and −13% in the near future, mid future, and far future, respectively, under RCP4.5, and by −11, −9, and −10% in the near future, mid future, and far future, respectively, under RCP8.5. Future changes in abstraction will be significant under all scenarios. The groundwater abstraction follows the same pattern of change in both wet and dry seasons under RCP4.5 and RCP8.5 scenarios. Urban area abstraction will increase in the twenty-first century under SSP1 and SSP3 scenarios, as shown in Figure 6. Compared to the base period, groundwater abstraction will increase by 80% and 61% in the near future, by 158 and 135% in the mid future, and by 193 and 203% in the far future under SSP1 and SSP3, respectively. Groundwater abstraction increases during the twenty-first century, stabilizing in the far future under SSP1, while in contrast, under SSP3, it shows a continuous increase. The projected decrease in groundwater abstraction in agricultural areas is associated with increasing future rainfall that will enhance water availability for agriculture, and groundwater recharge and pressure on groundwater resources will decrease. Unlike agricultural areas, the likely rise in future abstraction in urban areas can be associated with increasing population. Previous studies [2,71] suggest pressure on groundwater resources is dominated by anthropogenic pumping. Rising future temperatures and expansion in the built-up area due to imperviousness will exacerbate the pressure on groundwater resources.

3.7. MODFLOW Sensitivity Analysis

The study used a parameter estimation package (PEST) for model calibration and sensitivity analysis. Sensitivity analysis of the parameters such as; specific storage, specific yield, vertical and horizontal hydraulic conductivities, and recharge was performed to check their influence on groundwater level. The groundwater model was most sensitive to horizontal hydraulic conductivity, followed by recharge, specific storage, vertical hydraulic conductivity, and specific yield.

3.8. MODFLOW Calibration—Steady State

The study used the parameter estimation (PEST) module to calibrate the steady-state MODFLOW model. The observed piezometric water levels for the wet season in 2003 (Dec-03) were used as reference targets. The specific yield (Sy), horizontal (Kh), and vertical (Kv) hydraulic conductivities, specific storage (Ss), and recharge were used as calibration parameters. The optimized hydraulic conductivity of the total formation varies from 66.53 to 89.52 mm.day-1 and vertical hydraulic conductivity from 3.47 to 6.02 mm.day-1. Specific storage and specific yield vary from 1.09E-04 to 2.47E-04 and 0.17 to 0.22, respectively. The calculated heads are given in Figure 7. Results depict that the plotted points are close to the best fit line with a small degree of scattering (less than 2 m), as shown in Figure 7. The degree of scattering varies between 2 and −2 m with an average value of −0.077 m (0.25 ft), thereby depicting overall good agreement between simulated water levels. The coefficient of determination is very high (R² = 0.9869), showing a high level of correspondence between observed and simulated water levels.

3.9. MODFLOW Calibration—Transient

The transient model was calibrated for seven years; December 2003 to June 2009. The horizontal and vertical hydraulic conductivities were not considered because hydraulic conductivity is a non-time variant parameter. Therefore, recharge, specific yield, and specific storage were optimized using the PEST model. The optimized parameters are listed in Table S2 (supplementary material). The average values of specific storage and yield for the total formation vary from 1.09E-06 to 2.47E-03 and 0.17 to 0.22, respectively. The calibrated groundwater flow model performed well in simulating the water levels, as shown in Figure 8. The simulated hydrographs showed well match with observed hydrographs by effectively capturing the seasonal variations and overall trends in urban and agricultural zones. The average bias of model-simulated water levels in urban zones varies between 0.59 m and 1.23 m, and in the agriculture zones between 0.52 m and 0.81 m. The falling water levels during the calibration period show increasing pressure on groundwater resources in urban zones.

3.10. Groundwater Level Projections

This section presents the future evolution of groundwater levels projected using the MODFLOW model for four scenarios. The simulations were performed to predict the behavior of groundwater resources up to the year 2100. All three factors, climate, land use, and abstraction, pose a combined effect on groundwater resources. Therefore, the following sections present the results, projected under combined climate change, land-use change, and abstraction scenarios. Changes in groundwater levels are a measure of groundwater storage, and therefore, this study focuses on groundwater levels only. Groundwater levels of urban and agricultural zones were projected using the MODFLOW up to 2100, as shown in Figure 9a–c. The results elucidate that groundwater levels will decrease at a much faster rate in urban areas. The drawdown in urban areas will vary from 45 to 55 m under all scenarios by 2100. The surroundings of urban areas (peri-urban) will observe a drawdown of about 9 to 10 m. In agricultural zones, there will be a drawdown of more than 2 m under the RCP4.5R1S1 scenario.

The changes in groundwater levels for urban and peri-urban areas are significant at p = 0.05 under both scenarios. Changes in groundwater levels for the left and right zones will be insignificant. The fast decline in groundwater levels for urban areas is due to urbanization, causing a rise in abstraction and a fall in recharge. Unlike the increase in total annual precipitation under climate change scenarios, groundwater levels will decrease. The decrease in water levels is due to the dominant effect of land-use change and abstraction, both of which marginalize the climate change effect, causing groundwater levels to reduce at a much faster rate. In peri-urban areas, the rate of decline is relatively slower. Low population density and recharge from rainfall and irrigation water supply slows the rate of decline. In agricultural zones, the fluctuation in groundwater levels is apparent. Unlike urban and peri-urban areas, groundwater levels will increase. Besides rainfall recharging, the dense network of irrigation, link canals, and irrigation return flows contribute to the likely increase in water levels.

3.11. Future Impact on Groundwater Resources

The impact of climate change, land-use change, and abstraction on groundwater levels was estimated under combined scenarios by subtracting the groundwater levels for the base period and future. An increase in depth was assigned a negative sign and vice versa. Results show that compared to the base period, the degree of impact will be much higher in the future, as shown in Figure 10 and Figure 11. Unlike north and west, most areas will observe negative impacts under combined scenarios in the near future.

In the north, the effect will change from positive to negative in the mid and far future. Urban and southeastern areas will observe the highest negative impact due to groundwater levels, decreasing at a much faster rate from 23 m in the near future to over 55 m in the far future. The surrounding areas being positively impacted in the base period will observe negative impacts due to the fast-outward expansion of the depression front. The groundwater levels in these areas are likely to decrease from 2 m in the near future to over 23 m by 2100. The agriculture areas, already facing a negative impact due to a decrease in groundwater levels during the base period, are likely to observe a further increase in the magnitude of negative impact. The likely decrease in groundwater levels in agriculture areas varies from 2 m in the near future to over 11 m in the far future. The degree of impact will be higher under RCP8.5-R1S1-SSP3 scenarios than RCP4.5-R1S1-SSP1. Similarly, the combined R1S1 scenario project higher changes than the R2S2 scenarios. The areas of a high negative impact and low adaptive capacity are likely to be more vulnerable in the future under all scenarios, as shown in Figure 10 and Figure 11, including urban areas. The highly vulnerable urban areas will expand outward and triple spatially in the far future. Besides experiencing a negative impact during the base period and increasing in the future, agricultural areas will be less vulnerable in the future, under all scenarios. The surrounding areas of urban settlements are also likely to be less vulnerable in the future, although margins will decrease due to their proximity to the urban area. The degree and range of vulnerability will be higher and wider under RCP8.5-R1S1-SSP3 scenarios than RCP4.5-R1S1-SSP1. Similarly, the combined R1S1 scenarios project a larger change than R2S2 scenarios.

4. Policy Guidelines for the Adaptation of the Impact of Multiple Stresses on Groundwater Levels

The results of this study reveal that the groundwater resources in Lahore are vulnerable to the three dominant stressors: climate change, land-use change, and abstraction. The main reason is the increased deficit over time due to rising water demand and a decreasing supply (recharge) to the Lahore aquifer. The situation will become severe in the future if current groundwater practices continue. Therefore, groundwater resources need to be protected against the impact of climate change (SDG13), land-use change, and abstraction to ensure the continuous supply of water in ample quantities (SDG6). Protection means complying with the United Nations Sustainable Development Goals: SDG6 (“clean water and sanitation”) and SDG 13 (“climate action”). The said goal can be achieved by formulating practical adaptation measures.

The indicators used in this study can be grouped into physical and climatic. The results also reveal that population (density) has a dominant influence on groundwater, followed by recharge and water table depth, respectively. Irrigation supply and impervious areas have minimum effect compared to other factors.

The projected results show quantitively that the population (in urban areas) of Lahore will increase many-fold in the future, ultimately increasing abstraction. At present, the annual groundwater abstraction from the Lahore aquifer is dominated by domestic consumption (52.9%) [39] and is likely to exacerbate due to the future increase in population. An increase in population is associated with local birth rates and migration from other parts of the country due to socioeconomic factors [73]. Expansion in built-up (urban) areas involves many factors, defining the rate at which it overruns other land-use types [74,75]. These factors are biophysical, infrastructural, and socioeconomic. The socioeconomic (66.67%) and infrastructural (64.10%) factors contribute equally and twice the number of biophysical factors (33.33%) in Lahore [45]. The current study projects that built-up area expansion will be consequent to the increase in the future population. An increase in built-up areas depicts the contraction of other land use types, especially agriculture areas, thus reducing the recharge from rainfall and agricultural return flows. In order to control the rapid increase in abstraction and decrease in recharge, some workable adaptation policy actions should be implemented, supported by the results of this study, especially those presented in

Table 6. Proposed adaptation options to counter the negative impacts on groundwater resources in Lahore.

Sr. No.	Improve	Adaptation Options	Time-Based Effectiveness	Approximate Time to Observe Outcome
1	A	Population control	Slow	10–15 years
2		New economic zones	Slow	10–20 years
3		Regulation of abstraction and zoning	Fast	3–5 years
4	A/R	Supplemental supply of treated sewage	Fast	1–2 years
5	R	Building development laws	Slow	5–10 years
6		Rainwater and storm water harvesting	Fast	4–5 years
7		River ponding	Fast	3–4 years

Note. A = Abstraction; R = Recharge. The time-based effectiveness of each adaption measure is defined based on experts’ experience and judgment.

. The detailed description of adaptation options is provided in Table S3 (Supplementary material).

5. Assumptions, Limitations, and Future Work

Land-use change, especially in the built-up category, is identified to be posing a severe impact on groundwater in Lahore. However, this study relies only on land use projections, developed based on simplified assumptions (future built-up will extend only around the existing built-up areas). Future studies can prepare future land-use change maps based on robust methodologies and replicate the results of this study for Lahore. Our study proposes some adaptation options that could improve the quantitative situation of groundwater in Lahore. Further research may evaluate the potential for groundwater replenishment of these adaptation options: river ponding, rainwater, stormwater harvesting, and reuse of treated wastewater using modeling-based tools. Because of the lack of observation-based data on groundwater quality, this study focuses only on the quantitative aspect of groundwater resources without the quality component. The groundwater quantity, as well as quality aspects, may be considered in future work and management policy formulation as well.

6. Conclusions

This study uses a multi-model integrated approach to investigate the impacts of climate change, land-use change, and abstraction on groundwater resources in Lahore. Following are the conclusions.

Future annual temperatures will rise, with the Tmin increasing more than the Tmax. The dry season will be warmer than the wet season under climate change scenarios. Future annual rainfall will increase, while patterns of annual and seasonal rainfall will remain the same.

The annual rainfall will increase by 24% and 29% under RCP4.5 and RCP8.5 scenarios, respectively, in the far future. The changes will be more apparent under the RCP8.5 compared to the RCP4.5. The built-up area will increase in the future and dominate the agricultural area under land-use change scenarios. An increase in the built-up area will at the expense of the agriculture area. Future annual abstraction in urban areas will increase under shared socioeconomic scenarios. Future annual groundwater abstraction in agricultural areas will decrease but, seasonal groundwater abstraction in the wet and dry seasons will follow the patterns of annual groundwater abstraction under climate change scenarios. Annual groundwater abstraction in urban areas will increase with the highest increase under SSP3 and lowest under SSP1. Annual groundwater abstraction in the agricultural areas will decrease. The decrease will be fast under RCP4.5 than RCP8.5 scenarios.

Future groundwater levels in urban areas will decrease fast due to an increase in abstraction and built-up expansion. Surrounding areas will observe a decline in groundwater levels due to the outward expansion of the water level depression front. Groundwater levels in urban areas will decrease by 51 m and 56 m under RCP4.5R1S1SSP1 and RCP8.5R1S1SSP1 combined scenarios, respectively, in the far future. The groundwater levels in the surroundings of urban areas will decrease by 10.7 m and 11 m under RCP4.5R1S1SSP1 and RCP8.5R1S1SSP1 combined scenarios, respectively, in the far future. Agricultural areas will experience an insignificant change in groundwater levels due to rainfall and irrigation under combined scenarios in the future. Seven adaptation options to offset the negative effect of climate change and human development on groundwater resources have been proposed based on consultation with experts. However, their actual potential has yet to be determined.