Viscoelastic Strains of Palaeozoic Shales under the Burger’s Model Description

Abstract

This article presents the results of creep studies of Palaeozoic shales from the Baltic Basin in which the exploitation of shale gas in Poland was planned. Knowledge of instantaneous and long-term properties investigated in triaxial stress conditions is important from the point of view of exploitation techniques related to hydraulic fracturing. Rheological phenomena also play an important role in the analysis of the initial stress in shales, the knowledge of which is indispensable in the hydraulic fracturing process. The tests were carried out on samples representing four siltstone–claystone lithostratigraphic units occurring in the Baltic Basin. The studies and analyses were aimed at determining the character of creep in shales, selection of the appropriate rheological model for the analyzed rocks, and determination of the threshold of the linear creep under triaxial compression conditions. An original approach together with analysis results are presented here, which enable the separation and monitoring of shear and volume creep effects, and on this basis, the determination of the significance of the contribution of volume creep in the entire creep process. A relatively simple methodology for determination of the parameters of the Burgers model using this division is presented. The original value of the article is also due to the test results themselves and the parameter values of the analyzed model for triaxial creep of shales, which are not numerous in the literature. The investigations were performed at various loading levels in relation to the triaxial strength of the shales. Depending on the load, at its low values up to 0.7 (σ₁ − σ₃)_max, creep had a determined character and did not show features of progressive creep. The linear creep threshold was also analyzed in this range. The loading level of 0.7 (σ₁ − σ₃)_max was the limit of linear creep. Exceeding this load resulted in the loss of the linear character of creep, which in consequence lead to the subsequent third creep phase ending with rock damage. Parameters of the Burger’s model for gas shales from the Baltic Basin (northern Poland) were identified. There are significant differences in the behavior of shales depending on the lithostratigraphic unit from which the samples were collected. The mineral composition of the shales also influenced their behavior.

1. Introduction

All rocks display rheological features, which are emphasized to a larger degree in the case of soft rocks with high porosity rather than high strength rocks [1]. Various rock types have different creep predispositions [2]. Strain caused by creep may take place for several minutes or may last for several years. Thus, the study of creep is not only restricted to oil engineering but is also vital in long-term rock engineering and underground dynamics in mines, tunnels, and underground storage of nuclear waste, CO₂, and natural gas [3]. Salt represents a medium well recognized with regard to its rheology [4,5,6,7]; in addition to its exploitation (as potassium salt), salt deposits are used as magazines for gas or hydrocarbons. In the case of strong rocks, rheological phenomena are analyzed less frequently, but likewise in this case, due to the source and cause of creep, it is important to determine the time at which the rock is damaged or the moment when progressive creep begins. Gas shales are specific heterogeneous rocks [8,9,10], with a larger or smaller anisotropy [11,12,13]. Their rheological features have been analyzed with regard to the selection of the creep model [14,15], the influence of the mineral composition on creep [16,17], and the possibility of applying micro-indentation [18,19,20] to determine the rheological properties in shales. Due to the relatively complex description and complexity of tests, anisotropic creep of shales has been rarely studied [21]. Models taking into account creep effects are not frequently used in the simulations of hydrofracturing in shales or other rocks [15,22]; however, these phenomena may be observed during hydrofracturing [23,24]. Some attempts of the usage of rheological models in simulation fracturing may be found in [25,26], but such reports are rare, and their application is restricted. Rheological phenomena play an important role in the analysis of the initial stress state in shales, an issue that is indispensable in hydrofracturing [27,28]. An important but rarely studied aspect is the significance of volumetric change in the creep process [18,27,29,30]. This kind of change records information about the creep mechanism. This article analyses volumetric change related to creep, with shear and volumetric creep strains treated separately. Moreover, the significance of volumetric creep description is determined here for Paleozoic shales from the Baltic Basin. After creep character was initially recognized, the Burger’s model was selected for its description. The model was defined for triaxial loading conditions. The investigations and analyses were performed at different loading levels with regard to the triaxial shale strength. The linear limit of the creep process, i.e., the limit of application of the viscoelastic physical model, was determined for the analyzed rocks.

2. Characteristics of the Material

The investigated shale samples were collected from depth interval 3600–4000 m b.s.l. from the Baltic Basin (northern Poland). They represented various lithostratigraphic units from the Upper Ordovician (Sasino Claystone Formation) to the lower Silurian (Pelplin Claystone Formation, Pasłęk Claystone Formation, Jantar Bituminous Claystone Member). These units correspond to the interval between the Caradocian and the Wenlock [31]. The mineral composition of the analyzed rocks was determined using X-ray diffraction (XRD). Samples from the Pelplin and Pasłęk formations were characterized by a homogenous mineral composition.

