A Classification Tree for Modeling Ground Fractures from Subsidence

Abstract

This article presents a classification tree with predictive and prescriptive capabilities for the management of ground fractures in a crowded suburb in Tláhuac, a municipality of Mexico City. The tree is trained with observations of fractures parametrized with basic geotechnical and geological variables and specifications of the urban environment where they manifest. With the trained tree, the complexity of the scenario affected by the subsidence phenomenon is cleared because the relations between parameters can be viewed easily. With the trained tree, the influence of stratigraphic arrangements (geotechnical properties), geological conditions, size of roads and inhabited units, and location of water-pumping infrastructure in the appearance and exacerbation of cracks in soils is recognized, and this offers citizens and government administrators the possibility of anticipating damages and working on programs for improving structures and foundations, including relocation programs for communities at risk.

1. Introduction

In areas with rapid urbanization and demographic growth, the prolonged exploitation of groundwater is causing the land to sink with consequent dangerous impacts [1,2,3]. One of the responses of soil sinking is the collapse of the superficial layers, which affects fragile buildings and buried facilities that fail in blocs or by kilometers, resulting in huge economic losses and detriments to the quality of life of citizens. Superficial ground fracturing is one of the soil phenomena that Mexico City must face with extreme urgency and humanism, since it occurs in large portions of the deprived regions of the metropolis (Figure 1). Understanding what provokes a soil fracture is vital for the Mexican metropolis. Soil rupture, with opening and/or vertical displacement, is a geomechanical process that could be associated with the pumping of groundwater from unconsolidated sedimentary aquifer systems, and normally develops in arid or semiarid basins.

The generation and propagation of ruptures requires the development of failure conditions, that is, traction and/or shear stress not supported by shallow alluvial sediments. There are numerous analyses of stress fields that model the occurrence of ruptures, with locations in Mexico [4,5,6,7,8,9], China [10,11,12,13], and the United States [14,15,16,17,18,19,20,21,22,23] being among the most cited investigations. In CDMX, cracks associated with pumping have been studied from local and regional, timeless and evolutionary, and forensic and predictive perspectives [24,25,26,27,28,29,30,31,32], which have represented great contributions to the state-of-the-art research and the documentation of events. However, the discussion is still open in regard to defining the risk levels in specific properties (vulnerability and threats to the management of the effects) and recommending (preventing) certain structural solutions or accepting urban projects in localities that are affected by this phenomenon.

Tláhuac, a deprived mayoralty in southeastern Mexico City, is the subject of this research. In this study, area, cracks, fissures, and strong steps due to consolidation of soft layers on abrupt basement slopes, and differential subsidence in heterogeneous contacts, have been recognized [33]. The southern region of Mexico City has high levels of heterogeneity, and its typical stratigraphy is composed of superficial fillings (garbage and/or tuff strata) of poor or null capacity to resist load on extremely soft clays that were deposited on the steep slope of the perimeter of the basin (basement). These arrangements of materials manifest in complex behaviors whose prediction is a challenging open task. Identification of the processes is particularly difficult because the manifestations (cracks) are transfigured by interactions with inefficient foundations, inoperative communication routes, damaged buried infrastructure (leaks in water pipes), and the chaotic dynamics of the layers of debris (used to form a level when the urbanization of the area began) deforming under the stresses provoked by pumping wells.

When the best-known models that use groundwater flow equations and a complex geomechanical characterization are used to evaluate vertical deformations (associated with groundwater extraction) in the Tláhuac scenarios, the predictions fail because they entail difficulties and uncertainties, most of which are related to the parameterization of the environment. The transformation of specific information (often scarce) in a 3D space, with near and far borders, is a fragile point when predictions are faced with field observations. The disconnect between interpretations and reality has its most negative effect when ruptures occur without being able to alert government actors, developers, and the public themselves about the risk in their environment. The study of the behaviors and properties of soils has been analyzed with artificial intelligence and machine learning very profusely in recent years, showing the enormous potential of this technology to solve many geotechnical problems [34,35,36,37,38,39,40,41,42]; however, few investigations with sufficient solidity can be found on ground fracturing and subsidence. Using machine learning in the simulation of fracturing from subsidence makes it possible to generate concrete information about the aspects in the environment (natural and anthropic) that activate or exacerbate the breaks, while the processes that extend from the depths of the exploited aquifer to the superficial layers stayed immersed in the data.

