Finite Element Analysis of Expansive Bedrock Considering Electro-chemo-mechanical Coupling Phenomena in Crystal Layers of Clay Minerals and Internal Structural Degradation

This research proposes a novel model targeting expansive bedrock by combining an elastoplastic model considering electro-chemo-mechanical phenomena with the Cam-clay model, based on considerations of cementation and its degradation owing to plastic deformation. Additionally, tunnel excavation and swelling analyses are demonstrated by incorporating the proposed model into a finite element analysis code. Specifically, swelling analysis is conducted by decreasing the bedrock ion concentration and by considering the groundwater supply to the tunnel. The simulation results indicate that considerable swelling occurs in a limited area, such as the bottom part of the tunnel, even though changes in ion concentration are uniformly set in the entire analysis area; moreover, swelling phenomena develop instantaneously analogous to the observational records. The simulation reveals that shear stress increases while excavating an area and that softening with plastic expansion occurs in bedrock skeleton owing to decreasing ion concentration. In addition, the yielding of the bedrock skeleton triggers a rapid internal structural degradation and swelling of the clay mineral interlaminar layers based on simulations. Thus, this study clarifies the tunnel swelling deformation mechanism as an interaction problem between coupled electro-chemo-mechanical phenomena occurring in crystal layers and a plastic phenomenon involving the internal structural degradation of the bedrock skeleton. A model for expansive bedrock is proposed by considering electro-chemo-mechanical coupling phenomena in crystal layers of clay minerals and internal structural degradation. Tunnel excavation and swelling analyses are demonstrated by incorporating the proposed model into a finite element analysis code. Interactions between the expansion of interlaminar layers due to decreasing ion concentration and the plastic deformation of the rock skeleton causes swelling deformation. Significant tunnel swelling deformation occurs approximately 5 m below the tunnel base, because high-shear stress is applied subject to the influence of excavation. Yielding of the rock skeleton triggers rapid tunnel swelling deformation owing to the groundwater supply, and internal structural degradation accelerates the deformation. A model for expansive bedrock is proposed by considering electro-chemo-mechanical coupling phenomena in crystal layers of clay minerals and internal structural degradation. Tunnel excavation and swelling analyses are demonstrated by incorporating the proposed model into a finite element analysis code. Interactions between the expansion of interlaminar layers due to decreasing ion concentration and the plastic deformation of the rock skeleton causes swelling deformation. Significant tunnel swelling deformation occurs approximately 5 m below the tunnel base, because high-shear stress is applied subject to the influence of excavation. Yielding of the rock skeleton triggers rapid tunnel swelling deformation owing to the groundwater supply, and internal structural degradation accelerates the deformation.


