Mathematical Models for Stress–strain Behavior of Nano Magnesia-Cement-Reinforced Seashore Soft Soil

: The stress–strain behavior of nano magnesia-cement-reinforced seashore soft soil (Nmcs) under different circumstances exhibits various characteristics, e.g., strain-hardening behavior, falling behavior, S-type falling behavior, and strong softening behavior. This study therefore proposes a REP (reinforced exponential and power function)-based mathematical model to simulate the various stress–strain behaviors of Nmcs under varying conditions. Firstly, the mathematical characteristics of different constitutive behaviors of Nmcs are explicitly discussed. Secondly, the conventional mathematical models and their applicability for modeling stress–strain behavior of cemented soil are examined. Based on the mathematical characteristics of different stress–strain curves and the features of different conventional models, a simple mathematical REP model for simulating the hardening behavior, modified falling behavior and strong softening behavior is proposed. Moreover, a CEL (coupled exponential and linear) model improved from the REP model is also put forth for simulating the S-type stress–strain behavior of Nmcs. Comparisons between conventional models and the proposed REP-based models are made which verify the feasibility of the proposed models. The proposed REP-based models may facilitate researchers in the assessment and estimation of stress–strain constitutive behaviors of Nmcs subjected to different scenarios.


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][2][3]. In geotechnical engineering, deep mixing method is generally adopted to improve the strength of soil by adding cement to the soil [4][5][6][7][8]. The soft soil is reinforced to impart higher strength through a series of reactions between raw materials and curing agent [9][10][11][12][13][14][15][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][18][19][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][21][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.
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 different 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 different 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.

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.

Mathematical Characteristics of τ-δ Behavior
The shear stress-shear strain (τ-δ) 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.

Strain-Hardening τ-δ Behavior
The typical shear stress-shear strain (τ-δ) 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 τ-δ curve shows a straight line and the shear stress increases linearly with the shear strain at a gradient of initial elastic modulus E0. In the plastic stage, the tangent modulus Ei reduces gradually as strain accumulates, leading to a nonlinear τ-δ curve. In the failure stage, the τ-δ curve flattens out and the shear stress reaches the ultimate shear strength τp with a corresponding shear strain δ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 τ-δ curve.

Falling τ-δ 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 τ-δ curve after the elastic stage shows an evident softening behavior. As shown in Figure 2, the typical τ-δ 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 τ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 τ-δ 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 τ-δ 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. Based on the τ-δ 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 τ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.

Mathematical Characteristics of σ-δ 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 σp corresponding to a strain of ε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 E0. 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.

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

Conventional Shear Stress-Displacement (τ-δ) Models
The conventional shear stress-displacement models have three types: hyperbolic model, exponential model and power function model.

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 where τ is the shear stress (kPa), δ is the shear displacement (mm), E0 is the nominal elastic modulus (kPa/mm) and τp is the shear strength (kPa). According to the τ-δ 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 τ-δ 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.

Exponential Model
The exponential model converges faster, and it is suitable for the τ-δ 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 Cai et al. [41] investigated the mathematical defects of the exponential model using the halfvalue 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.

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 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 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 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 τ-δ curve under a highly corrosive environment. Therefore, the second derivative of power function cannot be equal to zero.

Mathematical REP Model for τ-δ 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 where a, b, λ, and k are parameters to be determined, and a ≥ 0, b > 0, λ > 0, k > 0. When k = 0, Equation degrades into exponential function.
Taking the limit and the first derivative of Equation (6) leads to The progressive limit of τ-δ curve is shear strength, i.e., τp. Hence, The first derivative of the new model with δ=0 is the initial modulus of elasticity, i.e., E0. Combining with Equation (8) gives According to Equation (10), (E0/τp − λk) > 0 if b > 0. Based on the above, the REP model for the hardening τ-δ curve can be expressed as 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 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 τ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.

Mathematical Models for Stress-Displacement (σ-δ) Behavior
Likewise, the conventional mathematical model for the constitutive relation of stressdisplacement (σ-δ) 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 (σ-δ) behavior, only power function, quadratic function, and the proposed REP function are discussed.
For the power function of σ-δ behavior, it can be easily modified from Equation (3), that is   where ε is the strain (%), σ is the stress (kPa) and σ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 where σp and ε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 σ-δ behavior is slightly different from that for the τ-δ behavior, i.e., Equation (11) 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 It can be observed that the first and second derivatives of REP model are affected 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.

Hardening τ-δ 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.  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 τ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 τδ 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.

Modified Falling τ-δ Behavior
After adding Nm, the τ-δ 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.  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 τ-δ curve. In contrast, it is easy for a double-"X" type discrepancy to appear in the conventional mathematical model simulation.

Strong Softening σ-δ 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 σ-δ 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 differs 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 effect 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 great engineering value in predicting the compressive strength. 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 convergent 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 different acid erosion concentrations, with relatively high fitting accuracy. The fitting effect 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.
where x, y, and z are the parameters to be determined; αs is the concentration of sulfuric acid solution.
In accordance with 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 Combining with Equations (15) and (19), the damage model of sulfuric acid erosion can be expressed as As shown in Equation (20), the damage model only contains parameters of sulfuric acid concentration (αs), cohesion (c) and internal friction angle (φ), making the new model more suitable for acid erosion environment.

Mathematical CEL Model for Modified S-type Falling τ-δ Behavior
As discussed above, the REP model has fairly good fitting effect which can be used for different types of stress-strain curves. However, for the S-type curve which has inflection points, it is timeconsuming 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.

Mathematical CEL Model
The CEL model is actually a coupled exponential and linear function model. Setting λ = −1, Equation (6) can be simplified as where α ≥ 0, b > 0 and k > 0; k is defined as the inflection factor. Taking the limit of Equation (21) gives Since the S-type curve eventually converge to the shear strength τp so p .
Taking the first derivative of Equation (21) leads to If δ = 0, then The first derivative of the new model is the initial elastic modulus when the displacement is 0, i.e., δ = 0. Combining with Equation (23)   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 τ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.

Conclusions
This study analyses the experimental data of τ-δ and σ-δ constitutive relationships of Nmcs, and examines its characteristics under various conditions. Two mathematical models, i.e., REP (reinforced exponential and power function) and CEL (coupled exponential and linear) models, are proposed to overcome the shortcomings of conventional models in the modeling of hardening, modified falling, S-type falling, and strong softening stress-strain behaviors. Some key findings are summarized as bellow:


The τ-δ curves with varying cement mixing ratio, Nm mixing ratio and acid erosion concentrations exhibit different behaviors. For example, hardening behavior occurs with low cement mixing ratio; falling behavior takes place after adding Nm and S-type happens with acid erosion environment.  The proposed REP model is able to satisfy the mathematical characteristics of the stress-strain curves with hardening and modified falling behaviors as well as strong softening behavior.
Compared with conventional hyperbolic model, exponential and power function models which easily produce double-"X" discrepancies, the REP model has an evidently higher fitting accuracy in modeling both hardening and modified falling curves. In the modeling of strong softening behavior, the REP model also performs the best. In addition, by introducing an acid erosion factor, the REP model can be further expressed as a sulfuric acid erosion damage model which includes only three typical parameters, i.e., sulfuric acid concentration, cohesion of silty clay, and internal friction angle.  In the modeling of S-type stress-strain behavior, the conventional models are unable to simulate the characteristic with inflection point and produce even bigger "X" discrepancy than that in the modeling of hardening and falling behavior. In contrast, the proposed CEL model performs much better especially in modeling the range around the inflection point at the early stage.
The applicability of the proposed method of this study in practice mainly include two parts. As introduced in the study, the seashore soft soil discussed herein was collected from coastal areas in Shaoxing, Zhejiang Province, China. Thus, for the regions with similar seashore soft soil, e.g., the Yangtze river delta region of China, the reported methods in this study could be used directly. On the other hand, for regions with different seashore soft soil (i.e., the soil specific gravity, LL, PL vary a lot), the proposed expressions are still applicable if the specific soil stress-strain behavior was exhibited, e.g., strong softening behavior. However, it should be noted that, for the latter scenario, the corresponding calculation parameters should be re-calibrated through laboratory tests.