Triaxial Deformation of the Goldwyer Gas Shale at In Situ Stress Conditions—Part II: Viscoelastic Creep/Relaxation and Frictional Failure

To understand the geomechanical implications of long-term creep (time-dependent deformation) response of gas shale, short-duration creep was recorded from laboratory triaxial tests on ten Goldwyer gas shale samples in the onshore Canning Basin at in situ stress conditions under constant differential axial stress. A simple power-law function captures primary creep behaviour involving elastic compliance constant B and time-dependent factor n. Experimental creep data revealed larger axial creep strain in clay and organic-rich rocks, than those dominated by carbonates. Anisotropic nature of creep was observed depending upon the direction of constant axial stress application (perpendicular or parallel to the bedding plane). Upon the application of linear viscoelastic theory on laboratory creep fitting coefficients, differential horizontal stress accumulation over a geological time scale was estimated from the viscoelastic stress relaxation concept. Further, this model was used to derive lithology-dependent least principal stress (Shmin) magnitude at depth for two vertical wells intersecting the Goldwyer gas shale formations. This newly proposed Shmin model was found to have a profound influence on designing hydraulic fracture simulation. Further, pore size distribution and specific surface area value SN2 were derived from low-pressure gas adsorption experiments. These physical properties along with weak mineral components were linked with creep constitutive parameters to understand the physical mechanisms of creep. A strong correlation was noted between SN2 and creep parameters n and B. Finally, an attempt was made to investigate how gas shale composition and failure frictional properties can influence shear fracturing. Bedding perpendicular gas shales show a higher creep magnitude than their bedding parallel orientation i.e., creep deformation is anisotropic in nature. A viscoelastic stress relaxation approach can explain lithology-dependent least principal stress Shmin profile estimated from 1-D power-law creep constitutive parameters (B—elastic compliance constant, n—time-dependent factor) and with a constrain of relative variation of in situ stress magnitudes ϕ at subsurface depth. Specific surface area SN2 is a first-order proxy for continuous creep constitutive parameters prediction as well as estimation of frictional failure properties of gas shales intersecting different stratigraphic layers. Bedding perpendicular gas shales show a higher creep magnitude than their bedding parallel orientation i.e., creep deformation is anisotropic in nature. A viscoelastic stress relaxation approach can explain lithology-dependent least principal stress Shmin profile estimated from 1-D power-law creep constitutive parameters (B—elastic compliance constant, n—time-dependent factor) and with a constrain of relative variation of in situ stress magnitudes ϕ at subsurface depth. Specific surface area SN2 is a first-order proxy for continuous creep constitutive parameters prediction as well as estimation of frictional failure properties of gas shales intersecting different stratigraphic layers.


List of Symbols
Tectonic strain in the direction of least principal stress ε H Tectonic strain in the direction of maximum principal stress

