Articles | Volume 382
Proc. IAHS, 382, 521–524, 2020
Proc. IAHS, 382, 521–524, 2020

Pre-conference publication 22 Apr 2020

Pre-conference publication | 22 Apr 2020

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

Regional subsidence at the former Texcoco Lake: numerical modelling and settlements prediction
Efrain Ovando-Shelley, Alexandra Ossa-Lopez, and Renata Gonzalez Efrain Ovando-Shelley et al.
  • Instituto de Ingeniería, Universidad Nacional Autónoma de México, Mexico City, 04510, Mexico

Correspondence: Alexandra Ossa (


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−1 whilst the settlement rates in its southwest sector reach 22 to 22 cm yr−1, average settlement rate has been estimated to be 13.2 cm yr−1 (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:


where, cve=k/(mveγw) is the consolidation coefficient associated to elastic deformations, mve=(κ/v0)/(σz) is the volumetric compressibility modulus along the elastic portion of the one-dimensional stress–strain curve, v0=(1+e0) represent the specific volume, k is the soil permeability and t0 is a reference time, set to 24 h or tEOP (the time at which primary consolidation takes place). Further, the ratios κv0 and λv0, are defined as the slopes of the swelling and normally compression lines in a εv vs. lnp plane. In the EVP model, κv0 and λv0, represent the instant time line and the reference time line, respectively. Finally, ψv0 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:

(4) u ( 0 , t ) = u ( 0 , 0 ) - V s × t u ( 2 H , t ) = u ( 2 H , 0 ) - V i × t

where, Vs and Vi 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:

(5) k x = 2 k p k q k p + k q .

Figure 1Finite 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 2Regional 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 3Geotechnical and piezometric conditions.


Table 1Geotechnical properties at the 7Pi site.

OCR = 1.3, t0=1 d.

Download Print Version | Download XLSX

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 4Measured and simulated pore pressures.


On the other hand, pore pressure depletion rates Vs and Vi 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 2Geotechnical properties at the site 7Pi.

Download Print Version | Download XLSX

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−1, 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.

Data availability

Data will be made available on request.

Author contributions

EOS defined the main goal of the investigation, AOL developed the model code, RG performed the simulations and prepared the first version of the manuscript with contributions from all co-authors. EOS and AOL revised and edited the document.

Competing interests

The authors declare that they have no conflict of interest.

Special issue statement

This article is part of the special issue “TISOLS: the Tenth International Symposium On Land Subsidence – living with subsidence”. It is a result of the Tenth International Symposium on Land Subsidence, Delft, the Netherlands, 17–21 May 2021.


The authors acknowledge the valuable support provided by the Grupo Aeroportuario de la Ciudad de Mexico during the development of this work.


II-UNAM: Investigaciones y estudios especiales relacionados con aspectos geotécnicos del Nuevo Aeropuerto Internacional de la Ciudad de México (NAICM) en el vaso del exlago de Texcoco, Zona Federal, Final Report presented to Grupo Aeroportuario de la Ciudad de México, unpublished data, Diciembre 2016. 

Nash, D. and Ryde, S.: Modelling consolidation accelerated by vertical drains in soils subject to creep, Géotechnique, 51, 257–273, 2001. 

Ovando, E., Ossa, A., and Santoyo, E.: Effects of regional subsidence and earthquakes on arquitectural monuments in Mexico City, Boletín de la Sociedad Geológica Mexicana, 65, 157–167, 2013. 

Ovando Shelley, E., Ossa López, A., and Romo, M. P.: The sinking of Mexico City its effects on soils properties and seismic response, Soil Dynam. Earthq. Eng., 27, 333–343, 2007. 

Yin, J. H. and Graham, J.: Elastic visco-plastic modelling of one dimensional consolidation, Geotechnique, 46, 515–527, 1996. 

Short summary
This paper paper descibes the application of an EVP one-dimensional constitutive soil model to analyse regional subsidence. This model 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.