This article presents a CART model (classification and regression tree) in classification mode that analyzes ground fracturing in a comprehensive way, eliminating subjective approaches and pointing out the spatial factors that most influence the adverse manifestations (cracks, fissures, and steps). The CART is trained with geotechnical and geological records that characterize the soils and rocks (composition and geometry), information on the buried infrastructure (pumping well location), and the arrangements of houses and streets in the most heavily fractured areas in Tláhuac. The tree makes it possible to define the areas that are most threatened by the phenomenon, the components (natural and anthropic) that exacerbate the manifestation, and the points on the surface that could potentially crack if their intrinsic characteristics and the water pumping patterns do not change. Its outstanding ability to predict breakpoints allows the use of machine learning to provide effective responses to city designers and government administrators for anticipating these kinds of events, work which is essential when the threatened areas have low resilience and a high level of poverty.

2. CART Basics

The CART model is a nonparametric learning technique that produces regression trees when the output is numerical and classification trees when the dependent variable is categorical. The general result of the CART algorithm [43] is a tree where the branches represent sets of decisions and each decision generates successive rules that continue the classification, also known as partition, thus, forming mutually exclusive homogeneous groups with respect to the variable discriminated. Trees are built using a recursive segmentation algorithm, ending once a stopping criterion is reached.

The method uses historical data to build the tree and, once completed, it can be used to classify new data. This means that $Y$ is a dependent variable and the $p$ predictor variables are $x_1,x_2,\dots ,x_n,$ where x is considered fixed and $Y$ d is a random variable. The statistical problem is to establish a relationship between $Y$ and x in such a way that it is possible to predict $Y$ based on the values of $x$. Mathematically it is required to study the conditional probability of the random variable $P [Y=y|x_1,x_2,\dots ,x_n]$ or a function of its probability, such as the conditional expectation [44].

2.1. Elements of a Tree

The tree in Figure 2 is made up of an initial node called the root node, formed from the attribute that creates the purest subset. It is divided in turn into two groups or decision nodes, to later apply the partition procedure separately to these groups. The divisions are selected in such a way that the purity of the decision nodes is greater than that of the root node. The goal is to partition the response into homogeneous groups while keeping the tree small enough. The recursive segmentation process continues until the tree is saturated, i.e., the subjects in the descendant nodes cannot be divided into an additional division since they do not fulfill the condition of the purity of these nodes being less than that of the node from which they come. These nodes that cannot continue to be divided are known as terminal nodes.

The methodology for developing regression and classification trees generally consists of three steps [45]. This section may be divided by subheadings. It should provide a concise and precise description of the experimental results, as well as their interpretation, and the experimental conclusions that can be drawn. This methodology is as follows:

• Construction of the tree;
• Pruning the tree;
• Selection of the optimal tree through a cross-validation procedure.

2.2. Node Impurity Function

This function is a measure that allows determining the quality of a node, expressed with i(t). The impurity measures that allow the analysis of different types of responses in the classification trees are available in [43,46].

The entropy index is defined as follows:

$$i\left(t\right)={{\displaystyle \sum }}_{j}p\left(j|t\right)lnp\left(j|t\right)$$

(1)

where

$p$ = relative frequency of class j in t/probability that an object is classified in a class

$j$ = class

$j$ = data set

The objective is to find the partition that maximizes $\Delta i\left(t\right)$ of the following equation:

$$\Delta i\left(t\right)=-{{\displaystyle \sum }}_{j=1}^{k}p\left(j|t\right)\ln p(j|t)$$

(2)

The Gini index defined as follows:

$$i\left(t\right)={{\displaystyle \sum }}_{i\ne j}p\left(j|t\right)p\left(i|t\right)$$

(3)

The following equation is used to find the partition that maximizes $\Delta i\left(t\right)$:

$$\Delta i\left(t\right)=-{{\displaystyle \sum }}_{j=1}^k{\left[p_j\left(t\right)\right]}^2$$

(4)

The most common algorithms for regression trees are based on the calculation of standard deviation, variance, and sum of squares.

2.3. Tree Pruning

In the first stages the tree obtained is over adjusted, so terminal nodes must be cut successively until the most suitable size is found. To solve this problem, an alternative is to look for a series of nested trees of decreasing size [47], each of which is the best of all trees of its size. These small trees are compared to determine the optimum. This comparison is based on a cost–complexity function, $R_{\alpha }\left(T\right)$.

For each T tree, this function [48] is defined as follows:

$$R_{\alpha }\left(T\right)=R\left(T\right)+\alpha \left|\tilde{T}\right|$$

(5)

where $R\left(T\right)$ is the average of the sum of squares between two nodes, which can be the total misclassification rate or the total residual sum of squares (depending on the type of tree); $\left|\tilde{T}\right|$ is the complexity of the tree, defined as the total number of nodes in the subtree, and $\alpha$ is the complexity parameter. The parameter α is a real number greater than or equal to zero; when $\alpha =0$ the result is the largest tree, and as it increases, the size of the tree decreases.

