Deformation Experiments on Bowland and Posidonia Shale—Part II: Creep Behavior at In Situ pc–T Conditions

To unravel their long-term creep properties at simulated reservoir conditions, we conducted constant stress deformation experiments at elevated confining pressures, pc = 50–115 MPa, and temperatures, T = 75–150 °C, on Posidonia (GER) and Bowland (UK) shale, which exhibit varying petrophysical and mechanical properties. Depending on applied pc–T conditions and sample composition, recorded creep curves exhibit either only a primary (decelerating) or additionally a secondary (quasi-steady state) and a tertiary (accelerating) creep phase during deformation. At high temperature and axial differential stress and low confining pressure, creep strain is enhanced and a transition from primary towards secondary and tertiary creep behavior is observable. Creep strain of Posidonia shale, which is rich in weak constituents (clay, mica, and organic content), is enhanced when compared to creep strain recorded during deformation of either carbonate- or quartz-rich Bowland shale. Electron microscopy observations revealed that creep strain is mainly accommodated by the deformation of weak minerals and local pore space reduction. In addition, microcrack growth occurred during secondary creep. An empirical correlation between creep strain and time based on a power law was used to describe the decelerating creep phase, also accounting for the influence of confining pressure, temperature, and axial differential stress. The results suggest that the primary creep strain can be correlated with mechanical properties determined from short-term constant strain rate experiments such as static Young’s modulus and triaxial compressive strength.


Introduction
Currently, energy consumption depends largely on conventional and unconventional resources, e.g., coal, oil, or gas (Kuchler 2017). The use of hydrocarbons from unconventional 1 3 reservoirs such as shale plays may help to bridge the transition from conventional towards renewable energy resources (Hausfather 2015; Zhang et al. 2016;Kuchler 2017). For example, shale gas production has increased tremendously over the past decades (Howarth and Ingraffea 2011;McGlade et al. 2013;Wu et al. 2017). Specifically, in North America, the gas production from shale reservoirs has had a major impact on the energy supply (McGlade et al. 2013;Kuchler 2017). In Europe, unconventional shale reservoirs (e.g., Posidonia, Bowland-Hodder formations) also show potential for economical hydrocarbon production (Andrews 2013;Kuchler 2017;Wiseall et al. 2018).
Since proppant embedment and fracture closure depend on the long-term creep behavior of shales, we conducted a series of deformation experiments at simulated reservoir conditions. In an effort to characterize the time-dependent deformation behavior of shale rocks, we conducted tests at a broad range of conditions, since only limited data on the creep properties of shale rocks exist until now (Sone and Zoback 2010, 2013bLi and Ghassemi 2012;Almasoodi et al. 2014;Rassouli and Zoback 2015;Rybacki et al. 2017). We established empirical relations between the stress-induced creep strain and time, which also account for the influence of confining pressure and temperature on the creep rate.
In addition, we compare our results to published data on an immature Lower Jurassic Posidonia shale. Sample material was recovered from a quarry in SW Germany, close to the village of Dotternhausen (DOT) (Rybacki et al. 2015). Details on the sample material and acronyms are given in Table 1.
Compositions of the specimens were determined using X-ray diffraction analysis (XRD). Samples contain a mixture of quartz (Qtz), clay (Cly), carbonates (Cb), mica (Mca), feldspar (Fsp), pyrite (Py), and organic matter (TOC) ( Table 1). Porosity (incl. micropores) was determined by Helium-pycnometry. For Posidonia and Bowland shale samples, porosity ranges between 1 and 12% ( Table 1). Note that, mineral data in Table 1 are given in vol%, as mechanical behavior of the samples largely depends on volumetric fraction and spatial distribution of the components. Posidonia shale is rich in mechanically weak phases (Cly + TOC + Mca + pores) and Upper Bowland shale samples are either quartz-or carbonate rich (Fig. 1, Table 1). For simplicity, we will use the attribute 'weak' in the following for Posidonia shale with respect to sample composition, where 'weak' represents Cly + TOC + Mca and pores. Outcrop material of Upper Bowland shale contains higher fractions of clay, mica, and pores compared to core-derived samples, likely due to weathering (Table 1). Microstructural observations reveal a very fine-grained matrix with a grain size ≤ 20 µm for Posidonia and Bowland shales (Fig. 5a, b). The bedding orientation is characterized by preferred alignment of organic matter and phyllosilicates (Fig. 5).
For deformation experiments, cylindrical samples of 10 mm diameter and 20 mm length were prepared with their axis-oriented perpendicular to bedding orientation. Although not fully representing in situ conditions, all specimens were dried at 50 °C for at least 48 h before starting deformations experiments. Typically, rocks at depth are at least partly saturated with formation fluids, which impacts their mechanical behavior (e.g., Ibanez and Kronenberg 1993;Rybacki et al. 2015). However, for comparison purposes, we dried all specimens prior to deformation, since we are not aware of the formation fluids present in the different shales and the degree of saturation under in situ conditions. In addition, the low permeability of the investigated shales favors testing under undrained conditions if they are fluidrich, which may not represent natural deformation conditions. Note that after drying at 50 °C, Posidonia and Bowland samples still contained a small water fraction, H 2 O , of H 2 O , (POS_HAR) = 0.6 vol% and H 2 O , (BOS_OC) = 1.6 vol%, respectively, which was determined by measuring the weight loss of some samples after drying them at a temperature of 110 °C.

