Failure analysis of the edge-notched beam test on fluid-exposed Berea Sandstone

Injection of fluids during wastewater disposal and geologic carbon sequestration causes induced stress and changes in rock’s elastic and failure properties, such as elastic moduli and fracture toughness. An accurate understanding of such changes in the response requires modeling and analysis of fluid-induced changes in rock’s stress and deformation states for which core-scale mechanical loading tests are often employed. Using experiments and simulations of Single Edge Notched Beam (SENB) test on water-exposed and supercritical (sc) CO2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_2$$\end{document}-exposed Berea sandstones, we quantify the effect of pore fluid on tensile crack propagation, spatial distributions of shear and tensile stresses, and fracture toughness and J-integral properties. We use Digital Image Correlation (DIC) to monitor the changes in local stress and deformation fields during crack initiation and propagation. A representative numerical simulation model of dynamic crack propagation under SENB loading is built using the extended finite element method. Results are analyzed using the stress path analysis to identify shear and tensile stress concentration regions, which provide precursory information for failure events during the test. Stress profiles along and across the growing crack are analyzed and compared with Irwin’s classical solutions to understand the impact of SENB loading. Evolution in stress intensity factor and crack length shows subcritical crack growth prior to Griffith-like failure at the structure scale.


Introduction
Long-term injection of CO 2 in deep saline aquifers for carbon sequestration [10,16,40,45,82], enhanced oil recovery [88,90], and enhanced methane recovery [36] causes substitution of in-situ pore fluids, such as water and hydrocarbons, with the supercritical (sc) CO 2 and a mixture of the scCO 2 and the ambient fluids. In general, exposure to scCO 2 alters the poromechanical response of the reservoir rock to hydromechanical loads of natural or anthropogenic origin, e.g. hydraulic fracturing [18], injection-induced stresses [88,90], earthquakes and ground deformation [56,79]. Multiple studies have investigated the effects of fluid substitution on the properties of the rock [1,5,13,15,17,25,28,34,65,75]. Water in the saline aquifer can react with the CO 2 to form acidic mixtures which can react with the minerals (calcite, dolomite, feldspar, quartz, and clay minerals like kaolinite) in the grains and in the cementation material between the grains to increase the rock's porosity [12,51,71] and decrease the rock's stiffness and strength [4,60]. It is also important to understand the mechanical effect of scCO 2 in a single fluid phase system, in the absence of significant water saturation, because scCO 2 and water have different properties and some of the CO 2 -water mixture properties, e.g., the mass density, may not follow from linear interpolation between the properties of the pure phases. The single-phase study is also relevant to understand the alteration in the nearwellbore region of many CO 2 injectors where in-situ fluids are almost completely displaced by the injected scCO 2 [23]. Gassmann's model, which assumes homogeneous isotropic mono-mineral rock saturated with chemically inert fluid under an undrained condition [14,29] provides a theoretical basis for estimating the change in elastic moduli (e.g. bulk modulus) resulting from the pore fluid substitution. Improved fluid substitution models have emerged since Gassmann's work, e.g. [44]. One study on Berea Sandstone [65] reported that, compared to dry specimens, water-saturated specimens show a higher dynamic shear modulus, and scCO 2 -saturated specimens show a lower shear modulus. Ao et al [7] reports a decrease in elastic moduli as well as in tensile and compressive strengths of their scCO 2 -treated shale specimens due to dissolution and pore pressure-induced strain. Reduction in static and dynamic elastic moduli has been reported in the presence of pore fluids such as water, kerosene, crude oil, and scCO 2 [7,44,85,86].
It is expected that fluid substitution leads to changes in rock's failure properties such as the fracture toughness [52]. Fracturing events near the injection well or the caprock-reservoir boundary, injection-induced seismic or aseismic events [24,89,90] on faults, and the resulting fluid leakage events [56,57] are important to detect during real-time storage monitoring so their impact can be controlled. The strength and toughness properties of a rock sample are related to each other due to micromechanical processes and granular texture. To estimate the rock strength at meso-scale, i.e., above the grain scale, there are mainly two kinds of measurement methods [30]: (1) static methods, which include uniaxial compression (ASTM D 2938-86), triaxial compression, Brazilian test (diametrical compression of disk-shaped rock specimens), bending test (ASTM E1820), torsion, and hollow cylinder test, and (2) dynamic methods, which include resonant bar method and ultrasonic pulse method [2,3,35]. To mimic the in-situ conditions, some of the tests are conducted on specimens saturated with a pressurized pore fluid. However, using a pressurized coreholder makes it difficult to measure the deformation field during the test using direct methods such as optical imaging. Some of the tension-based failure tests are hard to conduct safely in the lab for highly brittle specimens or for specimens with an unexpected brittle response. The Three Point Bending (3PB) test or Single Edge Notched Beam (SENB) test addresses some of these issues. SENB is a standard rock fracture toughness test that is used to estimate flexural and tensile failure properties of a rock under uniaxial compression. Several experimental studies have been conducted using the above-mentioned methods to quantify fluid-induced changes in elastic and failure properties of the rock [4,50]. The studies have been mostly complimentary. However, there have been some variations in their findings. For example, one study [67] found negligible changes in elastic properties of scCO 2 and CO 2 -water-saturated limestone specimens under in-situ conditions. Another study [80] showed a drop in the dynamic elastic modulus and strength of CO 2 -brine-saturated carbonate rocks.
These studies have raised the importance of analyzing stress distribution and evolution in fluid-saturated rocks undergoing mechanical failure, especially in standardized mechanical tests such as SENB. Induced stresses are important in the field because they dictate the timing, magnitude, and location of important geomechanical events such as fracturing, fault reactivation, and fluid leakage [56,77,89,90]. Stress distribution and stress trajectories of standardized tests, e.g., the Brazilian test method [21,48], can also inform us about fracture initiation and propagation characteristics. However, a measurement of the in-situ state of stress is difficult in practice due to the tensorial nature of the stress tensor; there are six independent stress components at each point in a 3D domain. The principal stress coordinate system also has six unknowns; three principal stresses and three principal directions, which may evolve in time with the near-wellbore flow and deformation processes and the far-field boundary conditions. Due to these issues, numerical simulation has emerged as the most promising tool for quantifying the state of stress during dynamic failure. The SENB test, which is the focus of this study, has been simulated using the microplane model [19] and the Lattice Discrete Particle Model [22,43]. However, those studies focused on reproducing the load vs. displacement curve and the crack pattern observed in the SENB experiments. A detailed analysis of the spatial distribution and temporal evolution of the SENB stress field does not exist in the literature.
Here, we present results from an integrated study that combines the experiment and simulation of SENB tests on fluid-exposed Berea Sandstone specimens with water and scCO 2 as the fluids. Digital Image Correlation (DIC) is utilized to quantify the displacement and strain fields, which are then analyzed to understand the effect of fluid type on crack propagation and toughness alteration processes. We build a high-resolution finite element model of the experiment, calibrate the model with experimentally measured displacement fields, and use the model-predicted stress field to analyze the evolution and mode of failure (tensile and shear) at different points in the domain that are representative of failure scenarios during CO 2 storage in real fields.