Introduction
Pore fluid composition affects the mechanical behavior of bedrock containing expansive clay minerals. Moreover, swelling of expansive bedrock often causes problems during tunnel excavation and earth cutting. For example, the deformation of the Mount (Mt.) Sakazuki Tunnel on the Yamagata Expressway in Japan led to bumps and cracks which emerged instantaneously on the road due to the swelling of expansive bedrock approximately 17 years after the initiation of the service. The road surface uplift was up to 494 mm, and it developed at a rapid rate of 10 mm/day from the onset of deformation. These tunnel deformations constitute a severe problem because of the subsequent significant costs for necessary tunnel modifications. Thus, it is necessary to elucidate the mechanism behind the emergence of bumps and cracks as well as quantitatively evaluate the deformation caused by swelling.
Electro-chemo-mechanical interactions between the crystal layers of expansive clay minerals are one of the predisposing factors for the swelling phenomenon in bedrock. Di Maio (1996) conducted one-dimensional compression tests on Ponza bentonite in various stress conditions by changing the cell fluid from distilled water to NaCl solution, and vice versa, and researched the osmotic compression and swelling behavior of expansive soils. Figure 1 shows that swelling and osmotic compression occur when the ion concentration of pore fluid decreases and increases, respectively; moreover, they depend on confining pressure. The arrows in Fig. 1 represent the moment at which the distilled water solution is replaced by the NaCl solution.
Various mechanical models focusing on the microscopic behavior at the surface of clay minerals have been proposed to explain consistently the characteristic behavior of expansive soils. Alonso et al. (1990) proposed the Barcelona expansive model that targeted unsaturated expansive clay. This model assumes that expansive clays have double-porous structures consisting of a macrostructure (soil skeleton) and a microstructure (mineral level), and that the behavior of the macrostructure can be expressed as a function of the effective stress for the saturated soil. Subsequently, Gens and Alonso (1992) extended the Barcelona basic model, and targeted nonexpansive soils. Musso et al. (2013) proposed a model for compacted active clay based on the concept of a double-porous structure, and described a microstructure using the van't Hoff equation (van't Hoff 1884). Kyokawa et al. (2020) proposed an elastoplastic constitutive model of expansive soils that focused not only microphenomena, such as electro-chemical phenomena on the surface of clay minerals but also on macrophenomena, using the Cam-clay model.
As such, various studies have been conducted on swelling clays. However, the swelling behavior of bedrock is more complicated because of various factors, such as high stiffness and internal structural degradation, which are absent in clays. Additionally, there are few studies on mechanical modeling of bedrock containing swelling clay minerals which is relevant to actual tunnel construction. Only a few models that describe the behavior of swelling bedrock simultaneously account for the electro-chemical microphenomena on mineral surfaces and macroscopic reduction in cementation caused by internal structural degradation.
Our study proposes a novel model targeting expansive bedrock that combines an elastoplastic model considering electro-chemo-mechanical phenomena (Kyokawa et al. 2020) in conjunction with the use of the Cam-clay model, and based on considerations of cementation and its degradation owing to plastic deformation (Yamada et al. 2022). In addition, an algorithm is constructed to implicitly update stresses and internal variables for numerical analyses using the proposed model. Finally, this research investigates the mechanism of swelling phenomena in bedrock that are responsible for problems in tunnel construction based on finite element analysis using the proposed model.  Maio (1996)

