**1. Introduction**

Soft soil is widely distributed in coastal areas with many defects such as large natural moisture content, excessive compression capacity, and poor bearing capacity [1–3]. In geotechnical engineering, deep mixing method is generally adopted to improve the strength of soil by adding cement to the soil [4–8]. The soft soil is reinforced to impart higher strength through a series of reactions between raw materials and curing agen<sup>t</sup> [9–16]. In addition to the traditional cement curing agent, researchers are constantly looking for some novel materials such as nano materials to improve the bearing capacity of the soft soil layer [2,3,14,15,17–20].

Nano cemented soil refers to the cement-soil mixture which is improved by adding nano materials as admixtures into the mixture of water, cement and soil [20–22]. At present, the nano materials for enhancing cemented soil mainly include nano titanium oxide, nano montmorillonite, nano magnesia, nano silicon, nano alumina, etc. Previous experimental results show that adding nano admixtures can improve the performances of cemented soil such as its soil strength and anticorrosive properties. Based on some recent studies [2,3], nanometer magnesia (Nm) can be added into cemented soil to improve its mechanical performance.

In the past few years, mathematical models have been adopted to study the stress–strain response of soils [23–31]. In term of the stress–strain behaviors of cemented soil, a large number of researches have been reported [2,3,5,6,13–18,32–37]. However, the reported mathematical models have some limitations particularly in the modeling of the constitutive behavior of nano magnesia-cement-reinforced seashore soft soil (Nmcs) which manifests different stress–strain behaviors under varying conditions [2–4,17,18,38,39].

This study aims to propose a REP (reinforced exponential and power function)-based mathematical model to simulate the various stress–strain behaviors of Nmcs under varying circumstances. The mathematical characteristics of di fferent constitutive behaviors are firstly examined. Then, the conventional mathematical models for stress–strain behavior of cemented soil are discussed. Based on the mathematical characteristics of various stress–strain curves and the features of di fferent conventional models, a new REP model for characterizing hardening behavior, modified falling behavior and strong softening behavior is proposed. Furthermore, a CEL (coupled exponential and linear) model improved from the REP model is also proposed which is able to simulate the S-type stress–strain behavior of Nmcs. Comparisons between conventional models and the proposed REP-based models are made which verifies the feasibility of the proposed models.

#### **2. Mathematical Characteristics of Stress–Strain Constitutive Relations of Nmcs**

#### *2.1. Materials and Samples*

The seashore soft soil discussed in this study was collected from coastal areas in Shaoxing, Zhejiang Province, China. Its specific gravity, liquid limit, and plastic limit were 2.6, 43.5%, and 30%, respectively. According to American Society of Testing Materials, ASTM (1994), it belongs to silty clay [2]. In direct shear test, the thickness of the shear plane is assumed to be zero which means the corresponding shear strain cannot be calculated; thus in this case shear displacement is applied to represent the strain behavior. In the following, a more general symbol δ is therefore used to represent shear deformation characteristics, i.e., shear displacement or shear strain.

The soil samples were prepared under a standard maintenance temperature of 20 ◦C with a relative humidity of 95%. The mechanical tests were performed after 28 days' standard curing time. For the tests considering sulfuric acid erosion, after the standard curing the samples were immersed in sulfuric acid solution for another 14 days.

#### *2.2. Mathematical Characteristics of* <sup>τ</sup>*-*δ *Behavior*

The shear stress–shear strain (<sup>τ</sup>-δ) behavior of Nmcs typically comprises four types [2,3,17,18]. They are strain-hardening behavior, falling behavior, S-type falling behavior and strain-softening behavior. According to Wang et al. [17] the strain-softening behavior can be well captured using a generalized mathematical model, so it will not be discussed in this study. In order to establish the constitutive model characterizing the shear stress-displacement curve, it is necessary to analyze the mathematical characteristics.

#### 2.2.1. Strain-Hardening τ-δ Behavior

The typical shear stress-shear strain (<sup>τ</sup>-δ) curve for strain-hardening behavior of Nmcs is shown in Figure 1. As can be seen, the strain-hardening process can be divided into three stages: elastic stage (OA), plastic stage (AB), and failure stage (BC). In the elastic stage, the <sup>τ</sup>-δ curve shows a straight line and the shear stress increases linearly with the shear strain at a gradient of initial elastic modulus *E*0. In the plastic stage, the tangent modulus *Ei* reduces gradually as strain accumulates, leading to a nonlinear <sup>τ</sup>-δ curve. In the failure stage, the <sup>τ</sup>-δ curve flattens out and the shear stress reaches the ultimate shear strength <sup>τ</sup>p with a corresponding shear strain <sup>δ</sup>p. In this stage, the shear stress is mainly contributed by the friction resistance at the failure surface of the soil sample. In sum, the

