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

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


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][4][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 welldesigned 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). Multi-component 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.

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.
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 2 Geofluids 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.

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, p c = p n − p w , 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 p c 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 80% 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 V denotes adsorbate volume, P denotes pressure, and P L and V L 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

Geofluids equation is selected to model the dynamic sorption balance between adsorption and desorption
where k a and k d are the adsorption and desorption coefficients, respectively, with the unit of S −1 , and V m denotes the saturated adsorption capacity. The gas concentration C 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 gc i denotes the distribution function in convectiondiffusion problems and the source term S; 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 S i denotes the adsorbed amount in the site of the ith direction and z denotes a coefficient balancing the units.

Multicomponent Ion Exchange
Modeling. The process of multicomponent ion exchange can be modeled as (using C a 2+ 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 ζðMI − XÞ term in the above equations represents the equivalent fraction of the cations on the exchanger, and MI denotes Na + , Ca 2+ , or Mg 2+ .
The wettability alteration is considered to result from ion adsorptions and corresponding surface charge change, which can be modeled as [13] where SO 2− 4 serves as a catalyst increasing Ca 2+ ion concentration. A similar formula can be applied to brine containing 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 HCO − 3 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, 4 Geofluids and structural force, all of which are functions of the wetting film thickness h. 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 p c is capillary pressure. The contact angel θ is given by [17] cos θ = 1 + 1 It can be referred from Equation [eqca] that if I > 0, the water film is unconditionally stable due to the constant contact angel θ = 0 ∘ . If the rock is water-wet, implying 90 >°θ > 0°or 0 < I < −σ ow , then the water film is meta-stable; otherwise, it is oil-wet with 180 >°θ > 90°or −σ ow < I < −σ ow .
The van der Waals force Π VDW in Equation [eqfo] can be calculated by [18] where ε 0 denotes the dielectric constant of vacuum, ε r denotes the relative dielectric constant of the aqueous medium, k B denotes the Boltzmann constant, T denotes the absolute temperature (K), and z 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 ψ ri in the above two equations is defined as the reduced potential of ith 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 ζ i is the zeta potential on the ith 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 C is the force coefficient relevant with the boundary conditions and h 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.

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 5 Geofluids 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.

Darcy's Scale Simulation.
A mass balance equation for immiscible incompressible oil-water two phase flow can be written as [25] ϕ ∂s ∂t where u o and u w 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 p o and p w represent the pressure of oil phase and water phase, k ro and k rw 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, p cow denotes the oil-water capillary pressure, θ is a scaling factor, and S orw denotes the residual oil saturation after waterflooding. The superscripts HS and LS 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 X c in the above equations, and the two thresholds of high salinity and low salinity are represented by the superscript HS and LS.

Geofluids
A dimensionless system to model the macroscopic flow of low salinity waterflooding can be written as [29,30] where the subscript D in the above system stands for dimensionless and s denotes the saturation. The initial condition is provided as reservoir saturation and formation water salinity (γ) before injection [31] t D = 0 : s = s I , γ = γ I : ð41Þ 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 threedimensional 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].

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 500 × 500 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.
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 porous media in shale Figure 2: A porous media generated using pore network model with a designed structure.
(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.

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 N 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 μ i ðnÞ is supposed to have two components μ convex i ðnÞ and μ concave i ðnÞ 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 1 , and C = ∂ðp 1 − p 2 Þ/∂V 1 . 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, Ψ + E should be sufficiently positive and E 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 ½250,700 K and overall molar density range ½0, 10000 mol/m 3 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 twophase 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.
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.
The effect of multicomponent ion exchange on relative permeability can be modeled as [39]

Geofluids
where β is the absorbed cations and the subscripts Ca and Mg represent the calcium and magnesium ion, respectively. The scaling function, F, 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 (I s ) calculation as where z i and m i denote the charge and molarity of the fluid species i, respectively. The thickness of double electric layer is then determined by where ε r and ε 0 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 waterwet in that study.

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 sub-surface 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, reducedspace 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.