The function $R_{\alpha }\left(T\right)$ will always be minimized by the largest tree; therefore, better estimates of the error are needed. From the sequence of nested trees, it is necessary to select the one that is optimal, and to achieve this, it is not effective to use comparison or a complexity penalty [48]; therefore, it is necessary to accurately estimate the prediction error, and in general, this estimation is made using a cross-validation procedure.

The objective is to find the optimal ratio between the misclassification rate (quotient between the misclassified observations and the total number of observations) and the complexity of the tree. The cross-validation procedure can be implemented in two ways, depending on the amount of data [49]. With sufficient data, the sample is split; half or less of the data are removed, and then the sequence of trees is constructed using the data that remain. Posteriorly, the predictions of each tree are calculated, and the error of the estimations is set to finally select the tree with the smallest prediction error. If there were not enough data, the k-partition cross-validation tool was applied (k-fold cross-validation).

2.4. J48 Algorithm

Quinlan’s C4.5 algorithm actualizes J48 to create a trimmed C4.5 decision tree. Every aspect of the information is split into minor subsets based on a decision. Here, J48 looks at the standardized data gain that results in the split of the information by choosing an attribute. To summarize, the attribute’s extreme standardized data gain is utilized. The minor subsets are returned by the algorithm. The split strategies stop if a subset has a place with a similar class in all instances. Then, J48 develops a decision node utilizing the expected estimations of the class. The J48 decision tree can deal with particular characteristics, lost or missing attribute estimations of the data, and varying attribute costs. Here, accuracy can be expanded by pruning [50].

The decision trees produced by J48 can be utilized for classification. At every node of the tree, J48 chooses the attribute of the data that most effectively splits its arrangement of tests into subsets improved in one class or the other. The splitting criterion is the standardized information gain (in contrast to entropy). The attribute with the highest standardized information gain is used for making the decision. The J48 algorithm at that point recurs on the smaller sub lists. The J48 tree develops a decision node utilizing the expected estimations of the class. The J48 decision tree can deal with characteristics, lost or missing attribute estimations of the data, and varying attribute costs [51]. Some successful CART applications in engineering and geosciences are [52,53,54].

3. Modeling of Ground Cracking

The investigation focuses on one of the municipalities of Mexico City that has experienced the strongest and most dangerous ruptures in the last 50 years, namely the region of Tláhuac. In the 85 km² that make up this municipality, very particular areas are recognized where this phenomenon is forcefully manifested, particularly the Del Mar neighborhood (approx. 2 km²) (Figure 3), where cracks of tens of centimeters with openings of 4 m deep and huge steps approximately 120 cm high have been recorded. Agreeing with previous studies [7,55,56], ground fracturing is caused by the interaction of different factors, between the most important, geological preexisting discontinuities (because of variations in the depositional environment), slope of the basement underlying the compressible layers, thermal cycles, heterogeneity in compressibility and permeability, and the extensive exploitation of aquifers. This last aspect is considered the trigger factor of subsidence due to the corresponding vertical and horizontal tensile stresses on certain soil strata [57,58].

When the water level of the aquifer system drops, it induces a gradual compaction of the sedimentary filling, causing superficial sinking. In Tláhuac, this sinking is uneven, so superficial cracks and scarps are formed. There are intrinsic (natural) factors, such as material heterogeneity and basement characteristics, that control the trace, shape, and size (Figure 4) [22,59]. It is also recognized that anthropic settings, such as foundations, buried installations, loads (imposed by houses and buildings), and covert anomalies (e.g., pre-Hispanic structures) are also responsible for the extent and spatial definition of these cracks [60].

For training the CART model, the selected physical variables (natural and anthropic) are integrated as a spatial database, to construct the training matrix. Natural information (geological and geotechnical) and man-made conditions (depletion of groundwater levels, location of dwelling arrangement, and communication routes) were compiled and related to cracks, fissures, and steps in the studied area.

3.1. Exploration and Expression of Variables

From geotechnical and geophysical campaigns (49 geotechnical boreholes, aerial photogrammetry of the area, and 31 measurements of environmental seismic noise), important parameters were defined. This information was mapped in 2D/flat units (squares of 121 m²) georeferenced to the center of each unit. The square is selected for the ease of translating its information into pixels and voxels. Water content W%, Atterberg limits (liquid limit, LL; plastic limit LP; and plasticity index, IP), number of blows of SPT NSPT (average for the first 35 m of the column), materials classification according to the USCS (Unified Soil Classification System), cohesion/friction values (average in the first 35 m), and the presence of a semirigid layer embedded in the clay matrix are the basic geotechnical inputs. Additionally, data concerning soft deposit thickness (clay strata), amplification, and soil period were included. The depth and slope of the basement that underlies the clay formation was also included. Some examples of these maps are shown in Figure 5, Figure 6, Figure 7 and Figure 8. The included anthropic parameters are the housing level (number of floors and type of structure), the size of the streets (number of lanes), the position relative to the pumping wells, and the groundwater level. A summary of the inputs is listed in