strain-hardening curve includes four mathematical features: through the origin, monotone increasing, convex to τ-axis, and infinite convergence.

**Figure 1.** Strain-hardening <sup>τ</sup>-δ curve.

#### 2.2.2. Falling <sup>τ</sup>-δ Behavior

In the direct shear test of seashore soft soil after adding Nm, the cohesive force was lost after failure and hence the shear stress decreased dramatically. The falling <sup>τ</sup>-δ curve after the elastic stage shows an evident softening behavior. As shown in Figure 2, the typical <sup>τ</sup>-δ curve can be divided into four stages. Stage 1 is the elastic stage (OA) which is similar to that of strain-hardening curve. In this stage, the shear stress increases with the gradient of initial elastic modulus until reaching the failure. The corresponding shear stress and strain at point A denote the failure displacement δp and shear strength <sup>τ</sup>p, respectively. Stage 2 is the falling stage (AB), which evidently shows the falling of shear stress after the soil fails. This falling of shear stress is attributed to the loss of cohesion and the tangent modulus of the <sup>τ</sup>-δ curve gives a negative value. Stage 3 is the plastic stage (BC) wherein the friction resistance starts to dominate after the loss of cohesion. In this stage, the tangent modulus of the curve approximates the initial elastic modulus at the beginning which subsequently decreases due to the occurrence of plastic shear stress. Stage 4 is the residual shear stress stage (CD) where the tangent modulus approaches zero and the shear stress is equal to residual shear stress. It should be noted that the <sup>τ</sup>-δ falling curve can be modified by ignoring Stages 2 and 3 (AB and BC), then it returns to the strain-hardening behavior. However, the shear strength in the falling curve refers to the shear stress at the end of Stage 1 which is different from its counterpart in the strain-hardening curve.

**Figure 2.** Falling <sup>τ</sup>-δ curve.

#### 2.2.3. S-Type Falling <sup>τ</sup>-δ Behavior

Based on the <sup>τ</sup>-δ response of Nmcs in the scenario with high corrosion, the typical S-type falling curve is shown in Figure 3. As can be seen, the S-type falling curve consists of five stages. Stage 1 refers to the erosion softening stage (OA) where the surface of the sample is eroded by highly corrosive sulfuric acid and becomes relatively soft. Initially, the increment in the shear stress is relatively small as the displacement increases, giving rise to a relatively small initial modulus. As the shearing develops and enters the hard soil, which is less corroded, the tangent modulus tends to increase until stabilizing when accessing Stage 2, i.e., linear elastic stage (AB). The shear strength <sup>τ</sup>p and corresponding displacement δp at failure occur at the end of Stage 2, i.e., point B. After reaching failure, the falling stage (BC), the frictional plastic stage (CD) and the residual shear stress stage (DE) emerge continuously. These three stages are similar to those in the above falling behavior. Thus, the S-type falling curve can be modified in a similar way as that for the above falling curve. However, the obtained modified curve is still unlike the strain-hardening curve which is convex to τ-axis; namely, an inflection point occurs on the curve, e.g., before this point, the curve is convex to δ-axis, while after this point, it is convex to τ-axis.

**Figure 3.** S-type Falling <sup>τ</sup>-δ curve.

#### *2.3. Mathematical Characteristics of* <sup>σ</sup>*-*δ *Behavior*

Based on the UCS (unconfined compression strength test) results, the stress–strain curve under uniaxial compression exhibits a strong softening and brittle failure behavior, as shown in Figure 4. As can be seen, the stress–strain curve is hump-shaped where the peak point of the curve gives the unconfined compressive strength <sup>σ</sup>p corresponding to a strain of <sup>ε</sup>p. This curve comprises three stages. In the hardening stage (OP), the stress initially increases linearly with strain; the gradient of which represents the initial elastic modulus *E*0. This gradient starts to decrease in the later phase of the hardening stage until arriving at the peak, i.e., point P where the brittle failure takes place and the tangent modulus equals to zero. After that the stress reduces monotonically with strain in the softening stage (PM) and remains unchanged when entering the residual stress stage (MN). In sum, the mathematical features of the curve include: through the origin, having extreme value point, and converging at residual stress.

**Figure 4.** Axial stress–strain curve.

#### **3. Established Mathematical Model**

The characteristics of modified falling <sup>τ</sup>-δ curve are almost consistent with those of strain-hardening <sup>τ</sup>-δ curve. Therefore, this study attempts to establish a simple, mathematical constitutive model which is applicable for modeling both strain-hardening and modified falling <sup>τ</sup>-δ behaviors of Nmcs.

