Land Subsidence Due to Groundwater Exploitation in Unconfined Aquifers: Experimental and Numerical Assessment with Computational Fluid Dynamics

Abstract

Land subsidence is a global challenge that enhances the vulnerability of aquifers where climate change and driving forces are occurring simultaneously. To comprehensively analyze this issue, integrated modeling tools are essential. This study advances the simulation of subsidence using Computational Fluid Dynamics (CFD); it assessed the effects of exploitation and recharge of groundwater on the vertical displacement of coarse and fine sands in a laboratory-scale aquifer. A model was developed by integrating the Navier–Stokes equations to study the groundwater flow and Terzaghi’s law for the vertical displacement of sands. The boundary conditions used were Dirichlet based on the changes in the hydraulic head over time. The specific storage coefficient was used to calibrate the model. The findings confirmed that subsidence occurs at slower rates in soil with fine sands with average particle diameters of 0.39 mm than in coarse sands with average particle diameters of 0.67 mm. The maximum discrepancy between the experimental and the numerical reaffirms that CFD platforms can be used to simulate subsidence dynamics and potentially allow the simultaneous simulation of other dynamics. Concluding remarks and recommendations are highlighted considering the up-to-date advances and future work to improve the research on subsidence in unconfined aquifers.

1. Introduction

The potential global exposure to land subsidence according to Herrera-García et al. [1] estimated a total of 31 countries with evidenced cases of subsidence in 2021 and up to 85 countries experiencing land subsidence in 2040. It is expected that by 2040 a total of 484 million inhabitants will be potentially threatened by land subsidence in an area of 12 million km² of the world, with a probability greater than 50%. This pressing situation asks for further investigations to prevent and manage land subsidence, especially when it is related to the exploitation of unconfined aquifers in coastal regions and threatens the sustainability of the freshwater supply.

Advances in Land Subsidence Studies in Aquifers

Subsidence due to the exploitation of fluids is the most studied type of subsidence. The first recollection of land subsidence due to fluid withdrawal was made by Poland and Davis [2], including the cases of subsidence due to the exploitation of oil, water, and gas. Important examples of land subsidence were assessed in the cities of Goose Creek, Texas [3]; Wilmington, California [4]; Lake Maracaibo, Venezuela [5]; Niigata, Japan [6]; Po Delta, Italy [7]; and lower Ganges-Brahmaputra Delta, Bangladesh [8], and cases of high exploitation of groundwater in Japan [9]; Mexico City (Mexico) [10]; and Arizona [11], Nevada [12], and California (USA) [13]. At that time, the San Joaquin Valley of California presented the highest levels of studied subsidence with up to 9 m from 1925 to 1977. The photographic evidence of the consequences of this rate of subsidence continues to bring attention to the matter nowadays, and it has become a benchmark for the great changes that subsidence can create to the landscape [14].

From 1992 onwards, radar images aquifers by SAR satellites have been used for InSAR studies as a supportive tool for land subsidence monitoring and to validate predictive models. This technique has shown effectiveness for the monitoring of subsidence in extensive regions [15,16]. Algorithms have also been developed to solve land subsidence problems assuming vertical compaction, such as MODFLOW SUB-WT and MODFLOW SUB [17,18], and with less frequency estimating numerical solutions with the application of Computational Fluid Dynamics (CFD). Both tools are based on the soil consolidation equations of Terzaghi’s law and Biot’s approach.

Major advances to assess and model land subsidence in coastal aquifers have occurred during the last decades in countries such as the USA, Japan, China, and throughout Europe [19,20,21,22,23,24]. Moreover, there is a variety of dynamics associated with land subsidence that require further study. The dynamics include the change of porosity and permeability due to soil compaction [25], subsidence in sandy aquifers [26,27,28], and the integration of subsidence and salinity intrusion in unconfined aquifers [29].

Table 1. Aquifers affected by land subsidence.

Aquifers and Study Area	Geological Characteristics a	Groundwater Exploitation Total or per Year	Data Collection	Type of Simulation	Maximum Subsidence	References
Chicot and Evangeline aquifer units (USA)	He, L, UNCO	5 m/year	Time series and InSAR	-	49 mm/year	[30]
Lower Bengal Delta (Bangladesh)	H, He, L, UNCO	5 m	Numerical simulation	Transient with MODFLOW	63 mm/year	[31]
Aguascalientes Valley (Mexico)	H, UNCO, Faults, R	3.5 m/year	SBAS InSAR	Fast Fourier Transform	120 mm/year	[10]
Morelia (Mexico)	H, UNCO, Faults, R	15 m/year	InSAR, settlement data	-	90 mm/year	[32]
Willcox Basin (Arizona)	He, L, CO	6 m/year Aprox.	InSAR, Hydraulic data	Storage loss estimation	140 mm/year	[11]
Wuxi City (China)	He, L, UNCO, CO	68 m	Extensometer	-	41.95 mm/year	[33]
Guangming Village (China)	He, L, UNCO, Faults, R	1 m	Experimental setup	Finite Element –Interfaced Elements	2–5 mm	[34]
Capo Colonna (Italy)	He, UNCO, Faults	-	InSAR, Geohazard Exploitation Platform	-	47 mm/year	[35]
Pingdu District (China)	H, UNCO	2400 mL/min	Experimental set up	-	0.708 mm	[36]
Chandigarh tri-city (India)	He, L, CO	0.2 m/year	InSAR, Field data	Neural network	8 mm/year	[37]
Bohai Bay (China)	He, L, CO	1.7 × 107 m3/year	InSAR, Field data	Neural network	80–150 mm/year	[23]
Fuhuayuan (FHY) deep foundation pit project	He, L, CO	0.24 m/year	Experimental set up	DEM–CFD	8.7 mm	[38]
Yangtze River Delta (China)	He, L, UNCO, CO	0.75 m/year	Experimental setup	-	7–10 mm	[39]
Xuwei area (China)	He, L, UNCO, CO	2.34 mm/month	Experimental setup	-	14.04 mm	[40]
Hypothetical aquifer	H, UNCO	Cycles of 27% exploitation and recharge	Experimental setup	CFD	2–4 mm	This study

Notes: ^a He: Heterogeneous aquifer, H: Homogeneous aquifer, L: Layered aquifer, Faults: Aquifer with evidenced faults, UNCO: Unconfined aquifer, CO: Confined aquifer, R: Analyzed recharge rates.

shows a summary of aquifers over the world that evidence land subsidence due to groundwater exploitation. It presents the various aquifer geological characteristics, groundwater exploitation, methods to collect subsidence data, type of simulation performed, and maximum evidenced subsidence.