Specimen properties and geometry
Berea Sandstone consists of mainly quartz mineral grains bonded by silica. The composition of our specimen is summarized in Table 1. The porosity of the specimens is in the range of 18-20%, permeability is in the range of 100-200 mD, and the dry bulk density is 2.1 g/cm 3 . The dimensions of all specimens are 50 mm by 25 mm by 10 mm (with an error of ±1 mm). The cores have been extracted from thick homogeneous isotropic sandstone beds, away from any bed boundaries. The long side of the specimen is parallel to the bedding orientation. At the length scale of the specimen, the specimens can be considered structurally and mechanically isotropic and homogeneous for our testing purposes. The notch is created using a diamond saw at least 48 hour before exposing the specimens to water or scCO 2 . The notch is located at the center of the longest dimension of the specimen, is 1 mm in depth (c), and extends across the entire thickness H of the specimen (Fig. 1). The two sides of the V-notch tip come together symmetrically at an angle of 30 degrees between the two sides to form the notch tip. Specimen's properties and geometry are summarized in Table 2. The faces and sides of the specimens are flattened using a Buehler MetaServ 250 Grinder-Polisher such that the opposite sides are parallel to minimize any stress nonuniformity during loading, which can alter the failure behavior of the specimen. The fineness of the grinding pad is approximately 15 to 18 micron.

Single edge notched beam (SENB) test
The Three Point Bending or Single Edge Notched Beam test is a type of mechanical loading test designed to measure the fracture toughness of solid materials, e.g. metal, ceramic, concrete, and rocks. It is especially useful to measure tensile strength and failure properties of highly brittle materials or materials displaying predominantly brittle responses under imposed boundary conditions. This holds true for our selected Berea sandstone specimens under uniaxial unconfined loading conditions. The specimen lays on two rigid rollers while axial loading is applied parallel to an existing notch [42]. Axial loading is applied to the rock specimen by a closed-loop hydraulic compressive machine (Instron system) with a static load cell of ± 30 kilonewton, 5 kg weight, and 50 Newton-meter bolt torque. The machine applies a downward load to the specimen via a rigid cylindrical pin attached to the load cell that moves at a prescribed displacement rate. The loaddisplacement data from Instron is output by the Bluehill data collection software. An illustration of the test setup is shown in Fig. 1, which also shows a schematic of the specimen with its load-displacement configuration. The . c A schematic of the specimen with the load-displacement configuration. P is the load applied by the loading pin at the middle of the top surface, and d is the load cell displacement. The specimen rests on two support pins. c is the depth of the notch into the specimen fracture develops from the notch along the specimen width axis.