Electro-chemo-mechanical Phenomena in Crystal Layers of Clay Minerals
An elastoplastic model considering the electro-chemomechanical phenomena in the crystal layers of clay minerals (Kyokawa et al. 2020) is used herein to investigate swelling bedrock. This section describes the model framework. This model assumes that swelling soils are composed of expansive clay minerals and nonexpansive soil particles, and possess a double-porous structure consisting of soil skeleton and interlaminar voids, as shown in Fig. 2.
Variables related to the interlaminar and soil skeleton aspects are indicated by a sub/superscripts "il" and "ss" , respectively; moreover, the variables representing volume are illustrated in Fig. 2. In this model, the void ratios (total void ratio e , interlaminar void ratio e il , and soil skeleton void ratio e ss ) are defined as When the volume changes from V 0 to V , the volumetric strain v (compression is positive) is expressed as Consequently, the interlaminar volumetric strain is defined as where d is equal to half the distance of the interlaminar space, d 0 is its initial value, S is the specific surface area of the crystals, n is the number of mineral layers, and N is the number of minerals in the soil element (Fig. 3).
Similarly, the soil skeleton volumetric strain ss v is defined as follows: From Eqs. (2)-(4), the total volumetric strain in the doubly porous structure is defined in terms of d as follows: where * is defined in terms of e ss 0 and e il 0 as follows: Therefore, the total strain using infinitesimal deformation theory is defined as the sum of the strain in the soil skeleton ss and isotropic strain with changing interlaminar distance, as follows: where is a second-order unit tensor, and its components are expressed by the Kronecker delta ij .
The swelling behavior is determined by the changes in d , and d is determined by the following equilibrium equation of the electro-chemo-mechanical forces between the crystal layers:  where c is the ion concentration of the pore fluid, f a is the Van der Waals force, f r is the osmotic pressure from the double diffusion layers, f h is the hydration force, and f e is the force from the (mean) effective stress . Hereafter, Eq. (8) is referred to as "the interlaminar equilibrium." The osmotic pressure, Van der Waals force, and hydration force are monotonically decreasing functions with respect to the interlaminar distance; moreover, the osmotic pressure changes with ion concentration. The osmotic repulsive pressure increases with decreasing ion concentration and vice versa. Therefore, increasing or decreasing the interlaminar distance can describe swelling or osmotic compression. The interlaminar equilibrium is a function of the interlaminar distance d , ion concentration c , and the effective stress . The characteristic mechanical behavior of swelling clay minerals, which depends on the pore fluid composition and confining pressure, can be expressed by adding the effective stress to the electro-chemical force. The osmotic repulsive pressure and Van der Waals force are based on the equations given by Komine and Ogata (1996), and the hydration force is based on the equations given by Afzal et al. (1984). Each of these forces and relevant other variables are defined as follows: and Boltzmann constants, respectively. These parameters are examined on a mineral-by-mineral basis or are physical constants.p = (tr )∕3 is the mean effective stress. Moreover, the shear stress q = √ 3∕2‖ ‖ with = − p and stress ratio = √ 3∕2‖ ‖ = q∕p with = ∕p are defined as such. In addition to the above equilibrium model, bedrock skeleton behavior is described with Cam-clay model considering cementation and its degradation in Sect. 2.2. Kyokawa et al. (2020) described the soil skeleton behavior with the Cam-clay model to target expansive soils. However, the target of this study is not soft expansive soil but swelling bedrock. Therefore, a method that considers the cementation effect and its degradation proposed by Yamada et al. (2022) is introduced into the Cam-clay model to describe the degradation of internal structure due to swelling. This model determines the tensile strength owing to the cementation effect of the material in soil skeletons, such as cementtreated soil and bedrock, by translating the yield surface to the tensile side, as shown in Fig. 4a. It can also express high cementation, because confining pressure virtually increases by treating translations of the yield surface as modifications of the effective stress. Moreover, it describes the transition from cement-treated soil and bedrock into normal soil by eliminating the deviation of the yield surface from the origin owing to progressive plastic deformation. Specifically, the model is expressed with modified effective stress , which is the effective stress plus the amount of isotopically parallel shifts of the yield surface Ψ:

Cementation and Internal Degradation of Swelling Bedrock
p , q , , , and , which are the derived quantities of , are defined in the same manner as those of . The yield surface with cementation passes through the origin in the modified effective stress space, as shown in Fig. 4b, such that the formulation of the model is common in many aspects to the standard Cam-clay model. Hereafter, quantities defined using the modified effective stresses considering cementation are denoted by an overline.
Ψ is referred to as cementation, and is assumed to asymptotically approach zero with plastic deformation. The evolution rule of Ψ is as follows: where d is the material constant that controls the decreasing rate of cementation and D = (̃−̃) (Mv ss 0 ) is the dilatancy coefficient. ̃ is the compression index, ̃ is the swelling index, M is the critical state constant, and v ss 0 is the initial value of a soil skeleton specific volume expressed as v ss = 1 + e ss .
This study applies the Cam-clay model by considering cementation and its degradation owing to plastic deformation in the soil skeleton of the elastoplastic model that satisfies the interlaminar equilibrium for a double-porous structure.

Formulation and Computational
Algorithm of the Proposed Model

Fundamental Equations of the Proposed Model
The fundamental equations of the modified Cam-clay model are given by and Λ is the plastic multiplier. In addition to these fundamental equations, the equations for the interlaminar equilibrium, double-porous structure, cementation, and its evolution rule mentioned in the previous section are used in this study Note that the ion concentration c of the interlayer pore fluid is a known condition for each incremental step.

Elastoplastic Computation Algorithm
In numerical calculations, the equations in the velocity form are approximated by backward difference as follows: (16) Strain rate:̇ ss =̇ e +̇ p , Interlaminar equilibrium: Modified effective stress: The evolution rule of cementation: ΔΨ = V(Ψ)‖̇ p ‖.
(26) Strain increment: Δ ss = Δ e + Δ p , Δ denotes that the quantity is an incremental amount. Variables of the n + 1 step are unknowns, and are computed by the updated algorithm described in the next section. In the formulation described above, p , , p c , and f use the value in the n + 1 step. Similarly, the subscripts which represent steps are omitted for variables other than those used to store history and for differential quantities.

