Abstract

Non-Darcy flow is observed in the shale gas reservoir because it is rich in organic nanopores. Generally, the permeability of shale gas reservoirs is modified because of non-Darcy flow. However, the viscosity is much less concerned. It has been verified that the viscosity of dilute gas depends on the size of the pore. In this paper, the viscosity of methane in organic slit nanopore is determined with equilibrium molecular dynamics (EMD) simulation. The result shows that the viscosity of bulk methane would decrease with dropping down pressure, while the confined effect would make the viscosity of methane in the organic slit nanopore lesser than that of the bulk phase, and it decreases severely at low pressure. The confined dense gas viscosity model is obtained by theoretical derivation. The EMD results were fitted with this model to obtain the viscosity correction method for dense methane in organic slit nanopores. The dimensionless viscosity (μeff/μb) would decrease sharply with the Knudsen number between 0.1 and 10. Unlike the confined effect on the dilute gas, the potential contribution of the dense gas and the wall also affects its viscosity. Because of the confined effect on the dense methane, the flow capacity of methane is enhanced 1.5 times at least with the pore being smaller than 10 nm and the pressure being lower than 5 MPa. It means that keeping a low reservoir pressure helps to improve the flow of shale gas. This work can improve the understanding of the importance of gas viscosity with the non-Darcy flow in shale gas reservoirs.

1. Introduction

The problem of modeling the transport of fluid in confined spaces has been attracting the attention of scientists and engineers for over a century because of its importance in a variety of applications and industrial interests [1]. Especially, the topic gained attention in the last decades in the petroleum development area as the shale gas reservoir was developed successfully. The flux of gas through the pore is composed of viscous donation and diffusion donation [2]. The transport of gas becomes more complex because of the presence of adsorption [1].

Generally, Darcy’s law is used to depict the viscous flow process of fluid in the pore. As shown in equation (1), the permeability and viscosity determine the velocity of the fluid in the pores at a specific pressure gradient [3]. Permeability is the ability of a porous material to allow fluids to pass through it, which depends on the porosity, tortuosity, connection, shape, and pore sizes of the porous material. The permeability (k) of the circular pore and the slit pore can be derived from the Navier–Stokes equation, which are and , individually. Viscosity reflects a fluid’s resistance to flow under a specific differential pressure [3], which is generally affected by the pressure, temperature, and fluid species.

Most of the shale gas reservoir pores are nanopores that are in the range of 2∼50 nm, with pressures between 2∼50 MPa and temperatures between 20∼140 °C [4, 5]. Darcy’s law should be modified in shale gas development because shale gas experiences multiscale flow at shale gas reservoir conditions. Thus, many methods are proposed to modify the permeability of a shale gas reservoir [59]. However, there are much fewer concerns about the viscosity of shale gas.

Gas is classified as dilute gas and dense gas by δ [10]. δ is defined as equation (2).

When , the gas is in the dilute regime. As molecular spacing decreases, the dense gas regime develops (). For dilute gas, the density of the gas is too small, and the kinetic-potential and potential-potential contributions are zero because the potential shear stress becomes zero. At this time, the viscosity of the dilute gas equals the kinetic-kinetic contributions according to the Chapman–Enskog formula, which is independent of pressure (or density) [11]. Because of the effect of potential, the viscosity of dense gas would increase with pressure (or density) [11]. Several methods for calculating the dense gas viscosity have been proposed [1215]. Based on the developed viscosity calculation theory [12], a semiempirical formula is proposed [14], which is consistent with the experimentally measured data from 37.8 to 171.2°C and 0.1013 to 55.158 MPa. The standard deviation of calculation is ±2.69%, and the maximum deviation is approximately 8.99% [14].

The existence of the Knudsen layer would affect the effective viscosity of the fluid. The free path area adjacent to a pore wall is referred to as the Knudsen layer, and the viscosity of the Knudsen layer is less than that in the bulk phase [16]. The effect of the Knudsen layer on the viscosity of a dilute gas was identified [1723]. The measured viscosity of the dilute gas in the Knudsen layer of nanopores is reported to be smaller than that measured in the bulk phase [2023], and the effective viscosity value decreases with increasing Knudsen number [1719]. Beskok found that unless the boundary condition was modified, the Navier–Stokes equation could not depict the velocity profile of the dilute gas in the pore when the Knudsen number is bigger than 0.1. He assumed that the viscosity could be expressed as the Bosanquet-type of approximation with a Knudsen number. Based on this thought, Beskok used the Navier–Stokes equation and general slip boundary condition to fit with the direct simulation Monte Carlo (DSMC) results to get the viscosity of the dilute gas in different shape pores [17]. NEMD simulations were used to mimic the Poiseuille flow of methane in the organic or inorganic pores of shale [2426]. NEMD showed that the gas flux was improved in the nanopore. Chen et al. concluded that the decrease of viscosity and slip effect led to the improvement of gas flux [25].

