Poroelastic Theory Research Articles

POROELASTICITY AND ROCK PHYSICS TEMPLATES

Two common rock physics templates (RPTs) that are used to identify geological facies and fluid trends derived from well log or pre-stack inverted seismic data involve cross-plotting VP/VS ratio against acoustic impedance, IP (the product of P-wave velocity, VP, and density) and mu-rho against lambda-rho, where mu and lambda are the Lam矣oefficients extracted from P and S-wave velocity, called the LMR, or LambdaMuRho, method. Using well log examples from an Alberta gas well and a deeper North Sea oil well, we show how to superimpose constant mu-rho and lambda-rho curves on a VP/VS versus IP cross-plot, which produce orthogonal trends of stiffness or porosity (from mu-rho) and fluid saturation (from lambda-rho) that match the data trends reasonably well. We then consider the related technique of KMR, or KMuRho, where K represents bulk modulus, and show how to superimpose constant K-rho trends on a VP/VS versus IP cross-plot. Again, we get orthogonal trends of stiffness or porosity (from mu-rho) and fluid saturation (from K-rho), but the fit to our real data examples is less accurate than with LMR. Using Biot-Gassmann poroelasticity theory, we generalize the LMR and KMR approaches into a technique we call FMR (Fluid-Mu-Rho). In the FMR techniques we introduce two new fitting parameters, f and d, where f is a fluid term derived from Biot-Gassmann theory, and d is the dry rock VP/VS ratio squared. When we superimpose fluid trends from the FMR technique on our two well log datasets, and adjust the f and d parameters, we achieve accurate fits to the datasets. Finally, we apply the FMR technique to a seismic case study from the Gulf of Mexico. In an appendix, we compare the FMR technique to the empirical Curved Pseudo-Elastic impedance, or CPEI, and Pseudo-Elastic Impedance for Lithology, or PEIL, methods.

The impact of boundary conditions on the acoustical behavior of lightweight granular particle stacks

The acoustical behavior of air-saturated granular stacks can differ significantly from that of more traditional sound absorbing materials such as fibrous webs and foams when tested in a standing wave tube. The latter materials can often be modeled as equivalent fluids by using one-dimensional theory, with the further assumption that their behavior does not depend on the input signal level. In contrast, level dependence of the response of lightweight granular materials such as glass bubbles has been consistently observed in previous studies. At low input levels, the absorption coefficient of glass bubble stacks shows solid-like behavior with multiple peaks associated with radial modes of the solid phase of the edge-constrained stack: i.e., the response is two-dimensional. In contrast, at high input levels, glass bubble stacks show one-dimensional fluid-like behavior, with the quarter wave depth resonance dominating the response. It is hypothesized here that at high input levels the particles slip axially at the circumferential boundary, thus showing a one-dimensional response, while at low levels, the particles are constrained to have zero axial velocity, thus creating a two-dimensional response. It is shown here that a two-dimensional finite difference implementation of the poro-elastic Biot theory can accurately reproduce the observed level-dependent behavior.

Modeling pore fluid pressure and deformation in unsaturated soils under gravitational compaction and cyclic loading

Cyclic loading of fluid-bearing soils often induces excess pore water pressure, leading to accumulation of strain in the solid matrix. Additionally, variations in material density and volumetric fraction due to gravitational effects generate momentum forces, which subsequently cause deformation in the soil. In this paper, we present a mathematical framework of poroelasticity that generalizes the simultaneous action of gravitational body forces and cyclic loading to analyze solid deformation and fluid pressure dissipation in partially-saturated porous media. To account for the effect of gravity, the gradient of the vertical displacement needs to be treated as an additional dependent variable alongside pore air and water pressures.We numerically solve a mixed strain-controlled and pressure-prescribed boundary-value problem using a finite difference scheme as a representative example to quantify the impact of gravitational forces. Our results show that gravitational compaction affects both total settlement and excess pore water pressure, with its influence being particularly significant at lower water saturations, regardless of excitation frequency, soil depth, or type. The effect of gravity is more pronounced in soils with higher compressibility but appears insensitive to excitation frequency. Notably, the gravitational effect correlates linearly with the depth of the soil deposit and, therefore, should be well integrated into the consolidation theory of poroelasticity.

Constraints on fracture connectivity from amplitude variations with incidence angle and azimuth analysis