Elastic Predictor Step
An elastic predictor step, which assumes ΔΛ = 0 and implies Δ p = , was performed at first. Moreover, the values obtained from this step are marked with the superscript "trial." In this step, the trial modified stress trial n+1 and trial interlaminar distance d trial n+1 , which simultaneously satisfy Eqs. (36)-(38), are obtained. Note that in this trial step, convergence calculation by the Newton-Raphson method is required, because the elastic modulus matrix ℂ e trial n+1 depends on stress Associated flow rule: Interlaminar equilibrium: Evolution rule of cementation: The first residual equation used in the Newton-Raphson method is obtained by substituting Eq. (36) into Eq. (37) as follows: where Δ il can be expressed as a function of d trial n+1 as follows: As Ψ does not change in the trial elastic calculation, effective stress trial n+1 can be expressed with modified effective stress trial n+1 as follows: , the interlaminar equilibrium is used as the second residual equation as follows: The corrections to stress and interlaminar distance in the Newton-Raphson method are obtained by solving the linearized equations of both residual equations as follows: , and C e, trial n+1 are given by The specific forms of f a ∕ d , f r ∕ d , and f h ∕ d are showcased in Appendix 7.1.
Upon convergence of the aforementioned iterative calculation, if the yield function value is negative, the loading judgment is designated as "unloading." Moreover, the internal variables are updated as n+1 = trial n+1 , d n+1 = d trial n+1 , and p v n+1 = p v n before terminating the update procedure. However, if the yield function value is positive, the loading judgment is designated as "loading," and a plastic corrector step is performed.

Plastic Corrector Step
The equations to be satisfied in the plastic corrector step are Eqs. (26)-(35). The number of unknowns in this step is reduced to three ( n+1 , ΔΛ , and d n+1 ) by eliminating variables from these equations. This step sets up three residual equations to find these unknowns by the Newton-Raphson method.
Combining Eqs. (26) and (32) gives the following: The first residual equation is obtained by substituting Eqs. (48) and (30)  where Δ il can be expressed as a function of d trial n+1 in the same way as in Eq. (40) as follows: Subsequently, the yield function, including the hardening rule, is defined as the second residual equation. As the plastic volumetric strain in n + 1 step can be expressed as Eq. (51) from the associated flow rule, the second residual equation is expressed by Eq. (52) as follows: Finally, the third residual equation is established. Substituting Eq. (30) into Eq. (34) gives the following: Note that Ψ n+1 can be expressed as a function of n+1 and ΔΛ . Therefore, f e ( n+1 ) is also a function of n+1 and ΔΛ , because n+1 can be expressed using Eq. (33) as follows: Consequently, the interlaminar equilibrium equation is designated as the third residual equation as follows: The corrections in the Newton-Raphson method are obtained by solving linearized equations of the three residual equations (Eqs. (49), (52), and (55)) as follows: where n+1 , C e n+1 , and f are given by The coefficients used in Eqs. (56) and (57) The modified stress, interlaminar distance, and plastic multiplier are determined by repeatedly correcting the solution until the values of the residual equations, i.e., Eqs. (56)-(58) yield values which are lower than the designated tolerances. After convergence, other hysteresis variables are updated using Eqs. (51), (53), and (54).

Consistent Tangent Modulus
A tangent modulus consistent with the aforementioned plastic corrector step is derived, such that the global Newton-Raphson method for satisfying the force equilibrium equation can converge to second order. Equations (49), (52), and (55) are linearized around current Δ . Herein, the coefficients are defined in the same manner as in the previous section The relations between d n+1 , d(ΔΛ) , and d(d) are obtained from these equations as follows: The following consistent tangent modulus for the plastic corrector step is obtained by substituting Eqs. (74) and (75) into Eq. (73) as follows:

Fig. 5 Calculation of unbalanced forces
A similar procedure is used to derive tangent coefficients consistent with the elastic predictor step. The following relations are obtained from linearizing Eqs. (39) and (42) around current Δ , as follows: The consistent tangent modulus for the elastic predictor step is obtained by substituting Eq. (79) into Eq. (78) as follows: The computational flow based on the above algorithm is summarized in Appendix 8. .

Tunnel Excavation and Swelling Analysis Methods
Tunnel excavation and swelling analyses were conducted by representing the model proposed in Sect. 2 with our inhouse finite element analysis code to achieve the main study objectives of (a) reproducing the qualitative tendency of the deformation in the Mt. Sakazuki Tunnel, and (b) elucidating the mechanism of rapidly developing swelling phenomena.
In this finite element analysis, the primary support was not modeled; only the swelling bedrock was analyzed. Figure 5 shows the considered case of excavating a circular hole in an area of ground. Moreover, first, it is assumed that there exists a horizontally stratified ground without a circular hole subjected to the action of overburden and dead weight with a known stress state satisfying the equilibrium equation. Here, the lateral pressure coefficient K 0 is also known. Unbalanced forces (= external nodal forces-internal nodal forces) act on the surface of the circular hole due to the presence of a circular hole in reality. These unbalanced forces are computed from the initial stress state and self-weight; moreover, tunnel excavation analysis is conducted by gradually decreasing these forces.

Swelling Analysis
Swelling analysis is conducted after the unbalanced forces are completely removed. Herein, the ion concentration c is lowered in the analysis area by assuming that groundwater is supplied from the surrounding ground toward the tunnel-cavity subject to the influence of excavation. When the ion concentration decreases, the osmotic pressure changes according to the interlaminar equilibrium, which induces swelling behavior. Accompanying this change, the state variables vary to satisfy the dominant equations, i.e., Eqs.
(16)-(25). This process is tracked by incremental analysis. In reality, the decrease in ion concentrations is thought to occur non-uniformly depending on the heterogeneity of groundwater flow within the bedrock. In this study, however, ion concentrations are reduced uniformly over the entire analysis region to distinguish the regions around the tunnel that are likely to be deformed, independent of the heterogeneity of the ion concentrations.  Figure 6 shows a finite element model (5040 nodes and 2400 elements). Note that 8-node hexahedral elements were used (Fig. 7). A tunnel was placed at the center of the analysis region (100 m × 100 m × 1 m), and calculations were performed by assuming planar strain conditions. Horizontal displacements were fixed on the left and right boundaries; furthermore, vertical displacements were fixed on the lower boundary. An equally distributed load of 1000 kPa, equivalent to the earth pressure of approximately 50 m, was applied to the top boundary to reproduce approximately 100 m of earth cover at the point where the deformation occurred in Mt. Sakazuki Tunnel.