Lv et al. got a new viscosity of the dilute gas in the cylindrical pore using the effective volume diffusion hydrodynamics method to fit with the results of the linearized Boltzmann equation and the BGK model of the Boltzmann equation [19]. Michalis et al. performed DSMC simulations and used the Green–Kubo relationship to calculate the viscosity of dilute nitrogen in the transition regime, the result of which showed that the value of viscosity depended on the Kn value [18]. On the basis of the Green–Kubo relationship, Fei et al. derived an expression of viscosity with the Knudsen number, which was usable when Kn was less than 0.5 [27]. According to the best knowledge of the authors, there is no such theory that derives an analytical equation to describe the viscosity of the dense fluid, not to mention the viscosity of the dense gas in nanopores [28]. Molecular dynamics (MD) simulation is a good method to study the viscosity of the fluid in the nanopores. Both equilibrium molecular dynamics (EMD) and nonequilibrium molecular dynamics (NEMD) can be utilized to calculate the viscosity of the dense fluid in nanopores [11]. The time-related correlation functions, such as the Green–Kubo relations, are usually applied to estimate the fluid (includes dilute gas and dense fluid) transport parameters, such as viscosity, diffusivity, and thermal conductivity, while the transport parameters are expressed as the integrals of time-correlated functions to the corresponding thermodynamic fluxes (stress tensor, velocity, and heat current) [29] in EMD [3032]. In the NEMD method, through a simulation of the dense gas flow in the pore, the velocity distribution and shear stress across the pore are generated, by which the viscosity of the fluid can be evaluated using Newton’s law of internal friction [3138]. The EMD and NEMD simulation methods can generate the same results [31, 39]. Zhang et al. [30] used the equilibrium molecular dynamics simulation (EMD) to study the viscosity of the methane liquid in silicate pores. The result showed that the viscosity of the methane liquid was affected by the temperature, density, and pore width. In this paper, the EMD method is chosen.

More than 90% composition of the shale gas is methane, which stays in the nanoporous shale in a dense state [40]. It is important to clarify the viscosity of the dense methane gas in the gas-bearing shale to describe the shale gas flow. Gas-bearing shale is rich in organic pores [1, 4]. Therefore, this paper is going to concentrate on the determination of the viscosity of methane at different states (dilute and dense states) in organic nanopores in shale reservoir conditions.

This paper is organized as follows: section 2 provides the calculation detail for the methane viscosity in equilibrium molecular dynamics simulations. In Section 3, the viscosity of methane in the bulk phase is calculated, and the results are compared by applying the Chapman–Enskog theory, the Lee function, and the EMD methods, firstly. After that, the effective viscosity of dense methane in the organic pore is determined using the EMD method. Next, a theory model is promoted to fit the result of the MD simulation. This method uses the experiential controlling coefficient, α, which reflects the effect of the Knudsen layer on gas viscosity. The value of the coefficient α is 2.75, which is a little bit bigger than the results of the dilute gas [10, 18, 20]. It may be caused by the interaction of gas and pore wall. The results show that the flow capacity of dense methane in the organic pore would be improved at low pressure, which indicates that reducing the reservoir pressure helps improve the flow of the shale gas. Finally, the key conclusions are presented in Section 4.

2. Molecular Dynamics (MD) Simulations

The LAMMPSTM (large-scale atomic/molecular massively parallel simulator) was used to perform the EMD simulations [41].

2.1. Force Field

