The city of Algiers (Algeria) is a highly seismic area, and therefore, soil liquefaction poses a major concern for structures resting on sandy soil. A campaign of 62 static penetration tests or cone penetrometer tests (CPT) was carried out on a site located in the commune of Dar El Beïda in Algiers. The soil Liquefaction Potential Index (LPI) values were assessed, for each borehole, based on the simplified procedure of Seed and Idriss. On the other hand, the geographic information system and geostatistical analysis were used to quantify the risk of soil liquefaction at the studied site. It is worth mentioning that the (LPI) was taken as a regionalized variable. In addition, the experimental variogram was modeled on the basis of a spherical model. Also, the interpolation of the LPI values in the unsampled locations was performed by the Kriging technique using both isotropic and anisotropic models. Kriging standard deviation maps were produced for both cases. The cross-validation showed that the anisotropic model exhibited a better fit for the interpolation of the values of the soil liquefaction potential. The results obtained indicated that a significant part of the soil is liable to liquefy, in particular in the northwestern region of the study area. The findings suggest that there is a proportional relationship between the risk of liquefaction and the increase or decrease in seismic acceleration.

#### Keywords

- Potential liquefaction
- geostatistics
- geographic information systems
- variogram

Soil liquefaction is a phenomenon of instability or loss of strength that can take place on a generally saturated or partially saturated granular medium. This event is characterized by a higher pore pressure linked to the contracting behavior of soil during the application of a rapid loading (earthquake, shocks, tidal wave, etc.). The occurrence of liquefaction in soils is often assessed using the originally simplified method proposed by Seed and Idriss (1971).^{[27]} This method is based on the Standard Penetration Test (SPT), Cone Penetration Test (CPT), Marchetti Dilatometer Test (DMT), Shear Wave Velocity Measurement Test, and other laboratory tests, such as the cyclic triaxial test. These methods are based on the calculation of the safety factor. The in-situ sand liquefaction behavior depends on several characteristics, namely the grain size, water table, average density of soil, layer thickness of loose sand, as well as the history of static and dynamic stresses.^{[1]}

Using the probabilistic approach (reliability methods) for the calculation of the risk of soil liquefaction made it possible to assess the effect of variability and uncertainty existing for each parameter of the CSR (cyclic stress ratio) and CRR (cyclic resistance ratio) functions on the estimation of liquefaction risk.^{[13,14,2]} However, the geostatistical approach makes it possible to analyze the spatial structure of soil Liquefaction Potential Index (LPI) for the zoning of the risk potential in a site.^{[21,22,29]}

Geostatistics, which was invented in 1962 by Georges Matheron, an engineer from the Corps des Mines in France, represents the set of probabilistic methods intended for the study of regionalized phenomena.^{[18]} They have been successfully applied to issues related to mining geology, hydrogeology, environment, and so on. Geostatistical analysis, which takes into account two important issues, differs from the rest of the classical interpolation methods. The first issue lies in identifying the spatial structure of the variable studied using the variogram, which is considered a basic tool of geostatistics; the second one concerns the use of this spatial structure with known measured values, for the optimal prediction at unmeasured points, using the Kriging technique. Another advantage of the geostatistical analysis, over traditional interpolation methods, is the production of standard deviation maps of the predicted values. This allows engineers to identify areas of uncertainty.^{[3]} Several studies on the use of geostatistical techniques (interpolation or simulation) in geotechnical applications have previously been carried out.^{[3,15,21]}

On the other hand, it is interesting to recall that the Geographic Information System (GIS) was developed by the Environmental Systems Research Institute (ESRI). It is a tool that is primarily intended to manage, visualize, map, query, and analyze all data possessing a spatial component.^{[10]} The geostatistical analysis extension of the GIS software package offers statistical tools for analyzing and representing continuous data on a map in order to generate surfaces. Exploratory spatial data analysis tools in the GIS software offer different perspectives on the data, such as their distribution, overall trends, and level of spatial autocorrelation.^{[8]}