Some of the mentioned studies and up-to-date advances consider surface deformation monitoring and field data to develop land subsidence models. Similarly, research about aquifers in Latin America includes field studies of the rates of land subsidence due to groundwater exploitation, and large-scale assessments [20,41,42]. However, comprehensive field campaigns involve extended periods of time and financial expenses. Mechanisms related to suffosion soil removal could be evidenced when assessing land subsidence in the field. It can manifest as a loss of volume or local erosion of the finer particles transported by seepage flow through the coarser particles or faults. It can result in soil instability and landslides [43]. Suffosion has not been linked to land subsidence in many studies; nevertheless, as mentioned by Najafi and Faghihmaleki, it could contribute to the enhancement of land subsidence in some specific cases [44]. Another area for enhancing land subsidence studies lies in the application of conventional numerical models. These models rely on assumptions related to the aquifer’s nature and primarily focus on groundwater exploitation as the driving force for land subsidence. This premise is applicable when land subsidence is the only evaluated threat for an aquifer. However, in the case of coastal aquifers where salinity intrusion poses an additional challenge, the application of conventional models is not enough. Various cases including Antonellini et al. [24], Eggleston and Pope [19], Yu and Michael [31], and Essink and Kooi [21] demonstrate the correlation between land subsidence and increased salinity concentrations in the aquifers. To analyze intrinsic factors that are affected by land subsidence and the changes that could affect salinity intrusion, the models require more robust capabilities such as the integration of various physics offered by CFD software.

To overcome the limitations of large-scale field studies and conventional numerical models, this study suggests the experimental recreation of land subsidence in a laboratory-scale setup subjected to exploitation and recharge. Lab-scale studies have proved to be adequate and cost-effective for the analysis of local dynamics and the study of small-scale subsidence dynamics. This study also considers CFD models to estimate aquifer compaction over time. It considers Terzaghi’s law as the general soil dynamics equation and couples it with the groundwater flow equation in terms of the aquifer storage coefficient. Applying CFD models allows the future coupling of land subsidence results with other dynamics such as salinity intrusion, consequently expanding the research scope of coastal aquifer numerical simulations. Their advantage over conventional models relies on the ability to integrate and solve partial differential equations representing simultaneous processes that occur in a specific control volume [45,46].

Two types of sand were analyzed for this research, namely coarse sand with an average particle diameter of 0.67 mm, also known as 12–20 grade sand, and fine sand with an average particle diameter of 0.39 mm, known as 30–40 grade sand, to analyze its vertical deformation after each cycle of exploitation and recharge. The experimental methodology followed the approach for the vertical deformation of sand suggested by Li et al. [28], where the vertical deformation was measured after each cycle of exploitation and recharge. Furthermore, the experimental findings were applied for the development of a numerical model coupling Navier–Stokes equations for groundwater flow and Terzaghi’s law.

This research provides a comprehensive theoretical background of land subsidence modeling, including mathematical approaches for simulating vertical displacement of the sands, groundwater dynamics, and changes in water pressure within the aquifer. The experimental setup and procedures are described in detail, highlighting the applied methodology and materials used. Additionally, the results and discussions from both experimental and numerical simulations are presented and analyzed. Finally, the study concludes with essential remarks and recommendations.

2. Theoretical Background

2.1. Mathematical Modeling of Land Subsidence

There are mainly two approaches in the mathematical modeling of land subsidence: Terzaghi’s law and Biot’s approach [47,48]. The former considers land subsidence as the vertical displacement resulting from compaction of the confining layers induced by water pressure variation, and the latter considers the land subsidence as a two- or three-dimensional problem, analyzing the consolidation of the soil and its vertical and horizontal movement due to water pressure changes and strain distribution.

Bear and Corapcioglu developed three mathematical models based on these approaches to simulate regional land subsidence induced by groundwater pumping [49,50,51]. They comprised equations for vertical displacements based on Terzaghi’s law [49], vertical and horizontal displacements [51], and equations for phreatic aquifers based on Biot’s approach [50]. Commonly used numerical models to simulate land subsidence are based on Terzaghi’s law.

Figure 1 illustrates the decline of the water table in an unconfined coastal aquifer, considering the effects of land subsidence and the interaction between the hydrostatic (p) and geostatic stresses (σ). Figure 1 also depicts that the water table is in a constant state and the hydrostatic stress and geostatic stress are not altered, showing a constant effective stress (σ’), and after the water exploitation occurs, and the hydrostatic stress is lowered, changing the geostatic stress, and increasing the effective stress.

Continuous cycles of water table decline and recharge can lead to a retained change in the effective stress even after the water table is replenished. This dynamic results in alterations to the aquifer soil and water dynamic. Consequent change to the effective stress strongly depends on the aquifer geology and the specific storage of the aquifer. The amount of subsidence is governed by various factors such as effective pressure, thickness, and compressibility of the aquifer, the length of time over which increased load is applied (for subsidence due to external loads), and the rate and type of stress [52]. The degree and rate of consolidation (vertical) and horizontal displacement depend on the stress-strain relationship of the solid matrix comprising the aquifer [9].

The mathematical premise states that a total load of soil and water (and everything that adds load at the ground surface, including atmospheric pressure) above the considered surface plane, is balanced by an effective intergranular stress (${\sigma }^{\prime }$) in the solid matrix of soil and by a pressure (p) in the water. This definition is described by Equation (1), which represents Terzaghi’s law, all terms with the magnitude of [ML⁻¹T⁻²].

$$\mathsf{\sigma }={\mathsf{\sigma }}^{\prime }+\mathrm{p}$$

(1)

To consider land subsidence in an aquifer, the solid matrix must be studied as compressible. The theoretical foundation to mathematically describe land subsidence includes two main stages: the fluid dynamic component and the structural deformation. The fluid dynamic component considers the mass balance of groundwater in the aquifer and expands the term for water accumulation. Equation (2) shows the abbreviated form of the groundwater momentum equation, with the fluid transport on the left side represented by Darcy’s law and the water accumulation term on the right side.

$$\nabla \left(\mathsf{\rho }\overrightarrow{\mathrm{v}}\right)=\frac{\partial \left(\mathsf{\rho }\mathsf{\phi }\right)}{\partial \mathrm{t}}$$

(2)