The available force field models for methane include TraPPE-UA [42], TraPPE-EH [43], Optimized Potentials for Liquid Simulations (OPLS) [44], and others. These models can accurately predict the critical parameters, the relations among the pressure-volume-temperature (PVT), viscosity, and the diffusion coefficients of methane. TraPPE-UA is a union atomic model that saves considerable time when compared to other models that calculate at the same scale, and therefore, it is used in this study. Two parallel graphene sheets are used to represent the organic slit pore of the gas-bearing shale [4547]. The width of the organic pore is adjusted by the distance of two graphene sheets. 12–6 Lennard–Jones potential function (Equation (3)) is used in the intermolecular interaction, and detailed parameters are presented in Table 1. The cutoff radius in the present study is assumed to be 5σ.

The interaction potential between the methane molecules and C atoms follows the Lorentz–Berthelot rules and is calculated by equations (4) and (5).

2.2. Viscosity Calculation

Using the kinetic theory concept, viscosity can be described in relation to the transfer of the momentum between the different layers of fluid moving at different velocities at the microscopic level [32]. The Green–Kubo function, developed by Green [48] and Kubo [49], can be used to calculate the fluid transportation properties. Shear viscosity can be calculated by the integration of an autocorrelation function of the off-diagonal elements in the stress tensor , as shown in equation (6) [29].

The microscopic expression of shear stress is presented in equation (7).where is the velocity in α direction of molecule i. The shear stress of fluid in the pore can be mathematically represented as equation (8).

For the bulk phase, all the off-diagonal elements in the stress tensor are equal. For the gas confined in the slit pore because of the interaction between the wall and the fluid atoms, μxz and μyz are enhanced compared with μxz. The viscosity of the gas confined in the slit pore can be written as equation (9) [32].

2.3. The Width of Slit Pore

Figure 1 shows the typical solid-fluid potential in a slit pore with the physical width L. Because of the repulsion between the gas and the wall, the gas molecule can only access the region with width L′. L′ is defined as the accessible pore width, as shown in equation (10), [50], which corresponds to the pore size from the experiment [51]. Therefore, L′ is going to be used to calculate the Knudsen number of simulation data.

L is the physical width of the pore, which is the distance between the center of the atoms at the surface of the two walls, Å; σf is the effective collision radius of the fluid, and Å; z0 is the distance at which the solid-fluid potential is zero. The potential distribution of methane in the slit pore is counted by moving a methane molecule in the pore, and the z0 is 3.03 Å for the pore with 10 Å width and 3.05 Å for pore with a width between 25∼400 Å.

2.4. Simulation Details

The workflow to calculate the viscosity of confined methane by molecular simulation is shown in Figure 2.

The viscosity of methane at the temperature of 343 K, the pressure from 1 MPa to 16 MPa, and the in-pore width from 1 nm to 40 nm was simulated. The length and width of the slit pore were 4.26∼42.6 nm and 4.92∼49.2 nm and were adjusted so that a sufficient number (>600) of methane molecules was contained in the box. Figure 3 shows a snapshot of the simulation to illustrate the geometry of the system. For the simulation of the methane molecules inside the channel, the confinement was only in the z-direction, perpendicular to the channel walls, and the periodic boundary conditions were applied in the x- and y-direction.

Firstly, the grand canonical ensemble Monte Carlo (GCMC) method is used to get a configuration of methane in a specific pore under a defined temperature and pressure. The GCMC method simulates gas molecules in the pore exchange with an imaginary ideal gas reservoir at the specified temperature and chemical potential (µ) [52]. It also attempts Monte Carlo (MC) moves (translations and molecule rotations) within the simulation cell or region. Eventually, the potential and temperature of the gas in the pores are equal to the imaginary ideal gas reservoir. The chemical potential is related to pressure (p) and fugacity coefficient (ϕ) as shown in Appendix A. The fugacity coefficient is obtained from the RK EOS in Appendix B. 0.1 to 1 million Monte Carlo move steps are run to get equilibrium, and then 1 million Monte Carlo move steps are run to the product equilibrium configuration. Secondly, the average configuration of GCMC in the production process is set as the initial configuration in the NVT ensemble molecular dynamic simulation. At the NVT ensemble, the number of molecular and system volume, the number of molecules, and the volume of the simulation box remain constant. The system temperature is controlled at a defined value by the Nose–Hoover method [32]. The equation of motion is integrated using velocity-Verlet scheme in every time step until the system achieves a steady-state with time running. In this study, the system is run to reach equilibrium in the last 4 ns, with each time step of 1 fs. Finally, the stress correlation function data is collected after equilibrium is acquired. The Green–Kubo equation is used to calculate the viscosity. The dump time (Td), which is the integration upper bond of the Green–Kubo equation, is a vital parameter to calculate viscosity. Theoretically, it is infinity. As emphasized by Stadler et al., the statistical error increases very rapidly [53]. Some authors proposed to choose dump time as the time for which the TAF first crosses the zero value [54]. In this paper, the dump time is correctly and suitably chosen by trying to assume the continuity of the functions and of their derivatives. Using this criterion, the dump time was chosen, which ranged from 3000 to 16,000 fs for different cases in this study. The sampling interval is 10 steps. The collection time is 400 times the dump time. To obtain accurate viscosity data, each data point was the average of 10 independent simulation outcome results [55].

