Impact of Climate Change Parameters on Groundwater Level: Implications for Two Subsidence Regions in Iran Using Geodetic Observations and Artificial Neural Networks (ANN)

Abstract

This study aims to investigate how changes in meteorological indicators affect groundwater resources, and hence to predict groundwater levels using these indicators, particularly in regions experiencing drought and subsidence. Precipitation, temperature, evapotranspiration and precipitable water vapor (PWV) are important meteorological parameters to predict groundwater levels. Two subsidence areas with different weather conditions were selected to conduct a comprehensive study on the effect of temperature and precipitation on groundwater level changes. The correct locations of the two subsidence areas were determined by analyzing Interferometric Synthetic Aperture Radar (InSAR) images of Sentinel-1A using the small baseline subset algorithm. The interferograms were processed to correct tropospheric effects using the advanced integration method. Specifying the exact locations of the two areas, the meteorological parameters were downscaled using the Statistical DownScaling Model (SDSM), synoptic observations, meteorological data, and the General Circulation Model (GCM). An Artificial Neural Network (ANN) was then employed to predict the groundwater level changes as a function of meteorological data, including Global Positioning System (GPS)-based PWV and the evapotranspiration index. The trained ANN, along with the downscaled meteorological indicators, was used to predict groundwater level changes over two time periods. In the first period, the prediction was performed over the current years to investigate the performance of the method using the available data, whereas in the second period, the prediction was performed for the coming years, up until 2030. The results confirmed the high performance of the prediction algorithm, and the importance of including PWV and evapotranspiration in groundwater level predictions. The Pearson correlation coefficient was used to check the relationship between groundwater level changes and meteorological variables. The statistical significance of these coefficients was tested at the significance level $\alpha =0.05$. In more than 80% of the cases, the correlation coefficients were statistically significant, reaching more than 0.70 in some of the months. It is also observed that an increase in the depth of groundwater level has an obvious relationship with an increase in temperature and a decrease in rainfall.

1. Introduction

Since the beginning of the Industrial Revolution, the increase in emissions of carbon dioxide, methane and other greenhouse gases has disturbed average weather conditions at global and local scales. Today, atmospheric levels of carbon dioxide have reached their highest levels in the last eight hundred thousand years [1]. Groundwater is considered an important resource, especially in arid and semi-arid regions, for reasons such as water potability, chemical composition, constant temperature, lower pollution coefficient and a higher level of reliability and sustainability in the water supply. The most important precipitation factors that affect nutrition, as a critical part of health and development, are the amount, distribution, intensity, duration and type of precipitation. For example, when the amount of precipitation exceeds the rate of evapotranspiration in a certain period, it is possible to feed the groundwater aquifer.

Many studies around the world have shown that temperatures are rising. Over the last 150 years, the average surface temperature in the Northern Hemisphere has increased by 0.76 degrees Celsius [2]. Since 1951, the average minimum and maximum temperatures in the world have almost doubled [3]. This increase in the average temperature has serious consequences for all biological and non-biological communities. Rising temperatures are also accompanied by a decrease in snow and ice cover, sea ice and mountain ice retreat, early spring in the Northern Hemisphere and an increase in the frequency of heavy precipitation events in the biological and physical indicators compatible with global warming [4]. Today, long-term changes in meteorological indicators studies are made possible by the output of GCMs and scenarios published by the Intergovernmental Panel. A review of past studies reveals the importance of these changes and their impact on water resources. Cohen examined the effect of precipitation and temperature change scenarios on precipitation and temperature and, consequently, changes in water levels in the Great Lakes of North America. The results showed that the water storage of large lakes will decrease in response to long-term changes in weather conditions [5]. Holman examined the effects of these changes and economic and social activities concerning groundwater recharge in the east of England and outlined solutions to its problems. The results showed that many factors such as changes in precipitation, changes in the thermal regime of coastal floods, urban development, forest lands and changes in cultivation and plowing can affect the future of groundwater resources [6]. Jykama and Sykes investigated the effects of weather changes on groundwater levels and of spatial changes on aquifer feeding in the basin Grand River in Ontario, Canada. The results showed that due to long-term changes in weather conditions, currents and feed volume of the aquifer increased [7]. Ng et al. reviewed the effect of these changes on groundwater recharge in the plains of northern Texas. The results showed that the amount of groundwater recharge will change between −75 to +35 percent [8]. Luoma and Okkonen studied the effects of future changes in weather conditions and rising Baltic Sea water levels on groundwater levels in southern Finland. The results showed that the groundwater recharge pattern has changed and that the seasonal effects of weather change on groundwater recharge will be in the form of floods in winter and droughts in summer [9]. Ertürk et al. evaluated the effects of current and long-term changes in weather conditions on groundwater resources in a small Mediterranean basin. According to the results, almost all components of the water balance have been reduced and water scarcity is expected to become an important issue in the future [10]. House et al. analyzed the impact of future weather changes on groundwater levels in a watershed in the United Kingdom. They concluded that due to changes in meteorological variables such as precipitation, temperature and evapotranspiration, the groundwater level decreases. In arid and semi-arid regions, groundwater is the main water resource for many applications. However, inappropriate use of groundwater reservoirs has caused the aquifers not to feed well. Therefore, it is important to study and predict the nutritional status of groundwater resources under the influence of weather change due to the fragility of these ecosystems [11]. Burgan checked the difference between ANN algorithms and Multiple Linear Regression for daily streamflow prediction [12]. Trabelsi et al. studied the effect of novel hybrid and benchmark machine learning algorithms to predict groundwater potentiality in drought areas [13].