where ρ [ML⁻³], $\overrightarrow{\mathrm{v}}$ [LT⁻¹], $\mathsf{\phi }$ [-], and t [T] are the fluid density, flux velocity, porosity, and time, respectively. The momentum equation exposes that the difference in the fluid flow entering and exiting the cavity is equal to the accumulated mass of fluid in the cavity (negative when it is being accumulated), expressed in Equation (3).

$$\nabla \left(\mathsf{\rho }\overrightarrow{\mathrm{v}}\right)=-\frac{\partial \left(\mathsf{\rho }\mathsf{\phi }\right)}{\partial \mathrm{t}}$$

(3)

The land subsidence formulation also considers the influence of the water pore pressure and the intergranular or effective stress, with the consideration of the specific storage [53]. Equation (4) shows the expansion of the groundwater accumulation term to consider those factors.

$$\frac{\partial \left(\mathsf{\rho }\mathsf{\phi }\right)}{\partial \mathrm{t}}=\mathsf{\rho }\frac{\partial \mathsf{\phi }}{\partial \mathrm{t}}+\mathsf{\phi }\frac{\partial \mathsf{\rho }}{\partial \mathrm{t}}=\left(\mathsf{\rho }\frac{\partial \mathsf{\phi }}{\partial \mathrm{p}}+\mathsf{\phi }\frac{\partial \mathsf{\rho }}{\partial \mathrm{p}}\right)\frac{\partial \mathrm{p}}{\partial \mathrm{t}}$$

(4)

where ρ [ML⁻³], $\mathsf{\phi }$ [-], t [T], and p [ML⁻¹T⁻²] are the fluid density, porosity, time, and fluid pressure, respectively. Each term of the parenthesis in Equation (4) relates to the aquifer’s elastic properties. The first term considers the change of the porosity over time, and the second term is the fluid density change. Both terms are related to the variation of the water pressure. The change in the porosity and the density can be further described using Equation (5) and Equation (6), respectively. The coefficients $\alpha$ and $\beta$ represent the solid matrix compressibility and the pressure changes in the voids normally filled with water, both with units of [ML⁻¹T⁻²]⁻¹.

$$\mathsf{\alpha }=\frac{1}{1-\mathsf{\phi }}\frac{\partial \mathsf{\phi }}{\partial \mathrm{p}}$$

(5)

and

$$\mathsf{\beta }=\frac{1}{\mathsf{\rho }}\frac{\partial \mathsf{\rho }}{\partial \mathrm{p}}$$

(6)

Replacing Equations (5) and (6) in Equation (4) results in Equation (7).

$$\frac{\partial \left(\mathsf{\rho }\mathsf{\phi }\right)}{\partial \mathrm{t}}=\mathsf{\rho }\left[\mathsf{\alpha }\left(1-\mathsf{\phi }\right)+\mathsf{\beta }\mathsf{\phi }\right]\frac{\partial \mathrm{p}}{\partial \mathrm{t}}$$

(7)

Equation (7) can also be simplified, defining the term in square brackets as the specific storage S [LT²M⁻¹] in terms of water pressure. From this first stage, the spatial gradient of the fluid pore pressure is estimated. Equation (7) is used in the structural stage to estimate the compaction and the resulting subsidence of the ground surface.

$$\frac{\partial \left(\mathsf{\rho }\mathsf{\phi }\right)}{\partial \mathrm{t}}=\mathsf{\rho }\mathrm{S}\frac{\partial \mathrm{p}}{\partial \mathrm{t}}$$

(8)

Replacing Equation (8) in Equation (2), the groundwater mass balance equation takes the form of Equation (9).

$$\mathsf{\rho }\mathrm{S}\frac{\partial \mathrm{p}}{\partial \mathrm{t}}+\nabla ·(\mathsf{\rho }\overrightarrow{\mathrm{v}})=0$$

(9)

Furthermore, the structural deformation component considers soil mechanics and calculates the total compaction of layers of the aquifer based on soil parameters such as the void ratio, e, and the soil porosity, $\mathsf{\phi }$. The void ratio is defined as the ratio of the pore volume to the grain volume, and the porosity is the ratio of the pore volume to the total volume. The following relationship holds:

$$\mathrm{e}=\mathrm{P}\mathrm{o}\mathrm{r}\mathrm{e} \mathrm{v}\mathrm{o}\mathrm{l}\mathrm{u}\mathrm{m}\mathrm{e}/(\mathrm{T}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l} \mathrm{v}\mathrm{o}\mathrm{l}\mathrm{u}\mathrm{m}\mathrm{e}-\mathrm{P}\mathrm{o}\mathrm{r}\mathrm{e} \mathrm{v}\mathrm{o}\mathrm{l}\mathrm{u}\mathrm{m}\mathrm{e})=\frac{\mathsf{\phi }}{(1-\mathsf{\phi })}$$

(10)

Considering that the deformation occurs only in the void spaces of the porous media and that the grains are incompressible compared to the total porous volume. The compaction ($\mathsf{\eta }$) can be expressed by Equation (11), as it is a consequence of the reduction of the initial void space caused by an increase in the effective stress.

$$\mathsf{\eta }=\mathrm{b}\frac{\mathsf{\Delta }\mathrm{e}}{1+{\mathrm{e}}_0}$$

(11)

where $\mathrm{b}$ [L] and ${\mathrm{e}}_0$ [-] are the initial thickness of the layer and the initial void ratio, respectively. Equation (10) can be further expanded to account for the coefficient of compressibility and to be in function of the hydraulic head difference. Assuming one vertical stress and strain, the vertical ground displacement resulting in vertical stress can be expressed in terms of a change in groundwater head, described by Equation (12).

$$\mathsf{\eta }\approx {\mathrm{S}}^{\mathrm{*}}∆\mathrm{h}\approx \mathrm{b}{\mathrm{S}}_{\mathrm{h}}∆\mathrm{h}$$

(12)

where $\mathsf{\eta }$ [L] is the compaction, ${\mathrm{S}}^{\mathrm{*}}$ [-] is the storage coefficient, and ${\mathrm{S}}_{\mathrm{h}}$ the specific storage [L⁻¹] in terms of the hydraulic head.

2.2. Groundwater Flow Equations

This numerical model for land subsidence caused by overexploitation includes the groundwater flow equations. Equations (9) and (12) are used to estimate the aquifer’s water pressure based on the drawdown, and the consequent compaction using Terzaghi’s approach. The flux velocity in Equation (9) is calculated following Darcy’s law, as stated in Equation (13).

$$\overrightarrow{\mathrm{v}}=-\frac{\mathrm{k}}{\mathsf{\mu }}\left(\nabla \mathrm{p}-\mathsf{\rho }\overrightarrow{\mathrm{g}}\right)$$