The analysis of seismic reflection amplitude variations with incidence angle and azimuth (AVAz) is a powerful tool for studying and characterizing fractured formations. There is evidence to suggest that wave-induced fluid pressure diffusion (FPD) within connected fractures can significantly affect the AVAz response. This, in turn, indicates that this seismic attribute may contain information on the connectivity degree of the fracture network, which represents a key hydraulic property of fractured formations. The assessment of the sensitivity of seismic reflectivity with respect to fracture connectivity accounting for FPD effects to date is limited to 2D models, which is mainly due to the complexity and the associated computational cost of modeling the effective properties of the corresponding generically anisotropic 3D media. To overcome these issues, we use a numerical upscaling procedure based on Biot’s theory of poroelasticity that allows us to obtain effective anisotropic representations of 3D fractured formations accounting for the FPD effects. Using these effective properties, we then perform a plane-wave analysis to compute the reflection coefficients at the top of the corresponding fractured formation. The results of our numerical analyses extend the previous 2D studies and demonstrate that seismic reflection coefficients are highly sensitive to the degree of connectivity of the fracture networks. The reflectivity responses of formations comprising connected and unconnected fractures indicate pronounced differences not only with respect to incidence angle but also with respect to azimuth. Quite importantly, our results also demonstrate that crossplots of commonly used AVAz coefficients permit us to identify regions with higher or lower fracture connectivity.

MPeat2D – a fully coupled mechanical–ecohydrological model of peatland development in two dimensions

Abstract. Higher-dimensional models of peatland development are required to analyse the influence of spatial heterogeneity and complex feedback mechanisms on peatland behaviour. However, the current models exclude the mechanical process that leads to uncertainties in simulating the spatial variability in the water table position, vegetation composition, and peat physical properties. Here, we propose MPeat2D, a peatland development model in two dimensions, which considers mechanical, ecological, and hydrological processes together with the essential feedback from spatial interactions. MPeat2D employs poroelasticity theory that couples fluid flow and solid deformation to model the influence of peat volume changes on peatland ecology and hydrology. To validate the poroelasticity formulation, the comparisons between numerical and analytical solutions of Mandel's problems for two-dimensional test cases are conducted. The application of MPeat2D is illustrated by simulating peatland growth over 5000 years above a flat and impermeable substrate with free-draining boundaries at the edges, using constant and variable climate. In both climatic scenarios, MPeat2D produces lateral variability in the water table depth, which results in the variation in the vegetation composition. Furthermore, the drop in the water table at the margin increases the compaction effect, leading to a higher value of bulk density and a lower value of active porosity and hydraulic conductivity. These spatial variations obtained from MPeat2D are consistent with the field observations, suggesting plausible outputs from the proposed model. By comparing the results of MPeat2D to a one-dimensional model and a two-dimensional model without the mechanical process, we argue that mechanical–ecohydrological feedbacks are important for analysing spatial heterogeneity, shape, carbon accumulation, and resilience of peatlands.

Poroelastic effects of chemical loading in consolidation of residual soils

An analytical chemo-hydromechanical model of residual soils was presented in this study. The modeling procedure utilized a unique phenomenological approach, which highlighted some of the macro-mechanical influences responsible for the behavior of residual soils. The macromechanical analysis yielded an extended poroelastic theory. The developed model was applied in solving a typical civil engineering problem of soil consolidation. For a chemical loading treatment, the problem was using requisite boundary conditions. Thereafter, the mathematical software Mathematica was used for solving the ensuing diffusion equations using physico-chemical properties typical of a residual soil sample. Even though the chemical loading was instantaneous, it took some time before its effect could be felt everywhere according to the poroelastic theory. From obtained results the evolution of the generated pore pressure showed increasing segments of the soil returning to their initial state after the passage of the pore pressure front. The vertical displacement showed a slight increase or elevation of the soil surface. The flux of water flow in the soil was initially positive before turning negative, which can be explained by initial successful penetration of the infiltrating liquid until it met increasingly tortuous paths as result of lower permeability. The volume (per unit area) of water leaving the soil after chemical loading demonstrated that a portion of the water flux associated with the infiltrating liquid ends up leaving the soil through the sides. The developed model was later compared with a similar porochemoelastic model found in literature with reasonable agreement.

Acoustic response of patchy-saturated porous media: Coupling Biot's poroelasticity equations for mono- and biphasic pore fluids.