The present study aims to calculate the soil Liquefaction Potential Index (LPI), which represents a measure of the total risk of soil liquefaction over a strip 20 m deep. For this, a series of 62 static penetration tests (Cone Penetration Tests) were carried out in the study area. In addition, a geostatistical analysis was performed to study the spatial variation of the liquefaction potential index and to map the areas at risk of liquefaction using the Geographic Information System (GIS). The soil susceptibility to liquefaction was checked for a magnitude of 6.8 (magnitude of Boumerdes earthquake in 2003).^{[2]} This value was used in the calculations, while taking as the Peak Ground Acceleration (PGA) value on the surface as a_{max} = 0.30 g (value provided by the Algerian Parasismic Regulations).^{[26]} The possible acceleration of the ground taken into account during earthquakes is between 0.2 g and 0.25 g

The site under study is located in the Commune of Dar El Beïda, 16 km east of Algiers, in the Great Plain of Mitidja that stretches between the mountains of the Tellian Atlas and the Mediterranean Sea. Figure 1 gives an overview of the location on the studied site.

The data used in this study originate from surveys carried out with a static penetrometer of Gouda type. The region under study is located in an area of high seismicity; it is classified in the seismic category Zone III, according to the Algerian earthquake regulation in force.^{[26]} The behavior of the water table under the site study ground exhibited significant variations of up to 5 m, depending on the season. It was therefore advisable to consider that the water table was flush.

The cyclic stress ratio (CSR) is calculated according to the following formula^{[27]}:
_{max}_{v0}^{′}_{v0}_{d}

To calculate the factor r_{d}, [31] proposed the following relationship:

With regard to the Magnitude Scaling Factor (MSF), the lower bound equation, which was suggested by [31], was used as follows:

The cyclic resistance ratio (CRR) is estimated from CPT data following the procedure in [23], [24] as:

The resistance to penetration of a standard cone may be calculated using the following equation:
_{c} is a correction factor that depends on I_{c} and defined by the equation of [23] as:
_{c1N} (CPT) was determined using the following formula:
_{at} is the reference pressure (100 kpa), q_{c} is the peak resistance, and σ_{v0} and σ^{′}_{v0} are the total and effective vertical stresses, respectively.

The soil behavior type index I_{c} is defined by the expression proposed by [23]:

The local factor of safety (FS) coefficient for the soil liquefaction is such that^{[12]}:
^{[12]} It is worth mentioning that the soil becomes more resistant to liquefaction as the safety factor increases.

The first step consists of determining the safety factor (FS), and the second one is to find the liquefaction potential index (LPI), which is a measure of the total risk of liquefaction of soil to a depth of 20 m. The liquefaction potential index is defined as follows^{[11,12]}:
_{i}(z) = 10 - 0.5z, and z denotes the depth in meters; Hi is the thickness of the discretized soil layers, _{Li} may be defined as [17]:

Table 1 presents an example about the evaluation of the liquefaction potential index for a typical borehole (BH34). Figure 2 shows some illustrative graphs to show the major steps involved in calculating the liquefaction potential at one CPT sounding location. The liquefaction potential index was calculated for each CPT profile for an earthquake scenario with Mw = 6.8 and having a peak horizontal ground acceleration a_{max} = 0.3 g (Table 2). The CPT positions were defined by the Universal Transverse Mercator (UTM 31 Zone 31N) coordinate system.

Computation of LPI for peak horizontal ground acceleration of 0.3 g corresponding to M_{w} = 6.8.

^{3}) |
_{V0} (Kpa) |
^{′}_{V0} (Kpa) |
_{c} (Kpa) |
_{s} (Kpa) |
_{d} |
_{s} |
||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|

2.50 | 19.58 | 48.95 | 23.95 | 2550.00 | 125.00 | 0.59 | 0.98 | 0.30 | 1.95 | 1.25 | 2.50 | 9.38 | 0.00 | 0.00 |

3.50 | 19.58 | 68.53 | 33.53 | 1140.00 | 41.00 | 0.19 | 0.98 | 0.30 | 0.62 | 3.00 | 1.00 | 8.50 | 0.38 | 3.21 |

8.50 | 19.58 | 166.43 | 81.43 | 2430.00 | 130.00 | 0.33 | 0.93 | 0.29 | 1.15 | 6.00 | 5.00 | 7.00 | 0.00 | 0.00 |

10.50 | 19.58 | 205.59 | 100.59 | 1450.00 | 51.00 | 0.13 | 0.89 | 0.28 | 0.48 | 9.50 | 2.00 | 5.25 | 0.52 | 5.43 |