Digital image correlation
Digital Image Correlation (DIC) is a method to quantitatively measure 2D displacement and strain fields on flat surfaces of solid bodies under dynamic loads using highresolution pictures of the surface at a relatively high frame rate [61,63,68,74]. The DIC method has been used to study the strain field in 3PB tests [49]. In this technique, reference speckles are created on a specimen surface. Before loading, the position field of all the speckles provides the reference or undeformed configuration of the specimen in terms of the material or Lagrangian position vectors. As the specimen deforms under the applied load, the speckles are displaced following the deformation, which gives us the current or deformed configuration of the specimen. Correlating the deformed configuration to the reference configuration allows one to quantify the displacement vector field (U x and U y in a Cartesian 2D space) and the strain tensor field (e xx , e yy , and e xy ), which evolve in time. We create the speckles by coating the specimen surface with flat white paint. The specimens are left to dry and then sprayed with a black flat paint until a speckling pattern is achieved. Our DIC system consists of a Grasshopper camera with a resolution of 5 MegaPixel and a Nikon Micro-Nikkor 60 mm f/2.8D lens. The image frames are stored by FlyCapture SDK software at a predetermined frame rate. Post-processing is performed using the Vic-2D correlation software, which is able to extract the correlation among images taken at successive time steps using one of the three criteria [81]: Squared Difference (SD), Normalized Squared Difference (NSD), and Zero-Normalized Squared Difference (ZNSD). The NSD criterion is insensitive (yields the same correlation coefficient) to scaling up or down of the brightness. However, it is sensitive (yields a different correlation coefficient) to linear shifts in brightness. We use the Normalized Square Difference (NSD) criterion and uniform lighting conditions to obtain the correlation coefficients in our experiments.

Experimental workflow
The tests are conducted on specimens with different pore fluid types: air, water at room pressure and temperature (25 C), and scCO 2 . We perform six tests on dry specimens, seven tests on water-exposed specimens, and four tests on scCO 2 -exposed specimens. Results from the tests on the dry specimens are used as the base case to analyze the effects of water and scCO 2 as pore fluids on the geomechanical properties of our specimens. The exposure protocol for water-exposed specimens is as follows. Dry specimens are immersed in a flask full of water at room pressure and temperature. The specimens are kept in water for different exposure durations ranging from 72 h to 240 h ( Table 3). The exposure protocol for scCO 2 specimens is as follows. An autoclave containing a stainless steel pressure chamber cell is used. The inside diameter of the chamber cell is 1.25 inch, wide enough to contain the specimen without touching the cell wall. The chamber is heated up to 50 C to reach scCO 2 conditions upon gas injection. We use a 99.99% pure CO 2 Coleman Grade cylinder at 830 psi as our CO 2 source. A Teledyne Isco pump Model 260D is connected to the chamber cell to raise the CO 2 pressure inside the chamber cell. To reach supercritical conditions inside the gas reservoir of the pump, we wrap the specimen with heating tape and allow it to heat up to 50 C for at least 1 h before injection. CO 2 is injected into the chamber cell until the chamber pressure increases to 3000 psi (AE50 psi). The same pressure is used for the four CO 2 -exposed specimens. We use 3000 psi because it is a typical subsurface pressure value during CO 2 injection in many projects on CO 2 storage and enhanced oil recovery, e.g., in the Morrow-B sandstone of the Farnsworth Unit CO 2 storage project in Texas [88]. The exposure durations for the scCO 2 specimens are given in Table 3. Prior to taking the specimens out of the chamber cell, the pressure is decreased at 100 psi/min until we reach room pressure to avoid possible cracking of the specimen due to sudden depressurization. The specimens are then taken to perform the SENB test with DIC. We perform the SENB tests at a constant loadcell displacement rate of 0.07 mm/min of the loading pin, which remains in contact with the specimen top surface. This displacement rate falls in the range of American Society for Testing and Materials (ASTM) standards for this test. The loading cell is adjusted so that the load axis is aligned with the notch axis. The DIC system is placed two feet away from the Instron machine, and it faces the front surface of the specimen as shown in Fig. 1b. The camera frame rate is set to be 2 frames/second, sufficient for the Vic-2D software to post-process the images. Analysis of the post-processed images indicates the path of crack propagation before the specimens fail catastrophically. Figure 2 shows our workflow in this study.

Deformation and elastic moduli
The flexural stress is calculated using the following equation [8]: where r is the stress at the outer surface in the load span region in MPa, P is the applied force in Newton, L is the support span in mm, B is the width of specimen in mm, and H is the thickness of specimen in mm. The maximum flexural strain is calculated using the following equation [8]: where e is the maximum strain at the outer surface, and d is the load-cell displacement at the mid-span L/2 corresponding to load P. The flexural modulus is where dP dd is the slope of the tangent to the straight line portion of the load-deflection curve; see Fig. 4 in the Results section. The reasoning behind using the straight line portion is that it constitutes the linear elastic region of the mechanical response of the tested specimens. This will allow us to obtain Young's modulus E of the specimens from their flexural modulus [64]. Table 4 shows the values obtained based on the specimen's pore fluid type.