#### *3.1. Conventional Shear Stress-Displacement (*<sup>τ</sup>*-*δ*) Models*

The conventional shear stress-displacement models have three types: hyperbolic model, exponential model and power function model.

#### 3.1.1. Hyperbolic Model

The most classical nonlinear model describing hardening curve is hyperbolic model, which is widely applied to the nonlinear modeling in various fields because of its simple expression, convenient simulation and easy determination of parameters. The hyperbolic model was proposed by Duncan and Chang [40] and its expression is

$$
\pi = \frac{\delta}{1/E\_0 + \delta/\tau\_\mathbf{p}} \tag{1}
$$

where τ is the shear stress (kPa), δ is the shear displacement (mm), *E*0 is the nominal elastic modulus (kPa/mm) and <sup>τ</sup>p is the shear strength (kPa).

According to the <sup>τ</sup>-δ curves of shear test on the soil-structure contact surface, the hardening curve is rigid and plastic, and the hyperbolic model is difficult to meet this condition. In contrast, the strain in the <sup>τ</sup>-δ curve of Nmcs is mostly elastic before failure and only a small amount of plastic strain occurs after failure. In addition, the hyperbolic curve produces a large error in fitting the softening curve, which causes the deviation between the fitting results and the measured data due to its mathematical characteristics. Therefore, the hyperbolic model may be not suitable for modeling both hardening and softening behaviors.

#### 3.1.2. Exponential Model

The exponential model converges faster, and it is suitable for the <sup>τ</sup>-δ curves with elastoplastic strain. Hence, it is superior to the hyperbolic model for hardening curves. The exponential model was firstly proposed based on the one-dimensional consolidation theory of Terzaghi, and was mainly applied to describe the stress–strain curve of soil. The mathematical constitutive function is

$$
\pi = \tau\_{\mathbb{P}} \left( 1 - e^{-E\_0 \delta / \tau\_{\mathbb{P}}} \right). \tag{2}
$$

Cai et al. [41] investigated the mathematical defects of the exponential model using the half-value strength index. The results showed that the half-value strength index of both the exponential and the hyperbolic function is a fixed value which is only related to the peak strength and the initial elastic modulus. In this regard, the simulated stress–strain curve is only applicable for specific cases. In contrast, the half-value strength index of the power function is a parameter with a non-fixed value. The shape of its curve varies with the variation of the parameters, and hence it is widely applicable.

#### 3.1.3. Power Function Model

Due to the mathematical characteristics of the power function model, its simulation effect is remarkable in modeling plastic curves. In addition, the power function model of the nonlinear constitutive model of clay has the parameter of tangent modulus index. The expression of the power function model is

$$\tau = \tau\_{\rm P} \left\{ 1 - \left[ 1 + \frac{(\theta - 1)E\_0}{\tau\_{\rm P}} \delta \right]^{\frac{1}{1-\theta}} \right\} \tag{3}$$

where the value of θ is larger than 1; when θ = 2, the power function becomes the hyperbolic shear stress-displacement constitutive equation, so the hyperbolic function is only a special case of the power function. The first derivative of Equation (3) is given as

$$\frac{d\tau}{d\delta} = E\_0 \left[ 1 + \frac{(\theta - 1)E\_0}{\tau\_\text{P}} \delta \right]^{\frac{\theta}{1-\theta}}.\tag{4}$$

Since θ > 1, the first derivative of the power function is always greater than 0 in the domain of definition, and the curve increases monotonically. As the parameter θ changes, it is more applicable to the hardening curve and modified falling curve.

The second derivative of Equation (3) leads to

$$\frac{d^2\pi}{d\delta^2} = -\frac{E\_0^{-2}\theta}{\tau\_\mathcal{P}} \left[ 1 + \frac{(\theta - 1)E\_0}{\tau\_\mathcal{P}} \delta \right]^{\frac{2\theta - 1}{1 - \theta}}.\tag{5}$$

It can be found that the second derivative of power function is always less than zero, meeting the characteristic of hardening and modified falling curves, i.e., convex to τ-axis. However, the inflection point occurs on the S-type <sup>τ</sup>-δ curve under a highly corrosive environment. Therefore, the second derivative of power function cannot be equal to zero.

#### *3.2. Mathematical REP Model for* <sup>τ</sup>*-*δ *Behavior*

Based on the above conventional mathematical models, a new REP (reinforced exponential and power function) mathematical model for modeling the shear stress-displacement behavior of Nmcs under acid erosion environment is proposed in this study. This model combines exponential and power functions, which is expressed as

