History Matching
The reservoir simulation module is designed to model depletion performance in tight unconventionals. The OPM Flow reservoir simulator, an opensource simulator using Eclipse 100 input syntax, is the engine of the tool.
1. Theory
1.1 Wellbox Model
The simulations are performed on a prebuilt model that represents the typical well in tight unconventionals, often referred to as a multifractured horizontal well (MFHW). The prebuilt model is called the wellbox model. A cartoon of the wellbox model is shown in Fig. 1.
Fig. 1: Wellbox model
The wellbox model will have dimensions that the user can change, but a set of default values will be given initially. The dimensions and default values are as follows:
Dimension  Symbol  Default Value 

Lateral Length  \(L_w\)  5280 ft 
Reservoir HalfExtent  \(x_e\)  800 ft 
No. of Fractures  \(N_f\)  100 
Fracture HalfLength  \(x_f\)  400 ft 
Reservoir and Fracture Height  \(h_f\)  150 ft 
As the wellbox model is symmetric, the numerical model will be a onequarter fracture, as highlighted by a redlined box in Fig. 1. The production from the numerical model is upscaled to by using a scaling factor. The scaling factor is defined as
where \(N_f\) is the number of hydraulic fractures (clusters) in the wellbox model. The dimensions of the numerical model are the same as the wellbox model, except for the ydimension, where instead of \(L_w\), we use \(y_e\) defined as
1.2 Enhanced Fracture Region (EFR) Option
It is possible to introduce a region of enhanced, effective permeability, \(k_{EFR}\), around each fracture with extent, \(y_{EFR}\). An example description for a singlefracture is illustrated in Fig. 2. The reservoir can extend beyond the fracture tip, and have contributions beyond the fracture height. The enhanced fracture region has the same height and length as the fracture, only difference being the extent in the y direction.
Fig. 2: Enhanced fracture region model
1.3 Reservoir Properties
1.3.1 GridBlock Properties
The reservoir properties will include some properties that can be changed by the user, and some properties that are locked. The changeable reservoir properties are given the following table:
Property  Symbol  Default Value 

Matrix Porosity  \(\phi_m\)  0.03 
Matrix Permeability  \(k_m\)  100 nd 
Matrix Initial HC Saturation (Fluid Dep.)  \(S_{HCi}\)  1 
Matrix Initial Water Saturation  \(S_{wi}\)  \(S_{wc}\) 
Initial Reservoir Pressure  \(p_{Ri}\)  7500 psia 
HydraulicFracture Dimensionless Conductivity  \(F_{CD}\)  100 
The locked reservoir properties are as follows
Property  Symbol  Default Value 

Fracture Porosity  \(\phi_f\)  0.1 
Fracture Width  \(w_f\)  0.1 ft 
Depth  \(D\)  10000 ft 
Note that depth is merely a cosmetic property as there is only one horizontal layer in the model, and no tubing calculations are made.
The dimensionless fracture conductivity is defined as
We calculate \(k_f\) from equation (3).
1.3.2 SaturationDependent Properties
1.3.2.1 TwoPhase Relative Permeability
The saturationdependent properties in this model are the relativepermeability curves. We use the Coreytype model to express twophase relative permeability. Baker's linear interpolation model is used to calculate \(k_{ro} = f(S_w,S_g)\) in threephase systems.
The Coreytype model calculates \(k_{rw} = f(S_w), k_{row} = f(S_w)\), \(k_{rog} = f(S_g)\), and \(k_{rg} = f(S_g)\) as follows:
where the following parameters must be provided
Property  Symbol  Default Value Matrix  Default Value Fracture 

