Understanding the effects of subsidence on unconfined aquifer parameters by integration of Lattice Boltzmann Method (LBM) and Genetic Algorithm (GA)

Excessive exploitation of groundwater has hitherto led to a significant land subsidence in a considerable number of plains in Iran. The compaction of aquifer layers ends up with changes in aquifer parameters, including hydraulic conductivity (Kx), specific yield (Sy), and compressibility (α). Accordingly, a precise estimation of aquifer parameters, Kx, Sy, and α seems essential for future water resources planning and management. In this study, an innovative inversion solution based on the combination of lattice Boltzmann method (LBM) and genetic algorithm (GA) was developed to determine the aquifer parameters, Kx, Sy, and α in Darab plain (in Fars province, Iran), which is highly subject to land subsidence. Herein, a newly developed LBM for unconfined groundwater flow was employed by incorporating the amount of subsidence measured by synthetic aperture radar interferometry (InSAR) spanning from 2010 to 2016. In order to optimize the aquifer parameters, the whole process of inverse modeling is replicated on the annual basis from 2010 to 2016 which leads to the temporal estimation of the aquifer parameters. Due to the compaction occurring in the aquifer system, a declining temporal trend is observed in the aquifer parameters in most parts of the plain. By fitting a function to time-dependent aquifer parameters, Kx, Sy, and α, their corresponding values and consequently the amount of subsidence in the near future, i.e., 2017, are predicted. The small average relative error (~ 3.5%) between the predicted land subsidence and the InSAR measurements demonstrates the high performance of the proposed inverse modeling approach.