The purpose of this paper is to investigate the relationship between meteorological parameters and groundwater depletion in subsidence regions. Iran is one of the countries that is severely facing the problem of declining groundwater levels. Iran, where groundwater is an important criterion for agricultural and urban-industrial development, is among the vulnerable countries in the field of water. A large number of shallow and deep wells have been drilled in the last few decades, resulting in the decline of the water table in aquifers. On the other hand, it has become apparent that groundwater prediction can help managers to manage groundwater resources. It is therefore important to study the causes of declining groundwater levels in different climate conditions within the country. In this research, we first use the InSAR technique to measure the subsidence in two study areas in Iran. Then, the meteorological indicators are downscaled and predicted on the subsidence center using the SDSM model and meteorological data. In the next step, the evapotranspiration and PWV parameters are computed based on the GPS observations. Unlike many studies, ignoring PWV and evapotranspiration, it is shown that these indicators have an important effect on a more accurate prediction of the groundwater level. The groundwater level is then predicted using an optimal ANN architecture. For this prediction, meteorological and groundwater level observations, along with PWV and evapotranspiration parameters, are considered as input to the training phase. The groundwater level prediction is performed over two time periods. In the first case, the prediction is performed for the current time so that the efficiency and accuracy of the method can fully be evaluated using the available observations. In the second case, the prediction will be made for the future. In both cases, the relationship between groundwater level changes and meteorological indicators will be investigated.

2. Background and Methodology

In the first step of processing, the subsidence points need to be carefully identified. Given that the weather conditions of the study area may result in a change in the subsidence signal position obtained from InSAR, a precise method should be used to eliminate the tropospheric effect in radar results and to obtain an accurate diagnosis of signals [14,15]. In the following, a detailed method to correct the tropospheric effect of InSAR measurements is described. After that, this section describes the downscaling and forecasting of climatic data as well as the neural network used to predict groundwater levels. We then show the required neural network data and how they are prepared. Finally, the results’ evaluation methods are presented.

2.1. Tropospheric Correction to InSAR Data

InSAR time series provides a powerful tool to detect surface displacement using a set of SAR acquisitions. The performance of this technique is affected by the accuracy and precision of the retrieved surface displacement and limited by the decorrelation of radar signals, phase-unwrapping error and atmospheric delay. Decorrelation is mainly caused by variations in the surface backscatter specifications over time and by the non-ideal acquisition strategy of SAR missions [16,17,18]. To overcome the problems related to the early SAR missions, including the large orbit separation between two acquisitions of the same track and the relatively long revisit time with non-regular images, two different types of InSAR time series techniques have been proposed: persistent scatterer and distributed scatterer algorithms. The persistent scatterer method uses phase-stable point scatterers. The distributed scatterer approach, which relaxes the strict limit on phase stability, includes regions that are affected by decorrelation through the exploitation of the redundant network of interferograms [19,20]. Based on the network of interferograms, the distributed scatterer methods can be divided into two categories: the first category uses the network of all possible interferograms [21] and the second category uses the network with small spatiotemporal baselines [22].

To reduce the tropospheric effect from the InSAR displacement field, both persistent scatterer and distributed scatterer algorithms traditionally rely on the spatiotemporal filtering of the phase by taking into account their different frequency specifications in the space and time domain [23]. Some of the previous studies used the empirical correlation between the tropospheric delay and topography [24]. Some new studies use atmospheric data, such as the Global Positioning System tropospheric delay and the Moderate Resolution Imaging Spectroradiometer [25,26], and some others use different approaches to estimate the tropospheric delay, such as 2D and 3D ray tracing as an advanced integration method based on meteorological data [27,28,29]. More information can be found in [30,31,32].

InSAR processes are performed using time delay and differential phase shift of radar acquisitions. When the radar signals propagate through the troposphere layers, their velocities are reduced, and the phase observations are affected due to spatial variations in tropospheric delays [33]. The tropospheric delay in the zenith direction is called the zenith tropospheric delay and is computed using the numerical integration as follows [33]:

$$\begin{array}{l}ZTD=10^{-6}\left\{\frac{k_1{R}_d}{g_m}P(z_0)+{\displaystyle \underset{z_0}{\overset{z}{\int }}((k_2-\frac{R_d}{R_v}{k}_1)\frac{e}{T}+k_3\frac{e}{T^2})dz}\right\}\hspace{0.17em}\hspace{0.17em}\\ k_1=0.776\hspace{0.17em}\mathrm{K}\mathrm{P}{\mathrm{a}}^{-}{}^1,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{k}_2=0.716\hspace{0.17em}\mathrm{K}\mathrm{P}{\mathrm{a}}^{-}{}^1,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{k}_3=3.75\times {10}^3\hspace{0.17em}{\mathrm{K}}^2\mathrm{P}{\mathrm{a}}^{-1}\end{array}$$

(1)