3. Results and Analysis

3.1. Bulk Gas Viscosity

Bulk methane viscosity at different pressures and 343 K was calculated using the Chapman–Enskog formula in Appendix C, the Lee equation in Appendix D, and the molecular simulation method. As Figure 4 shows, the methane viscosity varied slightly at low-pressure ranges but increased in high-pressure ranges. The results of the MD method matched well with the experimental data reported by NIST values [56], demonstrating that the MD simulation results were valid and reliable at shale gas reservoir conditions.

3.2. Gas Viscosity in Slit Nanopore

The effective gas viscosity for different pore widths at 1, 2, 4, 8, 12, and 16 MPa was simulated using the MD method described in Section 2, and the results are plotted in Figure 5.

Confine spaces may decrease the effective viscosity of methane. In the pore, the effective viscosity of gas decreased with the decreasing pore width, while in the smaller pores, the pore width has more effect on effective gas viscosity. Figure 6 reflects that the confined effect would make the effective viscosity of methane decrease at least 10% for pore size less than 40 nm when the pressure is 4 MPa. It means that the confined effect on viscosity should be paid more attention to. The viscosity change rate increases as pressure decreases. For small pores, the trend is even more pronounced.

There are two differences in effective viscosity between bulk gas and confined gas. Firstly, the effective viscosity of the confined gas decreases more rapidly than the bulk phase with decreasing pressure. Secondly, while the viscosity of the bulk gas is independent of pressure when the pressure is lower than 2 MPa, the viscosity of the confined gas decreases rapidly with decreasing pressure. The reasons for these differences are complex and will be further analyzed in the next section.

3.3. The Relationship between Viscosity and Knudsen Number (Kn)

The shear stress of fluid in the pore can be mathematically represented as equation (11).

There are three major components contributing to the above equation: kinetics contribution, potential contribution, and wall contribution.

In equation (12), which combines equation (11) with equation (6), the viscosity is divided into several different contribution groups.

In equation (12), there are six viscosity contributions that depended on temperature and density. When the density limits to 0, the effective viscosity of the dilute gas becomes . In this limit, the other terms on the right-hand side of equation (12) are zero because the potential shear stress becomes zero. and are the potential gas viscosities that are dependent on the potential shear stress between the gas molecules. and depend on the potential shear stress between the gas and the wall.

Because of the presence of the Knudsen layer, the effective viscosity of the dilute gas reduces and decreases with increasing the Knudsen number, which is expressed as a function of the Knudsen number equation (13). α is the viscosity coefficient reflecting the effect of the Knudsen effect on viscosity, which is affected by the shape of the pore and is taken as 2 for slit pores [18].

According to equation (12), the effective viscosity of dense gas can be expressed as equation (14). The detailed derivation is given in Appendix E.

Equation (16) shows that there is no direct functional relationship between the effective viscosity and the Knudsen number. Therefore, to clarify the effect of the Knudsen number on the effective viscosity, the relative viscosity is defined, μeff/μb. The relative viscosities of dilute and dense gases can be written as equation (15) and equation (16).

The semilog plot of relative viscosity obtained by the MD method and the Knudsen number is shown at Figure 7. The Knudsen number of methane can be acquired from Appendix F. As we know, when the value is below 0.001, it is a continuous flow, and when it falls between 0.001 and 0.1, a slip flow occurs. When the Knudsen number is greater than 0.1 but less than 10, it is a transition flow. A Knudsen number greater than 10 indicates that a free molecular flow is forming [10]. As shown in Figure 7, the results of relative viscosity with different width pores appear the same trend that the relative viscosity decrease with increase in Kn. Notably, the relative viscosity appears at different change rates with the Knudsen number in different regimes: when Kn is less than 0.1, relative viscosity slowly decreases with Kn increasing; however, effective viscosity decreases quickly with increase in Kn when Kn is between 0.1 and 10.