(13)

where $\mathrm{k}$ [L²], $\mathsf{\mu }$ [ML⁻¹T⁻¹], and $\overrightarrow{\mathrm{g}}$ [LT⁻¹] are the intrinsic permeability of the aquifer, the viscosity of the water, and the acceleration due to gravity, respectively. In terms of hydraulic head, the flux velocity can be expressed as Equation (14).

$$\overrightarrow{\mathrm{v}}=-\frac{\mathrm{k}{\mathsf{\rho }}_{\mathrm{f}}\mathrm{g}}{\mathsf{\mu }}\left(\nabla \mathrm{h}+\frac{\left(\mathsf{\rho }-{\mathsf{\rho }}_{\mathrm{f}}\right)}{{\mathsf{\rho }}_{\mathrm{f}}}\nabla \mathrm{Z}\right)$$

(14)

where

$$\mathrm{K}=\frac{\mathrm{k}\mathsf{\rho }\mathrm{g}}{\mathsf{\mu }}$$

(15)

Equation (15) represents the hydraulic conductivity, considering a constant density of the water and replacing the fluid and porous media characteristics. Equation (9) can be expressed in terms of the hydraulic head as denoted by Equation (16):

$$\nabla ·\left(-\mathrm{K}\left(\nabla \mathrm{h}+\frac{\left(\mathsf{\rho }-{\mathsf{\rho }}_{\mathrm{f}}\right)}{{\mathsf{\rho }}_{\mathrm{f}}}\nabla \mathrm{Z}\right)\right)={\mathrm{S}}_{\mathrm{h}}\frac{\partial \mathrm{h}}{\partial \mathrm{t}}$$

(16)

In addition, considering that the density in the groundwater remains the same, Equation (16) can be expressed as Equation (17), also known as groundwater flow or continuity equation.

$$\nabla ·\left(-\mathrm{K}\left(\nabla \mathrm{h}\right)\right)={\mathrm{S}}_{\mathrm{h}}\frac{\partial \mathrm{h}}{\partial \mathrm{t}}$$

(17)

Recharge and Discharge of Groundwater

The relationship between strain and stress is estimated with the measurements of the change on the hydraulic head after each cycle of recharge and exploitation of water. The stress-strain behavior depends on a variety of factors, including the sand properties, the amount of stress, and the number and magnitude of cycles. The variation of the effective stress can be estimated by averaging the additional stress according to the sand thickness above and below the water table, similar to the work of Li et al. [28]. Equation (18) describes the variation of the effective stress with the hydraulic head change.

$$∆{\mathsf{\sigma }}^{\prime }=({∆{\mathsf{\sigma }}^{\prime }}_1·∆\mathrm{h}+{∆{\mathsf{\sigma }}^{\prime }}_2·{\mathrm{h}}_2)/\left(∆\mathrm{h}+{\mathrm{h}}_2\right)$$

(18)

where $∆{\mathsf{\sigma }}^{\prime }$ [ML⁻¹T⁻²] is the effective stress variation and ${\mathrm{h}}_2$ [L] is the water table depth after withdrawal. Considering the exploitation and recharge cycles, Equation (18) can be modified for Equation (19), where γ [ML⁻¹T⁻²] is the unit weight of water.

$$\begin{gathered} ∆{\mathsf{\sigma }}^{\prime }=\frac{\mathsf{\gamma }·∆\mathrm{h}+{∆{\mathsf{\sigma }}^{\prime }}_2·{\mathrm{h}}_2}{\left(∆\mathrm{h}+{\mathrm{h}}_2\right)} \mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\mathrm{d}\mathrm{r}\mathrm{a}\mathrm{w}\mathrm{a}\mathrm{l} \\ ∆{\mathsf{\sigma }}^{\prime }=\frac{-\mathsf{\gamma }·∆\mathrm{h}+{∆{\mathsf{\sigma }}^{\prime }}_2·{\mathrm{h}}_2}{\left(∆\mathrm{h}+{\mathrm{h}}_2\right)} \mathrm{r}\mathrm{e}\mathrm{c}\mathrm{h}\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{e} \end{gathered}$$

(19)

3. Materials and Methods

The experimental setup is shown in Figure 2. It consists of a rectangular tank with 2.70 m (length) × 1.25 m (height) × 0.10 m (width) divided into three compartments (3, 7, 10) a central compartment containing the porous media (2.50 m long) and two lateral compartments (0.10 m long) for observing the water level. A fine screen mesh (4, 9) with holes of 0.074 mm (also known as mesh No. 200) separates the lateral compartments from the central compartment, preventing sand from flowing from the central into the side ones. The sand was silica sand, specifically quarzitic sand with 94.7% silica content. Fine and coarse sand with effective diameters of 0.39 mm and 0.67 mm, respectively, were used.

There are also two side chambers (2, 11) that control the water level. These chambers are connected to the tank through pipes and transparent flexible tubes. The tank is made of transparent Plexiglas with a thickness of 15 mm and is supported by a reinforced steel frame. Steel tension pins were also added to strengthen the tank and prevent side deformation. Three vertical displacement transducers were located at the top of the central compartment separated 1.25 m from each other (8). The displacement transducers (LVDT) were the Model Humbolt type HM-2310.04 range 0.4″ (10 mm) with a precision of 0.001 mm. Two side reservoirs (1, 14) with their respective pumps (5, 13) were included to store freshwater and supply the lateral chambers, the left tank with a volume of 250 L and the right tank with a volume of 500 L.

3.1. Experimental Procedure

The first step was to pack the sand simultaneously with the water on the lateral compartments. To guarantee the uniformity of the grain size and avoid the formation of air spaces and undesired strata, a packing technique was applied. It consisted of filling 10 cm of water prior to adding 10 cm of sand, this process was repeated until the amount of soil reached 1.10 cm. This process was applied for both coarse and fine sand.

The subsidence process was measured every 24-h using both analog and digital longitudinal vertical displacement transducers (LVDT) type HM-2310.04 Modelo Humbolt range 0.4″ (10 mm) and a precision of 0.001 mm. The exploitation was simulated by lowering the water level from 1.1 to 0.8 m in chambers (2) and (11), and the recharge was simulated by increasing the water level from 0.8 to 1.1 m. To initiate the exploitation cycles (lowering and increasing the water level), the hydraulic head was set to its initial height of 1.1 m. The surface level of the sand was registered at the beginning of the experiment, and it was considered the reference level to measure the vertical displacement of the sand. Hereafter, the hydraulic head was lowered to 0.8 m and this level was maintained for 24 h.