$$
\pi = a \Big[ 1 - e^{-b\delta} (1 + k\delta)^{-\lambda} \Big], \tag{6}
$$

where *a*, *b*, λ, and *k* are parameters to be determined, and *a* ≥ 0, *b* > 0, λ > 0, *k* > 0. When *k* = 0, Equation (6) degrades into exponential function.

Taking the limit and the first derivative of Equation (6) leads to

$$\begin{cases} \tau \big|\_{\delta=+\infty} = a \\ \frac{d\tau}{d\delta} = abe^{-b\delta} (1+k\delta)^{-\lambda} + a\lambda ke^{-b\delta} (1+k\delta)^{-\lambda-1} \end{cases} \tag{7}$$

τ =

The progressive limit of <sup>τ</sup>-δ curve is shear strength, i.e., <sup>τ</sup>p. Hence,

$$a = \tau\_{\mathbb{P}}.\tag{8}$$

The first derivative of the new model with δ=0 is the initial modulus of elasticity, i.e., *E*0. Combining with Equation (8) gives

$$
\tau\_\mathsf{P} b + \tau\_\mathsf{P} \lambda k = E\_{0,\prime} \tag{9}
$$

where

$$b = \frac{E\_0}{\tau\_\mathbb{P}} - \lambda k.\tag{10}$$

 as

> (11)

According to Equation (10), (*E*0/<sup>τ</sup>p − λ*k*) > 0 if *b* > 0. Basedontheabove,theREPmodelforthehardening<sup>τ</sup>-δcurvecanbe

 − Use, the REP model for the hardening \(\tau\text{-}\delta\) curve can be expressed as 
$$\tau = \tau\_{\rm P} \left[ 1 - e^{-\left(\mathbb{E}\_{0}/\tau\_{\rm P} - \lambda k\right)\delta} (1 + k\delta)^{-\lambda} \right].$$

In order to analyze whether REP model satisfies the mathematical characteristics of hardening curve, the zero point, the limit, the first derivative and the second derivative are discussed as below

$$\begin{cases} \pi|\_{\delta=0} = 0\\ \pi|\_{\delta=+\infty} = \pi\_{\mathsf{P}}\\ \frac{d\mathsf{T}}{d\delta} = E\_0 \big( \mathbb{E}\_0/\tau\_{\mathsf{P}} - \lambda k \big) e^{-(E\_0/\tau\_{\mathsf{P}} - \lambda k)\delta} (1+k\delta)^{-\lambda} + E\_0 \lambda k e^{-(E\_0/\tau\_{\mathsf{P}} - \lambda k)\delta} (1+k\delta)^{-\lambda-1} \\\ \frac{d^2\mathsf{T}}{d\delta^2} = -E\_0 \big( \mathbb{E}\_0/\tau\_{\mathsf{P}} - \lambda k \big)^2 e^{-(E\_0/\tau\_{\mathsf{P}} - \lambda k)\delta} (1+k\delta)^{-\lambda} \\\ -2E\_0 \big( E\_0/\tau\_{\mathsf{P}} - \lambda k \big) \lambda k e^{-(E\_0/\tau\_{\mathsf{P}} - \lambda k)\delta} (1+k\delta)^{-\lambda-1} - E\_0 \lambda (\lambda+1)k^2 (1+k\delta)^{-\lambda-1} \end{cases} \tag{12}$$

When δ = 0, the shear stress equals to 0 so the curve passes through the origin, satisfying the first characteristic of the hardening curve. The coefficients in the first-order derivative equation are all positive, and the first-order derivative is always greater than 0 in the domain of definition, which satisfies the characteristic of monotonic increase. When the displacement δ approaches infinity, shear stress gradually gets close to shear strength of <sup>τ</sup>p. Therefore, the new model goes through the origin and has both upper and lower bounds. The coefficients of the second derivative equation are all negative, and the second derivative of the new model is always less than 0. Hence, the REP model is theoretically suitable for the hardening and modified falling curves.

#### *3.3. Mathematical Models for Stress-Displacement (*<sup>σ</sup>*-*δ*) Behavior*

Likewise, the conventional mathematical model for the constitutive relation of stress- displacement (<sup>σ</sup>-δ) behavior also has many types such as hyperbolic, exponential function, power function, piecewise function and quadratic function. As aforementioned, the power function model has a better fitting effect than hyperbolic and exponential functions, while the piecewise function has some defects, e.g., it is troublesome to fit and has many parameters. Therefore, in the study of models for stress-displacement (<sup>σ</sup>-δ) behavior, only power function, quadratic function, and the proposed REP function are discussed.