Patchy saturation is a term used in the seismic prospecting literature to describe the state of a geological formation in which two immiscible pore fluids prevail in mesoscopic-scale clusters. If the pore fluids have contrasting compressibilities, wave-induced fluid pressure diffusion (FPD) processes may induce significant attenuation and velocity dispersion on seismic waves. Biot's monophasic poroelasticity theory is widely used to model the seismic response of rocks containing binary patches of two immiscible pore fluids. Even though effective fluid approximations may help to represent more realistic partially saturated patches using Biot's monophasic equations, the so inferred dissipation may not be representative of actual biphasic FPD phenomena. In this work, Biot's equations for mono- and biphasic fluids are combined to model FPD processes in porous media, comprising fully and partially saturated patches. An analytical solution for one-dimensional layered patchy-saturated media is presented which permits to explain some, as of yet enigmatic, experimentally observed characteristics such as increased seismic attenuation and stiffening effects occurring at low saturations. The results show that the existence of an additional diffusive wave mode within partially saturated patches may render conventional binary and effective fluid approaches incorrect and error prone.

Shear Bands Triggered by Solitary Porosity Waves in Deforming Fluid‐Saturated Porous Media

AbstractThe interplay between compaction‐driven fluid flow and plastic yielding within porous media is investigated through numerical modeling. We establish a framework for understanding the dynamics of fluid flow in deforming porous materials that corresponds to the equations describing solitary porosity wave propagation. A concise derivation of the coupled fluid flow and poro‐viscoelastoplastic matrix behavior is presented, revealing a connection to Biot's equations of poroelasticity and Gassmann's theory in the elastic limit. Our findings demonstrate that fluid overpressure resulting from channelized fluid flow initiates the formation of new shear zones. Through three‐dimensional simulations, we observe that the newly formed shear zones exhibit a parabolic shape. Furthermore, plasticity exerts a significant influence on both the velocity of fluid flow and the shape of fluid channels. Importantly, our study highlights the potential of spontaneous channeling of porous fluids to trigger seismic events by activating both new and pre‐existing faults.

Vibration and aeroelastic instability analysis in GPL-porous multi-layered beam with the rotation effect

This research studies the vibrations and instability of a rotating three-layer beam under aerodynamic forces consisting of two functionally graded (FG) composite surfaces and a porous intermediate layer. Linear poroelasticity theory is used for modeling the porous layer while Young’s modulus and density vary along the thickness. Equivalent coefficients of the graphene nanoplatelets-reinforced composite (GPLs) are extracted by modified Halpin-Tsai theory, according to five different configurations. The equations of motion in the GPL-porous multi-layered beam are derived using the Euler–Bernoulli beam (EBB), Timoshenko beam (TB), and Reddy’s higher-order shear deformation theory (HSDT) theories, as well as the energy method and Hamilton’s principle. The most important results show the effect of reinforcements and their different configurations, porosity and its distribution, speed of rotation on natural frequency, critical flutter point, and loss factor for a GPL-porous multi-layered beam. The unique properties of GPL-reinforced porous materials facilitate the design of components capable of effectively managing dynamic loads and resisting aeroelastic instabilities. This results in more efficient, reliable, and durable systems across various industries.

Phase-field modelling of dynamic hydraulic fracturing in porous media using a strain-based crack width formulation

This paper presents a comprehensive study on phase-field modelling in COMSOL MultiPhysics for simulating dynamic hydraulic fracturing in porous media based on Biot’s poro-elasticity theory. The focus is on addressing the challenges associated with crack width estimation in this context. A new strain-based crack width formulation is proposed, offering improved accuracy in predicting fracture permeability and ease of implementation in numerical approaches. The model’s capabilities are extended to consider dynamic crack propagation by incorporating the kinetic energy in the governing coupled hydro-mechanical-damage equations. The numerical implementations in COMSOL MultiPhysics are thoroughly explained, providing insights into the techniques used to solve the governing equations. Verification examples, including the benchmark KGD verification, are presented to demonstrate the model’s capabilities in simulating hydraulic fractures in porous media and validate its accuracy and reliability. A final numerical example focusing on the dynamics of crack propagation in a gravity dam is simulated, allowing for a comprehensive examination of the model’s performance. The proposed strain-based crack width formulation and consideration of dynamic crack propagation contribute to improved accuracy in predicting fracture permeability.

Shear-wave velocity prediction of tight reservoirs based on poroelasticity theory: A comparative study of deep neural network and rock physics model