where z₀ is the surface elevation, z is the highest altitude of the troposphere, R_v (461.495 Jkg⁻¹K⁻¹) and R_d (287.05 Jkg⁻¹K⁻¹) are the specific gas constants for water vapor and dry air, respectively, e is the water vapor pressure, P(z₀) is the surface pressure, T is the temperature in K and g_m is the averaged gravitational acceleration over the troposphere layer. To reduce the tropospheric effect from the radar phase, it is necessary to convert the zenith tropospheric delay to the line-of-sight direction. Based on the previous studies, this step is performed using a mapping function [29]. Using mapping functions can cause errors in the calculations and results. This error can even lead to a misdiagnosis of the displacement signal in the InSAR processing [29,34]. Therefore, in this study, an advanced but simple integration technique based on direct numerical integration on the line-of-sight direction is used to prevent errors related to the mapping function. As can be seen in Figure 1, this technique, like the simple integration method along the zenith direction, requires meteorological data in 3D space. The meteorological data on the signal path need to be determined using the interpolation approaches. In this study, the kriging and spline interpolations have been used in the horizontal and vertical directions, respectively [35]. More details about this technique can be found in references [29,34].

2.2. Downscaling and Long-Term Prediction

The SDSM model is used to downscale the meteorological data at a specific point in its current status and in future periods [36,37,38]. The input data of this model are time series of meteorological variables such as precipitation, minimum and maximum temperature and other atmospheric parameters. SDSM presents three different kinds of sub-models: monthly, seasonal and annual for the downscaling. The monthly sub-model is based on a special regression equation for each month, whereas the seasonal and annual sub-models generate one regression equation for each season and each year, respectively. In this study, the monthly sub-model type is preferred because of the large monthly variations in temperature and precipitation. In this mode, the SDSM generates different statistical parameters for each month. To produce and predict the data in the future, scenario A2 of the Hadley Centre coupled model-version 3 GCM [39], as well as ERA5 reanalysis data from the European Centre for Medium-range Weather Forecasts (ECMWF) along with the synoptic observations are used. Figure 2 shows the process of downscaling and production of the scenarios using SDSM.

The process of downscaling and predicting in this model consists of six steps. (1) Data gathering and quality control, (2) selection of appropriate predictor variables, (3) model training, (4) production of current meteorological data using observable variables, (5) statistical analysis of observed data, and (6) production of future meteorological data using the predictor variables [40].

It should be noted that there are three assumptions in downscaling and prediction of climate data:

      1-. Predictor and predictand data and their relationships are not affected by human-caused climate change and can be transferred to the next decades.

      2-. All downscaled parameters are obtained using the available data in the base period because it is assumed that, during the decades belonging to this time period, the data are of a higher quality and are more accessible than in other periods and, also, that the risk of instability in the relationship between predictor variables is small.

      3-. In climate simulations, it is assumed that the data follow the normal distribution.

2.3. ANN for Groundwater Level Prediction

The ANN is a powerful and intelligent network that is known as a model for solving complex problems in various fields such as modeling, optimization, simulation, prediction and others. Figure 3 presents the steps of groundwater level prediction using ANN. The ANN model structure consists of three different layers: the input layer which receives data, the output layer which presents calculated outputs, and one or more hidden layers that work as the connection between the input and output layer [41]. A neuron is a basic processing unit of the ANN and performs two various operations: the reception of the inputs and the presentation of the output. The input data are prepared by connection weights and then passed through an activation function to present the output [42]. A multilayer feedforward ANN has neurons structured into layers and can pass in one direction without feedback connection. Multi-Layer Perceptrons (MLPs) are the most popular and the simplest state of ANN.

This class of approaches is widely used to create connections between the input and output data. MLPs are feedforward networks that include one or more hidden layers. MLP is a user-friendly method that gives a chance to determine any relationship between input and output data using the training algorithms. The Levenberg–Marquardt algorithm along with the Broyden–Fletcher–Goldfarb–Shanno, Gradient Descent and Conjugate Gradient Fletcher–Reeves Update methods have shown their strength and success in different ANN applications [38,41,42,43,44]. Among them, we selected the reliable and widely used Levenberg–Marquardt algorithm in order to perform the ANN training stage optimally. This method is a modified version of the classic Newton algorithm for computing a proper solution to the optimization processes [43,44].

2.4. PWV as a Key Indicator

One of the most important indicators that are often ignored in the groundwater level prediction process is PWV. GPS observations have been first processed using the BERNESE software [45]. The Zenith Total Delay (ZTD) and the tropospheric gradients have been computed at intervals of 1 and 4 h, respectively. The ZTD is equal to the sum of the Zenith Hydrostatic Delay (ZHD) and the Zenith Wet Delay (ZWD) [46,47]. The ZHD can be calculated from the Saastamoinen equation with acceptable accuracy [46,47]. In this study, GPS observations were used to compute PWV based on the following relationships [46,48,49]:

$$ZHD=\frac{0.0022768P}{1-0.00266\cos (2\phi )-0.00028H}$$

(2)

$$ZWD=ZTD-ZHD$$

(3)

$$PWV=\Pi \times ZWD=\frac{{10}^6}{(\frac{K_3}{T_m}+{K_2}^{’})R_V}ZWD$$

(4)