For the power function of <sup>σ</sup>-δ behavior, it can be easily modified from Equation (3), that is

$$
\sigma = \sigma\_{\rm P} \left\{ 1 - \left[ 1 + \frac{(\theta - 1)E\_0}{\sigma\_{\rm P}} \varepsilon \right]^{\frac{1}{1-\theta}} \right\} \tag{13}
$$

where ε is the strain (%), σ is the stress (kPa) and <sup>σ</sup>p is the peak value of the stress–strain curve (UCS).

The quadratic model is able to simulate the stress–strain curve of compacted cement soil with an obvious peak value; its expression is

$$
\sigma = \sigma\_{\rm P} \left[ A \frac{\varepsilon}{\varepsilon\_{\rm P}} - B \left( \frac{\varepsilon}{\varepsilon\_{\rm P}} \right)^2 \right],
\tag{14}
$$

where <sup>σ</sup>p and <sup>ε</sup>p are the maximum stress and the corresponding strain, respectively. *A* and *B* are the fitting parameters to be determined.

When ε = 0, the stress of power and quadratic models is 0, which satisfies the characteristic of the stress–strain curve passing through the origin. As aforementioned, in the process of power function fitting, it is unable to converge in the failure stage, while the strong softening stress–strain curve will soften immediately after reaching the peak failure. Therefore, the convergence of the power function may be not timely, which may lead to a large deviation from the measured curve. Although the quadratic model is able to converge in time after the peak value; when the strain approaches infinity, the stress is also infinite, which theoretically does not meet the characteristics of infinite convergence of the stress–strain curve.

The expression of REP model for the <sup>σ</sup>-δ behavior is slightly di fferent from that for the <sup>τ</sup>-δ behavior, i.e., Equation (11)

$$
\sigma = a \Big[ 1 - e^{-h\varepsilon} (1 + k\varepsilon)^{-\lambda} \Big],
\tag{15}
$$

where *a*, *b*, λ*,* and *k* are parameters to be determined, and the value range of each parameter should be suitable for the hump curve. When ε = 0, the σ value of the REP model also equals to zero which satisfies the characteristic of passing through the origin.

The first and second derivatives of the stress–strain curve are as follows

$$\begin{cases} \frac{d\sigma}{d\delta} = ab\varepsilon^{-b\varepsilon}(1+k\varepsilon)^{-\lambda} + a\lambda k\varepsilon^{-b\varepsilon}(1+k\varepsilon)^{-\lambda-1} \\\frac{d^2\sigma}{d\delta^2} = -ab^2\varepsilon^{-(E\_0/\tau\_\mathbf{p}-\lambda k)\delta}(1+k\varepsilon)^{-\lambda} - 2ab\lambda k\varepsilon^{-b\varepsilon}(1+k\varepsilon)^{-\lambda-1} - a\lambda(\lambda+1)k^2(1+k\varepsilon)^{-\lambda-1} \end{cases} \tag{16}$$

It can be observed that the first and second derivatives of REP model are a ffected by the value of each parameter, and the positive and negative signs are uncertain which enables its adaptability to complex strain-softening curves. The specific value range of each parameter and the judgment of the positive and negative sign of the first derivative and the second derivative need further study.

#### *3.4. Application and Analysis*

⎪⎪⎪⎪⎪⎨

#### 3.4.1. Hardening <sup>τ</sup>-δ Behavior

Based on the measured data of direct shear test, comparisons of using hyperbolic, exponential and power function models as well as the proposed REP model are made. Two cases with 5% and 7% cement mixture ratio were examined. To determine the aforementioned model parameters, e.g., λ and *k*, four tests with vertical pressures of 100, 200, 300, and 400 kPa were carried out respectively. The comparison results for the cases under a vertical pressure of 400 kPa are shown in Figures 5 and 6.

**Figure 5.** <sup>τ</sup>-δ curve for mix ratio with 5% cement content.

**Figure 6.** <sup>τ</sup>-δ curve for mix ratio with 7% cement content.

As shown in Figure 5, on the whole, the four models have a good fitting effect, but the elastic stage of the measured curve is not smooth, and there is a prominent inflection point locally. In this stage, the four models have a large deviation from the measured value. The tangent modulus of the measured curve decreases gradually in the plastic stage, and the hyperbolic, exponential and power functions converge slowly, all of which appear below the curve while the fitting effect of REP model is very evident. When entering failure stage, both the hyperbolic and exponential fittings cannot converge to <sup>τ</sup>p, and lie in the upper part of the measured curve with big difference. In contrast, the power function converges well, and the fitting effect is better than the hyperbolic and exponential functions, although the end of the curve still appears above the measured curve. The REP model also has a good convergence effect and the whole fitted failure stage is very close to that of the measured curve.

