The purpose of this paper was to present the thermoporoelasticity model adapted for application in modelling processes, where phase transition may occur, such as during underground coal gasification (UCG). The mathematical model of the medium (soil/rock with pores filled with liquid/gas) in non-isothermal conditions is based on Biot's poroelasticity model. The poroelasticity model is expanded here by the influence of temperature and adjusted to the case where both liquid and highly compressible fluid are present in pores by using the gas laws. This requires considering temperature-dependent physical quantities such as pore fluid density, heat transfer coefficient and viscosity as functions of temperature. Based on the proposed mathematical model and the finite element method, a numerical model was built for the purpose of computing processes occurring in the vicinity of the UCG generator. The result of the authors’ work is a three-dimensional (3D) model, which was not only modified, but derived straight from the laws of thermodynamics, where fields of displacement, temperature and fluid flow are coupled. The model makes it possible to determine results significant to modelling of the UCG process, the reach of the gaseous phase's presence in pores, subsidence values, temperature distribution and directions and rate of seepage, without losing the simplicity and elegance of Biot's original concept. Next, the results of simulations for a hypothetical deposit to estimate the environmental impact of UCG are presented. After applying specific geometry and parameters, the model can be useful for verifying if the chosen technology of UCG in specific conditions will be safe for the environment and infrastructure.
- compressible pore fluid
- thermo-hydro-mechanical modelling
- underground coal gasification
The process of underground coal gasification, or UCG, is a technology for obtaining synthetic gas from deposits of hard coal or lignite
There are several UCG technologies; in general, they can be divided into the shaft (requiring human work underground) and shaftless (not requiring human work underground) methods. In shaftless methods (linked vertical wells, CRIP), there are usually drilled 2 boreholes (injection and production well), between which the UCG process is carried out. In modern methods, such CRIP (controlled retractable injection point) or P-CRIP (parallel CRIP), directional drilling is involved, and the point where the oxidant (air or oxygen with steam) is injected moves with the progress of the process, thus ensuring stable syngas quality. During UCG, processes take place that are analogous to those during traditional coal gasification over ground. Chemical reactions occur in the coal seam between coal and gases, and between the gases. Finally, the synthesis gas consists mainly of hydrogen, carbon oxides, methane and nitrogen. The proportions of these constituents will vary according to the seam's properties, the gasifying agent (oxidant), the amount of water used in the process and technology. A schematic of the UCG process is shown in figure 1.
In comparison to traditional deposit extraction methods, it is characterised by greater safety for the people employed during operation as well as less dust and noise pollution. It also makes it possible to reduce operating costs, including those generated by the transport and storage of coal, and does not require the storage of waste and ash . However, it poses a risk of polluting groundwater and air . Geotechnical issues involving analysis of the influence of the gasification process on the surrounding environment are an important part of the decision to undertake a project of this type, or whether to abandon activities in this respect. The influence of the installation on the environment must be carefully analysed in order to demonstrate that the given installation has no adverse impact on the natural environment .
In the gasification process, high temperatures are generated in the coal layer as a result of exothermic reactions, and this layer interacts strongly with the surrounding geological environment. The most important aspects from the perspective of environmental impact are subsidence, directions and rates of seepage and temperature distribution in the rock mass. The model used for simulations must account for coupled mechanical and thermal interactions as well as the possibility of phase changes in the fluid filling the medium's pores [10,11,12,13,14]. This article presents a coupled thermo-hydro-mechanical (THM) model based on Biot's poroelasticity model [15, 16].
The 3D consolidation theory, presented in isothermal form by Biot  and known as the poroelasticity theory, uses constitutive equations that make the strain of the medium and fluid discharge-dependent on stresses in both phases. It is applicable not only to soils or rocks, but also to porous media in general, including anisotropic media  such as concrete  or certain biological tissues [20, 21]. In the poroelasticity theory, the skeleton is treated as an elastic medium in which stress–strain relationships can be described by means of Hooke's law, and the pore pressure also has an influence on the skeleton's strain. The classical poroelasticity theory has been modified multiple times for the purpose of accounting for additional phenomena. Coussy [22, 23] uses changes in porosity in constitutive relationships instead of changes in fluid content. In addition, he distinguishes porosity in Lagrangian and Eulerian descriptions. Another form of the model is presented by Derski  and Strzelecki , who use the change in fluid volume instead of fluid content in the pores of the medium.
Using the poroelasticity model as a basis, there are also coupled THM models being built in various forms, depending on the phenomena analysed by the authors of these models. In Coussy's model , the state parameter used for creation of the model is porosity, but the final constitutive equations of the medium are transformed for using fluid content (depending on the volumetric strain of the skeleton, fluid pressure and temperature), so the model's structure is quite complicated. The model of Gawin et al.  focuses more on the flow and heat transfer of two immiscible fluids than on deformations. Coupling between the solid and fluid phase is implemented there by the Biot constant in the constitutive equation of the medium, which constant takes volumetric strain into account in both phases of the medium. In the heat conduction equations, stress is put on the effects related to the presence of 2 immiscible fluids; deformations of constituents are not taken into account there. The model of Bartlewska–Urban and Strzelecki  takes the rheological properties of the skeleton into account; it does not consider the presence of a gaseous phase in the pores of the medium. There is also a thermal seepage model  based on Biot's poroelasticity model, which takes volumetric strain into account. Liu et al.  presented the THM model with coupling implemented by the Biot constant for shale gas recovery, proving that gas recovery is strongly dependent on temperature.
In the case of elastic media, dependencies between temperature, strain and stress are described by Duhamel–Neumann relations. The method for deriving these relations was presented by Nowacki . The corresponding relations for poroelastic media can be derived analogously. In combination with the equation of fluid motion, they make up the thermoporoelasticity model. Numerical models applying different variants of the thermoporoelasticity model have been published, for example, in articles concerning geothermal springs , steam flooding of heavy oil reservoirs , the safety of storing radioactive waste , the influence of temperature changes on engineering structures , cement hydration  and ceramic filters . In this article, a proposal for a method of use for the thermoporoelasticity model, based on Derski's approach, for modelling coupled phenomena occurring in the area of coal gasification will be presented, using the example of a hypothetical brown coal deposit. The presented model concerns changes occurring in the vicinity of the reactor; process parameters such as gasification rate, temperature and pressure in the reactor are not modelled. Here, they constitute boundary conditions. A review of models used to simulate other processes occurring during UCG and their classification can be found elsewhere, for example, in the work of Rosen et al. .
Equations of flow continuity in the thermoporoelasticity model are the same as in the classical poroelasticity model. The medium's equation of continuity is in the form :
In the case of poroelastic media, mechanical coupling of phases in the conservation of momentum law is, for each of them, included in the form of kinetic energy and energy dissipation as a result of the fluid's flow through pores. The conservation of momentum law for the medium's solid phase is in the form:
Non-isothermal models of the biphasic medium are divided into two categories: LTE (local thermal equilibrium) – where the influence of the medium's mean temperature is accounted for in constitutive relations, and LTNE (local thermal non-equilibrium), where we consider the influence of the skeleton's and fluid's temperatures separately. In the event of an identical temperature for both phases, the results of both models are the same. The authors created both LTE and LTNE models, which were described and discussed in more detail in the works [15, 41, 42]. In this article, due to the scale of the phenomenon, the numerical model was based on the LTE approach. The procedure for deriving constitutive equations and heat transfer equations for the model are presented below in shortened form. Constitutive relations for the poroelastic medium under non-isothermal conditions can be obtained by accepting Helmholtz free energy
Heat transfer equations were obtained by using the entropy equation. The
As an example, we use one of deposits of brown coal, which is not being exploited due to geological conditions and environmental issues [45, 46]. It was assumed that part of the brown coal deposit serving as an example, with a thickness of 3 m, deposited at a depth of approximately 120 m below the ground level, would be gasified. Layers of clay below and above the layer of coal are insulating the deposit, and the model contains 8 geological layers, with boundaries whose ordinates are known at 16 points distributed over a 500-m mesh. The model's geometry is shown in figure 2.
A plane within the sandstones was set as the “0” ordinate, being the bottom of the model. It was assumed that the area in which gasification takes place is found in the central part of the model and has dimensions of 10 m * 100 m (figure 3).
Examples based on data from the literature concerning the deposit  and material tables containing the mechanical and thermal parameters of soils were adopted for simulations. After the end of the gasification process, it was assumed that the coal's mechanical parameters would decline. The parameters adopted for simulations are given in tables 1 and 2.
Material parameter values of individual geological layers.
initial values, which gradually fall to 0.1% of the initial value in the area of the gasification generator due to loss of mass of coal during the process
Thermal parameters of individual pore fluids and medium.
Volumetric expandability is given for fluids, linear for the skeleton.
Temperature-dependent value, initial value (for 0°C)
The medium's heat transfer coefficient and specific heat depend on the phase of the pore fluid and, in addition, the parameters of a given phase change with temperature. The permeability coefficient is dependent on viscosity: (
Similar to the case with viscosity, Sutherland's formula describes the dependency of gas thermal conductivity with respect to temperature:
The numerical model for FEM (finite element method) simulations was created using FlexPDE v.6 software. It was assumed that water is initially the pore fluid throughout the entire area, and as temperature grows, it undergoes phase transformation in certain areas. To determine the phase transformation temperature depending on the pressure in the medium, the Clausius–Clapeyron equation was used:
The use of different constitutive relations and parameters within a single model was made possible by the construction of a script using a coefficient defining the current state of matter, equal to 1 or 0, by which values pertaining to the given state of matter were multiplied (both values were computed), for example, in the case of fluid dilatation:
Model was run in 3 versions:
For the space, where gas is present in the pores, all mentioned modifications were applied: constitutive equations (10), conductivity and thermal parameters. Constitutive equations without modification (8), change in conductivities applied. No modifications in constitutive equations nor other parameters.
For the space, where gas is present in the pores, all mentioned modifications were applied: constitutive equations (10), conductivity and thermal parameters.
Constitutive equations without modification (8), change in conductivities applied.
No modifications in constitutive equations nor other parameters.
Our basic version is A. All the figures concern version A, unless version B or C is marked in the description. Versions B and C are used for comparison and identification of impact of given modifications.
Strain: The absence of displacements in the model's bottom plane and the absence of horizontal displacements in side planes was assumed. The condition of conformity of horizontal and vertical displacements was accepted on individual dividing surfaces between layers. Initial displacements are equal to zero, and initial stress results only from the mass of overburden.
Pressure: The absence of fluid flow through the bottom and side surfaces of the model was assumed, and the condition of conformity of fluid pressure on surfaces dividing layers was adopted.
Temperature: The known temperature at ground level and on the side and bottom planes was accepted, changing linearly depending on depth: from 10°
A gasification process scenario based on descriptions of coal gasification in the literature [7, 47] was adopted for the coal gasification generator. 2 years was accepted as the duration of the coal gasification period. Temperature increases gradually, as the UCG generator is accepted as a heat source. Then, as the gasification process progresses along the coal layer, heat is being dissipated. The evolution of temperature in exemplary points (x = 715, 745 and 775m) belonging to the gasified area (x = 700 to 800m) is shown in fig. 5. It was assumed that pressure reduced with respect to hydrostatic pressure will be maintained in the reactor over the course of the gasification process – this is one of approaches to UCG, applied to ensure that water is present in the process and to prevent the escape of gas and pollution [9, 48].
Displacements and pressure distribution of pore fluid resulting from self-weight were adopted as the initial condition, and temperature distribution changes linearly from 10°
Conducted numerical simulations made it possible to obtain fields of the medium's displacement, of stresses in the fluid and of temperature. This made it possible to determine the directions and rate of seepage, thermal flux and stresses in the rock mass. Thanks to the knowledge of the distribution of temperature and pressure in the rock mass, it was possible to determine the reach of gas presence in the medium's pores. In addition, numerical simulations made it possible to analyse temperature- and/or pressure-dependent parameters such as fluid density, conductivity coefficient and heat transfer coefficient. The FlexPDE program was chosen for the calculations because it enables any differential equation to be solved, and the number of simultaneously solved equations is not limited; it enables vector and array variables and is able to control timestep and automatically refines mesh – it monitors accuracy, finds areas with high error values and refines their mesh (user can also force mesh refining with a function – here on the front defined by the Clausius–Clapeyron equation). During simulations, the number of nodes ranged from 42,752 to 88,390 and the number of elements from 29,815 to 63,185 due to the mesh refinement in the area of interest.
The temperature distribution in the surroundings of the gasified deposit at certain instances of time is given in fig. 6a–f. As can be seen, the temperature maximum moves along the gasified layer of coal over time. The rock mass cooling process is slow, and elevated temperature around the deposit persists for many years after the process is completed.
Phase transition occurs as a result of the elevated temperature in the vicinity of the gasification area. The reach of phase transition is greatest at time 1.5 year – its reach is approximately 10 m outside the generator. It is possible to identify this region also in versions B and C, but in versions A and B, the temperature differences are small, and in version C, it is not necessary, because parameters do not change. Fig. 7a–f presents the area where water vapour is present at points in time corresponding to those given in fig. 6. A spatial view of the situation of the area in which gas (water vapour) is present in the medium's pores is shown in fig. 8 (example at t = 1 year).
The vector field of displacements resulting from elevated temperature is shown in fig. 9a–d. These displacements are the result of elevated temperature and changes in the coal layer's mechanical parameters due to gasification (table 1). Maximum displacements caused by elevated temperature amount to +0.14 m and occur near the roof of coal at time t = 1 month from the start of the gasification process, then subsidence starts due to the reduction of coal mechanical parameters.
Over time, small displacement values also appear at a greater distance from the gasified deposit, also at the ground level. Fig. 9a–d shows the vertical displacements in the cross-section over the course of the gasification process, initially positive as a result of thermal expansion, then negative due to subsidence after the end of the process, caused by the reduction of mechanical parameters. Final displacements of the ground level amounted to −5 cm in all versions (A, B and C). The cross-sections of displacements for model with modified constitutive equations (A) and constitutive equations of water (B) at selected instants above the deposit (y = 750, z = 58) is shown in fig 10a. The spatial distribution of displacements (A) on the background of the temperature field are shown in fig. 10b and 10c. Time instants t = 1 month, and the final instant of modelling, were selected. In figure 10c, besides displacements caused by thermal expansion, oppositely oriented displacements can be seen in the first section of the gasified deposit resulting from changes to the medium's mechanical parameters. On the cross-sections can be seen difference between behaviour of medium with modified constitutive equations (A – gas is present in the pores during the process of gasification) and with constitutive equations of water (B) – medium filled with gas deforms faster, but the final displacements are the same for both cases (−1.4 m). Subsidence evolves faster in the case A – 4 cm after 1.5 year (2 cm for version B).
Horizontal displacements are small compared to vertical displacements. Similar to the vertical displacements, the greatest values of displacements caused by temperature are reached in the layer of clays above the coal deposit. The distribution of horizontal displacement vectors in cross-section z = 51.5 at time t = 1.5 year (greatest displacements) is shown in fig. 11.
In the gas-filled space, conductivity is lower and seepage rate is higher due to its low viscosity with respect to water. According to Sutherland's formula, the viscosity of gas increases as temperature rises, while conductivity decreases.
Water seepage rate vectors at t = 1.5 year are shown in fig. 13 – their direction indicates that water reaches the gasification area according to the assumptions of the reduced-pressure method of coal gasification, providing the water vapour required for the process and also preventing migration of pollution. In case of no change in parameters (version C), seepage rates are smaller in the volume, where phase transition occurred, but change of rates outside this region is negligible.
The distribution of heat transfer coefficients at time t = 1 year is presented in figure 14. It shows a difference between the hot, gas-filled area and the cooler, liquid-filled area, which is due to lower heat transfer through gas and change in thermal conductivity along with temperature (the cross-section was made accounting for the shapes of finite elements). The difference in thermal conductivity is visible only in the closest proximity to the coal gasification generator, and its influence is relatively small, which coincides with the conclusions of Otto and Kempka .
However, in unsteady conditions, one must take into account also the specific heat. After the heat source is turned off, in case of water vapour in the pores, temperature drops relatively fast, while liquid water stores heat for many years until it is dissipated in the rock mass. The amount of heat needed to maintain temperature of the process during the process is also smaller after phase transition. Figure 15 shows temperature field in 6th year of simulation and the course of temperature in point 71,575,063 (10 m above the generator) in variant A (water vapour in pores) an variant C (no distinction between parameters of gas and liquid). The maximum temperature in variant C is 10°C higher due to thermal inertia and lack of insulating layer with gas in pores.
From the perspective of environmental impact, the most important issues are seepage (and migration of pollution) and subsidence. Final subsidence amounted to 5 cm in all three versions. Thanks to the application of reduced pressure with respect to the pressure in the rock mass, water flows into the generator, supplying it with water vapour and preventing propagation of pollution. The model, as presented, makes it possible to test different process parameters (such as temperature and pressure in the generator) for the purpose of determining the optimal conditions for the specific geology. Depending on the need, it is possible to use the model to analyse other issues that are not the subject of this article, such as plasticity potential in the rock mass and viscoelastic strain.
UCG is a problem that requires modelling of coupled fields of temperature, seepage and displacements. The proposed thermo-consolidation model enables modelling of such fields while simultaneously accounting for the presence of fluid in different states of matter within the medium's pores. Due to the different mechanical properties of both fluids – the linear stress/strain characteristic of the liquid and the non-linear characteristics of gases, there was a need for accounting gas laws in constitutive relations of the poroelastic medium.
In the proposed constitutive equations of the medium, successive terms correspond to the influence on stresses in the medium and in the fluid exerted by: strains of the medium, strains of the fluid phase, temperature changes in the solid phase and temperature changes in the fluid phase (or in the case where the medium is considered homogeneous with mean thermal parameters, the influence of the change in the medium's temperature). Terms of equations defining the influence of temperature changes on stresses use material constants of the poroelastic medium. The non-linear characteristic of gas also has an effect on heat transfer equations, in which compressibility and thermal expandability of both phases are also present.
The differences between the thermoporoelastic model (B) and proposed version with gas laws (A) are highest in vicinity of the UCG generator, so usefulness of such detailed model depends on scale – for modelling phenomena near the generator, the proposed model will be more accurate. However, time of computations is longer (about 270% compared to version B), so in modelling, only the phenomena at greater distance (such subsidence) model of thermoporoelasticity would be sufficient. In the UCG model, models for water and water vapour in the medium's pores were applied depending on the temperature and pressure. The area in which water vapour is present in pores of the medium was determined based on the Clausius–Clapeyron equation and, depending on the state of matter, constitutive equations and material parameters appropriate to gas were applied. For the assumed geology and gasification scenario, displacements as well as directions and rates of seepage were determined, which are of great importance to the safety of the process when considering migration of pollution. Under the assumed conditions, final vertical displacements of the ground's surface amounted to −5 cm, as a result of the drop in coal's strength parameters. Seepage was directed towards the generator, preventing migration of pollution in the rock mass and providing the water vapour required in the process. This model also made it possible to determine stresses in the skeleton and fluid, and to track changes of temperature-dependent parameters such as density, viscosity, heat transfer coefficient and hydraulic conductivity at selected points. This model can be modified (different geometry and material parameters can be defined in a manner that does not require high time input), and if additional data become available, it can be expanded by further material analyses.
Material parameter values of individual geological layers.
Thermal parameters of individual pore fluids and medium.