Experimental Methods
We used a Paterson-type deformation apparatus (Paterson 1970) to perform triaxial creep experiments at elevated confining pressures (p c = 50-115 MPa), temperatures (T = 75-150 °C), and axial differential stresses (σ = 9-453 MPa) ( Table 2). For simplicity, we use the term 'stress' instead of 'axial differential stress'. Argon gas was used as the confinement medium due to its inert properties. To prevent intrusion of the Argon gas into the samples, they were jacketed by thin (wall thickness ≈ 0.35 mm) copper   sleeves. Axial load was measured using an internal load cell, installed within the pressure vessel. Calibration runs were performed on pure copper samples to correct for copper jacket strength. Measured axial load was then converted to stress assuming constant volume deformation. Uncertainties of resulting peak stresses are < 4%. Recorded axial displacement, measured by a linear variable differential transformer (LVDT), was corrected for system compliance and converted to total axial strain, ε t . The error of calculated strains is < 6% of the strain value.
All experiments were performed with loading direction normal to bedding orientation. The triaxial compressive strength, σ TCS , and static Young's modulus, E, of investigated samples were determined in previous constant strain rate experiments performed at similar confining pressures and temperatures (Herrmann et al. 2018).
For creep testing, the axial load was increased to the desired level, which typically stabilized within ≈ 20 s. Subsequently, the load was held constant until the end of the experiment. Depending on deformation conditions, experiments were stopped either after sample shear failure or after a fixed time if no failure occurred. The total duration of creep tests was between 22 s and ≈ 15 days with maximum total axial strains of ε t ≈ 0.5-5.5% (Table 2).
Scanning electron microscope observations (SEM) were performed on mechanically polished thin sections using a Zeiss Ultra 55 Plus-microscope. High-resolution analysis of microstructures was performed using transmission electron microscopy (TEM, Fei Tecnai G2 F20 x-Twin) on focused ion beam (FIB, Fei FIB200TEM) prepared foils.

Results
We performed 49 triaxial creep experiments. During 42 experiments, we fixed two of the three parameters p c , T, and σ and changed the remaining third parameter within a given range to investigate the influence of this particular parameter on the resulting creep strain behavior (Table 2). In addition, in one triaxial creep experiment, axial stress was changed stepwise to investigate possible sample strain-hardening effects. The three different creep phases (primary, secondary, tertiary, cf., Fig. 3b) were distinguished in the creep curves (see Sect. 5.1).
After a tertiary creep phase, some samples failed along a single shear fracture inclined at φ ≈ 35 ± 2° with respect to the sample axis. From this, we estimated an apparent coefficient of internal friction of µ i ≈ 0.7 ± 0.05.

Effect of Loading History (Strain Hardening)
To investigate potential effects of loading history on the creep behavior of shale rocks, we performed one experiment at constant confining pressure of 50 MPa, 100 °C temperature and stepwise increased and decreased stresses between σ = 9 and 140 MPa (Fig. 2). The experiment was performed on Posidonia (HAD) shale, with a similar mineral composition as Posidonia (HAR) shale (cf., Fig. 1). During testing, the axial stress was first increased in two steps (σ 1 , σ 2 ), then decreased to σ 3, and subsequently increased up to σ 2 in four steps. Finally, the sample was loaded again by repeating stress steps σ 3 -σ 5 (Fig. 2). For the first two steps, creep strain increases with increasing stress, but for the next 4 steps, the creep strain is almost negligible, even when approaching the same stress, σ 1 , as applied in the first step. This was also found for the last three loading steps. Repeated loading at the highest stress level (σ 2 ) also reveals a considerably lower strain rate for the subsequent loading cycle (see lines in Fig. 2). This decrease in strain rate at similar stress demonstrates strain hardening of the shale, possibly caused by non-reversible pore collapse and closure of preexisting microcracks. Consequently, to exclude any influence of loading history on the creep behavior in the remaining tests, we performed creep experiments using a single loading step at predefined stress levels.

Effect of Composition
Triaxial creep experiments on Bowland (PH1) and Posidonia (HAR) shale were performed at 90 °C temperature and 75 MPa confining pressure, simulating p c -T conditions at about 3 km depth. Creep curves of Bowland shales recovered from borehole PH1 at different depths display varying creep strengths, which change with composition ( Fig. 3a). A clear trend accounting for varying amounts of intermediate strong (Cb) or strong (Qtz) sample constituents is not evident. Quartz-rich samples appear to be stronger than carbonate-rich samples and deform solely by primary creep at the applied stress of 442 ± 11 MPa. The latter exhibit only primary creep (BOS5) or in addition secondary (quasi-steady state) and tertiary (accelerated) creep stages (BOS1) until failure occurred (BOS3). In comparison, the Bowland (OC) shale and Posidonia (HAR) shale showed comparable creep strain values, but at considerably lower stress of 206 ± 4 MPa (Fig. 3b). This relatively low stress level was used, since Posidonia (HAR) would fail immediately when applying higher stresses (Fig. 3a). Under these conditions, quartz-rich Bowland shale collected from the outcrop was distinctively stronger than weak Posidonia (HAR) shale. Moreover, Bowland shale