Introduction
Hydraulic fracturing of tight gas shale reservoirs is required to increase permeability and thereby improve gas flow toward the production well.Organic-rich depth intervals with the lowest least principal stress magnitude are generally the most prospective and where fracturing operations are most effective in enhancing permeability and production while minimising injection energy.However, the rapid postfracturing decline in permeability and production is often observed, within 1-2 years (Hakso and Zoback 2019;Patzek et al. 2013).Therefore, economic gas extraction from such formations requires multiple fracturing stages along a given well to maintain production over time.
The root causes of the observed gas production decline in gas shales are the subject of ongoing research.One often invoked physical mechanism is the time-dependent fracture closure associated with shale creep deformation (Gupta and Mishra 2021;Herrmann et al. 2020;Liu et al. 2016;Rassouli and Zoback 2018;Rybacki et al. 2017;Sone and Zoback 2014a;Yang and Zoback 2016).It is, therefore, necessary to quantify the inelastic behaviour of the shale (i.e., creep and relaxation) to better understand its present-day properties and more robustly predict its elastic and mechanical behaviour for long-term production monitoring (using well logging and seismic monitoring), in situ stress magnitude estimation, and wellbore stability analysis (influenced by depletion and subsidence) (Hagin and Zoback 2004a, b;Sone and Zoback 2014a).
Gas shale reservoirs are often heterogeneous at multiple scales, from the nanoscale to the macroscale (i.e., formation scale).Their heterogeneity is governed by their microstructure, depositional history, diagenesis, the presence of organic matter and its maturation, and the content in clay minerals and their texture (Kohli and Zoback 2013;Sarout and Guéguen 2008a, b;Vernik and Nur 1992).Mineralogy, texture and fabric, pore size distribution, and pre-existing (natural) fractures/faults control the overall deformation characteristics and frictional properties of gas shales, therefore indirectly affecting the prevailing stress state and faulting regime at depth (Delle Piane et al. 2015;Kohli and Zoback 2013;Mandal et al. 2020b;Sone andZoback 2013a, b, 2014a;Yuan et al. 2019).
At in situ conditions, gas shales often exhibit lower creep rates compared to other sedimentary rocks, e.g., uncemented sands, or immature shales (Das and Zoback 2011;Hagin and Zoback 2004a, b;Herrmann et al. 2020;Rybacki et al. 2017;Sone and Zoback 2014a), although clay-rich gas shales tend to display higher elastic anisotropy and creep rates compared to other types of clastic shales (Sone and Zoback 2013a, b;Vernik and Nur 1992).However, robust estimation of the time-dependent deformation responses of shales is essential for (i) the estimation of the present-day (native, pre-fracturing) state of stress at depth resulting from viscoelastic strain accumulation over geological time scales; and (ii) the prediction of the post-fracturing behaviour of the reservoir in terms of changes in permeability, production, and stress state at depth, e.g., time-dependent post-stimulation fracture closure.Sone and Zoback (2014a) inferred creep constitutive parameters from primary creep data on samples from multiple gas shale formations in the USA, subjected to triaxial stresses representative of the native stress state at depth.They successfully used these creep parameters and the theory of viscoelasticity to estimate strain accumulation and the reciprocal stress relaxation taking place over geological timescales and causing the present-day (native) state of stress.Moreover, they successfully predicted the depth profile of accumulated differential horizontal stress and the resulting least principal stress magnitude with changing lithology in a vertical well of the Barnett shale formation.Yang and Zoback (2016), Rassouli and Zoback (2018), and Xu et al. (2019) followed analogous steps and reached similar conclusions for other gas shale formations in the USA.However, Rybacki et al. (2017) and Herrmann et al. (2020) mentioned that at elevated temperature and pressure conditions, primary creep data may not be sufficient to explain the observed production decline, and the presumed fracture closure and proppant embedment.They used the European Bowland and Posidonia shales to reach such conclusions.
Previous studies show that tectonic stress accumulation may not be governed solely by the elastic stiffness of the rock [e.g., Paris basin (Gunzburger and Cornet (2007), Barnett shale (Sone and Zoback 2014b;Yang et al. 2015)].This conclusion was based on the analysis of the variations with a depth of the least principal stress magnitude and on the inferred growth of hydraulic fractures beyond the intended depth interval.These facts (and intuition) suggest that the elastic properties of the rock formation are not sufficient to accurately predict the present-day stress state resulting from gravity and tectonic loadings over depositional/geological time scales.In other words, present-day in situ differential stress (S Hmax -S hmin ) in the horizontal plane results from tectonic loading over depositional (geological) time scales, which involves both the elastic (instantaneous or short-term) and the viscous stress and strain (time-dependent or longterm) response of the shale.Our approach here goes one step beyond elasticity by accounting for viscous creep/relaxation using laboratory data acquired over approximately 6 h (primary creep stage).In doing so, we neglect complex depositional processes such as burial, diagenesis, and organic matter maturation.
Successful applications of the theory of viscoelasticity to gas shale formations indicate that the rate of primary creep is essentially governed by the fraction of clay minerals (clay content), the direction of applied differential stress with respect to the shale bedding, the effect of water (i.e., saturation and salinity), and the specific microstructure of the shale (Sone andZoback 2014a, 2014b;Xu et al. 2019;Yang and Zoback 2016).However, quantitative assessments of the impact of these parameters remain scarce.Here, we contribute to bridging this gap by analysing the impact of Specific Surface Area (SSA) and Pore Size Distribution (PSD) on creep characteristics and frictional failure behaviour of gas shales based on laboratory data available for a suite of gas shales worldwide, i.e., Barnett, Haynesville, Eagle Ford, and Fort St. John in the USA (Sone 2012); Goldwyer in Australia (Mandal et al. 2021a).
The Goldwyer formation is a prospective shale gas resource located in the Canning basin, Western Australia.The mechanical and elastic (static and dynamic) properties of representative specimens of this formation recovered from a vertical well have been recently reported by (Mandal et al. 2022).We have reported creep data (primary stage) and frictional failure characteristics of 15 Goldwyer shale samples subjected to realistic triaxial stress conditions, i.e., along and across the bedding, at room temperature, and in dry conditions.In addition, complementary mineralogical, petrophysical, and petrological data were collected for these samples, i.e., X-ray diffraction (XRD) on powdered samples, Scanning Electron Microscopy (SEM) on thin sections, low-pressure N 2 gas adsorption (LPNA) on powder samples with grain size smaller than 450 µm, and Rock-Eval Pyrolysis on crushed chunks/powder sample.The strain evolution with time during creep is fitted with a power-law of time assuming linear viscoelasticity.A statistical analysis of the correlation between the recovered creep parameters and the mineralogical, petrophysical, and petrological properties of the studied shales is conducted.We further analysed the role of specific surface area and pore size distribution as measurable quantities representative of the microstructure on the time-dependent deformation behaviour of gas shale.Given the experimental creep data, field data from wireline logging operation, and correlation functions between laboratory-measured physical properties of rock, a new workflow is proposed to predict in situ S hmin profile at depth.For this purpose, the measured creep rates are used to derive the reciprocal stress relaxation rates and then integrate them over geological timescales to calculate differential horizontal stress along two vertical wells crossing multiple lithological layers within and beyond the Goldwyer formation.Thereafter, vertical stresses and differential horizontal stresses are combined with an assumption of faulting kinematics of the stress regime to build stress profiles.Finally, we investigated sample compositions, specific surface area, and frictional failure properties to obtain insights into shear slip mechanism of gas shale reservoirs.Mandal et al. (2022) provide a detailed geological setting of the Goldwyer formation.The following can be summarised from that study:  Mandal et al. (2022a) reported the composition of the various samples of the Goldwyer formation available for this study (Fig. 2).The key aspects of interest here can be summarised as follows:

Goldwyer Shale Samples
(i) Core plugs were extracted from 10 depth intervals of Theia-1 well covering G-I and G-III units.Total of 15 samples were prepared in a cylindrical shape with L:D ratio of 2:1.Ten samples are from the bedding perpendicular direction while the remaining five samples are from bedding parallel direction.Tables 1  and 2 of part-I study by Mandal et al. (2022) can be accessed to obtain information about sampling depth, sample orientation, petrophysical properties, mineral composition, maturity, and organic content.(ii) Mineralogical compositions of studied samples were grouped based on their mechanical strength into strong, intermediately strong, and weak phases, which are displayed in a ternary diagram (Fig. 2a).Compositions of clay material from XRD analysis (Table 1) revealed that illite is the major clay type, which is also non-swelling.Samples are grouped in two base categories-sub-group 1: Th3, Th6, Th7, Th9, Th10 are dominated by clay and organic content, and subgroup 2: Th1, Th2, Th4, Th5, Th8 are dominated by carbonate minerals.These subgroups will be used to report analysis and interpretation of the creep behaviour of currently studied samples whenever required.(iii) The correlation between bulk density and weak phase is rather trivial since the weak phase includes pores (see Fig. 2b).

Low-Pressure N 2 Gas Adsorption Experiment
The low-pressure nitrogen gas adsorption (LPNA) technique at room pressure and a temperature of 77 K were used to characterize the pore space in the available suite of Goldwyer shale specimens (Mandal et al. 2021a).Applied to degassed powder samples of shale, this technique yields the mesopore volume, the pore size distribution (PSD), and the specific surface area value (S N2 ) of pores in the size range of 2-100 nm (Kuila 2013).A Micromeritics (R) TriStar II PLUS instrument (Fig. 3b) is available to conduct these measurements on powdered specimens prepared following the ISO 9277:2010 standard, i.e., rock samples are crushed into powder with grain size less than 450 µm (< 40 µm mesh size).Prior to the measurement, 1 g of powder is degassed and dried at 105 °C over 12 h (degas setup in Fig. 3a) to remove adsorbed, capillary, and clay-bound water without modifying the internal clay structures (Kamruzzaman et al. 2019;Kuila and Prasad 2013;Yuan et al. 2019).Within the apparatus, the dry powdered shale sample is cooled down to 77 K, and the volume of N 2 gas adsorbed by the specimen is measured for a range of relative equilibrium pressures P/P 0 (0.1-0.99) along the corresponding isotherm, where P and P 0 are the absolute equilibrium pressure and the nitrogen condensation pressure, respectively (Zou et al. 2018).
In practice, the pressure in the test chamber is increased until the condensation peak pressure P 0 is reached.It is then gradually reduced to generate a desorption isotherm.Classification of isotherms (Type I-IV) according to the International Union of Pure and Applied Chemistry (IUPAC)  definition is described in great length by Rouquerol et al. (1998).For a purely mesoporous material, a hysteresis loop is expected (Kuila and Prasad 2013), which is associated with capillary condensation and evaporation mechanism occurring within mesopores (equivalent to Type IV isotherm).An important feature occasionally observed in the recorded hysteretic sorption/desorption curves is the forced closure (discontinuities) of the desorption curve at a P/P 0 value of ~ 0.4 (Iqbal et al. 2021;Kamruzzaman et al. 2019;Kuila and Prasad 2013) which is due to the instability of the hemispherical meniscus of nitrogen vapor during capillary evaporation in pores.These closure features are indicative of the presence of pores with a size < 4 nm (Kuila 2013).
The SSA and PSD are determined by inverting adsorption isotherm using the BET (Brunauer, Emmett, and Teller) theoretical model and BJH (Barrett, Joyner, and Halenda) model, respectively (refer to following papers to know more about the underlying principals and assumptions of BET and BJH inversion techniques-Brunauer et al. 1938;Rouquerol et al. 1998;Sing et al. 1985).The built-in software platform within Micromeritics (R) TriStar II PLUS instrument inverts the adsorbed isotherm data to supply SSA and PSD.
The PSD is derived from the inversion of the LPNA adsorption curve by assuming cylindrical pore geometry.Based on the model's assumption, pore size refers to the diameter of an equivalent cylindrical pore in this study.It is necessary to account for the dependency of adsorbed layer thickness on P/P 0 when inverting adsorbed data to obtain pore volumes.The algorithm inverts measured isotherm data using Kelvin's equation.The technique cannot evaluate pore diameter ranges below 2 nm since Kelvin's formula is invalid in micropores (Kuila and Prasad 2013).Since PSD is confined within mesopores range as per IUPAC definition (Micropore: pores < 2 nm; mesopores: pore between 2 and 50 nm; macropores: pore > 50 nm), it will be interesting to analyse what percentage of mesopores contribute to total porosity of the examined gas shale samples.
The semi-quantitative assessment of SSA derived from the LPNA corresponds to the amount of nitrogen molecules required to cover the specimen's mineral surfaces, including its external surfaces, and the walls of the macro-to mesopores (Kuila 2013;Zou 2019).

Creep Data Recording Under Triaxial Stress
Creep tests were carried out on cylindrical shale samples during the multistage triaxial tests (MSTs) using CSIRO's Autonomous Triaxial cell (Fig. 4a) at confining pressure (Table 2) consistent with in situ stress at the samples' native depth (Mandal et al. 2022a).The available samples were unpreserved (dry) since their recovery from the Theia-1 well.The MSTs were conducted at room temperature and in dry/drained conditions.No re-saturation of the samples or control of the pore pressure was attempted to (i) avoid shale swelling/shrinking associated with clay hydration and chemoporoelastic effects; (ii) avoid pore pressure build-up and (extremely) slow equilibration/ drainage times during triaxial loading; and (iii) focus on gas shales, which are often only partially saturated with formation water (as opposed to fully saturated sealing shales/caprocks).Following the recommendations of the International Society of Rock Mechanics (ISRM), cylindrical shale specimens were prepared with a length-todiameter ratio of two and flat and parallel end faces.They were cored either parallel (horizontal, 25 mm in diameter) or perpendicular (vertical, 38 mm in diameter) to the visible bedding plane (Fig. 4b).
Creep deformation was recorded during the fourth loading stage of each MST (Fig. 5), while the shale sample was held under a pre-defined and constant triaxial stress, i.e., differential stress of 30-40 MPa, depending on the native depth of each sample, which corresponds to the addition of mean effective stress p c defined as p c = (S v + S hmin + S Hmax )/3 − P p calculated from a 1-D geomechanical model of Theia-1 (Table 2) (Mandal et al. 2020b) and 40% value of peak stress.This triaxial state of stress was achieved by increasing the differential axial stress at a rate of 2 MPa/ min while maintaining the in situ confining pressure constant.In all creep tests, we observe that most of the strain occurs instantaneously upon application of the differential stress (elastic response).Once the target differential stress is achieved, it is maintained constant, while creep strains (axial and radial) are recorded for about 6 h (primary creep stage).The limitations of creep deformation are covered in great length by Mandal et al. (2021a).
As observed in Fig. 5b, the radial strain remains virtually unaffected during the axial creep of the sample over the six hours of monitoring.In other words, time-dependent deformation occurs in a direction orthogonal to the maximum principal stress direction, i.e., triaxial loading direction.This is observed for all tested samples, regardless of their orientation with respect to bedding, i.e., vertical, and horizontal.
For a few samples, we also conducted a stress relaxation test at the pre-defined triaxial stress, i.e., time-dependent stress change at zero displacement/constant strain.The aim is to assess the reciprocity between creep and relaxation parameters and validate the use of linear viscoelasticity for our shale samples.The relaxation test is achieved in a similar way to the creep test, except that once the target differential stress is achieved, the loading actuator is halted and locked in place so that the axial displacement is zero, while the evolution of the axial stress is recorded.
At the last (fifth) loading stage of each MST (confining pressure stage-1.25 × p c as outlined in Table 2), the sample is brought to failure (peak stress), and beyond to record the post-failure response up to a total axial strain of about 4% (Fig. 5a).The combined stress-strain data from the last stage of each MST was used to characterize the frictional failure properties of each sample, i.e., shear stress, effective normal stress, residual strength, and sliding friction coefficient µ s (covered in Appendix 3).Length to diameter ratio is kept under a 2:1 ratio per ISRM standard

Viscoelastic Creep and Stress Relaxation
Any material that exhibits time-dependent deformation with no irreversible strain is called viscoelastic, i.e., elastic characteristics with a viscosity factor.Since viscous materials behave in different ways, there are multiple idealized models, such as Maxwell, Voight, Standard linear, Burgers solid, etc.During viscoelastic deformation, viscous molecular rearrangement takes place in the material and gives rise to hysteresis in the stress-strain response upon cyclic loading.The extent of hysteresis reflects the amount of energy dissipated in this process.Stress and strain data from the Goldwyer shale samples obtained during quasistatic loading-unloading cycles around the in situ stress conditions are reasonably linear, with a mild to moderate hysteresis and relatively small irreversible strain (Fig. 5a).We, therefore, neglect the plastic contribution to the total recorded deformation and assume a viscoelastic behaviour of the Goldwyer shale (Hagin andZoback 2004a, 2004b;Jaeger et al. 2007;Sone and Zoback 2014a).Sone and Zoback (2014a) proved that idealized springs and dashpots model (Maxwell or Burgers) are not suitable for describing long-term deformation from shorter duration creep data from gas shales and unconsolidated sands.This is because neither of the studied samples contains a strain asymptote nor reaches a stable strain rate after some characteristic time constant.The same is true when we fitted Goldwyer gas shales with Burgers, Voight, and standard linear solid models.Power-law model, on the other hand, can explain creep deformation behaviour at a constant deceleration rate over a longer time duration.
The theory of viscoelasticity is covered in great length in Appendix 1.For Goldwyer shale samples, we determined the creep compliance function J(t) from the recorded creep data, then compute the corresponding stress relaxation modulus function E(t).Accounting for Eq.18 in Appendix 1, we subtract the instantaneous elastic strain achieved immediately after the application of the differential stress load from the recorded strain data and focus on the subsequent time-dependent strain recorded at constant stress (see Fig. 5).Dividing the record creep strain by the applied stress step yields the require creep strain compliance of the material J(t) = creep (t)∕ o (Sone  2012; Sone and Zoback 2014a; Xu et al. 2019; Yang and  Zoback 2016).
Several authors have shown that a simple power-law function of time is suitable to accurately describe the short-term, time-dependent deformation (primary creep) of sedimentary rocks in general and US gas shales in particular (Rassouli and Zoback 2018;Sone and Zoback 2014a;Xu et al. 2019;Yang and Zoback 2016).Following this approach for the Goldwyer shale samples, we postulate that the functional form of the creep strain compliance J(t) is a power-law of time, i.e., where B and n are fitting parameters often referred to as creep constitutive parameters, i.e., n is the time-dependent exponent, and B is the instantaneous elastic compliance in response to a unit stress-step loading (at t = 0 + ).The timedependent exponent n reflects the rate of creep, i.e., the rate at which strain accumulates over time after elastic deformation.For a given unit stress-step loading, the additional time-dependent strain compared to the instantaneous elastic strain can be written as where higher values of n reflect greater creep, while n = 0 means pure elastic and time-independent deformation.
Using the viscoelastic reciprocity principle (Refer to Appendix 1 for Eq.21), the Laplace transform e(s) of stress relaxation modulus function E(t) reads where n can be a complex number with a real part greater than − 1, and Γ is the Gamma function, i.e., an extension of (1) the factorial function to non-integer complex numbers.The Gamma function assumes the value of the factorial function for positive integers.
Applying the inverse Laplace transform to Eq. 3, and using Euler's reflection formula and the basic property of the Gamma function, yields E(t) in the time domain where n should not be an integer.Noting that lim X→0 sin X X = 1 , the stress relaxation modulus/function

E(t) can be approximated by
if n is a small (|n| ≪ 1) , non-integer, complex number with a real part greater than − 1.In this limit and considering the response of the rock in both the creep and relaxation formulations of viscoelasticity (Eq. 1 and approximation to Eq. 5), the instantaneous elastic response of the material implies that static Young's modulus E ≈ 1∕B.
In practice, we assess a posteriori the applicability of the approximate reciprocity relationship for our Goldwyer shale samples by checking whether: 1.The error associated with the assumption n ≪ 1 is small enough for the approximate reciprocity relationship to be usable 2. The laboratory data acquired at the target in situ stress support the equality E ≈ 1∕B for elastic deformation based on independent data sets, i.e., static Young's modulus E extracted from triaxial stress-strain data, and elastic compliance B inverted from subsequent creep data 3.The consistency between the viscoelastic parameters B and n inverted from either creep or relaxation test data independently.

Experimental Results
This section reports the data resulting from the low-pressure N 2 gas adsorption (LPNA) tests and the creep/relaxation tests of the Goldwyer shale samples, while an analysis of frictional failure behaviour of these samples is presented in Appendix 3. (4)

Pore Network Characterization
The N 2 sorption isotherm curves obtained for our suite Goldwyer shale samples are shown in Fig. 6 as the volume of adsorbed nitrogen per unit mass of test material versus the relative pressure P/P 0 in the test chamber.Overall, the clay and organic-rich samples with the highest TOC (subgroup 1: samples Th7, Th9, and Th10) (Fig. 6b) exhibit larger amounts of adsorbed gas than the carbonate-rich samples (sub-group 2: samples Th4, Th5, and Th8) (Fig. 6a).We observe a so-called forced-closure discontinuity mostly in the desorption isotherms of the clay and organic-rich samples at around 0.4 P/P 0 (Fig. 6b).This behaviour suggests that in these samples, the volume of pores < 4 nm is significant.At high relative pressures (P/P 0 ) in the range of 0.98-1, the sorption curves of clay and organic-rich samples indicate a higher volume of larger mesopores.
The overall outcomes of the LPNA technique applied to our Goldwyer shale samples are shown in Fig. 7 and summarised in Table 3 in terms of total porosity, pore size distribution, and specific surface area S N2 , where the pore diameter D is comprised in the range 2-100 nm.Total porosity is computed from the density porosity equation with grain density coming from laboratory helium gas pycnometer measurement, while bulk density from the volume-mass relationship of cylindrical sample and with an assumed fluid density of 1 g/cm 3 .
The pore size distribution is defined as the pore volume per unit mass of material versus (the decimal logarithm of) the pore diameter D (Fig. 7a).Overall, the pores exhibit a bi-modal size distribution with two noticeable peaks: (i) smaller mesopores narrowly centred around 2-3 nm; and larger mesopores more broadly centred around 20-30 nm.The total volume of pores in the investigated range of 2-100 nm is essentially composed of the larger type of mesopores (20-30 nm).As expected, we also observe higher overall porosity values and higher fractions of the smaller mesopores in the organic-rich samples (Th7, Th9, and Th10) compared to the organic-lean samples (Th4, Th5, and Th8).In Fig. 7b, the correlation of total porosity and mesopore volume indicates a linear positive trend.It is evaluated from Table 3 that on average more than 50% contribution in total pore space comes from interpreted mesopores in the Goldwyer gas shale formation.

Viscoelastic Response
Figure 8a shows a representative set of creep curves recorded for vertical and horizontal samples of the Goldwyer shale for six hours after the application of an average differential stress of 40% of peak stress at a constant confining pressure in the range of 14-17 MPa replicating the conditions at depth (Table 2).
Laboratory creep strain data were fitted with a linear regression function in the log 10 (t)-log 10 (J) space, corresponding to a power-law of time with the compliance B (y-intercept) and time-dependent exponent n (slope) as the constitutive parameters of creep to be identified (see example in Fig. 8b).The resulting creep parameters for all Goldwyer shale samples are shown in Fig. 9 and summarised in Table 4. Figure 9a plots B against n, for both vertical and horizontal samples, color-coded with the weak fraction, i.e., clay content, total organic carbon (TOC), and porosity fractions summed.Figure 9b plots the comparison (cross-plot) between the static Young's modulus (tangent from axial stress-strain curve in the range of 40-60% of differential peak stress) derived from the triaxial loading part of the test and 1/B derived from the fitting of the subsequent creep data.These results suggest that: (i) weaker samples (higher weak fraction, clay-and organic-rich-sub-group 1) tend to display higher values of B and n than carbonate-dominated samples (sub-group 2); and (ii) E and 1/B are nearly equal, with a high confidence coefficient of R 2 = 0.93.It is noticeable that vertical samples deformed orthogonal to bedding tend to display higher values of B and n compared to horizontal samples.Also, based on additional stress relaxation tests on a few samples, the viscoelastic parameters B and n derived independently from creep and relaxation data vary within 2%, which is consistent with prior results published on US gas shale formations by Sone (2012), and supports the hypothesis that the Goldwyer shale samples can be treated as viscoelastic under the in situ testing conditions replicated in the laboratory.These samples also exhibit anisotropic nature in their viscoelastic responses (Th5 sample in Fig. 8a; Table 4 and Fig. 9a vertical samples have relatively higher values of B and compared to horizontal samples,) equivalent to their observed elastic anisotropy (Mandal et al. 2021a, b).

Viscoelastic Numerical Modelling and Predictions
Viscoelastic modelling performed over an engineering timescale (Details in Appendix 2) shows that (i) viscoelasticity (power-law model) is applicable to the Goldwyer shale formation; (ii) the viscoelastic response of this shale is governed in part by its mineralogy and the fraction of weak phase (clays, organic matter, and porosity); (iii) the ratio of viscous-to-elastic strain/stress is significant.This implies that lithological layers with different compositions behave differently at the length and times scales of the hydraulic stimulation project, even if they have similar elastic (seismic) properties.Therefore, large errors can be introduced in the prediction of the stress state and deformation at depth if the time-dependent response of the shale is ignored.

Viscoelastic Stress Relaxation over Depositional Timescales
Existing methods for estimating the horizontal stress components S hmin and S Hmax at depth often refer to Eaton's approach (Eaton 1969;Thiercelin and Plumb 1994).This method predicts the horizontal stress components within a subsurface porous formation, a 3D volume with vertical boundaries subjected to a horizontal tectonic loading and vertical gravity loading.It accounts for the possible vertical transverse isotropy (VTI) of the formation, and poroelastic effects, i.e., where E v (E 33 ) and E h (= E 11 ) are the vertical and horizontal static young's modulus, respectively; ν v (= ν 31 = ν 32 ) and ν h (6)  Here the vertical direction is defined as X3, while the plane of symmetry (bedding) is horizontal (X1, X2).S v is the gravity-driven vertical stress; α is considered as isotropic Biot's consolidation coefficient for anisotropic material; P p is the formation pore pressure; ε h and ε H are the tectonic strains in the direction of S hmin and S Hmax , respectively.In the right-hand side of each equation (Eq.6), the first term describes the part of the horizontal stress generated by gravitational loading in the anisotropic formation assuming vertical transverse isotropy, as well as accounting for poroelastic effects.The second term describes the part of the horizontal stress generated by tectonic loading at the vertical and horizontal boundaries of the domain.
The viscoelastic parameters of the Goldwyer shale derived from laboratory triaxial stress-strain data are used to estimate present-day stress state at depth, accounting for (i) the short-term elastic response; (ii) the primary timedependent stress and strain response over approximately 6 h; (iii) available geological information about the basin, i.e., interpreted strain rates and duration of tectonic loading).However, we do not explicitly account for complex depositional processes, such as burial, diagenesis, and organic matter maturation.
Considering a constant-rate tectonic strain loading εo at the vertical boundaries of the subsurface domain of interest, the evolution of the stress with time is given by Eq. 19 (Appendix 1), which, in this specific scenario becomes Given the viscoelastic creep constitutive parameters B and n obtained from laboratory triaxial data, at a prescribed loading rate, the accumulated differential horizontal stress can be calculated over a given time duration t [s], within the described model assumptions, i.e., neglecting complex processes such as uplift and erosion.
(  For the Goldwyer formation, we use as input a tectonic loading duration of approximately 300 million years (Ma), which is two-third of the total geological age of the formation (440 Ma), during which maturation and diagenesis occurred, i.e., prior to any uplift or erosion event (Ghori and Haines 2006).It is indeed reasonable to assume that the processes governing the present-day mechanical properties of the shale result primarily from the maturation and diagenesis stages, rather than from the uplift and/or erosions episodes.Different gas shale formations may have different geological ages.However, Sone and Zoback (2014a) and Yang and Zoback (2016) have shown that the duration of the tectonic loading in these circumstances is irrelevant unless the geological timeframe of tectonic loading varies by an order of magnitude between reservoirs.
We focus here on gas shale prospects located within an intraplate region (Delle Piane et al. 2015;Ghori and Haines 2006).Relying on the rigidity of the plates over geological time scales provides an approximate upper bound of the intraplate tectonic strain rate of the order of 10 -18 s −1 (Haines 2004;Zoback and Townend 2001).The tectonic strain rate in the Goldwyer formation is derived from the computed 15 MPa horizontal stress difference (S Hmax -S hmin ) value of the vertical Theia-1 well in the Goldwyer gas shale (Mandal et al. 2020a).The value of S Hmax -S hmin at a specific depth point of G-III unit is derived from the field measured Leak-off test (LOT) and modelled S Hmax magnitude via the stress polygon method (Mandal et al. 2020a;Zoback 2010).Laboratory triaxial deformation experiments yield a range of horizontal static Young's modulus of Goldwyer gas shale between 20 and 35 GPa.The process required an average elastic strain rate of 4-7.5 × 10 -4 s −1 over the inferred geological time of 300 Ma to accommodate the evaluated horizontal stress difference if only elastic deformation is accounted for.Following this computation, we come up with a lower bound of strain rate which varies between 1 and 4 × 10 -19 s −1 .Given the bounds of intraplate strain rate, we used a strain rate of 2 x 10 -19 s −1 in the viscoelastic rheology-driven differential stress estimation, which is similar to outcomes from previous studies conducted on Barnett shale by Sone and Zoback (2014b) and Yang et al. (2015).
The horizontal differential stress (S Hmax -S hmin ) at depth accumulated over 300 Ma at an average intraplate tectonic strain rate of 2 x 10 -19 s −1 is shown in Fig. 10 as iso-stress contours in the (B, n) parametric space.The viscoelastic parameters B and n derived from the triaxial creep experiments on the vertical and horizontal Goldwyer shale samples are also displayed.
The viscoelastic rheology model predicts horizontal differential stress that reaches up to 15 MPa, which is in a close match with the observed value in the Goldwyer gas shale (Mandal et al. 2020a).The coloured contour lines represent the stress relaxation predictions after the application of a strain step perturbation.
The values shown on the top of the contours of Fig. 10 refer to pure elastic response of the gas shale when intersected horizontal axes at n = 0.It is resulting from the strain step perturbation.A higher magnitude elastic stress is expected for shale with lower values of elastic compliance B. The role of viscous stress relaxation comes into effect with the increment of n values along y-axis i.e., stress relaxation dominance of shale rock relative to the instantaneous elastic response.
It can be noticed that contours are closer to horizontal than the engineering timescale predicted values shown in Appendix 2. This emphasises the influence of creep constitutive component n on calculating current differential in situ stress magnitudes.In the case of n > 0.04, the amount of differential stress retained by rock is between 2 and 5 MPa irrespective of their elastic response.
Depending upon the rock's power-law exponent n value, it would be possible to guess the in situ differential stress magnitudes without knowing its elastic behaviour.This is because viscoelastic properties in gas shale reservoirs are mostly controlled by their composition and rock fabric.Therefore, these viscoelastic properties are indirectly related to lithology.

Well-Logs and Continuous Profiles of Rock Properties
Since obtaining a continuous profile of mechanical and elastic properties of rocks beyond reservoir intervals largely depend on sonic logs (compressional wave velocity V p and shear wave velocity V s ), quality control of these logs has utmost importance in mechanical characterization of subsurface rocks.These logs are required to extrapolate stress magnitudes, mechanical and elastic properties in carbonate-dominated Nita, Goldwyer-II, and Willara formations to observe their variations across lithologies.However, the experimental viscoelastic fitting parameters are mostly confined within G-I and G-III units.To overcome this limitation, the composition of each lithology obtained from Fourier Transform Infrared Spectroscopy (FTIR) technique (Mandal et al. 2020b) was plotted in a ternary diagram (Fig. 11a) along with investigated Goldwyer shale sample's mineral composition and organic quantity from pyrolysis and XRD analysis.
Lithologies of G-I and G-III units vary across wide mineral compositions spanning quartz, calcite, and clay with more than 50 vol% of clay-rich minerals.Moreover, the composition of the laboratory samples does not effectively cover the low-carbonate to carbonate-rich shale intervals of the Goldwyer formation.To fill that data gap, compositional representative US gas shale dataset is included for which creep constitutive parameters are acquired at similar experimental conditions and their mineralogy dataset are readily available in the literature (Sone 2012;Sone and Zoback 2014b).This additional data aided in obtaining a continuous profile of creep exponent parameter n where continuous depth profile of creep elastic compliance constant B is built from the inverse of advanced sonic log derived static Young's modulus.
Compressional wave and shear wave velocities were plotted from each lithological unit and laboratory-measured ultrasonic data at in situ condition of Goldwyer shale samples (Fig. 11b).The acoustic velocity ranges of the three limestone formations are within the velocity range of the Goldwyer shale formation.The laboratory ultrasonic velocity covers the whole spectrum of sonic logs.As a result, the mechanical and elastic properties of those limestone formations are within the studied samples from Goldwyer shale, which provides confidence for extrapolation beyond the reservoir intervals.Most of the velocity data Fig. 11 Mineral composition and acoustic log data from each stratigraphic unit of vertical Theia-1 well and comparison with laboratory measurement.Studied cores were recovered from Theia-1 well.a Ternary diagram of mineral compositions determined from FTIR technique and laboratory composition from XRD and pyrolysis.b Crossplot of acoustic V p versus V s and ultrasonic velocity of studied core plugs at in situ stress conditions.Solid and dashed curves refer to the empirical trend of brine-saturated shale, limestone, and dolomite, respectively (Castagna and Backus 1993).The pink line represents the dry organic-rich shale empirical trend (Omovie and Castagna 2019).Laboratory ultrasonic datapoint of clay and organic-rich shales follow the empirical dry organic pattern clusters are lies around and within brine saturated carbonate and shale empirical trend as reported by (Castagna and Backus 1993;Omovie and Castagna 2019).Datapoints lie in the left section of shale empirical trend line are relevant to the observed organic-rich gas shale trend.No obvious outlier was identified from the sonic log cross-plot, which gives further confidence in its usability for deriving continuous profiles of elastic and mechanical properties beyond the target zone.

Stress Profile Based on Viscous Stress Relaxation
The study conducted by Mandal et al. (2020b) reported that the stress state at the Goldwyer shale formation is hybrid (S v ≈ S Hmax > S hmin ) with a combination of normal faulting (NF)-strike-slip faulting (SSF) from conventional stress estimation method, while Bailey et al. ( 2021) proposed a depth-dependent stress regime in the Canning Basin.The differential horizontal stress accumulation can be predicted from Eq. 7 using a constant tectonic strain loading rate over geological timeframe along different lithological intervals of Goldwyer shale formations, although a constant strain rate may be questionable because deformation mostly happened at a non-constant strain rate over geologic time.However, Sone and Zoback (2014a) and Yang et al. (2015) have concluded that total strain (ε 0 ) is more important compared to loading history when evaluating current differential stress magnitudes.
To proceed, viscous stress relaxation and dynamic elastic properties (Young's modulus and Poisson's ratio) were evaluated from the quality-controlled cross-dipole sonic logs (Fig. 13b), followed by static-dynamic conversion available for Goldwyer shale (Mandal et al. 2020b) from the chosen vertical wells PE-1 and TH-1.The continuous profile of B was determined from the inverse of horizontal static Young's modulus (1/E h ), considering the requirement of differential horizontal stress computation which arises from tectonic horizontal strain difference ε diff = ε H -ε h .We have shown in Fig. 9b that static young's modulus is inversely proportional to elastic compliance constant B of laboratory creep data.Building a continuous profile of n along the depth, extrapolating beyond clay-rich intervals, necessitates sufficient laboratory creep measurement from carbonate-rich intervals.To mitigate data gaps and broaden the statistical significance of creep database and their usability, US gas shales (Haynesville and Eagle Ford) with similar experimental conditions was also included (Fig. 11a).From the empirical relationship (Fig. 12a) between creep constitutive parameters B and n as described by Eq. 8, profile of n was obtained as: (8) n = 0.0934 x B 0.344 , with R 2 = 0.74, for global gas shales Thereafter, the differential horizontal stress (S Hmax -S hmin ) accumulated due to tectonic loading was computed from Eq. 7. Finally, the least principal stress magnitude along depth was defined based upon constraining relative variation of in situ stress magnitudes ϕ with depth, where ϕ is defined using Eq. 9 provided by (Angelier 1979)  where S max , S int , and S min are maximum, intermediate, and minimum principal stress magnitude, respectively, and ϕ is a constant that varies between 0 and 1 when intermediate principal stress magnitude S int moves from S int = S min to S int = S max , respectively.In a physical sense, ϕ refers to the kinematics of slip along faults and can be used to define faulting categories (normal, strike-slip, combination of normal and strike-slip, and others).With the above assumption of a constant value of ϕ in a specific stress regime at depth, (9) Φ = S int − S min S max − S min , S hmin and S Hmax can be calculated from Eq. 10 with a viscous stress relaxation approach: Contours of constant ϕ are presented in Fig. 12b through a stress polygon on the defined stress regime based upon the presence of drilling-induced tensile fracture (DITF), rock strength properties, drilling data, and other components of principal stress magnitudes.It would be fair to consider that kinematics of the sedimentary formation remained uniform since different lithological units are part of it.Uniformity of ϕ with depth is in view of confirming regional crustal strain constraints (Sone and Zoback 2014b;Zoback 2010).The vertical stress S v is defined with a stress gradient of 23.1 MPa/km (Mandal et al. 2020b).
Despite our focus on viscous stress relaxation on two horizontal principal stress components, experimental creep data on vertical samples showed stress relaxation is also occurring under a vertical and a horizontal stress.The combined outcome of simultaneous differential stress relaxation among all principal stress components changes the stress state towards an isotropic stress condition (S v = S Hmax = S hmin ), which accords with a constant stress path shown in Fig. 12b.Since stress relaxation occurs along a constant ϕ, the effects of vertical stress relaxation do not impact spatial constraint on ϕ.From the constrained stress magnitudes at 1500 m depth through the stress polygon, we confine ϕ between 0.5 and 1 in a NF/SSF regime.Given the limitation of rock frictional strength, ϕ = 0.6 honours the upper limit of effective horizontal stress ratio within 3.1-4.3for frictional coefficients ranging between 0.6 and 0.8, assuming a strike-slip faulting scenario.The viscous stress relaxation model predicted stress profiles Fig. 14 Stress profile of TH-1 well estimated through viscous stress relaxation approach from well logs, laboratory data and relative in situ stress variation constraint.Tracks 1-4 display gamma ray (GR), creep power-law coefficients (B, n), differential horizontal stress for ϕ = 0.6 and 1-D in situ stress profile, which includes three principal stress components, directly measured LOT, regional S hmin , S Hmax bounds, and pore pressure, respectively.Black dashed lines represent the lithological boundaries interpreted from stratigraphy of Goldwyer gas shale formation are displayed in Figs. 13 and 14 where S hmin profile is calibrated with regional S hmin gradient of 18.5 MPa/km in onshore Canning basin and a nearby well leak-off test (LOT) data.A pretty good match is observed at depth point and stable regional gradient trend within clay-rich interval.The continuous profile of B and n honour the stratigraphic layering (track 2 in Figs. 13 and 14 for PE-1 and TH-1 well) and provides more confidence to the predicted model in shale intervals despite limited field measurements.The S Hmax profile is also confined within the pre-defined range of S Hmax magnitude obtained from the stress polygon.

Discussion
Compilation of different physical and rock mechanical deformation measurements such as creep, pore size distribution, specific surface area, mineral compositions, and frictional properties are required to interlink when analysing long-term effects of geomechanical impact and in situ stress of unconventional gas shale.Although viscoelastic modelling explains the consideration of time-dependent behaviour of rocks when dealing with gas shale production related activity from days to years, but extrapolating to depositional time scales requires better understanding of physical principle contributing to creep deformation.Rock samples of the studied gas shale formation have grain size in the nanometre range (Iqbal et al. 2021;Kohli and Zoback 2013;Passey et al. 2010;Rybacki et al. 2015) and have difficulty in describing deformation mechanisms only through qualitative microstructure analysis techniques, such as SEM, thin section, and X-ray CT scanning, from pre-and postdeformation tests.Moreover, an approximate order of axial creep strain magnitude of 1 × 10 -5 is reported in our triaxial tests, which makes it harder to identify direct microstructural changes owing to creep.Thus, we are going to discuss some obvious correlations that may provide insight for clarity on the physical mechanisms contributing to time-dependent deformation of gas shales.Sone and Zoback (2014a) demonstrated that the poroelastic effects had minimal contribution to the creep behaviour of studied gas shales when comparing time-dependent deformation of room-dry and oven-dried Haynesville samples under constant axial stress.Li and Ghassemi (2012), Sone andZoback (2014a), andRassouli andZoback (2018) did not find any dependence of creep on in situ confining pressure ranges.Rybacki et al. (2017) computed primary creep strain response of US shale under 30 MPa differential stress after 3 years using 1-D power-law equation and noticed that the compositional influence on primary creep is of similar order of magnitude as from the temperature and confining pressure influence on Posidonia shale.It is noteworthy that temperature and confining pressure conditions were tested beyond the Posidonia sample's original in situ stress conditions.The presence of water possibly enhances the amount of creep in the gas shale referred from the studies conducted by Sone and Zoback (2014a) and Almasoodi et al. (2014).The water effects could be effective through induced poroelasticity from pore pressure perturbation, capillary suction, swelling, and microcrack enhancement in highly stressed clay-rich rock with dominant swelling clay materials, such as Smectite and Illite-Smectite (I-S) mixture.In fine-grained mudrocks, the grain size is an order of magnitude smaller and the surface area is an order of magnitude larger compared to siliciclastic rocks (Passey et al. 2010).In the studied clay-rich rocks, which contains mostly Illite, the water Fig. 15 a Axial strain response of Savonnerie sample for different confining pressures where the axial strain was normalized by applied differential stress.b Regression analysis in log 10 J-log 10 t space to derive creep constitutive parameters.The initial instantaneous stress application grey colour coded section was discarded from the analysis retention capacity of clay surfaces played an important role in creep deformation.

Physical Significance of Creep
To understand the impact of grain rearrangement through pore collapse and the relative motion of grains, additional creep experiments were conducted using Berea sandstone and Savonnerie limestone as a proxy for strong and intermediately strong phase components.The experiment was conducted at room temperature and dry conditions at varying confining pressure conditions.The creep responses are presented in Fig. 15.Previously Sone and Zoback (2014a) and Rassouli and Zoback (2018) showed first-order influence of differential stress on axial creep strain.Our observations are no different from those of additional sample creep responses (Fig. 15a), as well as that from a clay-rich shale sample (Th10).In Fig. 15b, the creep compliance function of Savonnerie limestone can be described by 1-D power-law formulation with relatively small slope in log 10 J(t)-log 10 t space, i.e., lower value of power-law exponent n.Before the creep test and post creep deformation, porosity and permeability were also measured to observe the change in grain readjustment due to compaction of pore network of the sample from creep.Presence of large pore network in Savonnerie limestone allowed compaction during creep through reduction of pore volume by ~ 9%.Since the time-dependent behaviour of our studied samples mostly occurred in the axial direction followed by slight compaction, we expect the amount of pore volume (Micro to Meso pores) within and in between clay minerals to directly influence creep response through grain alternation and reduction of pore space contained inside weak phase component of gas shale.As reported in Table 3, most of the mesopore volume is confined within clay-rich samples and have contributed 30-70% of total pore volume.Further, we investigated SSA value S N2 and its correlation with weak mineral phase and mesopore volume.Most clay-rich rocks have higher SSA which is confirmed through its strong linear relationship (Eq.11 and Fig. 16a) with ClayTocPHI and therefore indirectly contributing to the interpreted mesopore volume (Fig. 16b).The correlation between S N2 and weak phase volume ClayTocPHI reads: Mesopore volume correlates with S N2 reflecting a higher correlation coefficient follows Sone and Zoback (2013b) and Kohli and Zoback (2013) emphasized load-bearing framework moved from clastic and carbonate supported to clay and organic matter supported when the amount of clay and TOC content reaches above ~ 30-40 vol%.Relying on that investigation, loadbearing capacity of most of the current investigated shale samples is expected to be controlled by weak phase constituents.Hence, the amount of weak phase becomes a proxy for deformation, especially the time-dependent component.

Impact of Specific Surface Area and Pore-Size Distribution
The surface of the grain is where most chemical reactions and rock mechanical deformations happen over depositional history.Hussaini and Dvorkin (2021) found a complex relationship between specific surface area and porosity in low to medium porosity sandstone and carbonate rocks.Also, Mavko et al. (2020) found the influence of SSA and porosity on the geometric effects of permeability as illustrated by the Kozeny-Carman equation.Although we are investigating the creep and failure frictional behaviour of ultra-low permeable gas shale, we believe SSA may act as a proxy to find any possible interrelationship with mechanical deformation.The gas shales reservoir rocks are predominantly clay rich and have a large specific surface area (Passey et al. 2010) depending upon clay type (Smectite has the largest total surface area of ~ 800 m 2 /g while Illite and Kaolinite have surface areas of the order of ~ 15-30 m 2 /g).We have tried to understand if indirectly derived S N2 can provide any empirical correlation (Eqs.13-14) with creep constitutive parameters as follows (see Fig. 16c, d).
It is clear from Fig. 16c, d and Eq.13-14 that the powerlaw exponent n and elastic compliance constant B are well correlated with S N2 .
Given the established relationship between S N2 and weak phase volume ClayTocPHI, followed by a strong correlation between creep parameters B, n with S N2 , it is evident to recognize how mineral composition and microstructure (here SSA describe the semi-quantitative measure of microstructure) possibly influence the creep deformation of clay and organic-rich shales.These interrelationships testify to the significance of the specific surface area of the studied gas shale reservoir rocks.Therefore, creep parameters n and B can be alternatively estimated from the studied rock formation's SSA without going through expensive experimental laboratory triaxial tests.

Limitation of Creep Measurement and Empirical Correlations
Creep deformation reported in this study both for bedding parallel and perpendicular shale samples of Goldwyer formation was acquired under 1-D consideration, neither reflecting actual far-field boundaries nor the formation pore pressure conditions.It is an indicator of the assumption of a drained creep measurement over a geological time frame owing to small tectonic loading rates.Gas shale reservoir rocks are mostly located under higher temperature and pressure conditions.Samples were received from the core library and then tested in the laboratory.Therefore, exact field conditions cannot be replicated.For example, (i) triaxial tests were conducted at room temperature, whereas reservoir temperature can be above ~ 100 °C; (ii) stress field at depth is anisotropic, whereas tested samples are exposed to axisymmetry conditions in which two horizontal stress components are equal (S Hmax = S hmin = S h ) and vertical stress is always greater than or equal to the horizontal stress (S v ≥ S h ); (iii) samples were tested in as-received condition at ambient pressure, while gas shales are partially saturated with water.Temperature largely influenced the creep behaviour of gas shales in a similar order as mineral compositions of gas shales reported in the literature (Herrmann et al. 2020;Rybacki et al. 2017).Further, creep constitutive parameters were empirically correlated with quantitatively measured petrophysical and elastic properties.In wellbore stability problems, subsurface rock strength properties, UCS, and friction coefficient (µ) are correlated with V p and porosity (Chang et al. 2006;Dewhurst et al. 2015).Although there is no physical resemblance between mechanical strength and compressional wave velocity, a strong correlation was reported.Here, we also refer to that strong correlation to help field engineers build subsurface stress profiles from data-driven correlation function as a proxy approach.For example, we found a relationship between creep constitutive parameters n and B with the semi-quantitative assessment of SSA and weak mineral phase fraction ClayTocPHI.To provide statistical significance of the proposed empiricalbased equations, we have tested the model's applicability by analysing r squared (R 2 ) value together with another statistical variable such as rank correlation coefficient and p-value.Spearman's rank correlation was reported along with their R 2 value and correlation category in Table 5.

Comparison of Stress Prediction Models and Hydraulic Fracturing Simulation
To make a comparison of the viscous stress relaxation model predicted horizontal stress profiles (S hmin , S Hmax ), two horizontal stress profiles (Figs. 18,19) of PE-1 and TH-1 wells were computed from Eaton's modified model (Eq.6).Anisotropic elastic properties (a vertical and horizontal component of Young's modulus and Poisson's ratio) were obtained from cross-dipole sonic logs (Fig. 17) with standard methods (covered in great length by Mandal et al. (2020b) in Appendix 1).We assume an isotropic value of Biot's coefficient as 1 for simplicity, which may not be an actual representation of the studied gas shale sample.Pore pressure was defined from the organic matter corrected total porosity log (Mandal et al. 2020b) with identified overpressure zones at G-I and G-III units.Horizontal tectonic strain components were kept constant with a magnitude of ε h = 3E−04 and ε H = 9E−04 adhering to the regional Canning basin stress study by Bailey et al. (2021).Two stress profiles derived from viscous stress relaxation and Eaton's modified model, respectively, are compared in Figs.18 and 19.Comparison of S hmin profiles from these two methods produced completely different stress profiles along depth and thus expected a variation of fracture gradient along lithological units as well.Eaton's modified model S hmin is mostly uniform while the other model showed a variable S hmin profile with depth from clay-rich to carbonate-rich formations.S hmin trends from both models show a relatively small difference in fracture gradients at the boundary between the Nita formation and G-I unit and a significant contrast in fracture gradients when the wells enter Willara limestone from the organicrich G-III unit.The viscous stress relaxation model predicts a decrease in S hmin magnitude.This observation indicates that Willara limestone may allow downward fracture propagation, even though this would be less likely using S hmin calculated from Eaton's modified model.To demonstrate the influence of stress layering on hydraulic fracture growth, we simulated the hydraulic fracturing process using a commercial software package (MFrac Suite, Baker Hughes).
A 3D planar hydraulic fracture model was used in the prospective G-III unit to quantify the impact of stress layering on its spatial propagation.Mandal et al (2021a, b) demonstrated the application of the viscous stress relaxation model in Perth Basin where direct field data followed lithology controlled S hmin profile.From the hydraulic fracturing simulation, the authors further showed that the downward propagation of simulated fracture from the perforated shale interval into the shaly-sandstone formation.
The key geomechanical parameters were derived from the available wireline logs, field data, and viscoelastic stress relaxation modelling reported earlier.A constant stress gradient was assumed in each layer of a four-layer geomechanical model (Fig. 20a).Average elastic properties such as Young's modulus and Poisson's ratio were computed from sonic measurement followed by dynamic to static conversion through empirical equations.Assuming a homogeneous rock property in each layer (Table 6), we simulated a vertical well configuration at the centre of the lower G-III unit with the following configurations: (i) slick water injection at a rate of 50 bbl/min for 2 h (ii) single perforation cluster with perforation diameter of 0.3 inch and 12 perforations per cluster (iii) proppant concentration of 1 ppg with 40/70 mesh size.Figure 20b displays the fracture geometry and proppant distribution 1 day after the shut-in of the hydraulic fracturing operation.
There are no continuous direct measurements available to validate the stress models expect a single depth LOT data in the Goldwyer shale formation.For US Barnett shale, it has been reported that the hydraulic fracture in the main reservoir interval propagates downwards into a limestone formation (Sone and Zoback 2014b;Yang and Zoback 2016).Propagation of a single fracture in an isotropic stress field shows uniform elliptical growth (Salimzadeh et al. 2019), while the application of stress layering controls vertical and lateral propagation of the fracture depending upon S hmin contrast at layer boundaries of the perforated interval (Singh et al. 2019).Our previous study in gas shale formations of Perth Basin (Mandal et al. 2021b) reiterated the fact that the fracture propagation and proppant distribution at subsurface intervals are driven by variation of least principal stress at depth.Our simulation shows that lateral hydraulic fracture growth is several times that of its vertical growth (Fig. 20).The fracture propagates in the upward direction where the least principal stress magnitude is lower.The lower layer Fig. 17 Stress profile estimated through Eaton's modified model from well logs and regional tectonic strain constraint.Tracks 1-4 display gamma ray (GR), static Young's modulus, Poisson's ratio, and 1-D in situ stress profile that includes three principal stress components, directly measured LOT, regional S hmin , S Hmax bounds, and pore pressure, respectively.Missing data points at the G-I unit of track 2 are due to poor data quality (Colour figure online) does not act as a full fracture barrier; it is rather a soft barrier, i.e., allows slight downward growth.Yang et al. (2015) showed that the viscous stress relaxation model derived S hmin profile can explain out-of-zone microseismicity, while Eaton's modified model failed to characterise the observed microseismicity beyond the perforated geologic formation.The uncertainty associated with Eaton's modified model is substantial, considering the contribution of elastic anisotropy parameters, Biot's coefficient, pore pressure, and arbitrarily defined horizontal tectonic strains.On the other hand, the advantages of the viscous stress relaxation model are several such as (i) direct utilization of laboratory creep deformation data, (ii) usage of physical rock properties from wireline logs, and (iii) honouring the kinematics condition of stress state within a lithologic interval.Further, Eaton's modified model suggests a hybrid-faulting regime, while the viscoelastic model confirms a consistent strike-slip faulting regime within the Goldwyer shale formation, similar to the present-day in situ stress variation study conducted in the onshore Canning Basin by Bailey et al. (2021).

Failure Frictional Properties, Specific Surface Area, and Weak Phase
We already covered how composition and microstructure can control overall mechanical deformation including creep.Further, we established a reliable correlation between S N2 and creep constitutive parameters.It was demonstrated that production from ultra-low permeable gas shale is feasible through inducing slip along pre-existing fracture networks by slick water injection.Favourably orientated natural fractures and faults experience instantaneous slip during hydraulic fracturing, but slow slip emanates from long-duration and long-period slip events from regional fault networks (Das and Zoback 2011;Kohli and Zoback 2013).Therefore, it is important to dig further to understand the influence of mineral composition and specific surface area on fault stability and frictional strength.Kohli and Zoback (2013) reported a strong influence of clay framework on the slip stability of several gas shale reservoirs (e.g., Barnett, Haynesville, and Eagle ford).Figure 21a shows a moderate negative correlation between failure frictional coefficient µ s and weak phase ClayTocPHI as shown in Eq. 15: (15) s = − 0.003xClayTocPHI + 0.86, with R 2 = 0.40, for Goldwyer gas shale.However, a rapid decrease in µ s (from an average of 0.8 to 0.6) value was observed when the weak phase fraction reached above 40 vol%.This behaviour indicates a possible shift from carbonate to a clay-rich load-bearing framework.The correlation between µ s and S N2 in Fig. 21b reads a negative trend with a correlation coefficient of 0.66.Equation 16presents the correlation as Again, this relationship points to the importance of the specific surface area in analysing intact mechanical properties as well as failure frictional properties of gas shale.A strong positive relationship between µ s and static Young's modulus (E) is reported separately for Goldwyer and US gas shales (Eq.17), equivalent to the observed correlation between mechanical compressive strength and Young's modulus (E) of intact global gas shales (Mandal et al. 2022a).( 16) s = −0.015xSN2 + 0.83, with R 2 = 0.66, for Goldwyer gas shale.
These empirically derived equations for Goldwyer and US gas shales allow direct prediction of continuous sliding friction coefficient profiles from cross-dipole sonic logs.Therefore, µ s profile indirectly helps in identifying possible intervals that are prone to shear slip following hydraulic fracturing operation.

Conclusions
Following acquisition and interpretation of short-duration primary creep, viscous stress relaxation modelling, analysis of failure frictional properties of Goldwyer gas shale, and interlinking the aforementioned weak mineral phase fraction, pore size distribution, specific surface area, and elastic properties, the following can be concluded: (i) Clay and organic-rich G-III gas shale reservoir unit accommodates more primary axial creep under constant differential stress as established from axial creep versus time plot of the studied samples and the increasing magnitude of creep power-law exponent n with higher clay content.(ii) Bedding perpendicular creep deformation magnitude and creep rate are higher than in those of bedding parallel direction, contributing to the anisotropic nature of creep.This observation is in line with the elastic anisotropy of gas shales.(iii) 1-D power-law creep constitutive parameters describing creep deformation were used successfully through the viscoelastic stress relaxation approach to build a continuous S hmin profile.The predicted stress profile is tied to the altering lithology of the studied formations.(iv) Clay and organic-rich shale samples have higher specific surface areas (S N2 ).From the reported relationship between ClayTocPHI and S N2 , it is imperative to consider that indirect assessment of gas shale's specific surface area may shed some light on long-term creep behaviour.(v) 1-D power-law creep parameters (n and B) can be estimated with some uncertainty from a semiquantitative assessment of gas shale's SSA value inverted from low-pressure nitrogen gas adsorption isotherms, despite the difficulty in providing a direct physical explanation.Alternatively, elastic compliance constant B can be calculated from the inverse of static Young's modulus built from recorded sonic logs.(vi) Lower stresses at the overlaying and underlying layer of the gas shale reservoir resulting from viscous stress relaxation driven S hmin profile allows the modelled artificial hydraulic fracture to propagate vertically into the carbonate-dominated layer.(vii) An empirical correlation between sliding friction coefficient µ s and specific surface area S N2 was observed, which could be used by a field engineer to create continuous frictional failure profiles at depth in the absence of experimental deformation data.

Appendix 1: Viscoelastic Theory
The theory of viscoelasticity implies that the total stress Σσ at any time t after application of a unit step strain load at t = 0 (Heaviside function) is given by the integral of the corresponding time-dependent strain ( ) over all previous times τ, i.e.
where E(t) is the stress relaxation modulus/function (See   1 3 during stress relaxation at constant strain (zero displacements) after application of the unit strain step at t = 0 reads It is well established that the reciprocity principle holds for viscoelastic materials, so that the strain changes at constant stress after a unit stress step is applied at t = 0 to a viscoelastic material can be written a priori as: where J(t) is the creep strain compliance function characterising the material (in Pa −1 ), reciprocal of the relaxation modulus E(t).This reciprocity between E(t) and J(t) in the time domain corresponds to an inverse relationship between their respective transform image e(s) and j(s) in the Laplace domain, i.e., where s is the complex Laplace variable.It is therefore necessary and sufficient to determine the functional form of either J(t) or E(t) to compute the other using the Laplace transform.

Appendix 2: Viscoelastic Response over Engineering Timescales
From linear viscoelasticity (covered in the modelling section), we estimate the accumulated creep strain and the residual stress after relaxation expected in the Goldwyer shale at depth, to quantify their impact on reservoir deformation and in situ stress changes at the time scale of the project (21) e(s) ⋅ j(s) = 1 s 2 , life cycle, i.e., responses 1 day and 1 year after application of the perturbation.This time scale is useful for (i) borehole stability analyses and predictions in the near-wellbore region (days) and for (ii) production forecast and reservoir stability analyses and predictions in the far-field region (years).
The viscoelastic analysis and predictions presented in Fig. 23 are based on stress or strain perturbations of 10 MPa or 1.5 × 10 -4 , respectively.These values are taken from the literature and correspond to typical drilling or fluid injection operations in the field (Sone and Zoback 2014a;Zoback and Kohli 2019).Note that this first-order viscoelastic analysis considers a rather simple 1D deformation field and neglects 3D boundary conditions and pore pressure effects.The aim is to assess whether accounting for viscoelasticity (beyond simplistic elasticity) in unconventional reservoirs improves reservoir and borehole predictions, before delving into an exhaustive 3D numerical analysis accounting for specific boundary conditions and actual project complexities, e.g., complex geological settings.Figure 23 shows the predicted accumulation and relaxation of viscoelastic strain and stress in the shale after 1 day (Fig. 23a, b), and after 1 year (Fig. 23c, d) for multiple combinations of B and n values.The amount of viscous strain change relative to the instantaneous elastic response is computed using Eq. 1.The time-dependent viscous stress relaxation relative to the instantaneous elastic response is computed using Eq. 5.The coloured contour lines displayed in Fig. 23a (or Fig. 23c) represent the total creep strain accumulated by the shale 1 day (or 1 year) after the application of a differential stress-step perturbation of 10 MPa.The reciprocal stress relaxation predictions 1 day (or 1 year) after the application of a strain step perturbation of 1.5 × 10 -4 is shown in Fig. 23b (or Fig. 23d).The response of the shale to a stress and strain perturbation after 1 day and after 1 year are qualitatively similar, although more viscous strain and stress develop over the longer time duration.
Note that in these graphs, the intersection of each contour line with the horizontal axis (n = 0) represents the purely elastic response of the shale, i.e., elastic strain resulting from the stress-step perturbation (Fig. 23a, c), and elastic stress resulting from the strain step perturbation (Fig. 23b, d).As expected, a higher elastic strain is achieved in response to a stress-step perturbation for shales with higher values of elastic compliance B. On the other hand, higher elastic stress is achieved in response to a strain step perturbation for shales with lower values of elastic compliance B. Larger values of n along the y-axis represent a larger contribution of viscous strain and stress to the total strain and stress undergone by the shale, i.e., creep strain or stress relaxation dominate the shale response relative to the instantaneous elastic contribution.For instance, when n > 0.02 the residual stress after viscous stress relaxation represents approximately two-thirds of the instantaneous elastic stress change.
The laboratory-derived viscoelastic parameters B and n for the vertical and horizontal Goldwyer shale samples are superimposed on the contour lines of the graphs in Fig. 23.The comparison between the predictions (contour lines) and the laboratory data suggests that the amount of viscous strain and stress in the Goldwyer shale is significant relative to the elastic response, i.e., n is larger than 0.02 for all samples but one (vertical sample).where σ A and p c stand for the differential axial stress and the effective confining pressure, respectively; β is the angle between the axial stress direction and the normal to the fault plane.Note that β is related to the internal friction angle ϴ as: We estimated the angle of the fault plane in each sample from the 3-D X-ray CT images obtained after triaxial tests, i.e., β varies between 50° to 70°, with an average value of 62° ± 4°, regardless of their sample and bedding orientation.The values of the angle β independently derived from the Mohr-Coulomb failure envelope lies within ± 5° of the directly measured fault angle values.Figure 24c illustrates the definition of the fault plane and its spatial orientation with respect to the direction of the differential stress, as well as the normal and shear stresses acting on the fault plane.The post-failure slip achieved along the fault plane by the end of each MST is in the range of 2-5 mm corresponding to an overall average shortening of 4%.
Figure 24d shows a representative 2D slice extracted from the complete 3D X-ray CT images of two contrasting samples of Goldwyer shale.The jagged and rough nature of the fault's surface is visible at this resolution.Besides the primary shear fault, multiple fracture-like features aligned with the bedding plane (planes of weakness) were also identified; they could have been induced by the post-test unloading of the sample at the end of the MST.Table 7 summarizes the frictional failure data for both vertical and horizontal samples.The data does not show an obvious directional dependency of the residual strength or the friction coefficient.The average value of µ s is 0.7 ± 0.1.tion, 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 holder.To view a copy of this licence, visit http:// creat iveco mmons.org/ licen ses/ by/4.0/.
(i) Target gas shale prospects originated from the deposited marine shales of the middle Ordovician Goldwyer formation in the onshore Broome Platform, Canning Basin.(ii) Vertical Theia-1 (TH-1) and Pictor East-1 (PE-1) wells intersected the Goldwyer formation of three stratigraphic subunits Goldwyer-I (G-I), Goldwyer-II (G-II), and Goldwyer-III (G-III) with a total depth of 1645 m and 1706 m, respectively.All necessary wireline logs are available from these two wells for the prediction of continuous least principal stress profiles through viscoelastic modelling.(iii) Clay and organic-rich G-I and G-III units are sandwiched between three carbonate-dominated intervals: Nita formation, G-II unit, and Willarra formation (Fig. 1).

Fig. 1
Fig. 1 Generalized stratigraphy of the onshore Canning Basin.Goldwyer formation is highlighted in red colour within Middle Ordovician [Modified from DMIRS (2019)].The right part shows the stepwise stratigraphy of the Middle Ordovician depositional age (Colour figure online)

Fig. 3 a
Fig. 3 a Micromeritics (R) VacPrep 061; b micromeritics (R) TriStar II PLUS instrument.Overall setup for sample preparation and low-pressure nitrogen gas adsorption

Fig. 4 a
Fig. 4 a Autonomous triaxial cell used for multistage deformation experiment including creep reported in this study.Axial and radial strains are acquired with axial linear differential displacement transducers (LVDTs) and radial strain cantilever as shown in this instru-

Fig. 5 a
Fig. 5 a Multistage triaxial experimental design over the five loading/unloading stages including stress relaxation and creep deformation during fourth stage.At the last stage, the sample is allowed to reach the failure point followed by residual strength recording.Loading-unloading cycles with mild hysteresis resulting from irreversible

Fig. 6
Fig. 6 Nitrogen gas adsorption experiment conducted on Goldwyer gas shale formation to quantitatively analyse pore structure and pore volume.a Carbonate dominated samples, b clay and organic-rich samples.A higher amount of gas adsorption of Th2 sample is an indication of heterogeneity of gas shale formation where the semi-quanti-

Fig. 7 a
Fig.7a Pore size distribution of Goldwyer gas shale formation from low-pressure gas adsorption isotherm curves.Th4, 5, and 8 are dominated by carbonates while Th7, 9, and 10 samples have a higher percentage of clay and organic matter.b Correlation between total porosity and mesopores volume interpreted from LPNA test.Samples are colour coded with their weak phase ClayTocPHI in vol% (Colour figure online)

Fig. 8 a
Fig. 8 a Few representative samples creep response under constant axial differential stress.V and H refer to vertical and horizontal samples.Vertical samples accommodate more time-dependent deformation than horizontal one (Th-5 V > Th-5H).b Linear regression of creep compliance data in log 10 t − log 10 J space for a tested sample to determine creep constitutive parameters.The initial section of the recorded data was discarded considering the pressure ramps up to reach the desired axial differential stress for creep data acquisition

Fig. 9 a
Fig. 9 a Creep constitutive parameters (B, n) obtained from linear regression analysis representing a 1-D power-law model of creep compliance data.All investigated samples from bedding parallel and perpendicular directions are colour coded with their total clay, TOC, and porosity vol%.b Relationship between static young's modulus E interpreted from stress-strain data and stiffness of the rock derived from the inverse of the elastic compliance constant B of the 1-D power-law model.The black line refers to 1 to 1 relationship between E and 1/B (Colour figure online)

Fig. 10
Fig. 10 Differential in situ stress accumulation over 300 Ma geological age at a constant intraplate strain rate of 2 × 10 -19 s −1 derived from 1-D viscoelastic rheology model.Values shown on the top of contour refer to estimated differential stress magnitudes (in MPa unit) are overlain on laboratory viscoelastic properties B and n of vertical and horizontal sample plugs of the Goldwyer shale formation (Colour figure online)

Fig. 12 a
Fig. 12 a Cross-plot of elastic compliance constant B and power-law exponent n.The Black dashed line refers to the established empirical relationship between B and n with the inclusion of Goldwyer and US gas shales.b Stress polygon to define in situ stress regime at 1.5 km depth with frictional strength limit constraint.Existence of DITFs and absence of borehole breakouts BOs are limiting the stress regime to the left side of tensile strength T 0 = 0 and below unconfined compressive strength C 0 = 60 MPa contour lines, respectively.The light blue area defines the allowable stress regime i.e., mostly in strike-slip faulting SSF but in close connection with the boundary of normalfaulting NF regime for the given condition of DITF presence and absence of BO.Constant contours of ϕ values are overlain within the defined stress regime (ϕ value can varies between 0.5 and 1 in the constrained stress regime) Fig. 13Stress profile of PE-1 well estimated through viscous stress relaxation approach from well logs, laboratory data, and relative in situ stress variation constrain.Tracks 1-4 display gamma ray (GR), creep power-law coefficients (B, n), differential horizontal stress for ϕ = 0.6 and 1-D in situ stress profile, which includes three princi-

Fig. 16
Fig. 16 Relationship between weak-phase ClayTocPHI, S N2 , mesopore volume and creep constitutive parameters B, n. a S N2 versus ClayTocPHI, b mesopore volume versus S N2 , c correlation of n versus S N2 , d correlation of elastic compliance B versus S N2 .A sample with

Fig. 18
Fig. 18 Comparison of minimum and maximum horizontal stress S hmin and S Hmax profiles of TH-1 well derived from viscous stress relaxation approach and Eaton's modified model.Tracks 1-3 present

Fig. 19
Fig. 19 Comparison of S hmin and S Hmax profiles of PE-1 well derived from viscous stress relaxation approach and Eaton's modified model.Tracks 1-3 present gamma ray, S hmin, and S Hmax profiles from viscoe-

Fig. 20 a
Fig. 20 a S hmin profile at depth for the four layers of the 1-D model used for hydraulic fracturing simulation.The red circle at the G-III unit pointed to the fracture initiation location.b Fracture geometry generated from hydraulic fracture simulation one day after shut-in.

Fig. 22 )
characterising the material (in Pa), and E o , o , and o are the limiting values of E(t), (t) , and (t) as time t approaches zero, respectively.Hence E o , o , and o correspond to the instantaneous elastic response of the viscoelastic material, i.e., o = E o o .Therefore, the stress changes

Fig. 21
Fig. 21 Cross-plots of sliding friction coefficient versus: a weak phase ClayTocPHI, b specific surface area (S N2 ) from nitrogen gas adsorption technique c static Young's modulus (Colour figure online)

Fig. 22
Fig. 22 Time-dependent deformation of a viscoelastic material is generally expressed as a creep strain or b stress relaxation

Fig. 23
Fig. 23 Creep strain and stress relaxation response over engineering time scale under a given instantaneous stress change of 10 MPa and strain perturbation of 0.0015 (1.5 milli-strain).Contours of stress (in MPa unit) and milli-strain are overlain on interpreted viscoelastic creep constitutive parameters B and n. a Total accumulated strain,

Table 1
Mineral composition of clay material at ten selected samples studied here I-S is an Illite-Smectite mixture.Total clay is calculated as a combination of clay minerals and mica

Table 3
Pore size distribution, specific surface area, and total porosity of ten shale intervals studied here Total porosity was estimated from helium gas pycnometer measurement *For Th1 sample, no data was acquired due to instrument malfunction at the time of experiment.No repetition of that sample was feasible

Table 4
Fitted creep parameters for all the studied Goldwyer shale samples Uncertainty of estimated parameters is derived from the error propagation method arising from measurement error

Table 6
Geomechanical properties used in the four-layer simulation