14.50 | 19.58 | 283.91 | 138.91 | 2860.00 | 112.00 | 0.18 | 0.78 | 0.24 | 0.73 | 12.50 | 4.00 | 3.75 | 0.27 | 4.02 |

Calculated values of liquefaction potential index (LPI).

BH1 | 519039.1 | 4062223 | 5.57 | BH32 | 517267.3 | 4060744 | 24.11 |

BH2 | 519247.8 | 4061957 | 0.00 | BH33 | 517582.2 | 4061421 | 16.97 |

BH3 | 519471.7 | 4061669 | 0.73 | BH34 | 516748.6 | 4060166 | 12.66 |

BH4 | 519697.4 | 4061414 | 0.00 | BH35 | 516205.9 | 4060353 | 26.85 |

BH5 | 519931.9 | 4061204 | 0.00 | BH36 | 516494 | 4060620 | 26.54 |

BH6 | 519925.8 | 4060761 | 0.00 | BH37 | 515785.8 | 4060475 | 21.12 |

BH7 | 519655.4 | 4061015 | 0.00 | BH38 | 515446.7 | 4061284 | 23.90 |

BH8 | 519438.6 | 4061314 | 0.00 | BH39 | 515676.1 | 4061850 | 30.48 |

BH9 | 519187.1 | 4061546 | 2.34 | BH40 | 515196.3 | 4061971 | 32.72 |

BH10 | 518939.2 | 4061757 | 5.73 | BH41 | 517075.7 | 4059634 | 27.74 |

BH11 | 518750.8 | 4062111 | 8.60 | BH42 | 516106.9 | 4061729 | 29.61 |

BH12 | 519674.6 | 4060860 | 8.75 | BH43 | 515283.1 | 4059841 | 13.44 |

BH13 | 519369.7 | 4060970 | 3.80 | BH44 | 515430 | 4060130 | 12.45 |

BH14 | 519124.5 | 4061169 | 0.00 | BH45 | 515829.7 | 4059942 | 23.93 |

BH15 | 518443.2 | 4061889 | 9.97 | BH46 | 518340.9 | 4061656 | 19.64 |

BH16 | 519199.3 | 4060493 | 8.61 | BH47 | 517981.9 | 4061644 | 6.18 |

BH17 | 518900 | 4060902 | 5.38 | BH48 | 517343.3 | 4058514 | 0.93 |

BH18 | 518591 | 4061257 | 16.25 | BH49 | 517268.7 | 4058736 | 0.00 |

BH19 | 518848 | 4060181 | 3.69 | BH50 | 517249.5 | 4058958 | 0.00 |

BH20 | 518474.9 | 4060857 | 19.59 | BH51 | 517536.8 | 4059125 | 0.21 |

BH21 | 518374.6 | 4060557 | 24.41 | BH52 | 517465.8 | 4058892 | 0.00 |

BH22 | 518288 | 4060102 | 26.39 | BH53 | 517511 | 4058648 | 3.24 |

BH23 | 518125.3 | 4059736 | 13.89 | BH54 | 517506.3 | 4058348 | 0.00 |

BH24 | 517678.6 | 4059724 | 24.30 | BH55 | 517628.5 | 4058870 | 1.83 |

BH25 | 517901.1 | 4060124 | 28.98 | BH56 | 517669.1 | 4058681 | 5.21 |

BH26 | 518152.4 | 4060413 | 25.75 | BH57 | 517711.5 | 4058493 | 1.61 |

BH27 | 518130.8 | 4060901 | 24.28 | BH58 | 515780.4 | 4058544 | 6.00 |

BH28 | 517660.5 | 4060689 | 13.15 | BH59 | 515724.8 | 4059143 | 3.29 |

BH29 | 517854.7 | 4060945 | 11.10 | BH60 | 515743.1 | 4059410 | 4.85 |

BH30 | 517444.4 | 4059335 | 0.00 | BH61 | 515869.5 | 4059687 | 11.92 |

BH31 | 517097.9 | 4061032 | 20.12 | BH62 | 515898.42 | 4059487.6 | 11.66 |

The variogram is a tool that allows measuring the dissimilarity between values based on their separations. It describes the spatial continuity of the regionalized variable.^{[8]} The theoretical variogram is defined as:

The experimental variogram, proposed by [18], is defined as follows^{[8]}:
_{i}_{i} + h)_{i})_{i}_{i} + h)_{i} + h)

Ordinary kriging is an interpolation procedure that provides, from the available data, an optimal, linear and unbiased estimate of the property under study and whose estimation error is minimized.^{[6]} The interpolated value at the point x_{0}, denoted _{0}), is given by:
_{i}^{[19, 7, 29]}:
_{i}_{j}_{i})_{j})_{i}_{j}; _{0}_{i}, x_{j})_{i}_{j}, x_{0})_{i}_{0}

Kriging also makes it possible to evaluate the uncertainty associated with the estimator in the form of the kriging variance, which is given as:

ArcGIS Geostatistical Analyst provides statistical tools that allow analyzing and mapping continuous data.^{[9]} Figure (3) illustrates the location of the data and the spatial distribution of the liquefaction index on the site under study. It is easy to notice that the distribution of boreholes is not regular and does not cover the entire study area. The statistical distribution of the Liquefaction Potential Index (LPI) is illustrated by the histogram in Table 3. The Kurtosis coefficient (1.79) is greater than 0, which implies that the distribution is more flattened than a normal distribution. The Skewness coefficient (0.47) indicates that this distribution is not symmetrical; it is an asymmetric distribution, with a spread towards the low values. Therefore, a logarithmic transformation could certainly help to make the distribution of the data more symmetrical. However, this transformation is not recommended by several authors.^{[4,20,25,30]}

Basic statistics of liquefaction potential index.

62 | 11.459 | 8.6759 | 10.456 | 0 | 32.718 | 0.472 | 1.793 |

The trend analysis tool provides a three-dimensional perspective of the data.^{[9]} The green line indicates the decreased value in the east-west direction, and the blue line shows a bell curve in the north-south direction. Liquefaction potential index values are higher in the northwest but there is not a sufficiently big trend in our data (Figure 4).

The variogram is used to characterize the spatial continuity of the liquefaction potential index (LPI).^{[21,22]} Figure 5 shows the experimental variogram, calculated for multiple distances, with a step of 300 meters, while assuming that this is an isotropic soil. It is worth noting that the fitting of several models, available in the ^{2} = 0.943 and minimum residual sum of squares (Table 4). The spherical model adopted in this analysis is defined by the following formula:
_{0} is the nugget, and C is the sill.

Properties of the fitted semivariograms of the liquefaction potential index.

_{0} |
_{0} |
_{0} |
^{2} |
|||
---|---|---|---|---|---|---|

Spherical | 0 | 120 | 1800 | 1 | 1728 | 0.943 |

Exponential | 0 | 135 | 1060 | 1 | 2633 | 0.885 |

Linear | 38 | 146 | 3440 | 0.735 | 7272 | 0.629 |

Gaussian | 5 | 125 | 1200 | 0.96 | 2205 | 0.900 |

Note: R^{2} is the coefficient of determination; RSS is the residual sum of squares, and c/c+c_{0} is nugget-sill ratio.

If the variogram actually depends only on the norm of h, it is said to be isotropic. On the other hand, if it depends on the h direction also, then it is anisotropic. Figure 6 shows directional variograms along four directions, that is, 0°, 45°, 90°, and 135°, with respect to the azimuthal direction. It is easy to see that the variograms have almost the same plateau, and the semivariogram in the 90° (east-west) direction grows more slowly than in the other directions. Therefore, one may say that there is a great spatial continuity of the liquefaction potential index in that direction, with a range of around 2600 meters. It can therefore be deduced that the geotechnical properties of the soil are more homogeneous in the east-west direction in comparison with the other directions. Semivariograms in the 0° (north-south) and 45° directions grow faster. Consequently, it can be said that there is a small spatial continuity of the liquefaction potential index along these 2 directions, with a range of about 1200 meters. It can hence be deduced that the geotechnical properties of the soil are rather heterogeneous.

Figure 7 depicts the theoretical model of a two-dimensional anisotropy, with a spherical model with a range of 2700 m, following the maximum direction of continuity of 90°, with a range of 1150 m. In the other minimum direction of 0° continuity, the plateau remains constant and equal to 120.