Effect of Confining Pressure, Temperature, and Stress
Selected creep curves of Posidonia (HAR) and Bowland (OC) shale at various p c , T, and σ conditions are shown in Fig. 4a, c, e and b, d, f, respectively. Creep strains of both shales at T = 90 °C decrease with increasing confining pressure ( Fig. 4a,  At similar stresses and constant confining pressure of 75 MPa, the creep strain of Posidonia (HAR) and Bowland (OC) increases with increasing temperature, but with a higher sensitivity in Posidonia (HAR) shale (Fig. 4c, d). Posidonia (HAR) shale loaded at high temperature deformed also by secondary creep, but Bowland (OC) shale showed only primary creep behavior in the investigated temperature range.
The effect of increasing stress at p c = 75 MPa and T = 90 °C on creep of Posidonia (HAR) and Bowland (OC) shale is shown in Fig. 4e, f. Samples subjected to high stress Fig. 2 Creep curve of Posidonia (HAD) shale recorded at 50 MPa confining pressure, 100 °C temperature and stepwise increased stress, σ i . Initially, creep strain is increasing with increasing stress. Sample strain hardening is evident after repeated loading at similar stress levels by a decrease of associated strain rate (e.g., α for σ 2 )  and Bowland (OC) shale (b, d, f) in relation to a variation of confining pressure (a, b), temperature (c, d) and stress (e, f). Creep curves display mainly primary creep. At high temperatures and stresses and at low confining pressures, the creep curves reveal in addition to primary also secondary and tertiary creep, leading to failure of some samples, as indicated (x). In general, with increasing confining pres-sure, the creep rate of both shales is reduced, whereas it is enhanced with increasing temperature and stress. Outliers of this overall trend are indicated by dashed lines. Weak Posidonia (HAR) shale is more sensitive to changes in confining pressure, temperature and stress than quartz-rich Bowland (OC) shale. For comparison, all curves were cut off at 5000 s, but the duration of most tests was much longer (Table 1). Deformation conditions are indicated display all three creep phases (primary, secondary, and tertiary) in contrast to samples deformed at low stress, which reveal only primary creep behavior. The transition from primary creep to secondary and tertiary creep of both shales appears to occur within a relatively narrow stress range of Δσ ≈ 15 MPa.
In summary, the creep behavior of both shales is sensitive to the applied p c , T, and σ-conditions. Quartz-rich Bowland (OC) shale is less sensitive to variations of p c , T, and σ-conditions than weak Posidonia (HAR) shale. Most creep experiments were finished within 1 day, either by manual termination or by sample failure (Table 1).

Microstructures
The microstructures of deformed samples of Posidonia (HAR) and Bowland (OC) shale were investigated using scanning electron microscopy (SEM, Fig. 5) and transmission electron microscopy (TEM, Fig. 6). The microstructures of samples of both shales that were deformed at low stresses in the primary creep regime (Fig. 5c, d) are hard to distinguish from those of undeformed samples ( Fig. 5a, b). Deformed Posidonia (HAR) shale sometimes reveals slightly deformed pyrite (Fig. 5c, bold black arrow) and intracrystalline fracturing of dolomite (Fig. 5c, open black arrows). Bowland (OC) shale displays minor bending of phyllosilicates (Fig. 5d, bold black arrow).  In a Bowland (OC) shale, sample deformed in tertiary creep, intracrystalline fractures of apatite (open black arrow), and bent phyllosilicates (bold black arrow) are observed (Fig. 5h).
High-resolution TEM micrographs of shales deformed by primary creep are shown in Fig. 6. The overview images of Posidonia (HAR) and Bowland (OC) shale show little evidence of brittle deformation features (Fig. 6a, b, respectively). Instead, the deformation of Posidonia (HAR) shale appears to be accommodated by local pore space reduction between calcite flakes (Fig. 6c, bold black arrow) and bending of weak phyllosilicates around strong quartz grains (Fig. 6c, open black arrow). Primary creep of Bowland (OC) shale is mainly accomplished by bending of weak phyllosilicates around stronger grains, e.g., apatite (Fig. 6d, bold black arrow).
In summary, our microstructure observations on deformed Posidonia and Bowland shale samples indicate that a range of deformation mechanisms operates including plastic dislocation activities and brittle microcracking.

Discussion
The time-dependent deformation (creep) behavior and deformation mechanisms of the investigated Posidonia (HAR) and Bowland (OC) shale rocks strongly depend on sample mineralogy and experimental conditions such as confining pressure, temperature, and stress, as recognized also for other shales Zoback 2013b, 2014;Rybacki et al. 2017). Here, we discuss the influence of these parameters on the creep behavior of Posidonia (HAR) and Bowland (OC) shale and compare the results to literature data that was obtained on other European and North American shale rocks (Sone and Zoback 2014;Rybacki et al. 2017).

