Energy extraction from conventional resources involves producing crude oil, natural gas, and its condensates from Copyright © 2020, for this paper by its authors. Use permitted under Creative Commons License Attribution 4.0 International (CCBY 4.0). rock formations that have high porosity and permeability (Bahadori, 2017) . These rock formations are found below an impermeable rock. However, energy extraction from unconventional hydrocarbon resources (Ahmmed and Meehan, 2016) involves using advanced drilling and stimulation techniques (e.g., long horizontal laterals and multi-stage hydraulic fracturing) to extract crude oil and natural gas that are trapped in the pores of relatively impermeable sedimentary rocks (e.g., shale, tight sandstones).

Typically, unconventional reservoirs have porosity in the range of 0.04-0.08 and matrix permeability on the order of nanodarcies (10-16-10-20 m2) (Rezaee, 2015; Belyadi et al., 2019) . Instead of the porous flow that dominates conventional reservoirs, fracture flow dominates unconventional reservoirs, with natural fractures dissecting the matrix and intersecting with the hydraulic fractures. As result, energy extraction is more difficult than conventional reservoirs. Model-based optimization of unconventional reservoirs is also challenging because due to the long horizontal laterals there is insufficient site data to inform high-fidelity physics models (Mohaghegn, 2017; Belyadi et al., 2019) . Despite these challenges and due to the abundance of unconventional resources, with reserves projected to last for many decades, energy extraction from these resources have gained prominence in recent years (Briefing, 2013 and Weijermars, 2014) . Current extraction efficiency from unconventional reservoir is very low (~5-10%) for tight oil and ~20% for shale gas (Sandrea, 2007; Muggeridge et al., 2014) compared to conventional reservoirs (~20-40%) (Zitha et al., 2008) . This is because the impact of resource development processes (e.g., slow drawdown or fast drawdown) and underlying physics that determine the energy extraction from the impervious rocks are poorly understood (Rezaee, 2015; Belyadi et. Al., 2019) . State-of-the-art workflows for unconventional reservoir management are data-driven, which perform poorly beyond their training regimes. Hence, innovative extraction strategies (e.g., pressure-drawdown management) coupled with advanced workflows (e.g., physics-informed machine learning) are needed to improve the hydrocarbon recovery efficiency (Seales et al., 2017; Lougheed et al., 2017; Mirani et al., 2018) . In this paper, we present a physics-informed machine learning (PIML) workflow (Fig.1) to address unconventional production for real-time reservoir management. One of the goals of this PIML workflow (Fig.2) is to develop fast and accurate ML-models grounded in physics for real-time history matching and production forecasting in a fracture shale gas reservoir.

Current workflows for unconventional reservoirs are predominantly based on production decline curve analysis and its extensions (Wu et al., 2013; Sun 2015) , data-driven machine learning (ML) approaches (Holdaway 2014; Chaki 2015; Puyang et al., 2015; Carvajal et al. 2017; Mohaghegn, 2017) , and/or extension of physics-based conventional reservoir workflows (Rezaee, 2015; Rajput and Thakur, 2016; Belyadi et al., 2019) . Decline curve analysis provides empirical models to forecast production data based on the past production history. This type of approach lacks physics and in-depth knowledge of the site behavior is not included in the forecasting models. The data-driven ML approaches perform poorly when faced with uncertain, missing, and sparse data – common problems with existing datasets related to unconventional reservoirs. Moreover, the data-driven ML analyses perform poorly in making forecasts outside of their training regimes, and the exploration of novel production strategies fundamentally requires extrapolation (where ML struggles) as opposed to interpolation (where ML excels). The physics-based workflows adopted for modeling conventional reservoirs use extensively available site-characterization data (which is acquired over months or years) for history-matching. Even though large unconventional reservoir data (e.g, fiber-optics) are sampled at sparse locations, simply combining all the data together would not improve the accuracy of real-time forecasting. This is because reservoir conditions change considerably from one basin to another basin. Therefore, physical constraints need to be incorporated in workflows. Existing workflows employ high-fidelity physics models to perform simulations, which are expensive to run. For example, it takes several days to months to run reservoir-scale model simulations (Rezaee, 2015; Rajput and Thakur, 2016; Belyadi et al., 2019) with degreesof-freedom in the (108) on state-of-the-art HPC machines. As a result, these conventional reservoir workflows are not ideal for usage in comprehensive uncertainty quantification studies, which require 1000s of forward model runs. To overcome the key gaps associated with existing workflows, we propose a PIML workflow (Fig.1 and Fig.2) to accelerate the development of ML-models while constraining it with physics. The aim is (1) to develop a library from a combination of site observations and physics-based models that is representative of unconventional reservoir site behavior, and (2) to develop fast and accurate ML-models for realtime history matching to refine key site parameters and forecast production quantities of interest (QoIs) with uncertainty estimates. Key production QoIs include, total cumulative production of hydrocarbons, gas and water production rates, stage-specific production, and spatial evolution of residual gas, reservoir pressure, temperature, and stress fields based on a user-defined pressure-drawdown strategy.

