**2 − 1. Study area and hydrogeological characteristics**

The Nahavand Plain lies within the upper watershed of the Karkheh River, which is fed by tributaries in the Zagros Mountains. This mountain range trends NW-SE from northwestern Iran to the Strait of Hormuz, spanning the length of the western and southwestern Iranian plateau and paralleling the coastline of the Persian Gulf. Elevation of the upper Kharkeh watershed elevations varies from 3400 meters above mean sea level (msl) in the southeast to less than 1400 above mean sea level (msl) in the northwest. The area of the Nahavand Plain aquifer is 513 km2 (Anonymous, 2015).

The mountain named as, Kuh-e Garin with peaks over 3600 m Elevation geologically is limestone (Miocene) and is structurally overlain to the north by a north-dipping thin thrust slice of Oligocene–Miocene sandstone, marl and limestone. SE of Dorud the Main Recent Fault is developed along steepened thrusts of the High Zagros Belt (Mohajjel and Fergusson, 2016).

The region falls within the Csa climate classification (hot-summer Mediterranean) in the Köppen-Geiger system (Climate-Data.org, 2022). Annual rainfall is estimated at 435.1 mm in the upper ranges and 366.5 mm in the plains. The estimated amount of free-surface evaporation is 1047 to 1211 mm per year.

Much of the upper Kharkeh watershed is underlain by carbonate rocks (primarily limestone), although volcanic rocks (andesitic tuff, other andesitic and basaltic rocks), metamorphic rocks (metavolcanics, schist, phyllite, slate, marble, and quartzite), and fluvial and piedmont conglomerates also occur (Fig. 1). The Nahavand Plain is marked by piedmont fan and valley terrace deposits. The alluvium reaches a maximum thickness of 70 to 80 m in the central parts of the plain. Carbonate rocks in the area attain a thickness of several kilometers and contain interbeds of shale, shaly limestone, marl and marly limestone, and dolomite. These rocks range from thickly bedded to massive and extremely fine grained to granular. High relief and active tectonics have promoted fracturing and development of karst features (Ghobadi et al., 2012), including enlarged joints with rounded edges, wavy ridges and pitted surfaces (Taheri Tizro, 2002). Infiltration through karstified carbonate rocks around the margins of the Nahavand Plain recharges the alluvial aquifer. It should be noted that karst formations and therefore aquifers are recharged by the infiltration of rainfall and snowfall on these formations. The percentage of rainfall that infiltrates through the carbonate rocks varies between 35 and 55% of the annual atmospheric precipitation. Transmissivity ranges from less than 100 m2/day along the margins of the plain to 100 to 300 m2/day for the alluvial aquifer.

**2–2. Selected model to study groundwater level changes in area**

2-2-1- Self-correlated model - integrated moving average ARIMA (p, d, q)

In a real-time series of hydrology and hydrogeology, most variables are unstable, and therefore the AR, MA, and ARMA models cannot be used in unstable processes. It has been observed that sometimes the difference between successive values of an unstable time series forms a static time series. Therefore, a second-order operator can use a differential operator to convert an unstable process to a static process. In short, it can be said that if Wt = ∇ d xt is static, where ∇d represents the Xt process difference, Wt is described by a static ARMA model of order (p, q). ARIMA model with order (p, d, q) can be called an integrated autoregressive model and moving average (Shahin et al., 1993)

The ARMA stochastic model is defined as Eq. 1:

$${\text{X}}_{\text{t}}= {\text{W}}_{\text{t}}+ {\text{B}\text{W}}_{\text{t}}+ {\text{B}}^{2}{\text{W}}_{\text{t}}+\dots$$

1

The ARMA model can be generalized from different directions. The models obtained from the ARMA model, if p = 0, the autoregressive model or (AR), if q = 0, the moving average model or (MA), if d = 0, the ARMA model. The ARIMA model is called an unstable model because it considers the instability of the data in the model. If the model is not seasonal and d = 0, then we have the static model ARMA (p, q).

2-2-2. Self-correlated model - average seasonally combined mobility SARIMA (p, d, q) (P, D, Q)s

If the remainder is unstable (seasonal trend), the initial difference can be used to eliminate the instability. This model is then defined as an integrated moving and autoregressive average (ARIMA) model, or seasonal ARIMA (SARIMA) with the difference of time series with the first series of the first-order operator itself as follows:

$${\text{r}}_{\text{t}}-{\text{r}}_{\text{t}-1}=(1-\text{B}){\text{r}}_{\text{t}}$$

2

\({\text{r}}_{\text{t}}\) The rest of the model after removing the seasonal dependence. B The backward operator is \(\text{B}{\text{r}}_{\text{t}}= {\text{r}}_{\text{t}-1}\)

To eliminate the seasonal nature of the data, we subtract the time series values from our previous values. For seasonal data, a seasonal difference operator is preferably used.

$${\text{r}}_{\text{t}}-{\text{r}}_{\text{t}-\text{s}}=(1-{\text{B}}^{\text{s}}){\text{r}}_{\text{t}}$$

3

s indicates the length of the seasonal period. For monthly data is s = 12 (Murthy et al., 2019) For non-negative integer d, Ds of rt series, the model is SARIMA (p,d,q) (P,Ds,Q) with period S.

If the differential series \({\text{Z}}_{\text{t}}={(1-\text{B})}^{\text{d}}({1-{\text{B}}^{\text{s}})}^{{\text{D}}_{\text{s}}} {\text{r}}_{\text{t}}\) is ARMA model, it is defined as follows

$${{\phi }}_{\text{p} }\left(\text{B}\right){{\Phi }}_{\text{p}}\left({\text{B}}^{\text{s}}\right){\text{r}}_{\text{h}}={{\theta }}_{\text{q}}\left(\text{B}\right){{\Theta }}_{\text{Q}}{\text{B}}^{\text{s}}{\text{e}}_{\text{t}}$$

4

For non-seasonal data:

$${{\phi }}_{\text{p} }\left(\text{B}\right)=\left(1-{{\phi }}_{1}\text{B}-{{\phi }}_{2}{\text{B}}^{2}-\dots -{{\phi }}_{\text{p}}{\text{B}}^{\text{p}}\right) , {{\theta }}_{\text{q}}\left(\text{B}\right)=(1+{{\theta }}_{1}\text{B}+{{\theta }}_{2}{\text{B}}^{2}+\dots +{{\theta }}_{\text{q}}{\text{B}}^{\text{q}})$$

5

For seasonal polynomials in B:

$${{\Phi }}_{\text{P}}\left(\text{B}\right)=\left(1-{{\Phi }}_{1}{\text{B}}^{\text{s}}-{{\Phi }}_{2}{\text{B}}^{2\text{s}}-\dots -{{\Phi }}_{\text{P}}{\text{B}}^{\text{P}\text{S}}\right) , {{\Theta }}_{\text{Q}}\left(\text{B}\right)=(1+{{\Theta }}_{1}{\text{B}}^{\text{s}}+{{\Theta }}_{2}{\text{B}}^{2\text{s}}+\dots +{{\Theta }}_{\text{Q}}{\text{B}}^{\text{Q}\text{s}})$$

6

p and q are the order of AR and MA, respectively. P and Q are the seasonal orders of AR and MA models, respectively. d is the number of non-seasonal differences and Ds is the number of seasonal differences (Fang and Lahdelma, 2016).

2–3. **Modeling and mapping of water-level changes**

Steps for modeling time series include model recognition, model fit, and model validation (Box and Jenkins, 1976). Time-series data in this study included water levels from 29 observation wells of the alluvial aquifer, recorded monthly from April 1997 through March 2017. To test for trend (mean instability), water-level data were linearly fitted using Minitab 17 software. Different models were fitted to the time series for each well and the best model was selected based on the lowest value of the Akaike information criterion (AIC) (Sharma et al., 2019) as well as root mean square error (RMSE), The MNSE range between 0 and 1 indicates an acceptable level of performance, while MNSE equal to 1 is considered the best and ∞ is the worst (Machiwal and Jha, 2015) and R2 values. To validate the models, simulated water-level data were plotted against observed water-level data for the period April 2012 to March 2014.

For each well, the selected model was fitted to the entire data series (April 1997–March 2017) and the water level was predicted for the next 5 years (April 2018–March 2022). Water-level changes were mapped annually for 2012 to 2019 and for the period 2018–2022. Maps were based on values in November because groundwater is pumped for agriculture in summer, when there is a lack of rainfall but not in late autumn. Ordinary kriging and inverse-distance weighting were used to interpolate modeled water-level data for mapping. Water levels predicted by the time-series models were divided into five intervals and the area for each of the defined intervals was calculated using ArcGIS10.4 software.