Deformation Mechanism
Black shales are usually highly anisotropic and deform by a combination of brittle and ductile mechanisms, since they consist of many different phases with different strength. Therefore, it is difficult to establish a constitutive equation that describes the macroscopic mechanical behavior based on microphysical deformation mechanisms. Commonly, the total axial strain is described by adding the elastic strain and the inelastic strain in each of the three different creep regimes separately (Gao et al. 2010;Brantut et al. 2012Brantut et al. , 2013Brantut et al. , 2014aDewers et al. 2017). This approach was successfully applied on monophase materials (e.g., metals, Chindam et al. 2013) and may serve as a first approach also to describe the mechanical response of polyphase shale rocks, although the semibrittle deformation behavior of shales suggests that several mechanisms operate in parallel, probably with a dominant process acting in each regime.
Depending on the applied deformation conditions (p c , T, σ) and composition of the investigated shale rocks, some samples showed only primary creep (group 1), while others deformed by primary and subsequently secondary creep and reached tertiary creep (group 2) (cf., Figs. 3, 4). It is conceivable that different deformation mechanisms may have been active in the samples, such as plastic dislocation activity and brittle microcracking. By definition, primary creep occurs at continuously decreasing axial creep strain rate, ̇c reep , whereas secondary and tertiary creep are characterized by a constant (minimum), ̇m in , and increasing strain rate, respectively (e.g., Brantut et al. 2014a, b;Rybacki et al. 2017). To distinguish samples belonging to group 1 or group 2, we calculated the strain rate as a function of time for every experiment. One example is given in Fig. 7, showing the change of strain rate with time for the two Posidonia (HAR) and Bowland (OC) samples (HAR43 and BOS16), as plotted in Fig. 3b. Obviously, the quartz-rich Bowland shale (BOS_OC-Qtz) deformed solely by primary creep (group 1), displaying a continuously decreasing strain rate over time with progressive hardening. In contrast, the weak Posidonia shale (POS_HAR-Cly) exhibits first a decrease of strain rate with time, but subsequently an increase after passing a minimum value, which is characteristic for the three phases of primary, secondary, and tertiary creep (group 2). For samples of group 2, we determined the minimum strain rate (Table 2).
Our microstructural observations of deformed Posidonia (HAR) and Bowland (OC) shales suggest that primary creep is mainly accommodated by deformation (e.g., bending) of weak material (Figs. 5f, 6b, d), and partly by local pore space reduction (Fig. 6c). We assume that plastic flow of organic matter and shearing of clay flakes as well as intracrystalline dislocation glide of phyllosilicates contributes to viscous creep. This may be combined with subcritical crack growth via chemical reactions (stress corrosion) during the secondary creep regime (Rybacki et al. 2017). Samples exhibiting tertiary creep show strongly localized deformation close to a macroscopic shear fracture. Here, deformation is mainly assisted by inter-and intracrystalline microcracking as well as crack coalescence.
Empirical models in the form of parabolic, exponential, logarithmic, or hyperbolic dependence of strain on time have frequently been used to describe the decelerating creep phase (group 1) of rocks, soils, and metals (Gupta 1975;Findley et al. 1976;Karato 2008;Paterson 2013). If the strain depends linearly on stress, e.g., in many cases for small stress levels, linear viscoelastic models can be used to describe the (primary) creep of a polyphase material Zoback 2013b, 2014). The basic elements of linear viscoelastic models (elastic spring and viscous dashpot, each of them representing the endmember phenomenon of elastic and viscous properties of a specific mineral) may be combined to describe the constitutive behavior of the rock (Li and Ghassemi 2012;Almasoodi et al. 2014). However, this requires knowledge of the elastic and viscous properties of each phase under the applied deformation (p c , T, σ) conditions that are typically not available. In addition, these models cannot capture the volumetric and spatial distribution of the different phases, or brittle fracturing of some constituents. For stress corrosion, Brantut et al. (2013) recently suggested a model using an exponential law to approximate crack-damage-related primary creep of rocks.
The secondary creep phase of rocks (group 2) at high confining pressure and temperature has been characterized by quasi-constant strain rate over time (indicating a balance between the rates of strain hardening and softening), where deformation may be accommodated by dislocation and diffusion activity or other ductile processes, which is commonly described by flow models based on a power law (Evans and Kohlstedt 1995;Dorner et al. 2014;Rybacki et al. 2017). Brantut et al. (2012) and Heap et al. (2009) established a similar power law creep relation between strain rate, stress, and confining pressure in the secondary creep regime, where they assume that the dominating deformation mechanism is subcritical microcrack growth (stress corrosion). The crack velocity of this thermally activated mechanism depends on several factors such as temperature, humidity, and stress intensity factor at the crack tip (Kranz 1980;Kranz et al. 1982;Atkinson 1987;Ciccotti 2009;Heap et al. 2009). Although our samples were dried at 50 °C for at least 48 h, fluids (clay-bound water) probably remain within unconnected pores (see Sample Material section above), potentially leading to time-dependent creep strain and static fatigue of the investigated samples due to subcritical crack growth by stress corrosion (Brantut et al. 2012(Brantut et al. , 2014a. Unfortunately, the Paterson deformation apparatus did not allow us to measure the load-normal lateral strain, and thus the volumetric strain, allowing to detect microcrack activity. Instead, we measured the porosity of samples before and after deformation. Figure 8 shows the calculated porosity reduction, ϕ reduction = ϕ prior deformation − ϕ after deformation , where ϕ reduction > 0 represents pore space reduction and ϕ reduction < 0 indicates pore space increase. For primary creep, Fig. 8 suggests that an increase of confining pressure yields porosity reduction and an increasing stress results in a porosity increase, likely by microcrack generation. Posidonia (HAR) (Fig. 8a, c, e) as well as Bowland (OC) shale (Fig. 8b, d, f) specimens belonging to group 2 (post primary creep) exhibit the largest porosity increase after sample deformation, independent of applied confining pressure (Fig. 8a,  b), temperature (Fig. 8c, d), and stress (Fig. 10e, f). An increasing ϕ reduction with increasing confining pressure may be explained by Goetze's criterion, which suggests that dilatancy and microcrack formation may occur as long as σ > p c (Evans and Kohlstedt 1995). These observations suggest that group 2 samples may have accommodated part of the deformation by subcritical crack growth due to stress corrosion. For porous sandstone deformed at room temperature, Heap et al. (2009) suggested the onset of dilatancy already at stresses of about 80-90% of the triaxial compressive strength leading to secondary and eventually tertiary creep. Since the applied stresses in this study are close to the short-term triaxial compressive strength obtained under similar p c -T conditions (Herrmann et al. 2018), we expect at least a partial contribution of this mechanism to creep of our samples.
Alternatively, secondary creep can be regarded as a (longlasting) transient creep phenomenon at minimum strain rate, where the deformation mechanisms, which are responsible for primary and tertiary creep balance by recovery or recrystallization over a certain period (Zener and Hollomon 1946). Brantut et al. (2013), noted that creep driven by irreversible  Fig. 3b) in double-logarithmic scale. Sample BOS_OC-Qtz deformed by primary creep, indicated by a continuously decreasing strain rate with increasing time. The Posidonia (HAR) shale sample shows in addition a secondary and tertiary creep phase, for which the strain rate reaches minimum value (secondary creep) before it increased with time (tertiary creep). Note that the strain rate data are smoothed because of numerical noise 1 3 crack growth is not a steady-state process and may be more likely an inflexion between the primary and tertiary creep phase, which exhibits a specific minimum creep strain rate.
The following discussion will focus on the influence of composition and boundary conditions (σ, T, p c ) on the primary creep phase of shale rocks, since most of the deformed shale specimens exhibited this creep behavior. An attempt to  c (a, b), T (c, d) and σ (e, f). ϕ reduction > 0 represents pore space reduction, whereas ϕ reduction < 0 displays increasing porosity due to deformation. Primary = samples exhibiting only primary creep behavior, post primary = samples exhibiting in addition to primary also secondary and partly tertiary creep behavior. Deformation conditions are indicated quantify the influence of σ, T, and p c on the secondary creep phase of the investigated shales is given in the appendix, revealing large uncertainties of the calculated parameters.

Effect of Confining Pressure, Temperature, and Stress on Primary Creep of Shale Rocks
To discuss the influence of confining pressure, temperature, and stress on the primary creep phase of Posidonia (HAR) and Bowland (OC) shales, we use a simple phenomenological power law approach: where ε t = total axial strain, ε el = elastic axial strain, ε inel = inelastic axial strain, σ = axial differential stress, E = static Young's modulus, t = experimental time, and 'a' and 'b' are constants. As suggested by Rybacki et al. (2017), the following empirical correlations were used to account for the influence of p c , T, and σ on the constants 'a' and 'b' used in Eq. (1): yielding where a 0 and b 0 are rock constants, n a,b = stress exponent, Q a,b = activation energy, V a,b = activation volume, R = molar gas constant, and T abs = absolute temperature. This approach was adopted from constitutive equations describing hightemperature creep of rocks, where the stress dependence is expressed by the stress exponent, the temperature sensitivity by the activation energy, and the confining pressure dependence by the activation volume. Note, however, that these correlations are empirical and do not claim to correctly express the underlying microphysical processes assumed here such as dislocation glide and bending and shearing of phyllosilicates, compaction, or granular flow.
The amount of elastic strain in Eq. (1) may be determined by calculating the ratio of applied stress and static Young's modulus. We used Young's modulus data from constant strain rate experiments previously performed on similar shale samples at similar p c and T conditions (Herrmann et al. 2018). Although the authors measured a confining pressure dependence of E, we assume no significant change of elastic axial creep strain with varying p c -T conditions in our experiments, since the changes in the absolute values of E are relatively low in the range of the confining pressures applied in the experiments reported by Herrmann et al. (2018). The inelastic axial creep strain of each experiment was calculated by subtracting the elastic strain from the total measured axial creep strain. Plotting the inelastic axial creep strain versus time and fitting the curve by non-linear regression yields the constants 'a' and 'b' for each experiment. The parameter 'a' varies between 0.0011 and 0.00639 s −b and 'b' between 0.03839 and 0.26815 (Table 2).
We also determined the parameters 'a' and 'b' of group 2 samples by fitting only the primary part of creep curves, which results in tremendously higher values compared to those obtained for group 1 samples. This indicates a difference in the relative contribution (from bending of minerals and pore closure towards generation and coalescence of intra-and inter-granular cracks) of the present deformation mechanisms, which act in the primary creep regime of both sample groups. Therefore, the approach given in Eq. (1) may not fully represent the time-dependent deformation behavior of shale rocks under all p c , T, and σ conditions, especially with respect to proppant embedment, where high stresses due to low contact areas between proppant agent and fracture surface are expected. However, since the lifetime of stimulated wells is typically only a couple of years, we assume that this approach is still useful to discuss the potential of a shale play for the economical extraction of hydrocarbons.
Taking the logarithm of Eqs. (2) and (3) yields the following correlations: which were used to determine the sensitivity of 'a' and 'b' to a change of p c , T, and σ. The corresponding plots to determine the parameters n, Q, and V for Posidonia (HAR) and Bowland (OC) shales are given in Fig. 9. In each graph, two of the three deformation conditions (p c , T, σ) are fixed to determine the sensitivity of parameters 'a' and 'b' to the remaining one. Linear regression fits for the parameters are shown by dashed lines in Fig. 9 and are summarized in Table 3, including the calculated constants a 0 and b 0 . They are discussed separately in the following three sections.
The equation suggested here to express the stress, temperature, and confining pressure dependence of primary creep of shale rocks is phenomenological. However, one may argue that the stress and temperature-sensitive creep parameter 'a' we obtained suggest that mainly viscous mechanisms

Effect of Axial Stress
The effect of stress on the primary creep behavior of Posidonia (HAR) and Bowland (OC) shale is expressed by the stress sensitivity n a and n b of the power law creep parameters 'a' and 'b', respectively. The stress exponent n a = 5.4 of weak Posidonia (POS_HAR) shale is almost three times higher than n a = 1.9 of quartz-rich (outcrop) Bowland (BOS_ OC) shale (Fig. 9a). Both values are distinctly higher then n a = 1, which is expected if the rocks would deform in a linear viscous manner. This non-linear viscous behavior may be caused by a range of mechanisms including dislocation glide or cracking. As the applied stresses of 160-316 MPa, are close to the triaxial compressive strength (Herrmann et al. 2018), we suggest that the generation of cracks and their interaction may also play a role at these conditions. Comparable results for Posidonia (DOT) shale (n a = 1.4), deformed at similar confining pressures and temperatures, were found by (Rybacki et al. 2017). In addition, Marcellus shale deformed at p c = 19 MPa, T = 20 °C and σ < 90 MPa exhibited non-Newtonian primary creep (Li and Ghassemi 2012). In contrast, other North American and Canadian shales such as Eagle Ford, Barnett, Haynesville and Fort St. John were found to deform in a linear viscous manner (Li and Ghassemi 2012;Zoback 2013a, b, 2014;Almasoodi et al. 2014;Rassouli and Zoback 2015). Note that in the latter studies, creep experiments were performed at low stresses (σ < 90 MPa, mostly < 45 MPa) and confining pressures (p c < 60 MPa), which are substantially lower than the applied stresses and resulting strains in our study. In addition, experiments reported in the literature were performed at ambient temperature only. Therefore, the expected deformation mechanisms are likely different for shales deformed at low and high stresses and temperatures, respectively. At lower stress conditions, the strain may be accommodated mainly by closing of preexisting microcracks and pores, whereas at high stress, generation (and subcritical growth) of microcracks, grain boundary sliding, and plastic deformation of weak mineral phases may have operated as well.
As observed by Rybacki et al. (2017) on Posidonia (DOT) shale, the stress dependence of the power law exponent 'b' (Eq. 6), n b , of Posidonia (HAR) as well as Bowland (OC) shale appears to be negligible (Fig. 9b).
A stepwise increasing stress irreversibly changes the microstructure of the samples within each single step, thereby affecting the deformation behavior of the subsequent creep steps (Fig. 2). Since most published creep experiments on shales were performed stepwise, we assume that the experimental protocol also affected the stress sensitivity of the investigated rocks (Brantut et al. 2013).

Effect of Temperature
With increasing temperature, we observed an increasing primary creep strain (cf. Fig. 4c, d). This may be due to weak phases with temperature-sensitive flow strength such as clays and organic matter (Mikhail and Guindy 1971). Least squares fitting of the temperature sensitivity of the creep parameter 'a' yields similar apparent activation energies of Q a = 14 ± 9 kJ/mol for Posidonia (HAR) shale and Q a = 13 ± 2 kJ/mol for Bowland (OC) shale (Fig. 9c). Regarding parameter 'b', however, shows no clear correlation with temperature (Q b ≡ 0 kJ/mol) for Posidonia (HAR) shale deformed at constant stress and confining pressure (Fig. 9d), suggesting minor temperature dependence. A 14.39 ± 8.49 kJ/mol -13.12 ± 1.55 kJ/mol - 1 3 best fit of the Bowland (OC) shale data (p c , σ = const.) gives an apparent negative activation energy of Q b ≈ − 5 kJ/ mol. Since this is counterintuitive, we also fixed Q b (BOS_ OC) ≡ 0 kJ/mol. This seems to be reasonable, since Bowland (OC) shale consists of ≈ 70 vol% quartz minerals, the strength of which is believed to be hardly sensitive to temperature over the temperatures and dry conditions tested. The calculated activation energies for primary creep of Posidonia (HAR) and Bowland (OC) shales are in relatively good agreement with estimates of Q a ≈ 3 kJ/mol and Q b ≈ 5 kJ/mol of other Posidonia (DOT) shales deformed in the primary regime (Rybacki et al. 2017). These values are lower than activation energies for steady-state creep of the weak phases (≈ 50-150 kJ/mol; cf., Rybacki et al. 2017;Herrmann et al. 2018) and for stress corrosion of sandstone and granite (≈ 30-50 kJ/mol; Brantut et al. 2012). The difference presumably results from the simultaneous operation of deformation mechanisms governing primary and secondary creep regime and from the multiphase composition of the investigated shales. Our determined apparent activation energies represent a bulk temperature sensitivity of primary creep of the shale samples.

Effect of Confining Pressure
With increasing confining pressure, we observed reduction in total axial strain (rate) for Posidonia (HAR) and Bowland (OC) shale deformed in the primary regime (Fig. 4a,  b), revealing material strengthening with increasing confining pressure, as expected for semibrittle deformation due to compaction of pores, closure of preexisting microcracks and enhanced friction between grains with increasing p c . This observation is in line with an increasing triaxial compressive strength with increasing confining pressures at constant strain rate as observed by (Herrmann et al. 2018) for the same types of shale rocks.
To empirically quantify the influence of confining pressure on the primary creep phase of Posidonia (HAR) and Bowland (OC) shales, we estimated an apparent activation volume, yielding V b (POS_HAR) = 31 ± 7 cm 3 /mol and V b (BOS_OC) = 12 ± 11 cm 3 /mol for the power law creep parameter 'b' (Fig. 9f). A higher value for Posidonia shale than for Bowland shale is likely related to the higher amount of weak sample constituents of Posidonia shale (Table 1). The power law parameter 'a' shows almost no pressure sensitivity and was, therefore, fixed to zero [V a (POS_ HAR) = V a (BOS_OC) ≡ 0 cm 3 /mol] (Fig. 9e), in agreement with results observed for Posidonia (DOT) shale (Rybacki et al. 2017). In contrast, Zoback 2013a, b, 2014) found no significant influence of confining pressures on the primary creep behavior of Barnett, Eagle Ford, Fort St. John and Haynesville shale. The authors observed also linear viscoelastic creep for their investigated shales, indicating that pressure-sensitive brittle deformation mechanisms were likely not the dominant mechanisms accommodating creep in their experiments. The low-pressure sensitivity may be explained by the relatively low confining pressures (p c = 10-60 MPa) and stresses (σ = 3-45 MPa) used by Zoback 2013a, b, 2014) which minimizes pore collapse and microcrack generation.