The absence of shear-wave (S-wave) velocity in well log data limits the effective seismic characterization of lithology, petrophysical properties, and fluid distribution of tight reservoirs in the Chang 7 shale oil formations of Ordos Basin, west China. However, conventional rock physics models (RPMs), e.g., the Xu-White model, struggle to accurately obtain the key parameters such as pore aspect ratio in the process of S-wave velocity prediction. This work performs a comparative study of S-wave velocity prediction based on deep neural network (DNN) and RPM for such reservoirs. By employing the poroelasticity equations which are capable of describing elastic wave propagation in the complex reservoirs saturated with viscous fluid, the relation between elastic and petrophysical properties is established with a DNN-based method. This approach quantifies the connection between P-/S-wave velocities and rock properties. On the other hand, we reformulate the traditional Xu-White model by incorporating the decoupled Kuster-Toksöz (K-T) model and the poroelasticity equations with viscoelastic fluid in the modeling process. The variation of pore aspect ratio with depth is also considered to characterize the complex pore structures, which is appropriate for improving the accuracy of prediction. Log data from the three wells are considered for verification, of which results show that the S-wave velocity prediction methods help improve the final result. Aided by the DNN-based method, the prediction exhibits a better performance if the training data is sufficient. In particular, provided a small amount of measured S-wave data (even ten data points), the trained DNN model can also provide reasonable predictions.

Modeling patchy saturation of fluids in nanoporous media probed by ultrasound and optics.

Nanoporous materials provide high surface area per unit mass and are capable of fluids adsorption. While the measurements of the overall amount of fluid adsorbed by a nanoporous sample are straightforward, probing the spatial distribution of fluids is nontrivial. We consider literature data on adsorption and desorption of fluids in nanoporous glasses reported along with the measurements of ultrasonic wave propagation. We analyze these using the so-called dynamic equivalent medium approach, which is based on Biot's theory of dynamic poroelasticity with coefficients that are continuous random functions of position. When further constrained by optical scattering data for similar systems, the calculations show that on adsorption the characteristic patch size is of the order of 100 pore diameters, while on desorption the patch size is comparable to the sample size. Our analysis suggests that one can employ ultrasound to probe the uniformity of the spatial distribution of fluids in nanoporous materials.

Estimation of Pore Structure for Heterogeneous Reservoirs Based on the Theory of Differential Poroelasticity

Estimation of Pore Structure for Heterogeneous Reservoirs Based on the Theory of Differential Poroelasticity

Coupled geomechanical analysis of irreversible compaction impact on CO2 storage in a depleted reservoir

The utilization of depleted gas reservoirs for carbon storage offers substantial advantages. However, it raises concerns related to the geomechanical effects of historic gas extraction on the CO2 sequestration operation. In this study, a coupled thermo-hydro-mechanical investigation is conducted on a multi-layered sedimentary system that includes sealed bounding faults. Our simulations exhibit various levels of reservoir compaction during gas extraction under different pre-consolidation conditions of the sediments, highlighting the pivotal role of reservoir compaction in geomechanical analysis. The results demonstrate the occurrence of strain-hardening compaction behavior in the reservoir during post-yield depletion. This compaction is accompanied by a significant reduction in porosity and permeability, as well as irreversible surface subsidence. Hysteresis in the stress state is induced by the aforementioned irreversible reservoir compaction through two primary mechanisms: poro-elastoplastic stressing and differential compaction on each side of the sealing fault. These mechanisms alter the magnitude and orientation of stress inside and outside the depleted reservoir. Moreover, caprock compaction is impeded and delayed by the irreversible reservoir compaction owing to poro-elastoplastic stressing. This implies that conventional methods relying on the poro-elasticity theory alone may overestimate pressure required to fracture the caprock by approximately 2 MPa. When considering the combined effect of poro-elastoplastic stressing and differential compaction, neglecting irreversible reservoir compaction may lead to underestimation of the critical pressure for inducing fault slip by up to 4.9 MPa. Additionally, regardless of whether plastic void compaction is considered, we recommend increased attention be focused on the sub-vertical faults to mitigate the risk of significant slip that potentially could also lead to upward CO2 leakage. In scenarios where a slip event occurs in a fault during reservoir depletion, the results show that a subsequent CO2 injection operation tends to stabilize the faults. Nevertheless, it is crucial to consider that fault stability could deteriorate rapidly over time, potentially leading to a second slip event during carbon sequestration.