Table 1. Input parameters description.

Parameter	Definition	Units
Geotechnical inputs
Water content W%	W%, also known as natural water content or natural moisture content, is the ratio of the weight of water to the weight of the solids in a mass of soil.	%
Liquid limit LL	LL is the water content where the soil starts to behave as a liquid.	%
Plastic limit PL	PL is the moisture content at which a fine-grained soil cannot be remolded without cracking.	%
Plasticity index IP	PI is a quantity which is the range of the water content over which the soil remains in plastic state.	%
NSPT	It is the number of blows required to penetrate 30 cm depth, known as the “standard penetration resistance SPT”, used to estimate some important geotechnical engineering properties of the soil.	Number of blows
Inclusion	Sufficiently thin semirigid layer	I–V, danger levels—function of stiffness, thickness, and depth
Superficial filling	Surface layer composed of low-quality fill (garbage, construction waste)	Meters
USCS	The world-accepted geotechnical soil classification.
Cohesion/friction values	Cohesion refers to shear strength under zero normal stress, or the intercept of a material’s failure envelope with shear stress axis in the shear stress–normal stress space. Friction angle is derived from the Mohr–Coulomb failure criterion, used to describe the friction shear resistance of soils (together with the normal effective stress).	kPa/°
Soft deposits thickness	Thickness of the soft clay layers; the materials susceptible to consolidation due to pumping in deep layers are here.	Meters
Amplification	It is the amplification factor used to account for site effects, i.e., how a site condition modifies the rock acceleration response spectra.	Dimensionless
Soil period	The period is the time interval required for one full cycle of a wave. It depends on the characteristics of the soil deposits.	Seconds
Depth of the basement	Measured from the surface, it refers to the depth at which the layer considered firm is found.	Meter
Slope of the basement	It refers to the slope of the layer considered firm.	Grades
Anthropic parameters
Housing level	It refers to the classification of houses, including their weight (related to build height) and structure solution (included foundation and superstructure).	Number of floors, from one to five and, the structuration as precarious, barely, and complete.
Size of the streets	It refers to the traffic load that goes from a minimum local flow to very heavy load motor transport circulation.	Number of lanes.
Pumping wells	It is expressed as the distance between the closest pumping-well(s) and the studied site.	X and Y, meters
Groundwater level	Referred to as the water table, this represents the top of the saturated zone; above the water table lies the unsaturated zone.	Meters

.

3.2. Tree Structure

The CART was trained using 3645 lines with 17 inputs and 1 output. The appropriate depth was determined by evaluating the tree on the held-out data set (700 lines) via k-cross validation. By resampling the data many times, splitting the data into training and validation folds, fitting trees with different sizes on the training folds, and looking at the classification accuracy on the validation test, the best tree depth (bias–variance trade-off) was found. The optimal tree (69 leaves and 137 terminal nodes) correctly classified 3338 instances (91.5775%), and 307 were incorrectly classified (8.4225%). Examples of the branches are shown in Figure 9. This tree predicts whether a site (square unit of 121 m²) cracks or not based on its intrinsic and extrinsic conditions, considering that consolidation (due to pumping) remains constant and acts on the same layers (those susceptible to consolidation) that drain toward sufficiently permeable strata.

Beyond its predictive power, this tree uncovers parametric relationships of paramount importance. The first partitioning of the data, the purest subset, is based on the characteristics (stiffness, thickness, and depth) of a stratum called inclusion. This clayey-sand layer is minor (its thickness varies between 50 and 400 cm) and ranges from a compact and fragile material to a soft, ductile soil mass. The depth of this stratum increases as it moves away from the steep gradient of the bedrock (underlying the clay layers), practically vanishing when the slope approaches zero (Figure 10). The behavior of sites with a superficial and fragile inclusion (labeled with ordinals from IV and V, with V being the more breakable) and of those with this layer in conditions of ductility and greater depth (labeled with ordinals from I to III, with I indicating increasing ductility) are clearly recognized by the CART. The next partition is based on the relative position to pumping wells, and the load imposed by houses and traffic is essential for the unit to crack.