Samples from the Pelplin Formation contain, ca., 46% of clay minerals, 40% quartz, feldspars and pyrite (QFP), and 14% carbonates, while samples from the Pasłęk For-mation are characterized by an elevated value of clay minerals (57.5%), and a lower content of both QFP minerals (36%) and carbonates (6.5%). Samples from the Jantar Member and Sasino Formation had a complex mineral composition. The Jantar Shales had 47% of clay minerals, 33% QFP, and 20% carbonates. The Sasino Shales are more siliceous, with an average to 45% clay minerals, 52% QFP, and only 3.5% carbonates. As indicated in [32], the mechanical and deformation properties of rocks from different units differ significantly. The results of the performed XRD analyses and petrophysical properties of the analyzed rocks were presented in detail in [32,33,34].

3. Laboratory Rheological Analyses of Paleozoic Rocks

3.1. Experimental Equipment

Creep analyses were performed in a servo-controlled strength press MTS-815 type equipped with a triaxial cell. Cylindrical samples cut out in a direction parallel to lamination (horizontal samples) with a diameter of 38 mm and height of 76 mm and height to diameter ratio of 2 were prepared for the laboratory tests (Figure 1a). Displacement of samples was measured using extensometers assembled on the samples; as a result, axial (ε_z) and circumference (ε_x,y) strain were determined (Figure 1b). All tests were performed at permanent temperature of 85 °C and a permanent confined pressure P_conf = 50 MPa. The confined pressure (P_conf = σ₂ = σ₃) was set with a permanent velocity of 10 MPa/min. An example of the sample CTC loading path is shown in Figure 1c. A detailed description of the research devices can be found in [34].

The testing procedure, according to which the samples were analyzed, is composed of three stages. At first, the temperature in the triaxial cell was increased to reach 85 °C and maintained at a constant level to the end of this stage. In the second stage, the confined pressure P_conf = σ₁ = σ₂ = σ₃ was increased with a permanent velocity of 10 MPa/min. to 50 MPa in three steps, to 10 MPa, 30 MPa, and finally 50 MPa; such pressure was maintained at a constant level until the end of the test. The confined pressure (P_conf = σ₂ = σ₃) was increased with a constant velocity of 10 MPa/min. The third essential loading stage usually comprised two differential loading levels, ca., 25% (σ₁ − σ₃)_max and, ca., 50% (σ₁ − σ₃)_max, during which shale creep was analyzed. The first level of, ca., 25% (σ₁ − σ₃)_max was usually maintained for 5 h, and the second from 5 to 20 h. In the case of two samples, creep was analyzed also at the level of, ca., 75% (σ₁ − σ₃)_max for 10 h.

In the last stage, the differential stress was increased to complete destruction of the samples with a constant velocity of axial strain at 10⁻⁵ (s⁻¹). The results presented in this report refer to samples and procedures determined by [34] as viscoelastic.

3.2. Initial Creep Analysis in Paleozoic Shales on 1D Models