Cross-validation consists of temporarily eliminating a point from the data set, then estimating its value by kriging using the remaining data and the variogram model that is to be adjusted.^{[5]} The clouds of correlations between the measured (LPI) values (on the abscissa) and the interpolated (LPI) values (on the ordinate) were compared for the two isotropic and anisotropic approaches. The results obtained showed a good value of the coefficient of determination R^{2}, while taking into account the anisotropy of the variogram depicted in Figure (8).

Ordinary kriging mapping of the soil liquefaction potential index over the entire study site was carried out. Maps for the ordinary kriging interpolation results as well as the corresponding kriging standard deviation maps for the two cases, isotropic and anisotropic, were presented. Figures 9a and 10a depict the spatial variations of the potential liquefaction index (LPI) of the soil in the study area. It was noted a large area was susceptible to liquefaction, particularly in the northwest of the region under study. Similarly, Figures 9b and 10b present the spatial

The coupling of GIS geographic information system techniques and geostatistical analysis made it possible to carry out a spatial modeling of the risk of soil liquefaction over the entire study site, using digital maps for the estimate of the probable spatial distribution of the liquefaction potential index (LPI) of soil. Moreover, the kriging standard deviation maps of the (LPI) predicted values are quite important tools for geostatistical analysis. These maps can be used to identify areas of uncertainty. The results obtained showed that the higher the standard deviation, the greater the uncertainty. The presence of a large zone susceptible to liquefaction, in particular in the northwestern part of the study area, was confirmed by the results obtained. Note that the kriging standard deviation did not depend on the measured (LPI) values, but depended only on the variogram model and the distributions of the measurement points. Furthermore, the findings indicated that there is a proportional relationship between the liquefaction risk and the increase or decrease in seismic acceleration.

#### Computation of LPI for peak horizontal ground acceleration of 0.3 g corresponding to Mw = 6.8.

^{3}) |
_{V0} (Kpa) |
^{′}_{V0} (Kpa) |
_{c} (Kpa) |
_{s} (Kpa) |
_{d} |
_{s} |
||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|

2.50 | 19.58 | 48.95 | 23.95 | 2550.00 | 125.00 | 0.59 | 0.98 | 0.30 | 1.95 | 1.25 | 2.50 | 9.38 | 0.00 | 0.00 |

3.50 | 19.58 | 68.53 | 33.53 | 1140.00 | 41.00 | 0.19 | 0.98 | 0.30 | 0.62 | 3.00 | 1.00 | 8.50 | 0.38 | 3.21 |

8.50 | 19.58 | 166.43 | 81.43 | 2430.00 | 130.00 | 0.33 | 0.93 | 0.29 | 1.15 | 6.00 | 5.00 | 7.00 | 0.00 | 0.00 |

10.50 | 19.58 | 205.59 | 100.59 | 1450.00 | 51.00 | 0.13 | 0.89 | 0.28 | 0.48 | 9.50 | 2.00 | 5.25 | 0.52 | 5.43 |

14.50 | 19.58 | 283.91 | 138.91 | 2860.00 | 112.00 | 0.18 | 0.78 | 0.24 | 0.73 | 12.50 | 4.00 | 3.75 | 0.27 | 4.02 |

#### Basic statistics of liquefaction potential index.

62 | 11.459 | 8.6759 | 10.456 | 0 | 32.718 | 0.472 | 1.793 |

#### Properties of the fitted semivariograms of the liquefaction potential index.

_{0} |
_{0} |
_{0} |
^{2} |
|||
---|---|---|---|---|---|---|

Spherical | 0 | 120 | 1800 | 1 | 1728 | 0.943 |

Exponential | 0 | 135 | 1060 | 1 | 2633 | 0.885 |

Linear | 38 | 146 | 3440 | 0.735 | 7272 | 0.629 |

Gaussian | 5 | 125 | 1200 | 0.96 | 2205 | 0.900 |

#### Calculated values of liquefaction potential index (LPI).

BH1 | 519039.1 | 4062223 | 5.57 | BH32 | 517267.3 | 4060744 | 24.11 |

BH2 | 519247.8 | 4061957 | 0.00 | BH33 | 517582.2 | 4061421 | 16.97 |