Fracture toughness
Many fracture modeling tools are based on the Irwin-Griffith brittle failure theory, according to which a crack extends when the strain energy release rate becomes equal to the surface energy required to create the two sides of the fracture segment [31,47]. When the energy available from   [47]. Irwin [41] made an important modification to the theory and proposed a critical stress intensity factor, normally referred to as the fracture toughness, which is defined as a measure of the ability of a material to resist an unstable fracture propagation. As per the Griffith-Irwin theory [41], a pre-existing crack of a certain length begins to grow rapidly when the stress intensity factor K at the crack tip exceeds the fracture toughness of the specimen, which for mode-I fracture is denoted by K Ic . K is related to the tiplocal stress tensor r ij , where ðr; hÞ denotes a point in polar coordinates near the crack tip, with origin (0, 0) at the tip and h measured from the crack axis, and g ij are trigonometric functions. For a mode-I fracture growing along the y axis [47], In the case of mechanical loading experiments where stress data is not available, the common practice is to calculate K and K Ic using the load vs. displacement data available at the point of loading, which we describe below under the energy method. In the case of a numerical model of SENB with simulated stress fields, we calculate K using the tip stress field.

Energy method
Rocks are considered quasi-brittle material due to the spatial heterogeneity created by the presence of grains and pore spaces. Failure in quasi-brittle materials is different from that in a brittle or ductile material primarily due to the presence of a relatively large fracture process zone which is comparable to the specimen size. Strain softening in the process zone leads to nonlinearity and deformation dynamics that are often negligible in metals and ceramics. To calculate the plane strain fracture toughness of our sandstone specimens, we use the energy method which is suggested by the International Union of Laboratories and Experts in Construction Materials, Systems, and Structures (RILEM). RILEM's fracture mechanics of concrete method (FMC 50) [55] uses the mechanical energy obtained from the load-displacement curve. The loading curve represents the derivative of the fracture energy [26]. Thus, the fracture energy can be calculated by integrating the load-displacement curve. To calculate the total fracture energy G F , we use the following equation: where W o is the area under the load-displacement curve shown in Fig. 4, m is the specimen mass between the supporting pins, which is calculated as m ¼ Total Specimen Mass Â Span=Specimen Length, g is the acceleration due to gravity, d o is the load-cell displacement at failure, and A f is the surface area of the fracture. The fracture surface area can be estimated from the cross-sectional area of a planar surface passing through the notch as illustrated in Fig. 1. Table 2 shows the values used in the calculations. The mode-I fracture toughness is then calculated for each specimen using the following equation:

J-integral analysis
Fracture initiation and growth processes in rocks and similar porous media are often nonlinear due to yielding at high stress, viscoelasticity and creep, sub-microscale plasticity, and intergranular interactions [37,69]. The small size of a test specimen can also invalidate some of the assumptions of LEFM. The concept of the critical stress intensity factor in Griffith's theory is limited to LEFM, which leads to crack-tip singularity and requires regularization techniques such as the idea of a process zone. For nonlinear elastic processes, which almost always precede the onset of plastic deformation [27], a popular method to characterize the tip-scale stress field is the J-integral method [66]. J-integral is a path-independent integral on a contour line around the crack tip, and it is defined as the rate of change in the mechanical energy (sum of elastic energy and cohesive energy) of the crack system per unit change in the crack length [47]. Neglecting body forces and kinetic energy, the J-integral for a crack along the ydirection is given as where i and j are dummy indices for the coordinate direction (x and y in the present case) and Einstein's summation rule is adopted, C is the contour around the tip, n ¼ ½n j is the outward unit normal on C, and d yj is the Kronecker delta function which is equal to 1 when j ¼ y.

Numerical modeling
To understand the effect of fluid type on the spatial distribution of the state of stress and how the stresses evolve over time, we construct a high-resolution simulation model of our SENB test that accounts for dynamic fracture propagation along the notch. The DIC data and load vs. displacement results discussed above cannot provide spatial distributions of the stresses, especially in high resolution. It requires a numerical model that must be representative of the SENB experiment in terms of specimen dimensions, loading conditions, and physical properties and processes that are active during the experiment. Once the model is calibrated with the experiment, we can analyze the simulated stress field in terms of the invariants and components of the stress tensor. This provides important insights into the failure processes during SENB. We build a high-resolution 3D finite element model of our specimen as shown in Fig. 3. We use Abaqus, which is a finite element simulator widely used in the computational solid mechanics community to solve dynamic equilibrium problems. We used the Extended Finite Element Method (XFEM) to model crack initiation and propagation from the notch, with a failure criterion based on the local maximum principal stress. The punch (loading pin) and support pins are modeled as discrete rigid bodies, and the specimen is modeled as a deformable and brittle material. We used a standard 3D stress element type (C3D8R) with Hex-shaped elements in a structured meshing routine. We mesh the domain using an average element size of 0.25 mm by 0.5 mm by 1.25 mm in the three Cartesian directions, respectively. To model the load, we apply a loading step with nonlinear large deformations where the initial time increment size is 1 s, and the minimum and maximum step sizes are 1 s and 180 s, respectively.
The material properties used in the 3-D numerical model are given in Table 5. The value of the maximum principal stress, which is related to the tensile strength of the fluidexposed specimen, is determined via model calibration, which is performed by comparing the load-displacement curve and the U x and U y fields between the experiment and the simulation.