where $\Pi$ is the conversion factor, R_V, ${K_2}^{’}$ and $K_3$ are constants with values of 461.50 (JK⁻¹kg⁻¹), 16.48 (KhPa⁻¹) and 3.776 × 105 (K²hPa⁻¹), respectively, T_m (K) is the atmospheric weighted mean temperature, P is the surface pressure (hPa), and $\phi$ and H represent the latitude and the ellipsoidal height (km) of the GPS station, respectively.

2.5. Evapotranspiration Estimation

Evapotranspiration is another important parameter that helps estimate the amount of water needed for crops; it is crucial to correct irrigation planning. This parameter plays an important role in groundwater level prediction; therefore, it is considered an important input in ANN for groundwater level prediction. Because this index is not directly available to the users, it is necessary to determine its values by the experimental relationships. The Food and Agricultural Organization-proposed Penman–Monteith (PM) equation is the most popular equation to estimate evapotranspiration [50]:

$${\mathrm{ET}}_{(PM)}=\frac{0.408\left(Rn-G\right)+\gamma \frac{900}{T_a+273}{U}_2\left(e_s-e_a\right)}{\Delta +\gamma \left(1+0.34U_2\right)}$$

(5)

where ET is evapotranspiration, G is soil heat flux density, Rn is net radiation at the crop surface, Ta is the mean daily air temperature at 2 m height, e_s is saturation vapor pressure, U₂ is the wind speed at 2 m height, e_a is actual vapor pressure, ∆ is slope vapor pressure curve and γ is psychrometric constant. It should be mentioned that this parameter is not equal to groundwater evapotranspiration. Due to the dependence of this relationship on a large number of meta-indicators that are often not published by databases or are accessible with different temporal and spatial resolutions, the use of this relationship is usually difficult. Unlike the PM equation, the TH method can estimate the amount of evapotranspiration based on air temperature [51].

$${\mathrm{ET}}_{(TH)}=\{\begin{array}{l}0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}T\le 0\\ 16K{(\frac{10T}{l})}^m\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}T>0\end{array}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}$$

(6)

where l and m are the heat factor and coefficient, respectively, and K is the calibrated coefficient of latitude and month. Due to the sensitivity and importance of this index and its role in predicting the groundwater level, it was necessary to calculate this index with the utmost possible accuracy. For this reason, a new method based on GPS observations has been used in this study [52]. In this method, the evaporation index value based on the PM method is considered as an accurate model. The aim of this algorithm is to model the difference between the PM and TH methods over time by using water vapor obtained from GPS, as well as temperature. Finally, the modeled difference should be added to the TH formula to correct it. DET can be modeled using a linear equation based on the temperature and water vapor:

$$DET=\hspace{0.17em}E{T}_{PM}-E{T}_{TH}$$

(7)

$$DET=\{\begin{array}{l}{a}_0\times PW{V}_1\hspace{0.17em}\hspace{0.17em}+\hspace{0.17em}\hspace{0.17em}{a}_1\times T\hspace{0.17em}\hspace{0.17em}+\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{a}_2\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}T\ge 0\\ b_0\times PW{V}_2\hspace{0.17em}\hspace{0.17em}+\hspace{0.17em}\hspace{0.17em}{b}_1\times T\hspace{0.17em}\hspace{0.17em}+\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{b}_2\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}T<0\end{array}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}$$

(8)

where a₀–a₂ and b₀–b₂ are the model coefficients, which can be determined based on the least squares method. PWV₁ and PWV₂ are the PWV at the moment of $T\ge 0$ °C and $T<0$ °C, respectively. Finally, the more accurate value of the evapotranspiration index can be calculated using the following equation:

$$E{T}_{Accurate}=DET+E{T}_{TH}$$

(9)

Therefore, the ET_Accurate can be computed at different times based on the differential model and without dependence on many meteorological indices. If there is a need to calculate this index as a two-dimensional map, it is also possible to enter the geographic location into the fitted model. Considering that it is necessary to calculate this index at the location of the station, the differential model is considered only in time. More information about how to use the differential model in temporal and spatial dimensions can be found in [52].

2.6. Effective Rainfall

Agriculture is one of the most important reasons for groundwater consumption and hence the reduction in its level. This variable should be considered in groundwater level prediction. In other words, to investigate the impact of meteorological parameters on the groundwater level, it is necessary to calculate the amount of water consumed through agriculture and consider it in the prediction process. The amount of consumed water is achievable according to the statistics provided by the water organization of each country. Based on the relationship between the amount of water consumed in agriculture and the effective rainfall (ER) index, it is possible to predict the amount of consumed water in the future. In this contribution, this value is removed from the prediction process. In this case, the final changes in the groundwater level can only be considered as a result of changes in meteorological parameters. The ER is estimated by different methods, such as the ratio method. The approximation is based on the properties of soil and runoff, and dryness also plays a role [53]:

$$ER=\frac{ET}{P}$$

(10)

where ET is evapotranspiration and P is rainfall. In this paper, the information provided by the Iranian Meteorological Organization (IMO) and the Iranian Geological Organization (IGO) has been used to correct the effect of consumed water.

2.7. Validation Methods

