1 Introduction

Hydraulic fracturing has been the critical procedure to exploit unconventional gas reservoirs, such as shale gas reservoirs (Saldungaray and Palisch 2013). However this technique has several major shortcomings (Cha et al. 2014). Firstly, the water used in this practice can cause clay swelling and thus the capillary fluid retention significantly reduces the permeability. Secondly, the required large amount of water consumption may be detrimental on the local water supply. Finally, the down hole injection of chemicals commonly required in hydraulic fracturing can contaminate the local underground water resource. On the contrary, cryogenic fracturing, such as liquid nitrogen cooling, is more promising due to its efficiency and it can avoid the defects of hydraulic fracturing.

When liquid nitrogen is injected, the heat from the rock quickly transfers into the liquid nitrogen, which causes rapid cooling of the rock. This rapid cooling can result in the contracting of rock surface and even fracturing when tension is sufficiently built up. In tension failure zones, the fractures change the permeability.

Liquid nitrogen stimulation has been increasing highlighted in recent year. Considerable amount of experiments, field tests and numerical studies have been performed. Finnie et al. (1979) presented a calculation model for a semi-infinite slab cooled on its surface using liquid nitrogen and conducted experiment modelling limestone injected with liquid nitrogen. Grundmann et al. (1998) treated a Devonian shale well with cryogenic nitrogen and observed that the initial production rate increases by 8 % compared with that of a nearby offset well using traditional hydraulic fracturing. Liquid nitrogen has also been used for fracturing reconstruction of reservoirs with favorable stimulation effects (Mcdaniel et al. 1997). Ren et al. (2013) studied the effects of liquid nitrogen on coal-rock samples in thermal shock experiments and ultrasonic tests. The experimental results show that micro-fractures can be induced inside the coal-rock matrix due to liquid nitrogen shock. The dramatic reduction of sound velocity and amplitude implies that the liquid nitrogen treatment effectively affects the internal structure of the coal-rock samples. Cai et al. (2014a, b) investigated the effects of liquid nitrogen cooling on rock pore structure. The results show that there are three major changes in the rock pore structure when the samples are cooled by liquid nitrogen: (1) a reduction in the number and volume of the pores; (2) an expansion of the micro-fissures (micro-pores); and (3) an increase in the pore scale.

As a promising stimulation method, the mechanism of liquid nitrogen fracturing, however, has not been fully understood. This is especially the case of the permeability evolution of rock under liquid stimulation, which is a key issue to understand the mechanism of liquid stimulation for engineering design. Interestingly, the permeability evolution has not been extensively studied in the literature. Zhang et al. (2015a) proposed a permeability evolution model subjected to liquid nitrogen cooling. However the permeability evolution in this model is dependent of the element size. This paper numerically models the cooling of a rock sample under liquid nitrogen. The enhanced permeability of the volumetric strain was defined, and the evolution of permeability is studied by relating to the volumetric change. The commercial package FLAC (ITASCA 2012) with FISH function is employed in this study. A semi-submerged test of liquid nitrogen was performed. The simulating results are proved to be right by the test results. This simulation results are discussed with the hope to provide some insight into understanding the nitrogen cooling practice.

2 Permeability evolution model

In this paper, the permeability is directly related to the volumetric change of pores and fractures of rock, which is in turn influenced by the thermal deformation, cooling fracturing by liquid nitrogen shock.

When water becomes ice, the volume increases by 9 %, which is much larger than any kind of thermally-induced expansion or shrinkage in geomaterials. However unsaturated rock-coal is subsistent and the effects of thermal shock is focused. So the freezing damage is ignored in this paper.

2.1 Temperature

By assuming that strain changes play a negligible role in influencing the temperature, the differential expression of the energy balance of the rock has the following form:

$$KKT_{,ij} + q_{V} = \rho C_{V} \frac{\partial T}{\partial t}$$
(1)

where KK is the thermal conductivity tensor, q V is the volumetric heat source intensity, C V is the specific heat at constant volume, \(\rho\) is the mass density of rock, T is temperature.

2.2 Mechanical coupling

The thermal strain and stress of rock due to temperature changes can be expressed as (Zhang et al. 2015a):

$$\Delta \varepsilon_{ij} = \alpha_{t} \Delta T\delta_{ij}$$
(2)
$$\Delta \sigma_{ij} = \alpha_{t} E\Delta T\delta_{ij}$$
(3)

where α t is the coefficient of linear thermal expansion, ∆T is the temperature change, δ ij is the Kronecher delta, E is Young’s modulus.

2.3 Elemental constitutive law

The constitutive relationship of the rock material can be simplified as Fig. 1, which includes three stages: linear elastic (OA or OD), perfect brittle (AB or DE) and plastic stage (BC or EF) under compression or tension, respectively (Zhang et al. 2015a).

Fig. 1
figure 1

Simplified deformation process for rock element (Zhang et al. 2015a)

Full size image

In linear elastic stage the components of stress are linear functions of the components of strain:

$$\Delta \sigma_{ij} = 2G\Delta \varepsilon_{ij} + \left( {K - 2/3} \right)G\Delta \varepsilon_{kk} \delta_{ij}$$
(4)

where K is the bulk modulus and G is the shear modulus. When the peak stress is researched, i.e. shear and tension failure begin, the perfect brittle degradation occurs. In this case, the non-linear Heok–Brown failure criterion (Hoek and Brown 1997) is appropriate to describe the non-linear behavior considering the confining pressure:

$$F = (\sigma_{1} - \sigma_{3} ) - \sqrt {m_{i} \sigma_{c} \sigma_{3} + s\sigma_{c}^{2} } = 0$$
(5)

where σ 1 and σ 3 are the major and minor effective principal stresses at failure, σ c is the uniaxial compressive strength of the intact rock material, m i and s are constants for the rock mass.

The tension failure criterion can be expressed as:

$$F = \sigma_{3} - \sigma_{t} = 0$$
(6)

where σ t is the tension strength of rock element.

The extension of inner micro fissures degrades the rock mechanical behavior. In case of tension failure, the residual tension strength, σ tr is maintained. For shear failure, the residual strength and modulus are dependent on confining pressure (Fang and Harrison 2001; Zhang et al. 2015a):

$$\sigma_{r} = \sigma_{11} \left( {1 - \frac{{r\delta \sigma_{u} }}{{\sigma_{c} }}} \right)$$
(7)
$$E_{r} = E_{0} - r\delta E_{0}$$
(8)

where δσ u is the strength degradation of the uniaxial case, σ is peak strength of the confining pressure of σ 3, E 0 is initial Young’s modulus, E r is the Young’s modulus during residual stage, δE 0 is the modulus degradation of the uniaxial case, r is the strength degradation index, and is expressed as (Fang and Harrison 2001; Zhang et al. 2015a) :

$$r = \exp ( - n_{d} \sigma_{3} )$$
(9)

where n d is a fitting coefficient.

The residual strength of rock element also follows Hoek–Brown failure criterion with the degradation of s and m i , which are written as (Fang and Harrison 2001; Zhang et al. 2015a):

$$s_{d} = \left( {\frac{{\sigma_{cd} }}{{\sigma_{c} }}} \right)^{2}$$
(10)
$$m_{d} = m_{i} \left( {\frac{{\sigma_{cd} }}{{\sigma_{c} }}} \right)^{9/14}$$
(11)
$$\sigma_{cd} = \sigma_{c} - r\delta \sigma_{u}$$
(12)

where s d and m d are the Heok–Brown material parameters for the residual strength.

The Mohr–Column failure criterion can also be employed to describe the shear failure:

$$F = \sigma_{1} - k\sigma_{3} - \sigma_{c} = 0$$
(13)

where k = (1 + sin φ)/(1−sin φ).

According to the relationship of material parameters between Heok–Brown failure criterion and Mohr–Column failure criterion, the equivalent inner friction angle and cohesion of Mohr–Column failure criterion can be expressed with m and s of Heok-Brown failure criterion (Fang and Harrison 2002; Zhang et al. 2015a):

$$\phi_{\text{h}} = 2{\text{arctan(}}\sqrt {k_{\text{h}} } - \pi/2)$$
(14)
$$c_{\text{h}} = \frac{{m\sigma_{3} + 2s\sigma_{c} }}{{4\sqrt {k_{\text{h}} } \sqrt {m\sigma_{3} /\sigma_{c} + s} }}$$
(15)

where subscript “h” represents expressions of Heok-Brown parameters, m and s are Hoek–Brown material parameters. A parameter k h equivalent to k in Eq. (13) can be obtained as:

$$k_{\text{h}} = 1+ \frac{m}{{2\sqrt {m\sigma_{3} /\sigma_{c} + s} }}$$
(16)

2.4 Dilatancy

In this paper, the dilatancy is considered to occur at the peak strength of the rock. This is a reasonable simplification as the dilatancy is minor when the stress is low (Besuelle et al. 2000; Arzúa and Alejano 2013).

Dilation angle usually is employed to predict the volumetric expansion. Dilation angle, which is closely related to confining pressure, can be described as (Yuan and Harrison 2004, 2005; Zhang et al. 2015b):

$$\psi_{p} = \arcsin \frac{{\tan (I_{d} \theta_{0} )}}{{2 + \tan (I_{d} \theta_{0} )}}$$
(17)
$$\theta_{0} = \arctan \left( {\frac{{2\sin \psi_{0} }}{{1 - \sin \psi_{0} }}} \right)$$
(18)

where ψ 0 is the dilation angle of unaxial compression state, I d is the dilatancy index, and is experimentally written as :

$$I_{d} = \left( {1 - \frac{{\sigma_{3} }}{{\sigma_{h} }}} \right)^{{k_{d} }}$$
(19)

where σ h is the lowest confining stress at which no dilatancy takes place. The parameter k d controls the curvature of the relations.

2.5 Evolution model of permeability

Permeability, describes how easily fluid can pass through a porous medium. This flow is related to the connected pore volume of a rock matrix and the hydraulic aperture of a fracture. As both of these attributes may be characterized in terms of either length units (e.g. aperture), or ratios of length units (e.g. connected pore volume as a proportion of total volume), it is reasonable to relate permeability to strain. Consequently, permeability of rock is considered as strain-dependent property in this paper.

When an element of rock undergoes dilatancy with degradation of both strength and stiffness due to changes of volume, its hydraulic properties change. Experimental results (Yu et al. 2013; Badrul Alam et al. 2014; Tan et al. 2014) indicate the relationship between the volumetric strain and the permeability takes on positive correlation. The enhanced permeability of the volumetric strain is defined as:

$$\lambda = \frac{{{\text{d}}\eta }}{{{\text{d}}\varepsilon_{V} }}$$
(20)

where η is the permeability, ε V is the volumetric strain.

The permeability may be expressed as:

$$\eta = \eta_{0} + \frac{{{\text{d}}\eta }}{{{\text{d}}\varepsilon_{V} }}\varepsilon_{V}$$
(21)

where η 0 is the permeability without the volumetric strain.

λ can be obtained by the tri-axial compression and seepage test(Zhang et al. 2015b). The common form is expressed as:

$$\lambda = a_{0} + a_{1} \varepsilon_{V}^{2} + a_{2} \varepsilon_{V}^{4} + \ldots$$
(22)

where, a 0, a 1, a 2… are fitting constants.

The permeability of some coal and sandstone has been tested, and the parameter of λ has been fitted by the test data (Zhang et al. 2015b).

2.6 Implementation

The above proposed evolution model of permeability in this paper could be implemented in either continuum or discontinuous numerical methods. The former is able to deal with fracture genesis, whereas the latter is primarily intended to handle a pre-existing discontinuous system. As the proposed model intends to simulate fracture and permeability evolution subject to liquid nitrogen cooling, it is appropriate to use a continuum method.

In this paper the explicit finite difference code FLAC (ITASCA 2012) has been selected for the implementation. In FLAC the mechanical deformation, temperature change and fluid flow can be coupled. FISH function approach is employed and the model in this paper is implemented as shown in the flow chart in Fig. 2.

Fig. 2
figure 2

Procedures for implementing the evolution of permeability

Full size image

3 Case study

An analysis example is presented using the implemented model in FLAC to study the effects of liquid nitrogen cooling on the permeability of a methane recovery well in Xinhua mine, Jilin province, China.

The mechanical and thermodynamic parameters in this paper were obtained through the mechanical and thermodynamic experiments as listed in Table 1. The depth of the well is 230 m, and the mean unit weight of overburden is 17.5 kN/m3.

Table 1 Parameters of case study
Full size table

Calculation domain is 20 m × 20 m. The diameter of the well is 0.48 m. The initial temperature of coal is 32 °C. The boiling point of liquid nitrogen at atmospheric pressure is −195.8 °C. The ratio of calculation domain length to radius of the well is more than 40, and the effects of boundary condition on results are ignored. The relationship between the permeability and the volumetric strain has been tested, and an approximative linear form is accepted. And through the fitting method the enhanced permeability of λ is 4.2 × 10−15 m2.The sketch of calculation model can be seen in Fig. 3.

Fig. 3
figure 3

Sketch of model

Full size image

Due to symmetry, the quarter of model is adopted, and 900 elements are generated.

Tension failure area and the permeability of elements of coal after five days under liquid nitrogen cooling can be seen in Figs. 4 and 5, respectively. After five days cooling, the tension failure zone with a band of 0.61 m, is observed in coal bed. In tension failure area cracks initiate and extend and the permeability of the element significantly increases. The maximum permeability for element is 4.4 × 10−12 m2, which is 1.83 × 105 times more than that before cooling. Apart from the bore hole, the thermal tensile stress leads to the growth of permeability at a rate of 1 %–16 %, far less than that in the tension failure area. The above conclusions have been proved to be right by the test results (Cha et al. 2014).

Fig. 4
figure 4

Tension failure area after five days of cooling coal

Full size image
Fig. 5
figure 5

Permeability distribution after five days of cooling coal (maximum is 4.4 × 10−12 m2)

Full size image

The permeability’s change of coal for varied cooling time can be seen in Fig. 6. It is inferred that with increasing cooling time, the area of tension failure expands. However, the permeability of the elements in failure area shows a slight reduction with time. The reason is that the volumetric stain of the element cuts down due to stress adjustment.

Fig. 6
figure 6

Permeability for varied cooling time

Full size image

The relationship between cooling time and failure area width can be seen in Fig. 6, too. It is found from Fig. 6 that tension failure area expands quickly at the beginning of cooling and slows down with time. At last, a 1.5 m band of cooling fracturing can be formed.

In order to verify the numerical model, a 50 mm × 50 mm × 100 mm coal sample was semi-submerged for 30 min. Figure 7 illustrates the coal sample scanned by SEM in its initial state and after being cooled in liquid nitrogen. It is found that new cracks appear and the natural cracks expand and widen, which will lead to an increase of permeability. So the above experimental results qualitatively verify the feasibility of the model in this paper.

Fig. 7
figure 7

Coal sample scanned by SEM (×100) a Initial state, b after being cooled in liquid nitrogen

Full size image

4 Concluding remarks

Based on the relationship between volumetric strain and permeability, the enhanced permeability of the volumetric strain is defined. And the permeability evolution stimulated by liquid nitrogen is numerically modelled. FISH function method in FLAC is employed to perform the model. A case study was conducted to model the behavior of a rock bed injected with liquid nitrogen. Besides a simple submerged test of liquid nitrogen cooling was performed. Our main conclusions include:

  1. (1)

    In tension failure area, crack’s initiation and extension lead to the significant growth of permeability.

  2. (2)

    The enhanced permeability of the volumetric strain was defined, and the permeability was directly evaluated based on element’s volumetric strain. Further a new method to model the permeability evolution of rocks subject to liquid nitrogen shock is given out.

  3. (3)

    At last a 1.5 m band of cooling fracture was found in this study.

  4. (4)

    A semi-submerged test of liquid nitrogen was performed in this paper. The method to model the permeability evolution of rocks subject to liquid nitrogen shock in this paper is proved to be right by the test results.

  5. (5)

    This study is hoped to provide a simple but reasonable description of the permeability evolution of rocks subject to liquid nitrogen cooling.