As shown in Figure 6, the linearity of the measured curve in the elastic stage is not obvious while the hyperbolic, exponential, and power function curves are linear, deviating slightly from the measured curve. In the plastic deformation stage, the three functions appear below the measured curve, and the difference between the power function and the measured curve is the smallest. Similar to that observed in Figure 5, in the failure stage, the three models generally appear at the top of the curve. In contrast, the REP model is very close to the measured curve and shows an evidently favorable fitting effect.

Among the conventional models, the power function model is relatively superior in modeling <sup>τ</sup>-δ curves, particularly in the plastic stage where the computed results are closer to the measured curve and failure, and in the failure stage showing a faster convergence speed. Compared with the hyperbolic model, the exponential model has a better simulation effect. Moreover, the REP model shows superior fitting accuracy than the conventional mathematical models.

#### 3.4.2. Modified Falling <sup>τ</sup>-δ Behavior

After adding Nm, the <sup>τ</sup>-δ curve exhibits modified falling type. In the comparison of different mathematical models, two cases with Nm mixing ratios of 10% and 20% under a vertical pressure of 400kPa were selected. The results are shown in Figures 7 and 8.

**Figure 7.** <sup>τ</sup>-δ curve for mix ratio with 10% Nm content.

**Figure 8.** <sup>τ</sup>-δ curve for mix ratio with 20% Nm content.

As can be seen, the shape of the two measured curves is similar. Before failure, the modified curve is mostly elastic stage and there is no plastic transition stage showing significant linearity. In this stage, the conventional models produce non-linear curves, while the measured curve approximates a horizontal line in the failure stage. The curves of conventional models locate below the measured curve in the former part and above the measured curve in the later part. In the two stages, the conventional model intersects with the measured curve, exhibiting two "X" shapes. In general, the difference between conventional models and measured curve is relatively big while the REP model fits fairly well with the measured curve. Therefore, the proposed REP model has significant advantage in modeling the modified <sup>τ</sup>-δ curve. In contrast, it is easy for a double-"X" type discrepancy to appear in the conventional mathematical model simulation.

#### 3.4.3. Strong Softening <sup>σ</sup>-δ Behavior

To compare the performance of the power function and quadratic function models as well as the REP model, UCS experimental results for cases with 0 mol/L and 0.08 mol/L H2SO4 erosion were adopted, Figure 9. As can be seen, the power function model behaves the worst in the modeling of the hump-shaped <sup>σ</sup>-δ curve. For instance, it is convex at the rising stage of the curve, which is contrary to the characteristic of concave in the measured curve. In addition, it cannot simulate the softening behavior as the measured curve although it converges to an almost constant value in the end. The quadratic function model performs better than the power function model. It is able to simulate the softening behavior and give a peak stress value. However, the shape of the quadratic curve is convex in all stages, which di ffers with the shape of measured curves. Worse still, there exists a large discrepancy in the maximum stress value and its corresponding strain value between the results of quadratic model and the measured results. This may raise an adverse e ffect in predicting UCS results. In contrast, the results of proposed REP model agree incredibly well with the measured results. The REP model not only shows consistent shapes, but also predicts a fairly close maximum stress value and a corresponding strain value. This is of grea<sup>t</sup> engineering value in predicting the compressive strength.

**Figure 9.** Comparison of stress–strain curves: (**a**) pure water condition (0 mol/L H2SO4); (**b**) 0.08 mol/L H2SO4.

To further examine the performance of the REP model, a series of comparisons for cases with varying erosion concentrations (i.e., 0.00 mol/L, 0.02 mol/L, 0.04 mol/L, 0.06 mol/L and 0.08 mol/L) is conducted, Figure 10. As can be seen, in the stress rising stage, the REP model fits fairly well with measured results; in the reducing stage, slight discrepancy occurs for the cases with 0.00 mol/L and 0.02 mol/L erosion concentrations. This could be attributed to the scattered data of the measured results. In the final residual stage, the REP model deviates from the measured data when the convergen<sup>t</sup> residual stress was relatively large, such as the cases in pure water (i.e., 0.00 mol/L) and 0.02 mol/L acid erosion environment. When the residual stress was relatively small in the environment with a high concentration of acid erosion, the REP model gives better performance in fitting. Therefore, it can be found that the REP model can adapt to the measured curve under di fferent acid erosion concentrations, with relatively high fitting accuracy. The fitting e ffect is more evident under a high concentration erosion environment. Table 1 provides the values of the four parameters in the REP model for these five cases.

**Figure 10.** Comparison of stress–strain curves of different H2SO4 erosion concentrations with REP model fittings.