The statistical indicators, namely, the Root Mean Square Error (RMSE), Nash–Sutcliffe efficiency (NSE) and the coefficient of determination (R²) have been used to compare and validate the results [54]. The Pearson correlation coefficient was used to determine the correlation between water level changes and meteorological indicators. In addition, the significance of the estimated correlation coefficient is determined based on the given confidence interval. The statistical tests have a null hypothesis. For most tests, the null hypothesis is that there is no relationship between your variables of interest or that there is no difference between groups.

Null hypothesis (H₀): there is no significant difference between the two sets of data

Alternative hypothesis (H₁): there is a significant difference between the two sets of data

The p-Value is a number, calculated from a statistical test, that describes how likely you are to have found a particular set of observations if the null hypothesis were true. p-Value is used in hypothesis testing to help decide whether to reject the null hypothesis. The smaller the p-Value is, the more likely you are to reject the null hypothesis. p-Value is most often used by researchers to say whether a certain model they have presented is statistically significant. Statistical significance is another way of saying that the p-Value of a statistical test is small enough to reject the null hypothesis of the test [54]. The threshold considered in this study is p-Value < 0.05. In other words, this test is performed at a significance level of α = 0.05.

3. Study Area and Data Set

Iran is a country with different weather conditions in different regions, and, in recent years, it has faced severe subsidence due to declining groundwater levels. Two regions with different weather conditions are selected in Iran. The first study area is part of the Qazvin plain with mountainous weather conditions and the second area is part of the Yazd plain with hot and dry weather. The location of the study areas can be seen in Figure 4. According to the statistics provided by the United Nations, Iran’s population is currently 86,022,843, which will increase by 3% in 2030 [55]. According to the information provided by the Statistical Center of Iran, the study areas of this research have the lowest population growth rate in Iran. These two areas are also outside urban areas, and increased personal consumption due to population increase has the least impact on the groundwater level. However, by employing the existing algorithm, the use of groundwater due to population changes has also been removed from the observations [56]. This includes, for example, the effects related to agricultural uses. Therefore, the results can be attributed to meteorological indicators.

The average precipitation and temperature in these two regions are different and this difference provides the basis for a detailed investigation. Figure 5 and

Table 1. Monthly average of precipitation and temperature from 1979 to 2021.

Parameter	Area	Jan	Feb	Mar	Apr	may	Jun	Jul	Aug	Sep	Oct	Nov	Dec
Average precipitation (cm)	Qazvin	0.50	0.62	0.78	0.75	0.41	0.09	0.05	0.05	0.05	0.33	0.56	0.49
	Yazd	0.28	0.26	0.37	0.30	0.16	0.01	0.01	0.00	0.00	0.05	0.17	0.22
Average temperature (k)	Qazvin	268.8	271.9	276.9	277.0	281.8	285.9	287.2	285.9	281.0	274.2	268.1	263.8
	Yazd	282.3	285.7	289.8	289.1	294.4	298.3	298.7	296.8	292.3	285.7	280.3	276.6

show the difference between precipitation and temperature in the center of these two regions. These two meteorological parameters show visible variations in the two study areas.

In this study, three types of meteorological data are used. The first type is the daily meteorological data from synoptic stations that are used for the training and evaluation step. The second type of meteorological information is the ERA5 reanalysis data from ECMWF that have been used as large-scale data for training. The tropospheric correction of the InSAR observations is also performed using the ERA5 reanalysis data from ECMWF. It presents values of the meteorological information on 37 pressure levels, and with a spatial resolution of about 31 km from 1950 to the present [57]. Based on the previous investigations, the reanalysis data are useful and practical in different applications such as geodesy, meteorology, remote sensing and geodynamics. The third type is the Hadley Centre coupled model-version 3 GCM information based on different scenarios that have been used for meteorological data prediction. This model presents different parameters such as temperature, humidity and wind speed. This model has 19 levels with a horizontal resolution of 3.75 degrees longitude and 2.5 degrees latitude. This is equivalent to the Earth’s surface resolution of about 417 km × 278 km at the Equator, reducing to 295 km × 278 km at 45 degrees latitude [39].

Based on the aim of this study, the InSAR technique has been used to identify the exact location of the subsidence signal. To implement the InSAR technique, a number of Sentinel-1A radar acquisitions taken from these two areas have been used. The radar images have been provided by the European Space Agency (ESA). The specifications of radar images can be seen in

Table 2. Specifications of the radar acquisitions using Sentinel-1A.

Area	Mission	Product	Track	Flight Direction	Beam Mode	Time Span
Qazvin	Sentinel-1A	S1A-IW-SLC	108	Descending	IW	2015–2020
Yazd	Sentinel-1A	S1A-IW-SLC	130	Ascending	IW	2015–2020

.

Groundwater level observations from the existing wells have also been used to calibrate and validate the ANN and, therefore, to predict groundwater level changes. The groundwater level data and needed information such as groundwater consumption have been prepared by Iranian regional water companies, IMO and IGO. Some of the statistical properties of groundwater level data can be seen in

Table 3. Statistical properties of groundwater level depth data.