Michalis et al. [18] used equation (15) to depict the effective viscosity of confined dilute Nitrogen gas and got α = 2. As shown in Figure 7, it is noteworthy that the relative viscosity of the dense gas in the organic slit pores obtained from the MD simulations in this paper is lower than that of the dilute gas obtained from Michalis et al. [18]. It means that the steeper drop in the relative viscosity of the dense gas occurs at the position where the Knudsen number is smaller. Equation (15) could fit the MD results well. 2.901 is for α, while R2 = 0.990. The last term on the right side of equation (16) has a good correlation with the Knudsen number and tends to increase and then decrease. It is assumed that this term has the following equation (17) functional form.

α is assumed to be the same with Michalis et al. [18]. The relative viscosity of dense methane in organic matter slit pores was fitted with equations (16) and (17), and the results of the fit are shown in the black solid line of Figure 6 with coefficients a = 1.563, b = 2.143, and deterministic coefficient R2 = 0.993. This difference is because of the fact that the relative viscosity of the dense gas has one more contribution from the dense gas potential and the wall potential than that of the dilute gas.

3.4. Analysis of the Flow Capacity Improvement

The flow capacity is defined as permeability divided by viscosity. The effect of viscosity on the flow capacity could be described by μb/μeff, which is shown in Figure 8. The flow capacity increases as the pressure decrease. When the pressure is 1 MPa, the flow capacity is 28.7 times that calculated by the bulk gas viscosity in 1 nm pore.

4. Summary and Conclusions

Gas-bearing shale is rich in organic nanopores. This paper determined the viscosity of methane in organic slit nanopores by EMD simulations at shale gas reservoir conditions. The result shows that the viscosity of the shale gas is affected by the presence of nanopores. We obtained the restricted dense gas viscosity model by theoretical analysis. We use this model to fit the EMD results to obtain the viscosity correction formula for dense methane in the organic slit nanopore.

The conclusions drawn are as follows:(1)Confine spaces may decrease the effective viscosity of methane. In the pore, the effective viscosity of gas decreased with decreasing pore width, while in the smaller pores, pore width has more effect on effective gas viscosity. The confined effect would make the effective viscosity of methane decrease at least 10% for pore size less than 40 nm when the pressure is 4 MPa, which means that the confined effect on viscosity should be paid more attention to.(2)There are two differences in effective viscosity between bulk gas and confined gas. Firstly, the effective viscosity of the confined gas decreases more rapidly than the bulk phase with decreasing pressure. Secondly, while the viscosity of the bulk gas is independent of pressure when the pressure is lower than 2 MPa, the viscosity of the confined gas decreases rapidly with decreasing pressure.(3)The viscosity of a confined dense gas devotes by kinetic contribution, the potential of gas contribution, and the potential of wall contribution. A viscosity model of confined dense gas is promoted, and specific model parameters are obtained by fitting EMD results, which can be used to characterize the effective viscosity of dense methane under shale gas reservoir conditions.(4)The confinement effect causes the methane flow capacity in the reservoir to increase. As this pressure decreases, the effect of the confinement effect on the methane flow capacity increases rapidly. Keeping a low pressure can decrease the viscosity of shale gas, which is beneficial to improving the flow capacity of shale gas.

The confinement effects on the viscosity of shale gas should not be negligible, and viscosity correction is recommended for shale gas flow description. In this paper, a viscosity correction method for dense methane in organic matter slit pores is provided. The proposed method still has shortcomings, which could be improved later, such as the diverse types and complex morphology of real shale pores, however, the target of this paper has been achieved, which is only for organic matter slit pores.

Appendix

A. Chemical Potential

The chemical potential of dense methane is written as (A.1) [51] follows:

B. Redlich–Kwong Equation of State

The Redlich–Kwong equation of state (RK EOS), which is shown as equations (B.1) and (B.2), is used in calculating the density and fugacity of light hydrocarbons with different polarities (nonpolar, middle polar, and highly polar) [57, 58].

The compressibility factor of RK EOS is determined by equations (B.2) and (B.4).

