## Abstract

The fractured carbonate rock reservoir is widespread in Sichuan Basin, and the characteristics of different areas are different. The development of natural fractures is varying degrees, lost circulations occur frequently, and the formation heterogeneity is strong, which causes the formation not sufficiently stimulated by acidizing. It may affect the effectiveness of reservoir stimulation. To advance the whole stimulation effect of the heterogeneous fractured carbonate reservoir, a new solution for determining the invasion radius during drilling in a fracture network reservoir is presented, which is based on solute transport and convection-diffusion equations. It can predict the invasion radius caused by mud loss and determines the range of mud loss invasion, which clarify the scope and degree of reservoir damage. The formation of skin factors polluted by mud loss was calculated. The experiments verified that the acidizing technology can remove reservoir damage and reduce the polluted formation of skin factors. The opening pressure of the nature fracture closed is calculated which can control the acidizing area. It is confirmed how many fractures in the carbonate reservoir can be opened under the wellhead pressure limit, which meets the construction conditions of acidizing fractured reservoirs. The framework of network-fracture deep acidizing technology was established, which can efficiently break through the detrimental zone caused by lost circulation, break down the natural fracture network, and decrease the formation of the skin. The restart pressure of natural fractures was calculated, and the design parameters such as pump pressure and displacement were optimized to quantify the scope of reservoir stimulation and the scale of acid fluid. The technique of network-fracture deep acidizing was applied for well A, the formation of skin after acidizing can be reduced to -4, and the testing production of well A was $58.87\xd7104\u2009m3/d$. The technique of network-fracture deep acidizing can quantify the acid scale and sweep area in acid fracturing design, which develops the fracturing efficiency and improves the fracturing engineering.

## 1. Introduction

Fractured reservoirs are widely distributed in Sichuan Basin. In recent years, exploration practice has confirmed that fractured reservoirs have great exploration and development potential [1–4]. The fractured reservoir has low porosity and permeability, strong heterogeneity, and locally developed natural fractures and vugs which cause drilling fluid leakage during drilling and completion; therefore, formation damage is serious. Given the development of natural fractures in this kind of reservoir, the phenomenon of lost circulation is common in the drilling process, and the reservoir is polluted and blocked to varying degrees, so the main goal should be to remove plugging acidification and dredge the natural fracture cave system [5–8].

Due to the development of natural fractures in the reservoir, drilling fluid leakage is common, and the drilling fluid invades the reservoir with the natural fracture system to form a large gel leakage damage zone; the conventional plugging removal acidification process is difficult to make the acid fluid break through the leakage pollution zone and cannot fully remove the reservoir pollution blockage [9–15]. Therefore, it is necessary to predict the depth of the damage zone formed by the invasion of drilling fluid into the reservoir, which can be achieved by measuring the changes in formation properties. Researchers have mainly studied the invasion model of drilling fluid filtrate to predict the invasion depth range and verified the model through field data fitting. At the same time, it can also be used to measure and monitor the properties of near-wellbore zones and reservoirs [16–18]. The drilling fluid intruded into the near-wellbore formation is displaced and mixed with the formation fluid, resulting in the formation of a certain depth of reservoir damage zone near the near-wellbore formation [19–22]. Civan proposed an improved migration equation of two-phase fluid in deformable porous media taking into account the diffusion, cake formation, and filtrate invasion of compressible/incompressible cake model with/without solid invasion [23–27]. Ramakrishnan and Wilkinson established a radial invasion model of water-based drilling fluid filtrate, which can determine the saturation of oil and filtrate in the formation water phase [28]. Pordel Shahri et al. simulated the convection-diffusion movement of drilling fluid filtrate in fractures with limited conductivity and the convection-diffusion movement facing the foundation block through the fracture wall [29]. The depth of damage zone formed by drilling fluid invading the reservoir is closely related to factors such as formation properties (porosity, permeability, gas saturation, and electricity), engineering parameters (drilling differential pressure and soaking time), and drilling fluid performance (density, water loss, and its system). For field applications, although too many parameters will make the prediction depth more accurate, the calculation is more complex, and the model calculation parameters may not be all available. Therefore, it is necessary to develop a new, fast, and accurate method to predict the depth of mud leakage invasion.

Another key factor is to determine how many closed natural fractures are fractured and treated by acidizing. It is also important to judge the fracture deformation behavior caused by acidizing. Hydromechanical coupling models were established based on the phase field method to investigate the effect of vugs on hydraulic fracture propagation, which were solved by finite element approximation and iterative method for the hydraulic fracturing problem [30–34]. Numerous kinds of literature have been carried out to analyze the influence of acid fracturing on formation rocks. Fu et al. presented an explicitly integrated, fully coupled discrete finite element approach for the simulation of hydraulic fracturing in arbitrary fracture networks, and their model was verified against the K–G–D closed-form solution for the propagation of a single hydraulic fracture and validated against laboratory testing results on the interaction between a propagating hydraulic fracture and an existing fracture [35]. Ghommem et al. validated a 3D core scale predictive model for carbonate acidizing and showed a good capability to qualitatively capture the experimentally observed wormhole morphology and quantitatively predict the pore volume to break through [36]. Ma et al. developed a three-dimensional thermal-hydro-chemical (T-H-C) coupled unified pipe-network method (UPM) to simulate the process of acidizing treatment in fractured carbonate rock, which considered thermal effects and incorporated multiphysical governing equations for explicitly represented fractures in the current method [37]. Guo et al. further carried out a series of hydraulic fracturing experiments and found that etching wormholes induced by the rock-acid reaction would change the pore structure of carbonate rocks, and gelled acid enables to increase both the width and the roughness of main fracture surfaces [38]. Gou et al. conducted true triaxle fracturing experiments on carbonate rocks and observed that viscoelastic surfactant acid made the hydraulic fracture more conducive to connecting the nature fractures [39]. Jin et al. established and verified a damage constitutive model considering the compression hardening process of the acid-corroded sandstone under uniaxial loading, and proposed a theoretical model for evaluating the brittleness index (BI) of the acid-corroded sandstone based on energy evolution theory and damage constitutive relation [36]. Cao et al. proposed a method to dissect the internal structure of fault-karst reservoirs that is helpful for $y$ targeted acid fracturing, which improved the prediction accuracy and the recovery efficiency of fault-karst reservoirs [41]. Xing et al. studied the uniaxial compressive strength (UCS) and uniaxial tensile strength (UTS) by five types of rock samples and obtained the fracture extension degree of each lithology via acid dissolution experiments [42].

By establishing the leakage and migration model of drilling fluid in the fracture network, the leakage and invasion depth of drilling fluid is quantified, the range of reservoir damage zone caused by well leakage is defined, and the natural fracture opening pressure is calculated. Based on not forming new fracturing fractures, the natural fracture system is fully blocked and dredged, the reservoir skin coefficient is reduced, and the network-fracture depth acidification technology based on leakage and invasion depth and natural fracture opening pressure is formed. Through field application, the scope of reservoir leakage damage zone of fractured formation in well A is predicted, and the opening pressure of natural fractures in the reservoir is calculated, to optimize the acidizing design parameters of the well, improve the acidizing efficiency, and achieve good acidizing transformation effect on site.

## 2. Characteristics of the Fractured Reservoir

### 2.1. Reservoir Physical Properties

Taking the fractured reservoir in Sichuan China, Western Sichuan as an example, its average depth is 6700 m, the average temperature is 150°C, the average formation pressure is 125 MPa, the average porosity is 0.75%, and the average permeability is 0.95 mD (Figure 1). The observation of the core and microthin section reveals that the reservoir space of the fractured formation is mostly fractures and corrosion holes associated with fractures. The relationship between core porosity and permeability is poor, and most samples have the characteristics of low porosity and relatively high permeability, indicating that fractures play a very important role in improving reservoir porosity and permeability (Figure 2).

### 2.2. Rock Mechanical Properties

Triaxial compression tests were carried out on rocks of the Maokou formation in Western Sichuan, Sichuan Basin. The experimental results show that the average compressive strength of the Maokou group rocks is 445.65 MPa and the average Young’s modulus is $3.60\xd7104\u2009MPa$, with an average Poisson’s ratio of 0.23 (Table 1). It can be seen from Figure 2 that there are visible natural fractures in the Maokou formation rock sample, which will inevitably lead to the reduction of the compressive strength of the rock sample.

The results of the differential strain test show that the maximum horizontal principal stress of the Maokou group rocks is 163.51 MPa, the minimum horizontal principal stress is 141.21 MPa, and the vertical principal stress is 190.66 MPa (Table 2). The horizontal two-dimensional stress difference is 22.3 MPa, and the horizontal stress difference coefficient is 0.16. Combined with the three-dimensional principal stress distribution law, it is considered that the formation of vertical fractures is the main form of fracturing in the Maokou formation. Combined with the paleomagnetic experiment and differential strain experiment (Table 2), it is concluded that the direction of horizontal maximum principal stress is NE125°.

## 3. Optimization of Network-Fracture Deep Acidizing Technology

### 3.1. Technical Connotation

Based on the geological and engineering characteristics of the reservoir, firstly, predict the depth of drilling fluid leakage and invasion, combined with the prediction model of reservoir skin factor, clarify the scope and degree of reservoir damage, and quantify the scope of reconstruction and the scale of acid fluid. Calculate the restart pressure of natural fractures, and optimize the design parameters such as pump pressure and displacement. Finally, improve the efficiency of reservoir stimulation of fractured reservoirs. The technique design process is shown in Figure 3.

### 3.2. Prediction of Drilling Fluid Leakage Damage Zone

For the fractured reservoir with natural fractures, the drilling fluid leakage is large and the reservoir damage is serious (Figure 4). The goal of reservoir stimulation is to remove the reservoir pollution and blockage and dredge the natural fracture cave system. Therefore, the key to reservoir stimulation is to break through the damage zone formed by the leakage of drilling fluid along the natural fracture system, unblock the natural fracture system, and reduce the reservoir skin coefficient to the greatest extent [43, 44].

#### 3.2.1. Mathematical Model Based on Drilling Fluid Solute Transport and Convection-Diffusion

Due to the leakage and invasion of drilling fluid, drilling fluid with different concentrations is distributed in the natural fractures of the formation. The fluid at the open end of the borehole wall fracture is pure drilling fluid, which is regarded as it is a dimensionless concentration of 1. The dimensionless concentration of drilling fluid is consumed with the radial extension of the fracture in the process of fracture leakage and flow and percolates and diffuses on the fracture wall. At the same time, it is also diluted by the formation fluid, which gradually reduces its concentration, until the concentration at the maximum invasion depth of drilling fluid is close to 0; as shown in Figure 5, the yellow curve represents the dimensionless concentration change curve of drilling fluid in the fracture. Therefore, the drilling fluid is regarded as a pure fluid different from the formation fluid. Based on the solute convection-diffusion theory, the dimensionless concentration distribution of drilling fluid leakage flow in the fracture network is simulated and studied, and the dimensionless concentration distribution profile of drilling fluid in the fracture is used to calculate the maximum depth of drilling fluid leakage invasion. Based on the fracture network geological model and the description of drilling fluid leakage and invasion, the leakage and migration model of drilling fluid in the fracture network is established.

^{3}), $Df$ is the diffusion coefficient of drilling fluid concentration in radial fracture (L

^{2}/T), $r$ is the penetration depth of drilling fluid leakage (L), $w$ is the average half-width of fracture (L), $q$ is drilling fluid leakage rate (L

^{3}/T), $q\u2019$ is drilling fluid leakage rate drilling fluid convection exchange flux between fracture and base block (L

^{3}/T), $v$ is steady flow velocity of drilling fluid in radial fractures (L/T), $re$ is the maximum penetration depth of drilling fluid leakage (L), and $t$ is leakage time (T).

^{3}/T), $rw$ is the $e$ radius of the wellbore (L), $r$ is the radial invasion radius of drilling fluid (L), and $h$ is the effective thickness of leakage layer (L).

^{2}/T) and $D\u2217$ is the molecular diffusion coefficient (L

^{2}/T).

^{2}/T), $nm$ is the volume fraction of the rock base block in the fracture network, $\varphi m$ is matrix porosity of rock base block, $\tau $ is the tortuosity of natural fracture, and $\u2202Cf/\u2202zi$ is the concentration gradient of drilling fluid on the contact surface between fracture and base block.

^{2}).

^{3}).

#### 3.2.2. Model Boundary Condition

*(1) Model Initial Condition*. When the drilling fluid does not leak into the fracture network, the dimensionless concentration of drilling fluid in the fracture network is

*(2) Model Boundary Condition*. The dimensionless concentration of drilling fluid at the borehole wall is

#### 3.2.3. Model Calculation

Formula (16) is to solve the dimensionless concentration distribution of drilling fluid with the invasion depth of drilling fluid leakage in the process of drilling fluid leakage.

Using mathematical software MATLAB, we can calculate and solve the matrix equations of drilling fluid leakage invasion depth distribution, and then, the distribution of drilling fluid leakage and invasion depth in natural fractures can be determined.

### 3.3. Calculation of Restart Pressure of Natural Fracture

#### 3.3.1. Mechanical Model of Fractured Formation with Natural Fracture

Firstly, two cases of natural fracture opening are discussed, and the model shown in Figure 6 needs to be established.

The tensile fracture of natural fractures needs to meet $\Delta P$$Pnet\u2265\Delta \sigma $.

It can be seen that the shear failure of natural fractures needs to meet $\Delta P$$Pnet\u2265\Delta \sigma +\sigma 0/\mu f$.

#### 3.3.2. Mechanical Model of Fractured Formation without Natural Fracture

For reservoirs without natural fractures, the above two criteria are not applicable, and new models and mechanical criteria need to be established. The mechanical model of fracture formation in a rock mass with undeveloped fractures is shown in Figure 7.

Elliptical fractures are formed in the plane wireless domain, the long half axis is $lf$, the short half axis is $w$, axial compressive stress at infinity is $\sigma H$, and compressive stress in the short axis direction is $\sigma h$, and there is uniform pressure in the seam $Pnet$.

The boundary condition of the model:

When $y=0,x\u2264lf,$$\sigma y=\u2212Pnet,\sigma xy=0,$ and $x2+y2\u27f6\u221e,$ we can take as $\sigma x\u27f6\sigma H,\sigma y\u27f6\sigma h,\sigma xy\u27f60.$

To break the rock body and produce new joints, it is necessary to meet $Pnet\u2009max\u2265\sigma H\u2014\sigma h+St$; that is, temporary plugging of natural fractures does not develop, and the reservoir forms new branch fractures in the rock mass need to meet $\Delta PPnet\u2265Sum$ of maximum and minimum horizontal stress difference and rock tensile strength.

#### 3.3.3. Critical Conditions for Closed Natural Fracture Opening

The stress of natural fractures depends on the total stress state of the formation and the orientation and dip angle of the fractures. Even if the natural fractures in the stratum are at the same depth, if they develop in different directions and dip angles, their normal effective stress is also different. For horizontal joints and vertical joints, they are only affected by vertical effective stress and horizontal effective stress, respectively. It is assumed that the principal stresses on the crack surface along the $X$, $Y$, and $Z$ axes are $\sigma x$, $\sigma y$, and $\sigma z$.

The opening of natural fractures depends on the effective stress of natural fractures, which is not constant. When the fluid pressure in the natural fracture increases, the effective stress acting on the natural fracture surface decreases, and the fracture opening will increase. As the crack opening changes, the flow capacity of the fracture also changes, showing stress sensitivity.

To judge whether the closed natural fracture can open under a certain acid injection displacement, a simplified physical model is established as follows: inject liquid into the wellbore at a certain flow rate. When the injection volume is greater than the outflow volume, the wellbore internal pressure will increase due to fluid compressibility. When the bottom hole pressure exceeds the critical pressure of the natural fracture opening, it is judged that the natural fracture can open and meet the conditions of network-fracture acidizing construction.

^{3}), $Vinj$ is the volume of injected fluid per unit time (m

^{3}) $Vout$ is the volume of fluid flowing out per unit time (m

^{3}), $Pi$ is the initial bottom hole pressure (MPa), and $\Delta P$ is the added extra pressure (MPa).

*V*is a certain value, according to the definition of fluid compressibility, $\Delta P$ can be expressed as

^{-1}).

### 3.4. Simulation and Calculation of Reservoir Skin Factor

^{3}), and $\varphi $ is the average porosity of the reservoir (%).

### 3.5. Experimental Verification and Evaluation

The reservoir core is used for artificial fracture. During the acid displacement experiment, by increasing the confining pressure, it is ensured that the cut fracture remains closed during the displacement experiment. The experimental results show that under the condition of not opening the crack, the acid liquid does not penetrate the core, and only intense dissolution occurs locally at the end face of the core (Figure 7). On the other hand, by appropriately reducing the confining pressure, it is ensured that the split fracture remains open during the displacement experiment. The experimental results show that when the fracture is opened, the acid can easily break through the pollution zone and form a certain depth of etching grooves on the fracture wall (Figure 8). The experimental simulation results also show that for the reservoir with closed natural fractures, the acid solution mainly dissolves along the acid injection end face, and no wormhole is formed. For reservoirs with open natural fractures, acid first flows along the natural fractures and erodes the wall of the natural fractures (Figure 9), and wormholes expand at the tip of the natural fractures (Figure 10).

## 4. Application and Discussion

### 4.1. Prediction of Leakage Intrusion Depth in Well A

#### 4.1.1. Model Parameter Assignment

The well section is 3330.4-3372.5 m. The gas reservoir in Section 1 is comprehensively interpreted, and the reservoir thickness is 35.9 m. The reservoir porosity is 1.1% ~6.9%, the average porosity is 2.5%, and the average permeability is 3.62 mD. Type I reservoir is 1.3 m, with an average porosity of 6.7%. Class II reservoir is 17.4 m, with an average porosity of 3.4%. Class III reservoir is 17.3 m, with an average porosity of 1.9%. The average number of fractures is 20/m, and the fracture width range is 0.17-0.41 mm. The basic parameters required by the model are shown in Table 3.

#### 4.1.2. Numerical Simulation Results

There was one lost circulation and one gas invasion in the drilling of the oil test section (3330.0~3358.0 m), with a cumulative loss of 629.4 m^{3} and a fracture width of 1.7~4.1 mm according to logging interpretation. The reservoir geological parameters and field leakage data are substituted into the prediction model, and the prediction value of drilling fluid leakage invasion depth is obtained through a finite difference solution. According to the prediction results of numerical simulation (Figure 11), the invasion depth of drilling fluid increases with the lost circulation time. When the lost circulation occurs for 50 hours and the loss is 629 m^{3}, the radius of the lost circulation damage zone of the drilling fluid is about 40 m.

### 4.2. Calculation of Natural Fracture Opening Net Pressure

According to the rock mechanics experiment in this well, the calculated cohesion is 24.1~92.4 mPa, with an average value of 69.3 mPa. The natural fractures with approach angle ≤60~80 are mainly tensile (Figure 12). According to the model calculation, the natural fracture opening PNET is 4.9~12.7 mPa (Figure 13), which provides a basis for the pump pressure in the stimulation construction.

### 4.3. Construction Pressure Calculation of Natural Fracture Opening

According to the formulas (27)–(33), the variation law of bottom hole pressure under different displacement conditions is calculated.

Calculate the construction friction, and calculate the pipe string concerning the friction: $tubing\u2009\Phi 88.9\u2009mm\u2009pipe\u2009wall\u2009thickness\u200912.09\u2009mm\xd71800\u2009m+tubing\u2009\Phi 88.9\u2009mm\u2009pipe\u2009wall\u2009thickness\u20099.53\u2009mm\xd74670\u2009m+tubing\u2009\Phi 73mm\u2009pipe\u2009wall\u2009thickness\u20095.51\u2009mm\xd7383\u2009m$.

