Application of finite difference method in determination of static shear stress influence on soil liquefaction

The aim of this study is to investigate the effect of static shear stress on the liquefaction potential of soil, and to propose a closed-form solution that takes this effect into account. Two-dimensional finite-difference-based numerical analyses, involving the use of various generic soil and earthquake combinations, were performed. The analyses include slopes with nonzero static shear stresses. It is concluded that the effect of the initial static shear stress on liquefaction resistance is highly dependent on the soil stiffness and on the initial on-site shear stress level, whereas nearly no effect is exerted by earthquake magnitude.


Introduction
Soil liquefaction potential is mainly estimated by Seed and Idriss's [1] "simplified procedure"; however, there are several other estimation methods that consider excess pore water pressure generation [2,3]. In the simplified procedure, the researchers proposed a cyclic resistance ratio factor (CRR) to define the resistance of soil to cyclic loadings, and the cyclic stress ratio (CSR) as the demand parameter. A detailed description of CSR calculation is given in Section 3. In the simplified procedure, the two parameters, CSR and CRR, are compared to each other. If CSR > CRR, it can be concluded that the soil under consideration can liquefy. If CRR is larger, then liquefaction is not probable. To make this comparison easier, charts were proposed in [1]. These charts were later on modified by various researchers [26,29]. However, the major flaw of this methodology is that CRR values and the corresponding resistance charts were developed for free-field site conditions only. These free-field sites are horizontal (level), and there is no shear stress at these sites before seismic excitation. Additionally, this methodology is valid for an atmospheric pressure of 1 atm, which is generally not the case for inclined surfaces. For example, shear stresses are present statically for inclined surfaces (slopes). Similarly, structures on level sites also create base shears, i.e., initial static shear stress. The effects of this initial (static) shear stress on liquefaction potential of the underlying soil is discussed in [4][5][6][7][8][9][10]. Before these studies, a still-valid correction factor (K α ) was proposed in [11], to take the effects of these initial (static) shear stresses into account. The CRR values for the stress state appearing at the site and the reference CRR values (found from "simplified procedure" and valid for zero initial static shear stresses) are multiplied by the correction factor K α . In addition to K α , there is another correction factor K s that corrects the CRR for the vertical effective stress. This is beyond the scope of this study but was included in the formulation for completeness. The CRR value is then calculated by Eq. (1) as defined in [26]: where CRR -cyclic resistance ratio at the current stress state CRR 1 -cyclic resistance ratio at the current stress state K s -effective confining stress correction K α -initial (static) shear stress correction.
The above procedure, which takes the (initial) static shear stresses and extra overburden into account in the liquefaction potential analysis, has been accepted by many researchers. However, there are some variations in the exact numerical values of K α as calculated by various researchers [12][13][14][15][16][17][18][22][23][24][25][26]. These research studies are discussed in detail in the following paragraph. The primary goal of this paper is to estimate the numerical value of K α by making use of simplicity of the numerical finite difference analysis. A closed-form solution, depending on the basic soil properties of relative density (D R ) and shear stress ratio (α), is proposed for this purpose.