Year	Area	Min (m)	Max (m)	Mean (m)	Median (m)
2013 2014 2015 2016 2017 2018 2019 2020	Qazvin (2 wells)	0.49 0.51 0.58 0.47 0.59 0.57 0.60 0.62	1.86 1.98 2.23 2.39 2.64 2.63 2.71 2.81	1.31 1.29 1.13 1.39 1.28 1.41 1.53 1.61	1.24 1.15 1.11 1.20 1.39 1.47 1.43 1.50
2013 2014 2015 2016 2017 2018 2019 2020	Yazd (2 wells)	0.50 0.48 0.59 0.63 0.58 0.61 0.67 0.71	2.54 2.64 2.87 2.86 3.12 3.25 3.39 3.31	1.42 1.49 1.69 1.87 1.88 2.01 2.25 2.19	1.31 1.28 1.54 1.63 2.01 2.15 2.11 2.40

. GPS measurements obtained from the existing stations in both study areas have been used to estimate PWV and ET.

4. Processing Results and Discussions

Figure 6 presents the processing steps to study the relationship between changes in meteorological indicators and groundwater levels.

In the first step, it is necessary to determine the subsidence area using InSAR. After determining the subsidence zone, SDSM is trained using large-scale and local-scale meteorological data. Then, using the local meteorological observations, the trained model is validated, and finally, the values of the temperature and precipitation are downscaled to the location of the desired stations. The meteorological parameters, in the future, will then be predicted using SDSM. In the next step, it is necessary to select the optimal ANN parameters to predict the groundwater level in the future. The meteorological data, along with groundwater observations, are then considered as input to the training phase. The trained ANN is evaluated using local observations, and the optimal ANN specifications are selected. The groundwater level is then predicted. Finally, the relation between the temperature and precipitation changes and the groundwater level can be studied.

4.1. Subsidence Detection

A number of Sentinel-1A radar acquisitions have been used to perform the InSAR analysis. In this process, the Amplitude Dispersion Index has been selected to be 0.42. First, the interferograms were selected according to their spatial coherence (Figure 7). In addition, a reference point has been selected with respect to the smallest phase variance and the maximum of coherent pixels. In order to compute the displacement velocity, disturbing effects must be removed from the obtained interferograms. After removing the disturbing parameters using the Digital Elevation Model and orbital data, the interferograms only contain tropospheric effects and the displacement phase. The tropospheric effect was estimated using the advanced integration method and applied to the interferograms. An example of the obtained tropospheric delay can be seen in Figure 8. Finally, the displacement velocity map was calculated using the small baseline subset method. Figure 9 shows the obtained displacement velocity for the two study areas.

The computed displacement velocity in the Yazd region indicates a subsidence signal at the center of this area. The amount of the subsidence reaches 150 mm/year, which is considered to be a high and dangerous rate. There is no significant displacement in other parts of the area. In the Qazvin area, the subsidence signal is clearly visible. The maximum subsidence rate is 100 mm/year. According to the location of the subsidence signals in both areas, the local groundwater observations, as well as the synoptic station, have been selected.

4.2. Downscaling and Prediction of Meteorological Parameters Using SDSM

The processing of this paper takes place over two periods of time. In the first period, the purpose is to predict the groundwater level between 2013 and 2020 so that the results and methods can be fully evaluated using the available observations. The second mode is making predictions for the future, between 2021 and 2030. Therefore, the prediction of meteorological indicators has also been conducted in two time frames. In the first case, the data between 1985 and 2010 have been used for training and the data between 2011 and 2012 were used to evaluate the model. In addition, the data between 1985 and 2018 have been used for training and the data from 2019 to 2020 were used to evaluate the predictions made for meteorological indicators. The comparison results between the downscaled mean monthly temperature and precipitation and their observed values, in terms of the above-mentioned statistical measures, are provided in

Table 4. Statistical measures of observed and downscaled mean monthly temperature during validation period from 2011 to 2012.

Area	Variable (°C)	RMSE (°C)	NSE	R2
Qazvin	Tmax	0.41	0.77	0.85
	Tmin	0.31	0.83	0.87
Yazd	Tmax	0.36	0.79	0.85
	Tmin	0.43	0.84	0.89

and

Table 5. Statistical measures of observed and downscaled mean monthly precipitation during validation period from 2011 to 2012.

Area	Variable (mm)	RMSE (mm)	NSE	R2
Qazvin	Precipitation	2.21	0.75	0.84
Yazd	Precipitation	2.02	0.71	0.82

.

Finally, using the downscaled data as well as the produced scenario, the values of meteorological parameters in the time span of 2013 to 2020 have been predicted. Figure 10 and Figure 11 show the monthly average of predicted meteorological variables in the time span of 2013 to 2020 compared to the base period. According to the prediction, the temperature in both regions will increase in the coming years. This increase is more noticeable in the warmer months of the year. In the Yazd region, the maximum temperature changes, compared to the base period, are about 3 degrees for the years 2013 to 2020. In the Qazvin region, the maximum rate of temperature increase, compared to the base period, is about 2.5 degrees. The results of precipitation prediction show that this parameter will decrease in the coming years and that this decrease is more significant in the cold months of the year. The biggest decreases for the Qazvin region in this period will occur in January, September and November by about 9 mm. After predicting the meteorological parameters and computing the ET and PWV, we can now predict the groundwater level changes using the ANN, which we will show in the next section.

4.3. Groundwater Level Prediction from 2013 to 2020

