## 2 Microstructure of Laminated Shale Specim and Micro-Constitutive Models
### Micro-Constitutive Model for Rock Matrix  
The simulation in this study is based on the interaction among adjacent discs under uniaxial compression. The motion of a disc follows the [PERSON]'s second law for the translation of mass center and the rotation around the mass center. The force-displacement law among adjacent discs is based on the contact mechanics. The contacts among the discs in a shale matrix are described through the parallel bond model (PBM), which is suitable for simulating the rock deformation and failure process [37]. Figure 2 shows the schematic diagram of this contact. Before the bond is broken, the force and moment at each contact will be updated through the following force-displacement relationship among adjacent discs:  
\[\overline{F}_{n}^{pb}=\overline{k}_{n}^{pb}A\ \Delta U_{n} \tag{1}\]  
\[\overline{F}_{s}^{pb}=\overline{k}_{s}^{pb}A\ \Delta U_{s} \tag{2}\]  
\[\overline{M}^{pb}=EIL^{-1}\theta_{b} \tag{3}\]  
where \(\overline{F}_{n}^{pb}\), \(\overline{F}_{s}^{pb}\) are the normal and shear components of contact force and \(\overline{M}^{pb}\) is the parallel-bond moment; \(\overline{k}_{n}^{pb}\), \(\overline{k}_{s}^{pb}\) represent the normal and shear components of contact stiffness, and the values of these two parameters are decided by the Young moduli of contact, the sum of adjacent disc radius and the stiffness ratio listed in Table 1. \(A\) is the cross-sectional area of contact, \(L\) is the length of elastic symmetrical beam, which is equal to the sum of the radii between two adjacent discs. and \(I\) is the moment of inertia of the cross-sectional area about the neutral axis of the beam. \(\Delta U_{n}\), \(\Delta U_{s}\) are the displacement increments in the normal and shear directions, respectively. They correspond to the extrusion/detachment and slip between adjacent discs. \(\theta_{b}\) is the relative bend-rotation increment, which corresponds to the rotation of discs.  
Figure 2: Force-displacement model for the interactions in rock matrix and beddings: (**a**) normal force versus normal displacement relationship; (**b**) tensile failure in matrix; (**c**) shear force versus shear displacement relationship; (**d**) shear failure in matrix (refer to Potyondy and Cundall [37]).  
Each contact may have two failure modes: tensile failure and shear failure. The failure is determined by the stress at the contact point. For the tensile failure, the tensile stress \(\overline{\sigma}\) is calculated by  
\[\overline{\sigma}=\frac{\overline{F}_{n}^{pb}}{A}+\beta\frac{\overline{M}^{pb} \overline{R}}{I} \tag{4}\]  
For the shear failure, the shear stress \(\overline{\tau}\) is calculated by  
\[\overline{\tau}=\frac{\overline{F}_{s}^{pb}}{A} \tag{5}\]  
where \(\overline{R}\) is the contact length in PFC2D, and \(\beta\) is the moment-contribution factor.  
A micro-crack is generated when the contact stress exceeds the strength parameters of the bonded contact. This implies that two failure modes may occur at each contact. For example, a micro-tensile crack of rock matrix (Pbt) occurs if the tensile stress \(\overline{\sigma}\) exceeds the tensile strength of rock matrix \(\overline{\sigma}_{c}^{pb}\). This tensile crack detaches the contact, forming an open micro-crack. A micro-shear crack of rock matrix (Pbs) occurs if the shear stress \(\overline{\tau}\) exceeds the shear strength of rock matrix \(\overline{\tau}_{c}^{pb}\). When shear failure occurs, the parallel-bond moment \(\overline{M}^{pb}\) and cohesion at the contact disappear and frictional sliding occurs.