Connate Water Saturation  \(S_{wc}\)  0.2  0 
Residual Oil Saturation (OW)  \(S_{orw}\)  0.2  0 
Residual Oil Saturation (OG)  \(S_{org}\)  0.2  0 
Critical Gas Saturation  \(S_{gc}\)  0.05  0 
Water Corey Exponent  \(n_w\)  2  1 
Oil Corey Exponent (OW)  \(n_{ow}\)  2  1 
Oil Corey Exponent (OG)  \(n_{og}\)  2  1 
Gas Corey Exponent  \(n_g\)  2  1 
Water Relative Permeability at \(S_{orw}\)  \(k_{rwro}\)  1  1 
Oil Relative Permeability at \(S_{wc}\)  \(k_{rocw}\)  1  1 
Gas Relative Permeability at \(S_{org}\)  \(k_{rgro}\)  1  1 
The relativepermeability values in the matrix will be open for change by the user. The fracture values are locked.
1.3.2.2 ThreePhase Relative Permeability
To achieve platform consistency between the different simulators (OPM Flow, Eclipse 100, Sensor, CMG IMEX), we have found that the Baker linear interpolation of \(k_{ro}\) is the only threephase relative permeability that is common for all simulators.
The Baker linear interpolation correlation is as follows
The testing of OPM Flow has been done by comparing OPM Flow to Sensor. The default threephase relative permeability correlation in Sensor is Stone 1 with Fayer's and Matthew's expression for \(S_{om}\). This option is not available in OPM Flow, which results in inconsistencies in predicted results for the two simulators. The default relative permeability correlation for OPM Flow (and E100) is the Baker linear interpolation. This correlation can also be activated in Sensor, making the two simulators consistent. However, Stone 1 with Fayer's and Matthew's expression for \(S_{om}\) remains the preferred choice, which was also the case for Keith Coats (see quote below)
The default threephase relativepermeability correlation for CMG IMEX is Stone 2. Hence, for consistency, our output files to CMG IMEX format will contain a keyword activating the Baker linear interpolation correlation.
The table below summarizes the threephase correlations available in the different simulators.
ThreePhase Rel.Perm. Correlation  Sensor  OPM FLow  Eclipse 100  CMG IMEX 

Stone 1  X  X  X  
Stone 1 with Fayer's and Matthew's  X  X  X  
Stone 2  X  X  X  X 
Baker  X  X  X  X 
Quote
"I would never use Baker kr. Stone 2 kr is pathological and Baker kr may be. To my knowledge Stone 1 – Fayers kr (sensor default) is the only kr logic free of pathology"
"Default kro is Stone 1 with Fayer’s modification. The ID option STONE 2 can be pathological, depending upon entered wateroil and gasoil krow and krog relperm data. The pathology is a stepfunction change from kro>0 to kro=0 for an epsilon change in saturations. The default Stone 1 is not pathological"
*Default refers to Sensor
1.3.3 MultiLayer Models
In whitson+, it is possible to have more than single layer in the model with fully or partially penetrating fracture.
The fractured layer properties are calculated by height weighted average of the \(n\) penetrated layers, as follows
1.4 Fluid Properties
whitsonʳˢ is designed to do simple and fast reservoir modeling to evaluate depletion performance. To make this tool as efficient as possible, the modified blackoil PVT model (outlined below) is used for the fluid treatment. An important feature of whitsonʳˢ is the direct link to whitsonᴾⱽᵀ. This makes the generation of blackoil PVT tables a backend feature in whitsonʳˢ, where the user only needs to provide key properties related to fluid PVT such as temperature, surfaceprocess information (separator pressure and temperature), and initial gas/oil ratio.
The key parameters of the modified blackoil PVT model are (details and definitions are found in Whitson and Brulé 2000) as follows:
Property  Symbol 

Oil Formation Volume Factor  \(B_{o}\) 
Solution Gas/Oil Ratio  \(R_s\) 
DryGas Formation Volume Factor  \(B_{gd}\) 
Solution Oil/Gas Ratio  \(r_s\) 
Oil Viscosity  \(\mu_o\) 
Gas Viscosity  \(\mu_g\) 
whitsonʳˢ will return the inplace volumes of the wellbox model when all model dimensions, reservoir properties, and fluid properties have been specified. The inplace numbers are calculated as
1.5 Water and Rock Properties
The following water and rock properties are used as default values for all reservoir simulators.
Property  Symbol  Default Value 