Parameter Setting
The analyzed parameters were set based on the report for the deformation at Mt. Sakazuki Tunnel (Okui et al. 2009). The geological profile of the Mt. Sakazuki Tunnel was dominated by Cenozoic Tertiary Miocene rhyolite and homogeneous tuff. Soft rhyolitic tuffs subjected to hydrothermal alteration were distributed in the deformed section and are known to contain large amounts of smectite. Based on these considerations, stiffness and initial specific volume in bedrock are assumed to be the soft rock properties of volcanic origin. Table 1 lists the material constants. The specific volume is set uniformly for all elements to easily conduct self-weight analyses. The reference specific volume on normal consolidation line N , critical state constant M , compression index ̃ , swelling index ̃ , and reference stress p r were set based on the material properties of London clay (Atkinson and Bransby 1978). The parameters of the Cam-clay model were substituted with the properties of the London clay, because it was difficult to collect samples locally. The soil particle density of expansive clay sw was set as the material property of bentonite sampled in Yamagata Prefecture in Japan (Konishi et al. 2010). The soil particle density s , Poisson's ratio , and lateral pressure coefficient K 0 were set as the general material properties of nonexpansive clay. Table 2 lists the physical constants, and Table 3 lists the parameters for osmotic pressure, Van der Waals force, and hydration force (Kyokawa et al. 2020). In this research, only sodium ions were considered. Table 4 lists the initial conditions and other parameters. The initial specific volume was calculated from the dry density of pyroclastic rock, classified as soft rock (Japan Rock Mechanics Committee 2002). The initial cementation was set to be consistent with the estimated uniaxial compressive strength of the soft rock (Kobayashi and Adachi 1982). Additionally, the uniaxial strength was the assumed strength of soft rocks classified as fresh-to-slightly weathered based on the engineering classification of bedrock (Japan Rock Mechanics Committee 2004, JGS 3811-2004). The calculation method for the initial values is explained in Appendix 9. The degradation index of cementation was set with reference to Yamada et al. (2022). Initial and critical ion concentrations (lower limit) were set with reference to Komine and Ogata (2003). The mass ratio of expansive to nonexpansive   clay is a parameter used to determine the initial interlaminar void rate, and their relationship is explained in Appendix 10.

Computational Results and Discussion
First, the proposed model was validated to determine whether it can represent the mechanical behavior of swelling rock. Figures 8 and 9 show the distributions of cementation Ψ and volumetric strain in an interlaminar void after excavation and swelling analyses. Hereafter, a volumetric strain of an interlaminar void is referred to as the "expansion amount." In the swelling analysis following the excavation analysis, cementation degradation and swelling progress within a limited area, such as that just below the bottom part of the tunnel, despite the fact that the ion concentration was reduced uniformly in the entire analysis area. Figure 10 shows distributions of vertical displacement and expansion amount in bedrock just below the tunnel bottom. Depth from the bottom is plotted on the vertical axis, and the cumulative values are plotted at a fixed point approximately 10 m from the bottom. Herein, the vertical displacement is only due to swelling. Note that the expansion and vertical displacement increase remarkably from approximately 5 m just below the bottom. These figures show that the proposed model can describe the mechanical behavior of swelling bedrock adequately, because this tendency in the vertical displacement distribution matches the observed results at Mt. Sakazuki Tunnel, which are shown in Fig. 11. Herein, Fig. 11 shows the displacement distribution after lining observed 17 years after service entry. According to calculations completed separately, it was confirmed that the expansion amount and vertical displacement distribution were not affected considerably by ground surface loading. Figure 12 shows the relationship between change in ion concentration and expansion amount of Element A. Pronounced swelling is observed, as marked in Fig. 9d. The dashed and solid lines show the excavation and swelling processes, respectively. The expansion amount increases rapidly when the ion concentration falls below a certain value. Considering the changing ion concentrations as equivalent to evolving time, this figure suggests that deformation due to swelling develops in a short time when the ion concentration falls below a certain value, because the ion concentration changes in actual bedrock due to a continuous supply of groundwater. In fact, deformations due to swelling were reported to be sudden and grew rapidly in the Mt. Sakazuki Tunnel (Okui et al. 2009).
In the next section, we investigate this remarkable swelling behavior that occurs approximately 5 m below the bottom part of the tunnel and its rapid growth after a certain point in time. Figure 13 shows the shear stress distribution just below the bottom part of the tunnel after excavation. The shear stress reaches a minimum at 7 m below the bottom and increases rapidly from approximately 5 m below the bottom part of the tunnel. This distribution of the shear stress is similar to that of the expansion amount after swelling analysis  Fig. 10. This correlation suggests that as the shear stress generated during tunnel excavation increases, the swelling increases when ion concentration decreases.
Considering the stress paths of the two elements shown in Fig. 9d, Element A is located at the bottom of the tunnel where swelling is considerable; moreover, Element B is located at increased depths in which the shear stress after excavation is minimal. Figures 14 and 15 show both the stress path and relationship between the ion concentration and volumetric strain in Elements A and B, respectively. Stress paths are plotted in the p − q plane considering the In the stress paths, Element A (Fig. 14a) achieves a stress state above the critical state line at the initiation of swelling analysis owing to the influence of the large shear stress caused by excavation. During swelling analysis, the effective stress increases because of interlayer spreading as the ion concentration decreases and the bedrock skeleton yields outcomes above the critical state line. In the Cam-clay model, the area above the critical state line is the softening zone with plastic expansion (Asaoka et al. 1994). In the stress path of Element A, the shear and mean effective stresses decrease after yielding. These indicate that softening occurs in the region with pronounced swelling. As the mean effective stress decreases due to softening, the distance between the crystalline layers also increases according to the mechanism described by Eq. (8), thus resulting in significant swelling just below the tunnel base. The proposed model describes the high stiffness of bedrock owing to the introduction of cementation, which also increases shear stress during the excavation and decreases ion concentration; in turn, these events lead to softening with plastic expansion. In addition, Figs. 8 and 9 show that as the degradation of the internal structure becomes more extensive, swelling becomes more pronounced. This means that the cementation reduction caused by degrading internal structure promotes swelling.
Moreover, it becomes clear that the rapid swelling begins when the bedrock skeleton reaches yield after comparison of Fig. 14a, b. This illustrates that the changing ion concentrations induce the bedrock skeleton to yield; moreover, the softening initiated by the yielding of the bedrock skeleton leads to significant swelling of the interlaminar layers. Consequently, the yielding of the bedrock skeleton triggers rapid expansion.
In contrast, Fig. 15a shows that Element B is in an almost isotropically stress state after excavation and continues to be compressed elastically during the swelling analysis. In addition, the volumetric expansion is neither sudden nor significant like in Element A, because the behavior of the bedrock skeleton has a minor influence on the behavior of the interlaminar in Element B. Therefore, swelling is suppressed in Element B.
The results outlined above conclude that the distribution of shear stress during excavation below the tunnel's bottom part is such that the minimum point exists at approximately 6-7 m (below the bottom part of the tunnel) and that significant swelling tends to occur in the high-shear stress region which is at lower depths than the depth at which the shear stress is minimized. Moreover, numerical simulations reveal that the yielding of the bedrock skeleton triggers considerable swelling deformation produced by the interaction of coupled electro-chemo-mechanical phenomena in the crystalline layers of swelling clay minerals and elastoplastic phenomena with internal structural degradation within the rock matrix.