Introduction
Various methods have been developed as essential tools to study water resources in recent decades. Some of these studies have addressed the effects of climate change (Ostad Ali Askari et al. 2019), whereas some others have dealt with the modeling of groundwater flow (Barati-Moghaddam et al. 2021;Hekmatzadeh et al. 2020;Ostad Ali Askari et al. 2017).
Modeling the fluid flow through porous media is of great importance quintessentially due to its applications in groundwater extraction, groundwater contamination and storage of nuclear waste (Fattahi et al. 2016;Cheng et al. 2019). The groundwater flow equation consists of partial differential equations as a function of aquifer parameters such as hydraulic conductivity and storage coefficient. An accurate knowledge of the aquifer parameters is of crucial importance for groundwater flow prediction. Determining the aquifer parameters based on hydraulic head as observations is known as 'inverse solution of the groundwater flow' (Garcia and Shigidi 2007).
The inversion methods based on the optimization tools have received special attention in recent years (Chan and Kaw 2020; Ghate 2020). Numerous optimization methods have been developed to improve the identification of groundwater parameters. Among them, the genetic algorithm (GA) has been extensively jointed with numerical solutions of differential equations (Xia et al. 2019) due to their specific advantages such as easy convergence to the global optimal solution as well as their high accuracy (Ayaz 2017). For instance, GA have been coupled with the advection-dispersion equation (ADE) (Han et al. 2020), groundwater equation (Aral et al. 2001), and MODFLOW and MT3DMS models (Sophia and Bhattacharjya 2020). Mahinthakumar and Sayeed (2005) used the real encoded genetic algorithm (RGA) to understand the location and release intensity of the two-dimensional pollution area. Significant endeavors have been made to solve the inverse method of groundwater equations and calibration process using GA, for aquifer parameters estimation (Sonnenborg et al. 2003;Geza et al. 2009;Abdelaziz and Merkel 2015;Lu et al. 2014). Gentry et al. (2003) employed the GA to solve the inverse problem of water movement from a shallow aquifer to a semi-confined aquifer. Finite element method (FEM) and finite difference method (FDM) were extensively employed in the inversion groundwater and pollution transport problems. For example, Seyedpour et al. (2021) used the GA coupled with FDM to optimize remediation process in groundwater flow. In this respect, more studies can be found in Kavousi et al. (2020), Li et al. (2019), Lorza et al. (2018), and Albu et al. (2020). Employing FEM or FDM joined with optimization tools for the inversion problems may be time dependent and consequently inefficient regarding excessive computational time (Hekmatzadeh et al. 2020).
Instability is another issue of time dependent FEM and FDM, leading to small time steps (Hekmatzadeh et al. 2020). In the last two decades, the lattice Boltzmann method (LBM) has emerged as a robust numerical tool to solve time-dependent partial differential equations. The LBM are known as fast computational methods with admirable stability in time dependent problems (Hekmatzadeh et al. 2019;Mishra et al. 2005;Servan-Camas and Tsai 2009). Other important advantages of LBM include but not limited to the ability to model complex boundary conditions, easy programing, and superior capability of parallel programming utilizing graphics processing unit (Safdari-Shadloo 2019; Zhang et al. 2019).
Groundwater flow as a diffusion equation has been solved using LBM by several scholars. In one of the pioneer studies, the applicability of the LBM to isotropic groundwater flow modeling was investigated in confined aquifers (Zhou 2007a, b). Budinski et al. (2015) used LBM with D2Q9 square lattice for modeling groundwater flow. Furthermore, LBM was extended to simulate fluid flow in porous media (Gharibi and Ashrafizaadeh 2020). With regard to the solution of groundwater flow using LBM, more studies have been reported in Zhou 2007a;Zhou 2011;Budinski et al. 2015;Arsyad et al. 2017;Hekmatzadeh et al. 2019Hekmatzadeh et al. , 2020, and especially multiphase flow in porous media (Anwar and Sukop 2008;Fattahi et al. 2016;Gharibi and Ashrafizaadeh 2020).
The above-mentioned studies were considered confined groundwater equation, in which the dependent variable is the water head. Recently, Yousefi et al. (2021) derived an innovative LB solution for the unconfined groundwater flow equation. In unconfined groundwater aquifer, the dependent variable is the water head to the power of two. Severe groundwater exploitation for irrigation aims has led to significant land subsidence in several plains of Iran (Rahnama and moafi 2009;Motagh et al. 2007;Dehghani et al. 2013;Yousefi et al. 2019). Increasing groundwater pumping in an aquifer system leads to decrease pore water pressure while increasing the effective stress which leads to aquifer compaction and surface deformation (Mahmoudpour et al. 2016). On the other hand, land subsidence, which is caused by excessive groundwater exploitation, influences the aquifer parameters (Guzy and Malinowska 2020;Yousefi et al. 2019). Accordingly, the incorporation of the amount of subsidence into the inverse solution of the groundwater flow is of substantial importance.
Due to the superior advantages of LBM for the solution of time dependent problem, an innovative LB solution of unconfined groundwater flow is coupled with GA to determine the model parameters. To the best of the researchers' knowledge, this is the first time that LBM is combined with optimization technique for a groundwater inversion problem. Afterward, an innovative framework is proposed to predict land subsidence. The proposed method was applied on the data of Darab plain, located in Fars province, Iran. Darab plain has faced considerable land subsidence due to over exploitation in the last decade because of fine-grained sediments compaction occurred in the aquifer (Yousefi et al. 2019). Land subsidence definitely affects the hydrology and hydrogeology properties of aquifer system. Hence, incorporation of the amount of subsidence into the inverse solution of the groundwater flow appears essential. The important goals of this study are inverse solution of groundwater flow equation and determination of unconfined aquifer parameters (K x , S y , and α) taking into account land subsidence measurements and finally prediction of land subsidence in near future.
The aquifer parameters including hydraulic conductivity, specific yield and compressibility are annually estimated using the proposed method of inverse solution and the annual measurements of the subsidence. The temporal trends of the aquifer parameters are further investigated to predict the values of the parameters associated with a desired time. This process is followed by the prediction of corresponding subsidence which is finally compared to the observed value.