The proposed approach to develop the PIML workflow consists of three key steps: (Step-1) Development of site behavior libraries based on fast and accurate physics, (Step-2) Development of ML-based inverse models to infer key site parameters, and (Step-3) Development of fast forward models that combine physical models with ML to forecast production and reservoir conditions.

PIML Workflow Step-1: Development of a site behavior library involves generating synthetic data for a range of possible site characteristics. This includes a large number of runs from a fast physics-based reduced-order model and a smaller number of runs from a high-fidelity, full physics model. The fast physics models (e.g., Patzek models) allow us to quickly model and build a site library on the evolving reservoir data. These models allow us to efficiently explore the parameter space. Moreover, combining fast physics models with ML allows us to identify the important physical processes or dominant mechanisms that must be represented during full physics simulations. The dominant mechanisms at different stages of production include • Early stages of production (<1 year): Primary fracture creation, geometry, and network connectivity, primary fracture behavior, propped fracture behavior, anisotropic permeability of fractures. • Middle stages of production (~1-5 years): Secondary fractures and their permeabilities, shear fracture geometries, and geochemical impacts of hydraulic fluids and formation water. • Late stages of production (~5-10 years): Matrix porosity and transport properties, water imbibition impacts, adsorbed gas properties, pore structure distribution.

These mechanisms are not accounted for in the fast physics models but are needed to describe the reservoir behavior at different stages of gas and water production. However, simulating these mechanisms using full physics models for entire parameter space is computationally intractable. As a result, a fast physics model library along with ML are used to guide and improve the development of the full physics model library (which consist of complex 3D simulations of matrix-fracture interactions) for characterizing site behavior.

PIML Workflow Step-2: The ML-inverse model is developed on this site behavior library using a transfer learning approach (Pan et al., 2009; Goodfellow et al., 2016; Yamada et al., 2018) where the numerous runs from the reduced order model are used to train an initial inverse model, then the smaller number of runs from the high-fidelity model are used to fine-tune the inverse model. This effectively represents a multi-fidelity approach to training the ML inverse model. This ML-inverse model provides capabilities for real-time history matching to update the key parameters (e.g., rock permeability, rock porosity, gas transport properties). Sensitivity analysis (e.g., Sobel indices, Random Forests) is performed to provide quantitative information on the key sensitive parameters (e.g, matrix permeability and porosity, fracture network parameters) that influence shale gas production rates.

PIML Workflow Step-3: The physics-based reduced order forward model uses these calibrated parameters for real-time forecasting. This model is trained on the site libraries along with evolving production data. Moreover, it allows for various operational decisions (e.g, slow drawdown vs. fast drawdown) to be evaluated relative to future outcomes. The ML-inverse and physics-based reduced order forward model can be combined to provide uncertainty estimates on the production quantities of interest (e.g., remaining hydrocarbon-in-place, spatial evolution of pressure and temperature, shale gas and water production as a function of time).

