Regional subsidence at the former Texcoco Lake: numerical modelling and settlements prediction

This paper describes the application of an EVP one-dimensional constitutive soil model to analyse regional subsidence. The model was modified assuming that the boundary conditions vary over time, as well as considering the stratification of the subsoil. It adequately simulates the consolidation process caused by the exploitation of the aquifers that underlie the Mexico basin, and is validated at a site in the former Texcoco Lake, some 14 km north of Mexico City, where geotechnical and piezometric information is available together with records of past subsidence. Settlement predictions were carried out for different periods of time, assuming that pore pressure depletion rates at the permeable borders of the compressible strata remain constant over time.

1 Introduction

Regional subsidence affecting the lacustrine zone of the Mexico Basin is the result of effective stress increments brought about by the depletion of pore pressure due to extensive extraction of water from the aquifers that underlie the compressible clay strata. Pore pressures have been gradually decreasing and trends in the evolution of pore pressures in the sands and sandy silts have been identified from records obtained from a number of piezometric stations installed throughout the city. High resolution levellings in the former Texcoco Lake area have shown that during the last decades the northern portion of the site is settling at about 6 to 8 cm yr⁻¹ whilst the settlement rates in its southwest sector reach 22 to 22 cm yr⁻¹, average settlement rate has been estimated to be 13.2 cm yr⁻¹ (II-UNAM, 2016). Water pumped from the aquifers provides about two-thirds of the city's supply and it is highly unlikely that pumping be stopped or even reduced in the future, given the trends of urban expansion observed over the last decades (Ovando Shelley et al., 2007).

Surface settlements are the result of changes that occur within the soil mass due to the consolidation process triggered by the depletion of pore pressures, as mentioned previously. These changes result from effective stress increments that have been taking place slowly and gradually for over 150 years in an on-going process of regional consolidation (Ovando Shelley et al., 2007; Ovando et al., 2013).

2 Simulation of excess pore pressure and settlements evolution

The evolution of excess pore pressure and settlements due to regional consolidation are studied using an EVP one-dimensional model proposed by Yin and Graham (1996) which overcomes some of the limitations of Terzaghi's theory. In this model, the total strain is the sum of the elastic and viscoplastic components. Hence, the soil is an elastic and viscoplastic material and consolidates according to the following equations:

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }{c}_{\mathrm{ve}}{\displaystyle \frac{{\partial }^{\mathrm{2}}u}{\partial z^{\mathrm{2}}}}={\displaystyle \frac{\partial u}{\partial t}}-{\displaystyle \frac{\mathrm{1}}{m_{\mathrm{ve}}}}g\left(u,{\mathit{\epsilon }}_z\right)\text{(2)}& {\displaystyle }& {\displaystyle \frac{\partial {\mathit{\epsilon }}_z}{\partial t}}=-m_{\mathrm{ve}}{\displaystyle \frac{\partial u}{\partial t}}+g\left(u,{\mathit{\epsilon }}_z\right)\text{(3)}& {\displaystyle }& {\displaystyle }g\left(u,{\mathit{\epsilon }}_z\right)={\displaystyle \frac{\mathit{\psi }/v_{\mathrm{0}}}{t_{\mathrm{0}}}}\exp \left[-\left({\mathit{\epsilon }}_z-{\mathit{\epsilon }}_{z_{\mathrm{0}}}^{\mathrm{vp}}\right){\displaystyle \frac{v_{\mathrm{0}}}{\mathit{\psi }}}\right]{\left({\displaystyle \frac{{\mathit{\sigma }}_z^{\prime }}{{\mathit{\sigma }}_{z_{\mathrm{0}}}^{\prime }}}\right)}^{\frac{\mathit{\lambda }}{\mathit{\psi }}}\end{array}$$

where, $c_{\mathrm{ve}}=k/\left(m_{\mathrm{ve}}{\mathit{\gamma }}_w\right)$ is the consolidation coefficient associated to elastic deformations, $m_{\mathrm{ve}}=(\mathit{\kappa }/v_{\mathrm{0}})/\left({\mathit{\sigma }}_z^{\prime }\right)$ is the volumetric compressibility modulus along the elastic portion of the one-dimensional stress–strain curve, $v_{\mathrm{0}}=(\mathrm{1}+e_{\mathrm{0}})$ represent the specific volume, k is the soil permeability and t₀ is a reference time, set to 24 h or t_EOP (the time at which primary consolidation takes place). Further, the ratios κ∕v₀ and λ∕v₀, are defined as the slopes of the swelling and normally compression lines in a ε_v vs. lnp^′ plane. In the EVP model, κ∕v₀ and λ∕v₀, represent the instant time line and the reference time line, respectively. Finally, ψ∕v₀ is the slope of the compression delay line.