The subsequent exploitation cycle began 24 h after the previous one, lowering the water table by 30 cm. A total of four scenarios of exploitations, recharge, and stabilization cycles were conducted, including two for fine sand and two scenarios for coarse sand. The selected drawdown was considered based on the evidenced depletion of unconfined aquifers subjected to overexploitation. Among the considered drawdowns are a 50% drawdown of the Lower Bengal Delta [31], a 43% drawdown of the Aguascalientes Valley, 42% at the Wuxi city aquifer, and an average drawdown of 27.66% of the total depth of the Arroyo Grande Aquifer, an unconfined sandy aquifer in Cartagena, Colombia. Therefore, for this study, the considered drawdown was 30 cm, also equivalent to a 27% drawdown of the experimental setup.

3.2. Numerical Modeling of Land Subsidence

The numerical model used the software COMSOL Multiphysics 6.0 under the assumptions of one-dimensional compaction and no additional sources of stress. Consequently, Equation (20) is obtained and applied for the estimation of the aquifer vertical displacement.

$$\frac{\partial }{\partial \mathrm{y}}\left(\mathrm{K}\frac{\partial \mathrm{h}}{\partial \mathrm{y}}\right)={\mathrm{S}}_{\mathrm{h}}\frac{\partial \mathrm{h}}{\partial \mathrm{t}}$$

(20)

where $\mathrm{K}$, $\frac{\partial \mathrm{h}}{\partial \mathrm{y}}$, ${\mathrm{S}}_{\mathrm{h}}$, and $\frac{\partial \mathrm{h}}{\partial \mathrm{t}}$, are the hydraulic conductivity in the vertical axis, the hydraulic head in the vertical axis, and the specific storage in terms of hydraulic head and the hydraulic head over time, respectively. Simulation with Terzaghi’s approach involves using the skeletal-specific storage or aquifer compressibility to calculate the vertical compaction $\eta$ (m) with Equation (21). This compaction is one-dimensional and follows a conventional flow model (Equation (9)). The results are used in post-processing to calculate vertical compaction based on Terzaghi’s theory [54].

$$\mathsf{\eta }={\mathrm{S}}_{\mathrm{h}}\mathrm{b}(-\mathrm{h})$$

(21)

where b is the standard notation for the vertical thickness of aquifer sediments [L] and $\mathrm{h}$ [L] is the hydraulic head.

3.3. Boundary Conditions

The selected boundary conditions for solving the set of equations using Terzaghi’s approach are as follows:

No-flow boundary conditions for the top and bottom boundaries are expressed by Equation (22).

$$-\mathsf{\rho }\overrightarrow{{\mathrm{v}}_0}=0$$

(22)

Lateral boundary conditions are head-controlled, with a change of the hydraulic head over time, as described by Equation (23).

$$\mathrm{p}={\mathrm{p}}_0-\mathsf{\rho }\mathrm{g}\mathrm{h}$$

(23)

where $\mathrm{h}$ is the hydraulic head at the boundary variating over time. Figure 3 shows the recreated scenarios; fine sand with the first scenario (A) and second scenario (B), and coarse sand with the first scenario (C) and second scenario (D). The variation of the hydraulic head over time represents the change of the water table level in the left- and right-side compartments used to recreate the cycles of exploitation and recharge.

3.4. Parameters Considered for the Simulation

The numerical simulation of the land subsidence experiment considered the initial parameters presented in

Table 2. Parameters for the numerical simulation.

Variable	Description	Value
$\rho$	Fluid density	1000 kg/m3
p	Fluid pressure	0 Pa
$h\left(0\right)$	Initial hydraulic head	1.1 m
$h\left(t\right)$	Hydraulic head over time	f(t)

, along with the physical and mechanical properties summarized in

Table 3. Physical and mechanical properties of the sands.

Sand Type	Specific Gravity, Gs	Effective Diameter (mm)	Porosity (%)	Hydraulic Conductivity, K (m/s)
Fine sand	2.65	0.39	43.3	2 × 10−4
Coarse sand	2.74	0.67	48.8	6.5 × 10−4

. The grain distribution and classification of the sand are depicted in Figure 4.

4. Results and Discussion

4.1. Experimental Results

This section presents the experimental results of the land subsidence laboratory scale simulation conducted to investigate the vertical displacement due to water exploitation and recharge in sands. The experiment involved subjecting the aquifer to recharge and exploitation cycles, with varying hydraulic heads. The measurements were obtained from three sensors, and the findings highlight the influence of cyclic exploitation on the pore pressure and subsequent changes in effective stress, in accordance with Terzaghi’s law.

4.1.1. Coarse Sand

Figure 5 depicts the vertical deformation of the coarse sand in two scenarios of water exploitation and recharge. The cycles are depicted with gray bars; the hydraulic head was changed with continued cycles of water tables of 1.1 m and 0.8 m. The first scenario depicted in Figure 5a showed an average maximum total displacement of 2.57 mm after 19 cycles and 692 h while the second scenario depicted in Figure 5b showed a maximum vertical deformation of 3.7 mm after 625.4 h and 15 cycles. The second scenario had a constant hydraulic head from 0 to 200 h; during that time, there was no vertical deformation. After the first cycle of exploitation was recreated, the vertical deformation increased, as observed after 200 h of the experiment.

The measurements were obtained from three sensors; during the first scenario, sensor 2 presented technical issues and therefore the measurements were neglected. For the second scenario, the issue was solved, and the three measurements were considered. The correlation coefficient for the measurements in the first and second scenarios was 0.99, suggesting a strong linear relationship.

Terzaghi’s law explains that an increase in effective stress within a soil or aquifer leads to enhanced compaction and deformation [55]. This is evidenced in the observed increase in vertical displacement during the cyclic exploitation of the coarse sand. As the sand underwent repeated cycles of recharge and exploitation, the change in the hydraulic head altered the pore pressure, affecting the effective stress on the sands. Consequently, compaction occurred, and it resulted in increased vertical displacement of the surface, as evidenced in Figure 5.

Figure 6 shows the continued variation of the effective stress generated with the cycles of exploitation and the resulting vertical strain for the coarse sand. The first scenario (Figure 6a) showed a continued increase and decrease in effective stress with the water table changes, this being compared to the constant increase in strain. The second scenario (Figure 6b) showed that the effective stress was constant, and the strain was low when the water table is stable, beginning to increase with the subsequent changes in the water table. This behavior represents the continued tendency of the coarse sand to vertically deform over time with the increase and decrease in the water table. This has been previously evidenced in laboratory scale experiments [56,57] and field data [58] where only a small proportion of compression was recovered in sandy aquifers after the groundwater level was replenished.

