A Fully Coupled Numerical Model for Microwave Heating Enhanced Shale Gas Recovery

Formation heat treatment (FHT) can be achieved by converting electromagnetic energy into heat energy (that is microwave heating or MWH). Experimental evidence shows that such FHT can significantly enhance oil and gas recovery. As relatively few research studies have been reported on microwave heating enhanced shale gas recovery (MWH-EGR), a fully coupled electromagnetic-thermo-hydro-mechanical (ETHM) model is developed for the MWH-EGR in the present study. In the ETHM model, a thermal-induced gas adsorption model is firstly proposed for shale gas adsorption and fitted by experimental data. This thermal-induced adsorption model considers the increase of matrix pore space due to the desorption of the adsorbed phase. Further, a thermal-induced fracture model in shale matrix is established and fitted by experimental data. Finally, this ETHM model is applied to a fractured shale gas reservoir to simulate gas production. Numerical results indicated that the thermal-induced fracturing and gas desorption make predominant contributions to the evolution of matrix porosity. The MWH can increase cumulative gas production by 44.9% after 31.7 years through promoting gas desorption and matrix diffusion. These outcomes can provide effective insights into shale gas recovery enhancement by microwave assistance.


Introduction
The exploitation of shale gas provides a powerful driving force for the development of oil and gas technology [1]. The shale gas reservoir is characterized by ultra-low matrix permeability, organic-rich deposition, and mineral-filled natural fractures. Current mainstream technologies available for commercial development of shale gas are horizontal drilling plus hydraulic fracturing, which can largely enhance reservoir permeability [2][3][4][5][6]. After hydraulic fracturing, the gas production rate is dramatically enhanced in the short term but this high production rate drops quickly in the longer term [7]. For conventional gas reservoirs, single hydraulic fracturing mainly generates great-scale fractures from the preexisting fractures. However, the reserves for tight gas reservoirs mainly exist in micro-scale pores [8]. Therefore, maintaining high gas production for the long term mainly depends on the release of gas stored in the matrix through some special techniques [9][10][11]. This study aims to explore an alternative method to improve both desorption and diffusion of gas in a shale matrix that combines a microwave heating treatment with hydraulic fracturing for the development of shale gas.
Formation heat treatment (FHT) is a process with great potential in enhancing gas production from unconventional oil/gas reservoirs [1,[12][13][14]. Different heating technologies, such as an electrical heater, hot water flooding, high-temperature steam, in-situ combustion, and electromagnetic heating (microwave/radio frequency) and so on, were applied to the oil and gas industry [4,6,[15][16][17]. It is asserted that the heat treatment of shale gas reservoirs can enhance permeability in a similar way for oil shales [3,13]. The combustion and pyrolysis can remove organic matter and decompose minerals from shale and further increase shale matrix permeability [4,6]. Jamaluddin et al. [12] confirmed an increment in permeability from 50% to 675% depending on the core type and treatment temperatures. The maximum temperature increment could reach around 650 • C. Gas-phase combustion for tight gas reservoirs may be an alternative for FHT. The temperature within the porous media can reach above 600 • C [12,13]. Meanwhile, the combustion enhancing recovery of shale gas was investigated through numerical simulations [3] and experiments [4,6]. The numerical simulations indicated that methane combustion did not induce a traveling wave of high temperature. Thus, methane combustion cannot be used for reservoir heating. Nevertheless, kerogen combustion (the oxidation of kerogen) can achieve a high temperature combustion wave. Sustaining a temperature wave that can reach 650 • C became possible to increase the permeability [3]. Chen et al. [6] investigated the impact of high temperature combustion on shale permeability enhancement. In their experimental study, the shale obtained from the Shanggu basin of Shangxi, China, was burned in an air environment at a temperature of 800 • C for 45 min to oxidize the organic matter and improve the physical properties of the shale. Kang et al. [1,8] investigated the change of porosity in the shale samples under a high temperature condition, up to the target temperature of 800 • C. They suggested that the application of FHT must be combined with multi-stage hydraulic fracturing in horizontal wells. The experience of Jamaluddin et al. [18] had great significance as a reference. Their patented downhole heater was a 60 kW heating system capable of inducing a mean temperature of 766 • C at a nitrogen flow rate of 3.5 m 3 /min [19]. In summary, FHT based methods have great potential for the recovery of oil and gas from oil or gas reservoirs.
In shale gas reservoirs, gases mainly include free gas in the porous space and adsorbed gas on the pore surfaces of the organic components. Some studies [20,21] estimate that the adsorbed gas shares are 20% to 85% of the total gas content. Therefore, enhancing gas desorption will inevitably increase shale gas recovery [9,10,22,23]. Existing experimental studies reported the effect of temperature on shale gas adsorption [24,25]. It is found that the adsorbed methane content in shale decreases with increasing temperature [25]. This indicates that thermal stimulation through changing the shale gas desorption behavior has the potential to enhance shale gas recovery, thus termed the thermal treatment technique. This thermal treatment technique enhances shale gas desorption and has been numerically investigated [9,10,22]. In the above numerical simulations, the fixed temperature boundary (over 200 • C) is generally used without considering the heat injection process. Wang et al. [23] used a 3D numerical model to explore the feasibility of an electrical heating style in enhancing shale gas desorption. In their model, the highest temperature reached 436 • C near the electrode when the electrical power was given as 500 kW. All these numerical studies only consider FHT induced gas desorption enhancement. No other mechanism in this FHT process is considered.
The FHT of shale gas reservoirs mainly induces the following mechanisms. Firstly, it accelerates the desorption of methane in the matrix, especially for organic-rich rocks [8][9][10]22,23]. Secondly, it prevents water blocking evaporating blocked water and dehydrating the clay structure [1,8]. Thirdly, it enhances matrix porosity and permeability or diffusivity in the microscale uniformly [1,8,26]. FHT is an auxiliary technique of hydraulic fracturing instead of refracturing. This can reduce the demand for more water, avoid the contamination of groundwater, and thus be a relatively eco-friendly technique [8]. As such, the heating method is the key to the success of FHT.
Microwave heating (MWH) can be a potential alternative thermal stimulation technique to FHT [1,15,27,28]. Microwave is a kind of electromagnetic wave with wavelength (1 mm-1 m) and frequency (300 MHz-300 GHz) [17,28]. It is a form of energy (not heat), but can be converted into heat in the medium. It can be transmitted, reflected, and absorbed by some materials. The polar molecules within shale reservoirs can absorb the microwave energy, which depends on the dielectric properties of the reservoirs. Dielectric materials exposed to microwave irradiation absorb and convert Energies 2018, 11, 1608 3 of 28 electromagnetic energy into heat energy. MWH does not rely only on convection or conduction, it can also heat objects internally regardless of the physical contact between the heating object and the microwave source [15,28].
Compared with conventional heating, MWH has the following advantages. Firstly, MWH is an energy conversion instead of heat transfer. Secondly, it is fast, non-contact, volumetric heating. Thirdly, it can start up and stop quickly. Finally, it has a higher level of safety and automation and is environmental friendly. A MWH technique based on a down-hole radiating antenna is less affected by formation geology and is able to distribute heat over a larger reservoir volume due to the propagation of electromagnetic energy in the medium [27]. MWH has been widely employed for oil shale [16,[29][30][31][32], oil sands [17], heavy oil [15,27,[33][34][35], coal-bed methane (CBM) recovery [28], tight sandstone gas [14,36], but little attention has been paid to shale gas recovery. In all these studies, the two-dimensional model for electromagnetic heating on heavy oil/CBM was solved using a simulator based on the finite element method [28,33,34]. It was noted that Li et al. [28] presented a fully coupled numerical model for microwave heating enhanced CBM recovery. They only focused on CBM desorption enhancement by MHW and did not consider thermal-induced coal fracturing. It is noted that gas flow and adsorption/desorption mechanisms have significant differences for CBM and shale gas, but few studies focus on the microwave heating enhanced shale gas recovery. Therefore, it is necessary to investigate the mechanisms of microwave heating-induced enhancement of shale gas recovery.
In order to investigate the feasibility of MWH-EGR, this paper proposes a fully coupled electromagnetic-thermo-hydro-mechanical model (termed ETHM). This model considers the interactions among microwave heating, reservoir deformation, gas transport in the shale matrix, fracture network, and hydraulic fracture. The thermal-induced interactions in reservoirs such as thermal expansion, temperature-dependent gas adsorption, and thermal-induced fracturing in the matrix are theoretically described and verified by some experimental data from the literature. The effects of the MWH stimulation are numerically evaluated and compared with the conventional drainage method through a 2D shale gas drainage example. These results are expected to provide a clearer understanding on the practical application of MWH to the enhancement of shale gas recovery. The outline of this paper is as follows: Section 2 discusses the interaction mechanisms and proposes a coupled ETHM model. Section 3 describes a numerical model for the coupled multi-physics process with multiple partial differential equations (PDEs). Section 4 quantifies the stimulation effects of MWH on shale gas recovery. The numerical simulations with the ETHM model are compared with those of the conventional drainage method. The findings and conclusions are presented in Section 5.

Electromagnetic Wave Excitation
The electromagnetic wave excitation is described by Maxwell's equations and boundary conditions. For a transient problem, the Maxwell's equations are Gauss law − magnetic : ∇·B = ρ Gauss law − electric : ∇·D E = 0 (4) where E is the electric intensity, D E is the electric displacement or electric flux density, H is the magnetic intensity, B is the magnetic flux density, J is the current density, and ρ is the electric charge density. Their constitutive relations are where ε 0 is the vacuum permittivity, ε r is the relative permittivity, ε r = ε + jε , ε is the dielectric constant that reflects the capacity of a material to store electromagnetic energy, and ε is the loss factor that dominates the conversion of electromagnetic energy into heat, µ 0 is the vacuum magnetic permeability, µ r is the relative magnetic permeability, and σ E is the electrical conductivity. After these constitutive models in Equation (5) are substituted into the governing equations, the Helmholtz vector equation for the time-harmonic electromagnetic field distributions is where k 0 = ω/c 0 is the free space wave number, j is the imaginary number, c 0 is the speed of light in vacuum, and ω = 2π f is the angular frequency. The microwave energy is weakened when propagating a lossy dielectric material and part of the microwave energy is converted into heat energy. This conversion is governed by the dissipated power P v (W/m 3 ) [28,34].

Energy Conservation Equation
If the inter-convertibility of thermal and mechanical energy is neglected, the energy conservation equation over an REV can be expressed as follows [37][38][39].
where T is the temperature, p f and ρ g f are the pressure and density of gas in natural fractures, respectively, ε v is the volumetric strain, K is the permeability tensor of natural fractures, µ is the dynamic viscosity of gas, K is the volumetric modulus, C eq and κ eq are the equivalent heat capacity and equivalent thermal conductivity. It is noted that the heat capacity of rock is at least 100 times larger than that of gas, and the formation porosity is normally small in shale reservoirs. Therefore, C eq and κ eq are dominated by the properties of formation rock [22]. K g and α g are the gas bulk modulus and the gas thermal expansion coefficient, respectively. As defined by Tong et al. [37], K g = p f and α g = 1/T, Q e denotes the total heat source. In the present study, this heat source is supplied by electromagnetic excitation through the electromagnetic loss. Therefore, the heat source is the sum of both resistive loss and electromagnetic loss (see Equation (9)), which are determined by the dielectric properties and magnetic permeability of reservoirs, respectively.
where E * and H * are the electric intensity and the magnetic intensity acting on the heated objects, respectively. Therefore, for the electromagnetic heating in porous media, Equation (8) becomes Energies 2018, 11, 1608

of 28
It is noted that the fourth term at the above right-hand side denotes magnetic loss. This loss is zero because shale is a non-magnetic material.

Reservoir Deformation Equation
Two categories of numerical models are usually used to consider the fracture deformation of poroelastic media: discrete and continuous approaches. Discrete approaches include discontinuities in the fracture modeling and continuous approaches consider fracture properties into the total deformation in the constitutive models [40][41][42]. For fractured reservoirs, the total elastic compliance tensor T e ijkl of cell e (where a cell refers to a matrix block) is contributed by both fractures and matrix as [7,43] T e ijkl = M e ijkl + C e ijkl (11) where M e ijkl is the isotropic elastic compliance tensor at the cell e, and C e ijkl is the anisotropic compliance tensor from the fracture.
where NFIE is the number of fractures intersecting the cell, K n and K s are the normal and shear stiffness of the fracture, respectively. δ ik is the Kronecker's delta, and l e is the trace length of the intersected fracture over this cell domain. F ij and F ijkl are non-dimensional tensors related to the fracture geometry, and ν is Poisson's ratio. Hence, the deformation equation is modified into where α B is the Biot's coefficient, p m the gas pressure in matrix, α T the volumetric thermal expansion coefficient, D the 4th order elasticity tensor which is the inverse of the total elastic compliance tensor considering fracture properties, I the identity tensor, u the displacement vector, F the body force.

Temperature-Dependent Adsorption Model
The shale gas adsorption depends on the content of organic carbon (TOC), organic matter type, thermal maturity, reservoir pressure, and temperature [9,23]. The Langmuir adsorption isotherm is widely applied to coalbed methane and shale gas [9]. This Langmuir isotherm can be extended to a temperature-dependent form based on the kinetic and thermodynamic theory [44].
K(T) is the temperature-dependent adsorption equilibrium constant.
where K 0 is a pre-exponential constant, E ads is the characteristic adsorption energy, R is the universal gas constant, ρ ga is the gas density in the standard condition, and V L is the Langmuir volume. The surface fitting of Equation (14) on available experimental data [24,25] is shown in Figure 1a for Longmaxi shale and in Figure 1b for Carynginia shale. These good fittings imply the applicability of Equation (14) to shale gas adsorption behavior under variable temperatures. In addition, the adsorbed phase occupies the following pore volume: where φ m is the matrix porosity, θ ads is the pore volume occupied by adsorbed phase, and ρ ads is the density of adsorbed phase.
K is a pre-exponential constant, ads E is the characteristic adsorption energy, R is the universal gas constant, ga  is the gas density in the standard condition, and L V is the Langmuir volume. The surface fitting of Equation (14) on available experimental data [24,25] is shown in Figure 1a for Longmaxi shale and in Figure 1b for Carynginia shale. These good fittings imply the applicability of Equation (14) to shale gas adsorption behavior under variable temperatures. In addition, the adsorbed phase occupies the following pore volume: (16) where m  is the matrix porosity, ads  is the pore volume occupied by adsorbed phase, and ads  is the density of adsorbed phase.  The temperature change can enlarge the porosity and increase the permeability and diffusivity of rock [1,6,8,26,39,45]. This may be induced by thermal-induced fracturing, which generates new pores and fractures. The description of this phenomenon is the key step in a thermal recovery model. It is found that the pore sizes in the shale matrix show significant fractal features [46,47]. These fractal features of the pores imply a relationship between the fractal dimension and the porosity of rock as [48]    The temperature change can enlarge the porosity and increase the permeability and diffusivity of rock [1,6,8,26,39,45]. This may be induced by thermal-induced fracturing, which generates new pores and fractures. The description of this phenomenon is the key step in a thermal recovery model. It is found that the pore sizes in the shale matrix show significant fractal features [46,47]. These fractal features of the pores imply a relationship between the fractal dimension and the porosity of rock as [48] where γ = d min p /d max p represents the ratio of maximum pore size to minimum one, and A is a constant. Therefore, the ratio of matrix porosity is expressed by where ς = γ/γ 0 , with the subscript "0" represents the initial state of corresponding variables.
In the present study, the field emission scanning electron microscopy (FE-SEM) images of the shale matrix under different high temperature treatments obtained by Kang et al. [1] were used to investigate the relationship between the temperature increment and fractal dimension. The FE-SEM images were firstly changed into binary images including matrix and pores only. Box-counting was then used to calculate the fractal dimension of pores in each image (see Figure 2). Finally, the relationship between the temperature increment and fractal dimension was fitted. A linear relation here is observed as (see Figure 3).
where T c is the room temperature, λ T is a sensitive coefficient of the fractal dimension to temperature, and D p0 is the fractal dimension of shale at the temperature T c .
where c T is the room temperature, Equation (20) shows that the change of porosity with the increase of temperature can be described by just one variable, in which the pore fractal dimension is only a middle variable. The sensitive Substituting Equation (19) into Equation (18), one can obtain Equation (20) shows that the change of porosity with the increase of temperature can be described by just one variable, in which the pore fractal dimension is only a middle variable. The sensitive coefficient, λ T , can be obtained by fitting available experimental data. Figure 4a,b present the fitting for the experimental data from Kang et al. [8] and Jha et al. [26], respectively. According to these fittings, the thermal-induced fracturing is closely related to the shale rock properties.
The general porosity response can be described as According to the constitutive relation including changes of temperature and pore pressure for the deformed shale reservoir, the volumetric strain is expressed as [38] ∆ε Substituting Equations (16), (18) and (22) into Equation (21) yields is the initial matrix porosity, and k s the solid bulk modulus of the grain.
Energies 2018, 11, 1608 9 of 28 Furthermore, the relationship between porosity and pore size is given by [22] φ m φ m0 where d p is the pore diameter within shale matrix and d p0 is the initial pore diameter.

K
Substituting Equations (16), (18) and (22) into Equation (21) yields  is the initial matrix porosity, and s k the solid bulk modulus of the grain. Furthermore, the relationship between porosity and pore size is given by [22]

Effect of Real Gas
The relation among gas density, gas pressure and temperature is described by the equation of state for real gas where g M is the molar mass of gas, g  and p respectively stand for gas density and gas pressure,

Effect of Real Gas
The relation among gas density, gas pressure and temperature is described by the equation of state for real gas where M g is the molar mass of gas, ρ g and p respectively stand for gas density and gas pressure, including gases in the shale matrix, natural fracture network and hydraulic fractures. The gas compressibility factor, gas bulk modulus, and viscosity change with temperature and pressure. The equation of state, such as the Peng-Robinson model, can be used to calculate the gas compressibility factor and viscosity. For simplicity, the present study calculates the gas compressibility factor by the following function of temperature and pressure [49]: where p r = p/p c and T r = T/T c are the relative pressure and the relative temperature of methane, respectively. p c and T c are the critical pressure and critical temperature of methane, respectively. The gas bulk modulus is then obtained as Gas viscosity is a function of both temperature and pressure as [50] where µ 1atm is the gas viscosity under atmospheric pressure and initial temperature. A 1 , A 2 and A 3 are constant. This study does not consider the gas fugacity.

Gas Flow in Matrix
The gas flow in the shale matrix is governed by a pseudo-steady diffusion equation as where ψ is the total gas mass in shale matrix. For a shale gas reservoir, this total gas mass consists of free gas and adsorbed gas. D app is the apparent diffusivity developed by our previous study [47,51], F s is the shape factor to reflect the geometry of matrix elements. It expresses the flow between matrix and fracture and depends on the form of multiple truncated fractures in the matrix element. The differentiation of φ m with respect to time is derived from Equation (23). Thus, the above governing equation can be obtained for gas diffusion in the matrix.

Gas Flow in a Natural Fracture Network
The governing equation for gas flow in a natural fracture network is written as where c g f is the bulk modulus of gas in natural fractures, φ f is the porosity of natural fractures, σ n and τ are the effective normal stress and shear stress of each fracture, respectively, and K ij is the permeability tensor in the whole domain, which is the integration of the K e ij of all cells in the domain.
where tensor P ij is related to the fracture geometry intersecting this cell in a 2D problem and is given as where w f is the fracture aperture, λ is a correction factor for the single parallel plate fracture model and is related to the fracture tortuosity and surface roughness. In the present study, λ = 1. l e w f /S e is a dimensionless factor for the revision of equivalent flux. The permeability tensor for the entire flow domain is obtained by integrating all permeability tensors in each cell. It is noted that the K e ij is overestimated by Equation (32) and should thus be revised. Further, the evolution of natural fracture aperture with stress should consider both normal closure and shear dilation [43].

Gas flow in Hydraulic Fracture
The gas flow in a hydraulic fracture is simulated using directional derivatives in COMSOL Multiphysics [52]. The governing equation of gas flow in a hydraulic fracture is where w F is the hydraulic fracture width, φ F the porosity of hydraulic fracture, ρ F gas density of gas in hydraulic fracture, and k F the proppant-pack permeability (i.e., permeability of hydraulic fracture). ∇ T denotes the gradient operator along the fracture, and Q m is the mass source term. p F is the gas pressure in hydraulic fracture. The same dependent variable of pressure is solved for Equation (30) in the fracture network.

Conceptual Design of the MWH-EGR Model
A conceptual design of the MWH-EGR model is presented in Figure 5a. This design is similar to those MWH for heavy oil and shale oil [1,27,32] and has the following three components: • A production well that is a multi-stage fractured horizontal well, whose completion scheme is specially designed in order to host microwave components and to allow electromagnetic irradiation. • An injection well whose completion scheme is specially designed in order to host the microwave components and to allow electromagnetic irradiation.

•
Microwave energy applicators which are composed of the surface unit with a high power microwave energy sources, the down-hole transmission lines and multiple bottom-hole waveguide antennas.

Numerical Model and Computational Parameters
A simplified 2D reservoir block is used to investigate MWH for shale gas recovery enhancement. A rectangular waveguide with TE10 mode or single propagating mode is used in this model. In this mode, if the model is in the xy-plane, the electric field only has a z-component, which furthermore is independent of the z-coordinate. This makes it possible to set up and solve the model in a 2D geometry with a single propagating mode. The main advantage with such a setup is that it is much faster and uses less memory. Figure 5b is a representative numerical model of a rectangular shale gas reservoir

Numerical Model and Computational Parameters
A simplified 2D reservoir block is used to investigate MWH for shale gas recovery enhancement. A rectangular waveguide with TE 10 mode or single propagating mode is used in this model. In this mode, if the model is in the xy-plane, the electric field only has a z-component, which furthermore is independent of the z-coordinate. This makes it possible to set up and solve the model in a 2D geometry Energies 2018, 11, 1608 13 of 28 with a single propagating mode. The main advantage with such a setup is that it is much faster and uses less memory. Figure 5b is a representative numerical model of a rectangular shale gas reservoir with dimensions of 100 m × 40 m. This numerical model contains a hydraulic fracture with a length of 60 m (the blue line). This reservoir is divided into ten matrix blocks with dimensions of 20 m × 20 m based on a previous study [7,43]. The total elastic compliance tensor and equivalent fracture permeability tensor in each matrix block are then obtained by Equations (11) and (31), respectively, after including the information of a discrete fracture in reservoir. Two microwave waveguides operated at the TE 10 mode are located at the two sides of the model. For the electromagnetic physics, port boundaries are applied on the waveguides for wave excitation, and the scattering conditions are applied on the remaining boundaries to make them transparent for microwaves [28]. The input power of the microwave is 100 W, and the microwave frequency is set as 2.45 GHz. The initial temperature of the reservoir is 350 K for the heat transfer process. No heat flux is given at any boundary. In the solid mechanics module, the left and lower boundaries are fixed in normal displacements and the in-situ stress states are acted on by the upper and right directions. For the gas flow, the initial gas pressures in the matrix, natural fracture, and hydraulic fracture are the same, p init . The external boundaries are set as no flow for gas transport, while the hydraulic fracture outlet is given a bottom-hole pressure. All computational parameters are taken from contemporary literature, experimental data or our estimations, and are summarized in Table 1. Table 2 presents the initial values and boundary conditions for all physical fields in this model.
Note: AB, BC, CD and DA denote the four boundaries of numerical model; EF and GH represent the waveguides; I is the intersection point between hydraulic fracture and horizontal well.

Coupled Multi-Physical Process
Equations (6), (10), (13), (29), (30), and (33) characterize the fully coupled ETHM model for MWH-EGR, including the Helmholtz equation, energy conservation equation, reservoir deformation equation, gas flow equations in the matrix, natural fracture, and hydraulic fracture. These equations are coupled through a suite of cross-coupling relations as shown in Figure 6. Note that the electromagnetic field provides the heat source for other fields and is only relevant to the heat transfer field. In fact, the temperature elevation may cause the changes in the dielectric properties of reservoir, which in turn affects the reservoir heating rate. Additionally, there may be also a coupling relationship between electromagnetic field and deformation field [54]. These couplings are not considered in this study. These cross-coupling relations link one process with others. The above governing equations are nonlinear partial differential equations (PDEs) and cannot be solved analytically. Alternatively, they are implemented into COMSOL Multiphysics, which provides the frequency-transient co-simulation solver as such that the electromagnetic field is solved in the frequency domain while the other physical fields remain in the time domain. It provides a powerful tool for the present study of MWH-EGR. Figure 7 shows the implementation procedure of the fully coupled electromagnetic-thermo-hydro-mechanical model in COMSOL with MATLAB. Firstly, the information of discrete fractures is generated based on the Monte-Carlo method. Secondly, the initial elastic compliance tensor and the improved permeability tensor are calculated in each cell. These two tensors consider both geometrical and mechanical properties of each fracture in this cell. Thirdly, the governing equations for the electromagnetic wave, heat transfer, reservoir deformation, the gas flow in the matrix and natural fractures, and the gas flow along the hydraulic fracture are simultaneously solved by the fully coupled solver within the platform of COMSOL software. Both elastic compliance tensor and improved permeability tensor are expressed by MATLAB code. The fracture aperture, normal stiffness, and shear stiffness (i.e., elastic compliance tensor and improved permeability tensor) are updated based on the numerical results at the previous time step. A new fully coupled process is restarted until the specified calculation time.

Coupled Multi-Physical Process
Equations (6), (10), (13), (29), (30), and (33) characterize the fully coupled ETHM model for MWH-EGR, including the Helmholtz equation, energy conservation equation, reservoir deformation equation, gas flow equations in the matrix, natural fracture, and hydraulic fracture. These equations are coupled through a suite of cross-coupling relations as shown in Figure 6. Note that the electromagnetic field provides the heat source for other fields and is only relevant to the heat transfer field. In fact, the temperature elevation may cause the changes in the dielectric properties of reservoir, which in turn affects the reservoir heating rate. Additionally, there may be also a coupling relationship between electromagnetic field and deformation field [54]. These couplings are not considered in this study. These cross-coupling relations link one process with others. The above governing equations are nonlinear partial differential equations (PDEs) and cannot be solved analytically. Alternatively, they are implemented into COMSOL Multiphysics, which provides the frequency-transient co-simulation solver as such that the electromagnetic field is solved in the frequency domain while the other physical fields remain in the time domain. It provides a powerful tool for the present study of MWH-EGR. Figure 7 shows the implementation procedure of the fully coupled electromagnetic-thermo-hydro-mechanical model in COMSOL with MATLAB. Firstly, the information of discrete fractures is generated based on the Monte-Carlo method. Secondly, the initial elastic compliance tensor and the improved permeability tensor are calculated in each cell. These two tensors consider both geometrical and mechanical properties of each fracture in this cell. Thirdly, the governing equations for the electromagnetic wave, heat transfer, reservoir deformation, the gas flow in the matrix and natural fractures, and the gas flow along the hydraulic fracture are simultaneously solved by the fully coupled solver within the platform of COMSOL software. Both elastic compliance tensor and improved permeability tensor are expressed by MATLAB code. The fracture aperture, normal stiffness, and shear stiffness (i.e., elastic compliance tensor and improved permeability tensor) are updated based on the numerical results at the previous time step. A new fully coupled process is restarted until the specified calculation time.

Cross-coupling processes for
Electromagnetic-Thermo-Hydro-Mechanical model  Additionally, the mesh size may influence simulation results to some extent. When the number of grids for the two waveguide ports is set as 16, 32, 64, and 128, the simulated minimum temperature is 124 °C , 127 °C , 128 °C , and 129 °C , respectively, and their respective maximum temperature is 298 °C , 309 °C , 315°C and 318°C . These results show that both the minimum temperature and the maximum temperature fall first and then rise, but gradually become stable when the number of the grid increases from 16 to 128. Here, the smallest grid size; i.e., 128-is used in this study to ensure the reliability of the simulation. In addition, the change in the mesh size of the overall model has little effect on the temperature for a given waveguide port mesh size (128 grids). Thus, an extra fine mesh is set in the computation.

Distribution of Electric Field and Magnetic Field
In the rectangular waveguide system, the transmission direction is along the z-axis, and the transmitted electromagnetic waves can be divided into three  Additionally, the mesh size may influence simulation results to some extent. When the number of grids for the two waveguide ports is set as 16, 32, 64, and 128, the simulated minimum temperature is 124 • C, 127 • C, 128 • C, and 129 • C, respectively, and their respective maximum temperature is 298 • C, 309 • C, 315 • C and 318 • C. These results show that both the minimum temperature and the maximum temperature fall first and then rise, but gradually become stable when the number of the grid increases from 16 to 128. Here, the smallest grid size; i.e., 128-is used in this study to ensure the reliability of the simulation. In addition, the change in the mesh size of the overall model has little effect on the temperature for a given waveguide port mesh size (128 grids). Thus, an extra fine mesh is set in the computation.

Distribution of Electric Field and Magnetic Field
In the rectangular waveguide system, the transmission direction is along the z-axis, and the transmitted electromagnetic waves can be divided into three categories according to the longitudinal components of the electric field E z and the magnetic field H z . If E z = 0 and H z = 0, E and H are in the cross-section of rectangular waveguides. This kind of electromagnetic wave is termed a transverse electromagnetic wave, which is simply termed the TEM wave. This wave cannot be solved by the longitudinal field method. If E z = 0 and H z = 0, only the electric field component is in the propagation direction. The magnetic field is only in the cross-section of rectangular waveguides and becomes a transverse magnetic wave (abbreviated as TM wave). If E z = 0 and H z = 0, only the magnetic field component is in the propagation direction, and the electric field becomes a transverse wave. This is abbreviated as the TE wave. The TE 10 module is termed the main wave form in the rectangular waveguide. The rectangular port is excited by a transverse electric (TE) wave that has no electric field component along with the direction of propagation. Given a microwave excitation frequency of 2.45 GHz, the TE 10 mode is only a propagating mode in the rectangular waveguide. Figure 8 shows the electric field distribution, E z , in two rectangular waveguide ports. It is noted that the propagation direction of electromagnetic wave is along x-axis here. The electric field distributions, which are similar to the sinusoidal shapes, are found for the left and right waveguide ports, respectively. The amplitude of the electric field is 9528.9 V/m. The varying electric field induces resistance loss and generates thermal energy. Figure 9a,b shows the magnetic field distribution (H x and H y ) in two rectangular waveguide ports, respectively. The maximum values of H x and H y are 0.84 A/M and 32.2 A/M, respectively. MWH includes two heating sources coming from the rapidly varying electric and magnetic fields: resistive loss and magnetic loss. The magnetic loss in shale reservoir does not contribute to temperature change because shale is regarded as a non-magnetic material and its relative magnetic permeability is unity. As such, the resistance loss dominates the whole MHW process. becomes a transverse wave. This is abbreviated as the TE wave. The TE10 module is termed the main wave form in the rectangular waveguide. The rectangular port is excited by a transverse electric (TE) wave that has no electric field component along with the direction of propagation. Given a microwave excitation frequency of 2.45 GHz, the TE10 mode is only a propagating mode in the rectangular waveguide. Figure 8 shows the electric field distribution, z E , in two rectangular waveguide ports. It is noted that the propagation direction of electromagnetic wave is along x-axis here. The electric field distributions, which are similar to the sinusoidal shapes, are found for the left and right waveguide ports, respectively. The amplitude of the electric field is 9528.9 V/m. The varying electric field induces resistance loss and generates thermal energy. Figure 9a The magnetic loss in shale reservoir does not contribute to temperature change because shale is regarded as a non-magnetic material and its relative magnetic permeability is unity. As such, the resistance loss dominates the whole MHW process.

Evolution of Temperature in the MWH Model
MWH is a process where high frequency electrical energy is transformed into heat energy by dielectric losses when an electromagnetic wave is radiated from the antennas into the shale gas reservoirs [33]. As the electromagnetic energy propagates into the formation, fluids and other reservoir materials prevent its passage by providing propagation resistance. As a result, the intensity of the propagating wave is reduced and the energy is converted into heat, thus increasing the formation temperature. This temperature change subsequently alters the gas desorption and the gas diffusion in reservoirs. When the electromagnetic wave propagates into the shale gas reservoir, a fixed portion of its residual energy at any point is delivered to a farther place in the reservoir. Therefore, most of the energy is deposited near the microwave source. Figure 10 depicts the energy deposition with the distance from the microwave source. The distance in Figure 10 is denoted as the "skin depth". It is the distance into a material the wave propagates before its field strength is reduced to 1 e ( 2.718  e ) of its initial value and the wave power is therefore reduced to about 37% of its original level [33,55,56]. Figure 11 presents the distribution of temperature at the production time of 9 1 10  s (approximately 31.7 years). As microwave heating time increases, the vicinity of the waveguides is first heated and the thermal energy is gradually transferred towards the front. After 31.7 years of MWH, the temperature close to the waveguides is greatly elevated to approximately 591 K (318 °C ), while the area away from the waveguide has a slight increase because of slow heat conduction and the minimum temperature finally reaches 402 K (129 °C ). The temperature difference between the maximum and minimum temperatures is about 189 °C . As MWH continues, the average temperature of the reservoir goes up while the temperature difference decreases. Figure 12  . They are at different distances from the right waveguide port. For point #1, the temperature starts to gradually elevate after three days and increases to 486 K (or 211 °C ) at the production time of 31.7 years. However, heat conduction is so slow that the temperature increment is confined to the vicinity of the two waveguide ports. The temperature at point #2 starts to increase after 200 days. The temperatures at points #3 and #4 increase slowly until the production time of 700 days. Therefore, how to uniformly elevate the temperature of the reservoir is the key issue to MWH during shale gas recovery. The ability that a material absorbs microwaves is mainly determined by its dielectric loss factor. Because the loss factor varies with material, MWH performs the characteristics of selective heating. An efficient method to improve the heating efficiency of microwaves is to increase the dielectric loss factor. The nanoparticles with high

Evolution of Temperature in the MWH Model
MWH is a process where high frequency electrical energy is transformed into heat energy by dielectric losses when an electromagnetic wave is radiated from the antennas into the shale gas reservoirs [33]. As the electromagnetic energy propagates into the formation, fluids and other reservoir materials prevent its passage by providing propagation resistance. As a result, the intensity of the propagating wave is reduced and the energy is converted into heat, thus increasing the formation temperature. This temperature change subsequently alters the gas desorption and the gas diffusion in reservoirs. When the electromagnetic wave propagates into the shale gas reservoir, a fixed portion of its residual energy at any point is delivered to a farther place in the reservoir. Therefore, most of the energy is deposited near the microwave source. Figure 10 depicts the energy deposition with the distance from the microwave source. The distance in Figure 10 is denoted as the "skin depth". It is the distance into a material the wave propagates before its field strength is reduced to 1/e (e = 2.718) of its initial value and the wave power is therefore reduced to about 37% of its original level [33,55,56]. Figure 11 presents the distribution of temperature at the production time of 1 × 10 9 s (approximately 31.7 years). As microwave heating time increases, the vicinity of the waveguides is first heated and the thermal energy is gradually transferred towards the front. After 31.7 years of MWH, the temperature close to the waveguides is greatly elevated to approximately 591 K (318 • C), while the area away from the waveguide has a slight increase because of slow heat conduction and the minimum temperature finally reaches 402 K (129 • C). The temperature difference between the maximum and minimum temperatures is about 189 • C. As MWH continues, the average temperature of the reservoir goes up while the temperature difference decreases. Figure 12 presents the evolution of the temperature at four specific points. Their coordinates are point #1 (98 m, 20 m), point #2 (90 m, 10 m), point #3 (70 m, 10 m), and point #4 (50 m, 10 m). They are at different distances from the right waveguide port. For point #1, the temperature starts to gradually elevate after three days and increases to 486 K (or 211 • C) at the production time of 31.7 years. However, heat conduction is so slow that the temperature increment is confined to the vicinity of the two waveguide ports. The temperature at point #2 starts to increase after 200 days. The temperatures at points #3 and #4 increase slowly until the production time of 700 days. Therefore, how to uniformly elevate the temperature of the reservoir is the key issue to MWH during shale gas recovery. The ability that a material absorbs microwaves is mainly determined by its dielectric loss factor. Because the loss factor varies with material, MWH performs the characteristics of selective heating. An efficient method to improve the heating efficiency of microwaves is to increase the dielectric loss factor. The nanoparticles with high ability to absorb microwaves, such as metal nanoparticles, iron oxide nanoparticles, can be taken into formation by the fracturing fluids during hydraulic fracturing [32,35]. The change of temperature impacts other physical processes. Thus, the contributions of each thermal effect (such as thermal expansion, thermal fracturing, and thermal-induced gas desorption) to the variation of porosity will be studied in the following section.
Energies 2018, 11, x FOR PEER REVIEW 21 of 30 ability to absorb microwaves, such as metal nanoparticles, iron oxide nanoparticles, can be taken into formation by the fracturing fluids during hydraulic fracturing [32,35]. The change of temperature impacts other physical processes. Thus, the contributions of each thermal effect (such as thermal expansion, thermal fracturing, and thermal-induced gas desorption) to the variation of porosity will be studied in the following section.  ability to absorb microwaves, such as metal nanoparticles, iron oxide nanoparticles, can be taken into formation by the fracturing fluids during hydraulic fracturing [32,35]. The change of temperature impacts other physical processes. Thus, the contributions of each thermal effect (such as thermal expansion, thermal fracturing, and thermal-induced gas desorption) to the variation of porosity will be studied in the following section.

Evolution of Matrix Porosity and Diffusivity
Our previous study [7] revealed that shale gas production mainly depends on matrix diffusivity and the matrix porosity dominates this change of matrix diffusivity process. Equation (23) shows that the porosity ratio represents the integrated effect of volumetric strain, gas pressure, gas desorption, thermal expansion, and thermal-induced fracturing. Figure 13a,b describe the evolutions of these components in contributing to the porosity variation in points #1 and #4, respectively. The contributions in the thermo-hydro-mechanical interactions vary distinctly at two different locations. The volumetric strain and pore pressure never significantly affect the porosity (see Figure 13a). When the adsorption layer is taken into account, the matrix porosity rapidly grows, indicating that the pore space released by the high temperature-induced gas desorption is an important factor to porosity enlargement. The immediate vicinity of the waveguide port (point #1) is characterized by a steep temperature gradient. The pore growth induced by thermal fracturing significantly increases the matrix porosity while thermal expansion slightly decreases the matrix porosity. For the points away from the waveguide port as shown in Figure 13b, the thermal effect exerts no influence on the porosity at the incipient stage of gas drainage. As time elapses, the temperature of point #4 progressively rises. The integrated factors induce a monotonic increase of porosity. Although thermal expansion is gradually enhanced, the thermal-induced fracturing and desorption enlarges considerably due to the high temperature change. As a result, the porosity constantly increases. The results indicated that the thermal-induced fracturing and gas desorption make the predominant contribution to the evolution of matrix porosity. Figure 14 shows the evolution of the matrix diffusivity ratio at four observation points. It is found that their matrix diffusivity increases considerably with the extraction time and the observation point near the waveguide exhibits larger change. The rapid variation in temperature can also induce the change of matrix porosity. The physical fields presented in the reservoir are strongly related and interact with each other.

Evolution of Matrix Porosity and Diffusivity
Our previous study [7] revealed that shale gas production mainly depends on matrix diffusivity and the matrix porosity dominates this change of matrix diffusivity process. Equation (23) shows that the porosity ratio represents the integrated effect of volumetric strain, gas pressure, gas desorption, thermal expansion, and thermal-induced fracturing. Figure 13a,b describe the evolutions of these components in contributing to the porosity variation in points #1 and #4, respectively. The contributions in the thermo-hydro-mechanical interactions vary distinctly at two different locations. The volumetric strain and pore pressure never significantly affect the porosity (see Figure 13a). When the adsorption layer is taken into account, the matrix porosity rapidly grows, indicating that the pore space released by the high temperature-induced gas desorption is an important factor to porosity enlargement. The immediate vicinity of the waveguide port (point #1) is characterized by a steep temperature gradient. The pore growth induced by thermal fracturing significantly increases the matrix porosity while thermal expansion slightly decreases the matrix porosity. For the points away from the waveguide port as shown in Figure 13b, the thermal effect exerts no influence on the porosity at the incipient stage of gas drainage. As time elapses, the temperature of point #4 progressively rises. The integrated factors induce a monotonic increase of porosity. Although thermal expansion is gradually enhanced, the thermal-induced fracturing and desorption enlarges considerably due to the high temperature change. As a result, the porosity constantly increases. The results indicated that the thermal-induced fracturing and gas desorption make the predominant contribution to the evolution of matrix porosity. Figure 14 shows the evolution of the matrix diffusivity ratio at four observation points. It is found that their matrix diffusivity increases considerably with the extraction time and the observation point near the waveguide exhibits larger change. The rapid variation in temperature can also induce the change of matrix porosity. The physical fields presented in the reservoir are strongly related and interact with each other.

Distribution of Gas Content within the Shale Matrix
The distribution of gas content within the shale matrix was investigated. Figure 15a,b are free gas content and adsorbed gas content with/without MWH, respectively. The adsorbed gas contents with MWH in each production time are significantly lower due to high desorption ratio under high temperature. However, the free gas contents with MWH in each production time are slightly higher, which result from rapid desorption and slow diffusion process. Without MWH, the adsorbed gas content in the shale matrix changes very slowly and this change depends on the shape factor of matrix block. When the reservoir is heated by microwave irradiation, the gas content decreases not only in the matrix block with a high shape factor but also in the temperature affected area. The larger the

Distribution of Gas Content within the Shale Matrix
The distribution of gas content within the shale matrix was investigated. Figure 15a,b are free gas content and adsorbed gas content with/without MWH, respectively. The adsorbed gas contents with MWH in each production time are significantly lower due to high desorption ratio under high temperature. However, the free gas contents with MWH in each production time are slightly higher, which result from rapid desorption and slow diffusion process. Without MWH, the adsorbed gas content in the shale matrix changes very slowly and this change depends on the shape factor of matrix block. When the reservoir is heated by microwave irradiation, the gas content decreases not only in the matrix block with a high shape factor but also in the temperature affected area. The larger the

Distribution of Gas Content within the Shale Matrix
The distribution of gas content within the shale matrix was investigated. Figure 15a,b are free gas content and adsorbed gas content with/without MWH, respectively. The adsorbed gas contents with MWH in each production time are significantly lower due to high desorption ratio under high temperature. However, the free gas contents with MWH in each production time are slightly higher, which result from rapid desorption and slow diffusion process. Without MWH, the adsorbed gas content in the shale matrix changes very slowly and this change depends on the shape factor of matrix block. When the reservoir is heated by microwave irradiation, the gas content decreases not only in the matrix block with a high shape factor but also in the temperature affected area. The larger the area in which the temperature changes, the larger the area in which the adsorbed gas content decreases.   Figure 16a,b show the evolution of the free gas and adsorbed gas at four specific points, respectively. The free gas at point #1 increases 6.5% of initial value firstly and decreases with time. This change is caused by the rapid increase in temperature near the right waveguide. No increase in free gas is observed at other observation points due to their slow rise in temperature. Therefore, the desorbed gas can flow into the natural fractures by diffusion. The desorption ratios are 44.4%, 34.2%, 25.0%, and 21.4% from point #1 to #4 after 31.7 years, respectively. Figure 17 compares the evolution of adsorbed gas and free gas with/without MWH at points #1 and #4, respectively. The difference in free gases with and without MWH at point #1 is gradually narrowing with time. The desorption ratio at Point #1 increases 434.9% due to MWH. The gap between free gases at point #4 with and without MWH is gradually enlarging with time. As the temperature rises, the gas continuously desorbs, and the increase in matrix diffusion is insufficient to deliver the gas to the natural fractures as soon as possible, which is the main reason for this phenomenon. Figure 18 indicates that the total gas content including adsorbed gas and free gas decreases significantly at point #1 due to a larger change in temperature. It is clear that the total gas content in the MWH-EGR case is lower than that in the routine shale gas recovery case, especially near two waveguide ports. The total gas contents at point #1 are 99.3%, 98.1%, 92.9%, and 84.8% of those in the routine shale gas recovery case for the drainage time of 1.3, 3.2, 12.6, and 31.7 years, respectively. The total gas contents at point #1 are 97.1%, 94.5%, 83.3%, and 66.9% of those in the MWH-EGR case for drainage times of 1.3, 3.2, 12.6 ,and 31.7 years, respectively. These results indicate that MWH enhances gas recovery significantly. The thermal-induced increase of matrix porosity accelerates the migration of desorbed gas from the matrix to natural fractures, which is preferred relative to gas desorption in MWH-EGR.  Figure 16a,b show the evolution of the free gas and adsorbed gas at four specific points, respectively. The free gas at point #1 increases 6.5% of initial value firstly and decreases with time. This change is caused by the rapid increase in temperature near the right waveguide. No increase in free gas is observed at other observation points due to their slow rise in temperature. Therefore, the desorbed gas can flow into the natural fractures by diffusion. The desorption ratios are 44.4%, 34.2%, 25.0%, and 21.4% from point #1 to #4 after 31.7 years, respectively. Figure 17 compares the evolution of adsorbed gas and free gas with/without MWH at points #1 and #4, respectively. The difference in free gases with and without MWH at point #1 is gradually narrowing with time. The desorption ratio at Point #1 increases 434.9% due to MWH. The gap between free gases at point #4 with and without MWH is gradually enlarging with time. As the temperature rises, the gas continuously desorbs, and the increase in matrix diffusion is insufficient to deliver the gas to the natural fractures as soon as possible, which is the main reason for this phenomenon. Figure 18 indicates that the total gas content including adsorbed gas and free gas decreases significantly at point #1 due to a larger change in temperature. It is clear that the total gas content in the MWH-EGR case is lower than that in the routine shale gas recovery case, especially near two waveguide ports. The total gas contents at point #1 are 99.3%, 98.1%, 92.9%, and 84.8% of those in the routine shale gas recovery case for the drainage time of 1.3, 3.2, 12.6, and 31.7 years, respectively. The total gas contents at point #1 are 97.1%, 94.5%, 83.3%, and 66.9% of those in the MWH-EGR case for drainage times of 1.3, 3.2, 12.6 ,and 31.7 years, respectively. These results indicate that MWH enhances gas recovery significantly. The thermalinduced increase of matrix porosity accelerates the migration of desorbed gas from the matrix to natural fractures, which is preferred relative to gas desorption in MWH-EGR.

Cumulative Production Enhancement Due to MWH
The effect of MWH on enhancing cumulative production is investigated. Figure 19 compares the cumulative productions with conventional shale gas recovery and MWH enhanced shale gas recovery. It clearly shows that MWH greatly enhances the cumulative gas production. To clearly describe the stimulation effect of MWH, a line for yield-increasing efficiency is also presented in Figure 19. MWH promotes the cumulative gas production by 11.5% after 3.2 years, 24.9% after 12.6 years, and 44.9% after 31.7 years. This yield-increasing efficiency significantly increases with production time. The MWH-EGR is mainly reflected by extending a mid-late period of gas production and effectively solving the low production induced by a low desorption ratio and diffusion capacity.

Cumulative Production Enhancement Due to MWH
The effect of MWH on enhancing cumulative production is investigated. Figure 19 compares the cumulative productions with conventional shale gas recovery and MWH enhanced shale gas recovery. It clearly shows that MWH greatly enhances the cumulative gas production. To clearly describe the stimulation effect of MWH, a line for yield-increasing efficiency is also presented in Figure 19. MWH promotes the cumulative gas production by 11.5% after 3.2 years, 24.9% after 12.6 years, and 44.9% after 31.7 years. This yield-increasing efficiency significantly increases with production time. The MWH-EGR is mainly reflected by extending a mid-late period of gas production and effectively solving the low production induced by a low desorption ratio and diffusion capacity. Cumulative gas production without MWH Cumulative gas production with MWH Yield-increasing efficiency Figure 19 Cumulative gas productions without and with MWH.

Conclusions
In the present study, a fully coupled ETHM model was proposed to include the interactions among the electromagnetic wave transmission, the mechanical deformation of reservoir, the thermal transfer, and the multi-scale gas flow within the shale reservoir. Based on the coupled model, a simplified 2D reservoir model was numerically studied on the interactions of different mechanisms

Conclusions
In the present study, a fully coupled ETHM model was proposed to include the interactions among the electromagnetic wave transmission, the mechanical deformation of reservoir, the thermal transfer, and the multi-scale gas flow within the shale reservoir. Based on the coupled model, a simplified 2D reservoir model was numerically studied on the interactions of different mechanisms and the effects of MWH-EGR. From these results of the present study, the following conclusions can be drawn: Firstly, the proposed ETHM model can better forecast the MWH effects on shale gas reservoirs. MWH can significantly enhance the cumulative gas production. The yield-increasing efficiency increases with production time. The present numerical simulations reveal that the MWH-EGR promotes the cumulative gas production by 11.5% after 3.2 years, 24.9% after 12.6 years, and 44.9% after 31.7 years, respectively. The simulated results can provide new insights into microwave-assisted shale gas recovery.
Secondly, with the MWH proceedings, the temperature increases preferentially in the vicinity of the waveguide ports. The change of temperature in the region away from the waveguide depends on reservoir properties, such as thermal conductivity, dielectric properties and so on. A significant overheating phenomenon occurs near the waveguide ports, which limits the further increase of temperature. Therefore, how to elevate temperature uniformly is the key issue of MWH-EGR.
Thirdly, the temperature elevation accelerates both gas desorption and gas diffusion. At a specific point near the waveguide port, high temperature enhanced gas desorption and released the pore space occupied by the original adsorbed gas. Thermal fracturing greatly increased the porosity and diffusivity of the shale matrix, but the thermal expansion of the shale matrix slightly decreased the porosity and diffusivity of the shale matrix. The increase in matrix porosity depends on heating rate.
Fourthly, the gas content decreases not only in the matrix block with high diffusivity but also in the thermal affected area due to the thermal-enhanced gas desorption. Gas desorption process is rapidly enhanced with the rapid temperature elevation in the vicinity of the waveguides. The free gas content rise first and then falls. The increase in temperature does not only enhance the shale gas desorption, but also improves the diffusivity of the shale matrix to a certain extent. The MWH-EGR mainly extends a mid-late period of gas production and can effectively solve the low production induced by a low desorption ratio and diffusion.
Finally, the numerical simulation analysis is done for a two-dimensional domain. This 2D model provides useful insights into MWH-EGR. Such a 2D model cannot describe the temperature distribution in the vertical direction of the reservoir. The heat loss in the vertical direction can occur in the MWH process through overburden and the underlying layers. In this respect, a 3D model is more accurate in the simulation of MWH-EGR. Besides, MWH-EGR may significantly remove the water blocking due to temperature elevation, thus enhancing gas recovery. The permittivity, thermal and electrical conductivity may also impact the MHW-EGR process dramatically. These should be studied further.