Water Formation Volume Factor at \(p_{ref}\) & \(T_R\)  \(B_w\)  1 RB/STB 
Water compressibility  \(c_w\)  3x10\(^{6}\) 1/psia 
Surface water density  \(\rho_w\)  62.37 lb/ft\(^3\) 
Water viscosity  \(\mu_w\)  0.5 cp 
Reference pressure for \(B_w\)  \(p_{ref}\)  14.696 psia 
Formation (rock) compressibility  \(c_r\)  4x10\(^{6}\) 1/psia 
1.6 Geomechanical Effects  PressureDependent Permeability
The matrix, EFR (if activated), and hydraulic fracture can be assigned a pressuredependent permeability on the form where \(\gamma\) is a property that be assigned in the frontend. There are two different gamma values, on for the hydraulic fracture \(\gamma_f\), and one for the matrix and EFR (if activated) \(\gamma_r\). To activate the pressuredependent permeability, toggle on "PressureDependent Permeability" under the "Well & Reservoir Data" card in the "History Matching" module.
1.7 Computing Original Gas in Place with Adsorption
Shale reservoirs possess a remarkable capacity to store methane through a process known as "adsorption." This involves gas molecules adhering to the organic material within the shale, creating a nearliquid state of condensed gas. A useful analogy is to envision a magnet sticking to a metal surface or lint clinging to a sweater, illustrating the concept of adsorption versus absorption, where one substance is trapped inside another like a sponge soaking up water. The adsorption process is reversible due to its reliance on weak attraction forces.
Compared to conventional reservoirs, shale reservoirs can store significantly more gas in the adsorbed state, even when compression is applied at pressures below 1000 psia. Notably, the gas stored in shale is predominantly held in the adsorbed state, as opposed to being free gas. This is because the volume of the cleat or fracture system within shale is relatively small compared to the overall reservoir volume. As a result, the pressure volume relationship is commonly described through the desorption isotherm alone.
Langmuir Isotherm
The release of adsorbed gas is commonly described by a pressure relationship called the Langmuir isotherm. The Langmuir adsorption isotherm assumes that the gas attaches to the surface of the shale, and covers the surface as a single layer of gas (a monolayer). At low pressures, this dense state allows greater volumes to be stored by sorption than is possible by compression.
The typical formulation of the Langmuir isotherm is:
Langmuir Volume, V_{LS}
The Langmuir volume is the maximum amount of gas that can be adsorbed by a shale at a "infinite pressure". The plot below of a Langmuir isotherm demonstrates that gas content asymptotically approaches the Langmuir volume as pressure increases to infinity.
Fig. 1: Langmuir Volume
Langmuir Pressure, p_{LS}
The Langmuir pressure, or critical desorption pressure, is the pressure at which one half of the Langmuir volume can be adsorbed. As seen in the figure below, it changes the curvature of the line and thus affects the shape of the isotherm.
Fig. 2: Langmuir Pressure
Computing Original Gas in Place (OGIP)
Ambrose et al.[1] argue that the free (OGIP) must be corrected for the adsorbed gas that occupies some of the hydrocarbon pore volume (HCPV). If porosity is estimated from core plugs where core preparation has removed the adsorbed gas, then the OGIP will be overestimated.
The convention of reporting OGIP when adsorption is included, is to report it in units of scf/ton (standard cubic feet per short ton). For a fluid system having only dry/wet gas and water, the total OGIP is
where
and with \(C_1\approx 5.7060\) and \(C_2\approx1.318\cdot10^{6}\) are unit conversion factors. The shale saturation is capped at \(1S_{wi}\), meaning the maximum amount of shale you can have is the entire HCPV. Hence, this will make the results independent of the adsorption numbers when violating this constraint.
The OGIP resolved in the ARTA and GFMB features yields the free OGIP. We can use this volume to compute the rock mass, which in turn can be used to calculate the adsorbed OGIP. Multiplying equation \eqref{eq:freegasperton} with the rock mass (\(G=\hat{G}_f m_r\)) gives the free OGIP in units of scf. Solving for the rock mass yields
This rock mass can then be multiplied by equation \eqref{eq:adsorbedgasperton} to get the adsorbed OGIP in scf, i.e.,
Total Compressibility when Adsorption is ON
In order to honor mass balance, it becomes necessary to modify total system compressibility when accounting for adsorption. Total system compressibility can be determined using the following equation:
The derivation of this equation is shown here and is contructed upon Ambrose et al.[1].
2. Numerical Model
2.1 Gridding
The numerical model will have a gridding as shown in Fig. 3. The gridding is logarithmic in the ydirection, and in the xdirection beyond the fracture tip. From the wellbore to the fracture tip the gridding is uniform in the xdirection.
Fig. 3: Gridding of Numerical Model
The gridding is calculated based on the dimensions of the model. The gridding of the numerical model will is a backend feature, where the user can choose between low, medium, and high grid refinement. Determining the number of grid blocks for the chosen level of grid refinement is discussed in the next subsection. The formulas used for gridding are
where the total number of grid blocks in the model are \(N_t = (N_{x,i} + N_{x,o})N_y = N_xN_y\).
2.2 Grid Refinement
The number of grid blocks in the model is controlled by the gridrefinement level set by the user. There are four levels of grid refinement which allow the user to decide the level of "compromise" between simulation run time and level of detail captured in the simulation results:
 Very Low (favours short simulation time)
 Low
 Medium
 High (favours capturing the details)