In order to check the effect of the use of PWV and ET and the statistical evaluation of all the mentioned processes, the prediction period is considered between 2013 and 2020. In this case, using the available observations, a comprehensive statistical evaluation of the method can be provided. In the first step, it is necessary to select an optimal ANN setting. The optimal ANN is organized to predict groundwater level changes. Based on the previous studies, the Levenberg–Marquardt algorithm has been selected as the most appropriate training function. The sigmoid log and hyperbolic tangent activation functions have been used for the hidden layers, and the linear activation function has been selected for the output layer. To determine the impact of PWV and ET on the prediction process, various input sets have been used to train the ANN.

Table 6. Different input sets of the ANN. D refers to groundwater level depth.

Set Number	Inputs
1	Precipitation, Tmin, Tmax, SR, D
2	Precipitation, Tmin, Tmax, SR, D, PWV
3	Precipitation, Tmin, Tmax, SR, D, ET
4	Precipitation, Tmin, Tmax, SR, D, ET, PWV

shows different types of input parameter sets for the ANN to predict groundwater levels. These parameters generally include precipitation, Tmin, Tmax, groundwater level depth, PWV, ET and Solar Radiation (SR). Based on the algorithm, the ANN has been launched using different sets of inputs. In this process, 80% of the data are used for training and the remaining 20% are considered for validation. Figure 12 shows the comparison between the observations and the simulated changes in the groundwater level from different input sets. RMSE of different modes obtained from the validation step can be seen in

Table 7. The statistical results of different input sets of ANN.

Area	Set Number	RMSE (m)	NSE	Correlation
Qazvin	1	0.989	0.71	0.79
Qazvin	2	0.615	0.75	0.82
Qazvin	3	0.361	0.81	0.89
Qazvin	4	0.189	0.86	0.93
Yazd	1	1.023	0.70	0.74
Yazd	2	0.084	0.73	0.81
Yazd	3	0.598	0.80	0.87
Yazd	4	0.241	0.84	0.91

.

The RMSEs show that data set No. 4 is the optimal ANN to predict groundwater level changes. These comparisons show that the use of PWV and ET obtained from GPS measurements leads to an increase in prediction accuracy. As it can be seen, when PWV and ET are not used (set 1), in addition to creating large errors, it also leads to a change in the trend of the graph in some parts. The biggest difference between the results obtained from data set 1 and the observations can be seen in the warmer months of the year. This is because, in the warmer months, PWV and ET play more important roles in the calculations than in the colder months. The results show that ANN with input set number 4 has a high ability to simulate the groundwater level in both areas and provides a reliable tool in predicting these variations in the coming years.

After checking the accuracy of the algorithm, in this step, the relationship between the time variations in temperature and precipitation and groundwater level changes is investigated. For this aim, the heatmap plot of the monthly average of groundwater level depth changes according to the baseline period is visible in Figure 13. Based on Figure 13 it can be concluded that most of the variations occurred in the warmer months of the year. After 2015, the process of increasing the groundwater level depth compared to the base period has become faster. This amount reached more than 50 cm in July and June 2020. This trend is evident in both study areas and it can be considered a result of a decrease in precipitation and an increase in temperature. The lowest number of changes are related to the beginning and end months of the year, which can be attributed to fewer changes in meteorological indicators compared to the warmer months of the year. The statistical results related to the comparison between the monthly averages of predicted groundwater level changes and temperature and precipitation changes can be seen in

Table 8. Correlation coefficient between changes in groundwater level and meteorological parameters. Green color refers to p-Value smaller than 0.05 and red color shows p-Value greater than 0.05.

Month		Jan	Feb	Mar	Apr	May	Jun	Jul	Aug	Sep	Oct	Nov	Dec
Precipitation	Yazd Well 1	0.39	0.56	0.49	0.54	0.58	0.61	0.68	0.51	0.51	0.68	0.73	0.62
	Yazd Well 2	0.40	0.35	0.37	0.48	0.55	0.63	0.69	0.59	0.55	0.63	0.69	0.66
	Qazvin Well 1	0.39	0.46	0.51	0.59	0.50	0.48	0.62	0.51	0.68	0.63	0.75	0.72
	Qazvin Well 2	0.43	0.40	0.49	0.51	0.57	0.51	0.60	0.63	0.71	0.61	0.68	0.66
Temperature	Yazd Well 1	−0.41	−0.37	−0.47	−0.55	−0.59	−0.69	−0.73	−0.75	−0.62	−0.54	−0.48	−0.44
	Yazd Well 2	−0.41	−0.42	−0.44	−0.57	−0.65	−0.63	−0.70	−0.67	−0.59	−0.58	−0.68	−0.69
	Qazvin Well 1	−0.36	−0.43	−0.59	−0.55	−0.63	−0.66	−0.62	−0.76	−0.71	−0.66	−0.59	−0.51
	Qazvin Well 2	−0.40	−0.51	−0.46	−0.57	−0.63	−0.59	−0.69	−0.74	−0.71	−0.68	−0.61	−0.59

.