Based on the results obtained with the optimal CART, the confusion matrix was constructed (

Table 2. Confusion matrix.

Total Population 3645	Actual Positive (AP) 309	Actual Negative (AN) 3336
Predicted Positive (PP) 166	True Positive (TP) 143	False Positive (FP) 23
Predicted Negative (PN) 3479	False Negative (FN) 166	True Negative (TN) 3313

). The total samples in the positive class are 309 and the number of samples in the negative class is 3336. The correct classifications are 143 for the positive class and 3313 for the negative class. Now, 166 samples that were expected to be of the positive class were classified as the negative class by the model (false negatives) while 23 samples were expected to be of the negative class but were classified as “positive” by the model (false positives).

With this confusion matrix, the accuracy of the tree was calculated using the following Equation (6):

$$Accuracy=\frac{TP+TN}{Total Population}=\frac{143+3313}{3645}=94.81\%$$

(6)

More than the 94% of samples were correctly classified out of all the samples present in the test set. To evaluate the worth of each attribute, the correlation (Pearson’s) between it and the output class was calculated (

Table 3. Correlation attribute level.

Parameter	Correlation Attribute Level
0.174244	Size of the streets
0.17366	Inclusion
0.14621	x
0.13615	Depth of the basement
0.07877	Water content %
0.02658	Groundwater level
0.02252	y
0.01159	Housing level
0.00996	Soft deposit thickness

). The nominal attributes were considered on a value-by-value basis by treating each value as an indicator, so, for a nominal attribute, the correlation was arrived at via a weighted average. The heterogeneity of the soil environment (represented with the inclusion) and the anthropic load of the streets are directly related to the occurrence of superficial cracks. For government administrators, it is essential to recognize that the distance in x that separates the inhabited sites from the pumping wells is very significant on the manifestation. Based on these findings, more threatened areas could be outlined, not just traces of breaks, but entire neighborhoods that will suffer the consequences of subsidence. The attention programs, and destination of the resources, can be more efficiently designed.

An illustration of the prediction process with CART is shown in Figure 11. On a section of the studied area, each patch (square unit) is filled with its corresponding input information to be able to follow the branches. When selecting a patch, for example, one with inclusion = V, the right branch of the tree must be followed to then question its position in space relative to the pumping wells (if x ≤ 13 or x > 13 and y ≤ 35 or y > 35). Because the unit is located at x = 11 and y = 32, then it is questioned about its housing level (>0 or ≤0 and >2 or ≤2). As the unit has a housing level = 3, the next node asks for the thickness of material with 10 < W% < 80%, and then its position in the final leaf categorizes it as “cracking not-activated”.

If, for example, the selected unit has an inclusion = V and its position in space is x <= 13 and y <= 35, the next suitable branch asks for the housing level (>2 or ≤2), then, due to its relative position to the pumping wells, it is categorized as “cracking activated”.

The application of CART to an expanded area is shown in Figure 12, and photographs of the evidence in the field are presented to point out the remarkable predictive capacity of the tree. It is important to mention that these cases are used for validation of the model, i.e., they are not contained in the training file.

4. Discussion

This article presents an alternative method to provide answers to citizens and governments that must build urban settlements and ensure their safety. The most sophisticated geomechanical models on fracturing have shown that their predictive capacity is very low due to subsidence, but what is most serious is that they do not allow for simple and unambiguous property-by-property evaluations. It is even practically impossible to include aspects of the city that exacerbate these breaks, such as the weight of the houses, heavy traffic, or the distance to pumping wells. The presented classification tree does not intend to explain the process of exploitation of groundwater but the connection between soils, urban scenarios, and fracturing.

This model could be built due to the enormous effort of the research team but also due to the government’s willingness to support field questionnaires and the execution of detailed geotechnical and geological tests. Until the 2017 earthquake, no campaign of this scale had been carried out in these lands. For this reason, the CART presented can be placed as one of the most detailed and best supported in this regard.

Streets, sidewalks, pipes, large buildings, and family homes are somehow affected by the “broken” condition of the environment, and idealization of the soils as a continuum (homogeneous layers) with regular rigid bases are strategies that have to be used when applying modeling tools that are not based on data. The CART presented here accepts inputs that describe heterogeneous stratigraphies (of very different materials, rigidities and continuities) that lie on a semi-firm and irregular, steep slope basement, but also, as in no other model, it accounts for houses and traffic loads, as well as the distance from pumping wells.

As this model was created using information contained in square pixel units (patches), expanding them to 3D (as voxels) would allow for refining the tree, so the manifestation is possible. The next step of this CART would be to generate metrics to qualify sites and to simulate subsidence evolution scenarios.