An internal logic is used to translate the level of grid refinement to a number of grid blocks. We use average gridblock dimensions in the \(x\) and \(y\) directions to set the refinement of the gridding. The following table contains the fixed numbers for \(\bar{\Delta} x\) and \(\bar{\Delta} y\) at the three gridrefinement levels.
Property  \(\bar{\Delta} x\) (ft)  \(\bar{\Delta} y\) (ft) 

Very Low  200  8 
Low  50  2 
Medium  25  1 
High  12.5  0.5 
The number of grid blocks in the \(x\) and \(y\) directions are then calculated as
The split of grid blocks in the \(x\)direction \(N_{x,i}/N_x\) is linear with \(x_f/x_e\). Because the number of grid blocks have to be integers, we round the numbers as follows:
The number of grid blocks outside the fracture tip are then calculated by \(N_{x,o} = N_x  N_{x,i}\).
2.3 Well Control
A well is always controlled by specifying a target rate. The most common rate is a surface product rate such as surface gas rate, surface oil rate, or surface water rate. Sometimes the target is a phase rate specified at flowing bottomhole conditions, either as a singlephase wellbore rate, or a combination of several wellbore phase rates, e.g. "reservoir" wellbore oil phase rate, "reservoir" wellbore liquid (oil phase + water phase ) phase rate or total wellbore fluid rate (gas phase + oil phase + water phase).
A zero target rate can also be specified, meaning that no production is seen at the surface. The well is "shutin" at the surface. Downhole within the wellbore, considerable flow rates usually are found during the shutin of a well with multiple well connections (to >1 grid cells). Some well connections show positive production rates and other well connections experience negative injection rates. The total injection rate of all components (surface gas, surface oil, surface water) must equal the sum of all well connection production product rates (to ensure no material balance loss in the wellbore during a shutin time step).
The well tries to satisfy the target rate during a time step, but only if the specified pressure constraint is not violated. Typical pressure constraints include flowing bottomhole pressure or flowing tubing pressure. If the pressure constraint is reached then the target rate can not be met and a rate is chosen that ensures the pressure constraint exists in the well.
During zero target rate and STOP (i.e. well closed at the surface), the pressure constraint will be honored. Layer crossflow may happen if the well is connected to multiple numerical layers with differential depletions. During zero target rate and SHUT (i.e. well closed at the sand face), the well is isolated, no well constraint is enforced.
A producing well can also have constraints on the well's produced surface product ratios GOR, WOR, OGR, surface water cut of total liquid rate or reservoir wellbore phase fractions such as water cut. These constraints are usually maximum values which, if reached, the wellbore pressure is adjusted to honor the constraint, and accordingly not allow the target rate to be met.
Well Control in whitson^{+}
There are six options currently available for controlling wells:
 Oil: the maximum surface oil production rate target or contraint.
 Gas : the maximum surface gas production rate target or constraint.
 Water : the maximum surface water production rate target or constraint.
 Liquid : the maximum surface liquid rate (water + oil) production rate target or constraint.
 BHP: the minimum bottomhole pressure taget or constraint
 BHP (synthetic): bottomhole pressure profile target or constraint, offering an option to include a maximum rate contraint. The BHP profile is set upon userdefined inputs such as (a) Minimum BHP, (b) Time to Minimum BHP, and (c) Total Simulation Days.