Rock creep is often described by 1D models, which in the case of both uniaxial and triaxial stresses show the dependency of axial strain on time. Such analyses allow the definition of the character of creep and the determination of the correctness of the initial assumptions in the creep model for the analyzed rocks. Figure 2 shows an example of creep strain obtained for a shale sample from the Pelplin Formation (W260A, see

Table 1. Parameters of the 3D Burger’s model for the Paleozoic shales of the Baltic Basin.

		Deviatoric Creep						Volumetric Creep
Sample	Sample	(σ1 − σ3)	GK	GM	${\mathsf{\eta }}_{\mathit{s}}^{\mathit{K}}$	${\mathsf{\eta }}_{\mathit{s}}^{\mathit{M}}$	R2	KK	KM	${\mathsf{\eta }}_{\mathit{m}}^{\mathit{K}}$	${\mathsf{\eta }}_{\mathit{m}}^{\mathit{M}}$	R2
Location	Number	${\overline{\mathsf{(}{\mathsf{\sigma }}_1-{\mathsf{\sigma }}_3\mathsf{)}}}_{\mathit{m}\mathit{a}\mathit{x}}$	[GPa]	[GPa]	[GPa/day]	[GPa/day]		[GPa]	[GPa]	[GPa/day]	[GPa/day]
Pelplin	W260A	0.23	473	25	8.0	292	0.89	420	40	38	48	0.87
		0.46	237	19	7.6	167	0.98	315	35	31	94	0.90
		0.69	81	16	3.0	63	0.99	221	39	16	65	0.89
	W264A	0.27	242	20	4.7	192	0.96	203	29	61	47	0.93
		0.55	166	20	3.0	63	0.99	413	46	38	70	0.90
	B211A	0.33	400	18	7.1	205	0.91	359	35	32	36	0.90
		0.65	133	18	2.2	86	0.99	231	62	57	68	0.91
	M30B	0.26	282	18	4.1	166	0.95	130	32	46	55	0.83
		0.52	118	15	2.7	60	0.98	281	35	6	67	0.87
	W260B	0.26	268	19	5.1	107	0.98	269	25	19	42	0.92
		0.52	128	17	6.1	145	0.98	174	37	7	59	0.89
Pasłęk	W278A	0.28	372	20	7.4	224	0.90	450	38	67	33	0.90
		0.56	162	17	4.7	191	0.97	178	28	18	66	0.92
		0.84	15	13	1.1	20	0.99	331	60	56	81	0.95
	W278B	0.34	185	16	3.7	194	0.88	439	38	28	37	0.85
		0.69	20	12	0.6	43	0.98	485	39	19	59	0.82
	B236A	0.20	711	24	9.1	195	0.96	363	27	27	114	0.92
		0.41	257	21	5.3	138	0.90	61	51	16	117	0.88
	W277B	0.51	68	17	1.1	71	0.92	173	39	19	81	0.90
Jantar	W305A	0.25	244	18	1.8	76	0.96	6	24	10	61	0.78
		0.49	123	14	1.3	32	0.98	226	34	10	45	0.89
	W306A	0.26	260	19	4.6	104	0.97	357	23	22	53	0.91
		0.52	110	14	3.3	80	0.99	182	37	22	73	0.94
		0.78	42	10	1.3	31	0.99	175	48	19	59	0.95
	B275B	0.26	246	15	4.1	155	0.88	287	20	21	30	0.90
		0.52	98	12	1.6	34	0.94	150	31	7	43	0.90
	B279D	0.23	329	19	5.7	119	0.96	226	26	44	164	0.73
		0.46	145	17	3.0	67	0.98	464	45	77	50	0.81
	M43BB	0.23	285	20	2.8	110	0.95	90	26	22	50	0.83
		0.46	142	16	3.0	60	0.90	261	33	6	62	0.78
Sasino	M60B	0.30	336	18	7.1	105	0.92	378	22	22	44	0.90
		0.60	169	16	3.1	62	0.98	191	28	28	12	0.92
	W320B	0.21	314	18	1.9	185	0.92	106	22	44	39	0.82
		0.42	161	15	3.2	144	0.98	87	29	12	43	0.84
	M62A	0.24	319	22	3.7	257	0.88	116	24	30	73	0.64
		0.48	304	20	7.4	219	0.94	254	35	1	47	0.66

). The sample was analyzed on three levels of differential load applied at 25, 50, and 75% of (σ₁ − σ₃)_max.

As shown in the case of sample W260A for all loading ranges, creep has classical characteristics, from initial (instantaneous) elastic strain, through primary creep, and ending with secondary creep (Figure 3) [35,36,37]. Notable is that the largest loading, assessed at the level of 75% (σ₁ − σ₃)_max before the experiments, did not cause a distinct transition to the third progressive creep in sample W260A. A similar case was observed in most tests, where creep ended at a constant creep velocity (secondary creep). Applying the most popular 1D creep models, i.e., standard, Burger’s and exponential models (Figure 4) for creep description, it can be initially assessed, which one corresponds well to the creep in Paleozoic rocks from the Baltic Basin in the analyzed loading range.

The standard model has a restricted application for the description of creep in the analyzed shales because creep in this model approaches a constant strain value, which usually is not observed in rock creep. This is also the case in the analyzed shales (Figure 2). In Burger’s model, creep approaches an asymptote described by the following equation:

$${\sigma }_0\left[\frac{1}{E_1}+\frac{1}{E_2}+\frac{t}{{\eta }_1}\right],$$

(1)

In advanced creep at a constant velocity—secondary creep—the inclination of the creep strain can be determined from the relationship $\frac{{\sigma }_0}{{\eta }_1}$.

This phenomenon is observed in most medium and high strength rocks [34,39,40]. It is also noted in analyzed shales (Figure 2). The advantage of the presented models is the rather simple physical interpretation of their constants, which may be used directly in the construction of constitutive relations of viscoelastic media. The exponential model used in rock creep description [37] is the most universal regarding the accuracy of creep description. Despite the fact that constants B and n may have a specific physical interpretation [27], they do not directly describe the constants of viscosity or elasticity as structural models (e.g., Burger’s). Due to good match to the obtained results, the exponential model is often used to determine the creep function J(t) and creep compliance below [15,27,37]:

$$J\left(t\right)=B{t}^n,$$

(2)

where B and n are constants.

