## 4 Development of a Mathematical Model of Pressure Distribution in a Structurally Complex Reservoir  
Prior to the development of any gas-bearing horizon, it is generally assumed that the pressure in the entire reservoir is conditionally equal. However, since the beginning of exploitation, pressure distribution has been non-uniform, and this process is influenced by many factors: development mode, production volume, properties of individual reservoir sections, heterogeneity of the structure, etc. All this suggests that the pressure in different points of the reservoir is formed as a single field distributed over space. Proceeding from the fact that the reservoir is an object with a heterogeneous structure, its description is built using a layered model, in which each individual section (layer) may have its own properties. Let us represent the gas-bearing horizon under study as a parallelepiped (Figure 2). Let us denote the main parameters of the modelling object. It is divided into 27 blocks in the three-dimensional domain. Geometric dimensions of the modelled system:  
\[x_{L}=210\ \text{m},y_{L}=210\ \text{m},z_{L}=120\ \text{m}.\]  
Each unit of the system has the following geometrical dimensions:  
\[x_{B}=70\ \text{m},y_{B}=70\ \text{m},z_{B}=30\ \text{m}.\]  
Figure 2: Layered model of a gas-bearing reservoir (\(x_{L}x,y_{L}y\) and \(z_{L}z\) are the geometric dimensions of the modelled system; 1–9 are the numbers of blocks). In the framework of the study, the system is considered to be non-deformable, thus the porosity of the blocks does not change with time. The system is also considered to be isotropic. When forming the mathematical model, the equation of pressure distribution in a gas-bearing deposit represented in three-dimensional pro-space (the equation of academician [PERSON]) was taken as a basis for the formation of the mathematical model [55]:  
\[\frac{\partial(\rho m)}{dt}=\frac{k}{\mu}\bigg{(}\frac{\partial^{2}p}{ \partial x^{2}}+\frac{\partial^{2}p}{\partial y^{2}}+\frac{\partial^{2}p}{ \partial z^{2}}\bigg{)}, \tag{1}\]  
where \(m\) is the porosity of the gas-bearing formation; \(\mu\)--absolute gas viscosity; \(\rho\)--gas density; \(p=p(x,y,z,t)\)--gas pressure at time \(t\); \(x,y,z\)--coordinates; \(x_{L},y_{L},z_{L}\)--final coordinates along the axes \(x,y,z;0<x<x_{L},0<y<y_{L},0<z<z_{L}\); and \(k\)--permeability of the gas-bearing formation. The Leibenzon function for a perfect gas is defined by the formula:  
\[p=\frac{\rho_{am}p^{2}}{2p_{am}}+C, \tag{2}\]  
where \(\rho_{am}\)--gas density at atmospheric pressure and \(p_{atm}\)--atmospheric pressure. Then, differentiating Equation (2) by coordinates two times, taking into account the equation of state of an ideal gas, we obtain for the left part of (1):  
\[\frac{\partial(\rho m)}{dt}=\frac{\rho_{am}m_{0}}{p_{am}}\frac{\partial p}{ dt},\]  
where \(m_{0}=m\)--is the porosity of the gas-bearing reservoir at \(p=p_{0}\). Considering this and the differentiated expression (2), substituting into (1), we obtain:  
\[\frac{\partial p}{dt}=\frac{k}{2\times\mu\times m}\bigg{(}\frac{\partial^{2} p^{2}}{\partial x^{2}}+\frac{\partial^{2}p^{2}}{\partial y^{2}}+\frac{ \partial^{2}p^{2}}{\partial z^{2}}\bigg{)}. \tag{3}\]  
The Leibenzon equation [56] can be written in another form by multiplying the right and left parts by \(p\) and substituting:  
\[p\frac{\partial p}{dt}=\frac{1}{2}\frac{\partial p^{2}}{dt}.\]  
Transform Equation (3) to the following form:  
\[\frac{\partial p^{2}}{dt}=\frac{k\times p_{0}}{\mu\times m}\bigg{(}\frac{ \partial^{2}p^{2}}{\partial x^{2}}+\frac{\partial^{2}p^{2}}{\partial y^{2}}+ \frac{\partial^{2}p^{2}}{\partial z^{2}}\bigg{)}, \tag{4}\]  
where \(p_{0}\) is the initial pressure of the gas-bearing reservoir.