Case study
Darab plain is in southern part of Iran, Fars province ( Fig. 1), at a distance of 250 km from Shiraz city, with an area of 2400 km 2 (Fig. 1) surrounded by high mountains with the maximum height of 3126 m above mean sea level (MSL). The minimum height of Darab plain is 1050 m. The average temperature ranges from 14 to 22 °C, while the average annual rainfall is recorded as 245 mm. The aquifer system in this plain is unconfined and that a total of 3444 pumping wells are in the area. Over a period of 23 years , the groundwater level in the plain has dropped up to 26 m. In addition, the annual evaporation, infiltration of rainfall, and groundwater discharge in this area are approximately 120 million cubic meter (MCM), 11 and 431.36 MCM, respectively (Regional Water Company of Fars 2016).
The highest amount of groundwater exploitation is due to farmlands irrigation. Due to the decline in groundwater level, Darab plain is subject to land subsidence which is mostly happening in the agricultural lands in the plains (Fig. 1).
Agriculture is the main activity in Darab plain (see Fig. 1) in which groundwater is the main water source for irrigation. It should be pointed out that 17 piezometric wells are located in the study area and its surrounding. The piezometric water heights at the boundaries were estimated using Kriging method. Land subsidence happens due to the fine-grained sediments compaction as a result of increasing effective stress. The compaction occurred in the fine-grained interbeds, will in turn influence the aquifer parameters including specific yield, compressibility, and hydraulic conductivity. The annual subsidence rate was estimated as 18 cm/year, based on measurements (Yousefi et al. 2019). A precise knowledge of the aquifer parameters and their relationship with the compaction mechanism will enable us to predict the subsidence occurrence.
The primary objective of this paper is to present a novel method based on a combination of LBM and GA to solve the inverse groundwater flow equation for 2D problem in order to estimate the aquifer parameters precisely. To do so, initial values of the parameters and other required quantities are first extracted from the existing reports on groundwater information. Figure 2 indicates the initial amounts of hydraulic conductivity which ranges from 0.32 to 18.31 (m/day) in the center and north of the area, respectively (Regional Water Company of Fars 2016).
Hydraulic head extracted from the piezometric information of 2010 is also shown in Fig. 3. According to this figure, the north and south parts of the area are inflow and outflow boundaries (perpendicular to contour lines, from high values to low values of hydraulic head), respectively, while there is no flow border in east and west direction (parallel to contour lines).
The head contours drawn in this figure are obtained using the data from all available piezometric wells (inside and outside the border (Fig. 1)) and subsequently the boundary conditions defined for the model represent the reality in the area under study. Moreover, the thickness of Darab aquifer system obtained from the piezometric information is illustrated in Fig. 4. The maximum thickness of 250 m belongs to the center of the plain which decreases to 20 m toward the area boundaries. The geological formations of this area include 4 main groups: (1) Permeable formations located in the northern slopes of the plain, (2) Semi-permeable formations which contains medium-grained sediments that are located in the continuation of the permeable formations, (3) Low permeable formations which are observed in the southeast and southwest parts of the plain, (4) Poorlypermeable formations which consist of fine-grained sediments that are mostly located in the center of the plain (Fars Regional Water Joint Stock Company, 1356). According to this division, it is clear that Darab plain alluvium has different grain sizes in different parts of the area. Therefore, the sediments become finer from the north to the south and from the east to the west, and as a result, the aquifer parameters vary in different areas of the plain. (Regional Water Company of Fars 2016).  In previous studies, the temporal behavior and spatial extent of the subsidence has been investigated by Yousefi et al. (2019) using synthetic aperture radar interferometry (InSAR). Two different SAR satellite images, Envisat ASAR and Sentinel-1A (spanning 2010 and 2015-2017), have been employed to extract the subsidence rate. As suggested by both datasets, there are three main subsidence lobes with different deformation rate as depicted in Fig. 5 (Yousefi et al. 2019). Figure 5 represents that land subsidence in Darab plain contains three main zones extending from northwest to southeast, approximately in the center of the plain that correspond to the size of sediment particles (fine-grained sediment) and the location of water pumping wells. It is supposed to apply the subsidence maps to infer the aquifer parameters by inverse modeling of groundwater flow expressed in the LBM framework. The rationale behind the LBM and its use in modeling the groundwater flow is initially explained below followed by the inverse modeling workflow.