Effect of Sample Composition and Mechanical Properties
To investigate the influence of shale composition on primary creep on mechanical behavior during reservoir extraction, we estimated the total primary creep strain after 3 year deformation, during which most of the production decline in shale gas reservoirs occurs. For comparison with literature data on Posidonia (DOT) shale (Rybacki et al. 2017) as well as on Barnett, Haynesville and Eagleford shale (Sone and Zoback 2014), we used a modified form of Eq. (1) by replacing the power law parameter 'a' by a = A * n a , yielding For the determination of ε t after 3 year deformation, we recalculated the values A and b for deformation conditions of p c = 20 MPa, T = 20 °C, and σ = 25 MPa to allow comparison with the experiments on N-American shales that were performed at low stresses and confining pressures and at ambient temperature. Results are given in Table 4 together with the approximate composition and porosity. Note that Sone and Zoback (2014) and Rybacki et al. (2017) fitted measured creep curves of total axial strain, ε t . Therefore, no elastic strain was added to the strain values determined from the literature. For comparison, we additionally calculated the elastic, ε el , and inelastic strain values, ε inel (Table 4), using the data provided by Sone and Zoback (2014).
The resulting total strains (in %) are shown in a tertiary diagram (Fig. 10a) to visualize the influence of sample mineralogy on the time-dependent creep behavior. As shown in Fig. 1, sample mineralogy is separated into strong (QFP), intermediate strong (Cb), and weak (Clay + TOC + Mica + ϕ) components. No correlation between ε t and the amount of QFP, carbonates, or weak phases is found (Fig. 10a, Table 4). There may be a negative correlation between ε inel and QFP (cf., Table 4), if shale rocks are recovered from different formations, although an inverse trend was found for shale samples acquired from the same formation. However, it is important to recall that the creep mechanisms of the various shales are probably different due to the different p c , σ, and T conditions applied in the different studies. Separating the influence of porosity on creep deformation also shows no trend (Table 4).