The key challenges to accelerate the proposed PIML workflow include: 1Q. How do we enhance the information content and fill the gaps in the limited unconventional site-data for PIML analyses? Specifically, how to use the short-time gas production data (e.g., 30-120 days) to forecast the long-term performance (e.g., 1-5 years) of an unconventional shale gas well? 2Q. What is the minimum number of high-fidelity physics model simulations (e.g., PFLOTRAN, FEHM, dfnWorks) needed to develop a gas production library that is representative of shale gas sites behavior? How do we improve the full physics models to accurately represent the site behavior? 3Q. How can we use the information learned from reducedorder physics models (e.g., Patzek model, graph-based models) to inform high-fidelity physics model simulations? 4Q. What are the key sensitive parameters (e.g., matrix and fracture properties) in high-fidelity physics models that influence short-term and long-term gas production rates? 5Q. What are the key operational parameters (e.g., pressure drawdown strategies) that can be used to inform decisions in real-time, leading to optimized production?

Our hypotheses/approaches to address the key challenges include: 1A. Augment limited unconventional reservoir data (e.g., PetroMehras) can be used to generate synthetic data for the short-term production data, well logs) with lab-scale experimental data, reduced-order physics simulation data, and high-fidelity physics simulation data. Short-term production data (e.g., 30-120 days) contains statistical information (e.g., matrix and fracture network properties) that can help us predict the long-term production data (e.g., 15 years). 2A. Reduced-order physics models provide information on the key sensitive parameters needed to inform high-fidelity physics model simulations. As we obtain more production data and/or site data, a ML-based active-feedback loop is used to improve full physics models. In this active feedback loop, ML-inverse model is used to update and constrain the key parameters based on newly available reservoir data. Based on these updated key parameters, site-behavior libraries are also updated to account for new production data.

ity physics. 3A. Transfer learning can be used to provide link and transfer information from reduced-order physics to high-fidel4A. Sensitivity analysis (e.g., Sobel indices, Random Forests) can provide quantitative information on the key sensitive parameters (e.g, matrix permeability, matrix lithology, fracture length and orientation, stage spacing, hydrocarbon-in-place) that influence short-term and long-term gas production rates through feature importance. 5A. Maximizing early production may not maximize total recovery efficiency. Optimal pressure management strategies are needed to enhance recovery efficiency. One of our hypotheses is that slower drawdown rates can lead to improved recovery efficiency in long-term shale gas production (e.g., 5-10 years).