Diffusiophoretic Fast Swelling of Chemically Responsive Hydrogels.

Acid-induced release of stored ions from polyacrylic acid hydrogels (with a free surface fully permeable to the ion and acid) was observed to increase the gel osmotic pressure that leads to rapid swelling faster than the characteristic solvent absorption rate of the gel. The subsequent equilibration of the diffusing ion concentration across the gel surface diminishes the osmotic pressure. Then, the swollen gel contracts, thereby completing one actuation cycle. We develop a continuum poroelastic theory that explains the experiments by introducing a "gel diffusiophoresis" mechanism: Steric repulsion between the gel polymers and released ions can induce a diffusio-osmotic solvent intake counteracted by the diffusiophoretic expansion of the gel network that ceases when the ion gradient vanishes. For applications ranging from drug delivery to soft robotics, engineering the gel diffusiophoresis may enable stimuli-responsive hydrogels with amplified strain rates and power output.

Semi-theoretical compression model of wet fiber web for the prediction of airflow rate in vacuum dewatering process of papermaking

This study expands the investigation of a mathematical method for predicting the airflow rate during the vacuum dewatering process of papermaking based on our previous research. By integrating poroelastic theory and initial saturation thickness, based on experimental data with a range of 20–60 kPa pressure difference and 20-200 g · m − 2 basis weight, a mathematical method with a semi-theoretical model for compression characteristics of the wet fiber web, the poroelastic porosity model (PEPM), was developed to support and expand on our previously established empirical model for porosity variations in wet fiber web, the variable porosity model (VPM). The accuracy of the mathematical method with the PEPM in predicting the airflow mass flux passing through the wet fiber web was validated using experimental data with a range of 0–60 kPa pressure difference and 20-200 g · m − 2 basis weight. The mean relative errors with the PEPM within the pressure difference range of 0-20 kPa and 0-60 kPa were 7.68% and 8.80%, respectively, compared to 17.44% and 11.67% for those with the VPM, demonstrating the improved applicability and generalization of the mathematical method, particularly in the zone of low pressure difference exceeding the range of the modeling experimental data. The PEPM effectively revealed the strain hardening behavior during the compression process and the diminishing returns of energy in the vacuum dewatering process with an increasing pressure difference. The developed mathematical method with the PEPM for the compression characteristics of the wet fiber web provides an effective and accurate way to predict the airflow rate during the vacuum dewatering process of papermaking.

How borehole wall elasticity affects wellbore storage capacity in periodic hydraulic tests

Characterizing fluid transport and pore pressure diffusion is key for monitoring and understanding many natural (e.g., seismically active zones and volcanic systems) and engineered environments (e.g., enhanced geothermal reservoirs and CO2 underground storage). Borehole hydraulic testing allows for inferring relevant properties of the probed subsurface volume, such as, for example, its transmissivity and diffusivity, assessing the governing flow regime as well as detecting the presence of hydraulic boundaries. Periodic hydraulic tests (PHT) achieve these objectives using an oscillatory fluid injection procedure while measuring the fluid pressure response in monitoring boreholes. The relevant information on the pressure diffusion process occurring in a probed formation is retrieved from the phase shifts and amplitude ratios between the flow rate and the interval pressure. In general, the interpretation of PHT data relies on the assumption that the pressure diffusion process is largely unaffected by the deformation of the probed rock volume. To assess this assumption, we present a conceptual PHT model based on the theory of poroelasticity to investigate the role played by hydromechanical coupling (HMC) effects during PHT and to assess whether and to what extent additional mechanical information can be extracted from these tests. These effects include the influence of the borehole wall deformation on the wellbore storage coefficient Sw, which quantifies the difference between the injected flow rates and those actually entering the porous formation. We show that, in homogeneous formations, the commonly taken approach based on the uncoupled diffusion solution for radial flow reproduces the results of our HMC model. However, neglecting the deformation effect of the borehole wall on the storage capacity of the system during injection can lead to underestimations of both transmissivity and diffusivity, particularly at shorter oscillation periods. We present analytical expressions for the borehole wall deformation contribution to Sw, which depends on the shear modulus and does not change with the oscillatory period. We show that it is possible to obtain the effective Sw along with the transmissivity and diffusivity using observations at various oscillatory periods. In the presence of hydraulic boundaries, such as no-flow or constant-pressure boundaries, at a given distance from the injection borehole, Sw becomes period-dependent, introducing apparent variations in the inferred shear modulus of the formation if the associated HMC effects are not accounted for. Our results thus highlight the importance of accounting for HMC effects in PHT data interpretation, not only to reliably and accurately characterize and monitor hydraulic properties but also to fully exploit the potential of PHT to provide additional information on the mechanical behavior of the probed rock volume.