Load vs. displacement
The load vs. displacement curves of the SENB experiments are plotted for dry and fluid-exposed specimens in Fig. 4.
The test data shows that the maximum flexural stress, which is the maximum stress corresponding to the peak load P o during the SENB test, is in the range of 6.11 to 7.26 MPa for the dry specimens (Fig. 5).
For the water-exposed specimens, Fig. 4 shows that the displacement values at peak load are between the respective values for scCO 2 and dry specimens. In the case of the specimens exposed to scCO 2 , the peak flexural stress values are similar to those for the dry specimens. However, the peak flexural strain values for scCO 2 specimens are the smallest. Figure 6 shows the load vs. displacement curve from the SENB simulation at a point inside the specimen  near the point of contact between the top pin and the specimen. Since a displacement boundary condition has been applied in the model, the load is calculated approximately using the cell-averaged stress values and the area of contact that increases with the loading pin displacement. We also plot the strain tensor components in Fig. 9, computed using the DIC displacement data. The DIC results show that the specimens undergo different stages of deformation and the pore fluid type exerts control on the onset and duration of those stages, which confirms our observations in the load-displacement analysis. f There are three main features of the deformation behavior: (1) the region below the loading pin and between the two support pins is displaced downward, (2) growth of the crack displaces the left side region of the notch leftward (negative x-displacement) and the right side region rightward (positive x-displacement), and (3) the two regions on the two sides of the notch rotate around the loading pin as the notch develops into a crack, causing rightward motion on the top left side of the specimen and leftward motion on the top right side. The regions outside of the two support pins have the least amount of displacement because of the zero-displacement boundary condition of these pins. The DIC results show the effect of the notch's presence in each specimen. In the case of U x , as we get closer to the notch, which has a tip position at Y ¼ 1 mm at t ¼ 0, U x Fig. 4 Load vs. displacement curves of dry (top), water-exposed (middle), and scCO 2 -exposed (bottom) specimens. The insets show horizontal displacement fields (U x ) at the highlighted time steps Fig. 5 Flexural stress and flexural strain of dry, water-exposed, and scCO 2 -exposed specimens are shown as a box plot. Dry and scCO 2 -exposed specimens display higher stress compared to the water-exposed specimens. For each fluid type, we have a range of values from the experiments; black lines indicate the minimum and maximum of the range, blue lines indicate 25 and 75 percentiles respectively, and the red line indicates the median value Fig. 6 Simulated load vs. displacement curve at a point below the contact point between the top pin and the specimen magnitude increases sharply due to stress concentration near the notch [49], which is a function of the tip radius [38]. We also notice in Fig. 7 that the discontinuity in the U x field exists close to the notch tip, which makes it the potential location of the onset of unstable crack propagation immediately before the peak load. Analysis of the U y displacement field in Fig. 8 shows that the emergence of a nascent crack decreases the stiffness and increases the compliance of the specimen as evidenced by large values of downward displacements in the region close to the notch.