suite, Landmark Nexus, MRST, BOAST, OPM) (see ref. site behavior library.

The reduced-order physics models (Patzek et al., 2013) to generate synthetic data are given by

− ̃ ( ) where ̃ is the dimensionless time, is the dimensionless distance, and

is the real gas pseudopressure. Eq.(1) corresponds to reduced-order physics of gas flow in fractured porous rock. These are given as follows ( ) =

+ (1 − ) where is the characteristic inference time, is the halfdistance between the hydrofractures, is the initial hydraulic diffusivity, is the fractured rock effective permeability dependent on reservoir pressure, is the fractured rock effective porosity,

is the fraction of pore space occupied by the gas, is the gas viscosity, is the isothermal compressibility of gas, is the differential equilibrium partitioning coefficient of gas, and is the compressibility factor of gas. These gas properties are dependent on evolving reservoir pressure and temperature. Eq. (1) is a nonlinear pressure diffusion equation, which is solved numerically. The cumulative production of gas mass is given as follows ( ̃) = ℳ =0

′ 0 ̃

where ℳ is the initial hydrocarbon-in-place.

The expensive full physics models to simulate gas flow and transport in fractured porous media (Rezaee, 2015; Salama et al., 2017; Belyadi et al., 2019) are given by , (2) (3) (4) (5) (6) where is the density of the gas which is dependent on the reservoir pressure, is the Darcy flux, is the gas concen ̃ = = =

, = 2 = 1 2 ( , , ) = 2 2

ℎ ( ) ( , ) ( , ) tration, is the tortuosity, and is the fracture rock effective diffusivity. Eq.(4) and Eq.(5) model the gas flow under varying reservoir and bottom hole pressures. The underlying assumptions include Darcy’s flow and Fick’s law. Adsorption and non-pore refinement effects on phase behavior are ignored. Eq.(6) models the gas transport from fractured stage to horizontal well based on the initial hydrocarbon-inplace. The amount of gas extracted from the pair of hydrofractures at a given section of the well is equal to . Different types of equation of state (EOS) models are used to evaluate the gas density. These include ideal gas law, exponential pressure-dependent model, and Redlich-Kwong-Soave model. The corresponding EOS model expressions are given by = Note that it is possible to incorporate nanopore confinement effects (e.g,shifts in bubble or dew points) into EOS models to account for density changes. For examp le, see ref. Islam et al., 2015 , Tan and Piri, 2015 , and Liu and Zhang, 2019 .

Through these multi-fidelity physics models, the reservoir physical behavior is captured accurately. Gas flow and transport mechanisms are accounted through conservation of mass and equation of state for real gases. State-of-the-art simulators (e.g., PFLOTRAN, dfnWorks) are used to develop high-fidelity simulation data. These simulators use finite volume methods and Newton-Raphson method to solve the discretized system of nonlinear equations given by Eq.(4)-(7). Moreover, these simulators account for accurate meshing of fractures, matrix, and upscaling of fracture network properties for reservoir-scale high-fidelity physics simulations.

Fig.3 shows the PIML workflow to create efficient inverse models to infer key site parameters from production information. First, we sample the relevant regions in parameter space. Two site behavior libraries are developed based on this sampling. One comes from running the reduced-order physics model with a large set of samples and the other comes from running the full physics model with a smaller set of samples. The first library (based on the reduced-order physics model) is used to train an initial ML-inverse model. The ML-inverse model is then fine-tuned with the library from the high-fidelity physics model using transfer learning, producing a final ML-inverse model. During this fine-tuning process the weights of the neural network are fine-tuned. The ML-inverse model takes past production as input and produces physical parameters as output. These physical parameters can then be fed into the reduced-order physics model. The loss function for this ML-inverse model is defined in terms of how well it works in combination with the reduced-order physics model at predicting future production. Note that the second library can also be augmented with field production data to improve the realism of the MLinverse model. This fine-tuning represents a multi-fidelity approach to machine learning where a large dataset is generated with a reduced-order physics model and a smaller dataset is generated with a high-fidelity, expensive full physics model. This multi-fidelity approach allows us to perform real-time history matching on new production data to determine critical site parameters that can be used to accurately predict future production. Fig.4 shows a DFN model of a single stage based on field data from the Marcellus Shale Energy and Environment Laboratory (MSEEL) shale gas site. This model was built using dfnWorks, which is a computationally expensive full physics software suite used to generate high-resolution representations of DFNs (Hyman et a l., 2015 ). This high-fidelity meshing of the fracture network is critical to accurately capture the physical processes in a fractured shale gas reservoir. To capture the matrix effects, we generate an octreerefined continuum mesh (grey color in Fig.4) based on the DFN. The original DFN model in Fig.4 consists of three hydraulic fractures and a swarm of natural fractures that are connected to hydraulic fractures. While generating the continuum mesh, the DFN model is simultaneously upscaled to account for matrix-fracture interactions, which results in accurate permeability and porosity values that are needed to simulate gas flow and transport in fractured shale. The final mesh contains approximately 500,000 mesh cells.

Fig.5 shows the flow and transport simulation using PFLOTRAN with a Barton-Bandis stress relationship (Barton and Bandis, 1990) . This figure shows the drainage over a period of 10 years. For gas flow simulations (left), the well is maintained at 12MPa and the reservoir initial pressure is at 21MPa. For gas transport simulations (right), the initial hydrocarbon-in-place is assumed to spread in the entire fracture stage and the figure shows the transport of hydrocarbon to the well over time at two different vertical heights. The main inference from these simulations is that the characteristic behavior of drainage is tied to hydraulic fractures. Our future work involves accounting for heterogeneity of hydrocarbon distribution in the matrix for production forecasts.

Fig.6 shows encouraging results of long-term production forecasts using ML-forward models. The red color represents the short-term production, which is used along with the site behavior library built on reduced-order physics models to infer key site parameters (e.g., initial hydrocarbon-inplace, hydraulic diffusivity). History-matching is performed on the short-term production data. That is, the ML-inverse model is trained on 90 days of production data (red color). Then, the resulting ML-model with site behavior library based on reduced-order physics is used to predict future production (green color). The markers represent the long-term production data and the solid green color line represents the ML-model forecasts. Quantitatively, the prediction accuracy based on R2-score is > 90%. From this figure, it is clear that ML-model predictions, which are combined with physics are able to accurately represent the long-term production data.

We note that if another site/formation (e.g., Woodford, Haynesville, Fayetteville, Barnett, Utica, EagleFord) has a similar parameter range, we can use a set of ML techniques derived from transfer learning to model this site/formation. Transfer learning helps us to transfer knowledge gained from one site (e.g., Marcellus) to another site (e.g., Woodford, Barnett). But the developed ML-models for one site (e.g., MSEEL) need fine-tuning (or minimal retraining the neural networks) to transfer knowledge across shale sites/formations. This is not burdensome when compared to developing a new site behavior library and ML-model for a different site altogether. The transfer learning approach is attractive for tasks where reusability of ML-models for similar types is of great importance.

If the production curve contains high frequency content or oscillations, which might need addressing, there are various ways to incorporate this in the ML analyses and physicsbased models. For example, from the production curve and bottom hole pressure, we can extract dominant frequencies or a range of frequencies we are interested in modeling through Fast Fourier Transformation (FFT). This FFT transformation of bottom hole pressure and production curve vs. time provides in-depth information on quantitative aspects of high-frequency oscillations to be incorporated in physics models (e.g., ℎ = 1sin 1 + 2sin 2 + … + −1sin −1 + sin ) when developing the site behavior library and pressure management strategies (e.g., drawdown frequencies).

In this paper, we have presented a PIML workflow for realtime history matching and forecasting of gas production QoIs. The workflow coupled the strengths of machine learning with the predictability of physics-based models for realtime history-matching and forecasting. The PIML workflow used short-term production data and site behavior libraries to perform real-time history matching. Site behavior libraries are developed based on many runs of a reduced-order physics models and a smaller number of runs of expensive full physics models. The initial ML-inverse model trained on the reduced physics site library and short-term production data provided us key site parameters, which are hydraulic diffusivity and initial hydrocarbon-in-place. This initial ML-inverse model is then fine-tuned with the library from the expensive full physics models using transfer learning. The expensive full physics model simulations were developed using dfnWorks and PFLOTRAN simulators. These high-fidelity simulations account for matrix-fracture interaction, which is needed to accurately simulate gas flow and transport in fractured shale reservoirs. Moreover, from these simulations we inferred that the characteristic behavior of drainage is tied to hydraulic fractures. This high-fidelity simulation library was used to fine-tune the ML-inverse model. Our ongoing work seeks to advance and test the PIML workflow, site behavior libraries, and ML-models for pressure management (e.g., slow drawdown vs. fast drawdown) to optimize recovery at MSEEL. Our preliminary results shown in this paper are encouraging for use of site behavior libraries with ML-inverse model to address this problem.

This work was funded by the U.S. Department of Energy (DOE) – Office of Fossil Energy (FE) through its Oil and Natural Gas Research Storage Program as implemented by the National Energy Technology Laboratory. We would specifically like to thank Jared Ciferno for his support and guidance. LANL is operated by Triad National Security, LLC, for the National Nuclear Security Administration of U.S. Department of Energy (Contract No. 89233218CNA000001). The opinions expressed in this paper are those of the authors and do not necessarily reflect that of the sponsors.

MSEEL: http://mseel.org/ PFLOTRAN: https://www.pflotran.org/

Figure-2: Details of PIML workflow. The workflow utilizes a library of data on site behavior to inform forecasting models. The scale of physical parameters used in the development of site behavior libraries are stage spacing (~100-200m), hydraulic fracture length (~100-150m), a stage may contain 3-4 hydraulic fractures and a swarm of natural fractures that are connected to hydraulic fractures, fracture rock reservoir permeability (~0.1-0.9mD), fractured rock porosity (~0.04-0.08), reservoir pressures (~2030MPa), well flowing pressures (~5-15MPa), and gas properties (e.g., density, saturation, viscosity, compressibility).