Table 4
Calculated primary creep strains after 3 years Subscripts I and II indicate clay and organic-rich (I) and organic-poor (II) rocks, respectively Ref. The correlation between extrapolated total strain after 3 year deformation and the static Young's modulus (Herrmann et al. 2018;Zoback 2013a, b, 2014) is shown in Fig. 10b. The axial total strain decreases with increasing static Young's moduli (Fig. 10b) and may be calculated using the correlation given in the figure (dashed line). Note that the correlation found in Fig. 10b is strongly influenced by the elastic strain (black line in Fig. 10b). Considering only the inelastic strain would yield a less pronounced trend (cf., Table 4). However, since the Young's modulus can be easily measured by wireline logs, this correlation may be especially interesting for practical purposes to estimate the long-term creep behavior of shale rocks. Note that this correlation is only applicable at depth levels representing the applied deformation conditions and will change at greater depth with high pressure and temperature.
Total strains are also smaller for strong samples with high triaxial compressive strength, σ TCS (Table 4). This is in line with the positive correlation between σ TCS and E of shale rocks (Herrmann et al. 2018).
In addition to the previously mentioned parameters (composition, E, σ TCS , and B), the angle between the maximum principal stress orientation and bedding orientation will also influence the primary creep strain due to the anisotropic character of shale rocks (Swan et al. 1989;Villamor Lora et al. 2016). If loaded parallel instead of perpendicular to bedding, samples would creep at slower rate yielding less creep strain after 3 years deformation (cf., Sone and Zoback 2014). In addition, the degree of anisotropy plays an important role in the deformation behavior of shale rocks, since bedding perpendicular and bedding parallel loading represent only endmembers, whereas in nature, an intermediate state (of a load-bearing framework of hard mineral phases or interconnected weak layers) is more likely, as pointed out by (Sone and Zoback 2013b).
The presence of water within samples also affects the creep behavior of shales in a way that creep is enhanced for samples with higher water content as shown for Haynesville shale by (Sone and Zoback 2014) and Alum shale by (Rybacki et al. 2017). This is presumably due to swelling of clay minerals or microcrack growth due to stress corrosion at higher stresses.