BH3 | 519471.7 | 4061669 | 0.73 | BH34 | 516748.6 | 4060166 | 12.66 |

BH4 | 519697.4 | 4061414 | 0.00 | BH35 | 516205.9 | 4060353 | 26.85 |

BH5 | 519931.9 | 4061204 | 0.00 | BH36 | 516494 | 4060620 | 26.54 |

BH6 | 519925.8 | 4060761 | 0.00 | BH37 | 515785.8 | 4060475 | 21.12 |

BH7 | 519655.4 | 4061015 | 0.00 | BH38 | 515446.7 | 4061284 | 23.90 |

BH8 | 519438.6 | 4061314 | 0.00 | BH39 | 515676.1 | 4061850 | 30.48 |

BH9 | 519187.1 | 4061546 | 2.34 | BH40 | 515196.3 | 4061971 | 32.72 |

BH10 | 518939.2 | 4061757 | 5.73 | BH41 | 517075.7 | 4059634 | 27.74 |

BH11 | 518750.8 | 4062111 | 8.60 | BH42 | 516106.9 | 4061729 | 29.61 |

BH12 | 519674.6 | 4060860 | 8.75 | BH43 | 515283.1 | 4059841 | 13.44 |

BH13 | 519369.7 | 4060970 | 3.80 | BH44 | 515430 | 4060130 | 12.45 |

BH14 | 519124.5 | 4061169 | 0.00 | BH45 | 515829.7 | 4059942 | 23.93 |

BH15 | 518443.2 | 4061889 | 9.97 | BH46 | 518340.9 | 4061656 | 19.64 |

BH16 | 519199.3 | 4060493 | 8.61 | BH47 | 517981.9 | 4061644 | 6.18 |

BH17 | 518900 | 4060902 | 5.38 | BH48 | 517343.3 | 4058514 | 0.93 |

BH18 | 518591 | 4061257 | 16.25 | BH49 | 517268.7 | 4058736 | 0.00 |

BH19 | 518848 | 4060181 | 3.69 | BH50 | 517249.5 | 4058958 | 0.00 |

BH20 | 518474.9 | 4060857 | 19.59 | BH51 | 517536.8 | 4059125 | 0.21 |

BH21 | 518374.6 | 4060557 | 24.41 | BH52 | 517465.8 | 4058892 | 0.00 |

BH22 | 518288 | 4060102 | 26.39 | BH53 | 517511 | 4058648 | 3.24 |

BH23 | 518125.3 | 4059736 | 13.89 | BH54 | 517506.3 | 4058348 | 0.00 |

BH24 | 517678.6 | 4059724 | 24.30 | BH55 | 517628.5 | 4058870 | 1.83 |

BH25 | 517901.1 | 4060124 | 28.98 | BH56 | 517669.1 | 4058681 | 5.21 |

BH26 | 518152.4 | 4060413 | 25.75 | BH57 | 517711.5 | 4058493 | 1.61 |

BH27 | 518130.8 | 4060901 | 24.28 | BH58 | 515780.4 | 4058544 | 6.00 |

BH28 | 517660.5 | 4060689 | 13.15 | BH59 | 515724.8 | 4059143 | 3.29 |

BH29 | 517854.7 | 4060945 | 11.10 | BH60 | 515743.1 | 4059410 | 4.85 |

BH30 | 517444.4 | 4059335 | 0.00 | BH61 | 515869.5 | 4059687 | 11.92 |

BH31 | 517097.9 | 4061032 | 20.12 | BH62 | 515898.42 | 4059487.6 | 11.66 |

Seismic bearing capacity of shallow strip footing embedded in slope resting on two-layered soil FEM modelling of the static behaviour of reinforced concrete beams considering the nonlinear behaviour of the concrete Evaluation of Tunnel Contour Quality Index on the Basis of Terrestrial Laser Scanning Data Efficiency of ventilated facades in terms of airflow in the air gap Characteristic parameters of soil failure criteria for plane strain conditions – experimental and semi-theoretical study A comprehensive approach to the optimization of design solutions for dry anti-flood reservoir dams Laboratory tests and analysis of CIPP epoxy resin internal liners used in pipelines – part II: comparative analysis with the use of the FEM and engineering algorithms Usefulness of the CPTU method in evaluating shear modulus G changes in the subsoil_{0}