Conclusions
This study proposed a novel model targeting expansive bedrock by combining an elastoplastic model considering electro-chemo-mechanical phenomena with the Cam-clay model based on cementation considerations and its degradation due to plastic deformation. The proposed model was represented by finite element analysis code. Moreover, tunnel excavation and swelling analyses were conducted to clarify the characteristics and mechanism responsible for tunnel deformation in swelling bedrock. Some of the key results and insights are the following: (1) Tunnel deformation caused by swelling in crystal layers tends to occur just below the tunnel base. Moreover, extensive swelling starts from approximately 5 m just below the bottom part of the tunnel. (2) There was extensive swelling in a limited area approximately 5 m below the bottom part of the tunnel, because high-shear stress was already acting in this region after the excavation; additionally, softening with plastic expansion occurred above the critical state line when the ion concentration decreased. (3) During excavation, the uplift of the tunnel bottom grew rapidly after a certain time owing to decreasing ion concentration in the bedrock caused by groundwater supply to the tunnel. The yielding of the bedrock skeleton triggered this rapid uplift. (4) The mechanism behind the aforementioned observation began with swelling in crystal layers owing to the influence of decreasing ion concentration that caused the yielding of the bedrock skeleton. Additionally, as high-stress ratios were acting owing to tunnel excavation at depths just below the bottom part of the tunnel, the bedrock skeleton yielded above the critical state line, and caused softening with plastic expansion. To elaborate, a decrease in confining pressure is associated with softening, which in turn causes expansion between crystal layers. The swelling then progressed rapidly, accompanied by reduced cementation owing to internal structural degradation that resulted in tunnel deformation. Therefore, the rapid swelling in the bottom part of the tunnel was attributed to the interactions between coupled electro-chemo-mechanical phenomena which occurred between the crystal layers as well as elastoplastic phenomena which involved internal structural degradation of the bedrock skeleton.