The strain slope for scenario 1 remains relatively constant over time, indicating consistent exploitation and recharge levels. However, a notable change in slope occurred after 315 h when the hydraulic head was reduced to 0.4 m (as shown in Figure 6a), increasing the strain slope. For scenario 2 (Figure 6b) the strain initially exhibited an average constant value of 1.8 × 10⁻⁴ up to 212 h, corresponding to a steady water level of 1 m. Upon the beginning of the first exploitation, the water level decreased to 0.8 m, causing the strain to increase to 1 × 10⁻³. The strain then remained constant as long as the water level was maintained and subsequently increased to 2 × 10⁻³ after the second cycle of exploitation. After the second cycle starting at 358 h, the strain slope remained stable, leading to progressive vertical deformation of the soil with each consequent cycle of recharge and discharge until reaching a maximum average strain of 3.4 × 10⁻³.

The first scenario of strain shows a constant increase in vertical deformation after approximately 150 h. The maximum strain for the second scenario was higher than the maximum strain for the first scenario due to an increase in deformation in the second scenario after 190 h. This behavior agrees with the constant rates of exploitation and recharges in the first scenario, different from the second scenario where the cycles were not constant, as depicted in Figure 6b. The change in the slope of the strain over time shows that the continued recharge and discharge of water in the coarse sand contributes to the densification of the sand, reducing the strain slope with the continued cycles, alike previous work from Li et al. [28].

The two repetitions of the coarse sand exhibited variations in strain, which can be attributed to the initial anisotropy of the samples caused by horizontally oriented grains or result from the sand packing technique [28]. The differences in the cycles of exploitation are a relevant factor that contributes to the strain difference. It was observed that the strain increased during water exploitation, which is also evident in Figure 5. Our experimental findings confirm the application of Terzaghi’s law to analyze the dynamic behavior of land subsidence in coarse sands. The findings with the experimental results of coarse sand highlight the significance of considering the effects of cyclic exploitation on effective stress and subsequent deformation.

4.1.2. Fine Sand

The fine sand scenarios showed a decrease in the rate of subsidence, with a vertical displacement that occurred mostly because of a delayed response to water table changes. The cycles of recharge and discharge for the fine sand were conducted in the same manner as with the coarse sand; however, the capillarity forces with the fine sand held more saturation than the coarse sand after the water table was lowered. Figure 7 shows the location of the capillarity fringe when the water table is fixed at 0.8 m for side A with the fine sand, and side B with the coarse sand. The capillary fringe occurs when a percentage of saturation is trapped in the micropores of the sand, depending on whether the porosity is high or low. Sands with a higher porosity have less capillary force than finer sands [53].

Figure 8 shows the vertical displacement in fine sand; it occurred with periods of no vertical displacement and periods of continued displacement, with a delayed response to the water table change, initiating vertical deformation after the second cycle of exploitation. This behavior agrees with the vertical deformation of dense sand found by Li et al. [28]. According to Li et al., for the first and second cycles the sand behaved as an elastic material, not evidencing vertical deformation. Consequently, with the increase in exploitation and recharge cycles, the vertical deformation increased. This indicates a non-recoverable plastic behavior of the sand, and continued densification of the soil until the deformation reduces with the number of cycles. The maximum vertical displacement for scenario 1 was 2.6 mm after 425.3 h, equivalent to 19 days of experiment, with a total of 12 cycles of recharge and discharge, and for scenario 2 the maximum displacement was 5.06 mm after three months of the experiment, with a total of 2404.8 h and 38 cycles of recharge and discharge. In scenario 2, at 424 h and 12 cycles, the maximum vertical displacement was 2 mm, evidencing a difference of 0.6 mm from the first scenario.

The correlation coefficient between the three measurements for scenario 1 was 0.97. During the second scenario, some of the measurements from the sensors were neglected due to errors in the data readings. These data reading issues were attributed to electricity failures in the facilities. After the data-cleaning process, it was observed that the three sensors followed the same tendency, maintaining a correlation coefficient of 0.96. This coefficient evidences a close agreement between the data measurements. The tendency of the fine sand shows that with the continued cycles the vertical deformation increases and sustains a constant slope.

The strain over time shows an increase after the cycles of exploitation, similar to the deformation with coarse sand. Figure 9 shows the strain and variation of effective stress over time, occurring after the cycles of exploitation of the fine sand. The strain progress over time shows a more subtle strain slope than the coarse sand scenarios, with gradual increases resulting from a delayed response to the water table change and lower compaction rates. Zhang et al. [57] obtained a similar experimental response to the simulation of fine sands subjected to intermittent exploitation, with a strain that increased with the exploitation rates and then reached a period of low deformation.

Figure 9a shows that for the first scenario of fine sand, the strain continued to increase when the stress was constant, and it did not rapidly increase with the change in effective stress. After the effective stress increased with the lowering of the hydraulic head at 156 h, and then recovered at 181.9 h, the change in strain reduced from 0.00109 to 0.0011. However, after this cycle and before the next exploitation, the strain continued increasing, evidencing the slower response of the vertical deformation after the exploitation cycles.

Figure 9 depicts the strain and the effective stress variation resulting from the changes in the hydraulic head presented in Figure 3. The relationship between effective stress variation and the hydraulic head changes are explained in Equation (19). It introduces that during the cycles of discharge, there is an increase in effective stress, and a decrease in effective stress during recharge. The increase in effective stress occurs due to the drainage of water from the thickness of the soil that was previously saturated, and that is no longer under the influence of pore water pressure. The opposite is evidenced by an increase in the water table, increasing the pore water pressure, and reducing the effective stress. This estimation has been previously applied in the works of Li et al. [28] and presented in the theory of land subsidence [53].

Figure 9b evidenced that for the second scenario of fine sand, the beginning of the vertical displacement was also delayed. After two cycles of exploitation and recharge the first readings of displacement were observed. Four moments of apparent recovery of the sand were evidenced at 703.5 h, 848.1 h, 919.5 h, and 1063 h, which show a lowering of the strain followed by an increase with the next cycle. This behavior has been observed in previous studies analyzing sandy strata, including observations of delayed land subsidence in aquifers subjected to artificial recharge and exploitation [39], the rebound of strain after the cycles of drainage and recharge in an experimental study [40], and the numerical simulation of groundwater exploitation and recharge [59]. It also shows an initial increase in strain after the effective stress was reduced, presenting a similar behavior to the first scenario of fine sand presented in Figure 9a.