**Table 1.** Four parameter values of REP model under different erosion concentrations.


In order to study the influence of erosion concentration on REP model parameters, an acid erosion factor was introduced

$$
\omega(\mathbf{x}, \mathbf{y}, z, a\_s) = \frac{\pi^2 \mathbf{x}}{2\pi y + 4za\_s},
\tag{17}
$$

where *x*, *y,* and *z* are the parameters to be determined; αs is the concentration of sulfuric acid solution.

In accordance with Table 1, substituting the acid erosion factor and acid erosion concentration αs into the four parameters of REP model leads to

$$\begin{cases} a = 8q \tan[\omega\_d(a\_s) + 1.84c] - 1.84c \\ b = -0.45c - 2.18q \cos[\omega\_b(a\_s) + 2.18q] \\ k = 0.11 \tan[\omega\_k(a\_s) + 0.16] - 0.16 \\ \lambda = -2.2c - a\_s - 2.1c \sin[\omega\_\lambda(a\_s) + 2.2c] \end{cases} \tag{18}$$

where *c* and ϕ are the cohesion and internal friction angle of silty clay in the experiment, respectively. Hence, the acid erosion factors of the four parameters can be obtained

$$\begin{cases} \boldsymbol{\omega}\_{\rm{h}}(\boldsymbol{a}\_{\rm{s}}) = \boldsymbol{\omega}(-7.96\boldsymbol{\varphi}, 0.881, 1.84\boldsymbol{c}, \boldsymbol{a}\_{\rm{s}})\\ \boldsymbol{\omega}\_{\rm{b}}(\boldsymbol{a}\_{\rm{s}}) = \boldsymbol{\omega}(2.18217\boldsymbol{\varphi}, -0.100816, 0.44586\boldsymbol{c}, \boldsymbol{a}\_{\rm{s}})\\ \boldsymbol{\omega}\_{\rm{k}}(\boldsymbol{a}\_{\rm{s}}) = \boldsymbol{\omega}(0.1111, -6.592\boldsymbol{\varphi}\times 10^{-6}, 0.158, \boldsymbol{a}\_{\rm{s}})\\ \boldsymbol{\omega}\_{\rm{h}}(\boldsymbol{a}\_{\rm{s}}) = \boldsymbol{\omega}(2.095\boldsymbol{c}, 2.1927\boldsymbol{c}, -105.507\boldsymbol{c}, \boldsymbol{a}\_{\rm{s}}) \end{cases} \tag{19}$$

Combining with Equations (15) and (19), the damage model of sulfuric acid erosion can be expressed as

$$\sigma = a\_{[\omega\_k(a\_\delta), a\_\delta]} \left[ 1 - e^{-b\_{[\omega\_k(a\_\delta), a\_\delta]} \ell} \left( 1 + k\_{[\omega\_k(a\_\delta), a\_\delta]} \varepsilon \right)^{-\lambda\_{[\omega\_\lambda(a\_\delta), a\_\delta]}} \right]. \tag{20}$$

As shown in Equation (20), the damage model only contains parameters of sulfuric acid concentration (<sup>α</sup>*s*), cohesion (*c*) and internal friction angle ( ϕ), making the new model more suitable for acid erosion environment.

#### **4. Mathematical CEL Model for Modified S-Type Falling** <sup>τ</sup>**-**δ **Behavior**

As discussed above, the REP model has fairly good fitting e ffect which can be used for di fferent types of stress–strain curves. However, for the S-type curve which has inflection points, it is time-consuming to calculate the zero point of the second derivative equation of REP model and determine the position of inflection points. To solve this, the REP is therefore simplified and a new CEL (coupled exponential and linear) model is put forth.

#### *4.1. Mathematical CEL Model*

The CEL model is actually a coupled exponential and linear function model. Setting λ = −1, Equation (6) can be simplified as

$$
\pi = a \Big[ 1 - e^{-h\delta} (1 + k\delta) \Big],
\tag{21}
$$

where α ≥ 0, *b* > 0 and *k* > 0; *k* is defined as the inflection factor.

Taking the limit of Equation (21) gives

$$\left. \tau \right|\_{\delta = +\infty} = a \left[ 1 - e^{-b\delta} (1 + k\delta) \right] \Big|\_{\delta = +\infty} = a(1 - 0) = a. \tag{22}$$

Since the S-type curve eventually converge to the shear strength <sup>τ</sup>p so

$$a = \tau\_{\mathbb{P}}.\tag{23}$$

Taking the first derivative of Equation (21) leads to

$$\frac{d\pi}{d\delta} = abe^{-b\delta}(1+k\delta) - ake^{-b\delta}.\tag{24}$$