During the historical period, the well is controlled on the userselected option. If the well is controlled by one of the surfacerate targets, then a BHP constraint of 100 psia is applied. If the well is controlled on a BHP target, then the software will set an unrealistic rate target with the desired BHP profile as a lower limit for producers. For the synthetic BHP profile, there has, since the Sept. 2023 release of whitson^{+}, been an option of setting a rate constraint.
During shutins period, the status of the well is set to STOP, indicating no fluid production to the surface and potential layer crossflow may occur if there are any open connections.
During the forecast period, the well is consistently controlled based on either a BHP or WHP (CHP or THP) target, with the added flexibility of applying a maximum rate constraint. If none of rate constraint is chosen, then the rate will be set to an unrealistic target.
It is important to note that whitson^{+} does not support "Reservoir" or bottomholerates targets. Additionally, constraining the well based on producing ratios (GOR, WOR, OGR, GLR, etc.) is also not supported. No economic criteria for production wells are currently being enforced. This implies that the well will only be shutin (denoted by a well status of STOP) if the well has zero rates and/or the pressure falls below the predefined minimum bottomhole pressure (whitson^{+} default is 100 psia).
Should your company require any of these constraint options, please feel free to contact us or use the Send Feedback button at the top right of whitson^{+}.
How does whitson^{+} recognize a shutin?
If Well Control is ,
 Oil: well is shutin if oil rate for time step is 0.
 Gas: well is shutin if gas rate for time step is 0.
 Water : well is shutin if water rate for time step is 0.
 Liquid : well is shutin if liquid rate (oil + water) for time step is 0.
 BHP: well is shutin if all rates (oil + gas + water) for time step is 0.
PS! Should all the rates for all time steps be uniformly zero, the software assumes that the user would like to control the model on a custom, uploaded BHP profile. Consequently, the well won't be considered as shutin.
3. What affects producing GOR profile and how?
3.1 Reservoir Parameters
We observe a widerange of producing GOR profiles in tight unconventionals. But which reservoir parameters change producing GOR when the bottomhole pressure is below the saturation pressure? And how?
Fig. 4: Effects of Reservoir Parameters on producing GOR
There are four broad categories of reservoir parameters that impact producing GOR
 Fracture conductivity
 Geomechanical effect
 Contributing volume beyond frac tips / height ("unstimulated rock")
 Relative permeability
This is summarized in the table above. Read more in this slide deck →→here←←.
3.2 Flow Regime
Fig. 5: Effects of flow regime and boundaries on producing GOR profile.
In addition to the reservoir parameters outlined above, flow regime matters. Everything else being equal producing GOR follows this relationship below saturation pressure:
Where IA refers to infinite acting, TF refers to transitional flow and BDF refers to boundary dominated flow.
4. Forecast
4.1 BHP Forecast Type
Fig. 6: Forecasting the numerical model
After one has finished matching the production history, the results can be used to forecast the production. In whitson+, the simulation forecast is controlled by pressure (i.e. BHP or WHP). Whitson+ provides two alternatives in creating BHP forecast:

Constant BHP throughout simulation time.

Parametric Decline to set the forecasting function between the input initial and final BHP as hyperbolic or linear based on the input Nominal Pressure Decline in %/yr or psia/d, respectively.