Comparison with Rheological Bodies
In addition to the empirical approach given in Sect. 5.2 to characterize the primary creep behavior of shale rocks, we also applied a linear viscoelastic model to fit the recorded primary creep data of shales, as has been suggested recently (Li and Ghassemi 2012;Almasoodi et al. 2014). These models typically consist of a combination of rheological bodies such as springs (elastic) and dashpots (viscous). Note that these models ignore effects of changing confining pressure and temperature. In addition, the models assume a linear correlation between strain and applied stress (n = 1), which was found to be different for the investigated shale rocks of these study. Irrespective of these limitations, we fitted a Zener model in series with a Kelvin-Voigt model, as proposed by (Almasoodi et al. 2014), to our recorded primary creep curves using the following equation: where µ = viscosity and subscripts 1 and 2 represent different sample constituents. This model was chosen, since it is in good agreement with the approach of separating shale Fig. 10 Calculated primary total strain (in %) after 3 years deformation at 20 MPa confining pressure, 20 °C temperature and 25 MPa stress versus composition (a) and static Young's modulus, E, (b). QFP quartz + feldspar + pyrite, Cb carbonates, TOC total organic carbon, ϕ porosity composition into compliant (cly, TOC, ϕ) and stiff (QFP, Cb, Mca) components, as suggested by Herrmann et al. (2018) for the same shale samples. Mca and Cb are regarded as stiff mineral phases here, since they display relatively large Young's moduli (Mavko et al. 2009). Stiff and compliant components are represented by subscripts 1 and 2, respectively. Fitting this model to the recorded creep curves, using the non-linear Levenberg-Marquardt algorithm yielded strongly unstable results for Young's moduli and viscosities both for Posidonia and Bowland shale. Therefore, the usage of such models based on springs and dashpots to characterize the time-dependent primary creep behavior of shale rocks may only be feasible at low stresses in the Newtonian viscous regime. For non-linear viscous creep or if brittle mechanisms such as the initiation and coalescence of microcracks contribute to the overall deformation, these models may not be appropriate.

Mechanical Equation of State
In addition to the phenomenological approaches given in Eqs. (4) and (8), the primary creep behavior of shale rocks may be described by a constitutive model based on a mechanical equation of state, as has been previously applied on metals, Carrara marble, and salt rocks (Zener and Hollomon 1946;Hollomon 1947;Hart 1970Hart , 1976Stone 1991;Covey-Crump 1994, 1998, 2001Stone and Plookphol 2004;Evans 2005). Assuming homogeneous deformation of an isotropic material with strain hardening, this approach aims to reliably extrapolate results of laboratory experiments to geological time scales (Hart 1970;Covey-Crump 1994).
Although not exactly representing the approach given by the authors, we like to introduce an exponential term, since this allows data fits at much higher quality. To describe our data, we used the following constitutive formulation to describe the strain(rate)-dependence of stress: where the parameter may be described as K = strength parameter, ε inel = inelastic axial strain, o = strain-hardening exponent, ̇0 = inelastic strain rate, and p = strain rate sensitivity. Again, in this equation, we assume an exponential dependence on strain, which fits much better to our data than the commonly assumed linear dependence (e.g., Karato 2008;Kassner and Smith 2014, cf., Eq. 10). Further assuming that the parameter p depends solely on stress to avoid overparameterization, we determined this parameter by plotting ln (σ) over ln ( ̇0 ) at given inelastic strains, ε inel , of those experiments, where p c and T were held constant and only σ was varied. This procedure yields p (POS_ HAR) = 0.027 ± 0.001 and p (BOS_OC) = 0.0395 ± 0.0025, respectively. The strength parameter K and strain-hardening  (Table 5). The parameters 'K' and 'o' appear to vary linearly with confining pressure and temperature ( Table 5). The corresponding contours are plotted in Fig. 11a, b and may be used to account for the influence of p c and T on these parameters of Posidonia (HAR) and Bowland (OC) shale at a certain depth using appropriate temperature and pressure gradients.
The inelastic strain rate as a function of inelastic strain calculated from Eq. (9) is compared in Fig. 11c, d with the values measured for two experiments (HAR33 and BOS26), suggesting relatively good agreement for both Posidonia (HAR) and Bowland (OC) shale. These samples were chosen, since they also represent experiments conducted under constant stress for the comparison of primary creep and constant strain rate behavior, as can be found in Sect. 5.2.7. Note that at small inelastic axial strains, corresponding to t < 30 min, calculated and measured strain rates do not fully coincide, possibly due deformation mechanism that are not covered by the used mechanical equation of state, such as pore collapse or closure and closing of pre-existing microfractures.