If δ = 0, then

$$
\frac{d\pi}{d\delta}|\_{\delta=0} = a(b-k). \tag{25}
$$

The first derivative of the new model is the initial elastic modulus when the displacement is 0, i.e., δ = 0. Combining with Equation (23) gives

$$
\pi\_\mathbb{P}(b-k) = E\_0. \tag{26}
$$

Hence, the parameter b can be derived

$$b = \frac{E\_0}{\tau\_\mathrm{P}} + k\_\mathrm{\prime} \tag{27}$$

and the CEL model can be expressed as

$$
\pi = E\_0 \left[ 1 - \varepsilon^{-(\frac{F\_0}{\mathbf{r}\_\mathbf{P}} + k)\delta} (1 + k\delta) \right]. \tag{28}
$$

⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨

In order to explore whether the CEL model satisfies the mathematical characteristics of the S-type curve, its zero point, limit, first derivative, and second derivative are discussed

$$\begin{cases} \tau \Big|\_{\delta=0} = E\_0[1 - 1 \cdot (1 - 0)] = 0\\ \tau \Big|\_{\delta \to +\infty} = \tau\_\mathrm{P} \\ \frac{d\tau}{d\delta} = \tau\_\mathrm{P} \Big|\_{\tau\_\mathrm{P}} + k \Big| e^{-(\frac{E\_0}{\tau\_\mathrm{P}} + k)\delta} \quad (1 + k\delta) - \tau\_\mathrm{p} k e^{-(\frac{E\_0}{\tau\_\mathrm{P}} + k)\delta} \\ \frac{d^2 \tau}{d\delta^2} = -\tau\_\mathrm{P} \Big| \frac{E\_0}{\tau\_\mathrm{P}} + k \Big|^2 e^{-(\frac{E\_0}{\tau\_\mathrm{P}} + k)\delta} \quad (1 + k\delta) + 2\tau\_\mathrm{P} \Big| \frac{E\_0}{\tau\_\mathrm{P}} + k \Big| e^{-(\frac{E\_0}{\tau\_\mathrm{P}} + k)\delta} \end{cases} \tag{29}$$

According to Equation (29), the CEL model goes through the origin and converges to <sup>τ</sup>p at infinity. Combining the similar items in the first derivative of Equation (29) gives

$$\frac{d\tau}{d\delta} = \left(E\_0 + E\_0 k \delta + \tau\_\mathrm{p} k^2 \delta\right) e^{-\left(\frac{E\_0}{\tau\_\mathrm{p}} + k\right)\delta} > 0. \tag{30}$$

Hence, the first derivative of CEL model is always greater than 0 in the domain of definition, and the curve increases monotonically.

The displacement at zero point of the second derivative of Equation (29) is obtained by

$$\delta\_{\mathbb{C}} = \frac{\tau\_{\mathbb{P}}k - E\_0}{\tau\_{\mathbb{P}}k^2 + E\_0k}. \tag{31}$$

To analyse the concavity of CEL model at zero point of second derivative, the relations are

$$\begin{cases} \frac{d^2 \tau}{d\delta^2} > 0, \left(\delta < \frac{\tau\_{\rm P}k - E\_0}{\tau\_{\rm P}k^2 + E\_0k}\right) \\\ d^2 \tau < 0, \left(\delta > \frac{\tau\_{\rm P}k - E\_0}{\tau\_{\rm P}k^2 + E\_0k}\right) \end{cases} \tag{32}$$

When δ < δc, τ" > 0 and the curve is convex to the δ-axis; when δ > δc, τ" < 0 and the curve is convex to the τ-axis.

#### *4.2. Application and Analysis*

Comparisons of hyperbolic, exponential and power function models and the above CEL model are carried out based on the experimental results of two cases with 0.09 mol/L H2SO4 erosion under vertical pressures of 300 kPa and 400 kPa, respectively. The results are shown in Figures 11 and 12.

**Figure 11.** <sup>τ</sup>*-*δ curve under 300 kPa vertical pressure.

**Figure 12.** <sup>τ</sup>-δ curve under 400 kPa vertical pressure.

As can be seen, the curves of conventional models are well banded which are above the measured curve before the inflection point and below the measured curve after the inflection curve. In contrast, the fitting effect of CEL model is obviously better than that of the other conventional models. It is basically consistent with the measure curve in the S-shaped range and converges faster and closer to <sup>τ</sup>p than the other models in the failure stage. However, at the initial stage, the CEL curve is slightly lower than the measure curve, while in the failure stage, the difference between the CEL curve and the measured curve is relatively evident.