Custom BHP provided by user in tablewise manner. By default, the first row will be the BHP at the last time step of history matching simulation run.
Forecasting Using WHP
More on WHP forecasting coming in this section soon!
5. Autofit
5.1 How to run Autofit
It is possible to automatically history match a numerical model to historical data with the "Autofit" option. To run the "Autofit",
 Click the AUTOFITbutton to the upper left.
 Select parameters you would like to alter (\(x_f\) and \(k\) are defaults).
 Provide their ranges for each parameter.
 Turn on/off weighting for relevant timeseries (global weighting switch).
 Click AUTOFIT to the lower right.
5.2 Set Weighting Points
In addition to the global weighting switch, you can set different weighting throughout time.
Each timestep can be set to either
 Weight = 1 (default)
 Weight = 10 (considered important data to match)
 Weight = 0 (considered less "bad data")
5.3 Weighting & Error Calculation
whitson+ facilitates regression on parameters to minimize the sum of squares of weighted residuals in the context of observed data and corresponding predictions. The objective function in is defined as:
Here, represents the RMS residual error between observed measurement and its corresponding prediction, while is the userassigned weighting factor for that residual. Ideally, these weighting factors should be inversely proportional to the standard deviations of the residuals. Minimizing in provides the maximum likelihood estimation of the model parameters, assuming independent and normally distributed residuals.
The default weighting factors are 1, while they can also be set to 0 (no importance), and 10 (high importance). Reported is a rootmeansquare (RMS) residual error () defined as:
This metric is related to but is more easily interpreted.
The residuals are calculated as relative percentages using the formula:
where is the historical production value, is the simulated value, and is the reference value for observation . Historical production values that is 0 is excluded from the error calculation, while historical production values that has a higher error than 100% is set to 100% error (to avoid that individual outliers dominate the whole error calculation)
5.4 Optimization Algorithm: NelderMead
NelderMead Simplex is a modified version of the downhill Simplex method used for nonlinear regression in reservoir simulation. It enhances convergence by applying constraints on parameters during the search, checking estimates against preset limits. After convergence, the routine restarts with a slight perturbation to avoid local minima, ensuring a global minimum. While it may require more function evaluations, its simplicity and lack of derivative requirements make it robust under various conditions.
5.5 Optimization Heuristics
In pursuit of practicality, our approach involves modeling 50 historical data points in each iteration, ensuring fast run times of each individual numerical run. he objective of the minimization function is to align production cumulatives (considering weighting) and bottomhole pressures over time.
6. Refrac Analysis
A ‘refrac’ or refracturing treatment refers to the process of performing another hydraulic fracturing operation on a well that has already undergone the process in the past.
Typically, you may want to evaluate a refrac that was recently performed on one of the wells. To do that, use the MFMB feature  The workflow is outlined in this section here.
Sometimes, you also want to be able to predict the initial rates and subsequent decline for a refrac that is planned on a well in the future. To do that, use numerical modeling as shown below.
6.1 Workflow for estimating uplift on planned refracs
This section outlines our recommendation for what the simplified workflow could be 
 History match the target well and resolve a numerical model with high confidence.
Then forecast this model to generate a baseline future expectation from the well, under current conditions as shown above in section 4 above.
Fig. 7: High confidence history matched numerical model

Copy the well (Wells page > Click the three dots to the right of the well name > Copy well)

On the copy of the well, modify the numerical model to reflect the expectation from the refrac, i.e. 
a. Change the completion  increase the number of fractures (Nf) or fracture halflength (xf) or both.
b. Increasing Nf will primarily accelerate depletion but yeild a higher intial rate.
c. Increasing xf (if we expect a bigger refrac completion) will add additional pore volume
d. Forecast this updated model.
Fig. 8: Updating the baseline numerical model

Delta between outcomes of Steps 1 and 2 is the refrac potential of this well.

Superpose these resulting rates at the time of refrac to get the uplift expected after refrac
Fig. 9: Forecasting the updated model and evaluating uplift
Let us know, if you'd like us to add another export option for your company. Contact us or use the Send Feedback button at the top right of whitson+.
References
[2] Open Porous Media. 2023. OPM Flow Reference Manual, Version 202304. Oslo, Norway: OPM.