When analyzing the characteristics of creep in the representative sample W260A on three loading levels (Figure 2b), it can be assumed that both the Burger’s and Norton’s models should correctly describe creep of the shale sample. Owing to the fact that the model constants can be interpreted physically, the Burger’s model was selected for further analyses. Axial and transverse shale creep was analyzed in the loading range of 25–75% (σ₁ − σ₃)_max. The model parameters were identified for the shales from the Baltic Basin, and then the model was verified on data from triaxial tests.

4. Theoretical Basics of Rock Rheology in Triaxial Test Conditions for the Burger’s Model

The presented analyses refer to the behavior of shale samples in a complex stress state, therefore the description of rock behavior should be linked with relationships in a triaxial stress field. Assuming isotropy of the medium, the stress and strain tensors may be represented by their axiator and deviator. The relationships of linear elasticity take in this case the following form:

$$S_{ij}=2G e_{ij}$$

(3)

$${\sigma }_m=3K{\epsilon }_m$$

(4)

where

• σ_m is axiator of the stress tensor σ_m = (σ₁ + σ₂ + σ₃)/3 = σ_okt,
• ε_m is axiator of the strain tensor ε_m = (ε₁ + ε₂ + ε₃)/3 = ε_okt,
• S_ij is deviatoric stress tensor,
• e_ij is deviatoric strain tensor.

Referring to the model of the linear theory of elasticity, its generalization to the differential form of the constitutive equations of linear viscoelastic media in a spatial state of stress has the following form [38,41]:

$$P_1\left\{S_{ij}\left(t\right)\right\}=2Q_1\left\{e_{ij}\left(t\right)\right\}$$

(5)

$$P_2\left\{{\sigma }_m\left(t\right)\right\}=Q_2\left\{e_m\left(t\right)\right\}$$

(6)

P₁, Q₁, P₂, Q₂ are linear differential operators, which are different for the law of shear change and different for the law of volume change. In an extended notation, the equations attain the following form:

$$\left(p_0^{\prime }+p_1^{\prime } \frac{𝜕}{𝜕t}+p_2^{\prime } \frac{{𝜕}^2}{𝜕t^2}+\dots +p_a^{\prime } \frac{{𝜕}^a}{𝜕t^a}\right)S_{ij}\left(t\right)=2\left(q_0^{\prime }+q_1^{\prime } \frac{𝜕}{𝜕t}+q_2^{\prime } \frac{{𝜕}^2}{𝜕t^2}+\dots +p_b^{\prime } \frac{{𝜕}^b}{𝜕t^b}\right)e_{ij}\left(t\right)$$

(7)

$$\left(p_0^{″}+p_1^{″} \frac{𝜕}{𝜕t}+p_2^{″} \frac{{𝜕}^2}{𝜕t^2}+\dots +p_a^{″} \frac{{𝜕}^a}{𝜕t^a}\right){\sigma }_m\left(t\right)=\left(q_0^{″}+q_1^{″} \frac{𝜕}{𝜕t}+q_2^{″} \frac{{𝜕}^2}{𝜕t^2}+\dots +p_b^{\prime \prime } \frac{{𝜕}^b}{𝜕t^b}\right)e_m\left(t\right)$$

(8)

They are linear differential equations, with constant coefficients pi and qi.

Using the Laplace transformation, the equations can be noted with sub-division into axiator and deviator changes in the following form [41]:

$$S_{ij}+\left(\frac{{\eta }_s^M}{G^M}+\frac{{\eta }_s^M+{\eta }_s^K}{G^K}\right){\dot{S}}_{ij}+\frac{{\eta }_s^M{\eta }_s^K}{G^M{G}^K}{\ddot{S}}_{ij}=2{\eta }_s^M{\dot{e}}_{ij}+\frac{{2\eta }_s^M{\eta }_s^K}{G^K}{\ddot{e}}_{ij}$$

(9)

$${\sigma }_m+\left(\frac{{\eta }_m^M}{K^M}+\frac{{\eta }_m^M+{\eta }_m^K}{K^K}\right){\dot{\sigma }}_m+\frac{{\eta }_m^M{\eta }_m^K}{K^M{K}^K}{\ddot{\sigma }}_m={\eta }_m^M{\dot{e}}_{vol}+\frac{{\eta }_m^M{\eta }_m^K}{K^K}{\ddot{e}}_{vol}$$

(10)

where

• $\dot{\sigma }$ and $\ddot{\sigma }$ are the first and second derivative of the average stress, respectively,
• $\dot{e}$_vol and $\ddot{e}$_vol are the first and second derivatives of the volumetric strain, respectively,
• $\dot{S}$_ij and $\ddot{S}$_ij are the first and second derivatives of deviatoric stress tensor, respectively,
• $\dot{e}$_ij and $\ddot{e}$_ij are the first and second derivatives of the deviatoric strain tensor, respectively.

The solution of both equations, after attaining a generalized notation follows:

$${\epsilon }_{ij}\left(t\right)=\frac{1}{3}{e}_{vol}\left(t\right){\delta }_{ij}+e_{ij}\left(t\right)$$

(11)

(δ_ij is the Kronecker’s delta) and can be presented in the following form [41]:

$${\epsilon }_{ij}\left(t\right)=\left(\frac{{\sigma }_m}{3K^M}+\frac{{\sigma }_m}{{3\eta }_m^M}t+\frac{{\sigma }_m}{3K^K}\left(1-e^{-\frac{K^K}{{\eta }_m^K}t}\right)\right){\delta }_{ij}+\frac{S_{ij}}{2G^M}+\frac{S_{ij}}{{2\eta }_s^M}t+\frac{S_{ij}}{2G^K}\left(1-e^{-\frac{G^K}{{\eta }_m^K}t}\right)$$

(12)

Due to the fact that in the process of sample loading in a conventional triaxial compression test, both the axiator and deviator of the stress and strain tensors are subjected to change (Figure 1c), determination of the material constants in Equation (12) using the test results is relatively difficult. Optimization techniques can be used in this case [41] or the parameters can be estimated separately for both equations describing the shear and volumetric changes after sub-division of the deviator and axiator part of the loading path [42]. The second possibility was used in this paper. Such an approach should also allow for a separate analysis of the volumetric and shear changes during the creep process.

It should be emphasized that using linear elasticity relationships (Equations (3) and (4)) and summing up the squared both sides of Equation (1), one obtains the following [43]:

$$S_{ij}{S}_{ij}=4G^2{e}_{ij}{e}_{ij}$$

(13)

Substituting respectively [43]:

$$S_{ij}{S}_{ij}=3{\tau }_{okt}^2\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{e}_{ij}{e}_{ij}=\frac{3}{4}{\gamma }_{okt}^2$$

(14)

One finally obtains:

$${\tau }_{okt}=G {\gamma }_{okt}$$

(15)

and

$${\sigma }_{okt}=3K{\epsilon }_{okt}$$

(16)

with

$${\tau }_{okt}=\frac{1}{3}\sqrt{{\left({\sigma }_1-{\sigma }_2\right)}^2+{\left({\sigma }_2-{\sigma }_3\right)}^2+{\left({\sigma }_3-{\sigma }_1\right)}^2}$$

(17)

$${\gamma }_{okt}=\frac{2}{3}\sqrt{{\left({\epsilon }_1-{\epsilon }_2\right)}^2+{\left({\epsilon }_2-{\epsilon }_3\right)}^2+{\left({\epsilon }_3-{\epsilon }_1\right)}^2}$$

(18)

Both equations describe the elastic relationships between the octahedral tangent and normal stress, and the corresponding octahedral strains. In a similar procedure as for the axiator and deviator of the stress and strain tensors, the linear viscoelastic relationships, which link octahedral stresses and strains can be obtained as below (Figure 5):

Figure 5. Burger’s model: (a) shear part, (b) volumetric part.

Figure 5. Burger’s model: (a) shear part, (b) volumetric part.

$${\gamma }_{okt}=\frac{{\tau }_{okt}}{G^M}+\frac{{\tau }_{okt}}{G^K}\left(1-exp\left(-\frac{G_K}{{\eta }_s^K}\right)t\right)+\frac{{\tau }_{okt}}{{\eta }_s^M}t$$

(19)

$${\epsilon }_{okt}=\frac{1}{3}\left[\frac{{\sigma }_{okt}}{K^M}+\frac{{\sigma }_{okt}}{K^K}\left(1-exp\left(-\frac{G_K}{{\eta }_m^K}\right)t\right)+\frac{{\sigma }_{okt}}{{\eta }_m^M}t\right]$$

(20)

Noting equation (12) in the form of separate volumetric and shear parts, and the relationships between the octahedral stresses and strains (19), (20), the material constants in the Burger’s model can be relatively easily determined based on the estimation of non-linear parameters of the equation, using for example Statistica software (Statistica v 13.3, Tibco Software Inc., Palo Alto, CA, USA) [44].

4.1. Parameters of Burger’s Model for the Paleozoic Shales from the Baltic Basin

The estimated parameters of the Burger’s model are presented in Table 1, Figure 6 and Figure 7. As commonly known, rock creep is sensitive to many factors, including stress state, temperature, or humidity [35,39]. The influence of the level of sample loading (σ₁ − σ₃)/(σ₁ − σ₃)_max, i.e., the stress state on creep in the case of the analyzed shales is clearly visible (Table 1, Figure 6).

The behavior of the analyzed shales is variable depending on the deviator or axiator changes. In the case of deviator loading, its increase causes a decrease of all parameters in Burger’s model, which means that the samples display lower stiffness and larger susceptibility to creep. Similar conclusions can be drawn from the analysis of viscosity coefficients. At higher values of deviator stress, creep velocity described as $\frac{{\tau }_{0kt}}{{\eta }^M}$ is higher than at lower levels of loading.

A different situation can be observed in the case of volumetric change (Figure 7). Sample stiffness increases with increased loading (σ₁ − σ₃)/(σ₁ − σ₃)_max, and the analyzed rock has a smaller tendency to creep. Creep velocity decreases with loading increase according to the relationship $\frac{{\sigma }_{0kt}}{{\eta }^M}$. Unfortunately, large dispersal of the results in the case of volumetric changes hampers their quantitative analysis.

The samples from particular lithostratigraphic units (Pelplin, Pasłęk, Jantar, and Sasino) are marked with different colors in Figure 6 and Figure 7. It is evident that creep in shales is different among the particular units. Results obtained for the Jantar Formation differ significantly from the samples from the Pelplin and Pasłęk Formations, which in turn behave similarly. Results for the Sasino Member are strongly dispersed and difficult to assign to any of the mentioned clusters.

The results are influenced by the mineral composition of the analyzed samples. As mentioned in Section 2, samples from the Pelplin and Pasłęk Formations have a similar mineral composition, while samples from the Jantar Formation and Sasino Member are characterized by a variable mineral composition, strongly influencing the observed results.

The determined shear G^M and bulk K^M moduli allow Young’s longitudinal elasticity modulus and Poisson’s coefficient to be defined. Results obtained by other methods are conformable with the results presented in Wilczyński et al. [34] in earlier analyses of the Paleozoic shales from the Baltic Basin.

In the drawings Figure 6 and Figure 7 one can notice a large scatter of the presented results. Due to the digital measurement of displacement and force (MTS system), the errors associated with the measurements themselves in experiments are relatively small. However, the greatest variability of the results is caused by the natural heterogeneity of the research material, because the shale samples came from the drill core. No preliminary selection of research material was carried out. This is clearly visible in Figure 6 and Figure 7, where a large spread of the obtained model parameters for individual formations is observed.

For each creep test, the coefficient of determination R² was calculated (Table 1), determining the quality of estimation of the parameters of the Burger’s model equations and the correctness of matching the model to the experimental test results. The obtained R² values indicate a good fit of the model to the creep results.

For comparison of shale creep from the Baltic Basin with literature data as on Barnett, Haynesville, and Eagleford shale [27], as well as on Posidonia shale [24] we used a power law model proposed by Sone and Zoback [27,45]:

$${\epsilon }_{creep}=B{t}^n$$

(21)

where B and n are constant.

The analyses were performed for the load level of samples corresponding to approximately (σ₁ − σ₃)/(σ₁ − σ₃)_max = 0.5, similar to the data of comparative literature tests. The results are given in

Table 2. Constant of the power law model for shales in the Baltic Basin.

Sample	Sample	(σ1 − σ3)	B	n	R2
Location	Number	${\overline{\mathsf{(}{\mathsf{\sigma }}_1-{\mathsf{\sigma }}_3\mathsf{)}}}_{\mathit{m}\mathit{a}\mathit{x}}$	[MPa−1] × 10−5	[-]
Pelplin	B-211-A	0.65	1.6	0.043	0.93
	M-30-B	0.52	2.0	0.044	0.94
	W-260-A	0.46	1.6	0.037	0.97
	W-260-B	0.52	1.7	0.045	0.97
	W-264-A	0.55	1.8	0.044	0.93
Average			1.7	0.043
Standard dev.			0.11	4.4 × 10−5
Pasłęk	B-236-A	0.41	1.7	0.023	0.94
	W-277-B	0.46	1.8	0.029	0.98
	W-278-A	0.56	1.9	0.032	0.98
	W-278-B	0.69	1.6	0.048	0.98
Average			1.8	0.033
Standard dev.			0.06	3.4 × 10−4
Jantar	B-275-B	0.52	2.6	0.040	0.99
	B-279-D	0.46	1.8	0.039	0.99
	M-43B-B	0.46	2.0	0.036	0.97
	W-305-A	0.49	2.0	0.048	0.76
	W-306-A	0.52	2.0	0.046	0.93
Average			2.1	0.042
Standard dev.			0.37	1.1 × 10−4
Sasino	M-60-B	0.60	1.9	0.036	0.96
	M-62-A	0.48	1.8	0.020	0.91
	W-320-B	0.42	2.4	0.027	0.82
Average			2.0	0.028
Standard dev.			0.17	1.3 × 10−4

. The literature data that were accepted for comparison, only horizontally oriented samples (bedding parallel) compared to those tested in this article, are presented in

Table 3. Parameters of the power-law creep model determined for shale rocks used for comparative analyses.

Shale	Ref.	B	n	QFP	Cb	Clay	σTCS	pc	Temp.
		[MPa−1] × 10−5	[-]	vol%	vol%	vol%	[MPa]	[MPa]	[oC]
Barnett-1 (Ba1)	a	2.0–2.6	0.012–0.021	48	2	50	210	<30	room
Barnett-2 (Ba2)	a	1.6–1.6	0.009–0.010	42	48	10	325	<30	room
Hayneville-1 (Ha1)	a	1.8–2.7	0.027–0.062	32	20	48	145	<30	room
Hayneville-2 (Ha2)	a	1.5–1.8	0.011–0.049	23	49	28	240	<30	room
Eagle Ford-1 (Ea1)	a	1.7–2.3	0.024–0.053	24	46	30	200	<30	room
Eagle Ford-2 (Ea2)	a	1.7–1.8	0.023–0.049	14	66	20	175	<30	room
Pos_Dot	b	2	0.05	14	42	43	175	20	20
Pelplin	c	1.7	0.043	40	14	46	215	50	85
Paslek	c	1.8	0.033	36	7	57	155	50	85
Jantar	c	2.1	0.042	31	22	46	215	50	85
Sasino	c	2	0.028	51	3	45	222	50	85

Ref.—reference, a—Sone and Zoback [27], b—Rybacki et al., [24], c—this study., B, n—constants, QFP—quartz + feldspar + pyrite, Cb—organic carbon, σ_TCS—triaxial-compressive strength, p_c—confining pressure, temp.—temperature.

. From the Baltic Basin, the average values of parameters B and n were taken for comparison (Table 3, Figure 8). In Figure 8, the power law model parameter values determined for American shales [27] are shown in the appropriate ranges consistent with the data in Table 3.

Considering that B essentially reflects the instantaneous elastic property of the rock [27], and that n determines the relative contribution of the time-dependent deformation to the total strain, a similar creep behavior of the compared shales can be observed. It is worth emphasizing that the testing conditions for the shales presented in our article were not exactly the same due to the different p_c, and T conditions. When comparing the results, differences resulting from the mineralogical composition of individual formations should also be taken into account.

5. Verification of Burger’s Model for Creep Description in Paleozoic Shales in Conventional Triaxial Loading Conditions

Using relationship (12), equations describing axial ε₁ (22) and lateral ε₂ (23) strain can be determined in the creep test in triaxial stress conditions (σ₂ = σ₃ = σ). The strains are the sum of elastic strain and viscous strain related with creep. The results of these analyses, verifying the correctness of assumptions in Burger’s model to describe creep in the analyzed shales are presented in Figure 9.

$${\epsilon }_1\left(t\right)=\frac{{\sigma }_1+2\sigma }{9K^M}+\frac{{\sigma }_1+2\sigma }{9{\eta }_m^M}·t+\frac{{\sigma }_1+2\sigma }{9K^K}\left(1-e^{-\frac{K^K}{{\eta }_m^K}·t}\right)+\frac{{\sigma }_1-\sigma }{3G^M}+\frac{{\sigma }_1-\sigma }{3{\eta }_s^M}·t+\frac{{\sigma }_1-\sigma }{3G^K}\left(1-e^{-\frac{G^K}{{\eta }_s^K}·t}\right)$$

(22)

$${\epsilon }_2\left(t\right)=\frac{{\sigma }_1+2\sigma }{9K^M}+\frac{{\sigma }_1+2\sigma }{9{\eta }_m^M}·t+\frac{{\sigma }_1+2\sigma }{9K^K}\left(1-e^{-\frac{K^K}{{\eta }_m^K}·t}\right)+\frac{\sigma -{\sigma }_1}{6G^M}+\frac{\sigma -{\sigma }_1}{6{\eta }_s^M}·t+\frac{\sigma -{\sigma }_1}{6G^K}\left(1-e^{-\frac{G^K}{{\eta }_s^K}·t}\right)$$

(23)

The experiments show that Burger’s model correctly describes creep in the representative sample W260A of Paleozoic shale in triaxial tests. In the presented dependencies, creep is described with regard to shear and volumetric strain. In this case, the contribution of particular creep components in the value of the entire strain is an interesting question.

6. Assessing the Significance of Volumetric Creep Description in Paleozoic Shales

In creep analysis it is often assumed that volumetric changes within transient creep are small and may be treated as elastic, while effects of creep are only related with shear change. In this case, Equation (22) attains the following form:

$${\epsilon }_1\left(t\right)=\frac{{\sigma }_1+2\sigma }{9K}+\frac{{\sigma }_1-\sigma }{3G^M}+\frac{{\sigma }_1-\sigma }{3{\eta }_s^M}·t+\frac{{\sigma }_1-\sigma }{3G^K}\left(1-e^{-\frac{G^K}{{\eta }_s^K}·t}\right)$$

(24)

Such assumptions in the creep calculations were accepted by some researchers [36,46] and discarded by others [41,47]. Including volumetric change in the creep description largely depends on the rock type and the mechanism of creep in that rock. The contribution of volumetric change in the creep of Paleozoic shales can be traced on the creep chart for sample W260A, where shear and volumetric changes are presented with regard to time (Figure 10).

For low levels of loading (σ₁ − σ₃)/(σ₁ − σ₃)_max, the values of volumetric and shear strain are similar and of the same order of magnitude; however for higher loading values, shear strain dominates, while volumetric strain remains at a similar, low level. These phenomena are related to changes of elasticity and viscosity parameters presented in Section 4.1, where rock deviatoric stiffness decreased with increased loading, and volumetric stiffness increased. The contribution of volumetric creep in Burger’s model is expressed by Equations (22) and (24), as well as Figure 10b. Worth noting is that volumetric strain has a very low influence on overall strain in the analyzed Paleozoic shales. To conclude, in the case of Paleozoic shales from the Baltic Basin the part related with volumetric strain may be ignored in creep description and this strain can be treated as elastic.

7. Linear Creep Limit in Paleozoic Shales

The determination of the linear creep limit is significant with regard to describing the loading range, in which the Boltzmann superposition is valid and linear viscoelastic relationships can be used.

The linear creep limit can be determined by analyzing creep velocity at various loading levels [39]. Creep velocity which is the tangent of the angle of inclination of the creep curve is constant in the secondary creep range (Figure 3 and Figure 11). Velocity increases proportionally (linearly) depending on the loading level. The linear limit is determined by the value of loading, for which the relationship of creep velocity loses its linear character.

Values of creep velocity at different loading levels for shales from different areas are compiled in Figure 12.

Based on the performed analyses (Figure 12) it can be assumed that the linear creep limit for the analyzed Paleozoic shales is 0.7 (σ₁ − σ₃)_max. It should be assumed that beyond this loading, creep in shale will be non-linear, and its description will exceed the Boltzmann superposition rule. Such creep can be described by the plasticity theory [48] or by failure mechanics [49] depending on the type of failure. However, in this case various creep tests should be made at loadings close to instantaneous strengths in order to recognize the failure mechanism accompanying creep.

8. Summary and Remarks on Creep in Shales with Regard to Loadings Exceeding Secondary Creep

The article presents the results of creep tests performed on Paleozoic shales in conventional triaxial stress conditions. The presented results and analyses were focused only on linear creep in the range of secondary creep, for the description of which viscoelastic equations can be applied. The character of creep in the accepted loading range shows that Burger’s model can be used for its description. An original approach together with analysis results are presented here, which allow the separation and monitoring of shear and volume creep effects, and on this basis, to determine the significance of the contribution of volume creep in the entire creep process. A relatively simple methodology to determine the parameters of the Burger’s model using this division is presented. The original value of the article is also the test results themselves and the parameter values of the analyzed model for triaxial creep of shales, which are not numerous in the literature. The obtained results show that in the case of the analyzed shales, creep related with volumetric change, which may be treated as elastic, can be omitted from creep description. This conclusion allows for further comments and assumptions related to the creep mechanism. The creep mechanism, in which the contribution of volumetric strain is low, is probably caused by a small share of crack propagation according to the growth of mode I subcritical cracks [2]. It may correspond on a micro-level to phenomena such as frictional slip, pressure solution, internal grain deformation, or grain sliding [18]. This indicates a ductile creep mechanism. In the studied shales, this behavior may be related to the high content of clay minerals. In the case of brittle or porous rocks, volumetric change during creep plays a significant role, and the creep mechanism is significantly related to development of mode I subcritical cracks [2,50,51,52], pore degradation, and stiffness reduction described by the damage mechanism. In this case the description of tertiary creep (Figure 3) may also be correlated with development of failure, and a scalar or tensor damage variable. Creep, whose mechanism corresponds to diffusion, may be described in the third range by plasticity, where the level of volumetric change may be regulated by the dilatation angle.