4.2. Numerical Results

After the experiments were developed, a numerical model was designed to explore the capabilities of CFD platforms to recreate subsidence experimental results in sands. To guarantee the numerical stability of the results a mesh convergence assessment was performed. The selected meshes with their respective degrees of freedom and number of elements are presented in

Table 4. Selected meshes for the mesh sensitivity assessment.

Mesh	DOF	Number of Elements
Coarse (fluid dynamics)	25,685	12,658
Normal (fluid dynamics)	170,799	28,192
Finer (fluid dynamics)	435,543	72,152
Extra fine (fluid dynamics)	689,835	343,972

. The estimation of compaction over time for the first scenario of coarse sand was selected for the mesh assessment and the results of the assessment are presented in Figure 10. The resulting compaction showed a close agreement with the results with each mesh. However, after 324 h, 336 h, and 360 h, the coarser and normal meshes showed a disagreement with the value of compaction. The more stable estimation of compaction was obtained with the finer to the extra fine mesh. The average standard deviation between the coarse mesh and the finer mesh was 0.35% and the finer mesh and the extra fine did not present deviation over time. For the estimation of compaction, the mesh with the best performance was the finer mesh with 72,152 elements and a maximum element size of 0.0308. The minimum element quality of the mesh was 0.085, the average element quality was 0.1224, and the geometry of the elements was triangular.

The Courant number was considered as a target to select the element size and the mesh. Due to the slower velocity and considering the mesh sensitivity assessment performed with the first scenario of coarse sand, the finer mesh with 72,152 elements was selected. In this research, the backward differentiation formula scheme was applied, considering that it has been shown to provide stability to problems that involve groundwater transport [60]. The free time stepping option was also employed, which allows the time step to be adjusted automatically based on local error estimates. For more detailed information about the time discretization scheme used in COMSOL Multiphysics 6.0, please refer to the COMSOL Multiphysics User’s Guide [60].

The hydraulic head was a key input parameter to define the variation of the water level in the numerical model. The numerical model considered all experimentally acquired data, such as initial porosity estimated using the Standard Test Methods For Specific Gravity Of Soil Solids By Water Pycnometer ASTM D854-00, the Hydraulic Conductivity estimated using a Permeameter ASTM D5084-16, and the soil grain size distribution estimated using the Unified Soil Classification System ASTM D2487-17e1.

The resulting vertical displacement obtained by the experimental results was used to compare and calibrate the subsidence results from the numerical model. COMSOL Multiphysics 6.0 was used to solve vertical displacement with Terzaghi’s law. The specific storage was selected as the calibration parameter and based on the range of specific storage suggested by Kuang et al. [61] for unconfined aquifers and the estimated average specific storage considering the average compaction per hydraulic head changes.

4.2.1. Coarse Sand

Using the finer mesh, the first and second scenarios of coarse sand were simulated, and the results are depicted in Figure 11. The numerical simulations were performed with Terzaghi’s Equation (12) and the groundwater flow equation in terms of storage coefficient (Equation (17)). The numerical results obtained a close agreement with the confidence interval obtained with the experimental results; this interval is depicted in gray shading, the selected specific storage for the numerical model that obtained the best agreement with the experiments is depicted with the black continued line, and the second closest is presented with the black hyphened line.

The specific storage was selected within a range of specific storages for coarse sands with low grade and hydraulic conductivity close to 0.00067 m/s. The evaluated range was 9 × 10⁻⁴ m⁻¹, 7 × 10⁻⁴ m⁻¹, 5 × 10⁻⁴ m⁻¹, 3.5 × 10⁻⁴ m⁻¹, 3 × 10⁻⁴ m⁻¹, and 1.91 × 10⁻⁴ m⁻¹ according to Kuang et al. [61].

The specific storage that was better adjusted to the first scenario was 3.5 × 10⁻⁴ m⁻¹ and for the second scenario was 9 × 10⁻⁴ m⁻¹. The average absolute error between the numerical model results and the media of the experimental results for the first scenario was 0.14 mm, with the highest discrepancy when estimating vertical deformation at 668.9 h and 692.35 h of 0.31 mm and 0.317 mm, respectively. For the second scenario, the average absolute error was 0.34 mm with the highest discrepancy at 309.8 h with 0.93 mm.

The first scenario of coarse sand (Figure 11a) showed the maximum vertical displacement obtained by the experiment was 2.8 mm by sensor two and 2.57 mm and 2.56 mm obtained by sensors one and three. The maximum displacement obtained by the model with the specific storage of 3 × 10⁻⁴ m⁻¹ was 2.01 mm and 2.6 mm with the storage of 3.5 × 10⁻⁴ m⁻¹. From the beginning of the experiment until 500 h, the specific storage of 3 × 10⁻⁴ m⁻¹ adjusted with the media of the experimental results and after 507.9 h the vertical displacement adjusted better with storage of 3.5 × 10⁻⁴ m⁻¹.

For the second case of coarse sand (Figure 11b), it is observed both in the experimental data and the numerical results that the initial displacement was delayed until the lowering of the water table started. The maximum vertical displacement according to the experimental confidence interval was 3.86 mm and the minimum was 3.61 mm, while the maximum value obtained by the numerical model was 4.19 mm. The range of specific storages that were evaluated for this scenario of coarse sand was 9 × 10⁻⁴ m⁻¹, 7 × 10⁻⁴ m⁻¹, and 1 × 10⁻³ m⁻¹. The storage that was better adjusted was the 9 × 10⁻⁴ m⁻¹ with an average absolute error to the experimental media of 0.34 mm with the highest discrepancy at 309.8 h with 0.93 mm. The specific storage of 1 × 10⁻³ m⁻¹ had a closer agreement to the media of the experimental results after 300 h and after 457 h the best agreement was obtained with 9 × 10⁻⁴ m⁻¹.

The numerical results and the experimental results for coarse sand showed a constant tendency of increasing vertical displacement; over time, it was evidenced that the slope started to decrease with the continued cycles of exploitation.

4.2.2. Fine Sand