Deformation and fracture propagation analysis
This behavior is further analyzed in Fig. 10 by comparing U x at different Y distances from the notch for a dry specimen (D4) and a water-exposed specimen (W4) before failure.
The line closest to the notch exhibits the largest displacement increment. As we move further away from the notch (Y ¼ 10 mm), the displacement increment shows an increasingly elastic response. The displacement gradient oU x ox j X¼notch along the Y ¼ 1:35 mm profile is approximately an order of magnitude larger compared to the gradient along the Y ¼ 10 mm profile. As a result, stress relaxation via crack propagation along a path parallel to the Y-axis is anticipated. Figure 11 shows U x and U y incremental displacements from the SENB simulation at three different timesteps during the test. We calibrate the model by achieving an acceptable agreement between the experimentally measured displacements (Figs. 7 and 8) and numerically simulated displacements (Fig. 11).
The overall displacement fields of the specimen appear consistent with the experimental results. The simulated horizontal displacement field shows blue color (U x \0) on Fig. 7 Time evolution of the deformation field U x for dry (left column), water-exposed (middle column), and scCO 2 -exposed (right column) specimens based on the DIC analysis. d is the load cell displacement that increases monotonically with time Fig. 8 Time evolution of the deformation field U y for dry (left column), water-exposed (middle column), and scCO 2 -exposed (right column) specimens based on the DIC analysis Acta Geotechnica (2023) 18:4035-4053 4043 the left and red color (U x [ 0) on the right side of the notch near the bottom edge of the specimen. This agrees with the experimental Fig. 7. Near the top edge of the specimen, the colors switch; U x [ 0 on the left side of the notch and U x \0 on the right side. This also agrees with the experiments; see the column for the water-exposed specimen in Fig. 7. The simulated vertical displacement field shows a yellow-colored region with downward displacement between the two support pins, which agrees with the experimental data in Fig. 8. We use the energy method, discussed above, to compute K Ic from the experimental data. The energy method is robust because it accounts for the entire load-displacement behavior, not just the peak load. Also, it is less sensitive to the specimen aspect ratio, which for our specimens is closer to 2:1. We used the flexural modulus E f to obtain values for Young's modulus E [64]. Young's Modulus of the water-exposed and scCO 2 -exposed specimens decreased by 11% and 5.5% compared to the dry specimen modulus, respectively. The corresponding shear modulus (G ¼ E=½2ð1 þ mÞ) of dry, water-exposed, and scCO 2 - Fig. 9 Three components of the plane strain tensor for Dry, Wet, and CO 2 -exposed specimens under SENB loading. The strain components are calculated by processing the DIC image data of the SENB experiments Fig. 10 Horizontal displacement (U x ) profiles at peak load along the horizontal direction (X-axis) for dry and water-exposed specimens at three different Y values from the bottom edge. The inset figures show the full spatial distribution of the U x displacement field in the two cases Fig. 11 Horizontal (U x ) and vertical (U y ) displacement results of a SENB test simulation. Top row: at early time, before crack initiation at the tip, the displacement fields look uniform. Middle row: Bending leads to tensile crack initiation at the tip of the notch. Bottom row: The crack has grown, and the displacement field is non-uniform exposed specimens are 113 MPa, 101 MPa, and 106 MPa, respectively. Table 3 and Fig. 12 show the fracture toughness values and their statistics for dry, water-exposed, and scCO 2 -exposed specimens. The values are within a range of 0

Stress results from the numerical model
Spatial distribution and temporal evolution of the stress components control the location and timing of micro-scale cracking events, which determine the catastrophic fracturing of the specimen during SENB loading. The cracknormal stress r xx is important because of its role in the onset of tensile failure at the notch tip. For quasi-brittle materials like our specimen, the magnitude of crack-parallel stress r yy determines whether microcracks in the process zone slip or lock under compression and how the slipping microcracks are oriented with respect to the notch direction [59]. This determines the width of the process zone and, thereby, the fracture toughness and fracture length. Given such importance of the stress solution and the difficulty in obtaining it from an experiment or the classical fracture mechanics theory, we use our high-resolution simulation model to obtain the stress components and analyze them to extract features that could be characteristic and predictive of failure during SENB loading.

Stress invariants
Since many constitutive models are expressed in terms of stress invariants, e.g., the first invariant I 1 of the stress tensor and the second invariant J 2 of the deviatoric stress tensor, we analyze the stress state in the invariant space. The first invariant is defined as is the volumetric or mean stress computed from the trace of the stress tensor. The second deviatoric invariant is defined as J 2 ¼ ð1=2Þ s : s ð Þ, where s ¼ r À r v I is the deviatoric stress tensor defined as the traceless part of the stress tensor. A commonly used quantity in those constitutive models is the von Mises stress r vm ¼ ffiffiffiffiffiffi ffi 3J 2 p . Shear and deviatoric stresses contributing to r vm lead to volume dilatancy of the specimen. When volume change during shear strain is restricted, the mean stress (r v ) increases. We analyze the stress path or trajectory of selected points within the domain to learn about the location and timing of micro-scale cracking events, which often precede the catastrophic tensile fracturing of the specimen. Understanding these precursory signals can benefit predictive modeling of rock fracturing during CO 2 injection and hydraulic fracturing. Figure 13 shows the change in the von Mises stress field and the volumetric stress field, at the same time steps as those of the displacement results shown above in Fig. 11.
To analyze how the state of stress at different points within the domain evolves toward mechanical failure or stability, we perform a stress path analysis in the invariant space of von Mises vs. volumetric stress (Fig. 14).
ffiffiffiffi ffi J 2 p controls inelastic deformation and shear-induced failure response of the specimen. I 1 controls the volumetric expansion or contraction response of the specimen. We see that the notch point evolves towards tensile failure because the change in volumetric stress at that point is tensile. Points near the top and left pins evolve towards shear failure, which is aligned with the DIC observation. We see that the stress trajectory at the left pin evolves farther than the trajectory at the top pin, which indicates that the left pin is closer to the Mohr-Coulomb failure envelope. The right pin, being symmetrically positioned, follows the same  stress path as the left pin. This analysis is valuable for identifying the failure onset time, failure location (relative to the pin positions), and failure mode (shear vs. tensile) at different points inside the domain. These observations also have practical implications for hydraulic fracturing in the field. It shows a mechanism for shear activation of natural or pre-existing fractures, especially those oriented favorably to the propagating hydraulic (tensile) fracture.

Stress components
We analyze the spatial distribution and time evolution of stress components in the vicinity of the notch, where the deformation magnitudes are the highest and most dynamic. Figures 15 and 16 show r xx and r yy fields at 12 successive time steps. See, e.g. t ¼ 75:89 and t ¼ 83:5 plots; crack propagation causes stress relaxation at the tip and increases the stress behind the tip. See t ¼ 123:5 and t ¼ 138:3 plots; crack propagation causes an increase in the stress behind the tip. This high-resolution analysis shows quantitatively how stress increase is followed by discrete crack initiation and propagation events, which are followed by either stress relaxation or stress increase depending on the crack length (nondimensionalized with the specimen width). Given the importance of the von Mises stress in determining the onset of failure at a small scale, we analyze its distribution and evolution in Fig. 17 within the same zoomin view as the normal stresses.
We see a region of low von Mises stress (blue-colored lobe at t ¼ 75:89; 83:5; 123:5) ahead of the tip, which promotes tensile cracking along the notch axis. Later, at t ¼ 170:4 and 180, two regions of high von Mises stress emerge at 45 angles from the crack path on the two sides of the crack, where shear failure may occur if a favorably oriented weak plane exists in the specimen at those locations.
An important question is how the stresses vary along and across the crack. In Fig. 18, we plot the stress profiles along the specimen length (X-axis) at selected distances from the bottom edge of the specimen. This is similar to the U x displacement profiles shown in Fig. 10, but now we do it for stress which holds information about deformation yet to occur.
It is insightful to compare the simulated stress field with Irwin's analytical stress, as shown in Fig. 19. Near the crack tip, away from the boundary effects, we expect the simulated and analytical solutions to show good agreement. As x increases and we move away from the tip, the effect of the support pins and finite dimensions of the specimen cause the simulated solution to differ from the analytical solution.

Stress intensity factor and J-integral
The mode-I toughness (K Ic ) calculated using Eq. (7) provides values at the specimen length scale. Also, these values are affected by the peak load P o , which marks the time when the specimen breaks due to an unstable Griffithlike fracture. Before this time, the experimental load vs. displacement data could be dominated by Subcritical Crack Growth (SCG) or static fatigue. SCG refers to pre-existing flaws and cracks growing quasi-statically (( 0:1 m/s) at stress intensity factor values less than the critical stress intensity factor required for unstable, instantaneous rupture [9], i.e., dc=dt [ 0 for K Iscc \K\K Ic , where K Iscc is the critical stress intensity for SCC to occur. The corrosive action of water on strained silicon-oxygen bonds can create SCG conditions [6,46]. The presence of CO 2 can modulate the SCG characteristics, e.g., the subcritical crack growth index in Charles' law relating the crack velocity with the stress intensity factor [50]. We discuss SCG in more detail in the discussion section.
Given the high resolution of our simulation model, we can investigate the early-time crack dynamics in terms of the stress field and the crack length available from the simulation. We can calculate stress intensity K near the tip using the Irwin crack-tip solutions, shown in Eqs. (4)-(5), which can be plotted along with the crack length to test the SCG hypothesis. We plot the results in Fig. 20 for a point ðr; hÞ ¼ ð0:1 mm; 0Þ on the center line of the notch where g ij ¼ 1. Based on this result, many crack propagation events indeed seem to correspond to SCG and lead to a drop in K, which implies stress relaxation and is expected from the fracture mechanics theory. Also, these K values are smaller than the experimental K Ic values calculated above. This is expected because the experimental value is at the length scale of the specimen and based on the load boundary condition P instead of the local stress tensor r ij .
In Fig. 21, we show J-integral values using the DIC displacement data. The water-exposed specimen has the lowest J-integral values for the reasons discussed above. This result is in agreement with the smallest fracture toughness values of the water-wet specimens in Fig. 12. The CO 2 -exposed specimen has the largest J-integral value at the peak load (at time = approximately 662 s), which is different from Fig. 12 showing the dry specimen's toughness as the largest of three types of specimens. The reason is that the J-integral method accounts for nonlinear tip processes such as microcracking, plastic yielding, and grain sliding at the tip. These nonlinear processes occur in  both water and CO 2 -exposed specimens, and the processes can differ for the two fluid types. We also use our simulation results to compute J-integral values using Eq. 8. We perform this computation at multiple timesteps during the quasi-static period of the crack growth process, i.e., before the failure becomes unstable. As shown in Fig. 22, the J-integral values show an increasing trend with the crack length, which suggests that the process zone may be evolving.

Discussion
We find that the total fracture energy values for the waterexposed and scCO 2 -exposed specimens are smaller compared to those of the dry specimen. Also, water exposure leads to a 16-20% drop in the Berea Sandstone's plane strain fracture toughness value compared to scCO 2 -exposed and dry specimens. There are a number of mechanisms proposed in the literature that can explain these observations: mineral dissolution and fines migration, scCO 2 adsorption on mineral surface [70], and poor consolidation due to chemo-mechanical action of the pressurized fluid. The effect of moisture in reducing the tensile strength and fracture toughness of sandstone [32,58,62,78,92] and compacted clay [33,84] is known in geotechnical engineering. Barton's classic paper [11], and references therein, suggest a reduction in the tensile and shear strengths of a crack due to strength corrosion, especially at high pore pressure, which causes a reduction in the effective intergranular compression as per Terzaghi's relation [83]: r 0 ij ¼ r ij À apd ij , where r 0 ij is the effective stress tensor (compression positive), r ij is the total stress tensor (from external loading and surrounding rock), p is the pore pressure, and a is the isotropic Biot-Willis parameter. The difference in CO 2 and water's chemical action can cause different effective stresses during CO 2 and water exposure because a decreases with an increase in the effective stress [39,72]. Specimens experiencing higher pore pressure (lower effective compression) can display weaker tensile strength due to the cracking of the cementation material and sliding and rolling of the grains.
Another mechanism by which pore fluid affects failure properties is stress corrosion cracking (SCC), i.e., water molecules reacting with silicon-oxygen bonds in quartz to convert them into mechanically weaker silicon-hydroxyl bonds held together by hydrogen bonds. Based on the Charles-Hillig theory [20], SCC refers to changes in the crack tip radius, and, thereby, stress concentration [38], due to the formation of the reaction byproduct. SCC leads to subcritical crack growth (SCG) in rocks [6,9,46], as mentioned above. Such moisture-induced weakening, or time-dependent failure, leads to crack extension and subsequent coalescence of several cracks to form a macroscopic failure surface, which results in a drop in the macroscopic toughness [73]. SCG can be important for near-wellbore permeability and caprock integrity during CO 2 sequestration because the pore fluids at reservoir temperature and injection-induced stress levels can lead to stress corrosion cracking and dissolution at crack tips.
Khazanehdari and Sothcott [44] finds that viscosity and molecular polarity of the pore fluid lead to hardening and softening of the porous specimens, which show an increase and decrease in bulk and shear moduli, especially at low values of the effective compressive stress (r 0 v ) when the cracks and voids are open and provide the additional surface area for the fluid to interact with the solid. Clay content, dry pore surface area, differential stress, loading rate, fluid compressibility, and pH in the presence of background formation water are other important parameters that affect rock's mechanical properties. For example, water, which has a higher polarity Fig. 17 Von Mises stress at the same timesteps as above than hydrocarbons such as kerosene and crude oil, leads to a decrease in the shear modulus, whereas hydrocarbons lead to an increase in the shear modulus at low differential stresses. CO 2 is not dipolar like water (CO 2 is quadrupolar), and it is not organic like hydrocarbons. Because of these differences, the effect of CO 2 on rock mechanical properties is expected to be different from that of either water or hydrocarbons.
Our Berea Sandstone specimens have less than 5% clay minerals. However, even a relatively small fraction of clay, in the presence of pore fluids, can lead to noticeable changes in the host rock's mechanical properties because clays have a tremendous surface area that is electrically charged and allows diffusion of the fluid such as scCO 2 [65]. The effect of clay on the host rock's properties also depends on whether the clay is present as laminae, structural clasts, or dispersed as grain coating [87].

Conclusions
We demonstrated the potential of integrating experimental and simulation methods to obtain spatial distribution and temporal evolution of the stress and deformation fields during Single Edge Notch Bending (SENB) tests. We examined the failure processes and properties of Berea Sandstone specimens exposed to different pore fluids-water and CO 2 -separately. We built a numerical model of our SENB experiment with dynamic crack propagation to obtain the stress field and its invariants, which have not been analyzed before in such high resolution. Using the stress intensity factor vs. crack length analysis, we find evidence for subcritical crack growth in fluid-saturated specimens. Using the stress path analysis, we show how two different stressed regions-one ahead of the growing crack and the other along the line joining the loading pin and a support pin-are subject to different modes of failure. To account for nonlinear crack tip processes, we computed J-integral values by post-processing the DIC displacement data. The water-exposed specimen showed the lowest Jintegral values, whereas CO 2 exposure leads to an increase in the J-integral values, which suggests a decrease in brittleness, which is aligned with observations at CO 2 storage sites such as Decatur, Illinois and Farnsworth, Texas.
The stress state can be used to predict the timing, location, magnitude, and mode of failure events and their impact on hydraulic/mechanical properties [76,77]. The integrated analysis can be used to estimate changes in rock's fracturing potential, which is important for modeling near-wellbore hydraulic fracturing [18,53,54], injectioninduced caprock fracturing [91], and leakage of fluids [56] during hydrocarbon recovery and CO 2 storage.
Author contributions BJ designed the study. RD and BCJ conducted the experiments. RD and BJ conducted the simulations. RD and BJ analyzed the results. BJ and RD wrote the manuscript.
Funding Open access funding provided by SCELC, Statewide California Electronic Library Consortium.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright Fig. 20 Evolution of the stress intensity factor (SIF) K I ahead of the crack tip as the crack length evolves. Crack propagation events are marked by a drop in K I ahead of the tip. The crack length is measured along the loading axis straight line. Crack lengths measured along the tortuous crack path (not shown) are approximately 1-5% longer Fig. 21 J-integral values for dry, water-exposed (wet), and CO 2exposed specimens during SENB experiments. The DIC data is used to compute the J-integral values Fig. 22 Evolution of the J-integral in the SENB simulation during the quasi-static crack growth period. The solid line is the best-fit line and is shown as a visual aid holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.

Data availability
The data used in this study are provided in the manuscript and at this URL: https://figshare.com/s/ b9ddf72b6c5abcc14d38