Correlation Between Primary Creep and Constant Strain Rate Behavior
As suggested by Kassner and Smith (2014), we used a slightly modified form of Eq. (9) to correlate between results obtained from constant stress and constant strain rate experiments: The strain rate sensitivity, p, was set p ≡ 0.01 based on experiments deformed at varying strain rates in the range of 5 × 10 −6 to 5 × 10 −4 s −1 (Herrmann et al. 2018 (OC) shale yields under the assumption of ε inel = 0.01 and σ = 314 MPa (taken from constant stress experiment BOS26) the following rates: ̇0 = 6.7 × 10 −8 s −1 for K, o determined from constant strain rate experiment, and ̇0 = 1.0 × 10 −7 s −1 for K, o determined from constant stress experiment. For Posidonia shale (ε inel = 0.01, σ = 194 MPa), we obtained ̇0 = 4.9 × 10 −9 s −1 and = 8.4 × 10 −6 s −1 , respectively. Surprisingly, the variations are small for Bowland shale but rather large for Posidonia shale. This may be caused by the slight difference of the p c -T deformation conditions for the selected Posidonia (HAR) samples deformed under constant strain rate and constant stress, or by sample-to-sample variations.
Another approach to compare the constant strain rate tests with our constant stress tests may be given by first rearranging Eq. (10) and replacing ̇0 by inel t that yields Integrating Eq. (11) and assuming σ = const. (creep experiment) results in the following equation: The similarity of Eq. (12) and Eq. (1) allows to derive the creep parameters, 'a' and 'b', from constant strain rate experiments using the following equations: Setting σ = 194 MPa (from HAR33 creep experiment) for Posidonia (HAR) shale, we obtained a (POS_HAR) = (1.92 ± 0.08) × 10 −3 and b (POS_ HAR) = 0.137 ± 0.002 from the performed constant strain rate test using Eqs. (13) and (14). Extrapolated to strain achieved after 3 year deformation based on Eq. (1) yields a total strain of ε t = 0.045 ± 0.002. This value is somewhat lower than the total strain of ε t = 0.08, extrapolated from the corresponding creep experiment.
For Bowland shale (σ = 314 MPa, BOS26), we calculated the creep parameters a (BOS_OC) = (1.74 ± 0.07) ×10 −3 and b (BOS_OC) = 0.172 ± 0.003. Using these values, we obtain a total axial creep strain of ε t = 0.064 ± 0.004 after 3 year deformation, which is relatively close to the extrapolated total strain of ε t = 0.04 based on creep testing. The differences between extrapolated total axial strains from constant strain rate and constant stress tests may arise from slightly different deformation conditions or sample-to-sample variations. In addition, the approach is based on a mechanical equation of state, which assumes that the stress required for inelastic deformation depends only on the instantaneous values of strain, strain rate, temperature, and confining pressure and not on the loading history, which is in contrast to our observations (Fig. 2). However, the obtained results suggest that quantitative estimates of the (long-term) creep strain of shale rocks from short-term constant strain rate experiments may be feasible, but only as a first order approximation.

Conclusions
Constant stress experiments performed on Posidonia (HAR) and Bowland (OC) shale reveal time-dependent creep behavior characteristic of semibrittle deformation, influenced by confining pressure, temperature, and applied stress. At high stress and temperature, and at low confinement samples displayed primary creep also secondary and tertiary creep prior to sample failure. Primary creep strain is mainly accommodated by the deformation of weak sample constituents (clays, TOC, and mica) and local pore space reduction. At the onset of tertiary creep, samples displayed localized deformation by formation of a single macrocrack, which is composed of coalesced subparallel microcracks for Bowland (OC) shale. Even in these high stress samples, the dislocation activity in quartz and carbonates is low. The effect of deformation conditions on the creep behavior is more pronounced for weak Posidonia (HAR) shale than for strong Bowland (OC) shale. Empirical power law correlations are applicable to account for the influence of stress, confining pressure and temperature on the primary creep phase of shale rocks.
The influence of sample composition on the creep behavior is in the same order of magnitude as the effect of applied deformation conditions. A negative correlation between primary creep strain and sample strength and static Young's modulus is evident. These results are useful to assess the potential of a reservoir for economical and extended extraction of hydrocarbons, since the Young's modulus can be easily and quickly estimated from wireline logs.
A correlation between the primary creep and constant strain rate behavior of shale rocks was found to some extent if calculations are based on a mechanical equation of state. This may be useful to estimate the long-term creep behavior of shale rocks from their short-term properties.
With respect to proppant embedment, our results suggest a lower fracture closure rate of strong Upper Bowland (OC) shales compared to weak Posidonia (HAR) shale constant strain rate represents the range, which was used by the authors (Fig. 5 in Brantut et al. 2014a, b) to calculate the stress deficit, ∆Q = σ creep test − σ const. strain rate test , for every single inelastic strain value. In our tests, the stress-inelastic axial strain curves did not intersect each other twice, since the sample deformed at constant stress failed earlier (Fig. 12a). This may be due to the slightly different temperatures, since axial strain increases with increasing temperature (cf. , Fig 4c, d). Therefore, for the calculation of ∆Q, we used the area between the strain at the first intersection of recorded curves and the strain, where tertiary creep started (dashed grey area in Fig. 12a). Brantut et al. (2014a, b) proposed the following correlation between creep strain rate, ̇c reep , at any given value of inelastic strain during constant stress deformation and ∆Q: where ̇0 = applied constant strain rate during constant strain rate test and σ* = characteristic activation stress, which is a combination of parameters describing subcritical crack growth (Brantut et al. 2014a, b). Figure 12b shows the natural logarithm of the creep strain rate normalized by the applied constant strain rate as a function of ∆Q. In line with observations made by (Brantut et al. 2014a, b), the resulting curve has different slopes for the decreasing ̇c reep (↓) and the increasing ̇c reep (↑) with a minimum ∆Q at minimum ̇c reep . Fitting Eq. (17) to Posidonia (HAR) shale, data yielded σ* = 6.4 ± 2.7 MPa (Fig. 12b), which is larger than what was found for different sandstones (σ* = 1-2.7 MPa) by (Brantut et al. 2014a, (17) creeṗ 0 ≈ e ΔQ * , b). The contrasting results may originate from the different deformation conditions in our experiments and those performed by (Brantut et al. 2014a, b) and the different rock types. In addition, during the deformation of Posidonia (HAR) shale, other deformation mechanisms may be active, but are not considered in the formalism above.
Interestingly, this approach was not applicable on Bowland (OC) shale, although it contains ≈ 70 vol% QFP and, therefore, should be better comparable to sandstones than Posidonia (HAR) shale. The applied constant stress, which is necessary to reach the secondary creep phase, is slightly larger than the recorded stress during performed constant strain rate test. This may be due to the relatively high applied constant strain rate, at which the shale sample cannot accommodate deformation over a certain time period, leading to lower strength due to increased pore pressures as a result of pore collapse, as all experiments have been performed at undrained conditions. This may be more pronounced for porous Bowland (OC) than for low porous Posidonia (HAR) shale (cf., Table 1). Therefore, no ∆Q could be calculated.
However, the obtained results may be useful to subsequently assess specific rock properties of some shales, such as the stress intensity factor, time to failure as well as creep strain rates by only performing one constant stress and one constant strain rate test at the same p c -T conditions. Fig. 12 a Stress-inelastic axial strain curve of two Posidonia (HAR) samples deformed at either constant strain rate, ̇0 , or (quasi) constant stress, σ, recorded at slightly different temperatures. Due to the assumption of constant volume deformation (yielding an increasing sample diameter with increasing strain), σ of constant stress test is slightly decreasing, since a constant load was applied. Grey dashed area was used to calculate stress deficit, ∆Q. b ln of normalized creep strain rate versus stress deficit of the two Posidonia (HAR) samples shown in a. σ* = activation stress. Deformation conditions are indicated