Simulation of dynamic pulsing fracking in poroelastic media by a hydro-damage-mechanical coupled cohesive phase field model

This study develops a hydro-damage-mechanical fully coupled numerical method capable of modelling complex fluid-driven transient dynamic crack propagation in quasi-brittle poroelastic media. In this method, the fluid flow in both fractures and porous media is described by a fluid continuity equation with the modified Darcy-Poiseuille law based on the Biot's poroelastic theory. The fluid pressure and the inertial force of solids are coupled by governing equations of the mesh-insensitive phase-field regularized cohesive zone model that can simulate quasi-brittle multi-crack initiation, propagation, branching and merging without remeshing or crack tracking. The resultant displacement-pressure-damage coupled multiphysics system of equations is solved using an alternative minimization Newton-Raphson iterative algorithm with an implicit Newmark integration scheme within the finite element framework. The new method was first validated by a few 2D problems, including crack branching and deflecting in solids, fracking in a concrete cube, and consolidation and stress wave propagation in poroelastic media, subjected to various impulsive loadings. It was then applied to pressure pulsing fracking of a 40 m granite rock reservoir with extensive parametric studies of fluid viscosity and pulsing injection rate, mode and period. It was found that pulsing injection with higher fluid rates and lower fluid viscosities resulted in more developed crack patterns, and in particular, there existed an optimal pulsing injection period that could promote fracking under relatively low injection pressures. 3D horizontal well problems with multiple non-planar crack propagation and bifurcation were also successfully simulated to demonstrate the capacity and potential of the new method for engineering design and optimization of pressure pulsing fracking.

Structural relaxation in a Schiff base‐crosslinked aldehyde modified xanthan gum/gelatin hydrogel: Experimental and numerical validation using linear poroelasticity model

AbstractHydrogels hold immense promise in biomedical applications due to their biocompatibility, high water content, and versatile fabrication. This study focuses on the mechanical behavior of a novel polysaccharide/protein hybrid hydrogel, GEL‐AXG, synthesized via a Schiff base reaction between aldehyde‐modified xanthan gum (AXG) and gelatin (Gel). Hydrogel samples with varying AXG‐to‐Gel ratios were subjected to unconfined compression tests to assess their mechanical properties. The observed stress‐relaxation mechanism in deformed hydrogels primarily involves water migration. To quantify these mechanical properties, we applied the linear poroelasticity theory. Our results highlight that Gel‐AXG hydrogels with a 2:1 AXG‐to‐Gel ratio exhibit significantly higher peak and equilibrium stresses. This enhancement can be attributed to increased crosslink density and reduced dangling chain presence. Moreover, the linear poroelasticity formulation yielded a shear modulus of G = 44.91 ± 0.25 kPa for Gel‐AXG hydrogels with a 1:2 AXG‐to‐Gel volume ratio, which we identified as our optimized hydrogel.

Numerical study of multiple hydraulic fractures propagation in poroelastic media based on energy decomposition phase field methods

This paper proposes a novel phase-filed model for simulating hydraulic fracturing in poroelastic media under complex stress conditions. The main theoretical contribution lies in the fact that a generalized strain energy density decomposition is incorporated to hydraulic fracturing phase field model, which unifies hydrostatic-deviatoric and spectral decompositions and predicts material failure through classical Mohr-Coulomb failure criterion. Based on Biot’s poroelastic theory, an indicator function is used to distinguish between fluid flow in fractures and in pores. The fully coupled equations are solved by combining staggered schemes with Newton-Raphson methods. The proposed model is validated by comparing laboratory experimental results with numerical results. In order to demonstrate the influence of introducing hydrostatic-deviatoric and spectral energy decompositions, we simulated propagation of multiple hydraulic fractures and coalescence of hydraulic fractures with natural fractures. Finally, the effects of stress ratios (σx/σy) on multiple hydraulic fractures are analysed. It is found that when the stress ratio is 2, a more complex hydraulic fracture network is formed in the porous media.