Figure 12 depicts the vertical displacement over time for the fine sand. The evaluated range of specific storages for the fine sands included 4 × 10⁻⁴ m⁻¹, 4.5 × 10⁻⁴ m⁻¹, 5 × 10⁻⁴ m⁻¹, and 6 × 10⁻⁴ m⁻¹. For the first scenario of fine sand (Figure 12a), the best agreement was obtained with the storage coefficient of 5 × 10⁻⁴ m⁻¹ and for the second scenario was 4.5 × 10⁻⁴ m⁻¹. The average absolute error for the first scenario was 0.33 mm with respect to the media of the experimental results. From the beginning of the experiment up to 279 h the best agreement was with 5 × 10⁻⁴ m⁻¹ and after 279 h the best agreement was with the coefficient of 6 × 10⁻⁴ m⁻¹ with an average error of 0.34 mm. The maximum vertical deformation obtained with 6 × 10⁻⁴ m⁻¹ was 2.5 mm, a difference of 0.18 mm with respect to the experiment. The average absolute error for the second scenario (Figure 12b), was 0.29 mm with respect to the media of the experiments. The agreement between the experimental and numerical results for the second scenario shows the potential of integrative tools such as CFD to simulate land subsidence. This scenario lasted more than three months and the simulation recreated the tendency of vertical deformation observed in the experiments.

Along the different periods, the vertical displacement showed a close increase in subsidence over time. The simulation shows a constant slope of compaction, but for the experimental case, periods with low to no compaction were observed. This was evidenced in the first scenario after 200 h and at 300 h again, with constant cycles of exploitation and recharge, and in the second scenario after 250 h and at 500 h. The late response of finer sediments has been evidenced by various studies [39,40,57] and one of its consequences is delayed land subsidence after the water head declines.

From all experimental and simulated results, it is observed that the tendency of land subsidence in sands is to increase with each cycle of exploitation and recharge. The strain slope reduces over time and, according to the literature, it continues to deform until it reaches a maximum deformation and reduction of the specific storage [62]. In this study, the maximum vertical displacement for all scenarios was not reached, but a vertical displacement of sand of up to 5 mm was recreated with the longest experiment, evidencing that high displacement can occur in sands subjected to continued exploitation rates. The highest slopes for vertical deformation were obtained with the coarse sand, with an estimated deformation after 600 h of 4.19 mm for the simulation and an average of 3.7 mm with the experimental results. The tendency of the results agrees with Wu et al. [62] and Lv et al. [27], reaffirming that sands can rearrange their pores with each cycle of exploitation and recharge, consequently generating vertical displacement of the surface.

The specific storage of an aquifer indicates the percentage of water storage alteration corresponding to a unit change in the hydraulic head within the aquifer. Considering the extensive collection of data and references for unconfined aquifers presented by Kuang et al. [61], the specific storage for unconfined aquifers is expected to be distributed in the range of 10⁵–10⁻³ m⁻¹, with the highest frequency for unconfined aquifers of 10⁻⁵–10⁻⁴ m⁻¹. The specific storages evaluated and considered to adjust to the numerical model agree with the range provided by Kuang et al. [61], with storage coefficients of 3 × 10⁻⁴ m⁻¹, 9 × 10⁻⁴ m⁻¹, 3.9 × 10⁻⁴ m⁻¹, and 5 × 10⁻⁴ m⁻¹ for coarse and fine sand, respectively.

5. Conclusions

The understanding of land subsidence impact in unconfined aquifers is key to improving the estimation of the driving forces that act upon aquifers and threaten their sustainability. Among those forces, salinity intrusion and the transport of anthropogenic pollutants are threats that can be altered due to specific storage changes and land subsidence. This study incorporates its novelty into the estimation of land subsidence based on Terzaghi’s approach using Computational Fluid Dynamics, and the simulation of scenarios of land subsidence for up to three months in coarse and fine sands using an experimental setup. Relevant remarks are highlighted in the following statements:

• During the continued cycles of recharge and exploitation, both sands showed continued compaction that kept increasing over time. However, the deformation of coarse sand occurred at higher rates than for the fine sand. At times when the hydraulic head was maintained constant, the deformation for the coarse sand was significantly reduced, and for fine sands, it evidenced a delayed response;
• The estimation of land subsidence following Terzaghi’s approach agreed with the vertical displacement behavior observed in laboratory experiments. Nevertheless, it required the adjustment of the model considering the aquifer’s specific storage and detailed characterization of its physical characteristics;
• The average absolute error of the numerical model to recreate the experimental results of land subsidence was 0.14 mm for the first scenario of coarse sand, 0.34 mm for the second scenario of coarse sand, 0.33 mm first the first scenario of fine sand, and 0.29 mm for the second scenario of fine sand. The maximum discrepancy of 0.34 mm was obtained for the second scenario of coarse sand where the maximum experimental subsidence was 3.86 mm;
• The variation of effective stress and strain in sand subjected to withdrawal-recharging cycles revealed different patterns of deformation behavior with grain sizes and distributions. Both types of sand showed mostly inelastic deformation that did not rebound, except the fine sand, which evidenced three rebounds after a month of cycles of exploitation and recharge;
• The experimental results evidenced that in addition to the specific storage that contributes to deformation, fine and coarse sands show a different response to capillarity effects due to their different effective diameter and porosity. Therefore, unconfined aquifers with fine sands are able to maintain more saturation within the micropores after the water table is lowered, contributing to a slower or delayed response to the water table changes;
• The variation of the specific storage in unconfined aquifers is a pressing topic that requires further study. The consequent reduction in storage coefficient could alter the groundwater flow and other simultaneous dynamics occurring in the aquifer.

Large-scale simulation of land subsidence in unconfined aquifers based on satellite and in situ data is suggested to test more robust subsidence models with the application of CFD. Countries with abundant coastal aquifers and groundwater resources such as Colombia require further research regarding land subsidence, especially in coastal zones where it has been evidenced and where there is a potential for groundwater supply with the coastal groundwater reserves.

Future work could oversee the microscopical rearrangement of pores in unconfined aquifers due to constant cycles of recharge and discharge in fine sands, including its relationship to the porosity and possible increase in the aquifer’s capillarity. Assessments should also consider the land subsidence distribution in highly heterogeneous and layered fine-sand aquifers. The relationship between land subsidence and other mechanisms such as suffosion is an interesting subject for future study and will be considered by the authors during future work that will be undertaken.

Additionally, research could consider the change of storage coefficient in unconfined aquifers due to land subsidence and its relationship with the changes in groundwater flow in the aquifer. Other approaches could explore large-scale land subsidence with CFD considering Biot’s approach and include the estimation of geotechnical properties of the soil, such as Young’s modulus and the Poisson’s ratio. The authors will continue advancing this topic to reduce the existing knowledge gaps.