Lattice Boltzmann method (LBM) for modeling the groundwater flow
General form of the 2D groundwater flow equation in unconfined aquifer is defined as (Budinski et al. 2015): where K x and S y indicate the hydraulic conductivity and specific yield, h and R are the groundwater head and the recharge function, respectively, and x i shows two directions in the horizontal plane.
The lattice Boltzmann equation is stated as Eq. 2. It is proved that the groundwater equation is the same as LB equation. Basic form of 2D lattice Boltzmann method (Eq. 2) is used to achieve groundwater flow in unconfined aquifer by Yousefi et al. (2021): where f k and f k eq represent the particle distribution function and the local equilibrium distribution function, respectively, Δt is the time step, x is a vector in xy coordinate system, c k is te particle velocity vector, τ is the single relaxation time, and a is the number of the particle. LBM consists of two main steps: streaming and collision steps. The particles move to the neighboring lattice points in their directions of their velocities, in the streaming step, and the arriving particles interact with one another and change their velocity directions based on the scattering rules in collision step (Zhou 2004).
Indeed, the type of equilibrium distribution function (f k eq ) and the relaxation time (τ) changes Eq. 2 to other differential equations. Different types of Lattice, including D2Q4, D2Q5 and D2Q9, can be used for streaming in two dimensions (Hekmatzadeh et al. (  (Yousefi et al. 2021). In this method, it is assumed that the particles can move around on specific paths as shown in Fig. 7. Accordingly, a = 9 and c k, indicates the velocity and is described as follows; (c x = Δx/ Δt) (Zhou 2004): f k is also specified as: where h(x, t) and w are groundwater hydraulic head and weight factor, respectively, for each direction according to the next terms. The weight factors are indicated as: and hydraulic head for D2Q9 is calculated by Eq. 6: Assuming open boundary for the east boundary (unknown boundary) in Figs. 6, 7 and 8 is employed for calculating the distribution functions: Dirichlet boundary condition at the west boundary (known boundary) in Fig. 6 suggests Eq. 9 to measure the distribution functions: Regarding the case of Neumann Boundary condition (north boundary-inflow boundary in Fig. 6), the distribution function is estimated by Eq. 10 as follows: Apart from the boundary conditions, the following equations have been already introduced as the equilibrium distribution function for unconfined aquifers (Yousefi et al. 2021): Eq 11 is characterized by the following relations: Modeling the study area using LBM Application of boundary condition By considering Eqs. 3-14, the groundwater flow can be modeled using Eq. 2 which is a function of aquifer parameters, i.e., K x and S y as well as the aquifer compressibility (α). The initial values of K x , S y and α are first estimated using the piezometric measurements followed by calibrating the relaxation time. The optimized aquifer parameters are then estimated using GA by freezing the amount of relation time in the LBM equation. The whole process will be explained in more details in the following section.

Proposed methodology to estimate the aquifer parameters
The whole workflow applied for inverse modeling of groundwater flow for the determination of the aquifer parameters of the study area is represented in Fig. 8.
At first, the initial values of the parameters, K x , S y and α, should be determined using the piezometric information of the study area. Afterward, these values are extracted using interpolation for the whole area. For modeling of LBM, the first parameter to be calibrated is the relaxation time with a little range of variation. An approximate amount of aquifer parameters, K x , S y and α, would suffice to estimate the relaxation time. The first attempt, hence, is to estimate the compressibility, α, using the relationship between the groundwater level and the land subsidence (aquifer compaction) at piezometric wells location. Knowing about the compressibility would help approximate S y and the sediment types. In addition, a rough estimate of K x is already available from the existing reports. These values are further employed to solve the inverse groundwater flow equation using LBM which is combined with GA in order to estimate the precise values of aquifer parameters, K x , S y and α. In the following, the values of the parameters are determined annually, and by using the extracted temporal trend, the values are extrapolated for the near future. Consequently, the subsidence is predicted for the same time period of interest which is finally evaluated by the InSAR measurements of the subsidence. Different parts of the proposed method are explained in more detail in the following sections.

Estimation of approximate values of aquifer parameters using the relationship between land subsidence and piezometric information
In this study, in order to generate the proper model for Darab aquifer, the land subsidence obtained from SAR interferometry time series analysis from 2010 to 2016 (Yousefi et al. 2019) is considered as an effective factor. In order to incorporate the amount of subsidence measured by InSAR with aquifer parameters, the relationship between land subsidence and groundwater level drop in an unconfined aquifer is defined using one dimensional compressibility equation, (Eq. 14) (Hoffmann et al. 2003), which shows the ratio of the relative changes in the thickness of the aquifer to the changes in effective stress: where α is compressibility, db is the amount of subsidence, b is the thickness of aquifer and dσꞌ is the effective stress increase. By replacing the amount of effective stress (Eq. 16) in Eq. 15, the new equation for compressibility is obtained as follows (Eq. 17): where ρ is the water density, g is the gravity acceleration and dh is the hydraulic head difference which can be replaced by the difference between two piezometric measurements in a specific time period (Eq. 18).
By replacing Eq. 18 in Eq. 17 (Eq. 19), and rewriting it (Eq. 19) based on the head, we will have: According to Eq. 20, the hydraulic head is a function of compressibility where the compressibility is also affected by land subsidence (db is the thickness changes or amount of subsidence). To achieve the compressibility behavior of Darab aquifer system approximately, Eq. 20 is applied by incorporating the amount of subsidence and the variation of hydraulic head at piezometric wells with annual time step. The land subsidence is obtained from the results of InSAR (Yousefi et al. 2019). These results include the land subsidence time series from 2010 to 2017 which specifies the annual amount of subsidence. Based on the approximate values of compressibility, rough estimates of sediment types (S.T) and S y are also obtained (Table 1). S y is a dimensionless quantity, ranging between 0 and 0.4. Moreover, the maximum amount of compressibility occurs in the main subsidence zone (wells 3, 5, 9, 11), and consequently, the sediment type is mainly fine-grained as clay. On the other hand, despite the significant hydraulic head decrease at some piezometric wells (wells 1, 7, 10, 12), no remarkable subsidence occurs due to existence of coarse-grained sediment type such as gravel.

Inverse solution of groundwater flow equation based on LBM combined with GA
Equation 20 (the relation between land subsidence (db) and piezometric head (h)) is utilized for inverse modeling of the groundwater flow using LBM. Regarding this equation, h 2 is assumed to be computed using LBM by Eq. 6, while h 1 is the initial head. By substitution of h 2 as a final solution of groundwater flow equation based on LBM, Eq. 21 can be written as: In other words, the final solution of groundwater flow equation based on LBM should satisfy the above relationship.
Since the study area is large, the subsidence information was presented in a square grid with a cell size of 500 m (Yousefi et al. 2019). For this reason, the same cell size was used for the lattice size. The grid size of 250 m was also examined aiming to check the discretization accuracy. The estimated head values resulted (year 2010-2011) from cell size of 250 m are nearly identical to those stemmed from the grid size of 500 m ( Table 2). The computational time on inversion problem increases when an optimization tool is combined with a numerical solution. The higher the grid size, the faster optimization of model parameters is obtained. Of note, the discretization error is usually a matter of concern in the advection-diffusion equation especially when the advection term is dominated (Hekmatzadeh et al. 2019;Servan-Camas and Tsai 2009). In order to obtain the inverse solution of the groundwater flow equation, a gridded layer with a cell size of 500 m and thickness of the aquifer (Fig. 4), is generated. The relaxation time and the time step required for modeling are to be considered as 1.35 and 0.002 day, respectively. The relaxation time is estimated using approximate values of the aquifer parameters. Based on the hydraulic head contour lines of Fig. 3, the north part of the aquifer is inflow boundary with known hydraulic head (dirichlet boundary). The south boundary is, however, the outflow one whose hydraulic head is not known (open boundary), whereas the east and west borders are no flow boundaries (solid boundary). In addition to the recharge and discharge of the aquifer for the time period of 2010-2016 which is determined from the existing reports, the annual amount of hydraulic head extracted from the piezometric information are used in the process as well (Regional Water Company of Fars 2016). The piezometric head (h), corresponding to 2010 is considered as the initial head, h 1 . Moreover, the highest and lowest values of rainfall are 400 and 100 mm, at high altitudes and in the plains, respectively. The annual evaporation and annual infiltration of rainfall in the study area are considered to be about 67 and 30% which are equivalent to 120 and 11 MCM, respectively. According to the existing reports, there were 3444 active pumping wells in the study area in 2010, which corresponded to the annual groundwater discharge volume of about 431.36 MCM which was mostly utilized for agricultural irrigation.
Inverse modeling of groundwater flow with three unknown aquifer parameters including K, S y and α is performed with the incorporation of LBM and GA optimization. The objective function to be minimized in GA is the difference between the calculated hydraulic head from LBM (h (x, t), Eq. 6) and the one computed using the subsidence value (the right hand of Eq. 20, (h 2 )). In fact, the purpose of using the GA in this study is to determine the best values of the aquifer parameters, K , S y and α, in such a way that the value of the objective function is minimized. The values of the GA parameters used in this study are listed as follows: the maximum number of generation is 300, the size of population is 50, the migration fraction is 0.2, the crossover fraction is 0.8, the number of elites equal to 5% of the initial population and also other parameters are defaults. The lower and upper bounds of the aquifer parameters extracted from Table 1 and Fig. 2 are illustrated in Table 3.

Results and discussion
According to the flowchart of Fig. 8, the aquifer parameters, K , S y and α, are initially estimated based on the piezometric information. There are 17 piezometric wells available in the study area. These values are exploited in order to estimate the relaxation time which 1 3 is further used in the inverse solution of groundwater flow equation. By integration of LBM and GA and application of the amount of subsidence measurements as the innovative model of numerical optimization model to solve the inverse problem, K x, S y and α are optimized on the annual basis from 2010 to 2016 (Figs. 9, 10, 11) As observed in Figs. 9 and 10, the maximum and minimum amounts of K x and S y occur in the north (permeable formations with larger values of K x and S y ) and center of Darab plain (Poorly-permeable formations with smaller values of K x and S y ), respectively, that shows the agreement of the results with the geological features of the area. On the other hand, the maximum value of α is detected in the main zone of the subsidence which are consistent with the geological setting of the area (Fig. 11). The temporal decreasing trends of K , S y and α are clearly observed due to the progressive aquifer compaction. The optimized values of K , S y and α can be used effectively to calculate the hydraulic head at any location of interest. In order to evaluate the accuracy of the estimated aquifer parameters, the calculated hydraulic head is compared with that measured at the locations of piezometric wells. The annual relative errors between the calculated hydraulic head and the one observed at the locations of piezometric wells are indicated in Table 4.
According to Table 4, the maximum and minimum amounts of relative error are 0.0008 and 0.00001 which occur in wells 2, 5, 9, and 11, respectively. The small range of the relative error confirms that the results of the proposed model are strictly consistent with the observed measurements which is an indication of the accurate estimation of the aquifer parameters.
Since the main reason for the subsidence occurrence is the groundwater level decline, it is expected that any changes in the aquifer parameters K x , S y , and α might be related to the changes of the groundwater level. This is investigated at the location of piezometric wells. In Figs. 12, 13 and 14, the relationship between the piezometric head (h) and hydraulic conductivity (K x ), specific yield (S y ) and compressibility (α), are illustrated. These Figs contain RMSE values as the averaged model misfit.
As shown in Figs. 12, 13 and 14, the values of the parameters (K x , S y , and α) decrease with the reduction of the hydraulic head, in most parts of Darab plain that makes perfect sense. Water level decline changes the values of aquifer parameters and mainly decreases K x , S y . The decrease in water level is typically associated with the compression of the aquifer system. As the water level continues to decrease, its compressibility decreases as well.
In addition to the groundwater level decline, the subsidence occurrence is highly dependent on the aquifer system properties, specifically its hydrological and hydrogeological parameters. Therefore, one of the main objectives to optimize the aquifer parameters on the annual basis is to model their temporal behavior in order to be able to predict the amount of subsidence in the near future. This may simply enable us to study the subsidence behavior in the absence of satellite measurements. To this end, the temporal behavior of K x , S y and α at piezometric wells is first modeled as a function of time. The modeling results at the locations of piezometric wells are illustrated in Figs. 15, 16 and 17. It should be mentioned that the function employed for the temporal modeling is mainly based on the logarithmic family. Figure 15, 16 and 17 indicate that the temporal trends of K x , S y , and α, are mostly decreasing between 2010 and 2016. The values of K x and S y are, however, constant outside the subsidence zones (wells 4, 6, 10, and 12). The compressibility (α) at wells 2, 5, and 12 increases with time while K x and S y are constant or have a decreasing trend. 1 3 The slope of the line fitted to these parameters is also decreasing, which can be understood by considering the capacity of soil compressibility. The averaged model misfit is expressed as the root mean square error (RMSE) for each piezometric wells. The small values of RMSE are an indication of the proper curves fitted to the aquifer parameters time series.
By extrapolating the aquifer parameters, K x , S y , and α using the generated temporal model, we are able to estimate their values corresponding to the near future, i.e., 2017. This estimation is performed for every grid cell in the study area. This process is followed by the estimation of the subsidence in 2017 using LBM (Fig. 18). Meanwhile, the subsidence map corresponding to 2017 is produced by InSAR measurements (Yousefi et al. 2019) which will be served as the assessment criteria. Figure 19 shows the relative error as the model misfit calculated between the predicted subsidence and the one obtained from InSAR measurements.
The maximum relative error is estimated as 15%, which belongs to only 1% of all pixels at boundaries. The reason for the high relative error in the borders is the small amount of subsidence in these areas (since the absolute error of these areas is insignificant, the relative error would be high. It should be noted that the relative error is the ratio of the absolute error of a measurement to the actual measurement being taken). In the subsidence area, however, the amount of the model misfit is considerably smaller. The averaged relative error of the whole study area is estimated as 3.5% which is an indication of the high performance of the proposed subsidence prediction model. Moreover, investigation of the temporal behavior of the aquifer parameters enables us to predict the subsidence in near future which is one of the novel outcomes of the present study.

Conclusions
As an effective and novel framework for inverse solution of unconfined groundwater equation, LBM and GA optimization were jointly employed to determine the aquifer parameters including hydraulic conductivity (K x ), specific yield (S y ) and compressibility (α) considering land subsidence in Darab plain.
GA was used to determine the best values of the aquifer parameters with minimize the difference between the calculated hydraulic head by LBM and the one computed using the subsidence value as the objective function. Using the compressibility relation in unconfined aquifer systems, the relationship between land subsidence and groundwater equation was defined. The land subsidence has been already measured by InSAR. By considering the possible ranges of the aquifer parameters, the LBM combined with GA optimization was applied to achieve the aquifer parameters on the annual basis precisely. The results show that the values of hydraulic conductivity (K x ) and specific yield (S y ) decrease from the north to south and the maximum values of compressibility (α) are observed in the center of the plain. A model fitted to the temporal behavior of the aquifer parameters allowed us to figure out the hydrological and hydrogeological properties of the Drarb aquifer system and predict the subsidence in the near future.  The temporal trends of K x , S y , and α, between 2010 and 2016 suggest that the aquifer system is highly affected by the land subsidence. Furthermore, the small relative error between the predicted subsidence and the InSAR measurements demonstrates the   total storage capacity of the aquifer system which most probably leads to large surface runoff after a precipitation event. It means that the aquifer can no longer store water from the atmospheric precipitation. There is really a cause for severe concern specifically during drought events where people mostly rely on the groundwater resources. As a significant outcome of this research, precise information on the aquifer system properties and storability plays an important role for a successful groundwater management which is of important significance.