A.1 Specific Forms of Linearization Terms in the Elastic Predictor Step
This section presents the specific forms of linearization terms used in the elastic predictor step in Sect. 3.2. First, the linearized forms of the residual equations in the elastic predictor step are Additionally, each linearization term is Partial differentiation of osmotic repulsive pressure, Van der Waals force, and hydration force by interlaminar distance are

C.1 Calculation of Initial Specific Volume
This section describes the method used to determine the initial specific volume. In this research, the initial specific volume was calculated from the dry density of a pyroclastic rock. Dry density d and initial bedrock skeleton void e ss 0 are related as follows: Therefore, from the dry density d and soil particle density s , the initial specific volume of the bedrock excluding interlaminar void is obtained as follows: The method of determination of initial interlaminar void is described in Appendix 10.

C.2 Calculation of Initial Cementation
This section describes the method used to determine the initial cementation. The uniaxial tensile strength of a rock is in the range of 5-20% of the uniaxial compressive strength (Kobayashi and Adachi 1982). In this research, the initial cementation is determined from the uniaxial tensile strength q u calculated from the rock's uniaxial compressive strength. Figure 16 shows the relationship between the stress path in a uniaxial tensile test and yield function considering a parallel shift. Herein, a uniaxial tensile strength test is assumed to perform with no confining pressure. In Fig. 16, the initial cementation Ψ 0 is equal to the deviation from the origin of the yield surface based on cementation considerations if the yield surface considering the parallel shift by cementation and the stress path of a uniaxial tensile test intersect at the point satisfying q = q u on the tensile side. Therefore, the initial cementation is determined by finding Ψ 0 , such that the two curves intersect at q = q u .

C.3 Calculation of the Size of the Initial Yield Surface
This section describes the method used to calculate the size of the initial yield surface. When the initial yield surface and initial stress state are as shown in Fig. 17, the initial (104) v ss 0 = 1 + e ss 0 = s d . and yielding states must lie on the same elastic wall in the p − q − v ss space. In this case, on p − v ss plane, Point A (initial state) and Point B (foot of the initial yield surface) have the positional relationship shown in Fig. 18. As points A and B lie on the same swelling line, it follows that: Moreover, as point B lies on the normal consolidation line  where v ss c0 is the soil skeleton specific volume on the NCL at p = p c0 .
The relationship obtained for the size of the initial yield surface by eliminating v ss and v ss c0 from Eqs. (105)- (107) is As showcased in Eq. (108), the size of initial yield surface p c0 is determined with material constants, initial stress state, and initial soil skeleton specific volume.

Appendix D: Relationship Between Mass Ratio and Initial Interlaminar Void
This section explains the relationship between mass ratio and initial interlaminar void. The simplified double structure in the initial period is shown in Fig. 19. A is the bottom area, a is the soil particle height of nonexpansive clay, b is the soil particle height of expansive soil, and c is the representative height of interlaminar void. Additionally, d 0 is the initial interlaminar distance determined from the interlaminar equilibrium with initial ion concentration and stress state, and t is the thickness of expansive clay minerals.
Herein, the initial interlaminar void e il 0 is described using the initial interlaminar volume V il v , soil particle volume of nonexpansive clay V s s , and soil particle volume of expansive clay V sw s as follows: (108) p c0 = p 0 exp N − v ss 0 −̃ln(p 0 p r ) −̃. (109) The relationship between the initial interlaminar distance and thickness of expansive clay mineral is Additionally, the mass ratio of nonexpansive clay to swelling clay r for soil particle densities sw and s is Using the mass ratio, the ratio of a and b is set as Finally, the initial interlaminar void is set by substituting Eqs. (109) and (110) into Eq. (112) as follows: