Geofluids
Volume 2020 (2020), Article ID 8876153, 12 pages
https://doi.org/10.1155/2020/8876153

Research Article

A Digital Twin for Unconventional Reservoirs: A Multiscale Modeling and Algorithm to Investigate Complex Mechanisms

Tao Zhang,¹ Yiteng Li,¹ Jianchao Cai,2,3 Qingbang Meng,⁴ Shuyu Sun,1,2 and Chenguang Li⁵

¹Computational Transport Phenomena Laboratory, Division of Physical Science and Engineering, King Abdullah University of Science and Technology, Thuwal 23955-6900, Saudi Arabia
²Institute of Geophysics and Geomatics, China University of Geosciences, Wuhan 430074, China
³State Key Laboratory of Petroleum Resources and Prospecting, China University of Petroleum, Beijing 102249, China
⁴Key Laboratory of Tectonics and Petroleum Resources, Ministry of Education, China University of Geosciences, Wuhan 430074, China
⁵China National Oil and Gas Exploration and Development Company Limited, China

Correspondence should be addressed to Shuyu Sun; shuyu.sun@kaust.edu.sa

Received 9 July 2020; Revised 8 September 2020; Accepted 15 September 2020; Published 2 November 2020

Academic Editor: Jinze Xu

Copyright © 2020 Tao Zhang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

The special mechanisms underneath the flow and transport behaviors in unconventional reservoirs are still challenging an accurate and reliable production estimation. As an emerging approach in intelligent manufacturing, the concept of digital twin has attracted increasing attentions due to its capability of monitoring engineering processes based on modeling and simulation in digital space. The application potential is highly expected especially for problems with complex mechanisms and high data dimensions, because the utilized platform in the digital twin can be easily extended to cover more mechanisms and solve highly complicated problems with strong nonlinearity compared with experimental studies in physical space. In this paper, a digital twin is designed to numerically model the representative mechanisms that affect the production unconventional reservoirs, such as capillarity, dynamic sorption, and injection salinity, and it incorporates multiscale algorithms to simulate and illustrate the effect of these mechanisms on flow and transport phenomena. The preservation of physical laws among different scales is always the first priority, and simulation results are analyzed to verify the robustness of proposed multiscale algorithms.

1. Introduction

The successes in the commercial exploitation of unconventional resources, such as shale gas and tight oil, in North America have already changed the current world energy market, and the growing public concerns on the depletion of conventional oil and gas resources in the foreseeable future also stimulates more efforts in both academia and industry to investigate unconventional reservoirs [1, 2]. A large number of technical studies have been carried on for a better characterization of oil and gas storage in either or both adsorbed or free states and a clearer description of the rock properties including porosity and permeability [3–5]. On the other hand, environmental issues have challenged the development of the unconventional oil and gas resources. As a consequence of hydraulic fracturing and shale gas production, groundwater pollution has become a serious issue that haunts oil companies. Earthquake and gas explosion, the other two problems often criticized due to unconventional reservoir production, are increasingly arising public concerns. In order to achieve a better balance between recovery efficiency and environmental impacts, we need to pay more efforts for a thorough understanding of the special mechanisms that control the storage, flow, and transport of unconventional resources in subsurface reservoir in order to meet the growing global energy demands in an environmentally friendly and economical approach. For instance, the classical viewpoint that transient linear flow dominates the flow regime for multifractured horizontal wells is challenged by the anomalous diffusion behavior, which also enlightens us that our understanding of the complex mechanisms and flow behaviors in unconventional reservoirs still needs to be improved [6].

Numerical modeling and simulation have become a popular approach in the study of unconventional reservoirs, mainly due to their significant superiority on efficiency and flexibility [3]. In practice, experimental studies are limited to time and space scale and also restricted by field or laboratory conditions. However, these constraints can be avoided by well-designed numerical simulation. Moreover, numerical schemes can be further verified through experimental analysis, and some key parameters are tuned from experimental results. After continual improvement and optimization, such mathematical model is believed to be “realistic” and is reliable to be applied to solve practical problems in much larger time and space scale in order to guide or suggest engineering practice. As one of the cutting-edge techniques, digital twins have been extensively used in automatic production, predictive maintenance, and complete-cycle management. Digital twins, also known as DT, can be defined as a simulation process or simulation-based system that integrates multidiscipline, multiphysics, and multiscale numerical methods to make full use of physical models, sensor data, operation history, and other information and to complete mapping in virtual space so as to reflect the whole life cycle of a corresponding process or equipment in physical space [7]. The concept of digital twins was first proposed in the field of advanced manufacturing, along with the application of state-of-the-art information technologies in industrial processes [8, 9]. With the arrival of the big data era, the entire product life cycle produces plenty of data in the aspect of designing, manufacturing, marketing, and service. These data can be transferred into the digital models for simulation and analysis, and in turn, the numerical results can provide supports to improve practical manufacturing. The data obtained in physical reality is often fragmented and isolated, but well-designed numerical schemes in virtual space can describe the underneath physical or chemical correlations among the data. Thus, the intelligence and efficiency of industrial manufacturing are constantly improved to achieve intelligent production and management. The similar idea can be applied to reservoir simulation in which data of the reservoir geology, fluid properties, and environmental and operation conditions can all feed into the model in virtual space.

In this paper, we will combine several promising numerical models and algorithms that describe the special mechanisms behind the flow and transport behaviors in unconventional reservoirs in different scales as an exploratory investigation to construct the digital twin. A thermodynamically consistent flash calculation scheme is designed to consider the effect of capillary pressure on phase equilibria, and the dynamic sorption is included in the particle distribution function to establish a delicate Lattice Bhatnagar-Gross-Krook (LBGK) scheme to simulate the shale gas flow and transport using Lattice Boltzmann Method (LBM). Multicomponent ion exchange and double-layer expansion, both directly relevant to fluid salinity, are also modeled in different scales. Simulation results are presented to show the effectiveness of the proposed algorithms for the investigated mechanisms. The remainder of this paper is organized as follows. In Section 2, several important mechanisms in unconventional reservoirs are modeled, and the corresponding simulation algorithms and results are illustrated in Section 3. Conclusive remarks are provided in Section 4.

2. Design of the Digital Twin

A digital twin is expected to realize interaction and integration between physical space and virtual world due to its feasibility on real-time synchronization, real mapping, and high fidelity. In order to create the virtual model of physical entity in unconventional reservoirs, the fluid flow and transport process in realistic reservoir conditions are mathematically described in the digital space, while data is transferred into the twin and information is fed back from the twin. The procedures and processes in physical entity are expected to be enabled or expanded with new capabilities based on the feedback, which may be obtained via various approaches, such as data fusion analysis and decision iteration optimization. In this section, a digital twin for unconventional reservoirs is designed to bridge the physical space and digital world, with focus on several selected mechanisms which affect the complete cycle of evaluation and production. As shown in Figure 1, some processes are extracted first from the physical space, such as reserve prediction, recovery evaluation, hydraulic fracturing, and environmental effect. To provide more effective, real-time, and intelligent improvement and optimization schemes to these processes, the digital twin is designed with three parts: models on reservoir geometry, Darcy’s scale fluid flow, and pore scale fluid flow. Each part is decomposed further with representative techniques and research topics. It should be noted that phase equilibrium calculations are conventionally investigated in pore scale, but in [10], Darcy’s scale phase equilibrium was studied.

Figure 1: Design of a digital twin for unconventional reservoirs.

In this paper, we will start to model capillarity and dynamic sorption, both of which play significant roles in unconventional reservoirs and affect many engineering processes in practice. The widespread nanopores are often considered as the main cause of significant capillary effect and confinement effect. The strong interactions between gas molecules and rock surface result into the dynamic sorption, and the related concepts of Knudsen diffusion, Knudsen layer, and Knudsen number are introduced. The large amount of water needed for hydraulic fracturing is often questioned by the local communities, and seawater or produced water after proper treatment is expected as an alternative to save fresh water. Water salinity, one of the key factors affecting injection and fracturing performance, has also attracted numerous attentions. In this section, thermodynamic equilibrium modeling for multicomponent fluid mixtures in unconventional reservoirs is used to describe phase behaviors and properties. A delicate LBGK model is constructed to take into account the effect of dynamic sorption on particle distribution functions. For the two processes directly relevant to salinity, the electrical interaction process is represented by multicomponent ion exchange equations and the electrostatic forces in a double-layer expansion process needed for numerical simulation are calculated using the DLVO-type theory, short for Derjaguin–Landau–Verwey–Overbeek theory.

2.1. Capillarity Modeling

The presence of capillary pressure could significantly deviate the phase behaviors of reservoir fluids in unconventional reservoirs from their bulk properties. Therefore, the capillary effect has to be taken into account in order to accurately estimate phase amounts and compositions of unconventional reservoir fluids at geological conditions. Defined as the pressure difference between wetting phase and nonwetting phase, , the effect of capillary pressure can be concluded as moving the interface toward the nonwetting phase due to a positive capillary pressure and conversely due to a negative capillary pressure. The work done by capillary pressure can be formulated as

where is assumed to be constant along the interface between nonwetting and wetting fluids. Here, the Young-Laplace equation is used to calculate the capillary pressure as

with the interfacial tension , in the unit of N/M, being estimated by the Weinaug-Katz correlation

2.2. Dynamic Sorption Modeling

Usually, shale gas can exist as three states, including dissolved gas, adsorbed gas, and free gas, and the dominant reserve is made up of the adsorbed gas, which has been reported in [11] that the adsorbed gas covers up to of the total gas reservation. The dynamic balance between adsorption and desorption, as a result of shale structures and fluidity, can provide helpful knowledge to design and optimize fracturing and recovery processes. Such critical information is of significant importance to evaluate the reservoir and predict the well production. A large number of studies have been reported to analyze and model the dynamic sorption in unconventional reservoirs. Molecular accumulation, as a consequence of surface energy minimization, is often considered the main cause of adsorption on shale surface, while the van der Waals force is leading the physical sorption in potential theory. The sorption capacity is dependent on temperature, pressure, and other geochemical properties, such as the TOC content, also known as total organic carbon. Generally, if there are more organic matters in shale, a higher gas adsorption amount can be detected together with a higher surface area, total pore volume, and porosity. Moreover, permeability, which is the key factor relevant to flow and transport in porous media, is also affected by the adsorption and desorption process in gas production. For example, pressure drop facilitates gas desorption from kerogen, and on the other hand, the free gas production further decreases the pore pressure. As a result, pressure difference between the pores and bulk matrix will reinforce the desorption on the matrix surface. Plenty of isotherm models have been developed to mathematically describe the sorption mechanism, among which the most commonly used one is the Langmuir’s model

where denotes adsorbate volume, denotes pressure, and and denote the Langmuir pressure and Langmuir volume, respectively. Other isotherm models, including Freundlich model, D-R model, BET model, and Toth model, are proposed later so that the estimations of these models get closer to the experimental results. In this paper, the following equation is selected to model the dynamic sorption balance between adsorption and desorption

where and are the adsorption and desorption coefficients, respectively, with the unit of , and denotes the saturated adsorption capacity. The gas concentration can be calculated from the deformed advection-diffusion LB scheme:

In the above equation, the total gas concentration at certain location can be calculated as the summation of distribution functions in all the directions:

where denotes the distribution function in convection-diffusion problems and the source term ; in distribution balance equation, see Equation [eq: eq (6)], which represents the effect of dynamic sorption which is determined by

Correspondingly, the macroscopic velocity can be formulated as

where denotes the adsorbed amount in the site of the th direction and denotes a coefficient balancing the units.

2.3. Multicomponent Ion Exchange Modeling

The process of multicomponent ion exchange can be modeled as (using as an example of divalent cations)

A generalized model has been proposed in [12] for different ions, e.g., sodium, calcium, and magnesium cations, which are commonly seen in reservoir brine and rocks.

Two coefficients are presented in their study to model the ion exchange selectivity as

The term in the above equations represents the equivalent fraction of the cations on the exchanger, and denotes , , or .

The wettability alteration is considered to result from ion adsorptions and corresponding surface charge change, which can be modeled as [13]

where serves as a catalyst increasing ion concentration. A similar formula can be applied to brine containing ions

In high salinity injection, pH reversal may occur and the ion exchange can be modeled as

Dissolution and precipitation can happen with the interaction of injected brine and rock

The in Equation [eqim] is also obtained from the interaction in the aqueous phase as

2.4. Double-Layer Expansion Modeling

The modified DLVO theory, which describes electrostatic forces, is commonly used to model and calculate the double-layer expansion process [14, 15], and the results can be used in numerical simulation of the low salinity waterflooding process [16]. Force contributions can be modeled as

where the disjoining pressure is a summation of contributions from London van der Waals force, electrostatic force, and structural force, all of which are functions of the wetting film thickness . This formulation can be used in wettability calculation, as the following augmented Young-Laplace equation uses disjoining pressure to express the interaction equilibrium in oil/brine/rock system [15]:

where is capillary pressure. The contact angel is given by [17]

It can be referred from Equation [eqca] that if , the water film is unconditionally stable due to the constant contact angel . If the rock is water-wet, implying or , then the water film is meta-stable; otherwise, it is oil-wet with or .

The van der Waals force in Equation [eqfo] can be calculated by [18]

where is the Hamaker constant. The electrostatic double-layer force in Equation [eqfo] can be modeled using a bounded estimation of the Poisson-Boltzmann equation [19]

where denotes the dielectric constant of vacuum, denotes the relative dielectric constant of the aqueous medium, denotes the Boltzmann constant, denotes the absolute temperature (K), and denotes the ion valence. An upper limit of the estimation can be established with a constant potential assumption that attraction between the plates with the same charge but unequal potential always exists

On the contrary, if a repulsion force is assumed to exist between the plates,

The term in the above two equations is defined as the reduced potential of th plate, and it can be calculated by

The above approximation model is called the analytical compression approximation (CA) model, which is suitable for cases with low to intermediate electrostatic potential. A linear superposition approximation (LSA) is proposed, also called weak overlap approximation (WOA), as a correct answer between the two extremes as [20]

where is the zeta potential on the th surface of each plate. It is proved in [19] that the LSA model is more favorable to fit the force measurement on surface experiments. in the above models are defined as the Debye-Hückel reciprocal length and determined by the following expression:

The hydration force can be calculated by [21]

where is the force coefficient relevant with the boundary conditions and is the clay length. The hydrophobic force near the surface can be calculated using

A more general form has been proposed in [22] to calculate the structural forces for all the cases.

3. Simulation and Results

All the aforementioned mechanisms, such as capillarity, sorption, and salinity, have challenged the conventional continuum modeling and simulation using Navier-Stokes equations. Consequently, modifications and improvements have to be introduced to account for these mechanisms so that the governing systems are generated in order to obey the physical laws as well as realistic conditions. Thermodynamic equilibrium schemes are reconstructed to account for capillary effect and to meet the first and second laws of thermodynamics. The calculated equilibrium solutions can ensure a good initial estimate for multiphase multicomponent flow simulation, and the new energy and entropy balance formulations can lead to a stable convergence while tolerating a relatively much larger time step. Mesoscopic numerical approaches, representative of LBM [23] and Pore Network Method [24], are widely employed in the direct simulation of flow and transport in porous media, and some simple but effective terms can be added to take such special mechanisms into account. The new mass and momentum conservation properties can be proved rigorously, which further promote such delicate numerical approaches. However, there is still not a comprehensive model or simulation approach covering all the effect of capillarity, sorption, and salinity. In another words, more investigations are expected to develop a reliable and practical scheme to simulate engineering processes, and the concept of a digital twin can be a potential platform to cover as much mechanisms as we want. With conservation equations governing fluid flow through porous media, Darcy’s scale simulation can be performed to model the migration and transport of oil and gas and predict the oil production. Pore scale simulation is conducted to investigate the detailed mechanisms of surface interactions and to show the correlation between the microscopic details in a single pore, thermal equilibrium conditions, and macroscopic flow and transport properties. Mesoscopic simulations, like LBM and pore network modeling, are used as a bridge to link between micromechanism and macrophenomena.

3.1. Darcy’s Scale Simulation

A mass balance equation for immiscible incompressible oil-water two phase flow can be written as [25]

where and denote the velocity of oil phase and water phase, respectively. The extended Darcy’s law for aqueous and oil phases can be modeled as follows if gravity is ignored:

where and represent the pressure of oil phase and water phase, and are the oil and water relative permeability, respectively. A commonly used model to account for the effect of salinity on relative permeability and capillary pressure has been proposed in [26] as

where

In the above system, denotes the oil-water capillary pressure, is a scaling factor, and denotes the residual oil saturation after waterflooding. The superscripts and denote the high salinity water and low salinity water, respectively. It is a simple model capable of predicting the oil recovery using low salinity waterflooding at field-scale studies or single-well tests. An obvious dependency of capillary pressure and relative permeability on injection salinity is expressed in the formulas, regarding the salt species in brine only as an additional single lumped component in the aqueous phase. A balance between the two extreme conditions, lower salinity limit and upper salinity limit, is conducted using the scaling factor.

A more comprehensive model is proposed in [27] which starts from the classical Corey equations [28]

where

The first modification is assuming the residual oil saturation is a function of salinity only.

An additional modification introduces the salinity effect on the end-point water relative permeability

Next, the effect of salinity is also applied to the exponent of oil relative permeability

The salinity is denoted by in the above equations, and the two thresholds of high salinity and low salinity are represented by the superscript and .

A dimensionless system to model the macroscopic flow of low salinity waterflooding can be written as [29, 30]

where the subscript in the above system stands for dimensionless and denotes the saturation. The initial condition is provided as reservoir saturation and formation water salinity () before injection [31]

Several popular industrial software has been developed to calculate the mechanisms related with injection salinity. PHREEQC is an industry-standard geochemistry software which has been successfully applied in the study of low salinity waterflooding, with emphasis on the electrostatic reaction, ion exchange, and mineral dissolution [32, 33], and used to verify new models and approaches [12]. UTCHEM simulator is another widely accepted software in petroleum industry, developed by the University of Texas at Austin, to predict the effect of injected brine with various ion compositions, and the injected low salinity water is described using the integrated tool, UTCHEM-IPHREEQC [34]. IPHREEQC is also a state-of-art geochemical engine, and UTCHEM-IPHREEQC is an accurate, robust, and flexible tool that enables to model low salinity waterflooding and many other enhanced oil recovery techniques with respect to geochemistry. Later, another three-dimensional equation-of-state-based compositional simulator, also developed by the University of Texas at Austin, UTCOMP, is coupled with IPHREEQC and the effect of hydrocarbon components soluble in the aqueous phase on the pH buffering and other related reactions in the oil/brine/rock system [35].

3.2. Mesoscopic Simulation

In addition to LBM mentioned in Section 2, another representative mesoscopic approach, pore network modeling, has also shown promising potentials in simulating the flow and transport behaviors in unconventional reservoirs. By constructing a porous network in which pore bodies are connected through pore throats, such a model could represent highly irregular structure from the perspective of topology and geometry. After selecting certain distribution functions and key parameters to control the size of pore body and pore throat, the network is connected and a double permeability media is then constructed for further investigation. As explained in [11], the two structures, pore body and pore throat, can be treated as fracture and matrix, respectively, and this body-throat connection can be easily extended to carry on the streaming and collision process of distribution functions.

A network constituting of pore bodies is constructed as shown in Figure 2, where the black band represents matrix and white band represent fracture. It can be easily referred that this porous structure is generated by two sets of pore parameters, and the corresponding porosity and fluidity is different in these layers. The parameters of two types of media are listed in Table 1, and it can be stated that the three layers of Media 1 contain more matrix compared with the two layers of Media 2. Furthermore, a better mobility is expected in the two layers of Media 2, and more resistance may occur in the three layers of Media 1.

Figure 2: A porous media generated using pore network model with a designed structure.

Table 1: Parameters for pore network modeling.

Parameter
Media 1
Media 2
Unit

Min. pore body inscribed radius
0.0372
0.0625
mm

Max. pore body inscribed radius
0.254
0.366
mm

Mean. pore body inscribed radius
0.125
0.246
mm

Standard deviation
0.128
0.187
mm

After constructing the porous media using pore network modeling, the detailed mesoscopic algorithm can be described as follows:

(1)Generate the optimized LBGK scheme with sorption coefficients, weight matrix, medium structure, and flow scenario. Determine the inclusive parameters(2)Apply the free flow distribution function in fractures and transport distribution function in matrix. The dynamic sorption is then calculated, while at the first iteration, this adsorption is set to be zero. The free flow distribution function for fractured porous media reads as

and the distribution function at equilibrium state is

(3)Use this calculated adsorption amount to update the free flow simulation, and further calculate the diffusion and transport process

This two-scale LBM can be easily recovered back into Navier-Stokes equation and advection-diffusion equation, respectively, for fracture and matrix scale LBGK scheme by Chapman-Enskog expansion [36]. The effectiveness of this algorithm with proper modifications on the general advection-diffusion LBGK scheme and the coupling of scales using the free flow velocity can be verified with the following example with constant gas injection on the left boundary of Figure 2. The adsorption distribution at different time step is illustrated in Figure 3. It can be seen free flow is much faster in Media 2 with more fracture, while adsorption amount is much larger in Media 1 with more matrix. The black color in Figure 3 corresponds to the “fracture” structure. The result is reasonable in both scales, and it can be concluded that more adsorption in the matrix may be due to higher saturation sorption amount or adsorption coefficient (referred from Langmuir-type isothermal models) and can lead to smaller free flow velocity in fractures. On the contrary, the result of our scheme is reasonable in both media and the effect of dynamic sorption on free flow region is illustrated. The effect of porosity in both two scales, fracture and matrix, and the effect of sorption parameters in a Langmuir-type isothermal sorption model are all tested and analyzed. Generally, the increasing adsorbed amount in the matrix due to the higher adsorption coefficient or saturation sorption amount will result in a slower velocity in the free flow scale. However, if the increase of adsorbed gas amount is the result of larger matrix porosity, then the free flow velocity could be accelerated as the total resistance in the media has been decreased.

Figure 3: Adsorption distribution at time steps 12000, 22000, 76000, and 130000.

3.3. Pore Scale Simulation

Phase equilibrium calculation is essentially needed for generating a physical-meaningful initial phase distribution benefiting further multiphase flow simulation. The following fugacity equation can be used to establish thermodynamic equilibrium, which can be calculated based on various equations of state (EOS), for example, Peng-Robinson EOS [37, 38]:

Generally, volume constraint is also needed for a complete conservation relation, which could be described as

where denotes the mole density, denotes the mass density, and is the constitutive equation. Unconditional stability, the most essential property determining algorithm robustness and applicability in practice, is preserved by using the convex-concave splitting scheme on chemical potential. In details, the chemical potential is supposed to have two components and and the counterpart splitting can be written as

The unconditional stability of the above semi-implicit scheme has been proved in details in [37]. The evolution equations of mole and volume can be formulated based on Onsager’s reciprocal principle as

The computational efficiency and reliability require an energy stable system consistent with the second law of thermodynamics. Regarding the Onsager coefficient matrix , it can be divided into 4 submatrices, shown as below

Here, , , and . It is essential to ensure the positive definition of the Onsager coefficient matrix; otherwise, a modified Cholesky factorization will be introduced to preserve its positive definiteness. Generally, should be sufficiently positive and is added as a diagonal matrix with suitable positive entries. This positive definite property can keep the continuous increasing of entropy in the iterations, which will ensure to reach the local maximum using the Newton-Raphson method. The effect of capillarity can be illustrated by the difference in tangent-plane distance (TPD) function and phase envelope of fluid mixtures in porous media with various pore radius. As shown in Figure 4, the TPD function with respect to temperature range and overall molar density range is plotted for an EagleFord oil in two cases either with or without capillary effect. Within the specified molar density and temperature intervals, there is a single two-phase region and the phase boundary between single-phase and two-phase states is drawn in red and blue for the case with and without capillary effect, respectively. It can be seen that capillary pressure can significantly reshape the bulk phase envelope by suppressing the bubble point curve and meanwhile expanding the dew point curve.

Figure 4: Phase envelope for an Eagle Ford oil considering capillarity at different pore radii.

The effect of pore radius on the work done by capillary pressure can be explained in Figure 5. If capillary pressure is taken into account, the dew point pressure will be decreased in the lower branch and the suppression of dew point pressure becomes significant as pore radius decreases. Moreover, dew point pressure increases in the upper branch and deviates more significantly from the dew point curve of the bulk phase where the capillary effect is negligible. The overall effect of the dew point pressure increasing in the upper branch and decreasing in the lower branch enlarges the vapor-liquid region compared with bulk phase envelope.

Figure 5: Global minimum of TPD function, as well as corresponding phase envelope, as a function of the overall molar density and temperature for an Eagle Ford oil with (red line) or without (blue line) capillarity.

The effect of multicomponent ion exchange on relative permeability can be modeled as [39]

where is the absorbed cations and the subscripts and represent the calcium and magnesium ion, respectively. The scaling function, , is dependent on the divalent ion adsorption conditions in the precipitation and dissolution processes.

The microscopic displacement efficiency as a function of trapping number is proposed in [40], using ionic strength () calculation as

where and denote the charge and molarity of the fluid species , respectively. The thickness of double electric layer is then determined by

where and denote the relative permeability and the permeability of the free phase, respectively. A clay mineral model is built in [41] in terms of composition, structure, and charge density on the clay surface. Various salinity conditions, ranging from freshwater to seawater, have been considered, and system wetness can be either oil-wet or water-wet in that study.

4. Conclusion and Remarks

Due to the tight formations and related special mechanisms often met in unconventional reservoirs, there are still nonnegligible public concerns on the economic efficiency, production safety, and environmental friendship. As an effective approach to describe the flow and transport behaviors in subsurface reservoirs, a digital twin is designed in this paper to cover the purposes of media construction, mechanism investigation, and production estimation in physical entity. Representative mechanisms, such as capillarity, sorption, and injection salinity, have been mathematically characterized in details, and multiscale algorithms are developed to simulate the effect of these mechanisms. Physical reservations and equivalence between various scales can be preserved in the generated schemes, and several results are illustrated to prove the reliability and robustness.

More models and algorithms are expected to be included in the digital twin in the future to construct a more comprehensive numerical platform capable of simulating realistic engineering cases efficiently and accurately. Extensions to a wider range of applications are easy to perform as long as the physical correlations can be described numerically in the virtual space. Field scale studies can be enabled by the usage of parallel computing, bound-preserving, reduced-space methods and many other scale coupling techniques [42]. Molecular dynamics simulation and Monte Carlo simulation can also be added into this twin to lay a more solid theoretical foundation on fundamental microscopic mechanisms [43]. Accelerated flash calculations using deep learning algorithms have also been carried out by many researchers [44, 45], while the pore-scale flash calculation schemes have been extended to solve related engineering problems including carbon dioxide sequestration [46]. Hydraulic fracturing and rock properties are directly relevant with the exploitation of unconventional reservoirs, where numerous models have been developed to simulate the observations [47, 48]. Gathering and transportation are the connection between reservoir exploitation and market utilization, where the flow and transport in pipelines can also be modeled to resolve the engineering problem including scaling and corrosion [49, 50]. The next step is to find the proper connectors that link the many aspects of mechanisms, models, and algorithms to establish a comprehensive digital twin.

Data Availability

Data are available on request via email sent to tao.zhang.1@kaust.edu.sa.

Conflicts of Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

The work of Tao Zhang, Yiteng Li, and Shuyu Sun was supported by funding from the National Natural Scientific Foundation of China (Grants Nos. 51874262 and 51936001) and King Abdullah University of Science and Technology (KAUST) through the Grant no. BAS/1/1351-01-01.

References

Y. Shen, H. Ge, X. Zhang, L. Chang, D. Liu, and J. Liu, “Impact of fracturing liquid absorption on the production and water-block unlocking for shale gas reservoir,” Advances in Geo-Energy Research, vol. 2, no. 2, pp. 163–172, 2018.

J. Cai, D. Lin, H. Singh, W. Wei, and S. Zhou, “Shale gas transport model in 3D fractal porous media with variable pore sizes,” Marine and Petroleum Geology, vol. 98, pp. 437–447, 2018.

H. Wang, L. Chen, Z. Qu et al., “Modeling of multi-scale transport phenomena in shale gas production—a critical review,” Applied Energy, vol. 262, article 114575, 2020.

D. Wang, X. Wang, H. Ge, D. Sun, and B. Yu, “Insights into the effect of spontaneous fluid imbibition on the formation mechanism of fracture networks in brittle shale: an experimental investigation,” ACS Omega, vol. 5, no. 15, pp. 8847–8857, 2020.

G.-p. Zhu, J. Yao, H. Sun et al., “The numerical simulation of thermal recovery based on hydraulic fracture heating technology in shale gas reservoir,” Journal of Natural Gas Science and Engineering, vol. 28, pp. 305–316, 2016.

R. Christopher, “Anomalous diffusion or classical diffusion in an anomalous reservoir? Evaluation of the impact of multi-phase flow on reservoir signatures in unconventional reservoirs,” in Unconventional Resources Technology Conference, pp. 22–24, Denver, Colorado, July 2019.

S. Sun and T. Zhang, “A 6M digital twin for modeling and simulation in subsurface reservoirs,” Advances in Geo-Energy Research, vol. 4, no. 4, pp. 349–351, 2020.

F. Tao, J. Cheng, Q. Qi, M. Zhang, H. Zhang, and F. Sui, “Digital twin-driven product design, manufacturing and service with big data,” The International Journal of Advanced Manufacturing Technology, vol. 94, no. 9-12, pp. 3563–3576, 2018.

S. Boschert and R. Rosen, “Digital twin—the simulation aspect,” in Mechatronic futures, pp. 59–74, Springer International Publishing, 2016.

S. Sun, “Darcy-scale phase equilibrium modeling with gravity and capillarity,” Journal of Computational Physics, vol. 399, article 108908, 2019.

T. Zhang, S. Sun, and H. Song, “Flow mechanism and simulation approaches for shale gas reservoirs: a review,” Transport in Porous Media, vol. 126, no. 3, pp. 655–681, 2019.

C. Dang, L. Nghiem, N. Nguyen, Z. Chen, and Q. Nguyen, “Mechanistic modeling of low salinity water flooding,” Journal of Petroleum Science and Engineering, vol. 146, pp. 191–209, 2016.

P. Zhang, M. T. Tweheyo, and T. Austad, “Wettability alteration and improved oil recovery by spontaneous imbibition of seawater into chalk: impact of the potential determining ions Ca²⁺, Mg²⁺, and SO₄²⁻,” Colloids and Surfaces A: Physicochemical and Engineering Aspects, vol. 301, no. 1-3, pp. 199–208, 2007.

H. Ding and S. Rahman, “Experimental and theoretical study of wettability alteration during low salinity water flooding-an state of the art review,” Colloids and Surfaces A: Physicochemical and Engineering Aspects, vol. 520, pp. 622–639, 2017.

Z. Jalili and V. A. Tabrizy, “Mechanistic study of the wettability modification in carbonate and sandstone reservoirs during water/low salinity water flooding,” Energy and Environment Research, vol. 4, no. 3, 2014.

P. C. Myint and A. Firoozabadi, “Thin liquid films in improved oil recovery from low-salinity brine,” Current Opinion in Colloid & Interface Science, vol. 20, no. 2, pp. 105–114, 2015.

S. Basu and M. M. Sharma, “Measurement of critical disjoining pressure for dewetting of solid surfaces,” Journal of Colloid and Interface Science, vol. 181, no. 2, pp. 443–455, 1996.

H. C. Hamaker, “The London—van der Waals attraction between spherical particles,” physica, vol. 4, no. 10, pp. 1058–1072, 1937.

J. N. Israelachvili, Intermolecular and surface forces, Academic press, 2015.

J. Gregory, “Interaction of unequal double layers at constant charge,” Journal of Colloid and Interface Science, vol. 51, no. 1, pp. 44–51, 1975.

J. A. Molina-Bolivar and J. L. Ortega-Vinuesa, “How proteins stabilize colloidal particles by means of hydration forces,” Langmuir, vol. 15, no. 8, pp. 2644–2653, 1999.

G. J. Hirasaki, “Wettability: fundamentals and surface forces,” SPE Formation Evaluation, vol. 6, no. 2, pp. 217–226, 1991.

H. Wang, X. Yuan, H. Liang, Z. Chai, and B. Shi, “A brief review of the phase-field-based lattice Boltzmann method for multiphase flows,” Capillarity, vol. 2, no. 3, pp. 33–52, 2019.

A. Golparvar, Y. Zhou, K. Wu, J. Ma, and Z. Yu, “A comprehensive review of pore scale modeling methodologies for multiphase flow in porous media,” Advances in Geo-Energy Research, vol. 2, no. 4, pp. 418–440, 2018.

L. W. Lake, R. Johns, W. Rossen, and G. Pope, Fundamentals of Enhanced Oil Recovery, Prentice Hall, Englewood Cliffs, NJ, 2014.

G. R. Jerauld, K. J. Webb, C.-Y. Lin, and J. Seccombe, “Modeling low-salinity waterflooding,” in SPE Annual Technical Conference and Exhibition, pp. 24–27, San Antonio, Texas, USA, 2006.

I. Tripathi and K. K. Mohanty, “Instability due to wettability alteration in displacements through porous media,” Chemical Engineering Science, vol. 63, no. 21, pp. 5366–5374, 2008.

M. Delshad and G. A. Pope, “Comparison of the three-phase oil relative permeability models,” Transport in Porous Media, vol. 4, no. 1, pp. 59–83, 1989.

A. Zeinijahromi, T. K. P. Nguyen, and P. Bedrikovetsky, “Mathematical model for fines-migration-assisted waterflooding with induced formation damage,” SPE Journal, vol. 18, no. 3, pp. 518–533, 2013.

F. Hussain, A. Zeinijahromi, P. Bedrikovetsky, A. Badalyan, T. Carageorgos, and Y. Cinar, “An experimental study of improved oil recovery through fines-assisted waterflooding,” Journal of Petroleum Science and Engineering, vol. 109, pp. 187–197, 2013.

S. Borazjani, A. Behr, L. Genolet, A. van der Net, and P. Bedrikovetsky, “Effects of fines migration on low-salinity waterflooding: analytical modelling,” Transport in Porous Media, vol. 116, no. 1, pp. 213–249, 2017.

D. A. Afekare and M. Radonjic, “From mineral surfaces and coreflood experiments to reservoir implementations: comprehensive review of low-salinity water flooding (LSWF),” Energy & Fuels, vol. 31, no. 12, pp. 13043–13062, 2017.

A. Kazemi Nia Korrani, G. R. Jerauld, and K. Sepehrnoori, “Mechanistic modeling of low-salinity waterflooding through coupling a geochemical package with a compositional reservoir simulator,” SPE Reservoir Evaluation & Engineering, vol. 19, no. 1, pp. 142–162, 2016.

K. N. Korrani, K. S. Aboulghasem, and M. Delshad, “A novel mechanistic approach for modeling low salinity water injection,” in SPE Annual Technical Conference and Exhibition, New Orleans, Louisiana, USA, 2013.

A. K. N. Korrani, G. R. Jerauld, and K. Sepehrnoori, “Coupled geochemical-based modeling of low salinity waterflooding,” in SPE Improved Oil Recovery Symposium, Tulsa, Oklahoma, USA, 2014.

T. Zhang and S. Sun, “A coupled Lattice Boltzmann approach to simulate gas flow and transport in shale reservoirs with dynamic sorption,” Fuel, vol. 246, pp. 196–203, 2019.

J. Kou and S. Sun, “Thermodynamically consistent simulation of nonisothermal diffuse-interface two-phase flow with Peng–Robinson equation of state,” Journal of Computational Physics, vol. 371, pp. 581–605, 2018.

J. Kou and S. Sun, “A stable algorithm for calculating phase equilibria with capillarity at specified moles, volume and temperature using a dynamic model,” Fluid Phase Equilibria, vol. 456, pp. 7–24, 2018.

A. V. Omekeh, H. A. Friis, I. Fjelde, and S. Evje, “Modeling of ion-exchange and solubility in low salinity water flooding,” in SPE Improved Oil Recovery Symposium, Tulsa, Oklahoma, USA, 2012.

E. W. Al-Shalabi, K. Sepehrnoori, G. Pope, and K. Mohanty, “A fundamental model for predicting oil recovery due to low salinity water injection in carbonate rocks,” in SPE Energy Resources Conference, Port of Spain, Trinidad and Tobago, 2014.

T. Underwood, V. Erastova, P. Cubillas, and H. C. Greenwell, “Molecular dynamic simulations of montmorillonite–organic interactions under varying salinity: an insight into enhanced oil recovery,” The Journal of Physical Chemistry C, vol. 119, no. 13, pp. 7282–7294, 2015.

H. Yang, S. Sun, Y. Li, and C. Yang, “A scalable fully implicit framework for reservoir simulation on parallel computers,” Computer Methods in Applied Mechanics and Engineering, vol. 330, pp. 334–350, 2018.

Y. Yang, A. K. Narayanan Nair, and S. Sun, “Adsorption and diffusion of carbon dioxide, methane, and their mixture in carbon nanotubes in the presence of water,” The Journal of Physical Chemistry C, vol. 124, no. 30, pp. 16478–16487, 2020.

K. Wang, J. Luo, Y. Wei, K. Wu, J. Li, and Z. Chen, “Artificial neural network assisted two-phase flash calculations in isothermal and thermal compositional simulations,” Fluid Phase Equilibria, vol. 486, pp. 59–79, 2019.

T. Zhang, Y. Li, Y. Li, S. Sun, and X. Gao, “A self-adaptive deep learning algorithm for accelerating multi-component flash calculation,” Computer Methods in Applied Mechanics and Engineering, vol. 369, article 113207, 2020.

Y. Li, Z. Qiao, S. Sun, and T. Zhang, “Thermodynamic modeling ofCO2solubility in saline water using NVT flash with the cubic-Plus-association equation of state,” Fluid Phase Equilibria, vol. 520, article 112657, 2020.

S. Al-Rbeawi, “The performance of complex-structure fractured reservoirs considering natural and induced matrix block size, shape, and distribution,” Journal of Natural Gas Science and Engineering, vol. 81, article 103400, 2020.

S. Al-Rbeawi and J. F. Owayed, “Fluid flux throughout matrix-fracture interface: discretizing hydraulic fractures for coupling matrix Darcy flow and fractures non-Darcy flow,” Journal of Natural Gas Science and Engineering, vol. 73, article 103061, 2020.

H. Bai, “Mechanism analysis, anti-corrosion techniques and numerical modeling of corrosion in energy industry,” Oil & Gas Science and Technology–Revue d’IFP Energies nouvelles, vol. 75, 2020.

T. Zhang, Y. Li, C. Li, and S. Sun, “Effect of salinity on oil production: review on low salinity waterflooding mechanisms and exploratory study on pipeline scaling,” Oil & Gas Science and Technology–Revue d’IFP Energies nouvelles, vol. 75, 2020.