Literature review
The undrained cycling response of soil is mainly affected by the (initial) static shear stresses. As mentioned previously, to consider this effect, Seed [11] introduced the K α correction factor, which accounts for these initial static shear stresses for inclined ground conditions. K α was found to be a shear stress factor before an earthquake (i.e., initial static shear stress) and relative density of soil. This (initial) static shear stress is measured in terms of a dimensionless parameter a, which can be found by dividing the (initial) static shear stress by vertical effective stress, as proposed in [11] and as presented in Eq. (2): ( This dimensionless parameter is the main feature for calculation of the initial stress correction factor (K α ), which is also highly dependent on relative density. When calculating K α , the most important parameters are the failure criteria, the confining stress, and the laboratory test procedure (such as the cyclic simple shear and anisotropically consolidated cyclic triaxial tests). As mentioned previously, these shear stresses that are determined at the site before cyclic loading (i.e., earthquakes) also affect sand resistance during cyclic loading. When initial static shear stresses exist, the K α values are higher than unity for soils with relative densities higher than 50 % (D R > 50 %). This results in an increase in the cyclic resistance ratio (CRR). However, for looser sands (D R < 45 %), the situation is just the opposite: the CRR values decrease when there is an initial static shear stress, i.e., K α < 1. Seed [11] studied dense soils (soils with D R > 50 %) and found that the (initial) static shear stresses increase the cyclic resistance, thus resulting in a K α correction factor higher than unity. Similarly, studies by [12][13][14][15][16][17][18] gave the same results. They also concluded that the cyclic resistance ratio of loose, sandy soil decreases with the initial static shear stresses, i.e., K α is less than 1.
In a more recent research [19], the authors conducted a number of undrained cyclic torsional shear tests using Toyoura sand. It was found that static shear stresses may have a positive or a negative effect on the liquefaction potential of soils, which is also dependent on the failure behavior and loading pattern. Similarly, based on undrained cyclic triaxial tests conducted at varying levels of initial deviator stresses, which enables us to see the combined effect of cyclic and static shear stresses on the undrained GRAĐEVINAR 71 (2019) 12, 1143-1151 Application of finite difference method in determination of static shear stress influence on soil liquefaction cyclic response of saturated sands, it is concluded in [20] that the relationship between the normalized residual pore pressure ratio and the normalized number of loading cycles is not affected by the amplitude of cyclic stress, but is significantly influenced by static shear stresses. A series of undrained cyclic triaxial tests was performed in [21] under various initial states on silty sands and it was concluded that the a concept can be extended to silty sands. This means that the static shear stress can be beneficial or detrimental to the liquefaction resistance depending on the initial state of the samples and the degree of stress reversal. K α correction factors proposed in [22,23] are also consistent with the previously mentioned studies for a range of relative densities, and the researchers verified the results of these studies. A practical guideline was proposed by [24] in which the effects of the relative density and confining stress as well as the relative state parameter index proposed by [25] (x R ) are used simultaneously. A summary of the recommendations by [26] for the K α correction factor is presented in Figure 1.
As apparent from this figure, K α varies between 0.3 and 1.7. Similar to the findings of the above mentioned researchers, the effect of initial static shear stresses on cyclic resistance of soils is positive for dense soils (D R > 50 %) and negative for looser sands (D R < 45 %). This results in a K α correction factor value greater than 1 for D R > 50 % and less than 1 for D R < 45 %.

Numerical simulations
A series of analyses were performed using the finite difference commercial software FLAC [27] to assess the effect of the initial shear stress on the soil liquefaction potential. The analyses were carried out in two dimensions (2D) as this is the universally accepted procedure for slope stability analysis. In the numerical assessment scheme, a 2D static and dynamics analysis of an inclined surface was carried out. To be consistent with recommendations of the numerical analysis software, and in order to prevent any numerical distortions, the sizes of the meshes were adjusted so that they were less than 10 % of the wavelength of the highest frequency component of the input wave. The boundary conditions were selected as "free field, " which enables a reduction in model size and, in turn, a reduction in computational effort. A typical mesh used in numerical simulations of slope is presented in Figure 2.

Figure 2. Two-dimensional finite difference model
In the numerical analysis, the material was selected as clean sand with varying shear wave velocities and, in turn, with varying stiffness values. The shear wave velocities of the soils ranged from 50 m/s to 400 m/s. The hysteretic degradation and damping curves [28] for PI = 0 were used to implement the nonlinear behavior of the model. The Mohr-Coulomb failure criteria were implemented in the model. The properties related to the soils used in the finite difference analyses are listed in Table 1. As can be seen in this  [31]. The shear and bulk modules are utilized in the numerical analysis. The shear modulus is calculated using the well-known equation G = r × V s 2 , where r is the unit mass of the soil. K is calculated using elastic relationships-K = E/(3 × (1 -2ν)) and G = E/(2 + (1 + ν))-where E is the elasticity (Young's) modulus, and ν is the Poisson's ratio. Point A in Figure 2 indicates the location where initial (static) shear stresses are greater than zero. Point F in the same figure represents the level surface response. That is the reason for such a skewed mesh, where a larger portion is located on the right side of the mesh. An attempt was made to obtain both the slope and level surface responses from the same analysis. Point F is the location at which there are no initial (static) shear stresses, i.e. it is the free-field point. Figure 3 shows shear stress values at the end of the static analysis. As seen in this figure, the shear stress value at Point F is nearly zero, whereas shear stresses are present if we cut the slope from Point A. To compare the effects of this initial shear stress on liquefaction potential, points on a vertical line below Points A and F are considered at every 1 m.

Figure 3. Shear stresses at the end of static analysis
Using the mesh in Figure 2 (10 to 20 m thick, 50 m wide), five generic homogeneous soil profiles with different relative densities ranging from 8 to 80 % (whose properties are listed in Table 1) were shaken by nine different earthquakes, from records: The details for the input motions used are presented in Table 2.

Evaluation of results
The results of the numerical analyses were evaluated using a simplified procedure [1]. As described in Section 1, in the simplified procedure, the CSR and CRR values are compared and, if CSR is larger than CRR, then the soil is said to be liquefiable.
The CSR values are calculated by dividing the average shear stress by the static effective stress. This is given in Eq. (3) below, as proposed in [1]: In this equation, CSR eq (z) is the CSR value at a depth of z, a max is the peak ground acceleration, g is the gravitational acceleration, g n is the unit weight of the soil layer, and s' ν (z) and r d (z) are the vertical effective stress and mass participation factor at depth z, respectively. In this equation, g n is multiplied by the depth to find the mass of the soil column (also known as the total vertical stress). Further, this mass is multiplied by a max to find the maximum shear stress that occurs during the earthquake. The factor 0.65 is the conversion of this maximum stress to the average shear stress. This average shear stress is then divided by the effective vertical stress, which is calculated by subtracting the pore water pressure from the total vertical stress at the depth considered. The mass participation factor r d (z) was calculated as proposed by NCEER [26] for noncritical projects, and details on this are given in their study. This factor is proposed to be a function of depth only, and it assumes the value of 1 at the surface. Then it decreases while going deeper, and assumes the value of 0.504 at a depth of 30 m. The next step is to calculate the CSR value at the free-field state, where the effects of the slope (initial shear stresses) can be neglected. The CSR value is also calculated using Eq. (3) for the free-field point (Point F in Figure 2), where there is no static shear stress. These two CSR values are then divided by each other to calculate the resulting K α correction factor, as shown in Eq. (4).
The CSR F and CSR A values are corrected for an overburden stress K s . This correction is made because, although they seem to be at the same level, the overburden stresses they encounter are different as their depths from the ground surface are not the same. In the calculation of K α , the steps below were followed for a specific example: V s = 50 m/s during the 1980 Loma Prieta earthquake (BRN). The data are presented in Table 3: -After application of strong ground motion data from the base of the model: -The maximum shear stress and maximum acceleration at point F are read (Column 12 in Table 3 right). -The maximum shear stress and maximum acceleration at point A are read (Column 3 in Table 3left).
-The values of α are calculated by dividing the shear stress values by the effective vertical stress at the points considered (Columns 7 and 16 in Table 3 for Points A and F, respectively).   When all cases, including all depths, earthquakes, and relative densities, are treated together, the result is not that meaningful. In the literature [11,12,24,26], the K α correction factors were calculated according to the stiffness, i.e., relative density, and were presented in graphs that were plotted for different relative densities. Thus, in this study, this factor is classified according to the relative densities. Figs. 5 through 9 in Section 4 show the NCEER recommendations as well as the results obtained from probabilistic analyses. These figures present the K α correction factors for relative densities (D R ) of 8 %, 15 %, 45 %, 60 %, and 80 %, respectively. In addition to the values obtained from the finite difference analyses, the recommendations of NCEER [26] are also embedded in these figures. NCEER [26] recommendations are mainly developed using the results of cyclic simple shear, cyclic ring torsional shear, and anisotropically consolidated cyclic triaxial tests that investigate the effects of static shear stresses on cyclic resistance. The results of the numerical analyses seem to fit properly with the NCEER [26] recommendations. These results could also be used to obtain a formulation to find the values of K α . For this purpose, a probabilistic analysis was performed to relate the initial shear stress ratio (α) and the stiffness of the soil (defined as the relative density D R of the soil profile in this study) with K α .

Probabilistic analyses
Although technological advancements and considerable progress have been made over the last several decades, dynamic analyses are still lengthy and complex. Therefore, dynamic analyses are not preferred at the preliminary design stage. This study aims to define a formulation for K α that can be used in preliminary analyses, as it is thought that the conduct of dynamic analyses is still difficult and time consuming, and that it requires a wide variety of data.

Selection of descriptive variables
The values of K α are predicted from the results of numerical analyses using the methodology explained in previous sections. Afterwards, a simple and user-friendly relationship that uses the results of these numerical analysis for predicting the K α value, is estimated. The maximum-likelihood methodology is selected as probabilistic tool for this purpose. For a proper probabilistic analysis, the descriptive parameters that dominate the liquefaction potential of the soils in the existence of initial static shear stresses must first be described. The important parameters defining K α are the initial shear stress ratio (α) and the stiffness of soil (defined as the relative density D R of the soil profile in this study). Figure 4 shows the variation of K α values with different peak ground acceleration (PGA) values of the earthquakes. As can be seen in this graph, the PGA of an earthquake does not mean anything from the point of view of K α . For example, when PGA = 0.2 g, the K α values range from 0 to 2.5 and are the same when PGA = 0.9 g. For this reason, PGA is not selected to be a variable for predicting K α in probabilistic analyses.

Relationship proposed for estimation of K α
After testing a series of alternatives as descriptive variables, D R and α are selected to be the main variables in the relationship proposed for the estimation of K α . The limit state function of the equation that best fits the results obtained from the numerical analysis is presented in Eq. (5). As is clear from this equation, the trends of numerical values of the initial shear stress correction factor (K α ) depend on the functional forms of D R and α. These functional forms are estimated separately using the maximumlikelihood methodology. The details of this methodology can be found in various resources such as [29].
There is a random model correction term (e) in the proposed model. The reasons for such an error term can be summarized as the insufficiency of the mathematical model proposed (it may not have the ideal form) and because some descriptive parameters may be missing that may also affect the soil liquefaction potential when (initial) static shear stresses exist. This random model correction term is assumed to have a normal distribution with a zero mean in an unbiased model. The standard deviation of this term is referred to as se and GRAĐEVINAR 71 (2019) 12, 1143-1151 Application of finite difference method in determination of static shear stress influence on soil liquefaction must be estimated. Therefore, we have a set of unknown parameters Θ which includes both the θ values (θ1-θ5) in Eq. (5) and se. The values of θ (θ1-θ5) are estimated so that the likelihood function in Eq. (5) assumes its maximum value. All values of D R and the corresponding a values obtained as a result of the numerical analyses are listed for this purpose. Then, the K α values, which were calculated as presented in Table 3 for each depth and scenario, are added to this list and are denoted as K α,measured . Then, using Eq. (5), K α values are calculated and denoted as K α,calculated . θ values for the maximum-likelihood function are then calculated. The θ values obtained as part of the maximum-likelihood methodology are presented in Table 4. Eq. (5) assumes the form presented in Eq. (6): The final form of Eq. (6) after inserting the model coefficients is: The proposed and calculated K α values are presented comparatively in Figs. 5-9. There are some deviations in these figures. However, it can be said that results of numerical analyses are in agreement with the results obtained from the formulation (Eq. 5). This equation is the best fit among many others when considering all the data. The differences in these figures are the results of possible missing parameters in the mathematical model. It is impossible to cover all the parameters in such a model, but this limited number of parameters is sufficient to obtain an appropriate value for K α . There may be other factors that affect K α . These are assumed to be included in the error term.

Conclusion
Given the confines of this study, 2D finite-difference-based slope stability analyses were performed. These analyses considered not only the static but also the dynamic stress states and performances. Using the results of these numerical analyses, a simplified formulation to obtain K α was proposed, as presented in Eq. (7). The results obtained from the numerical analyses were consistent with the values given in the literature. As stated by NCEER [26], the K α values increase with an increase in α for dense soils and decrease with an increase in α for loose soils. A maximum threshold value of 2.15 was proposed based on the results of numerical analyses. The results of probabilistic analyses showed that, from the liquefaction point of view, the most important parameter affecting the behavior of the slopes during cyclic loading is the stiffness of soil. Although it is a well-known fact that the liquefaction resistance increases with the soil stiffness at level sites, the existence of the initial shear stresses increases or decreases the liquefaction resistance of the soil when compared to level sites. If the soil is stiff (D R > 35 %), the presence of the initial shear stress is beneficial for liquefaction resistance. However, if the soil is loose (D R < 35 %), everything is opposite. Increasing the initial (static) shear stress, and thus increasing α, decreases the liquefaction resistance of soils. The results of the numerical analyses led to the same result. For soils with D R values higher than 35 %, the existence of the initial shear stress benefits the soil for liquefaction. The other important factor is the initial stress ratio (α). A change in α will clearly change the liquefaction resistance. Although not included in the formulation, the shape of the slope (height, obliquity, etc.) is one of the main factors. The shape of the slope is an indirect indicator of the shear stresses developed in the area. It is denoted by α and is also affected by the stiffness of the soil and potential failure surface. The formulation that estimates K α , which was obtained as a result of numerical analyses and is presented in this paper, should be used with caution. As the number of cases is limited and K α does not depend on the case history or laboratory test data, the formula can be used for comparison purposes in a preliminary analysis. For a final accurate result, a dynamic numerical model should be implemented, and its results should be used.