Modelling the behaviour of a soft rock using a kinematic hardening constitutive model

Soft rocks are of great geological and geotechnical interest, due to their mechanical and structural characteristics falling between those of a soil and those of a rock. This presents several challenges for constitutive and numerical modelling. Here we propose a calibration strategy of an advanced kinematic hardening model with structure degradation against experimental data on soft rocks, accounting for their features at the micro- and macro-scales. The experimental results show a highly structured material, with strength and stiffness being controlled by inter-granular bonding as well as by mean effective stress. The clay content of soft rock samples seems to be a dominant factor: the specimens show lower stiffness and strength when the content is higher. Accordingly, the samples have been classified into two different types (A and B). The constitutive model reproduces well both the mechanically strong and stiff behaviour of Type A as well as the more ductile mechanical response of Type B specimens observed during drained triaxial compression tests. With increasing confining pressures, the transition from brittle and dilatant to more ductile and contractant behaviour observed for both types of rocks is well captured by the adopted constitutive model.


Introduction
The application of engineering science to the study of soils and rocks requires conceptual and mathematical models, able to describe the real behaviour of these materials [1]. Soft rocks are borderline cases between soft soils and hard rocks, with a stiffer and stronger structure compared to soft soils, and are not influenced by large-scale discontinuities, unlike hard rocks [2].
Experimental observations have shown similarities in the behaviour of soft rocks and natural clays, even though their structure may have originated from different processes [1,[3][4][5][6][7]. A key point in the experimental and modelling investigation of this class of materials is the characterisation of their preand post-gross yield behaviour, where the term "gross yield" indicates a state in the effective stress space outside the elastic domain at which a significant discontinuity in stress-strain response can be easily identified [8]. Behaviour before gross yield is generally stiff, but not necessarily elastic, as structure degradation may start with stress changes within the structure locus [1]. After gross yield, the structure is progressively lost with increasing plastic strains, causing the mechanical response to converge towards that of the reconstituted material [9].
Several constitutive models have been developed specifically for soft rocks [10][11][12][13][14][15] using a single surface to represent the in-situ structure inside which the behaviour is elastic, or with an additional surface to describe the intrinsic state. As an alternative, the multi-surface constitutive model proposed by Rouainia and Muir Wood [16] has been adopted in this work to simulate the mechanical behaviour IOP Conf. Series: Earth and Environmental Science 833 (2021) 012116 IOP Publishing doi:10.1088/1755-1315/833/1/012116 2 of a specific soft rock. The model was initially developed for structured clays and successfully used in many geotechnical engineering applications [17][18][19]. The model offers the advantage in defining a small elastic bubble that moves inside the structure surface, which acts as a bounding surface [20], according to a kinematic hardening law. This implies that plastic strains may develop due to a stress change within the structure surface as a function of the distance between the current stress state and its conjugate stress on the bounding surface. Kinematic hardening models, originally proposed by Mroz [21], have also the advantage to account for the effect of recent stress history. Furthermore, the formulation of the model adopted in this work is based on the intrinsic characteristics of the material, through the definition of a so-called reference surface. Indeed, it has been demonstrated that the intrinsic properties of a natural soil can provide a robust frame of reference against which the in-situ state of the soil, its structure and subsequent degradation under loading can be effectively assessed [9].

Experimental results
The experimental results presented in this study refer to a calcareous mudstone, which has been extensively investigated at the micro-and macro-scale by Simpson [22]. These young soft rocks were likely deposited within a low energy, but chemically dynamic hypersaline environment, and are predominantly biogenic/chemical in origin. They are characterised by weak cementitious bridges that may precipitate within the voids very rapidly, introducing a certain degree of structure [5,23]. The testing program carried out to study these mudstones include the use of X-Ray Powder Diffraction (XRPD) and Scanning Electron Microscopy (SEM), along with standard soil and rock characterisation tests to gain insight into the mineralogy of the rocks and their cementing agents and to assess their physical and mechanical properties. The samples, on which the tests were performed, arose from two different site investigations of similar superficial deposits [22]. From the geological point of view, the two sites are close to one another, but the undisturbed samples for triaxial tests were taken from different boreholes, some at several kilometers of distance from each other.
Two distinct subtypes, namely Type A and Type B, have been identified within the calcareous mudstone, according to their characteristics and behaviour. XRPD analysis has shown that both subtypes consist of the same minerals but in different percentages. Type A samples have overall higher percentage fractions of Dolomite (with an average of 88.5%) and lower percentages of Palygorskite, a fibrous clay mineral (on average equal to 6.9%), whilst having similar amounts of Quartz and Halite to the Type B samples. The Type B specimens, instead, have significantly lower percentages of Dolomite (average of 78.0%), but higher amounts of Palygorskite (on average equal to 16.3%). These trends correspond well to the findings of the PSD analysis, where Type B samples have a greater percentage of clay sized grains than Type A. In terms of Atterberg limits, Type A is characterised as being of intermediate to low plasticity (plasticity index of 17%) in contrast with Type B (plasticity index of 25%).  The increase in plasticity is interpreted as a result of high concentrations of Palygorskite in Type B samples, as confirmed by the SEM analysis ( fig. 1). Straight fibers at random orientations can be observed within the void spaces of Type B samples, while in Type A specimens the grains are tightly bonded due to repeated precipitation and dissolution of Dolomite over time.
The mechanical behaviour of the calcareous mudstone has been investigated by conducting a set of four consolidated drained triaxial (CDT) tests per subtype, at increasing confining pressures ranging from 235 kPa up to 1300 kPa. Samples are indicated with the letter A or B, depending on the subtype, followed by a number from 1 to 4 for increasing cell pressures (e.g. A1 refers to the triaxial test result on Type A sample sheared at p = 235 kPa). Table 1 summarises the initial properties of the samples subjected to CDT tests. The stress-strain behaviour shows that all tests on Type A samples are characterised by higher stiffness and deviatoric peak stress than the Type B specimens, sheared at the same mean effective pressure. In both cases, there is a transition from brittle (rock-like) to ductile (soil-like) behaviour with increasing confining pressure, as also observed by Pellegrino [24] on calcarenite and tuff samples. This transition is seen to occur at significantly lower values of mean pressure for Type B samples compared to Type A. Volumetric behaviour observed during the 235, 470 and 940 kPa Type A tests is seen to shift from being initially contractant to dilatant as the shear strains increase. Only the A4 sample is noted to be purely contractant throughout the test, indicating that plastic strains may have caused a degradation of the material bonding during the test. In Type B specimens, the transition occurs starting from the test conducted at 940 kPa.

Constitutive model
In this work the kinematic hardening model with elements of bounding surface plasticity formulated by Rouainia and Muir Wood [16], is adopted. The model is characterised by three surfaces in the stress space, with an elliptical section in the p-q plane ( fig. 2). The structure surface holds information about the current magnitude and anisotropy of the structure. As plastic strains occur, the structure surface tends to collapse towards the reference surface, which represents the behaviour of the reconstituted soil. The bubble surface encloses the elastic domain and moves within the structure surface according to a kinematic hardening rule. A translational rule ensures that the elastic bubble does not intersect the structure surface. When the stress path moves towards the structure surface, the behaviour is elastic if the current stress state is within the bubble surface, while it is elasto-plastic if the current stress state lies on the bubble surface. The bubble is carried along by the stress increment according to the hardening and translation rule. The magnitude of plastic strains depends on the plastic modulus, which consists of two terms: the first is a function of the distance between the current stress state and its image on the structure surface, ensuring a smooth variation of the stiffness as the bubble approaches the structure surface, while the second term is derived from the condition of monotonic loading when the bubble and the structure surface are in contact at the current stress point. The model is also characterised by a volumetric hardening law, so that as plastic volumetric strains occur, all three surfaces change in size, as well as by a destructuration rule describing the monotonic reduction of the degree of structure as a function of both volumetric and deviatoric plastic strains. The shape of the three surfaces remains the same and can change in size due to volumetric hardening and structure degradation.

Calibration of the constitutive model
The constitutive model has been implemented in the finite element code PLAXIS [25], used in this work for the simulation of triaxial tests and the calibration of the model parameters, based on experimental results, typical values available in the literature and a sensitivity analysis. For a few parameters, the same values for both Type A and Type B samples have been assumed. These include the slope of the intrinsic normal compression line λ * (equal to 0.03), the Poisson's ratio ν (0.25), the destructuration parameter A * (0.5) and the stiffness interpolation exponent ψ (1.15). Since the available data do not give any information on a possible initial anisotropy of the material, it has been considered appropriate to set η0 to zero. The model simulates the destructuration process with a monotonic reduction with increasing plastic strains of the parameter r, which represents the degree of structure, at a rate given by the factor k. The destructuration parameter A * can take values between 0 (when the destructuration depends entirely on the plastic volumetric strains) and 1 (if it depends entirely on the plastic deviatoric strains). In this paper, A * has been chosen such that the destructuration process is dependent on both deviatoric and volumetric plastic strains, in equal proportion. The rate of destructuration k has been estimated based on a sensitivity analysis, where different values have been adopted within the same subtype, as the destructuration process is somehow faster for the Type A sample sheared at the lowest confining pressure (p = 235 kPa). The parameter k is equal to 2.2 for A1 and 0.3 for the remaining Type A samples, while it is equal to 0.03 for all Type B samples, since a slower degradation of structure seems to take place in these samples.
The values of the stiffness interpolation parameters B and ψ, adopted in the model for the definition of the plastic modulus at current stress, have also been estimated by means of sensitivity analyses. While ψ has been assumed to be the same for all samples, B is different for each subtype (33 for Type A and 15 for Type B).
The elastic behaviour is described by the slope of the swelling line in a volumetric strain-logarithmic mean stress compression plane, κ * , and by the Poisson's ratio ν, used to calculate the pressure dependent bulk modulus K and shear modulus G. The higher stiffness of Type A samples has been obtained with a κ * value of 0.002, while a value of 0.007 has been used for Type B. Based on the results plotted in the q/p vs εs chart ( fig.3) Figure 3. Experimental stress ratio curves for Type A (left) and Type B specimens (right), at increasing confining pressures.
The size of the bubble has been considered relatively large compared to the corresponding reference surface: the ratio of the axis of the elastic bubble with respect to axis of the reference surface, R, is equal to 0.5 and 0.7 for Type A and Type B, respectively. As an expression of the larger peak deviatoric stresses observed in the results of Type A samples, the initial degree of structure r0 is larger and equal to 2.3 compared to the value of 1.8 assumed for Type B specimens.
The overconsolidation ratio should not be intended as a geological OCR, but it is used to initialise the reference surface, and consequently the structure and bubble surfaces. Its value allows to calculate the intrinsic pre-consolidation pressure (2pc), i.e. the axis of the reference surface along p, as the initial mean pressure multiplied by OCR. For both subtypes, the reference surface has the smallest size in the case of the sample sheared starting from a cell pressure of 235 kPa, as the A1 specimen shows a different behaviour from the other samples. For the Type A material, a distinction has also been made between the A2 sample, retrieved from a borehole quite distant from the locations where A3 and A4 samples were retrieved. For Type B, a different OCR has been specified for B2 and B3 compared to B4. A summary of the OCR values specified for each test is reported in Table 2, together with the corresponding sizes of the three surfaces. Note that the overconsolidation ratio assumed for each test refers to an initial isotropic state of 235 kPa and that the initial void ratio is not required for the model initialisation, since its effect is taken into account through the normalisation of λ * and κ * . The representation of the surfaces in the p-q plane is shown in fig.4 for each subtype.

Results
The triaxial tests simulations consisted of an isotropic compression phase, starting from the initial cell pressure of 235 kPa up to the specific pressure reached in each experiment, followed by a straincontrolled shearing phase in drained conditions.
The stress-strain results and the volumetric behaviour of each subtype can be seen in fig.5 and 6, respectively. The results have been presented up to a deviatoric strain value preceding the onset of strain localisation in the sample, since its effects can only be captured using special numerical algorithms beyond the scope of this work. Different scales for each subtype have been used to show the volumetric strains, as their values are much larger in the case of Type B samples, characterised by a more contractant behaviour at lower confining pressures compared to Type A specimens.

Discussion and conclusions
The simulation results can be considered in good agreement with the experimental data of the soft rocks under investigation. The model can capture both the evolution from a contractant to a dilative behaviour during the tests performed at lower confining pressures in both subtypes, and the purely contractant behaviour exhibited by B3 and B4 samples. The simulations are also successful in capturing the magnitude of the observed peak deviatoric stress for both subtypes.
The model can reproduce well the stiffer behaviour of Type A in contrast to Type B samples, although a good agreement in terms of stiffness between the experiment data and the simulations is observed, especially in the case of Type B specimens. It can be noted that the shear stiffness observed in the experimental results of Type A samples does not increase with higher confining pressure, as expected, but tends to be the same for the A1 and A2 tests. Instead, the constitutive model predicts a stiffness which increases linearly with p. The use of a different formulation for the elastic behaviour, such as the small-strain stiffness equation proposed by Viggiani and Atkinson [26] or the hyperelastic formulation proposed by Houlsby et al [27], may improve the model predictions and will be considered in further developments of this research.
Concluding, soft calcareous mudstones should be thought of as highly structured/bonded/cemented soils, as they show a similar behaviour to natural clays in terms of additional strength and stiffness due to the presence of a micro-structure. Deviations in behaviour from structured soils occur in terms of higher shear stiffness at lower mean effective stress, likely caused by the presence of a well-developed bonding/cementation. Within this perspective, the work presents a first attempt to use an advanced soil constitutive model, accounting for structure degradation, to predict the mechanical response of a calcareous mudstone. The results seem to indicate that the shear behaviour of the investigated soft rock can be adequately interpreted within the elasto-plastic framework developed for cemented soils.