For methane gas, TC is 190.6 K, pC is 4.599 MPa, molecular mass is 16.04 g/mol, and ω is 0.012 [59].

The fugacity coefficient of RK EOS is determined by equation (B.7) [60].

C. Chapman–Enskog Formula

There are numerous pieces of research on the viscosity of dilute gas, which contains the methodologies of bulk gas and gas in the pore [18, 61]. A formula to estimate the dilute gas viscosity at low pressure was proposed and shown in equation (C.1) [18].

D. Lee Function

The correlation about viscosity from Lee is a common method, which is shown in equations (D.1) to (D.2) [14].

E. The Effective Viscosity of Dense Gas in Slit Pore

By bringing equation (13) into (12) and splitting the second and third terms on the right-hand side of equation (12), equation (12) can be rewritten as equation (E.3).

Let

Then, the effective viscosity can be further rewritten as follows:

F. Knudsen Number

The Knudsen number is defined as equation (F.1) [10].

The characteristic length (L) of a pore that is used in the criterion to classify flow regimes is the radius (rp) for a circular pore and the width () for a slit pore.

The density of methane is calculated by equation (F.3). For a dense gas, the compressibility factor is identified by RK EOS. However, the compressibility factor is 1 for dilute gas.

Nomenclature

Notations
a:A parameter of the Redlich–Kwong equation of state
b:A parameter of the Redlich–Kwong equation of state
:Mean velocity of molecule in the system
d:Average intermolecular distance, Å/fs
h:Planck constant, 6.62607004 × 10−34 m2kg/s or J·s
k:Permeability of pore, m2
m:Molecular mass, g/molecule
n:The number density of gas, molecules/m−3
n0:The number density of gas at standard condition, molecules/m−3
p:Pressure of fluid, Pa
r:The vector between two molecules
rp:The radius of circular pore, m
t:Time, s
u:The speed of the stream layer, m/s
:The average flow rate in pore when it is in Darcy’s law, m/s
:The velocity of atom when it is in the Green–Kubo function, Å/fs
:The molar volume when it is in the Redlich–Kwong equation of state, m3/mol
:The width of slit pore, m
x:The position at flow direction, m
y:The location of the stream layer, m
A, B:The coefficients for equation of state
A:The cross-area of pore, m
F:The intermolecular forces, N
L:The physical width of the pore, m
L′:The effective width of pore, m
M:The total number of wall molecules
N:The total number of fluid molecules in the system
T:System temperature, K
U:The van der Waals potential between two atoms, Kcal/mol
V:System volume, m3
Z:A compressibility factor of real gas
Kn:The Knudsen number
:Relative molecular mass, g/mol
:Averaging for the system
Greek symbol
α:The coefficient of viscosity
δ:Rarefied gas criterion. When δ ≥ 7, the gas is in the dilute regime
ε:The potential well of Lennard–Jones potential model
λ:The mean free path, m
μ:The viscosity of fluid, Pa s
ρ:Density, g/cm3
σ:The molecular effective collision diameter, m
σT:Tangential accommodation coefficient
τ:Shear stress, N/m2
ϕ:Fugacity coefficient, fraction
ω:The acentric factor of fluid
Λ:de Broglie wavelength, m
Subscripts
b:The parameter of bulk phase
c:The parameter under the effect of confinement
C:The parameter at the critical point
i, j:Molecule number
p:The parameter of gas in pore
r:The reduced parameter
eff:The effective parameter
α, β:Coordinate direction
0:The parameter when Kn approaches 0
∞:The parameter when Kn approaches infinity
Constant
kB:Boltzmann constant, 1.381 × 10−23 m2·kg·s−2·K−1
R:Ideal gas constant, 8.314 J/(mol K).

Data Availability

The raw data required to reproduce these findings cannot be shared at this time as the data are also a part of an ongoing study. If someone is interested in our work, the result can be repeated by the details of the model, and the methods of data processing are described in the paper.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This work was supported by China Postdoctoral Science Foundation (No. 2019M650965), China National Major Projects (Nos. 2016ZX05010-002-002 and 2016ZX05048-002), and The Fund of SKL of Oil & Gas Reservoir Geology and Exploitation Engineering (No. PLN2019020). The first author would also like to acknowledge the China Scholarship Council (CSC) for providing financial support for him to study at the Faculty of Engineering and Applied Science of the University of Regina as a coeducation Ph.D. student.