In Table 6, two items are considered: first, the Pearson correlation coefficient, and second, the significance of the computed correlation using p-Value. It should be noted that this test has been conducted at a significance level of 0.05. As can be seen, the correlation coefficient between water level depth changes and meteorological indicators reaches 76%. As expected, this correlation is positive for the precipitation index and negative for the temperature index. In fact, with the increase in temperature, the groundwater level depth decreases and, with the increase in precipitation, the groundwater level depth increases. The p-Value test shows that the correlation coefficients were significant in 83% of cases. It takes at least one month to establish a meaningful relationship between indicators. For this reason, at the beginning of the prediction process, the correlations obtained from the p-Value test were not accepted.

4.4. Groundwater Level Prediction from 2021 to 2030

As mentioned in the previous sections, in this study, we tried to take into account all indicators that affect water level changes. They include meteorological changes, agricultural uses and population increases. It was observed that the algorithm has the ability to model water level changes with appropriate accuracy. At this stage, we will predict the changes in groundwater level in the coming years, between 2021 and 2030. The processing process is exactly the same as before. The only difference is the time frame of the prediction of meteorological indicators as stated in Section 4.2. In this case, more accuracy is expected from the final results because the duration of the training period for meteorological data prediction has increased. The groundwater level changes over the period from 2021 to 2030 have been predicted. Figure 14 represents the heatmap plot of the variations in groundwater level depth from 2021 to 2030 compared to the baseline period. The results obtained for the correlation coefficient and its significance test can be seen in

Table 9. Correlation coefficient between changes in groundwater level and meteorological parameters. Green color refers to p-Value smaller than 0.05 and red shows p-Value greater than 0.05.

Month		Jan	Feb	Mar	Apr	May	Jun	Jul	Aug	Sep	Oct	Nov	Dec
Precipitation	Yazd Well 1	0.31	0.39	0.41	0.61	0.50	0.63	0.69	0.51	0.43	0.65	0.62	0.68
	Yazd Well 2	0.31	0.35	0.37	0.47	0.51	0.55	0.50	0.61	0.52	0.58	0.63	0.68
	Qazvin Well 1	0.31	0.33	0.43	0.59	0.68	0.69	0.73	0.76	0.72	0.62	0.69	0.61
	Qazvin Well 2	0.31	0.40	0.51	0.59	0.58	0.63	0.66	0.67	0.71	0.70	0.62	0.53
Temperature	Yazd Well 1	−0.37	−0.42	−0.53	−0.57	−0.61	−0.52	−0.69	−0.61	−0.59	−0.55	−0.51	−0.61
	Yazd Well 2	−0.40	−0.33	−0.61	−0.55	−0.51	−0.61	−0.68	−0.53	−0.48	−0.54	−0.49	−0.57
	Qazvin Well 1	−0.37	−0.39	−0.44	−0.51	−0.59	−0.52	−0.46	−0.41	−0.44	−0.58	−0.51	−0.58
	Qazvin Well 2	−0.44	−0.45	−0.51	−0.58	−0.64	−0.67	−0.63	0.51	−0.69	−0.66	0.53	−0.49

.

Like the previous processing, Figure 14 shows that the most changes occurred in the warmer months of the year. These changes will reach more than 75 cm in 2030. The groundwater level changes for the Yazd region are more than Qazvin. This issue can be attributed to the rapid decrease in precipitation and the increase in temperature in Yazd compared to Qazvin. The most changes have happened in June and July, which is due to more significant changes in temperature and precipitation in these months.

The calculated coefficients in Table 7 show that the correlation between groundwater level changes and meteorological indicators reaches 76%. In 86% of the cases, the relationship between water level changes and meteorological indicators is significant at the significance level of 0.05. As mentioned before, it takes at least one month to establish a meaningful relationship between the variables. For this reason, the p-Value test was rejected in the first months.

5. Summary and Conclusions

This paper focused on two issues: first, providing a specific and reliable algorithm to predict groundwater level changes, and second, investigating the correlation between these changes and meteorological indicator variations. Two regions in Iran with different weather conditions were considered. Using Sentinel-1A acquisitions, the Interferometric Synthetic Aperture Radar (InSAR) technique and advanced integration methods, the exact location of the subsidence in both areas was determined. The meteorological parameters were downscaled and predicted in these areas using the Statistical DownScaling Model (SDSM), meteorological data, synoptic observations and the General Circulation Model (GCM). The groundwater level prediction was conducted using meteorological data, groundwater level observations, Perceptible Water Vapor (PWV) from Global Positioning System (GPS) measurements and the evapotranspiration index based on the Artificial Neural Network (ANN). Predictions were conducted in two different modes. In the first mode, the prediction was made for the current years to validate the proposed algorithm using observations. In the second mode, the prediction was conducted for future time spans. Comparing the prediction with the observations showed that water vapor and evaporation are very important in predicting surface changes. In addition, the mentioned algorithm has a strong ability to model groundwater level changes due to variations in meteorological parameters if additional effects are removed. After making predictions, the relationship between groundwater level changes and meteorological indicators was studied using the Pearson correlation coefficient. The significance test was performed at the significance level of 0.05 using the p-Value test. The results showed that there is a strong relationship between groundwater level and meteorological indicators. The value of the multiple correlation coefficient between groundwater level variations and meteorological data reaches more than 70%. In more than 80% of cases, the relationship between groundwater level changes and meteorological parameters is significant.