*/*m

^{3}), $q$ is the pump displacement (m

^{3}/min),

*L*is the well depth (m), and $D$ is the tube diameter (m).

According to formula (39), take three points in the stable section of the curve to calculate the average friction factor $F$, and then, bring the calculated friction factor into formula (39) to calculate the friction and friction coefficient under different displacements. The calculation results are shown in Figure 14.

When the acid injection displacement is greater than 4.0m^{3}/min, the conditions for natural fracture opening are met, and the displacement at this time is the minimum displacement required for natural fracture opening in formation. When the acid injection displacement reaches 5.0m^{3}/min, the wellhead pressure is close to the pressure limit (125 MPa), and the displacement can no longer be increased. At this time, the bottom hole pressure is 160 MPa, and the displacement at this time is the maximum displacement that can be achieved during construction (Figure 15). Under deep conditions, the minimum fracture opening pressure is 143 MPa (0° inclination), 149.8 MPa (30° inclination), 155.4 MPa (60° inclination), and 165 MPa (90° inclination). Therefore, under this condition, 0°, 30°, and 60° natural fractures can open, while 90° natural fractures cannot open. It is confirmed that most of the fractures in the reservoir of this well can be opened under the condition of a wellhead pressure limit, which meets the construction conditions of acidizing fractured reservoirs.

### 4.4. Acidizing Design and Construction Optimization Analysis

#### 4.4.1. Optimization of Acidizing Construction Parameters

According to the results of the prediction model of drilling fluid leakage invasion depth established above, it is concluded that the invasion depth of drilling fluid is approximately 30 m, so the acid corrosion fracture and acid fluid sweep range should be greater than 30 m in the acidizing design, FracPT software is used to simulate that when the scale of acid fluid is 320m^{3}, and the corresponding length of acid corrosion fracture is 35 m, which can remove the leakage pollution of drilling fluid, as shown in Figures 16 and 17. Through the simulation of the reservoir skin coefficient, it is obtained that the reservoir skin coefficient after 320m^{3} gel acid acidification is -4 (reservoir depth is about 20 cm), which is shown in Figure 18. Combined with the prediction model of drilling fluid leakage and invasion depth, it shows that reservoir acidification in this scale can reduce the skin coefficient of the fractured reservoir.

#### 4.4.2. Analysis of Acidizing Construction Effect

When the oil pressure of the acid solution in the well is blocked, the oil pressure of the acid solution in the well is gradually reduced, which means that the oil pressure of the acid solution in the well is blocked, and the oil leakage of the reservoir is gradually reduced in the acid injection stage, which means that the acid solution in the well is blocked, and the oil pressure of the reservoir is gradually reduced in the acid injection stage, as shown in Figure 19. The phenomenon of oil pressure drop in acidizing construction curve is analyzed, which shows that the acid solution relieves the reservoir pollution and blockage caused by lost drilling fluid and dredges the natural fracture cave system. After acidizing, the well obtained a high test production of $58.87\xd7104\u2009m3/d$. The field application proves that the prediction model of drilling fluid leakage and invasion can optimize the acidizing design of reservoir stimulation and improve the efficiency of reservoir stimulation.

## 5. Conclusions

For the fractured reservoir with serious leakage pollution, first, remove the pollution blockage near the well, reduce the reservoir skin coefficient, break through the leakage blockage zone, realize deep acidification, expand the transformation scope to increase the discharge area, and improve the efficiency of conventional plugging removal and acidification transformation.

Based on the drilling fluid leakage invasion model of the upper fracture network, the scope of the leakage damage zone is defined. Combined with the calculation of the restart pressure of natural fractures, the opening pressure of nature fracture closed is calculated which can control the acidizing area. It is confirmed how many fractures in the carbonate reservoir can be opened under the wellhead pressure limit, which meets the construction conditions of acidizing fractured reservoirs.

The formation of skin factors polluted by mud loss was calculated. The experiments verified that the acidizing technology can remove reservoir damage and reduce the polluted formation of skin factors.

The network-fracture depth acidification design parameters are further optimized, which provides a theoretical basis for field application.

The framework of network-fracture deep acidizing technology was established, which can efficiently break through the detrimental zone caused by lost circulation, break down the natural fracture network, and decrease the formation of the skin. The restart pressure of natural fractures was calculated, and the design parameters such as pump pressure and displacement were optimized to quantify the scope of reservoir stimulation and the scale of acid fluid. It has been applied to gas wells and achieved a remarkable acidizing effect.

## Data Availability

Some or all data, models, and code generated or used during the study appear in the submitted article. All data, models, or code generated or used during the study are available from the corresponding author by request.

## Conflicts of Interest

The authors declare no conflict of interest.

## Authors’ Contributions

Mingwei Wang was responsible for the conceptualization, funding acquisition, project administration, resources, and software. Song Li was responsible for the data curation, methodology, formal analysis, methodology, and writing of the original draft and wrote, reviewed, and edited the manuscript. Wen Wu was responsible for the investigation, methodology, software, and visualization. Tao Li wrote, reviewed, and edited the manuscript. Gensheng Ni was responsible for the visualization and supervision. Yu Fu and Wen Zhou were responsible for the funding acquisition and wrote, reviewed, and edited the manuscript.

## Acknowledgments

This work was supported by the National Natural Science Foundation of China: Study on Dynamic Characteristics of Methane/Carbon Dioxide in Shale Heterogeneous Reservoir under Multi-Field Coupling (41772150).