Equations (1) to (3) comprise a system of non-linear differential equations that can be solved through the finite difference method using the implicit Crank–Nicholson scheme. This system must be solved under boundary conditions that simulate the evolution of the pore pressures in the permeable materials that confine the compressible clay strata. Initial conditions represent the distribution of pore pressures at a given site at the beginning of the study period. To specify the boundary conditions, it is necessary to know pore pressure depletion rates in the upper and lower ends of the confined stratum. Once these rates are determined, the functions defining the boundary conditions of the clay stratum can be specified from simple linear relationships such as:

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& u(\mathrm{0},t)=u(\mathrm{0},\mathrm{0})-V_s\times t\\ & u(\mathrm{2}H,t)=u(\mathrm{2}H,\mathrm{0})-V_i\times t\end{array}\end{array}$$

where, V_s and V_i are the pore pressure depletion rates at the upper and lower boundaries of the clayey stratum respectively and 2H is thickness of the clay stratum and t is the time. Conceivably pore pressure depletion rates can in general be expressed as polynomial functions of time.

To describe more realistically the consolidation process, soil permeability is assumed to vary with depth and time, considering the sketch shown in Fig. 1 (Nash and Ryde, 2001), where the value of the permeability at point X located at the interface between strata p and q is expressed as the ratio between their permeabilities as indicated below:

$$\begin{array}{}\text{(5)}& k_x^{\prime }={\displaystyle \frac{\mathrm{2}{k}_p{k}_q}{k_p+k_q}}.\end{array}$$

Figure 1 Finite difference scheme.

2.1 Validation of the EVP model

To validate the EVP model, the evolution of pore pressure depletion and settlement rates were simulated at site 7Pi located inside the area of the former Texcoco Lake as shown in Fig. 2. This site offered a wealth of historical, geotechnical, and piezometric data that resulted from three exploration surveys carried out between 2001 and 2014 as well as from records of regional subsidence rates obtained from direct high precision topographic surveys at 90 points during October 2013 and April 2014 (II-UNAM, 2016).

Figure 2 Regional subsidence rates and location of site 7Pi (II-UNAM, 2016).

As indicated in Fig. 3, the stratigraphy of site 7Pi consists of two compressible strata: the upper clay formation (FAS) and the lower clay formation (FAI), both interspersed with lenses of sand and silt. The superficial crust and the hard layer represent the upper and lower borders of the FAS, respectively, while the hard layer is constituted by a sequence of hard and soft lenses, and the deep deposits represent the lower border of the FAI. The deep deposits were considered as an incompressible layer. For numerical modelling, the compressible soil deposit were divided into 20 sub-layers, having the geotechnical properties listed in Table 1.

Figure 3 Geotechnical and piezometric conditions.

Table 1 Geotechnical properties at the 7Pi site.

OCR = 1.3, t₀=1 d.

2.1.1 Definition of the initial and boundary conditions

The evolution of the pore pressure was modeled to simulate the piezometric behavior of the site 7Pi. The initial pore pressure condition was fixed according to piezometric data collected at station EPA-1 during 2001 (see Fig. 3). The results obtained after five years were compared with data collected in 2013 from station 6-PZA, located near EPA-1.

Figure 4 Measured and simulated pore pressures.

On the other hand, pore pressure depletion rates V_s and V_i in the upper and lower boundaries of the compressible clay series were calculated from the piezometric records of the EPA-1 and EPA-3 stations from 2001 to 2013, and are presented in Table 2.

Table 2 Geotechnical properties at the site 7Pi.

2.1.2 Evolution of the piezometric conditions

Figure 4 shows the evolution of the pore pressure distribution calculated using the EVP model for the period 2001 and 2013 and compared with the data recorded at the EPA-1 and EPA-3 piezometric stations during the same period. As seen in this figure, there is quite a good agreement between the measured data and the pore pressures calculated with the EVP model.

Our results also indicated that settlements rates calculated at the top of the superficial crust closely matched the settlement rate measured at site 7Pi, 12 cm yr⁻¹, as indicated in Fig. 1.

2.2 Evolution of pore pressure and settlements

Having calibrated the model, a complementary analysis was performed to estimate the evolution of pore pressures and surface settlements at the site between 2013 and 2061.

As indicated in Fig. 5, the total settlement, estimated between 2013 and 2061 was 4.59 m. Settlements in the FAS will be 56 % of the total settlement; the remaining 44 % is associated to the compression of the FAI.

Figure 5 (a) Pore pressure evolution and (b) accumulated settlement, site 7Pi.

3 Conclusions

The elasto-viscoplastic model EVP can be readily applied to model soil consolidation as shown in this paper. Furthermore, the model can be adapted to estimate the evolution of pore pressures and settlements in soil strata subjected to a water extraction process by pumping.

Predicted settlements due to the regional subsidence at site 7Pi are significant and should be accounted for in the design of the structures in and around it.