University of Southern California Dissertations and Theses

Integration of multi-physics data into dynamic reservoir model with ensemble Kalman filter

(USC Thesis Other)

Integration of multi-physics data into dynamic reservoir model with ensemble Kalman filter

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content Integration of Multi-Physics Data into Dynamic Reservoir Model with
Ensemble Kalman Filter

by
Siavash Hakim Elahi

B.Sc., Chemical Engineering
University of Tehran, 2010

M.S., Petroleum Engineering
University of Tehran, 2012

Ph.D. Dissertation

Doctor of Philosophy
(Petroleum Engineering)

Faculty of the USC Graduate School
University of Southern California
May 2018

Dissertation Committee

Dr. Behnam Jafarpour
Dr. Iraj Ershaghi
Dr. Charles Sammis
ii
© Copyright by Siavash Hakim Elahi, May 2018
All Rights Reserved

3
Table of Contents
Abstract ................................................................................................................................5
Introduction ..........................................................................................................................7
Chapter 1: Literature Review .............................................................................................10
1.1 Characterization of Hydraulic Fracturing in Unconventional Reservoirs ......11
1.1.1 Fracture Models
1.1.2 Data Types
1.1.3 Uncertain Reservoir and Fracture Properties
1.2 Continuous Model Calibration Methods…………………………………… 21
1.2.1 Deterministic Modeling and Gradient Based Methods
1.2.2 Stochastic Modeling and Ensemble Based Methods
1.3. Description of Rock Facies Distributions Using Sparse Data………………25
1.3.1 Direct Conditional Facies Simulation
1.3.2 Discrete Regularization
1.3.3 Continuous Parameterization
Chapter 2: Dynamic Fracture Characterization………………………….………………..33
2.1 Fracture Models and Parameters……………………….……………………..33
2.1.1 Planar Model
2.1.2 Discrete Fracture Network
2.2 Data Assimilation with Ensemble Kalman Filter …………………………….39
2.3 Results and Discussion ……………………………………………………….43
Chapter 3: Continuous Parameterization of Discrete Geological Facies………………….52
3.1 Problem Description…………………………………………………………….53
3.2 Data Assimilation Technique…………………………………………………...55
3.2.1 Distance Transform
3.2.2 SVD-Parameterized Distance Transform
3.2.3 Implementation Steps
3.3 Experiments and Results…………………………………………………..........62
3.3.1 Estimation of Facies Distribution

4
3.3.2 Estimation of Facies Distribution and Properties
3.4 Discussion and Future Work ..............................................................................73
Chapter 4: Integration of Discrete Microseismic Data .......................................................76
4.1 Coupled Flow ad Geomechanical Model ...........................................................77
4.2 Failure Criterion ................................................................................................79
4.3 Source of Uncertainty in Potential Microseismic Events ...................................80
4.4 Microseismic Data Simulation ...........................................................................81
4.4.1 Deterministic MEQ Simulation Model
4.4.2 Stochastic MEQ Simulation Model
4.5 Sensitivity Analysis ............................................................................................84
4.5.1 Deterministic Approach
4.5.2 Stochastic Approach
4.6 Proposed Data Assimilation Techniques ............................................................86
Kernel Density Function
4.7 Results ................................................................................................................87
4.7.1 Estimation of Facies Distribution
4.7.2 Estimation of Facies Distribution and Rock Properties
4.8 Discussion and Future Work ..............................................................................90
Chapter 5: Summary and Conclusion .................................................................................94
Appendix A: Flow Regimes in Production from Fractured Models ..................................99
Appendix B: Nomenclature .............................................................................................100
References ........................................................................................................................104
Tables ...............................................................................................................................122
Figures..............................................................................................................................128

5
ABSTRACT
Construction of predictive subsurface flow models involves subjective interpretation and
interpolation of spatially limited data, often using imperfect modeling assumptions. Hence, the
process can introduce significant uncertainty and bias in predicting the flow and transport
behavior of these systems. The uncertainty in the distribution of discrete fracture network and
rock facies distribution in complex geologic environments along with their flow and mechanical
properties can be consequential for forecasting the dynamic response of these systems to
perturbations due to fluid injection and development activities. In this work, tracer test,
performance and microseismic data are used for dynamic characterization of important
parameters of dynamic reservoir models, including (1) matrix fluid and mechanical properties,
(2) rock facies distribution (3) planar fracture half-length and hydraulic conductivity, (4) discrete
fracture networks density and conductivity, and (5) fracture closing (conductivity decline) rate
during production. Conventional model calibration techniques are designed to update continuous
model parameters and have difficulty in estimating discrete parameters such as rock facies
distribution (or distribution of discrete fracture network). Here, we present a distance transform
approach for converting discrete parameters to continuous parameters that can be updated, using
continuous model calibration methods, from flow and transport-related data. Distance-based
transforms are widely used in discrete image processing, where the discrete values in each pixel
are replaced with their distance (i.e., a continuous variable) to the nearest cell with a different
value (i.e., nearest boundary cell). The resulting continuous distance map can then be back-
transformed to retrieve the original discrete image. For facies model calibration, the facies maps
are first transformed into continuous spatial maps of distances to facilitate model calibration with
continuous inversion methods. For facies with large-scale connectivity, truncated singular value
decomposition (SVD) parametrization can be used to represent the distance maps with low rank

6
parameters. In this work, a variant of the ensemble smoother is used to update the final
continuous parameters (either distance maps or their SVD coefficients if used). The resulting
updated parameters are back-transformed to obtain calibrated models. The distance transform
method addresses an important problem in facies model calibration where model updating can
result in losing facies connectivity and discreteness.
In general, high pressure fluid injection into geologic formations can trigger rock deformation
and failure and lead to microseismic events. Injection induced seismicity has been proposed as a
monitoring technique that can be used to constrain the description of rock flow and mechanical
properties. To establish the correlations between rock properties and microseismic clouds, the
developed framework combines coupled flow and geomechanics simulation outputs with Mohr-
Coulomb’s failure criterion to describe the spatiotemporal distribution of seismicity potential in
the formation. This model is then used with an ensemble smoother with multiple data
assimilation, ES-MDA to estimate rock properties form observed microseismic data. Since the
ensemble smoother is designed for assimilation of continuous data, kernel density estimation
(KDE) is applied to convert discrete microseismic events to a continuous map of seismicity
density. We present the developed formulation and discuss its application to inversion of
microseismic data for characterization of reservoir flow and mechanical properties. Several
examples are presented to evaluate the performance of the method for estimation of rock
properties from microseismic data.

7
INTRODUCTION
This work focuses mainly on the subsurface model calibration and estimation of dynamic
reservoir properties including matrix and fracture properties using different types of monitoring
and performance data. The usual method to improve the productivity in the geothermal and
unconventional reservoir is to stimulate reservoir with hydraulic fracturing. There are several
monitoring tools and data for characterization of fractured reservoir. Microseismic, tracer and
production data are primary data for this purpose. In general, microseismic data have limited
signal to noise ratio and provide information about fracture location and geometric attribute,
which may not be sufficient for understanding fracture connectivity, flow properties, and
contribution to production or SRV. Tracer tests are useful tools that have several applications,
including understanding reservoir connectivity, verifying flow barriers, and characterizing flood
patterns and inferring reservoir heterogeneity. For fracture characterization, tracer test is
performed using a single well, by injecting fracturing fluid with tracers and measuring its
concentration in the flow back. Tracer tests have advantages over other monitoring tools such as
production logging. It can be used to infer fracture transport properties, the contribution of each
stage to production, optimize production, and improve well completion design. Tracer test
response data for different fracturing stages can be used to potentially estimate fracture
conductivity, length, and other geometric attributes along with rock matrix properties. In
addition to tracer test, production data can be included in characterization of dynamic reservoirs.
Subsurface flow model calibration involves solving an inverse problem in which the mismatch
between observed and model-predicted data is judiciously reduced. Physical rock and fluid
properties, initial conditions, boundary conditions, sources or sinks terms and a mixture of the
above are examples of unknown parameters. While the model calibration process is used to

8
ensure that a given model reproduces the historical measurements in the field, the main objective
is to improve the model to accurately represent the flow and transport processes, and hence
predict the future performance of the system. There are two main subsurface flow model
calibration: Deterministic and stochastic modeling. In this work, we mainly focus on stochastic
modeling.
Fracture complexity can vary significantly depending on several factors including geological
condition, preexisting natural fracture, and anisotropy in principle stress direction. There are two
main models for fractures, planar model, and complex discrete fracture network model. In this
research, we are interested to characterize both. In the first set of our experiments, we used a
simplified subsurface model, and we integrated tracer and production data to characterize
fractured reservoir by implementing variation of the Ensemble Kalman Filter (EnKF). And for
estimation of discrete parameters such as rock facies distribution in the matrix and discrete
distribution of fracture elements in the reservoir, we introduced distance transform to relax the
limitation of conventional history matching techniques such as EnKF for integration of non-
Gaussian parameters.
To have more complex and representative model for a subsurface reservoir, a coupled flow and
geomechanical simulation model was built using CMG-GEM module. The purpose is to have a
more realistic model, and to estimate hydraulic and geomechanical properties of the rock such as
Young’s modulus, Poisson’s ratio, friction angle, cohesion, and permeability by integrating
discrete microseismic data. Therefore, we have combined the complex coupled flow and
geomechanical simulation model with Mohr-Coulomb rock mechanical failure to predict
potential rock failure map and eventually microseismic cloud. However, Microseismic data are
discrete data and cannot be used as an input to conventional data assimilation techniques such as

9
EnKF and we need to modify our data assimilation technique to be adaptive to discrete data.
Therefore, WEpropose a workflow for integration of discrete micro-seismic data to estimate rock
flow and mechanical properties. This workflow includes the Kernel density estimation (KDE) for
transformation of discrete micro-seismic data to continuous observation
There are several important aspects of proposed workflow and method to characterize fractured
reservoirs. 1). Integrating different types of data such as tracer, production and microseismic, 2).
Having more representative model for subsurface reservoir by coupling flow and geomechanical
model, 3). Addressing the existing issues of conventional data assimilation techniques, extending
the application of Ensemble Kalman Filter for non-Gaussian reservoir parameters (Discrete
facies map) and introducing techniques for integration of discrete microseismic data.

10
CHAPTER 1
Literature Review
The growth in the implementation of optimization techniques in reservoir management has
increased demands on the application of history matching to make models that not only
reproduce the historical production behavior, but also preserve the geological structure and
assess forecast uncertainty. The inverse theory is the fine art of estimating model parameters from
data. In conventional reservoirs, the most important model parameters are rock flow properties
such as permeability and porosity of the field. These properties mainly control flow performance
of the subsurface reservoir. However, there are other parameters, geomechanical rock properties
that play a key role in risk assessment and mitigation of environmental issues associated with
high-pressure fluid injection into the reservoir. There are also different types of data such as
production rate, chemical and radioactive tracer, pressure, the temperature which can be used to
characterize subsurface reservoirs.

Researchers have used various history matching techniques to estimate reservoir flow properties
such as permeability and porosity in the conventional reservoirs. However, there is a limited
number of publications on development and implementation of stochastic data assimilation
techniques for characterization of unconventional reservoirs and simultaneous estimation of flow
and geomechanical properties. People have used traditional manual history matching long time
ago to estimate reservoir parameters. In one of the real case studies, it has been reported that
manual history matching of the bottom-hole pressure of producer was partially accomplished by
modifying the transmissibility within a zone surrounding the well, which results in a geologically
unrealistic permeability field. Also, the water cut history of producer, and the fluctuating gas oil
ratio data cannot be easily matched using manual history matching [Zhang et al., 2009].

11
Moreover, due to the complex relationship between data and reservoir parameters, it is generally
difficult to achieve well-by-well manual history matches that are geologically feasible. In
general, there are two main automated history matching techniques to integrate data,
deterministic technique and stochastic technique. Both of them will be discussed briefly in the
upcoming sections.

There are a numerous amount of field examples in the literature in which people implemented
various data assimilation techniques to diagnose the reservoir heterogeneity, flow barriers, and
reservoir connectivity. Evensen introduced a thorough workflow to update reservoir parameters
in North Sea field [Seiler et al., 2010]. They showed that updating relative permeability
properties by production data assimilation allows accountig for geological heterogeneity effects.
Haugen et al. presented a successful study for a North Sea field case, where real production data
have been integrated using stochastic technique [Haugen et al., 2006]. They improved the
estimated model parameters (permeability and porosity) by assimilating data such as bottom-hole
pressure and production rate of oil, gas and water. Oliver et al. presented a field example in their
paper and compared the application and feasibility of manual and automated history matching
technique in their work [Zhang et al., 2009]. In other work, Al Reynolds et al. implemented one
of the well-known history matching techniques, Ensemble Kalman Filter (EnKF) with
covariance localization to history match permeability field by assimilating production data
[Emerick et al., 2011].

1.1 Characterization of Hydraulic Fracturing in Unconventional Reservoirs
The efficiency of hydraulic fracturing combined with multi-well pad horizontal drilling has
caused a rapid growth in production from shale oil and gas reservoirs in the past decade [Boyer,

12
2006; Drake, 2007; King, 2010]. Hydraulic fracturing in tight reservoirs in order to increase the
permeability has become a standard industry practice. With fractures initiating and growing
thousands of feet below the surface, understanding fracture distribution, geometry, and
production response are quite challenging. Moreover, based on the field experience, hydraulic
fracturing is proving to be vastly more complex and unpredictable than initially thought. In some
cases, different production performance has been reported for identical treatments within the
same formation. While fracture treatments continue to be designed using the best tools and
techniques available, characterization of fracture geometry and hydraulic properties is a
nontrivial task. Numerous fracture diagnostic techniques are proposed to achieve this goal and to
improve our understanding of hydraulic fracture behavior [Cipolla and Wright 2000].

Characterization of fracture distribution and growth is one of the main challenges in hydraulic
fracturing. Fracture complexity can vary significantly depending on some factors, including
geologic conditions, anisotropy in the principal stress directions, and pre-existing natural
fractures. Several fracture models have been proposed in the literatures, including planar model,
wire mesh model, and complex fracture network model [Crockett, 1989; Xu, 2009; Cipolla,
2010]. The planar growth of hydraulic fracturing is the primary model that has been used to
describe hydraulic fracturing. The planar model is valid when the difference between the
maximum and minimum horizontal stress in the field is large [Cipolla et al., 2011]. When natural
fractures exist in the field, they can give rise to the development of complex fracture networks,
which are more difficult to characterize, usually requiring stochastic representations to account
for uncertainty in fracture distribution and connectivity.

13
1.1.1 Fracture Models
Fractured reservoirs can be divided into four different types depending on their permeability and
porosity [Nelson, 2001]: Type-I, fractures are the main reservoir porosity and permeability;
Type-II, fractures provide the main permeability while the matrix has the main porosity; Type-
III, matrix permeability is high, so fractures only increase flow capacity; and Type-IV, fractures
are filled with minerals and provide no additional porosity or permeability. Fracture
characterization is a critical step in developing a fracture reservoir simulation model, in
particular parameters such as fracture length, orientation, spacing, porosity, connectivity,
aperture, and permeability can have a detrimental effect in describing fracture models and
predicting their flow behavior.
Two general approaches for simulating flow in fractured systems involve discrete fracture model
(DFM), also known as discrete fracture network (DFN) model, and an equivalent continuum
model. In the DFM approach, flow is explicitly modeled in each fracture and in the matrix using
Darcy’s law, which can be computationally very expensive if fine-scale geological models are
used. Alternatively, in the equivalent continuum approach, cell properties are assigned equivalent
properties that represent the combined effects of fractures and matrix. Dual-porosity and dual-
porosity/dual-permeability models are examples of the continuum modeling approach and are
generally more approximate and computationally more efficient than the DFM, which is why
they are applied more widely to fractured reservoir simulations.
Representing and simulating the physics involved in the production from hydraulically fractured
reservoirs is complex and typically involves coupled flow and geomechanical simulation.
Moreover, capturing the dynamic deformation and fracture closing behavior during production

14
from hydraulically fractured reservoirs requires coupled flow and geomechanical models,
although more simplified approaches have been proposed in the literature.
A simple way of predicting production from hydraulically fractured reservoirs is to represent
fracture effects by adjusting the hydraulic properties of the cells that contain fractures. Appendix
B provides a brief discussion on how these adjustments are performed. This modeling approach
has been adopted and studied in several previous publications [Cox, et al. 2008; Lolon and
Cipolla, 2009; Cipolla et al., 2010; Li et al., 2011; Yu et al. 2014]. In the Appendix A of their
paper, Cipolla et al. (2010) discuss the validity of this modeling approach even when relatively
large cell sizes (~ 2ft) are used to represent fractures. Li et al. (2011) also use this modeling
approach and show that with proper adjustment, larger cell sizes (up to 10 ft in their examples)
can be used as long as proper adjustments are made. In our study, we use this simplified
approach with adjusted fracture properties to model and predict production from fractured
reservoirs. Elahi and Jafarpour (2017) showed that predictions generated using logarithmic grid
refinement with relatively large values of the fracture cell width, and properly adjusted hydraulic
properties for fracture cells, were quite similar.
While accurate simulation and predictive models are generally desirable, in practice
computational constraints motivate the use of more simplified simulation techniques. In
particular, in ensemble-based forecasting for data assimilation and uncertainty quantification,
where a large number of simulations with different input parameters are performed to explore the
uncertainty space, fast proxy models provide an attractive option, motivating the need for the
development of more efficient models that capture the general behavior of an exact simulation
model. In this context, there is a ubiquitous trade-off between using a few realizations with more
accurate simulations or a large number of simulation runs with fast approximate simulation

15
techniques. In general, for large-dimensional parameter spaces, an extensive search of the
parameter space is essential and is only feasible using proxy models.
For ensemble data assimilation, a detailed model of flow and transport, and potentially coupled
with geomechanics to include dynamic effects related to fracture closing, could be deemed
computationally too demanding for this application (Tarrahi et al., 2015). In case of the EnKF,
the simulations are performed to compute sample statistics, including covariances and Kalman
gain matrices, which are often considered to be approximate due to limited sample sizes (Tarrahi
et al., 2016).

1.1.2 Data Types
There are various diagnostics tools and methods such as well testing, fracture modeling, borehole
image logging, production logging, surface and downhole tilt fracture mapping, micrsoeismic
mapping, tracer test and production data analysis in the literature that people have used in real
cases to diagnose the fracture and characterize its properties. Combination of these tools and
method can be used to verify the estimation of the fracture properties from fracture models and
improve fracture treatment design.
Sudden energy release due to mechanical rock failure propagates outward through elastic waves
and resembles a small-scale earthquake event that can be recorded by geophones, and
subsequently processed and analyzed to locate the event source and estimate its intensity. The
resulting distribution of events in space forms a cloud of micro-earthquake (MEQ) or
microseismic data that carry important information about physical rock properties and possible
fracturing mechanisms and fracture propagation behavior [King, 2010; Cipolla, 2010]. However,
several challenges exist that can affect the reliability of microseismic data. In particular, the

16
seismic signal received by sensor arrays is typically very weak and has a low signal-to-noise
ratio, introducing significant uncertainty about the location of the events, typically calculated by
time separation of P-wave and S-wave [Peyret et al. 2012]. Furthermore, microseismic data is
mainly used for imaging fracture distribution and geometric attributes, which is not sufficient to
distinguish propped open fractures that contribute to production from closed (non-conductive)
fractures. Figure 3-a demonstrates this concept, where for Stage 3, the MEQ data overestimate
fracture length (the real fracture length is shown with rectangle). Also, Figure 3-b shows how
MEQ data can mislead geometry interpretation of DFN (Brown rectangle represents conductive
fracture, white one indicates nonconductive natural fracture). On the other hand, flow and
transport related data, such as tracer test and production data, may be used to infer open
(conductive) fracture length as well as fracture hydraulic conductivity [Ghods and Zhang, 2012;
Moreno et al., 2014; Elahi and Jafarpour, 2015].
Flow back and fracture load recovery reveal useful information about hydraulic fracturing. By
adding chemical and radioactive isotope tracers to fracturing fluid throughout treatment, we can
trace fluid flow back [King 2010]. Tracers have been used in the oil field for several decades.
Reservoir characterization including the identification of trends of fluid movement, layering,
flow barriers and saturations is the main purpose of a tracer test [Zemel, 1995]. In a single-well
test, the tracer is injected and produced from the same well. Tracer concentration is measured in
the fluid flow back. Residual oil saturation, heterogeneity [Descant et al., 1989] and wettability
[Ferreira et al., 1992] can be estimated by using a single-well tracer test.
Tracer tests are performed to obtain important information about the stimulated formation,
including evaluating the performance of the fracturing treatment (that is, the ability of the
stimulated formation to return fluids), understanding fluid breakout (that is, measuring flowback

17
efficiency of each fracture stage directly), and marking fracture points [Silber, 2003; Woodroof,
2003]. In some cases, fracturing treatments are not able to recover the injected fluids and tracers
completely, partly because of the reduction in fracture conductivity and fracture closing [Veach,
1982; Barree, 1995]. Hence, the amount of recovered tracer can be related to fracturing
performance, indicating that the tracer tests can be used to estimate fracture geometry and
conductivity. Since tracer tests are performed during a short period of time, they tend to reveal
fractures short-lived (or new well) flow and transport response. Once the production phase
begins, additional performance data (e.g., production flow rate) become available and can be
used to infer fracture and matrix properties.
There are other conventional and advanced technologies to better understand unconventional
reservoir, hydraulic fracturing and improve treatment design. One of them is surface and
downhole tilt fracture mapping; it has a very simple principle in fracture mapping. Created
hydraulic fracture causes characteristic deformation of the rock surrounding the fracture. Then,
by measuring hydraulic fracture tilt of the earth at the surface and downhole using extremely
accurate carpenter’s level, we can estimate fracture orientation and geometry.
Other techniques such as radioactive trace technology, temperature logging, production logging,
borehole image logging, downhole video and caliper logging are known as direct near-wellbore
fracture diagnostic techniques. Not being able to provide information for the fracture growing
more than 1-2 feet from the wellbore is the main limitation of these techniques. Radioactive
tracers are added to the proppant and fracturing fluid. After the fracture treatment, gamma ray
logging is used to measure radioactivity within the neighborhood of the wellbore. Temperature
logging can be used to measure the cooling effect caused by fracturing fluid during the treatment.
Then, we can look for irregular heat up effect as the formation returns to thermal equilibrium.

18
Intervals with the most fracturing fluids show the slowest heat up. One of the most popular
production logging used for fracture diagnostics is noise logging which uses the created sound by
fluid entry into the wellbore to identify fluid entry point. The other production logging is spinner
survey which is the combination of several sensors including flowmeter (spinner), temperature,
pressure, fluid density and the gamma ray to allow a comprehensive evaluation of amount and
type of fluid entry into the wellbore from each perforated zone. Image borehole logging provides
orientation of induced and natural fractures in wellbore proximity. Interpretation of fracture
orientation can be related to the direction of maximum stress and fracture azimuth. The
downhole video provides information about the zones contributing to the production. Caliper
logs which are performed before fracture treatment. They can measure the orientation of
wellbore breakout and ellipticity which can be related to the maximum principle stress.
Pressure transient test (Well test), and production data analysis are indirect fracture diagnostic
techniques. They can be used to estimate fracture dimension and conductivity by using indirect
measurements of pressure response during fracture treatment and pressure and flowrate during
the production.
Fiber optic distributed temperature sensing (DTS) monitoring is a real-time monitoring tool for
hydraulic fracture characterization. It has many application including estimation of fracture
initiation depth, vertical coverage and number of generated fractures and characterizing the
undesired flow behind casing. However, the desired application of this tool in our study is to
assess the effectiveness of the stimulated treatment and fracturing fluid placement. Fiber
conveyed with coil tubing inside the flow path or cemented behind casing. The location of fiber
has a major effect on the recorded temperature depending on the thermal conductivity between
the fiber cable and injected fluid path. People also have recommended using induced hot/cold

19
thermal tracers along with temperature modeling to calculate flowrate distribution along the
perforated intervals. DTS can be run during the actual hydraulic fracture treatment and also
during any shut-in periods. Contrary to the temperature logs which give a snapshot of
temperature during shut-in, DTS generates a complete dynamic temperature profile. Temperate
logs run after a pre-treatment and multiple wireline temperature logs are run over time, and it
creates a temporal temperature profile of each logged depth in contrast to DTS which gives us a
complete wellbore temperature profile at the same time [Sierra et al., 2008]. Therefore, DTS can
be used to quantify the flow distribution over the perforated intervals during hydraulic fracture
treatment.
There are several case histories in the literature that people used the combination of the above
technologies to characterize fracture properties. Although it is ideal to apply all the technologies
on fracture treatment, it is not neither practical nor economic [Cipolla, 2000]. By making some
assumption about fracturing process, we can make a reservoir and fracture model. This model
later can be used to match observed data on the field. However, the associated uncertainty with
reservoir and fracture properties will complicate this process. In the next section, we discuss the
possible uncertainty in the reservoir properties, and after that, we present an overview of data
assimilation techniques which can be used to reduce the uncertainty in the reservoir properties.

1.1.3 Uncertain Reservoir and Fracture Properties
Lack of complete knowledge about the factors contributing to the uncertainty or the way to
quantify them can cause uncertainty underestimation. In the oil industry, optimum sampling and
observations are always limited which results in high degree of uncertainty when building the

20
reservoir model. Stochastic modeling is a useful approach for quantifying uncertainty by
generating several reservoir property models (See next section).
In conventional reservoir, the most important uncertain parameters are permeability and porosity
distribution and also reservoir connectivity in complex channel systems. Researchers have
implemented different data assimilation techniques to estimate these parameters and reduce the
associated uncertainty [Jafarpour et al., 2009; Jafarpour et al., 2011; Khaninezhad et al., 2012].
However, in unconventional reservoirs, due to the existence of induced hydraulic fracture, there
are additional fracture properties and their associated uncertainty. The most important properties
of unconventional reservoirs can be categorized into fracture geometry, fluid and transport
fracture properties, and geomechanical properties of matrix and fracture.
Fracture complexity can vary significantly depending on several factors including geological
conditions, stress anisotropy, and preexisting natural fractures. There are two main fracture
models in the literature, planar model and complex discrete fracture network model. Depending
on which model is suitable for a particular reservoir, description of fracture geometry and
corresponding properties can be different. The planar model can be described deterministically
and the most important and uncertain fracture parameters are fracture half-length, height,
conductivity, azimuth, width. On the other hand, for DFN, the complexity of the fracture
networks makes it difficult to describe those using deterministic techniques, and a stochastic
description is more appropriate. In this case, some of the important parameters of the fracture
network include regional fracture density, the fraction of open fractures that contribute to
production and the estimated SRV, and effective fracture conductivity.
By improving the descriptive model of the reservoir and including the geomechanical properties
of the reservoir, we are adding more source of uncertainty in our data assimilation. These

21
additional uncertain properties are geomechanical properties such as Young’s modulus,
Poisson’s ratio cohesion and friction angle. Initial stress state of the reservoir or initial stress
state applied on the preexisting natural fracture plane which can result in different interaction
between induced fracture and natural fractures such as crossing natural fracture, or opening
nature fractures are other possible source of uncertainty.

1.2 Continuous Model Calibration Methods
Subsurface flow model calibration involves solving an inverse problem in which the mismatch
between observed and model-predicted data is judiciously reduced. The unknown parameters can
include physical rock and fluid properties, initial conditions, boundary conditions, sources or
sinks terms and a mixture of the above. While the model calibration process is used to ensure
that a given model reproduces the historical measurements in the field, the main objective is to
improve the model’s ability to accurately represent the flow and transport processes, and hence
predict the future performance of the system and mitigate its associated environmental issues.
There are two main subsurface flow model calibration which will be discussed with more details
in the upcoming sections: Deterministic and stochastic modeling.

1.2.1 Deterministic Modeling and Gradient-Based Methods
Subsurface model calibration is a nonlinear inverse problem. The assumption of the perfect
dynamical model and no model errors lead to a cost function which is the distance between the
model solution and the observation. In this problem, we are trying to find model parameters that
minimize the squared norm of the data mismatch (cost function)
 
2
obs
) (
d
m m J
C
g d   (1)

22
where m is model variables, g(m) is the output of the model (i.e. the fluid flow simulation results
such as pressure and rate) and
obs
d is the vector of observed data. The norm in the above
equation is the weighted l2 norm. A common choice to weight the squared data mismatch in
automatic subsurface model calibration is to use the inversion of noise covariance matrix.
      ) ( ) (
obs obs
m m m J
d
T
g d C g d    (2)

where C
d
is the covariance of noise in observed data. The most important difficulties associated
with this optimization problems are the problems of non-uniqueness. Non-uniqueness is present
when more than one set of parameters leads to minima of the objective function.
When the number of model parameters is large (greater than the number of independent data),
the problem of minimization is ill-posed. To overcome this problem, one approach is to reduce
the number of parameters to be estimated without losing the ability to reproduce spatial
variability. This approach is called parameterization. Another approach is to add more
constraints by adding a regularization term to the objective function. More details can be found
in literatures [Tikhonov and Arsenin, 1977; Carrera et al., 2005; Hunt et al., 2006; Oliver et al.,
2008].

There are different types of descent methods which can be used to solve the above minimization
problem. Typical methods involve the calculation of gradient of predicted observation (g(m))
concerning reservoir parameters (m). The gradient may be evaluated by solving the adjoint
model [Evensen 2006]. The solution obtained is a single realization which hopefully will be the
global minimum of the J(m) as a cost function. However, the probability of converging to a local

23
minimum is high and in most realistic applications the global minimum is never found. With a
single realization as a solution, there is no error estimate available. Since we are focusing mainly
on stochastic modeling and uncertainty assessment of the reservoir parameters in this work, we
will not further discuss this type of model calibrations and gradient based methods here.

1.2.2 Stochastic Modeling and Ensemble Based Methods
Stochastic modeling is a useful approach for quantifying uncertainty by generating multiple
reservoir property models. Data assimilation techniques can be formulated by using Bayes’
theorem, which states that the posterior pdf of a random variable conditioned on data is
proportional to the prior pdf for the random variable multiplies with the likelihood function
which describes the uncertainty in the data.

An alternative to the proposed standard minimization problem in the last section is sampling
methods which applies Monte Carlo sampling of the solution space. Sampling methods increase
the likelihood of finding the global minimum, but it causes a high numerical cost associated with
a large number of simulations when the dimension of the parameter space is large. If the global
minimum can be found, the sampling methods, such as genetic algorithms [Goldberg, 1989], also
provide an estimate of the posterior pdf around the global minimum with an error estimate.

Evensen proposed an alternative method for sequential processing of measurement starting from
Baye’s theorem with two assumptions. The first assumption is that the simulator is Markov
model which means that the solution at the current time is only dependent on the solution at the
previous time instant. Secondly, measurement errors are uncorrelated in time. In another word,
error covariance matrix is diagonal. The assumption of independent data in time allows for

24
Bayes’ theorem to be written recursively where data are processed sequentially in time, and it
results in a sequence of inverse problems. Each of the inverse problems can be solved using the
desired minimization method. However, because the process is recursive, we require having
access to both the prior estimate and its uncertainty.

The particle filter method is proposed for solving the sequence of inverse problems [Doucet et
al., 2001]. In this method, pdf for the solution is approximated by a large number of particles
(model states). Each of the particles is integrated in time based on the model equations with
stochastic terms to represent model errors. At each time, when data are available a Bayesian
update is computed by re-sampling of the posterior distribution. A solution is a large number of
realizations which represents the posterior pdf. A large number of realizations is required for the
particle filter to converge, and it is only applicable to low dimensional problems.

One of the most popular stochastic methods is Ensemble Kalman Filter (EnKF). The EnKF is
marginally similar to the particle filter. However, it is assumed that the predicted prior pdf is well
approximated by only the first and second order moments of the pdf. Therefore each realization
is updated according to the standard Kalman filter update equation, but the measurement errors.
The EnKF is very simple to implement and can be developed for a black box model, without
requiring any changes to the simulator. It has shown promising results for data assimilation in
nonlinear dynamical systems in various applications [Evensen and van Leeuwen 1996, Reichle et
al., 2002; Naevdal et al., 2005; Chen and Zhang, 2006; Aanonsen et al., 2009; Jafarpour and
McLaughlin, 2009; Jafarpour and Tarrahi, 2011]. While the filter has well-known limitations,
including sampling errors and sub-optimality for nonlinear non-Gaussian problems, it is

25
recognized as a practical estimation method in applications with weakly to moderately nonlinear
models and nearly multi-Gaussian problems.

1.3 Description of Rock Facies Distributions Using Sparse Data
Variability in rock facies distribution and within facies heterogeneity is another source of
uncertainty in predicted performance and monitoring data. A complicating factor in accurate
description and prediction of fluid flow and transport processes is the complex heterogeneity in
rock physical properties that directly impact the fluid flow and transport behavior [Carle et al.,
1998; Zinn and Harvey, 2003; de Marsily et al., 2005]. Heterogeneity exhibits itself across
multiple scales, from pore to basin scale. In the face of data scarcity, identification of the
heterogeneity in hydraulic properties of geologic formations at all scales presents a formidable
challenge to the subsurface flow modeling community. Geostatistical modeling techniques have
traditionally been adopted to represent the heterogeneity in geologic formations [Deutsch and
Journel, 1998]. Classical geostatistical modeling approaches use simplifying assumptions, which
are often driven by the ease of mathematical treatment (e.g., use of variogram or covariance
functions), to describe the natural variability in geologic formations. In some cases, the adopted
statistical/mathematical models can lead to over-simplification of the existing geologic features
and their connectivity. Together with data scarcity, imperfect modeling assumptions introduce
significant uncertainty into fluid flow and transport predictions [Carle et al., 1998; Zinn and
Harvey, 2003; de Marsily et al., 2005; Hill and Tiedeman, 2007; Oliver et al., 2008].
A common practice to improve the predictive power of subsurface flow models is their
calibration against observed monitoring and performance data [e.g., Carrera and Neuman, 1986a-
b; Dai and Samper, 2004; Hill and Tiedeman, 2007; Yeh and Khaleel, 2008; Oliver et al., 2008;

26
Zhou et al., 2014]. In general, the resolution of dynamic flow and transport data seldom lends
itself to characterization of heterogeneity at small scales, often necessitating a coarse-scale
description of rock properties for model calibration. Examples of coarse-scale model parameter
representations include upscaling/zonation, low-rank parameterization and regularization
techniques (see e.g., [Jacquard and Jain, 1965; Jahns, 1966; Tikhonov and Arsenin, 1977;
Gavalas et al., 1976; Shah et al., 1978; Lawson and Hanson, 1995; Weiss and Smith, 1998;
Chavent and Bissell, 1998; Grimstad et al.., 2004; Linden et al., 2005; Tonkin and Doherty,
2005; Doherty and Skahill, 2006; Jafarpour and McLaughlin, 2009]). One of the key issues in
formulating model calibration problems is reconciling the discrepancies between model and data
resolutions so that the parameters of interest properly represent the salient and relevant
information (to fluid flow and transport processes) while enabling the estimation of parameters
from available (often low-resolution) pressure and flow rate measurements. While choosing
optimal parameter resolution is important and challenging from an inverse modeling perspective,
the main objective of model calibration is to improve the prediction of flow and transport
behavior.
An important property in predicting the fluid flow and transport behavior in a geologic formation
is the connectivity of the extreme (low and high) values in the hydraulic conductivity field. The
connectivity in high-conductivity rocks can form a preferential flow and transport pathway (with
minimum resistance to fluid flow) while the connected regions with low conductivity can act as a
flow barrier that impede fluid flow. One way to characterize connectivity in subsurface rocks is
through characterization of rock facies types and their spatial distribution [Deutsch and Journel
1998; Carle et al., 1998, Western et al., 2001; Zinn and Harvey, 2003; Ritzi et al., 2004, 2006; de
Marsily et al., 2005]. The shape and connectivity of rock facies directly affects the flow and

27
transport behavior. Facies connectivity is particularly important when the hydraulic properties of
different facies types have sharp contrasts; such is the case in fluvial systems where high
permeability sand channels are embedded in low-permeability mud or shale rocks [Guardiano
and Srivastava, 1993; Deutsch and Wang, 1996; Holden et al., 1998; Strebelle, 2002; Zinn and
Harvey, 2003; Caers and Zhang, 2004] .
In practice, significant uncertainty exists in describing rock facies distributions, especially away
from observations points, which typically consist of scattered injection, extraction, and
monitoring wells. A common way to reduce the uncertainty in facies distribution is through the
integration of dynamic response data that often encode important information about interwell
flow connectivity [Hill and Tiedeman, 2007; Oliver et al., 2008; Oliver and Chen, 2011; Zhou et
al., 2014]. Within-facies heterogeneity presents another form of variability in describing
heterogeneous rock properties. However, in many cases, within facies variability is less
impactful than cross-facies heterogeneity [Ritzi et al., 2004, 2006]. While presence of small-
scale heterogeneity is important for capturing the transport behavior in more detail, the resolution
and content of scattered dynamic response data generally do not afford the ability to resolve such
level of detailed variability in flow models.
To constrain a large number of spatially distributed parameters using sparse data, additional
assumptions are generally made about the parameter field. Examples include Tikhonov
regularizations [Tikhonov, 1963; Tikhonov and Arsenin, 1977] and geostatistical methods
[Kitanidis and Vomvoris, 1983; Dagan, 1985; Carrera and Neuman, 1986; Kitanidis, 1995; Yeh
et al., 1995; Zhang and Yeh, 1997]. Geostatistical methods assume that parameter fields can be
modeled as realizations of a spatially correlated random field and incorporate information related
to the expected spatial correlation (e.g., variogram function) to penalize departure from those

28
models. For more structured random fields, in which the main component of parameter
variability is between geologic facies, variogram-based methods are not suitable. This type of
parameter distribution is not uncommon and data from wells and other measurements have
shown sharp contrasts in parameter values at geologic facies boundaries with relatively little
intrafacies variability [e.g., Fienen et al., 2004]. Also, in hydrologic characterization at the
watershed-to-basin scale [Robinson et al., 2008; Abu-el-Sha'r and Rihani, 2007], identification of
the boundaries between distinct geologic facies is of primary importance.
One approach for imaging of structured parameter fields, such as facies distribution, is to reduce
the problem to the estimation of a few meta-parameters that control the shape and location of the
structures. Included in this family are zonation [Jacquard and Jain, 1965; Jahns, 1966],
multiscale methods, which successively subdivide an initially homogeneous parameter field into
smaller homogeneous regions in order to improve data fit, using predefined refinement methods
(e.g., halving along one dimension) [Yeh and Yoon, 1981; Liu, 1993; Ben Ameur et al., 2002;
Grimstad et al., 2003], control point methods that reproduce the unknown parameter field using a
set of points at which parameter values are estimated and define an interpolation method to be
used between these points [Sun and Yeh, 1985; Tsai et al., 2003; Tung and Tan, 2005; Tsai and
Yeh, 2004], postprocessing routines to detect zones within images generated by other inverse
methods [Hyndman et al., 1994; Eppstein and Dougherty, 1996; Jarvis and Knight, 2002;
Tronicke et al., 2004; Fienen et al., 2008].
Calibration of rock facies distributions against dynamic flow and transport response
measurements presents several challenges. Most notably, facies maps are described as discrete
variables. While this limits the possible outcomes (facies type) that each grid block can assume,
it introduces a hurdle in using well-established model calibration methods that are designed for

29
continuous parameters. Furthermore, the discrete nature of these facies implies that their
distribution is far from multi-Gaussian and that their treatment within the classical stochastic
inverse modeling formulations is not convenient. Because of these inherent properties of facies
models, direct application of continuous model calibration methods to identify facies distribution
is not feasible. This difficulty has led to several efforts to manipulate model calibration
formulations to enable facies model identification, including direct conditional facies simulation,
discrete regularization, and continuous parameterization.

1.3.1 Direct Conditional Facies Simulation
A more direct approach for calibrating facies models is through conditioning on flow data using
geostatistical simulation techniques, instead of post-simulation updating. Direct conditional
simulation has the inherent advantage of generating samples that are, by construction, consistent
with the prior model used to generate them. Recent examples of these techniques are the
probability perturbation method (PPM) [Caers, 2003; Hoffman et al., 2005; Hu, 2008], the
probability conditioning method (PCM) [Jafarpour and Khodabakhshi, 2011] and a stochastic
inverse method proposed by [Harp et al., 2008] for inversion of a transition probability model in
an analytical framework, where statistical structural parameters such as mean length are used to
update transition probability model. The PPM uses a parameterization of the random seed of the
stochastic (SNESIM) facies simulation to generate a new realization that improves the data
match obtained by the current solution. This process is repeated by changing the random seed
and generating a new iterate until no further improvement in data match is observed, or the
maximum number of iterations/simulations is reached. The PCM, on the other hand, uses the
flow data to estimate a facies probability map that provides important information about facies

30
distribution. The estimated facies probability map is then used as input into the SNESIM
algorithm to simulate updated facies models from the training image. Therefore, PCA indirectly
incorporates the flow data (by using them to estimate the facies probability maps) to generate
conditional facies realizations. In a recent paper, Ma and Jafarpour [2017] use the concept of
pilot points method to introduce a new approach for conditioning discrete facies simulation on
nonlinear flow data. Similar to the PCM, the pilot-point-based method uses the flow data to
identify facies type at pilot points, which are then included as hard data to generate conditional
facies models (using the SNESIM algorithm). Other probabilistic parametrization methods for
geologic facies include alternative forms of the pluri-Gaussian simulation, in which continuous
Gaussian random fields are projected onto discrete facies maps using a thresholding approach
[Liu and Oliver, 2005; Sebacher et al., 2013].

1.3.2 Discrete Regularization
Discrete regularization methods can also be used to promote discrete solutions to facies model
calibration problems. Regularization methods are often implemented as a (generally soft) penalty
to encourage solutions with the desired property [Tikhonov and Arsenin, 1977; Carrera et al.,
2005; Hunt et al., 2006; Oliver et al., 2008]. A popular example is the use of total variation (TV)
regularization penalty [Rudin et al., 1992; Lee and Kitanidis, 2013] that allows local
discontinuity features (e.g., facies boundaries) to be detected by penalizing the l1-norm of the
spatial gradient of the solution (unlike the harsher l2-norm penalty that is used in Tikhonov
regularization [Tikhonov and Arsenin, 1977] to obtain smooth solutions). Sparse reconstruction
techniques, especially those based on l1-norm minimization (formalized under the compressed
sensing paradigm) [Donoho, 2006, Elad, 2010] have also recently been applied to estimate

31
geologic facies [Jafarpour et al., 2010, Li and Jafarpour, 2010; Khaninezhad, et al., 2012a-b].
Discrete regularization techniques, which are widely used in discrete tomography, provide
another set of methods for estimating discrete geologic facies. These methods introduce discrete
penalty functions into least-squares model calibration formulations to promote discrete solutions.

1.3.3 Continuous Parameterization
A general class of continuous parameterization methods for discrete geologic facies is low-rank
linear expansion techniques. Examples of these transform-domain parameterization techniques
include the principal component analysis (PCA) [Jolliffe, 2005; He et al., 2013] and its kernel
version known as K-PCA [Sarma et al., 2008], the discrete cosine transform (DCT) [Rao and
Yip, 1990; Jafarpour and McLaughlin, 2009], grid connectivity transform (GCT) [ Bhark, et al.,
2011], and the discrete wavelet transform (DWT) [Lu and Horne, 2000; Mallat, 2008; Li and
Jafarpour, 2010; Jafarpour, 2011; Kovacevic et al., 2013]. Low-rank parameterization with these
methods leads to close approximation of the original model while substantially reducing the
number of parameters. The basis images and their corresponding expansion coefficients in these
transform-domain methods are continuous. As a result, continuous model calibration methods
can be conveniently used with these parameterization techniques. However, the final calibrated
model is generally continuous and violates the discrete nature of the facies. A common property
of these parameterization methods is that they tend to provide a smooth large-scale
approximation, which helps with capturing the global facies connectivity. Furthermore, the
resulting transform-domain parameters tend to be near-Gaussian, a property that is desirable in
employing classical estimation methods that use multi-Gaussianity assumption, e.g., Kalman
filtering [Kalman, 1960] and its extensions.

32
Other transform-domain representations are also used to impart multi-Gaussian properties on the
parameters (without reducing the number of parameters) in order to improve the performance of
model calibration methods that use multi-Gaussian assumptions. A recent example of such
methods is the use of CDF transformation (e.g., the normal score transform) before calibration of
non-Gaussian groundwater model parameters with the ensemble Kalman filter (EnKF) [Zhou et
al., 2011]. Gaussian Mixture Models [Dovera, 2010] have also been introduced to approximate
multimodal priors and use a modified implementation of the EnKF to update Gaussian mixture
(GM) distributions. In this method, non-Gaussian distributions are described as weighted sums of
Gaussian PDFs. Another set of parameterization methods for discrete models that also involve
transformation to continuous domains are those that provide exact discrete representations after
back transformation. A popular approach in this category is the level-set method [Osher and
Sethian, 1988; Santosa, 1996; Cardiff and Kitanidis, 2009; Chang et al., 2010, Moreno et al.,
2011, Lorentzen et al., 2012], in which facies boundaries are described using continuous
parameter fields. In this case, during model calibration, the continuous parameters that control
the location and shape of facies boundaries (through the zero values of the underlying signed
distance function) are updated using the flow data.

33
CHAPTER 2
Dynamic Fracture Characterization
In this chapter, WEexamine the applicability of the proposed EnKF-based inversion framework
for characterizing hydraulic and natural fracturing parameters and matrix properties by
integrating tracer and production data in series of numerical experiments (with planar and DFN
models). The workflow begins with generating an ensemble of model parameters. Each
realization of these parameters has different fracture and matrix properties such as fracture half-
length, fracture density, matrix permeability and matrix porosity and is used to predict the
dynamic response of the reservoir based on the specified parameter realization. These responses
are then integrated with the observed response data, using the EnKF update equation, to update
the initial parameters. The main steps of the framework are schematically shown in Figure 1.

2.1. Fracture Models and Parameters
In this work, we are interested in characterizing both planar fractures and complex fracture
networks from flow-related data. Figure 2 shows a simple schematic of planar (2a) and DFN
(2b) models. Each of these models and their underlying assumptions and important parameters
are described in this section.

2.1.1 Planar model
Planar fractures are typically modeled with rectangular or oval shapes and can be parameterized
using their simple geometric attributes, e.g., the length and height of the fracture. The
conductivity of the fracture is assumed to be constant in the entire fracture plane. In our

34
numerical simulations, the fracture is modeled as very high permeability cells. Since the
minimum dimension in our model (3 ft) is larger than realistic fracture apertures (10-2-10-1ft), in
our model the permeability is adjusted to account for the difference between the realistic and
modeled fracture apertures. To compensate for this discrepancy, the fracture permeability in the
model is reduced according to the following equation
𝐹𝐶𝐷 = 𝐾 1
𝑊 1
= 𝐾 2
𝑊 2
(3)

where K and W represent fracture permeability and aperture, respectively. Based on this relation,
for fracture aperture and permeability of 0.10 ft and 120 D (Crain, 2010), respectively, the
corresponding adjusted fracture permeability (for a 3ft aperture) is 4 D.

Another important effect to model is the loss of fracture conductivity due to fracture closing,
which can have an important impact on production performance. To model this effect, we
assume that fracture conductivity changes with pressure. In the literature, the permeability of
unconventional reservoirs is reported to change during the production (Clarkson et al. 2012).
Different mechanisms, including fracture closure under high effective stress and proppant
embedment and crushing due to stress changes, could change the fracture conductivity. The
impact of stress-dependent fracture conductivity variations has been investigated recently
(Clarkson et al. 2012). In Haynesville Shale, both matrix permeability and fracture conductivity
have been reported to reduce during production (Diaz de Souza et al. 2012). Incorporating this
effect in modeling the dynamic behavior of fractured reservoirs is necessary for reliable
identification of fracture properties and prediction of production and optimal re-fracturing time.

35
A prevailing approach to account for permeability changes in conventional simulators is
introduced by Yilmaz and Nur as follows (Yilmaz et al. 1991):
𝐾 = 𝐾 𝑖 𝑒 B
𝑝 (𝑃 −𝑃 𝑖 )
(4)

where 𝐵 𝑝 is the permeability modulus, K is the absolute permeability, and P is pressure. In
several studies, this equation has been used to model permeability changes in unconventional
reservoirs (Franquet et al. 2004, Thompson et al. 2010, Okouma et al. 2011). Figure 3a (first
row) shows permeability changes as a function of time and pressure. Because permeability is
typically measured as a function of stress (Clarkson et al. 2012), the relation in Eq. (4) can be
expressed regarding stress variations as follows (Bachman et al. 2011):
K = K
i
e
(-B
ơ
Δơ)
(5)

Here, B
ơ
is the permeability modulus, and ơ is the average effective poroelastic stress. The
average effective stress can be converted to pressure changes using a uniaxial-strain state
(Batchman et al. 2011, Settari et al. 2005)
Δơ = -α
c
1+ν
3(1-ν)
ΔP (6)

where, 𝛼 𝑐 is the Biot’s constant and ν is Poisson’s ratio. By plugging Eq. (6) in Eq. (5) and
comparing it with Eq. (4), the permeability module, for the case of uniaxial stress state, can be
calculated as
B
p
= B
ơ
α
c
( 1-
2
3
η
α
c
) , η =
α
c
(1-2ν)
(1-ν)
(7)

36
The permeability modulus (B
𝑝 ) and initial permeability (Ki) are the two additional important
uncertain parameters in Eq. (4) which should be updated. We studied the effect of fracture
permeability changes on pressure and cumulative production (Figure 3b, second row) and, as
shown, found that fracture permeability reduction has a significant effect on cumulative
production. Therefore, being able to predict changes in fracture permeability during the life of
the fracture is imperative to improving the hydraulic fracturing design and enhancing the
effectiveness of re-fracturing.

Table 1 summarizes a list of uncertain fracture and reservoir properties along with their typical
range. Tables 2 and 3 provide a summary of the fracture and reservoir properties for the
reference model. The fracture half-length, conductivity, decline rate, as well as the matrix
permeability and matrix porosity are assumed uncertain and described with probability
distributions, which are used to generate an initial ensemble for these parameters. The
corresponding ensemble of fracture parameters and reservoir properties are provided as input
into the 3D model to perform Monte-Carlo simulation and predict the ensemble of observed data
for assimilation.

2.1.2 Discrete fracture network
For the discrete fracture network model, we used a stochastic object-based method to generate
3D fracture networks in the entire reservoir. The goal of model calibration for stochastic
fractures is to find the open fracture density (that is, fractures that contribute to flow and
production) along with induced and natural fracture hydraulic conductivity values, as well as the
matrix permeability and porosity. The dynamic decline in both natural and hydraulic fracture

37
conductivities with decreasing pressure is also considered, as described in the previous section.
The parameters of our stochastic simulations are fracture density (location) and size (length)
(Cacas et al. 1990). We assumed a Poisson distribution for fracture density and a log-normal
distribution for fracture length. First, the coordinates of fracture centers are generated using the
Poisson's distribution and a random number generator. The firing rate (λ) in Poisson distribution
is related to the fracture density. For instance, if reported fracture density is n per cubic feet and
the size of the reservoir model is Nx×Ny×Nz cubic feet, the firing rate of Poisson distribution
equals to n×Nx×Ny×Nz. In other words, Poisson's firing rate represents the expected fracture
quantity in the reservoir model. In the next step, fracture length is randomly generated using a
log-normal distribution with the mean of 10 ft and standard deviation of 4 ft. After generating the
spatial distribution of natural fracture elements, hydraulic fracture half-length, hydraulic and
natural fracture conductivities, hydraulic and natural fracture decline rates, matrix permeability
and porosity are randomly generated using a uniform distribution with the specific range
summarized in Table 1. Table 4 is a summary of true fracture parameters for discrete fracture
network model. Figure 2(b) illustrates a fracture network generated using this stochastic
simulation process and dynamic change of hydraulic and natural fracture conductivities over
production time is shown in Figure 4.

For simulation and data assimilation purpose, we consider two stages consisting of tracer test (a
transport model), and prediction production (flow model). The tracer test is modeled using the
governing transport equations. During the hydraulic fracturing treatment, tracer is injected into
the formation of the fracturing fluid. After treatment, during flow back, the tracer concentration
levels are recorded for each stage of hydraulic fracturing. The governing equation for an

38
environmental tracer in a single phase fluid flow system is given in Eq. (8) (see Appendix A for
notations):
d
dt
(
VSC
B
) +
d
dt
(V ρ
r
C
a
1-ϕ
ϕ
) = ∑ [
Tk
r
Bµ
( δP- ρgD
z
)C + DFD
c
S δC] + QC-V
S
B
λC (8)

where S, P, C, and 𝐶 𝑎 represent, the host phase saturation, pressure, flowing tracer concentration
and adsorbed tracer concentration, respectively. The notations V, Q, and D
z
stand for the block
pore volume, host phase production rate, and center depth, respectively. Terms 𝜌 𝑟 , 𝜌 , 𝜙 , T, 𝑘 𝑟 ,
B, µ, 𝐷 𝑐 , DF, 𝜆 and 𝑔 respective;y denote mass density of the rock and host phase, porosity,
transmissibility, relative permeability, formation volume factor, viscosity, tracer diffusion
coefficient, diffusivity, tracer decay constant and gravity acceleration, respectively.
Tracer decay is modeled as a simple half-life process, which implies that the tracer starts to
decay as it is injected into the reservoir according to its specified half-life. The diffusive flow of
tracer from cell i to a connected cell j is expressed in Eq. (9).
F
t
= DF. D
c
. S. δC = DF. D
c
. S(C
ci
-C
cj
) (9)

The diffusivity (DF) has the form DF =
Aϕ
d
where A, d, and ϕ denote, respectively, the cell
interface area, the distance between cell centers, and porosity. We solve the diffusivity equation
using a finite difference-based commercial fluid flow simulator (Eclipse, 2011). Assuming that
convection is the main transport mechanism for the tracer (i.e., neglecting diffusion, decay, and
adsorption), Eq. (8) can be simplified to:

39
𝑑 𝑑𝑡
(
𝑉𝑆𝐶
𝐵 ) = ∑ [
𝑇 𝑘 𝑟 𝐵 µ
(𝛿𝑃 − 𝜌𝑔 𝐷 𝑧 )𝐶 ] + 𝑄𝐶 (10)

Production prediction is also modeled using the continuum approach with adjusted fracture
properties as discussed in Appendix B. Adjustments are applied to account for the
approximations involved in representing fractures as porous media with unrealistically large cell
sizes for fracture aperture. The simulations in this study are performed using Eclipse 2011 flow
simulation software.

2.2Data Assimilation with EnKF
For dynamic characterization of fracture and matrix properties from tracer test and production
data, we use the ensemble Kalman filter (EnKF), which is a nonlinear approximation of the
classical Kalman filter (KF) (Kalman, 1960). The KF is a sequential state estimation method that
characterizes the first and second statistical moments (mean and covariance) of the states
posterior distribution. For a linear state-space system with jointly Gaussian states and
observations, KF presents a full characterization of the multi-Gaussian posterior distribution
(Kalman, 1960; Gelb, 1974). For nonlinear models (or non-Gaussian priors), EnKF provides an
approximation to the KF, which has found widespread application in practice (Evensen and van
Leeuwen 1996; Evensen, 2007; Evensen, 2009). The main appeal of the EnKF is related to its
simple and gradient-free implementation, which does not require changes to the simulator. The
filter has shown promise for data assimilation in nonlinear dynamical systems in various
applications (Evensen and van Leeuwen 1996, Reichle et al., 2002; Chen and Zhang, 2006) and
more recently for history matching of oil reservoirs (Naevdal et al., 2005; Aanonsen et al., 2009;

40
Jafarpour and McLaughlin, 2009; Jafarpour and Tarrahi, 2011). While the filter also has its well-
known limitations, including sampling errors and sub-optimality for nonlinear non-Gaussian
problems, it is recognized as a practical approach in applications with weakly to moderately
nonlinear models and nearly multi-Gaussian statistical characteristics.
The EnKF consists of a sequence of forecast and update steps. The forecast step uses Monte-
Carlo simulation to approximate, through an ensemble of model inputs, the predicted
observations of the system using nonlinear forward simulations. This step can generally be
represented as
x
t|t-1
j
= f(x
t-1|t-1
j
, m
j
) , y
t|t-1
j
= h (x
t|t-1
j
) + ε
j
(11)

where the first function represents a nonlinear state transition operator while the second function
is a generally nonlinear measurement operator that related current system states to observations.
In Eq. (11), x, y, 𝘀 and h represent state vector, perturbed observation vector, uncorrelated
Gaussian observation error, and observation operator, respectively (see Appendix A), the
subscript (t|s) in x
t|s
j
refers to state variables at time t that are conditioned on all measurements
up to time s and the superscript j denotes the jth sample in the ensemble. The ensemble of
predicted states and observations is used to compute the necessary covariance and cross-
covariance information for the analysis step.

In the analysis step, the predicted ensemble of states and observations (and model parameters) is
updated using the observed dynamic responses from the actual systems (field data). In history

41
matching applications, it is common to estimate the parameters of the reservoir and use them to
predict the state variables (e.g., pressure and tracer concentration distributions). This will ensure
that the reservoir states and parameters remain consistent and honor the governing physical laws
of flow and transport in porous media. In this case, the augmented states of the system that are
updated at the analysis step are defined as z
j
= [
m
j
h
j
]
T
z
j
= [
m
j
h
j
]
T
, where m and g
represent model parameters and predicted observations. The EnKF analysis step can then be
compactly expressed as
z
u
j
= z
f
j
+ G(y
j
-Hz
f
j
) , 𝐾 = 𝐶 𝑧 𝑒 𝐻 𝑇 (𝐻 𝐶 𝑧 𝑒 𝐻 𝑇 + 𝐶 𝑦 )
−1
(12)

where, G, H, 𝐶 𝑧 𝑒 , and 𝐶 𝑦 are Klaman Gain, measurement, state sample covariance, and
observation covariance matrices, respectively. The subscript (u) and (f) refer to the forecast and
updated augmented states.

In our data assimilation experiments, the data consist of time-dependent tracer concentration
observations and oil flowrate for each stage (the parameters in each example are defined in the
corresponding description of the experiments). The forecast and updates steps in Eqs. (11) and
(12) are sequentially implemented until all observations are assimilated.
In practice, computational resources limit the ensemble size that can be used and result in
sampling errors, which can lead to erroneous updates due to unphysical (spurious) correlations.
Additionally, an important effect of the spurious correlations is the underestimation of the
ensemble variance and potential filter divergence. Two approaches are proposed to alleviate the

42
impact of the spurious updates: ensemble inflation and localization (or local analysis) (Evenson
2007). In our EnKF implementation, we used a damping factor to reduce the effect of the
Kalman gain in Eq. (12). The damping factor is a parameter in the range of [0,1] that controls the
weight of observations, relative to the prior model, in the update equation. In effect, the damping
factor can be viewed as inflation that is applied to the observation error variance to avoid sudden
and large updates to the prior models. The use of a damping factor was introduced by Anderson
(2007) and was later included in the ensemble smoother with multiple data assimilation (ES-
MDA) update scheme that was used by Emerick and Reynolds (2013). To compute the damping
factor in our implementation, we follow the proposed approach in (Anderson 2007). In our
implementation, the forecast step of the filter remains unchanged while the update step is
modified to:
𝑥 𝑢 𝑗 = 𝑥 𝑓 𝑗 + 𝛼𝐺 (𝑦 𝑗 − 𝐻 𝑥 𝑓 𝑗 ) (13)

in which 𝛼 is the damping factor. The observations are perturbed using uncorrelated Gaussian
noise with covariance matrix C
d
, where the standard deviation of the observation error is equal to
20% of the observed values (𝑦 𝑜𝑏𝑠 ). That is, the realization j of the perturbed observation vector
can be expressed as
𝑦 𝑗 = 𝑦 𝑜𝑏𝑠 + 𝘀 𝑗 , 𝑗 = 1,2, … , 𝑁 𝑒 , 𝘀 𝑗 ~𝑁 (0, 𝑜 2
) (14)

43
In this paper, we assume an uncorrelated Gaussian observation error with zero mean and the
standard deviation defined above. Next, we present the numerical experiments and the results
obtained using the above methods, followed by a discussion of our observations.

2.3 Results and Discussion
We consider a 3D reservoir model of size 600 ft×800 ft×90 ft. As illustrated in Figure 2, four
stages of hydraulic fracturing with different fracture properties are assumed. In the first set of
experiments, only planar fracture models are considered (Figure 2a) while in the second set of
simulations a DFN model is adopted (Figure 2b). For the planar model, each fracture is modeled
as a rectangular plane with fracture length, conductivity and conductivity decline rate as
parameters. The fractures are assumed to cover the total depth of the shale rock formation and
are perpendicular to the horizontal well. For complex DFN model, the hydraulic and natural
fractures are modeled with rectangular shapes. The natural fractures are assumed to be
perpendicular to the main hydraulic fractures.

Experiment 1: Planar Model and Tracer Data
In Experiment 1, we use a planar fracture model as a commonly adopted simplified description
of hydraulically induced fractures. The model assumes large stress anisotropy with the minimum
and maximum stresses aligned in the x and y directions, respectively. The experimental setup is
designed to simulate a horizontal injection well with completion in the y-direction (perpendicular
to the minimum stress direction). A single phase fluid (water) with tracer is injected into the
formation with no-flow boundary conditions. The reservoir properties for this example are

44
reported in Table 3. In this example, we assume that the hydraulic fracture design is the same for
all four stages and that water is injected with tracer to induce the fracture.

The unknown parameters are hydraulic conductivity and half-length of fractures, and
permeability and porosity of the rock matrix. This experiment is designed to evaluate the use of
tracer data in estimating fracture parameters simultaneously. Tracer data are recorded for a
minimum of 24 hours and a realistic maximum of 60 days. Recently, tracer data has been used to
assess the performance of the fracturing process; however, as tracers are attached to the proppant
or fracturing fluid, they can provide important information about fracture propagation as well as
fracture geometry and conductivity. It is important to note that tracer test and production
responses follow different physics (transport versus flow, respectively) and take place at
different times and over different time periods. Therefore, in our experiments, we first assimilate
the tracer test data using a transport model to reduce the initial uncertainty in each parameter. In
this case, we simulate the tracer data for 60 days, followed by production from the fracture for
360 days. Furthermore, the decline rate of fracture conductivity is only considered during the
production phase as tracer tests typically occur over a short period when production (pressure
drop) has not started; therefore, during tracer tests, fracture conductivity is not expected to
decline significantly.

Before data assimilation, a sensitivity analysis is performed to evaluate the impact of fracture
properties on observed responses of the model. The conductivity and half-length of planar
hydraulic fractures are estimated from tracer concentration data. Figure 5 depicts the tracer
concentration travel path. In the case of the planar model, the sensitivity of the tracer test and

45
production data to fracture half-length, conductivity, matrix permeability and matrix porosity is
evaluated and shown in Figure 6. When the fractures at different stages have the same
conductivity, but different half-lengths, the sensitivity analysis results suggest that as fracture
half-length increases so does the arrival time of a fixed tracer concentration. Intuitively, tracers
injected into longer fractures take longer time (as they travel farther) before they are recovered
during flowback. The effect of fracture conductivity on tracer concentration is shown in Figure
6b (middle), where the fractures at different stages have the same half-length, but different
conductivity values. As fracture conductivity increases, the arrival time of a fixed tracer
concentration decreases, suggesting that fractures with higher conductivity can lead to faster
tracer recovery during the backflow. The last row in Figure 6 shows that the production is more
sensitive to fracture half-length than it is to fracture conductivity.

Figure 7 summarizes the sensitivity of tracer concentration and production data to fracture and
matrix properties. Here, 10-40% increase is applied to each parameter and the corresponding
cumulative oil production and mean tracer concentration over time have been plotted. The results
in Figure 7a indicate that tracer data are mostly sensitive to fracture conductivity. Matrix
permeability and porosity do not show a strong correlation with tracer response, mainly because
the short duration of the tracer test may not allow much interaction between the DFN fractures
and hydraulically induced fractures.

For data assimilation, the ensemble of model realizations is generated by assigning permeability
and porosity values at fracture and matrix grid blocks. The initial realizations have different
(stochastically described) half-length, hydraulic conductivity, matrix permeability and porosity.

46
The ensemble size is 100, and a uniform PDF is used to generate the initial realizations for these
parameters. The results of simultaneously estimating these parameters, from tracer data only, are
shown using the black (empty) histograms in Figure 8a (the results are shown for Stage 1,
similar results are obtained for other three stages). The black box shows the theoretical uniform
distributions that are used to generate the initial parameter ensembles. The histograms of the
updated parameter ensembles (after integration of tracer data only) are shown in black. The
updated histogram of fracture conductivity shows smaller variance around the reference values
(red line). However, the updated histogram for fracture half-length shows a large variance, which
is consistent with the results of the sensitivity analysis. The updated matrix permeability and
porosity histograms also show large variability due to the lack of sensitivity of tracer data to
these parameters. The predicted tracer concentration profiles, for the initial (cyan) and final
(grey) ensembles, are shown in Figure 8b (left). The results show that the reduction in the
variance of fracture hydraulic conductivity has had a noticeable impact on reducing the ensemble
spread. However, the final ensemble of predicted tracer concentration profile still shows some
uncertainty. The results from our first experiment show that while tracer data can be used to
reliably estimate fracture hydraulic conductivity, complementary data are needed to constrain
fracture half-length and matrix permeability and porosity.

Experiment 2: Planar Model, Dynamic Fracture Conductivity, and Production Data
In this experiment, we use the results from the assimilation of tracer test data in Experiment 1 as
input (initial ensemble) for integration of production data. Also, we include the dynamic fracture
conductivity decline rate as a new parameter. The results from sensitivity analysis for production
data are shown in Figure 7b and suggest strong sensitivity to matrix permeability and some

47
sensitivity to fracture half-length, due to its impact on the stimulated reservoir volume, and
matrix porosity and fracture decline rate. For the same fracture length, fracture conductivity has
little impact on cumulative production; the only exception being the early time production, which
is dominated by production from the fractured zones and their immediate vicinity. For
permeability values higher than 20mD, up to 9000mD, known as infinite conductivity, this effect
is not noticeable. Figure 7c shows the sensitivity of pressure response to fracture and matrix
properties, where only matrix permeability shows a significant effect on the pressure response.
In Experiments 2, it is assumed that as reservoir undergoes production pressure in the reservoir
depletes and the effective stress increases, leading to the gradual closing of the fractures. In this
example, we dynamically update the permeability of the planar fracture by uniformly updating
the permeability values within fracture cells using updates based on production data.

The final results from the sequential integration of tracer and production data are also shown in
Figure 8. In Figure 8a, the final histograms of rock matrix and fracture properties are shown in
blue. The updated fracture half-length shows less uncertainty after integration of production data
while the histogram of fracture half-length does not change after integrating the production data.
Integration of production data has also resulted in reduced uncertainty about matrix permeability
and porosity, and the fracture conductivity decline rate. In particular, the spread in matrix
permeability is significantly reduced, suggesting that stronger update is applied to it due to its
correlation with production data (as shown in Figure 7b). The corresponding predicted
cumulative production profiles are shown in Figure 8b (second and third column). The results
show a significant reduction in production uncertainty after data assimilation.

48
The results from Experiments 1 and 2 demonstrate that tracer test and production data provide
complementary information about rock and fracture properties for planar models. Tracer test data
are primarily sensitive to fracture conductivity and less sensitive to fracture half-length and
matrix properties (permeability and porosity). On the other hand, production data show strong
sensitivity to parameters that directly affect the produce volumes, mainly matrix permeability
and porosity and fracture half-length (which controls the SRV). Production data does not show
much sensitivity to fracture conductivity, except at very early stages where production is mainly
from the fracture cells and their immediate neighborhood. The data assimilation results are
consistent with the outcome of our sensitivity analysis and show that the EnKF has promising
performance in updating parameters that show sensitivity to each respective data type. In the
next experiment, we explore the role of tracer test and production data in estimating the
parameters of a DFN model and examine the performance of EnKF in characterizing hydraulic
fracturing of formations that contain natural discrete fractures.

Experiment 3/4: Stochastic Characterization of DFN Using Tracer and Production Data
In our next experiments, we consider hydraulic fracturing of a formation that contains complex
natural fractures. We apply our sensitivity analysis and data assimilation approach to
characterize the stochastic parameters of the discrete fracture network model along with the
planar fracture and matrix properties. We first evaluate the sensitivity of tracer, pressure and
production data to fracture and matrix properties. The list of parameters that are included in our
sensitivity analysis can be found in Table 1.

49
In these experiments, the reservoir domain is divided into four regions, each corresponding to a
fracturing stage. Tracer, production and pressure data are used to estimate fracture density
(number of fracture per unit volume), along with other fracture and matrix properties (Table 1).
Table 4 lists the reference fracture density other fracture properties for each region. It is assumed
that the natural fractures are all propped open and contribute to production. In our
characterization of the DFN, we have not attempted to estimate individual fracture attributes
(e.g., location, length, conductivity, etc.) as the monitoring and production data do not afford the
detailed resolution needed to resolve natural fractures individually; instead, we have focused on
estimating the parameters of the stochastic model which is used to generate the DFN.

Figure 9 summarizes the results of our sensitivity analysis, which shows that tracer response
data are mostly sensitive to hydraulic fracture conductivity (Figure 9a), similar to the planar
model, and fracture density. The lack of sensitivity to other parameters is mainly related to
limited interactions between the tracer and the matrix and the fracture network. The sensitivity of
tracer data to fracture density can be explained by noting that a higher fracture density increases
the probability of having fractures closer to the tracer injection point, and more interaction
between fractures and the injected tracers. On the other hand, production data show sensitivity to
fracture density, half-length, matrix permeability, hydraulic fracture decline rate (Figure 9b). It
also shows slight sensitivity to hydraulic fracture conductivity and matrix porosity. Figure 10
illustrates the lack of sensitivity of production data to hydraulic and natural fracture
conductivities. To show this, we have used an identical ensemble of all model parameters except
or the parameter under study. For this parameter, there (low, medium and high) values are
assigned. For the entire ensemble, the histogram is plotted for the difference between the

50
forecasts forecast of the high and medium values and low and medium values. For both hydraulic
and fracture conductivities, Figures 10a and 10b, respectively, these differences are very small,
suggesting that over all values of other parameters (randomly assigned in generating the
ensemble), the values of hydraulic and fracture conductivity do not affect the production
forecast. For comparison, we have shown the same results for fracture density (Figure 10c). The
impact of changes in fracture density on production forecast is far more significant. The strong
sensitivity to fracture density is due to the direct relationship it has with the stimulated reservoir
volume and cumulative production from the reservoir. As a result, in our data assimilation
experiment, we expect that the production data can strongly constrain the fracture density for
different fracturing stages.

For data assimilation, we followed the same sequential approach, integrating first the tracer data
and using the resulting ensemble of model parameters to initialize the integration of production
data. The ensemble size is 100, and all the parameters are drawn from uniform PDFs, with a
specified interval for each stage (Table 1). The estimation results for fracture and matrix
properties after assimilating tracer and production data are shown in Figure 11. In Figure 11a,
the initial (theoretical) distributions of the parameters are shown with a solid box. The
histograms of the updated ensembles after integrating the tracer and production data are shown in
black and blue colors, respectively. The results are shown in Stage 1 fracturing (similar results
are obtained for other three stages). We note that the conductivity decline rate is only included in
assimilating the production data. The updated ensemble after tracer data assimilation is improved
and shows significant variance reduction only for fracture conductivity and natural fracture
density. These outcomes are consistent with the results of our sensitivity analysis. The final

51
histogram after integration of production data, however, shows a significant reduction in
uncertainty and closer estimates to the reference case for all parameters. The corresponding
predicted tracer profile and production data, for the initial and final ensembles, also show a
significant reduction in uncertainty after data assimilation. We note that the stochastic nature of
the DFN results in an inherent level of variance, which cannot be reduced during data
assimilation and is reflected in the production forecasts. These results indicate the
complementary nature of tracer and production data in estimating fracture and matrix properties
in the case of DFN and the promise of EnKF for characterization of hydraulic fracturing in
formations with natural fractures.

52
CHAPTER 3
Continuous Parameterization of Discrete Geological Facies
In this chapter, WEintroduce a novel distance transform approach as a continuous
parameterization method, which are subsequently updated using ensemble-based data
assimilation. Parameterization methods for discrete models involve transformation to continuous
domains and provide exact discrete representations after back transformation. In addition to
converting the parameters from discrete to continuous, another goal of some parameterization
methods is to create nearly multi-Gaussian parameters that can be updated using effective multi-
Gaussian estimation methods. The ensemble Kalman filter (EnKF) [Evensen, 1994, 2007] is an
especially popular example of such methods. The EnKF is a Monte-Carlo-based approximation
of the classical Kalman filter, which is an optimal filtering method for linear-Gaussian state
space models. The EnKF is simple to implement, treats the forward model as a black box, and
has been widely used for data assimilation in nonlinear dynamical systems [Evensen and van
Leeuwen, 1996, Houtekamer et al., 2001; Reichle et al., 2002; Bertino et al., 2003; Moradkhani
et al., 2005; Naevdal et al., 2005; Chen and Zhang, 2006; Aanonsen et al., 2009; Jafarpour and
McLaughlin, 2009; Jafarpour and Tarrahi, 2011]. While the EnKF has its well-documented
limitations, including sampling errors and sub-optimality for nonlinear non-Gaussian problems,
it is recognized as a practical estimation method in applications with weakly to moderately
nonlinear models and nearly multi-Gaussian problems. We also show the result of proposed
method for characterization of the rock facies map with two and multiple facies types (e.g. shale
and sand) and in the next chapter, WEpresent the updated discrete map for rock geomechanical
properties using distance transform technique.

53
3.1 Problem Description
Effective development of subsurface hydrological and energy resources hinges on accurate
description and prediction of fluid flow and transport processes that take place in complex
geologic environments. A complicating factor in predicting these processes is the complex
heterogeneity in rock physical properties that directly impact the fluid flow and transport
behavior [Carle et al., 1998; Zinn and Harvey, 2003; de Marsily et al., 2005]. Heterogeneity
exhibits itself across multiple scales, from pore to basin scale. In the face of data scarcity,
identification of the heterogeneity in hydraulic properties of geologic formations at all scales
presents a formidable challenge to subsurface flow modeling community. Geostatistical
modeling techniques have traditionally been adopted to represent the heterogeneity in geologic
formations [Deutsch and Journel, 1998]. Classical geostatistical modeling approaches use
simplifying assumptions, which are often driven by the ease of mathematical treatment (e.g., use
of variogram or covariance functions), to describe the natural variability in geologic formations.
In some cases, the adopted statistical/mathematical models can lead to over-simplification of the
existing geologic features and their connectivity. Together with data scarcity, imperfect
modeling assumptions introduce significant uncertainty into fluid flow and transport predictions
[Carle et al., 1998; Zinn and Harvey, 2003; de Marsily et al., 2005; Hill and Tiedeman, 2007;
Oliver et al., 2008].

A common practice to improve the predictive power of subsurface flow models is their
calibration against observed monitoring and performance data [e.g., Carrera and Neuman, 1986a-
b; Dai and Samper, 2004; Hill and Tiedeman, 2007; Yeh and Khaleel, 2008; Oliver et al., 2008;
Zhou et al., 2014]. In general, the resolution of dynamic flow and transport data seldom lends

54
itself to characterization of heterogeneity at small scales, often necessitating a coarse-scale
description of rock properties for model calibration. Examples of coarse-scale model parameter
representations include upscaling/zonation, low-rank parameterization and regularization
techniques (see e.g., [Jacquard and Jain, 1965; Jahns, 1966; Tikhonov and Arsenin, 1977;
Gavalas et al., 1976; Shah et al., 1978; Lawson and Hanson, 1995; Weiss and Smith, 1998;
Chavent and Bissell, 1998; Grimstad et al.., 2004; Linden et al., 2005; Tonkin and Doherty,
2005; Doherty and Skahill, 2006; Jafarpour and McLaughlin, 2009]). One of the key issues in
formulating model calibration problems is reconciling the discrepancies between model and data
resolutions so that the parameters of interest properly represent the salient and relevant
information (to fluid flow and transport processes) while enabling the estimation of parameters
from available (often low-resolution) pressure and flow rate measurements. While choosing
optimal parameter resolution is important and challenging from an inverse modeling perspective,
the main objective of model calibration is to improve the prediction of flow and transport
behavior.
An important property in predicting the fluid flow and transport behavior in a geologic formation
is the connectivity of the extreme (low and high) values in the hydraulic conductivity field. In
practice, significant uncertainty exists in describing rock facies distributions, especially away
from observations points, which typically consist of scattered injection, extraction, and
monitoring wells. A common way to reduce the uncertainty in facies distribution is through
integration of dynamic response data that often encode important information about interwell
flow connectivity [Hill and Tiedeman, 2007; Oliver et al., 2008; Oliver and Chen, 2011; Zhou et
al., 2014].

55
3.2 Data Assimilation Technique (Ensemble Smoother)
The classical Kalman filter [Kalman, 1960] is a sequential state estimation method for
characterization of the first and second order statistical moments of the state variables’ posterior
distribution. Therefore, for linear and multi-Gaussian state-space systems, Kalman filter
completely characterizes the posterior distribution [Kalman, 1960; Gelb, 1974]. For nonlinear
models and/or non-Gaussian systems, the ensemble version of Kalman filter (i.e., the EnKF)
provides an efficient approximation for practical applications [Evensen, 2007, 2009]. The EnKF
consists of a forecast step in which the mean and covariance of the states are predicted by
implementing a linear state propagation model; and an analysis step that updates the mean and
covariance of the states using the dynamic observations and the mean and covariance of the
forecast states.
The Kalman smoother is also used to update states variables based on future (retrospectively) as
well as past and current observations. The ensemble smoother (ES) was also introduced for
practical data assimilation in 1996 [van Leeuwen and Evensen, 1996]. Recently, an iterative
form of the ES, known as the ensemble smoother with multiple data assimilation (ES-MDA),
was introduced by Emerick et al. [2013] and applied to subsurface flow model calibration
problem. The ES-MDA is implemented with a damping mechanism, which is equivalent to
observation covariance inflation, to avoid sudden updates and ensemble variance
underestimation. In the EnKF, the updates at different times are performed sequentially, whereas
in the ES-MDA all the forecasts are collected in a long vector to perform a single update step. In
an iterative ES-MDA formulation, a single update is repeatedly performed in the hope of
improving the solution. The ES-MDA formulation can be summarized as (see Appendix A for
notation):

56
𝑥 𝑗 = 𝑓 (𝑥 𝑗 , 𝑚 𝑗 ) , 𝑦 𝑗 = 𝑔 (𝑥 𝑗 ) + 𝘀 𝑗 (15)
𝑧 𝑗 = [
𝑚 𝑗 𝑔 𝑗 ] , 𝑧 𝑢 𝑗 = 𝑧 𝑓 𝑗 + 𝛼𝐾 (𝑦 𝑗 − 𝐻 𝑧 𝑓 𝑗 ) 𝐾 = 𝐶 𝑧 𝑒 𝐻 𝑇 (𝐻 𝐶 𝑧 𝑒 𝐻 𝑇 + 𝐶 𝑦 )
−1
(16)

where in the analysis step (Eq. 16), the update vector consists of the states (or augmented states)
forecasts at all observation times. The forecast and updates are repeated multiple times, typically
until the data match falls below a certain threshold or a maximum number of iterations is
reached. The main steps of the ES-MDA, can be summarized as follow:
1. Choose number of data assimilations (Nd) and the corresponding damping factor (α).
2. Run the ensemble forecast to predict all observations with the current parameters.
3. For each realization, the observation is perturbed as follows
𝑦 𝑗 = 𝑦 + √𝛼 𝑖 𝐶 𝑦 1/2
𝑍 𝑑 , 𝑗 = 1,2, … , 𝑁 𝑒 , 𝑖 = 1,2, …, Nd, 𝑍 𝑑 ~𝑁 (0, 𝐼 𝑁𝑑
)
4. Update the ensemble using Eq. 4 after replacing 𝐶 𝑦 with 𝛼 𝐶 𝑦
It is important to note that a difficulty in ensemble data assimilation is due to the finite ensemble
size in practical settings, which leads to sampling errors and inaccurate (spurious) updates. In
general, smaller ensemble sizes are known to result in spurious (overestimated) numerical
correlations that can significantly underestimate the ensemble variance. Two general approaches
to alleviate the impact of spurious updates are covariance inflation and localization (and local
analysis) [Evenson, 2007; Anderson et al., 2007]. A third approach is to increase the ensemble
size by using fast surrogate models. However, the latter approach introduces errors associated
with forecast approximations, which could be difficult to quantify. Recently, [Tarrahi et al.,
2016] introduced a fast linear (approximate) forecast scheme that enables the use of large

57
ensemble sizes to reduce the effect of spurious correlations and to improve the EnKF
performance. In the examples of the next section, observations are perturbed using an
uncorrelated Gaussian noise model with a diagonal covariance matrix C
d
. The standard deviation
of the noise is chosen to be 20% of the observed values.

3.2.1 Distance Transform
In this study, the unknown of interest is the spatial distribution of different facies types in the
formation. It is assumed that the hydraulic conductivity and porosity for each facies type are
constant and known a priori and the unknown parameter is the distribution of different facies
types in the model. Therefore, when the facies type at each cell location is identified, the related
values of permeability and porosity are assigned to that cell. We use a distance transform to
characterize discrete facies distribution as continuous map of distances. Formally, the distance
transform of a function 𝑓 defined on a regular grid 𝒢 , that is
𝑓 : 𝒢 → ℝ, is defined as
𝐷 𝑓 (𝑝 ) = min
𝑞 ∈ 𝒢
(𝑑 (𝑝 , 𝑞 ) + 𝑓 (𝑞 )) (17)

where 𝑑 (𝑝 , 𝑞 ) is some measure of the distance between 𝑝 and 𝑞 . For each point 𝑝 , the above
distance transform definition finds a point 𝑞 that is close to 𝑝 , and for which 𝑓 (𝑞 ) is small. The
above definition is closely related to the traditional distance transform of a set of points on a grid
𝑃 ⊆ 𝒢 , that is
𝐷 𝑃 (𝑝 ) = min
𝑞 ∈ 𝑃 𝑑 (𝑝 , 𝑞 ) (18)

58
where for each location on the grid 𝐷 𝑃 (𝑝 ) defines the distance to nearest point on the grid (see
[Felzenszwalb et al., 2004] for a detailed discussion). In our application, a distance transform of
a binary facies image at each pixel specifies the distance from that pixel to the nearest pixel with
a different facies type. The proposed transform is relatively simple, easy, and fast to implement
[Felzenszwalb et al., 2004]. The generated distance map has only positive values as they
represent Euclidean distances (see Figure 12). A back-transform to recover the original image is
carried out by specifying a threshold distance that locates the facies boundary. This distance can
be adjusted based on the transformed representation of the prior model realizations as discussed
in the sections discussing the numerical experiments. From Figure 12, it can be seen that the
back transformation can reconstruct the original map perfectly. This indicates that, during model
calibration, if the updated distance map detects the presence of a facies feature, the back
transformation should be able to recover that feature without any errors.

Multifacies maps can be considered as a combination of several two-facies maps. To describe the
facies distribution in distance domain, we use multiple distance functions, each corresponding to
two facies. For instance, for three facies, two distance functions are needed (in general for n
facies types, 2
𝑛 −2
distance functions are required). For three facies case, we first combine facies
1 and 2 as a new facies A, and then combine facies 2 and 3 as another facies B. Hence, we have
two sets of two facies maps which are facies A and 3, facies B and 1. The variables D1 and D2
represent the first and second distance functions, respectively while T1 and T2 denotes the
thresholds for the inverse transform of D1 and D2, respectively. At each grid block, if D1 < T1,
then 1 is assigned to that grid block and the P1 binary map is generated. If D2 < T2, 1 is assigned
to the grid block and we generate a P2 map.

59
The three facies types are then designated using the following definitions:
{
𝑃 1 = 1 𝑎𝑛𝑑 𝑃 2 = 1 𝐹𝑎𝑐𝑖𝑒𝑠 1
𝑃 1 = 1 𝑎𝑛𝑑 𝑃 2 = 0 𝐹𝑎𝑐𝑖𝑒𝑠 2
𝑃 1 = 0 𝐹𝑎𝑐𝑖𝑒𝑠 3

If P1=1, the point belongs to set A (facies 1 or 2). If P2 is zero at a point, then that point belongs
to set B (facies 2 or 3). Hence, the intersection of these two sets identifies the desired facies,
which is facies 2.

3.2.2 SVD-Parameterized Distance Transform
To better preserve the prior connectivity in the distance maps, we use a low-rank singular value
decomposition (SVD) parameterization of the prior model realizations. The SVD
parameterization uses a finite number of the leading left singular vectors of the matrix that has as
its columns the vectorized form of each prior model realization [Lawson and Hanson, 1995;
Tonkin and Doherty, 2005; Tavakoli and Reynolds, 2009]. During model calibration, instead of
updating the distance values, the coefficients of the SVD parameterization corresponding to the
distance map for each realization are updated. Due to orthogonality of the SVD basis, the
updated distance maps can be easily recovered through multiplication by the transpose of the
reduced SVD basis. The workflow for implementing the distance transformed followed by SVD
parameterization is shown in Figure 13. In the case of three facies, first two binary images are
generated based on the distribution of the three facies. The two distance maps indicate the
distances from medium and high conductivity facies boundaries. The resulting distance maps are
then parameterized using the SVD parameterization method. In the next step, the model
calibration method (the ensemble smoother with multiple data assimilation steps) is used to
update SVD coefficients of the two distance maps. The calibrated facies models are then

60
obtained through an inverse SVD transformation followed by a back transformation of the
distance maps to facies domain. When SVD is used for low-rank representation in the examples
of this paper, the truncated SVD basis contains S=70 leading left singular vectors, which are
generated using 5000 initial realizations of the respective training data (facies models or distance
maps).

3.2.3 Implementation Steps
Figure 13 shows the workflow for implementing the distance transform for model calibration.
As discussed above, the proposed framework can be applied with and without SVD
parameterization, which is also shown in Figure 13. The SVD parameterization is useful when
the facies models are expected to exhibit smooth and large-scale connectivity patterns. When the
facies do not exhibit smooth and continuous boundaries, the use of distance transform without
SVD parameterization is more appropriate. The calibration methods with and without SVD
parameterization consists of the following main steps (see Figure 13):
Without SVD parameterization of distance maps:
1) Transform discrete geological facies into continuous distance maps as new parameters
2) Use data assimilation (e.g., ES-MDA) to update the ensemble of distance maps
3) Perform inverse distance transform to convert the updated the ensemble of distance
maps to an ensemble of calibrated facies models

61
With SVD parameterization of distance maps:
1) Transform discrete geological facies into continuous distance maps
2) Apply SVD parameterization to the distance maps to obtain near-Gaussian SVD
coefficients as new parameters
3) Use data assimilation (e.g., ES-MDA) to update the SVD coefficients from Step 2
4) Back transform the updated SVD coefficients to updated ensemble of distance maps
5) Perform inverse distance transform to convert the updated distance maps to calibrated
facies models

An important implication of using a low-rank SVD basis to parameterize the distance maps is the
introduction of approximation error due to the truncation of insignificant basis elements. The
resulting error in the distance maps may affect the facies description after back transformation.
More specifically, the approximation error in the distance map may lead to displacement of
facies boundaries and results in local errors in their detection. One way to compensate, at least
partly, for this error is to use the prior model realizations and their distance transforms as training
data to learn this error. A simple approach that we have used in this paper is to use the prior
training data to learn an appropriate threshold that is used to back transform the distance maps to
facies distribution. To this end, using 5000 prior training facies, for different number of retained
SVD basis elements, different thresholds are applied to their SVD-based distance transform
approximation. For each threshold, the back transformations are performed and the resulting
absolute errors (magnitude of the difference between original facies maps and the corresponding
reconstructed map) for all 5000 realizations are computed. Figure 14 displays the average error
resulting from the SVD approximation with different number of leading left singular vectors and

62
for different thresholds used in inverse distance transformation. For different number of leading
SVD basis elements used, the corresponding threshold with the lowest absolute error can be
selected from this figure.

3.3 Experiments and Results
In the first example, the model domain is discretized into a 64 × 64 × 1 uniform grid system and
each grid cell has a dimension of 10 m × 10 m × 1 m. Initially, the contaminant is spread within
the center of the domain (Figure 15, middle column). To clean up the contamination, water is
injected into the surrounding injection wells at a constant rate of 5 m3/day and for a period of
800 days. A commercial simulator [Eclipse 2014] is used to solve the governing equation for a
standard single-phase transient flow equation in water-saturated aquifers. The contaminant
transport is modeled by solving the following governing equation:
𝑑 𝑑𝑡
(
𝑉𝑆𝐶
𝐵 ) = ∑ [
𝑇 𝑘 𝑟 𝐵 µ
(𝜌𝑔𝛿 ℎ − 𝜌𝑔 𝐷 𝑧 )𝐶 ] + 𝑄𝐶 (19)

where advection is assumed to be the main transport mechanism (see Appendix for notations).
Figure 15 shows two well configurations, in which I stands for Injector and P represents a
pumping well (top and bottom rows). The top row refers to Cases 1 and 2 while the bottom row
is used in Cases 3 and 4, all related to Example 1. The only difference in these two
configurations is that in Cases 1 and 2 (top row in Figure 15), the well pattern is placed such that
Well-I4 falls inside the isolated high-permeability facies pattern on the top-left corner. For Cases
3 and 4 (bottom row in Figure 15) the well-pattern is kept the same, but it is shifted upward so

63
that Well-I4 falls outside the disconnected high-permeability facies pattern on the top-left corner.
The main distinction between these two well configurations is that in Cases 3 and 4, in which
Well-I4 does not intersect the high-permeability facies pattern, the responses from all the wells
show no sensitivity to the presence or absence of the disconnected facies in that region. This was
confirmed by removing the isolated high-permeability facies from the model and examining the
responses of all the wells (results not shown). Therefore, these two configurations correspond to
when observed data is sensitive to presence of the isolated facies (Cases 1 and 2), and when the
observed data is not sensitive to presence of the isolated facies (Cases 3 and 4). Another
important consideration in solving model calibration problems is the quality and
representativeness of the prior information. In Example 1, we use the prior model (initial
ensemble), that is statistically consistent with the reference model (the reference map comes
from the same distribution). However, we study two sets of initial ensembles, one that is
conditioned on the hard data (facies types) at the 9 well locations, and another without
conditioning on hard data. Taken together, in Example 1, the following four cases are considered
for model calibration:

Case 1: High Data Sensitivity and Conditional Prior: In this case, the well configurations
are such that Well-I4 is placed in the high-permeability facies and the data
collected from it show sensitivity to the disconnected facies feature on the top-left
corner. The initial ensemble is conditioned on hard data and has the correct
connectivity structure.

64
Case 2: High Data Sensitivity and Unconditional Prior: This case is similar to Case 1,
except that the initial ensemble is not conditioned on hard data, and hence is not
informed about the isolated facies on the top left.
Case 3: Low Data Sensitivity and Conditional Prior: For this case, Well-I4 is not placed
in the high-permeability facies on the top-left corner. The data recorded from all
the wells do not show sensitivity to the presence of the high-permeability facies
on the top-left corner. In this case, the initial ensemble is conditioned on the hard
data, which does not include the isolated high-permeability facies.
Case 4: Low Data Sensitivity and Unconditional Prior: This case is similar to Case 3 with
the difference that the initial ensemble is not conditioned on hard data. In that
case, some of the initial realizations may include the disconnected high-
permeability facies on the top-left corner.

In all the experiments presented in this paper, the spatial distribution of facies is considered as
the uncertain parameter. The permeability and porosity of each facies type are assumed known
and used to distinguish between the facies (permeability = 100 mD and porosity=0.30 for high-
permeability facies; permeability = 50 mD and porosity=0.22 for medium –permeability facies;
and permeability = 0.10 mD and porosity=0.15 for low-permeability facies). For two facies
examples, the high and low values are considered. Measurements of pressure head and pollutant
concentration at each observation well are recorded and used for model calibration. In all
experiments, the following four data assimilation schemes are considered:

65
(i) ES algorithm (using the ES-MDA implementation) to update the log-permeability
(ii) ES algorithm to update continuous distance maps (ES-DT)
(iii) ES algorithm to update the SVD parameterization of log-permeability (ES-SVD)
(iv) ES algorithm to update the SVD parameterization of distance maps (ES-SVD-DT)

Finally, to evaluate the performance of each data assimilation approach, we use ensemble-based
Root-Mean-Squared-Error (RMSE), which is defined for vector variable y as:
𝑅𝑀𝑆𝐸 (𝑦 ) = √
1
𝑁 𝑒 ∑ (𝑦 𝑗 − 𝑦 𝑡𝑟𝑢𝑒 )
𝑁 𝑒 𝑗 =1
2
(20)

3.3.1 Estimation of Facies Distribution
Results (Example 1)
The facies estimation results for all four methods (Columns 2-5) and in Cases 1-4 (Rows 1-4) are
displayed in Figure 16, using a random sample (top plots) and the mean (bottom plots) of the
ensemble. It can be seen from Figure 16 that only for Cases 1 and 2, where data sensitivity to the
high-permeability facies in the top-left corner is strong, the estimation results recover this facies
pattern. Note that none of the methods in Cases 3 and 4 is able to detect this facies feature,
suggesting that regardless of the method used the information in the data does not allow for
estimation this feature. When the facies feature is detected in Cases 1 and 2, the distance is able
to recover it in its discrete form. For Case 1, the isolated small-scale facies pattern on the top-left
corner can be detected and the distance transform can reconstruct it. However, in Case 2, where
conditioning on hard data is not performed, the isolated facies is detected, but it shows

66
connectivity to the facies pattern on the lower-left part of the domain. An interesting observation
is that in Case 3, where Well-I4 is outside the facies, conditioning on the hard data can be
misleading as it prevents almost all of the initial ensemble members (see the mean plot in Figure
16) to include the isolated facies feature. This highlights the complexity of the estimating facies
distribution where sudden changes in facies types can occur within a short distance.

The results from using the regular ES method show the general trend and location of high-
permeability features without reproducing the expected connectivity and the discrete nature of
facies maps. The estimation results with distance transform produce discrete maps and the
obtained solutions tend to exhibit discontinuity and small-scale variability, which is consistent
with the initial ensemble. This indicates that the direct use of distance transform with ensemble
smoother tends to preserve the variability in the prior model. We note that the SVD
parameterization has two important impacts: first, the resulting parameters are better
approximated with a Gaussian distribution; second, by removing high-frequency spectral
components, the parameterization tends to better represent large-scale features and hence
improve facies connectivity. This can be seen in the results from Figure 16. In all 4 cases, when
the SVD parameterization is applied to the facies models directly, the solution shows a smooth
behavior, which is expected. This smoothness is also observed in the estimation results using the
ES-DT-SVD scheme. In summary, the results in Figure 5 show that when data does not show
sensitivity to facies distribution (e.g., a well does not intersect a facies feature), none of the
implemented methods was able to capture the facies feature of interest in that region. However,
when the data has important information about the presence of the facies pattern, the discrete
transform-based methods are better able to identify and represent those features. The results also

67
show that the initial ensemble plays an important role in increasing the accuracy of the updates
and identifying the lack of connectivity in the isolated facies.

Figure 17 summarizes the RMSE plots for both log-permeability distributions (top) and the data
predictions (bottom) for all cases. The data prediction plots include both pressure and
concentration data. The RMSE plots are plotted against the number of ES-MDA iterations. The
results show that in all cases, the ES-DT-SVD scheme provides the lowest RMSE for facies
identification and data prediction. Figures 18a and 18b display the pressure head and pollutant
concentration plots, respectively, before and after updating with the ES-DT-SVD method for all
the wells and in all four Cases. In this figure, the red line represents the predictions with the
reference model while the grey and cyan curves show the predictions before and after updates.
The prediction ensemble means with the initial and updated models are shown with a dark black
and dark blue colors, respectively. As it can be seen, after data assimilation the predicted
ensemble of pressure head and concentration are closer to the observed data and the uncertainty
in these quantities is significantly reduced. In general, the updated models provide a less
uncertain and more unbiased forecast of the pressure heads and contaminant concentrations.

3.2 Example 2: Three-Facies Model
Experimental Setup
In Example 2, we consider a three facies model consisting of facies with low, medium, and high
permeability values. The reference model and well configuration for this example are displayed

68
in Figure 19 (top row). In this case, the estimation is performed to update two different maps
that correspond to distances from two facies boundaries as discussed in the methodology section.
The reference model includes three facies types with their major continuity aligned diagonally in
the north-west and south-east direction. The facies with highest permeability is located on the top
right while those with medium and low permeability are in the center and lower-left parts of the
domain, respectively. The initial and final pollutant concentrations are shown in Figure 19 (top
row, left).

Results (Example 2)
The plots on the bottom of Figure 19 contain the estimation results for the log-permeability
(facies) distribution. Two sample realization, the ensemble mean and variance are shown in the
first to fourth rows. The initial facies maps are shown in the first column and indicate that, while
having consistent orientation, the three facies types are located in different areas of the domain.
The updated maps from different data assimilation techniques are shown in columns 2-5 of
Figure 19 (bottom). The mean maps of the update models from the ES and ES-DT methods
show that the calibrated models can capture the locations of facies with low and medium
permeability, and to some extent the facies with high permeability, in the reference model.
However, the results from the ES method are continuous and do not show the presence of three
facies types. Furthermore, the variance of the updated facies maps from the ES method seems to
be severely underestimated. While applying the SVD parameterization before the updates
improve connectivity of geologic features, the ES-SVD method still results in continuous
solutions that do not reveal the three facies types. On the other hand, estimation with the ES-DT-

69
SVD method provides updates that include three facies types, which are located in areas
consistent with those in the reference model. Figure 20 depicts the ensemble prediction of
pressure (top) and concentration (middle) data before and after the update for all wells using the
ES-DT-SVD method. As expected, the results show improved data match and forecast after the
update. The log-permeability and data prediction RMSE plots are also shown on the bottom in
Figure 9. The two RMSE maps reveal the lowest error in updating the facies and predicting the
data to be associated with the ES-DT-SVD method.

3.3.2 Estimation of Facies Distribution and Properties
3.3 Example 3: Estimation of Facies Distribution and Hydraulic Conductivity
In our previous examples, it was assumed that the values of facies hydraulic conductivity are
known and we are only interested in identifying the facies distribution in space. The underlying
assumption in those examples is that the uncertainty and impact of the spatial distribution of rock
types are far greater than those of facies values. Assuming that rock types have very different
(and disjoint) ranges of hydraulic conductivity values and that well log and core analysis data
can provide some information about the roc types, it is possible to approximately estimate the
value of hydraulic conductivity for each rock type. Furthermore, our sensitivity analysis (Figure
21), also indicates that, when the facies boundaries are correctly specified, a modest change in
the facies hydraulic conductivity does not have a significant impact on the flow response.
However, we observed some sensitivity to the transport response for large changes in facies
hydraulic conductivity.

70
A logical extension of the proposed method is to assume that the facies hydraulic conductivity
values are unknown and need to be estimated from dynamic flow and concentration data. With
the EnKF estimation framework, adding the hydraulic conductivity of facies to the unknown
parameters is straightforward and can be achieved by augmenting the unknown parameters with
the facies values, as discussed in the following section.

Experimental Setup
Example 3 is identical to Example 2 with three facies types. However, in this case, it is assumed
that the facies hydraulic conductivity values are uncertain. The estimation is performed to
simultaneously update two different spatial maps that correspond to distances from two facies
boundaries as discussed in the methodology section, and three permeability values each
corresponding to a facies type. The results of sensitivity analysis of pressure head and pollutant
concentration with respect to the value of permeability for each facies type are depicted in
Figure 21. Figure 21a and 21b show the logarithms of the mean pollutant concentration and
pressure head over time, respectively. Each height of each colored box indicates the variation in
model response as a result of the changes made to facies permeability. We introduced 10%, 20%
and 40% changes in the value of permeability for each facies type. And the corresponding
differences between the observed responses with those from the reference map are shown in the
vertical axis for all three facies (low, medium and high permeability). The maximum change in
the value of permeability is chosen to introduce significant changes without changing the facies
type. For example, a reduction of 50% in the value of facies with high permeability results in a
permeability value of 50 mD, which corresponds to the facies type with medium permeability.

71
The height of each box in Figure 10 represents the well response for each case while the change
in the height of each box indicates the sensitivity of the well response to that parameter. As can
be seen from these figures, pressure head in all wells does not show sensitivity to the value of
permeability (i.e., the boxes with different colors have similar heights); however, the pollutant
concentration observed in some wells shows sensitivity to the value of permeability. This is
expected as a significant change in the permeability value affects the velocity values and the
resulting pollutant arrival times. In general, the observed sensitivities are related to the location
of the initial pollutant source, well configuration, and facies distribution between injection and
pumping wells. For example, the concentrations recorded in Well P6 show sensitivity to the
value of permeability for facies with low and high permeability as P6 and I2 are located within
the facies with high permeability (Figure 10). Also, the concentrations in Well P1 is sensitive to
the value of permeability of facies with low permeability as I1 and P1 are both located within
facies with low permeability (Figure 10). Moreover, Well P1 shows sensitivity to the value of
permeability for the facies with high permeability only after a significant change of 40% which
is because of the effect of I2 and location of well P1 close the boundary of facies.

Results (Example 3)
The plots on the Figure 22 contain the estimation results for the facies distribution as described
in Section 3.2 for Example 2. The updated maps from different data assimilation techniques are
shown in Columns 2-5 of Figure 22. The mean maps of the update models from the ES shows
that the calibrated model can capture the location of facies with low and medium permeability
and to some extent the facies with high permeability. However, the results from the ES method

72
are continuous and do not show the presence of three distinct facies types. Furthermore, the
variance of the updated facies maps from the ES method seems to be severely underestimated,
suggesting that the ensemble has collapsed. The calibrated model using the ES-DT method can
capture the location of each facies type in the reference model and doesn’t suffer from an
ensemble collapse issue. While applying the SVD parameterization before the updates improves
the connectivity of geologic features, the ES-SVD method still results in continuous solutions
that do not reveal the three facies types. On the other hand, estimation with the ES-DT-SVD
method provides updates that include three facies types that located in areas consistent with those
in the reference model. Figure 23 shows the updated value of permeability for each facies type.
In this figure, the red crosses represent the reference value and the cyan bars represent the
ensemble variance. As it was observed in our sensitivity analysis, some wells are only sensitive
to the value of high and low permeability (the extreme values), which partly explains why the
value of permeability in the medium-permeability facies doesn’t change significantly. However,
the estimated permeability values in facies with low and high permeability exhibit marked
improvements.

Figure 24 depicts the ensemble prediction of pressure (top) and concentration (middle) data
before and after the update for all wells using the ES-DT-SVD method. The results show
improved data match and forecast after the update step. The log-permeability and data prediction
RMSE plots are also shown in Figure 24 (bottom). The two RMSE maps correspond to the
facies permeability and the predicted data and show that the error associated with the ES-DT-
SVD method is the smallest.

73
3.4 Discussion and Future Work
Estimation of facies models offers a way to identify the first order approximation of rock
property distributions without including fine details that represent within facies heterogeneity.
When the hydraulic properties of different rock types exhibit sharp contrasts and the within
facies variability is negligible, the facies dominate the parameter field tend to dominate the flow
and transport behavior. In those circumstances, locating the facies boundaries becomes the main
problem for model calibration as identification of the facies spatial distribution contribute
significantly to constraining the flow and transport behavior in the geologic formation. In
addition to the importance of facies spatial distribution in describing fluid movement in geologic
formations, the resolution of flow and transport data from scattered well locations limits the
resolution and details that can be included in a numerical model. Therefore, one could view the
estimation of facies as an additional constraint in model calibration where fine-scale detailed has
been eliminated. In fact, regularization methods such as the total variation regularization promote
piecewise constant distributions, which have been shown to be useful for estimating facies
boundaries [Lee and Kitanidis, 2013].

Estimation of facies distribution is partly complicated by the discrete indicators that are assigned
to each rock type. The proposed method in this paper offers a simple transformation that can
parameterize discrete facies types into continuous distance maps that represent, for each grid
cell, the distance from the nearest facies boundaries. The resulting distance maps can then be
updated based on dynamic flow data and using methods that are developed for continuous
parameter estimation (i.e., the ensemble smoother in this paper). In estimating facies distribution
one may assume that once a facies type is identified the corresponding hydraulic property of the

74
facies are known. While observations of facies types and their hydraulic properties may be
available at well locations, assuming known hydraulic properties for facies may introduce an
approximation, which could be small depending on the available data and the contrast in
hydraulic properties of different rocks. When information about the hydraulic properties of
different rock types is not available, those properties should be considered as unknown
parameters that need to be estimated along with facies distribution. Our Kalman filter based
model updating approach offers a convenient way to augment the facies distribution (i.e.,
unknown parameter vector) with their hydraulic properties during the data integration step.

It is important to point out that when the contrast in hydraulic properties of different rocks is not
significant the facies boundaries play a less important role in governing the flow and transport
behavior. Therefore, in those circumstances, using facies distributions to describe a model may
not be effective as the available data is likely to offer little sensitivity to the locations of facies
boundaries. Furthermore, within facies heterogeneity can become as important as the facies
distribution and their boundaries when the contrast in hydraulic properties of different rocks is
not very significant, which also implies that a facies-based description is not likely to be
effective. As in other facies model calibration techniques, such as the level set method, when the
number of facies types increases, so does the complexity of the proposed method. Specifically, in
our method, the number of facies types dictates the number of parameters that are estimated (i.e.,
distance maps that are generated through transformation). For instance, when three facies types
are present, two distance maps are generated after transformation, and need to be estimated
during inversion. In general, for F facies types (F-1) distance maps of the same dimension as the
size of the numerical model should be generated. However, the SVD-based implementation in

75
this paper (i.e., ES-DT-SVD method) can be used to provide a low-rank representation of the
resulting distance maps.
Another important consideration is the back-transformation from the distance domain to the
original facies domain, which involves thresholding to identify the boundaries of the discrete
facies. In this paper, we learned this thresholding by using the initial ensemble as a training
dataset. More sophisticated methods may be developed to identify what threshold need to be use,
especially when training data is not available to learn a threshold.

76
CHAPTER 4
Integration of Discrete Microseismic Data
The objective of this chapter is to introduce a novel geomechanically-driven modeling approach
for accurate representation and inversion of discrete microseismic monitoring data, as a key
monitoring technology, during bulk fluid (Gas and liquid) injection into geologic formations.
Here, we present the developed approach for application to CO2 injection. To establish the
correlation between rock mechanical properties and microseismic cloud, the developed
framework combines a coupled flow and geomechanical model with Mohr-Coulomb’s failure
criterion to describe the seismicity potential in the reservoir. The seismicity potential map (i.e.,
probability of geomechanically-induced seismic events at different times and locations) will
depend on the stress conditions and mechanical strength of the formation rock. It is then used to
simulate the predicted microseismic clouds for any input model of rock property distributions.

In general, there are two major sources of uncertainty in interpreted rock failure (potential source
of MEQ events), Coupled physics modeling uncertainty and seismic modeling uncertainty. The
first category (coupled physics modeling uncertainty) consists of rock facies distribution,
magnitude of physical rock properties and within facies heterogeneity. To account for the first
source of uncertainty in the estimated failure map and microseismic cloud, we considered
multiple realizations and implemented Monte Carlo based data assimilation technique. Here, we
designed three case studies to analyze the sensitivity of microcseimic data to each of these
variables and we implemented the proposed geomechanics-based approach to reduce the
associated uncertainty. Second source of uncertainty in interpreted microseismic map includes
velocity model and signal to noise ratio (SNR), in general, uncertainty in MEQ event detection.

77
We considered an uncertain region around the event location to include the effect of them in our
analysis.

If the simulated microseismic clouds differ from those observed in the field, the flow and
geomechanical model input parameters are updated to match the observed microseismic clouds.
An ensemble-based method is used to assimilate the observed discrete microseismic events and
estimate rock flow and mechanical properties. we also combined the proposed ensemble-based
method with distance transform (introduced in chapter 3) to estimate the distribution of rock
properties including Young’s modulus, Poisson’s ratio, friction angle and cohesion along with
the hydraulic properties.

4.1 Coupled Flow ad Geomechanical Model
Prior to development of our microseismic data integration method, we should build a forward
model to simulate coupled flow and geomechanical effects. To develop the proposed data
assimilation formulation for inversion of discrete microseismic data, we set up a 3D fully
coupled fluid flow and geomechanics simulation model, which also involves geochemical
reactions to simulate reactions between the formation rock and the injection CO2 (Table 5). The
3D model was developed using the CMG simulation package. The model has 191*191*6 grid
blocks, with a uniform cell size of 20m*20m*30m. Figure 25 shows the model and along with
its dimensions and the assumed boundary conditions.
We have used this 3D model with approximated distribution of fracture network (DFN) as
depicted on Figure 26. In this proxy model, we replaced each cluster of natural fracture elements
with a lumped irregular-shape region. The discrete fracture network is generated by model

78
described in chapter 2 and we used connected volume algorithm in PETREL software to
approximate each group of fracture network with a 3D object (Figure 26-b). As it can be seen in
Figure 26-c, clusters of fracture elements are shown with various colors to distinguish them as a
separate objects in geologic formation. Therefore, the fractured geologic formation can be
approximated with distribution of irregular shape objects with their specific hydraulic and
geomechanical properties. In other word, we assumed a simplified model for dynamic fracture
opening and propagation during the injection process. Since shale rock represents a rock type
with higher brittleness and higher chance to be fractured, we defined the generated irregular
shape objects as a shale rock type (or weak region). Therefore, we can describe the geologic
formation with two rock facies map, shale and sand. In section 4.2, we will show that the proxy
model statistically represents the same distribution of microsesimic cloud as discrete fracture
network in terms of b-value in power’s law.
We considered homogeneous and heterogeneous descriptions of rock property distributions to
perform sensitivity analysis. For homogeneous distribution of rock properties, we have assumed
that all rock properties are constant. For the heterogeneous case, we considered two rock facies
types with spatially stochastic distribution within the formation. We have also considered the
effect of small-scale within-facies heterogeneity for hydraulic properties. For geomechanical
properties, we have assumed facies-dependent constant property values (without including the
variability within each facies). While these assumptions can be easily relaxed if necessary, the
results of our sensitivity analysis are shown in next sections. Figure 27 illustrates a sample
realization of the rock facies types (top) along with two realizations of heterogeneous
petrophysical rock properties for the same facies model. We have generated multiple realizations
of geomechanical rock properties by probabilistically assigning different values of

79
geomechanical properties to each facies type. For petrophysical properties, we have used the
Sequential Gaussian Simulation (SGS) algorithm in PETREL to generate several realizations that
represent the uncertainty in within-facies variability (two samples are sown in Figure 27). Table
6 summarizes the range of values for rock fluid and geomechanical properties, which we have
assigned based on the published studies in the literature [].

The result of our forward simulation model for a 3D geologic formation with two rock facies
type and homogeneous properties is shown in Figure 28. As you can see, the distribution of rock
displacement and stress can be affected by the location of lump objects (shale rock types with
high brittleness shown with dark blue color). In the following sessions, we used one layer of 3D
model to perform sensitivity analysis and develop our proposed data assimilation technique.

4.2 Failure Criterion
To study the sensitivity of potential microseismic data (MEQ) to different rock properties, we
have assumed that microseismic events are related to rock failure. Geomechanical rock failure
can be described by predicting the stress variation due to CO2 injection and establishing the
Morh-Coulomb’s failure criterion as follows:
𝐹 =
𝜎 1
+𝜎 3
2
sin𝜙 𝑓 −
𝜎 1
−𝜎 3
2
+ 𝑐 𝑓 cos𝜙 𝑓 (21)
in which σ 1 and σ3 represent maximum and minimum stresses, respectively, and ϕf and cf denote
rock frication angle and cohesion coefficient. These parameters are shown in Figure 29.

80
The Mohr Coulomb failure criterion depends on important rock properties. Some of these
properties, such as cohesion and friction angle, explicitly appear in the definition of failure
criterion while other properties such as the Young’s Modulus, Poisson’s ratio and permeability
and porosity are included implicitly as they have important impact on the stress state. Assuming
a homogeneous formation with constant rock flow and mechanical properties, Figure 30a, and
30b illustrate the correlation between rock properties and failure map (Middle layer of 3D
model). The result of coupled flow and geomechanical simulation runs for two sample
permeability maps (with high and low permeability values, respectively) are shown in this figure.
The corresponding vertical displacements, minimum/maximum stress, and failure maps are
displayed for each case. In Figure 30, negative values on the failure map indicate potential
locations for microseismic events. As can be seen when the permeability is low (Figure 30b),
the failure map shows a large area around the injection point where failure has occurred.

4.3 Source of Uncertainty in Potential Microseismic Events
Recorded microseismic data are believed to indicate rock failure in the geologic formation.
Microseismic data acquisition can be done by both surface and downhole receiver arrays as
illustrated in Figure 31. Recorded data from surface arrays often have lower quality and signal to
noise ratio (SNR), but provide a higher angle coverage. On the other hand, downhole receivers
can capture smaller microseismic events while having narrower angle coverage. The locations of
microseismic events can be interpreted by inverting the P-wave and S-wave differential arrival
times, using a velocity model. Two major sources of uncertainty that should be accounted in our
modeling work are the uncertainties associated with coupled physics modeling and those related
to seismic modeling. The first category consists of various inputs used for simulating coupled

81
flow and geomechanics (including the description of rock and fluid properties and their
distributions in space, boundary and initial conditions, and assumption used in deriving the
governing equations). In this project, we will focus on rock property distribution as a dominating
source of uncertainty. Figure 32 shows the effects of some of these parameters on the
distribution of pore pressure, minimum and maximum stresses, and the failure map. We observe
that, overall, within facies heterogeneity is not as significant as the heterogeneity in facies
distribution and the value of rock properties in each facies type. For this project, these sources of
uncertainty will be incorporated in our coupled physics simulation.

To account for the both source of uncertainty (Coupled physics modeling and seicmic
uncertainty) in the estimated failure map and microseismic cloud, we propose a probabilistic
approach to simulate microseismic cloud using failure map. This approach will be discussed with
more detail in the next section.

4.4 Microseismic Data Simulation
Before performing the sensitivity analysis, two MEQ simulation models will be introduced,
deterministic and stochastic models, to predict the location of MEQ events (MEQ cloud) using
described failure map in previous sections.
4.4.1 Deterministic MEQ Simulation Model
In a deterministic view, at each grid block, negative F values on the failure map indicates 100%
chance of observing an event while positive F values imply 0% chance of observing events.
However, the estimated Mohr circle and failure envelope can be affected by the uncertainty in
rock properties. For example, if the estimated cohesion (Cf) is higher than what is depicted in

82
Figure 29, the failure envelope will move up and does not intersect with the Mohr circle. In
other words, uncertainty in cohesion and friction angle explicitly affects the rock failure point
through failure line (Figure 33-a). On the other hand, the other rock properties, such as
permeability, porosity, Young’s Modulus and Poison’s ratio, can influence rock failure implicitly
by changing minimum and maximum stresses on the Mohr circle (Figure 33-b).

As a result, the uncertainty in rock flow and geomechanical properties can result in
misinterpretation of microseismic clouds. Therefore, the uncertainty in the estimated failure map
can be modeled stochastically as it will be described below.

4.4.2 Stochastic MEQ Simulation Model
To account for the uncertainty in the estimated failure map and microseismic cloud, we propose
a probabilistic approach to simulate microseismic clouds from the rock failure map. We used
sigmoid function to relate the failure map to the probability of failure, as follows
P =
1
1+exp (𝐹 )
(22)
in which, F is failure value at each grid block and P is the corresponding probability. Figure 34
(left) depicts the shape of the sigmoid function for different values of failures. As it can be seen,
there is a 50% chance of failure when F is zero (the intersection blue point on Figure 2) and
100% chance of failure for very large negative value of F and 0% chance of failure for very large
positive value of F (i.e., no intersection between Mohr Circle and Failure envelope). Figure 34
(right) also illustrates how the uncertainty in rock properties and the calculated stress state on
Mohr-Coulomb’s failure diagram can be interpreted as probabilistically using the sigmoid
function.

83

We implemented the proposed stochastic method on several examples to study its properties. In
the first example, we first consider a simple 2D homogeneous aquifer with constant rock
property (Figure 35). The estimated probability map and its corresponding microseismic cloud
using the above stochastic model indicate that the probability of observing microseismic events
will be decreased by increasing distance from injection point. Accordingly, the microseismic
cloud is more concentrated around the well. In the second example, we considered a more
complex geological structure with two rock facies types. We compared the simulated
microseismic cloud using both a deterministic and stochastic approach. In deterministic
approach, we use F≤0 as a condition for generating microseismicity. In the stochastic approach,
first we generate the probability map as described above and then we sample from the Bernoulli
distribution using the calculated probability map. Each sample from the Bernoulli distribution
represents the distribution of microseismic events. In other words, the stochastic method
generates multiple possible microseismic clouds rather than one single map.

To study the effect of failure map on the simulated microseismic cloud, we show two extreme
cases, one with weaker (i.e., low cohesion) rock, shown as blue regions, and another with a
stronger (high cohesion) rock (Figure 36). When the cohesion of the weak region is very low
(top figure), the generated microseismic maps using both deterministic and stochastic approach
reveal the shape of the weak region. On the other hand, for the case with higher cohesion value
in the weak region, the microseismic response, by both determinist and stochastic methods, does
not show the exact shape of the weak rock. This example shows how the uncertainty in rock

84
properties can impact the simulated microseismic events. However, the examples do not
incorporate the uncertainty in seismic source inversion for locating microseismic events.
The main sources of uncertainty in the interpreted microseismic events distribution include the
velocity model and signal to noise ratio (SNR) of the recorded seismic data. The uncertainty in
both the velocity model and the signal to noise ratio often manifest themselves as uncertainty in
the time and/or location of the events (notice that the velocity model related time to distance). To
model the uncertainty in the event location, we consider an ellipsoid (ellipse in 2D) centered at
the location of the predicted microseismic event. The principal axes of the ellipsoid are
determined based on the level of uncertainty in the velocity model and the recorded seismic data.
Figure 37-a illustrates this approach by placing ellipses around sample microseismic event
locations. In Figure 37-b depicts the simulated microseismic clouds for the two cases in Figure
36. As it can be seen, accounting for the seismic modeling uncertainty results in dispersed
microseismic clouds that do not represent the exact shape and location of the failed regions.

4.5 Sensitivity Analysis
To study the sensitivity of the failure map (i.e., microseismic response) to different rock
properties, we considered 10%, 30% and 50% increase in the values of each rock property
(Table 6). The microseismic response is quantified as the number of cells in the model that meet
the Mohr-Coulomb failure criterion using stochastic approach (The same observation was
reported by implementing deterministic approach). The results of our sensitivity analysis for the
homogeneous case are summarized in Figure 38. In this figure, C, E, F, poi, K and pro stand for
the cohesion, Young’s modulus, friction angle, Poisson’s ratio, permeability and porosity,
respectively. In general, microseismic data show more sensitivity to geomechanical properties

85
compared to rock flow properties. The results show that decreasing the friction angle and
cohesion, increases the number of potential microseismic events. This can be explained by
noting that lower values of cohesion and friction angle are indications of a downward shift in the
Mohr-Coulomb failure line, which increase the chance of intersection the failure line with the
Mohr circle. The number of failure events show positive correlation with the Young’s modulus
and a weak negative correlation with the Poison’s ratio. Moreover, MEQ data are less sensitive
to rock flow properties (permeability and porosity), which is why they are shown in linear scale
(not log-scale). By decreasing the value of permeability and porosity, the number of events
increases primarily due to pressure buildup in the aquifer (increase in pore pressure reduces rock
effective stress).

In the second experiment, we assumed a model with two rock facies type (discrete heterogeneous
map) as depicted on Figure 26, and both rock flow and mechanical properties are constant
within each rock facies (ignoring within-facies variability). Similar to our previous experiment,
we apply 10%, 30% and 50% increases to each rock property value (one at a time) to study their
effect on the simulated failure map (microseismic data). The results of our sensitivity analysis
are shown in Figure 39. In this figure, Es, pois, Ks and pros represent Young’s modulus,
Poisson’s ratio, permeability and porosity for sand, respectively; and C, E, F, poi, K and pro
stand for cohesion, Young’s modulus, friction angle, Poisson’s ratio, permeability and porosity
for shale, respectively. The main observations form this experiment are in agreement with the
results from the previous experiment with homogeneous rock properties.

86
4.6 Proposed Data Assimilation Technique
Standard inverse modeling methods such as ensemble Kalman filter (EnKF) that are designed for
continuous data cannot be used with discrete microseismic events. Therefore, a new inverse
modeling paradigm that respects the discrete nature of microseismic data is needed to extract
rock hydraulic and mechanical properties from these observations. In this work, we propose
transformation techniques to generate the continuous representation of discrete microseismic
cloud based on the location of events and associated uncertainty with it. Then, the continuous
map can be used with ensemble smoother (ES-MDA) to estimate rock properties. The proposed
data assimilation workflow is shown on Figure 40.
Kernel Density Estimation:
Kernel density estimation (KDE) is a method to interpret MEQ events as a continuous
measurement. In this method, spatial MEQ events are replaced with a Gaussian kernel function
centered at the event location. Then, by adding Kernel function in the space, we construct a
continuous function that represents the spatial density of MEQ events. The continuous seismicity
density map can be written as:
𝐷 𝑠 (𝑙 ) =
1
𝑛 𝑒 ∑ 𝐾 𝑖 (𝑙 )
𝑛 𝑒 𝑖 =1
(23)
The 𝐾 𝑖 = 𝑁 (𝑙 , 𝛴 )is a Gaussian kernel, 𝑛 𝑤 is the number of events at each time step. u indicates
the location coordinate of MEQ event. Σ is covariance matrix of the Gaussian kernel. 𝐷 𝑠 (𝑙 ) is a
continuous function which represents the seismicity density at all locations in the geologic
formation. The Gaussian kernel can be represented as:

87
𝐾 (𝑙 ) = 𝑁 (𝑙 𝑀𝐸𝑄 , 𝛴 ) =
1
2𝛱 |𝛴 |
1
2
(−
1
2
(𝑙 − 𝑙 𝑀𝐸𝑄 )
𝑇 𝛴 −1
(𝑙 − 𝑙 𝑀𝐸𝑄 )) (24)
The covariance matrix for the Kernel function (Σ) determines the shape, size and orientation of
the Gaussian ellipsoid centered at the MEQ event location. The value of bandwidth for Kernel
function can be determined by uncertainty in locating seismic source through seismic source
inversion methods. Figure 41 illustrates the continuous representation of discrete MEQ data
using KDE method for one dimensional problem [Tarrahi, 2015]. In this figure, blue crosses
represents the location of discrete data, the transformation of each point is shown with black
Gaussian curve, and red curve represent the sum of kernel function in space (black curves).
Continuous representation of a two-dimensional microseismic cloud using KDE approach is
depicted in Figure 42. As you can see, continuous map shows higher value in the region with
larger number of MEQ events.

4.7 Results
The results of implementing the proposed workflow will be presented with three experiments. In
the first experiment, we assumed that the distribution and boundary of failed (or weak) regions is
known, but there is uncertainty in the value of hydraulic and geomechanical properties. In the
second experiment, the distribution of the weak regions is considered as the main source of
uncertainty. And in the last experiment, the distribution of the weak regions along with the value
of hydraulic and geomechanical properties are estimated.
Experiment 1
In this experiment, we use a two-dimensional geological model with two regions (this model
corresponds to the fourth layer of the 3D model in Figure 28). The reference model is depicted in

88
Figure 43 (top). The weak regions are shown in blue color and are assumed to have high brittleness.
We assume that the spatial distributions of these regions are known, but the value of hydraulic and
geomechanical properties in each facies is uncertain. The assumed range for hydraulic and
geomechanical properties can be found in Table 6. The results of updating rock properties with the
described workflow are summarized in Figure 43 (middle and bottom rows). In this figure, the initial
and updated ensemble of parameters is shown with grey and blue histograms, respectively. The true
value of parameter in the reference model is shown with the red line. We also plotted the rock
properties for two samples to illustrate how they are updated. The initial and updated parameters for
Sample 1 are shown with black asterisk and cyan circle, respectively. For Sample 2, the red asterisk
and green circle illustrate the initial and updated parameters. After update, the variance of the updated
rock properties is reduced, and the updated ensemble gets closer to the reference value (red line).
However, a higher variance is observed for some of the updated parameters, including sand
permeability, sand friction angle and shale Poisson’s ratio compared to other parameters. This
observation can be explained by the results of our sensitivity analysis (not shown). The microseismic
clouds along with their KDE transformation before and after updates are shown in Figure 44. After
updating the rock properties, the distribution of microseismic events is more consistent with those
predicted by the reference model.

Experiment 2
Here, we assume the magnitude of rock properties within each region is constant and estimate
the distribution of the discrete weak regions from microseismic data. We use the distance
transform for continuous parameterization of discrete geological regions along with the KDE
approach for generating continuous representation of the discrete microseismic cloud. The main

89
steps of proposed workflow illustrated in Figure 45. In Figure 45, the top part shows the
continuous parameterization of one initial discrete geological map. The middle part illustrates the
use of KDE for converting discrete microseismic data to continuous observations, and the bottom
part depicts the updated distance map by integrating microseismic data and applying the inverse
distance transformation to retrieve the discrete distribution of geological facies. The results from
applying this workflow to a 2D model are summarized in Figure 46. The comparison between
the initial and updated models with the reference case show improved results in identifying the
weak regions after data assimilation. The resulting predictions of microseismic events
distribution and their corresponding KDE transformation with the discrete geological facies
before and after data assimilation are shown in Figure 47. These results also demonstrate the
improved performance in predicting microseismic event distributions after data assimilation.

Experiment 3
In this experiment, we estimated the distribution of the discrete weak regions and magnitude of
rock properties simultaneously by integration of discrete microseismic data. In fact, this
experiment is the combination of first and second experiments. The results are summarized in
Figure 48 and Figure 49. The comparison between the initial and updated models with the
reference case show improved results in identifying the weak regions and rock properties after
data assimilation. The resulting predictions of microseismic events distribution and their
corresponding KDE transformation before and after data assimilation are shown in Figure 50 as
well. These results also indicate the improved performance in predicting microseismic event
distributions after data assimilation.

90
4.8 Discussion and Future Work
Accurate representation of the distribution and attributes of discrete microseismic monitoring
events is critical for characterization of rock flow and mechanical properties. These properties, in
turn, determine the changes in the subsurface stress and strain distributions due to fluid injection.
The proposed workflow establishes physical correlations among rock mechanical properties and
discrete microseismic events, a fundamental requirement for formulating and solving inverse
problems to estimate rock physical properties from observed microseismic monitoring data. In
this work, hydraulic and geomechanical rock properties were estimated in various scales. The
proposed distance based kernel density estimation approach is used to determine the first order
approximation of rock properties distribution by integration of discrete microseismic data.
Results of the sensitivity analysis and model calibration showed that the effect of fine details that
indicate within facies (or failed region) heterogeneity on distribution of microseismic events is
not significant. Specifically, when the properties of different rock types exhibit sharp contrasts,
the within weak region variability is negligible. In those cases, distribution of weak regions (or
failed region) dominates the stress and strain distribution along with the flow and transport
behavior in the geologic formation.

It is important to point out the role of two major sources of uncertainty in interpreted failure map
(potential source of MEQ events) and relationship among them. Coupled physics modeling
uncertainty consists of rock facies distribution and properties. The uncertainty in the flow and
geomechanical properties can result in misinterpretation of MEQ cloud. Therefore, error in the
estimated failure map can be modeled stochastically using sigmoid function. This function
converts the failure map to the probability map based on the magnitude of rock failure at each

91
grid block. Second source of uncertainty in detection and processing of microseismic map
includes velocity model and signal to noise ratio (SNR). Velocity models are generated by using
sedimentary layers derived from geology. Therefore, existing uncertainty in geological structure
of the reservoir (First source of uncertainty) can also lead to uncertainty in velocity model and
interpreted location of MEQ events. On the other hand, signal to noise ratio (SNR) of recorded
wave depends on the survey design and data acquisition. The uncertainty in both the velocity
model and the signal to noise ratio often manifest themselves as uncertainty in the time and/or
location of the events (notice that the velocity model related time to distance). To model the
uncertainty in the event location, we consider an ellipsoid (ellipse in 2D) centered at the location
of the predicted microseismic event. The principal axes of the ellipsoid are determined based on
the level of uncertainty in the velocity model and the recorded seismic data. By combining
sigmoid function with the random allocation of event location within its vicinity, we developed a
stochastic method to model distribution of microseismic data in the geologic formation.
Characterization of hydraulic and geomechanical rock properties using microseimic data is
complicated due to the assigned discrete indictors to the location of microseismic events. The
proposed model calibration workflow offers a simple transformation that quantify discrete
microseismic events along with their associated uncertainty. Our ensemble-based model
updating approach offers a convenient way to augment the rock facies distribution with their
hydraulic and geomechanical properties during the data integration step with discrete
microseismic data. However, transformation intrinsically involve information loss and
representativeness errors. To enhance the description of discrete microseismic events and
correlate them to rock physical properties, direct filtering approach without transformation of
discrete data to continuous domain is required. Therefore, in the next step of this work, a

92
stochastic point process filtering paradigm can be developed for accurate representation and
inversion of discrete microseismic events. A space-time point process is a random process,
whose realizations define the spatial location and occurrence time of the events. Point processes
are well studied in statistics and probability theory (Kallenberg, O., 1986; Daley, D.J, 1988) and
have many applications in various disciplines, including computational neuroscience,
seismology, materials science, and telecommunications. The theory of point processes is covered
in details in (Daley and Vere-Jones, 2003; Karr, 1986). Vere-Jones (1970) and Kagan 919730
were among the first to study seismicity. While point process modeling of the binary events is
quite natural, it is not typically used in modeling microseismic data. A point process for a
discrete signal is a series of random 0 and 1 events that occur in continuous time (Daley and
Vere-Jones, 1988; Truccolo et al., 2005). At any given location, the 1’s indicate the time point at
which a microseismic event has occurred. A stochastic point process is completely specified by
its conditional intensity (or instantaneous rate) function. Given an observation interval (0,T ] ,
let’s denote by n(t) the number of events that occur in the interval (0, t] for t (0,T ] . The
instantaneous rate function is defined as (Daley and Vere-Jones, 1988).
λ(𝑡 |H (t)) = lim
𝛥 →0
P (𝑛 (𝑡 + Δ ) − n(t) = 1|H (t))
𝛥

where P(. | .) denotes a conditional probability and H (t) denotes all covariates that affect the
event’s firing rate. For microseismic data, these covariates correspond to rock mechanical
properties and failure mechanisms (criteria). Therefore, the probability of having an event in a
small time interval (t, t +Δ] is approximately λ (t | H (t))Δ . This step of point process modeling
and filtering paradigm is similar to proposed stochastic description of microseismic distribution

93
to account for coupled physics and seismic modeling uncertainty. The next step of point process
modeling and filtering paradigm is derivation of recursive Bayesian inversion framework which
involve a sequence of forecast and update steps, similar to the Kalman filter, with the forecast
equation given by the coupled flow and geomechanics simulation to compute the predicted
states. Then, the update equation of the point process filter is used to obtain updated state
variables.

94
CHAPTER 5
Summary and Conclusion
Characterization of dynamic reservoirs, in particular, fractured reservoirs in unconventional oil
and gas development can significantly improve production efficiency and reduce the related
environmental impacts. In this work, we use a coupled flow and geomechanical model with the
Ensemble-based data assimilation workflow to infer matrix and fracture properties and
heterogeneity within rock facies from tracer test, production and microsesimic data. In the first
series of numerical experiments, we considered planar fracture models and described them
deterministically. We first perform a sensitivity analysis to quantify the effect of fracture half-
length, conductivity, and conductivity decline rate, along with matrix permeability and porosity,
on chemical tracer transport and production response of the reservoir. The data assimilation was
performed sequentially, first using the tracer concentration data and a transport forecast model to
update fracture and matrix properties. The outputs from tracer data assimilation were then used
in the second experiment to assimilate production data and further constrain these properties. In
the third experiment, we used a discrete fracture network model and included additional
properties related to the stochastic DFN description, including natural fracture density,
conductivity, and conductivity decline rate.
The results from our sensitivity analysis and data assimilation experiments show that tracer and
production data have complementary information. The tracer data follow a transport model over
a short period of time and tends to show sensitivity to parameters that are involved in the short-
time response of the model (mainly fracture conductivity and natural fracture density). For
instance, matrix properties show little sensitivity to tracer response due to lack of interaction
between the tracer and the matrix. On the other hand, the production data follows a fluid flow
model over a long period, which tends to show strong sensitivity to fracture half-length and

95
natural fracture density (which control the SRV) and matrix permeability and porosity (which
control per-volume production). When tracer and production data are both integrated, the
parameters of the model are all updated and the uncertainty in future predictions is significantly
reduced. We incorporated the impact of dynamic changes in bot hydraulic and natural fracture
conductivity due to fracture closing during production, using a simple model to relate the
changes in stress (pressure) to fracture conductivity. The model provided a fracture conductivity
decline rate as a new parameter that was successfully updated by production data.
For the DFN model, we considered a stochastic description of the fractures and used the
production data to update the fracture density in space without estimating the exact locations of
the fractures. This assumption was adopted to make sure that the data contains the information to
update the resulting parameter (fracture density and conductivity). We note that since our
parameterization of the natural fractures did not include spatial resolution (other than allowing
different density for each fracturing stage), the same fracture density could lead to very different
natural fracture distributions. However, the cumulative production response of such models
would be more sensitive to fracture density and connectivity (i.e., SRV) than it would be to
fracture distribution in space.
The results from our study suggest that the EnKF can be used to successfully characterize
fracture and matrix properties from dynamic tracer test and production data. However, typically
matrix properties are not homogeneous and there is variability in rock properties in various
scales. Therefore, in the second part of this study, we limit our study to matrix without fractures
and we propose a distance transform method to characterize the rock facies distribution.

96
Estimating discrete facies models from dynamic flow and transport data is complicated by the
categorical (discontinuous) nature of the parameters (facies values). Existing model calibration
techniques are designed primarily for estimating continuous model parameters and, hence, are
not suitable for inferring discrete model parameters from flow data. We develop a distance
transform method that converts discrete facies models into continuous parameters, namely
distances from facies boundaries, to facilitate model calibration. The resulting continuous
distance maps tend to inherit similar connectivity patterns as the facies models. To better
preserve the implied connectivity patterns in the prior models of facies and to improve the
estimation of connected distance maps, we propose the use of singular value decomposition
(SVD) to parametrize the distance maps before model calibration. The SVD parameterization is
an optional step that is useful when large-scale facies are considered. When this parameterization
is used for model calibration, the ensemble smoother (ES) with multiple data assimilation
method is used to update the SVD weights of the distance maps. The updated SVD coefficients
are then back transformed to obtain the updated distances, followed by an inverse distance
transform to convert the updated distance maps to calibrated discrete facies models. Application
of the distance transform without SVD parameterization only involves data assimilation to
directly update the initial ensemble of distance maps, followed by an inverse distance
transformation to recover the updated facies. The results from our experiments with two and
three facies models show that the distance-transform method can be applied to ensure that the
estimated facies maps during model calibration remain discrete. In addition, the introduction of
the SVD parameterization of distance maps is justified when the facies exhibit large-scale
connectivity. For improved implementation of the SVD-based distance transform, the prior
ensemble of facies models is used to develop a relation between inverse distance transform

97
threshold and the SVD parameterization level. We also studied a case in which the proposed
method was used to successfully estimate the facies distribution and their associated permeability
values simultaneously.
By applying the discrete distance transform to convert categorical facies maps to continuous
distance functions, the presented work enables the calibration of discrete geologic facies models
through standard model calibration techniques for updating continuous rock properties. The use
of ensemble-based data assimilation is motivated by the promising performance of these
techniques as reported in the recent literature. The results show that combining distance
transforms with SVD parameterization in the distance domain can lead to better preservation of
facies connectivity. However, care must be exercised in using the SVD parameterization when
the expected facies patterns have short correlation lengths, which may not be properly
represented or updated using SVD parameterization. Lastly, the SVD method is not the only
parameterization that can be applied to improve facies connectivity. Other parameterization
methods can also enhance solution continuity and lead to improvements in capturing the
expected geologic connectivity patterns. An alternative approach to calibration of discrete
geologic facies is learning the salient geologic features from prior facies realizations and using
the learned features in reconstructing solutions to model calibration problems.
While the forward model of our data assimilation study for these experiments was based on fluid
flow and transport equations, a more in-depth analysis would require the use of coupled physics
models (e.g., flow and geomechanical simulations) to better describe fracture and rock
properties. There are also other types of data such as microseimic data which can provide useful
insight about the location of failure and geomechanical properties of the reservoir. To have more
complex and representative model for subsurface reservoir, we built a coupled flow and

98
gemechanical simulation model using CMG-GEM module. The purpose is to have a more
realistic model and to estimate hydraulic and geomechanical properties of the rock such as
Young’s modulus, Poisson’s ratio, friction angle, cohesion and permeability by integrating
discrete microseismic data. Therefore, we have combined the complex coupled flow and
geomechanical simulation model with Mohr-Coulomb rock mechanical failure to predict
potential rock failure map and eventually microseismic cloud. Since Microseismic data are
discrete data and cannot be used as an input to conventional data assimilation techniques such as
EnKF, we propose a stochastic modeling and filtering paradigm to be adaptive to discrete data.
In this study, the main aspects of proposed workflow for characterization of fractured reservoirs
can be summarized as: 1). Integrating different types of data such as tracer, production and
microseismic data, 2). Having more complex and representative coupled flow ad geomechanical
model for subsurface reservoir, 3). Addressing some of the existing issues associated with
conventional data assimilation techniques, 4). Extending the application of Ensemble Kalman
Filter for characterization of discrete parameter distribution (such as rock facies map) and
introducing techniques for integration of discrete microseismic data.

99
Appendix A: Nomenclature
S: host phase saturation
C: flowing tracer concentration
𝐶 𝑎 : adsorbed tracer concentration
𝜌 𝑟 : mass density of the rock formation
ϕ: porosity
ρ: host phase density (here is water)
𝜇 : host phase viscosity (water)
𝐷 𝑧 :center depth
B: host phase formation volume factor
T: transmissibility
𝑘 𝑟 : host phase relative permeability
V: block pore volume
Q: host phase production rate
P: pressure
g: gravity acceleration
K: absolute permeability
ơ: effective poroeastic stress
𝐵 𝑝 𝑎𝑛𝑑 𝐵 ơ
: permeability modulus
𝛼 𝑐 : Biot’s constant
ν: Poisson’s ratio
𝐷 𝑐 : the tracer diffusion coefficient
DF: diffusivity
λ: the tracer decay constant
𝐺 : Kalman gain
𝐻 : measurement matrix
𝐶 𝑥 𝑒 : state sample covariance matrix
𝐶 𝑦 : observation covariance matrix
𝑗 : realization index (superscript)
𝑡 : time index (subscript)
x: state vector
𝑧 : augmented state vector
𝑦 𝑜𝑏𝑠 : perturbed observation vector
𝑚 : fracture parameter vector
𝑓 : nonlinear state propagation operator
ℎ: measurement operator
𝘀 : uncorrelated Gaussian observation error
𝛼 : damping factor

100
Appendix B: Flow Regimes in Production from Fractured Models
In this appendix we investigate the validity of the forward modeling approach used in this paper.
This is done by examining the adjustment to fracture properties when large fracture apertures are
used and by illustrating the resulting production behavior. Finally, we compare the forward
simulation result of our approach with a more direct approach to simulating discrete fracture
models, which is implemented in the MRST simulation package. To reduce the error associated
with representing fractures as porous matrix with large aperture sizes, the cells that contain
fractures have to be treated with care. This is typically done by assigning small grid sizes
through local grid refinement around the fractures and adjusting cell properties to ensure that the
effective fracture properties (e.g., conductivity and porosity) are preserved. A description of
these changes in our modeling approach is presented below.

B1) Adjustment to Fracture Properties
Transmissibility within fractures
Considering 1D flow in a horizontal fracture across neighboring fracture cells, shown in Figure
B-1a, and the fracture permeability and conductivity 𝑘 𝑓 and 𝐶 𝑓 , respectively, the transmissibility
is expressed as
𝑇 𝑦 = 0.001127 𝑘 𝑓 . 𝐴 /𝐿 = .001127 𝑘 𝑓 . 𝑑 𝑥 . 𝑑𝑧 /𝐿 = .001127 𝐶 𝑓 . 𝑑𝑧 /𝐿 (B1)
where 𝑘 𝑓 𝑊 = 𝐶 𝑓 is the fracture conductivity, which stays constant after adjusting the
permeability corresponding to the fracture aperture (grid size). Therefore, for a fixed width, the

101
transmissibility does not change along the fracture length after adjusting the fracture
permeability. The flow rate of phase i across neighboring blocks is expressed as:
q
i
=
.001127k
f
A
L
Kr
i
µ
i
(P1-P2) =
.001127C
f
dz
L
Kr
i
µ
i
(P1-P2) (B2)
where with 𝐴 = 𝑊𝑑𝑧 is used. The resulting equation shows that the adjusted fracture
permeability does not affect the flowrate and productivity index.

Transmissibility across fracture and matrix
In the case of transmissibility across fracture and matrix (Figure B-1b), we have
𝑇 𝑥 =
0.001127 𝐴 𝐵 =
0.001127 𝐴 𝑑𝑥 𝑓 𝑘 𝑓 +
𝑑𝑥 𝑚 𝑘 𝑚 2

Assuming 𝑑𝑥 𝑓 1
= 0.01 and 𝑘 𝑓 1
=120 D, and the adjusted 𝑑𝑥 𝑓 2
= 3 ft and 𝑘 𝑓 2
= 0.4 D, we get
dx
f1
k
f1
=
0.01
120
= 8.3e
-5
(𝑏𝑒𝑓𝑜𝑟𝑒 𝑎𝑑𝑗𝑢𝑠𝑖𝑛𝑔 𝑝𝑒𝑟𝑚𝑒𝑎𝑏𝑖𝑙𝑖𝑡𝑦 )
dx
f2
k
f2
=
3
0.4
= 7.5 (𝑎𝑓𝑡𝑒𝑟 𝑎𝑑𝑗𝑢𝑠𝑖𝑛𝑔 𝑝𝑒𝑟𝑚𝑒𝑎𝑏𝑖𝑙𝑖𝑡𝑦 )
dx
m
k
m
=
30
1e-8
= 3.0e
9
(𝑚𝑎𝑡𝑟𝑖𝑥 )
By adjusting the fracture permeability, the ratio
𝑑𝑥
𝑓 𝑘 𝑓 may be overestimated, but when the
transmissibility between fracture and matrix is calculated, this ratio is negligible compared to the

102
same ratio for the matrix
𝑑𝑥
𝑚 𝑘 𝑚 . Hence, for this example, the error in estimation of transmissibility
in x direction is negligible:
𝑇 𝑥 =
0.001127 𝐴 8.3𝑒 −5
+ 3.0𝑒 9
2
= 7.5𝑒 −13
𝐴 (𝑏𝑒𝑓𝑜𝑟𝑒 𝑎𝑑𝑗𝑢𝑠𝑖𝑛𝑔 𝑝𝑒𝑟𝑚𝑒𝑎𝑏𝑖𝑙𝑖𝑡𝑦 )
𝑇 𝑥 =
0.001127 𝐴 7.5 + 3.0𝑒 9
2
= 7.5𝑒 −13
𝐴 (𝑎𝑓𝑡𝑒𝑟 𝑎𝑑𝑗𝑢𝑠𝑖𝑛𝑔 𝑝𝑒𝑟𝑚𝑒𝑎𝑏𝑖𝑙𝑖𝑡𝑦 )
Figure B2 shows the expected flow regimes over a wide range of time-scales. These flow
regimes include pseudo linear flow (slope =0.5), pseudo steady sate flow (slope =1), compound
linear flow (slope =0.5), pseudo radial flow (slope =0), and the boundary behavior (slope =1).
Porosity:
The fracture porosity can affect the amount of hydrocarbon stored in the fracture and the
computed pressures as follows (assuming a single phase flow). The fracture pore volume (which
is modeled as a porous medium) can be calculated as Vp = length × aperture × height × porosity.
In order to preserve the pore volume, fracture porosity also should be adjusted based on the
assumed fracture aperture.

B2) Validity of Forward Modeling Approach
Several authors have investigated the accuracy of finite difference models with adjusted fracture
properties for production prediction from unconventional reservoirs (Cipolla et al, 2010; Le et
al., 2011). In this part, we also show the sensitivity of production to fracture cell sizes using
different grid refinements. Figure B3 shows the sensitivity analysis by sequentially refining the
fracture cell width (and adjusting the permeability and porosity accordingly) and comparing the

103
simulation outputs for each case. The results clearly show that, if the fracture properties are
adjusted properly, increasing the fracture aperture from 0.1ft to 3ft does not result in a noticeable
change in the simulation output. We also compared the results from our fracture modeling
approach with a direct discrete fracture modeling framework that is implemented in MRST
(MRST-2015b). Figure B-4 shows that the pressure map and production profile generated by
DFM model is in good agreement with those from grid refinement approach used in our study.
Finally, we note that for data assimilation with DFN models, the uncertainty in the distribution
and connectivity of the fractures tends to dominate the prediction errors. Therefore, the
approximations involved in the simulation tend to be overshadowed by the dominating effect of
uncertainty in natural fracture distribution.

104
REFERENCES
Aanonsen, S., Naevdal, G. Oliver, D. and et al. 2009.The ensemble Kalman filter in reservoir
engineering- a review, SPE J., 14(3):393–412.
Alcolea, A., Hidalgo, J. J., Carrera, J., Alcolea, A., & Medina, A. 2005. Inverse problem in
hydrogeology Inverse problem in hydrogeology, http://doi.org/10.1007/s10040-004-0404-7
Anderson J.L. 2007. Exploring the need for localization in ensemble data assimilation using a
hierarchical ensemble filter, Physica D: Nonlinear Phenomena. 230(1–2):99-111. ISSN
0167-2789, http://dx.doi.org/10.1016/j.physd.2006.02.011.
Bachman, R. C., Sen, V., Khalmanova, D., et al. 2011. Examining the Effects of Stress
Dependent Reservoir Permeability on Stimulated Horizontal Montney Gas Wells. Society
of Petroleum Engineers. doi:10.2118/149331-MS
Barree, R. D. and Kukherjee, H. 1995. Engineering Criteria for Fracture Flowback Procedures.
SPE Rocky Mountain Regional/Low Permeability Reservoir Symposium, Denver, CO,
March 20-22.Cooke, C. E., Jr.: “Effect of Fracturing Fluids on Conductivity,” Transactions
of AIME. SPE 29600. http://dx.doi.org/10.2118/29600-MS
Bertino L, G. Evensen, H. Wackernagel (2003). Sequential data assimilation techniques in
oceanography. Int Stat Rev. 71(2):223–41. http://dx.doi.org/10.1111/j.1751-
5823.2003.tb00194.x.
Bhark, E. W., B. Jafarpour, and A. Datta-Gupta (2011). A generalized grid connectivity-based
parameterization for subsurface flow model calibration. Water Resources Research, 47(6).
http://doi.org/10.1029/2010WR009982
Borgefors G. (1986). Distance transformations in digital images. Computer Vision, Graphics and
Image Processing, 34(3):344–371

105
Boyer, C., Kieschnick, J., and Suarez-Rivera, L. R., Waters, G. 2006. Producing Gas from its
Source, Oilfield review, 36-49.
Bradford, S. A., Yates, S. R., Bettahar, M., & Simunek, J. (2002). Physical factors affecting the
transport and fate of colloids in saturated porous media, 38(12), 1–12.
http://doi.org/10.1029/2002WR001340
Cacas, M. C., Ledoux, E., De Marsily, G. and et al. 1990. Modeling Fracture Flow with A
Stochastic Discrete Fracture Network: Calibration and Validation: 1. The Flow Model.
Water Resources Research. 26(3):479-489.
Caers J., Zhang T. (2004). Multiple-point geostatistics: a quantitative vehicle for integrating
geological analogs into multiple reservoir models. AAPG Mem 80:383-394.
Caers J. (2003). Efficient gradual deformation using a streamline-based proxy method. J Pet Sci
Eng, 0920-4105 2003;39(1–2):57–83. http://dx.doi.org/10.1016/S0920-4105(03)00040-8.
Cardiff M., and P.K. Kitanidis (2009), Bayesian inversion for facies detection: An extensible
level set framework, Water Resources Research, Vol. 45(10), W10416, DOI:
10.1029/2008WR007675.
Carle, S.F., and Fogg, G.E. (1997). Modeling spatial variability with one and multi-dimensional
continuous Markov chains. Mathematical Geology 29(7), 891–918.
Carle, S.F., LaBolle, E.M., Weissman, G.S., van Brocklin, D., and Fogg, G.E. (1998).
Conditional simulation of hydrofacies architecture: A transition probability/Markov chain
approach. In Hydrogeologic Models of Sedimentary Aquifers, G.S. Fraser, and J.M. Davis,
eds., Concepts in Hydrogeology No. 1, SEPM (Society for Sedimentary Geology) Special
Publication, pp. 147–170.

106
Carrera, J., A. Alcolea, A. Medina, J. Hidalgo, and L. J. Slooten (2005), Inverse problem in
hydrogeology, Hydrogeol. J., 13(1), 206–222
Carrera, J. and S.P. Neuman (1986a) Estimation of aquifer parameters under steady-state and
transient condition: 2. Uniqueness, stability, and solution algorithms, Water Resour. Res.,
22(2), 211-227.
Carrera J, S.P. Neuman (1986b). Estimation of aquifer parameters under transient and steady-
state conditions. 1. Maximum likelihood method incorporating prior information. Water
Resour Res 22(2):199–210.
Chang, H., D. Zhang, Z. Lu (2010). History matching of facies distribution with the EnKF and
level set parameterization, Journal of Computational Physics, Volume 229, Issue 20, 1.
Pages 8011-8030, ISSN 0021-9991, http://dx.doi.org/10.1016/j.jcp.2010.07.005.
Chen, Y., and Zhang, D. 2006. Data assimilation for transient flow in geologic formations via
ensemble Kalman filter. Adv. Water Resour. 29: 1107–1122.
Chen, Y., and D. Zhang (2006), Data assimilation for transient flow in geologic formations via
ensemble Kalman filter, Adv. Water Resour., 29, 1107–1122.
Cipolla, C.L., Wemg, X., Mack, M.G. et al. 2011. Integrating Microseismic Mapping a Complex
Fracture Modeling to Characterize hydraulic Fracture Complexity. SPE Hydraulic
Fracturing Technology Conference, The Wood Land, Texas, USA, 24-26 January. SPE
140185.
Cipolla, C. L., Williams, M. J., Weng, X., et al. 2010. Hydraulic Fracture Monitoring to
Reservoir Simulation: Maximizing Value. SPE. doi:10.2118/133877-MS.
Cipolla, C.L., Lolon E.P., Erdle J.C., Rubin B. (2010), Reservoir modeling in shale-gas
reservoirs. SPE Reservoir Evaluation & Engineering Journal. SPE-125530-PA.

107
Cipolla, C. L. and Wright, C. A. 2000. State-of-the-Art in Hydraulic Fracture Diagnostics. SPE.
doi:10.2118/64434-MS.
Clarkson, C. R., Qanbari, F., Nobakht, M. and et al. 2012. Incorporating Geomechanical and
Dynamic Hydraulic Fracture Property Changes into Rate-Transient Analysis: Example
from the Haynesville Shale. SPE. doi:10.2118/162526-MS
Cox, S. A., Cook,D. M. and et al.(2008). Unconventional Resource Play Evaluation: A Look at
the Bakken Shale Play of North Dakota. SPE Unconventional Reservoirs Conference, 10-
12 February, Keystone, Colorado, USA. SPE-114171-MS.
http://dx.doi.org/10.2118/114171-MS
Crain E. R. 2010. Crain’s Petrophysical Handbook, Rocky Mountain House.
Crestani, E. 2013. Tracer Test Data Assimilation for the Assessment of Local Hydraulic
Properties in Heterogeneous Aquifers. Doctoral thesis. Universit_a degli Studi di Padova
Crockett, A. R., Willis, R. M., and Cleary, M. P. 1989. Improvement of Hydraulic Fracture
Predictions by Real-Time History Matching on Observed Pressures. SPE.
doi:10.2118/15264-PA.
Dai, Z. and Samper, J. (2004), Inverse problem of multicomponent reactive chemical transport in
porous media: Formulation and Applications, Water Resour. Res., Vol 40, W07407,
doi:10.1029/2004WR003248.
Daley, D.J, Vere-Jones, D. (1988). An Introduction to the Theory of Point Processes. Springer,
New York. ISBN 0-387-96666-8, MR 950166.
De Marsily, G., Delay, F., Goncalves, J., Renard, P.H., Teles, V., and Violette, S. (2005).
Dealing with spatial heterogeneity. Hydrogeology Journal 13(1), 161–183.

108
Descant, F. J., Blackwell, R., Pope, GA. and et al. 1989. The Use of Single Well Tracer Testing
to Estimate Heterogeneity. Richardson, TX (1989). SPE 20303.
Deutsch C.V., Journel A.G., (1998). GSLIB - Geostatistical Software Library and User’s Guide.
Oxford University Press, New York - Oxford.
Deutsch, C.V., and L.Wang (1996). Hierarchical object-based geostatistical modeling of fluvial
reservoirs, paper SPE 36514 presented at the 1996 SPE Annual Technical Conference and
Exhibition, Denver, Oct 6–9.
Diaz De Souza, O. C., Sharp, A., Martinez, R. C. and et al. I. 2012. Integrated Unconventional
Shale Gas Reservoir Modeling: A Worked Example From the Haynesville Shale, De Soto
Parish, North Louisiana. 2012, January 1. Society of Petroleum Engineers.
doi:10.2118/154692-MS
Doherty, J., and B. E. Skahill (2006). An advanced regularization methodology for use in
watershed model calibration, J. Hydrol., 327, 564–577, doi:10.1016/j.jhydrol.2005.11.058.
Donoho, D. L. (2006), Compressed sensing, IEEE Trans. Inf. Theory, 52(4), 1289–1306
Dovera, Laura, D.R. Ernesto (2010), Multimodal ensemble Kalman filtering using Gaussian
mixture models, Computational Geosciences. 15: 307. doi:10.1007/s10596-010-9205-3
Drake, S. 2007. Unconventional Gas Plays, AAPL, Southwest Land Institute.
Eclipse (2014), Eclipse reference manual and technical description, Schlumberger-SEES,
University of Southern California.
Eisner, L., Hulsey, B. J., Duncan, P. and et al . 2010. Comparison of surface and borehole
locations of induced seismicity. Geophysical Prospecting. 58(5):809–820.
http://doi.org/10.1111/j.1365-2478.2010.00867.x

109
Eisner, L., Thornton, M., and Griffin, J. 2011. Challenges for microseismic monitoring. SEG
Annual Meeting, San Antonio, Texas, 1519–1523. http://doi.org/10.1190/1.3627491
Elad, M., (2010), Sparse and Redundant Representations: From Theory to Application in Signal
and Image Processing, Springer.
Elahi, S. H., and Jafarpour, B. 2015. Characterization of Fracture Length and Conductivity From
Tracer Test and Production Data With Ensemble Kalman Filter. Unconventional Resources
Technology Conference. doi:10.15530/URTEC-2015-2174307
Emerick, A. A. and Reynolds A. C. 2013. Ensemble smoother with multiple data assimilations,
Computers& Geosciences, 55, 3–15. http://dx.doi.org/10.1016/j.cageo.2012.03.011
Evensen, G. 2009. The ensemble Kalman filter for combined state and parameter estimation:
Monte-Carlo techniques for data assimilation in large systems. IEEE Control Syst. Mag.,
29(3) 83–104. doi:10.1109/ MCS.2009.932223.
Evensen, G. 2007. Data Assimilation: The Ensemble Kalman Filter, Springer, New York.
Evensen, G., and P. J. van Leeuwen (1996), Assimilation of Geosat altimeter data for the
Agulhas current using the ensemble Kalman filter with a quasigeostrophic model, Mon.
Weather Rev., 124, 85–96.
Evensen, G., and van Leeuwen P. J. 1996. Assimilation of Geosat altimeter data for the Agulhas
current using the ensemble Kalman filter with a quasigeostrophic model, Mon. Weather
Rev., 124: 85–96.
Eclipse. 2011. Eclipse reference manual and technical description, Schlumberger-SEES,
University of Southern California.
Felzenszwalb, P. F. Huttenlocher, D.P (2004), Distance Transforms of Sampled Functions.

110
Ferreira, L. E. A., Descant, F., Delshad, M., and et al. 1992. A Single-Well Tracer Test To
Estimate Wettability. SPE/DOE Eighth Symposium on Enhanced Oil Recovery, Tulsa, 22-
24 April 1992. SPE/DOE 24136
Franquet, M., Ibrahim, M., Wattenbarger, R. A., and et al. 2004. Effect of Pressure-Dependent
Permeability in Tight Gas Reservoirs, Transient Radial Flow. Petroleum Society of
Canada. 1 January 2004. doi:10.2118/2004-089
Gardien, C. J., Pope, G. A., & Hill, A. D. 1996. Hydraulic Fracture Diagnosis Using Chemical
Tracers. Society of Petroleum Engineers. 1 January 1996. doi:10.2118/36675-MS
Gavalas G.R., Shah P.C., Seinfeld J.H. (1976).Reservoir history matching by Bayesian
estimation. Soc. Petrol. Eng. J., 16(6), 337–350.
Gelb, A. 1974. Applied Optimal Estimation, MIT Press, Cambridge, Mass.
Ghods, P., & Zhang, D. 2012. Automatic Estimation of Fracture Properties in Multistage
Fractured Shale Gas Horizontal Wells for Reservoir Modeling. Society of Petroleum
Engineers. 1 January 2012. doi:10.2118/153913-MS
Grimstad A. A., Mannseth T., Nævdal G., Urkedal H. (2003). Adaptive multiscale permeability
estimation, Computational Geosciences 7, 1-25.
Grob M. , Mirko van der Baan (2011). ”Inferring in-situ stress changes by statistical analysis of
microseismic event characteristics.” The Leading Edge, 30(11), 1296-1301.doi:
10.1190/1.3663403
Guardiano, F., and R. Srivastava (1993). Multivariate geostatistics: Beyond bivariate moments,
Geostatistics-Troia, Kluwer Academic Publications, Dordrecht, 1,133–144.

111
Harp D., Z. Dai, A. Wolfsberg, J. Vrugt, B. Robinson, V. Vesselinov (2008), Aquifer structure
identification using stochastic inversion, Geophysical Research Letters, 35, L08404,
doi:10.1029/2008GL033585.
He J, P. Sarma , L. Durlofsky (2012). Reduced-order flow modeling and geological
parameterization for ensemble-based data assimilation. Comput Geosci. 55:54–69.
http://dx.doi.org/10.1016/j.cageo.2012.03.027.
Hill MC, Tiedeman CR. (2006). Effective groundwater model calibration: with analysis of data,
sensitivities, predictions, and uncertainty. John Wiley & Sons.

Hoffman BT, J. Caers (2005). Regional probability perturbations for history matching. J Pet Sci
Eng, 0920-4105 2005;46(1–2):53http://dx.doi.org/10.1016/j.petrol.2004.11.001.–71.
Holditch, S. A., Holcomb, D. L., & Rahim, Z. 1993. Using Tracers To Evaluate Propped
Fracture Width. Society of Petroleum Engineers. 1 January 1993. doi:10.2118/26922-MS
Holzer, S., S. Hinterstoisser, S. Ilic, and N. Navab (2009). Distance transform templates for
object detection and pose estimation. In Computer Vision and Pattern Recognition, 2009.
CVPR 2009. IEEE Conference on (pp. 1177-1184). IEEE
Houtekamer PL, HL. Mitchell (2001). A sequential ensemble Kalman filter for atmospheric data
assimilation. Mon Weather Rev, 129(1):123–37. http://dx.doi.org/10.1175/1520-
0493(2001)129<0123:ASEKFF>2.0.CO;2.
Hu LY. (2008). Extended probability perturbation method for calibrating stochastic reservoir
models. Math Geosci, 1874-8961. 40(8):875–85. http:// dx.doi.org/10.1007/s11004-008-
9158

112
Hunt RJ, J. Doherty, MJ. Tonkin (2007). Are models too simple? Arguments for increased
parameterization. Ground Water. 45(3):254–62. http:// dx.doi.org/10.1111/j.1745-
6584.2007.00316.x.
Huseby,O.K., Valestrand, R., Naevdal, G., and et al. 2011. Natural and Conventional Tracers for
Improving Reservoir Models Using the EnKF Approach. EUROPEC/EAGE Conference
and Exhibition, 8-11 June, Amsterdam, The Netherlands. SPE-121190-MS.
http://dx.doi.org/10.2118/121190-MS
Huttenlocher D, G. Klanderman, and W. Rucklidge. (1993), Comparing images using the
hausdorff distance. IEEE Transactions on Pattern Analysis and Machine Intelligence,
15(9):850–863.
Jacquard P., Jain C. (1965). Permeability distribution from field pressure data. Soc. Pet. Eng.
Journal, 281-294.
Jafarpour, B. (2011), Wavelet reconstruction of geologic facies from non- linear dynamic flow
measurements, IEEE Trans. Geosci. Remote Sens., 49, 1520–1535,
doi:10.1109/TGRS.2010.2089464.
Jafarpour, B., and D. McLaughlin (2009), Estimating channelized reservoir permeabilities with
the ensemble Kalman filter: The importance of the ensemble design, SPE J., 14(2), 374–
388.
Jafarpour, B., and D. McLaughlin (2009), History matching with an ensemble Kalman filter and
discrete cosine parameterization, Computational Geosciences. 12(2), 1573-1499.
Jafarpour B., Goyal V.K., McLaughlin D.B., Freeman W.T (2010): Compressed history
matching: Exploiting transform-domain sparsity for regularization of nonlinear dynamic
data integration problems, Mathematical Geosciences, 42(1), PP 1-27.

113
Jafarpour, B., and M. Tarrahi (2011), Assessing the performance of the ensemble Kalman filter
for subsurface flow data integration under variogram uncertainty, Water Resour. Res., 47,
W05537, doi:10.1029/2010WR 009090.
Jafarpour, B., M. Khodabakhshi (2011), A Probability Conditioning Method (PCM) for
Nonlinear Flow Data Integration into Multipoint Statistical Facies Simulation. Math.
Geosci. 43, 133–164 .
Jahns HO. (1966). A rapid method for obtaining a two-dimensional reservoir description from
well pressure response data, Soc. Pet. Eng. J. 6(4), 315–327.
Jolliffe, I. (2005), Principal Component Analysis, Springer.

Kalman, R. E. 1960. A new approach to linear filtering and prediction problems, J. Basic Eng.
82: 35–45.
King, G. E. 2010. Thirty Years of Gas Shale Fracturing: What Have We Learned? SPE.
doi:10.2118/133456-MS.
Khaninezhad, M. M., and B. Jafarpour (2014). Sparse Randomized Maximum Likelihood
(SpRML) for subsurface flow model calibration and uncertainty quantification. Advances
in Water Resources, 69, 23–37. http://doi.org/10.1016/j.advwatres.2014.02.005
Khaninezhad, M. M., B. Jafarpour and L. Li (2012a). Sparse geologic dictionaries for subsurface
flow model calibration: Part I. Inversion formulation. Advances in Water Resources, 39,
106–121.
Khaninezhad, M. M., B. Jafarpour and L. Li (2012b). Sparse geologic dictionaries for subsurface
flow model calibration: Part II. Robustness to uncertainty. Advances in Water Resources,
39, 122–136.

114
Kovacevic, J., V. K. Goyal, and M. Vetterli (2013), Fourier and Wavelet Signal Processing,
Cambridge Univ. Press.
Lange, A., Bourbiaux, B. J., & Bouzian, J. 2005. Tracer-Test Simulation on Discrete Fracture
Network Models for the Characterization of Fractured Reservoirs. Society of Petroleum
Engineers. 1 January 2005. doi:10.2118/94344-MS

Lawson, C. L., and J. R. Hanson (1995), Solving Least Squares Problems, Classics Appl. Math.,
vol. 15, 337 pp., Soc. for Ind. and Appl. Math., Philadelphia, Pa.
Lee, J., and P.K. Kitanidis (2013), Bayesian inversion with total variation prior for discrete
geologic structure identification, Vol 49(11), PP 7658–7669.
Li J, Du C,M, Zhang Z. 2011. Critical evaluation of shale gas reservoir simulation approaches:
single-porosity and dual-porosity modeling. SPE Middle East Unconventional Gas
Conference and Exhibition, Muscat, Oman, 31 January - 2 February2011.
Li L., Jafarpour B. (2010), Effective solution of nonlinear subsurface flow inverse problems in
sparse bases, Inverse Problems, 26(10), 105-016 (24pp).
Liu, N., D.S. Oliver (2005). Ensemble Kalman filter for automatic history matching of geologic
facies. J. Pet. Sci. Eng. 47(3–4), 147–161
Lolon E., Cipolla, C.L. and Weijers, L. and et al. (2009). Evaluating Horizontal Well Placement
and Hydraulic Fracture Spacing/Conductivity in the Bakken Formation, North Dakota. SPE
Annual Technical Conference and Exhibition, 4-7 October, New Orleans, Louisiana. SPE-
124905-MS. http://dx.doi.org/10.2118/124905-MS

115
Lorentzen, R. J., G. Nævdal and A. Shafieirad (2012), Estimating Facies Fields by Use of the
Ensemble Kalman Filter and Distance Functions--Applied to Shallow-Marine
Environments. Society of Petroleum Engineers. doi:10.2118/143031-PA.
Lu, P. B., and R. N. Horne (2000), A multiresolution approach to reservoir parameter estimation
using wavelet analysis, paper presented at SPE Annual Technical Conference and
Exhibition, Soc. of Pet. Eng., Dallas, Tex., 1–4 Oct.
Lu P, Horne RN (2000). A multiresolution approach to reservoir parameter estimation using
wavelet analysis, paper presented SPE annual technical conference, Soc. of Pet. Eng.
Ma W, Jafarpour B (2017). Conditioning Multiple-Point Geostatistical Facies Simulation on
Nonlinear Flow Data Using Pilot Points Method, Paper SPE-185629-MS, presented at the
SPE Western Regional Meeting held in Bakersfield, CA, USA. 23-27 April, 2017.
Mallat, S. (2008), A Wavelet Tour of Signal Processing: The Sparse Way, 3rd ed., Academic
Press.
McDaniel, R. R., Holmes, D. V., Borges, J., Bajoie, B. J., Peeples, C., & Gardner, R. 2009.
Determining Propped Fracture Width from a New Tracer Technology. Society of
Petroleum Engineers. 1 January 2009. doi:10.2118/119545-MS

Moradkhani H, S. Sorooshian, HV. Gupta H, PR. Houser (2005). Dual state-
parametervestimation of hydrological models using ensemble Kalman filter. Adv Water
Resour. 28(2):135–47. http://dx.doi.org/10.1016/ j.advwatres.2004.09.002.
Marquart, G., Vogt, C., Klein, C., and et al. 2013. Estimation of Geothermal Reservoir Properties
Using the Ensemble Kalman Filter. Energy Procedia. 40:117-126.
http://dx.doi.org/10.1016/j.egypro.2013.08.015

116
Moreno, J., Tarrahi, M., Gildin, E., & Gonzales, S. 2014. Real-Time Estimation of Hydraulic
Fracture Characteristics From Production Data. Unconventional Resources Technology
Conference. 25 August. doi:10.15530/URTEC-2014-1923687
Moreno, D.L. and S.I. Aanonsen (2011). Continuous Facies Updating Using the Ensemble
Kalman Filter and the Level Set Method. Mathematical Geosci. 43 (8), 951–970.
http://dx.doi.org/10.1007/s11004-011-9347-4.
Naevdal, G., Johnsen, L., Aanonsen, S., and et al. 2005. Reservoir monitoring and continuous
model updating using ensemble Kalman filter. SPE J. 10(1): 66–74.
Nelson, R. 2001. Geologic Analysis of Naturally Fractured Reservoirs. Gulf Professional
Publishing. Butterworth–Heinemann
Ogata, Y. 1999. Seismicity analysis through Point-process modeling: A Review. Pure and
Applied Geophysics.
Okouma Mangha, V., Guillot, F., Sarfare, M., et al. 2011. Estimated Ultimate Recovery (EUR)
as a Function of Production Practices in the Haynesville Shale. Society of Petroleum
Engineers. 1 January 2011. doi:10.2118/147623-MS
Oliver D.S., Chen Y.(2011).Recent progress on reservoir history matching: a review.
Computational Geosciences, Volume 15, Number 1. 185-221, DOI: 10.1007/s10596-010-
9194-2
Oliver, D. S., A. C. Reynolds, and N. Liu (2008), Inverse Theory for Petroleum Reservoir
Characterization and History Matching, Cambridge Univ. Press
Parney, R., Cladouhos, T., La Pointe, P., Dershowitz, W., & Curran, B. 2000. Fracture and
Production Data Integration Using Discrete Fracture Network Models for Carbonate

117
Reservoir Management, South Oregon Basin Field, Wyoming. Society of Petroleum
Engineers. 1 January 2000. doi:10.2118/60306-MS
Peyret, O., Drew, J., Mack, M., and et al. 2012. Subsurface To Surface Microseismic Monitoring
for Hydraulic Fracturing. Society of Petroleum Engineers. doi:10.2118/159670-MS
Rao, K. R., and P. Yip (1990), Discrete Cosine Transform: Algorithms, Advantages,
Applications, Academic Press.
Reichle, R., D. McLaughlin, and D. Entekhabi (2002), Hydrologic data assimilation with the
ensemble Kalman filter, Mon. Weather Rev., 130(1), 103–114.
Ritzi, R., Z. Dai, D. Dominic, and Y. Rubin (2004), Spatial correlation of permeability in cross-
stratified sediment with hierarchical architecture, Water Resour. Res., 40,
doi:10.1029/2003WR002420.
Ritzi R., Z. Dai, D. Dominic, Y. Rubin (2006), Reply to comment by Shlomo P. Neuman on
"Spatial correlation of permeability in cross-stratified sediment with hierarchical
architecture", Water Resour. Res., 42, W05602, doi:10.1029/2005WR004402.
Rudin, L. I., S. Osher, and E. Fatemi (1992), Nonlinear total variation based noise removal
algorithms, Physica D, 60(1), 259–268
Sarma P, LJ. Durlofsky, K. Aziz (2008). Kernel principal component analysis for efficient,
differentiable parameterization of multipoint geostatistics. Math Geosci. 40(1):3–32.
http://dx.doi.org/10.1007/s11004-007-9131-7.
Sebacher, B., R. Hanea, A. Heemink (2013). A probabilistic parametrization for geological
uncertainty estimation using the ensemble Kalman filter (EnKF). Computational
Geosciences. doi:10.1007/s10596-013-9357-z

118
Settari, A. T., Bachman, R. C., and Walters, D. A. 2005. How to Approximate Effects of
Geomechanics in Conventional Reservoir Simulation. Society of Petroleum Engineers. 1
January 2005. doi:10.2118/97155-MS
Shah P.C., Gavalas G.R., Seinfeld J.H. (1978). Error Analysis in History Matching: The Optimal
Level of Parametrization”. Soc. Petroleum EnginemJ., June 1978, 219-228.
Silber, R., Martin, J., Willis, S., and et al. 2003. Comparing Fracture Simulation Design to
Radioactive Tracer Field Results: A Case History. Society of Petroleum Engineers. 1
January 2003. doi:10.2118/84842-MS
Strebelle S. (2002). Conditional simulation of complex geological structures using multiple-point
statistics. Mathematical Geology. 34(1), 1-21.

Smith, L., Mase, C. W., & Schwartz, F. W. 1987. Estimation Of Fracture Aperture Using
Hydraulic And Tracer Tests. American Rock Mechanics Association. 1 January 1987.
Tarrahi, Mohammadali., S.H. Elahi, B. Jafarpour (2016). Fast linearized forecasts for subsurface
flow data assimilation with ensemble Kalman filter. Computational Geosciences. 1420-
0597
Tarrahi, M. Jafarpour, B., and Ghasemi, A. 2015. Integration of microseismic monitoring data
into coupled flow and geomechanical models with ensemble Kalman filter. Water
Resources Research. 51(7): 5177–5197. 10.1002/2014WR016264
Tavakoli, R., A.C. Reynolds (2009). History matching with parameterization based on the
singular value decomposition of a dimensionless sensitivity matrix (SPE-118952). In: SPE
Reservoir Simulation Symposium.

119
Thompson, J. M., Nobakht, M., and Anderson, D. M. 2010. Modeling Well Performance Data
from Overpressured Shale Gas Reservoirs. Society of Petroleum Engineers.
doi:10.2118/137755-MS
Tikhonov, A. N., and V. Y. Arsenin (1977), Solutions of Ill-Posed Problems, 330 pp., WH
Winston, Washington, D. C.
Tonkin, M., and J. Doherty (2009). Calibration-constrained Monte Carlo analysis of highly
parameterized models using subspace techniques, Water Resour. Res., 45, W00B10,
doi:10.1029/2007WR006678.
Tonkin, M. J., and J. Doherty (2005). A hybrid regularized inversion methodology for highly
parameterized environmental models. Water Resources Research, 41(10), 1–16.
http://doi.org/10.1029/2005WR003995.
Van Leeuwen, P.J., G. Evensen (1996). Data assimilation and inverse methods in terms of a
probabilistic formulation. Monthly Weather Review 124, 2898–2913.
Veatch Jr., R. W. and Crowell, R. F. 1982. Joint Research/Operations Program Accelerates
Massive Hydraulic Fracturing Technology. Journal of Petroleum Technology. 34(12):
2763-2375. SPE-9337-PA. http://dx.doi.org/10.2118/9337-PA
Weiss, R., and L. Smith (1998). Parameter space methods in joint parameter estimation for
groundwater flow models, Water Resour. Res., 34(4), 647– 661, doi:10.1029/97WR03467.
Woodroof, R. A., Asadi, M., and Warren, M. N. 2003. Monitoring Fracturing Fluid Flowback
and Optimizing Fracturing Fluid Cleanup Using Chemical Frac Tracers. Society of
Petroleum Engineers. doi:10.2118/82221-MS

120
Wright, C. a., Davis, E. J., Weijers, L., and et al. 1997a. Horizontal Hydraulic Fractures: Oddball
Occurrences or Practical Engineering Concern? Proceedings of SPE Western Regional
Meeting. http://doi.org/10.2118/38324-MS
Wright, C. A., Minner, W. A., Weijers, L., and et al. 1997b. Wellbore-to-Fracture
Communication Problems Pose Challenges in California Diatomite Horizontal Wells. SPE
Annual Technical Conference and Exhibition. Society of Petroleum Engineers.
http://doi.org/10.2118/38632-MS
Wright, C., Weijers, L., Davis, E., and et al. 1999. Understanding Hydraulic Fracture Growth :
Tricky but Not Hopeless. SPE. 56724, 10. http://doi.org/10.2523/56724-MS
Xu, W., Thiercelin, M. and Walton, I. 2009. Characterization of Hydraulic Induced Shale
Fracture Network Using a Semi-Analytical Model. SPE Annual Technical Conference and
Exhibition, New Orleans, LA, USA, 4-7 October. SPE 124697
Yilmaz, O., Nur, A., and Nolen-Hoeksema, R. 1991. Pore Pressure Profiles in Fractured and
Compliant Rooks. SPE.
Ye, M., and R. Khaleel (2008). A Markov chain model for characterizing medium heterogeneity
and sediment layering structure, Water Resources Research, 44, W09427,
doi:10.1029/2008WR006924.
Yu, W., Luo Z., and et al. 2014. Sensitivity analysis of hydraulic fracture geometry in shale gas
reservoirs. J Pet Sci and Eng. 113:1-7. http://dx.doi.org/10.1016/j.petrol.2013.12.005
Zemel, B. 1995. Tracers in the Oil Field, Elsevier Science B. V., T
Zhou, H., Gómez-hernández, J. J., and Li, L. (2014). Advances in Water Resources Inverse
methods in hydrogeology : Evolution and recent trends. Advances in Water Resources, 63,
22–37. http://doi.org/10.1016/j.advwatres.2013.10.014

121
Zhou H, JJ. Gómez-Hernández, H.J. Hendricks Franssen, L. Li (2011). An approach to handling
Non-Gaussianity of parameters and state variables in ensemble Kalman filter. Adv Water
Resour. 34(7):844–64. http://dx.doi.org/10.1016/j.advwatres.2011.04.014.
Zinn, B., and Harvey, C.F. (2003). When good statistical models of aquifer heterogeneity go bad:
A comparison of flow, dispersion, and mass transfer in connected and multivariate
Gaussian hydraulic conductivity fields. Water Resources Research 39(3), 1051,
doi:10.1029/2001WR001146.

122
Tables

Table 1. List of uncertain parameters
Planar Model Discrete Fracture Network
Parameter Range Parameter Range
Half-length (ft) 16-320 Fracture density 0-4
Conductivity (md) 2-10000 Half-length (ft) 80-370
Decline rate (1/Psi) 1.5e-4 – 0.01 Conductivity HF (md) 7-10000
Matrix permeability (nd) 0.3-130 Conductivity NF (md) 2-10000
Matrix porosity 0.03-0.1 Decline rate HF (1/Psi) 1.5e-4 – 0.01
Decline rate NF(1/Psi) 2.8e-4 – 0.1
Matrix permeability (nd) 0.3-130
Matrix porosity 0.03-0.1

123

Table 2. True fracture parameters (Experiment 1, 2)
Stage
#
Fracture Parameter
Half Length (ft) Initial Permeability
(md)
Decline Rate
(1/Psi)
1 96 113 0.0015
2
360 1425 0.0012
3 220 6762 0.0011
4 136 9839 0.001

124

Table 3. Simulation Rock/fluid properties
Fluid density (kg/m
3
) 1000
Fluid viscosity (Pa.s) 10
-3

Matrix Porosity 0.06
Fracture Porosity 0.2
Matrix Permeability (nd) 10
Production BHP (Psi) 500
Poisson’s ratio 0.3-0.5
Biot’s constant 0.1-0.3
Permeability modulus (Psi
-1
) 10
−3
− 10
−2

125

Table 4. True fracture parameters (Experiment 3, 4)
Stage
#
Fracture Parameter
Fracture
Density
Hydraulic Fracture
Conductivity (md)
Half-
length (ft)
Natural Fracture
Conductivity (md)
Decline rate
HF (1/Psi)
Decline rate
NF (1/Psi)
1 0.1 12 124 1819 0.0015 0.001
2 0.5 224 172 1161 0.0012 0.0001
3 1.2 917 268 9 0.0011 0.01
4 3 144 79 274 0.001 0.0015

126
Table 5. Geochemical Reactions-mineral dissolution-precipitation reactions
Reaction No. Reaction Log-Equilibrium-Constant
(log(Keq))
1 H2O = H+ + OH-

-13.2631
2 CO2 + H2O = H+ + HCO3—

-6.3221
3 CO2 + H2O = 2H+ + CO3—

-16.5563
4 CALCITE + H+ = (Ca++) + (HCO3-)

1.3560
5 KAOLINITE + 6(H+) = 5(H2O) + 2SiO2(aq) + 2(Al+++)

5.4706
6 ANORTHITE + 8H+ = 4(H2O) + (Ca++) + 2[SiO2(aq)] + 2(Al++)

23.0603

127

Table 6. Rock geomechanical properties
Rock property Facies 1 (yellow, sand) Facies 2 (grey, shale)
Young’s Modulus (kpa) 4.2e6 - 4.2e7 1.4e5 - 1.4e6
Poisson’s ratio 0.1 - 0.18 0.15 - 0.2
Cohesion (kpa) 2.1e3 - 4.9e3 50 - 600
Friction angle 20 - 40 15 - 25
Permeability (mD) 20 - 300 0.1 5
Porosity 0.25 - 0.45 0.25 - 0.35

128
Figures

Figure 1. Workflow relating fracture parameters to observed tracer and production data.

Hydraulic Fractures
Updated Fracture
Parameter
Economides et al., 2011
Woodroof et al., 2003
Tracer Concentration
Profile
Production Data
Observed
Data
Data
Assimilation
Data Acquisition
Data
Prediction
per Stage
Tracer Data Production
Data
Predicted
Data
Modeling and predictions

129

(a) Plane model.

(b) Discrete fracture network model.

Figure 2. Shale reservoir model with four stages of hydraulic fracturing: (a) Planar Model, (b) DFN

130

(a) Fracture conductivity change

(b) Pressure and cumulative production

Figure 3. a. Fracture conductivity changes for one single fracture element (planar model) with pressure (first
row, left), changes with time (first row, right). b. Effect of facture conductivity reduction on pressure (second
row, left) and on cumulative production (second row, right)

Log-Conductivity
Time (month)
Log-Conductivity
Pressure (psi)
Time (month)
Cumulative Oil Pro (STB)
Time (month)
Pressure (Psi)

131

(a) Initial permeability map (b) Final permeability map (12 month)

Figure 4. Map of hydraulic and natural fracture conductivity changes with time.

132

Figure 5. Tracer concentration travel distance (left), showing the schematic representation of tracer profile
(natural logarithm of tracer concentration is shown); concentration-distance profile for different stages (right).

(a) Tracer Travel
Distance
(b) Concentration
Profile
Tracer Concentration
Distance (ft)

133

(a) Fracture Model

Different Half-Length Different Conductivity Different Length and Conductivity

(b) Tracer Concentration

Normalized Concentration

Time (hr)

Time (hr)

Time (hr)

a. Cumulative Production

Cumulative Oil Pro (STB)

Time (day)

Time (day)

Time (day)
Figure 6. Sensitivity of predicted tracer concentration and production for four fracturing stages with different half-length
(left), different conductivity (middle), and different half-length and conductivity (right).

(a) Tracer Response

134

(b) Production Response

(c) Pressure Response

Figure 7. Planar model: Sensitivity of Tracer , production and pressure data to fracture and matrix
properties including hydraulic fracture conductivity (HFK), Half-length (HFL), Decline rate of conductivity
(a), Matrix permeability (mK) and matrix porosity (mphi).

135
(a) Parameters Histograms
Freq.

Half-Length (ft) Log K_HF (md) Decline Rate (×1000, 1/Psi)
Freq.

Log K_Matrix (md) Matrix porosity

(b) Well Responses
Tracer concentration

Cum oil pro(STB)

Pressure (Psi)

Time (hr) Time (day) Time (day)
Figure 8. Evolution of logarithm of fracture conductivity, half-length, decline rate, matrix permeability and matrix
porosity and corresponding predicted data by assimilation of the tracer, production and pressure data. The initial
ensemble of parameters are shown with black box, and updated ensemble by consecutive assimilating tracer data
and production (and pressure) data are shown with grey and blue histograms. The predictions with the initial and
updated ensemble are shown in grey and cyan (with ensemble means shown in black and blue, respectively);
reference predictions are shown in red.

136

(a) Tracer Response

(b) Production Response

(c) Pressure Response

Figure 9. Discrete fracture network model: Sensitivity of Tracer, production and pressure data to fracture
and matrix properties including fracture density, hydraulic fracture conductivity (HFK), Half-length (HFL),
natural fracture conductivity (FK), Decline rate of hydraulic fracture conductivity (aHF), Decline rate of
natural fracture conductivity (aF), Matrix permeability (mK) and matrix porosity (mphi).

137

(a) HF Conductivity (b) NF Conductivity (c) NF Density
Freq.

Freq.

Freq.

Production difference (H-M) Production difference (H-M) Production difference (H-M)

Freq.

Freq.

Freq.

Production difference (L-M) Production difference (L-M) Production difference (L-M)

Figure 10. Sensitivity of production data to hydraulic and natural fracture conductivity compared to fracture
density. “H”, “M”, and “L” refer to high, medium and low values of respective parameters.

138

(a) Parameters Histograms

Freq.

NF Density HF Half-Length (ft) Log K_HF (md) Log K_NF (md)
Freq.

HF decline (×1000, 1/Psi) NF Decline (×100, 1/Psi) Log K _Matrix (md) Matrix porosity

(b) Well Responses
Tracer concentration

Cum oil pro(STB)

Pressure (Psi)

Time (hr) Time (day) Time (day)

Figure 11. Evolution of fracture density, half-length, logarithm of hydraulic and natural fracture conductivity,
hydraulic fracture and natural fracture decline rates, logarithm of matrix permeability and matrix porosity and
corresponding predicted data by consecutive assimilation of the production and pressure data. The initial ensemble
of parameters are shown with black box, and updated ensemble by assimilating tracer data and production (and
pressure) data are shown with grey and blue histograms. The predictions with the initial and updated ensemble are
shown in grey and cyan (with ensemble means shown in black and blue, respectively); reference predictions are
shown in red.

139

Figure 12. Illustration of forward and inverse distance transform for a binary image. The example on the top
illustrates the concept while the one on the bottom shows application to a facies map.

140

Figure 13. Schematic of the ES-DT and ES-DT-SVD implementation workflows for models with three facies;
dashed box is only included for the ES-DT-SVD method where the SVD coefficients of the distance maps are
updated using the ES algorithm.

141

Figure 14. The error associated with back transformation threshold and the number of SVD basis elements
used for distance transformation. The error is reduced with increasing number of SVD basis elements and
decreasing threshold value.

142

Figure 15. Reference map, well configurations and pollutant concentration for Experiments 1; top row:
reference model and well configuration (left), the corresponding initial pollutant concentrations (middle), and
final concentrations (right) for Cases 1 and 2; bottom row: reference model and well configuration for Cases
3 and 4.

143

Figure 16. Updated log-permeability maps for Experiment 1 with two facies types using different data
assimilation techniques; the top set of plots show the updates to a sample realization for Cases 1-4 (rows);
Columns 1-5 show the a sample from the initial, ES, ES-DT, ES-SVD, and ES-DT-SVD updated realization.
The plot on the bottom show similar results for the ensemble mean.

144

Figure 17. Evolution of ensemble-based RMSE of the log-permeability maps (top) and pressure head and
concentrations data for Case 1-4 in Example 1.

145

Figure 18. Ensemble-based predictions of pressure and concentration for Cases 1-4 using different
assimilation schemes in Example 1: (a) pressure predictions corresponding to, clockwise from top left, ES,
ES-DT, ES-SVD, and ES-DT-SVD, respectively; (b) concentration predictions corresponding to, top to
bottom, ES, ES-DT, ES-SVD, and ES-DT-SVD, respectively; the predictions with the initial and updated
ensemble are shown in grey and cyan (with ensemble means shown in black and blue, respectively); reference
predictions are shown in red.

146

Figure 19. Experimental setup and log-permeability estimation results for Example 2 with three facies types;
top plots: reference log-permeability map with well configuration (left), and the corresponding initial
(middle) and final (right) pollutant concentration maps; bottom plots: two sample realizations (rows 1 and 2),
mean (row 3), and variance (row 4) of the initial, ES, ES-DT, ES-SVD-, and ES-DT-SVD updated ensembles
(Columns 1-5, respectively).

147

Figure 20. Data prediction results for Example 2 with three facies; top: ensemble-based predictions of
pressure data at all wells using the ES-DT-SVD method; middle: ensemble-based prediction of the
concentration at extraction wells using the ES-DT-SVD method; bottom: evolution of ensemble-based RMSE
of log-permeability (left) and data predictions (right) using ES, ES-DT, ES-SVD, and ES-DT-SVD methods.

148

Figure 21. Sensitivity of pressure head and concentration to the facies permeability for each well.

True facies map

149
Initial ES ES-DT ES-SVD ES-DT-SVD
Sample

Sample

Mean

Variance

Figure 22. log-permeability estimation results for Example 3 with three facies types: two sample realizations
(rows 1 and 2), mean (row 3), and variance (row 4) of the initial, ES, ES-DT, ES-SVD-, and ES-DT-SVD
updated ensembles (Columns 1-5, respectively).

Figure 23. Estimated value of permeability for facies with low permeability (left), with medium permeability
(middle) and high permeability (right) using ES-DT-SVD method.

150

Figure 24. Data prediction results for Example 3 with three facies; top: ensemble-based predictions of
pressure data at all wells using the ES-DT-SVD method; middle: ensemble-based prediction of the
concentration at extraction wells using the ES-DT-SVD method; bottom: evolution of ensemble-based RMSE
of log-permeability (left) and data predictions (right) using ES, ES-DT, ES-SVD, and ES-DT-SVD methods.
Facies map
map
Well responses
Iteration Iteration
RMSE

151

Figure 25. 3D coupled flow and geomechanical model of an aquifer and its boundary conditions.

152

Figure 26. Approximated 3D object based model with two rock facies type for a discrete fracture network
(DFN). a). DFN model, b). proxy lump object based model. c). 3D connected lump objects.

153

Figure 27. Rock facies and petrophysical properties map.

154
a). Rock facies distribution b). Pressure (*10^4 kpa)

c). Displacement in Z direction (m) d). Minimum stress (kpa)

Figure 28. 3D Forward simulation model. a). rock facies distribution, b). Pressure, c). Displacement in y
direction, d). Stress map.

155

Figure 29. Simple illustration of the Mohr-Coulomb’s failure criterion.

156
(a) High-permeability rock

(b) Low-permeability rock

Figure 30. Coupled flow and geomechanical model response for high (a) and low (b) permeability rocks. The
following quantities are shown: permeability (top-left), Young’s modulus (top-middle), vertical displacement
(top-right), min stress (bottom-left), max stress (bottom-middle), and failure map (bottom-right).

157

Figure 31. Illustration of microseismic data acquisition (Surface and downhole receiver array).

158
Permeability Pressure Min stress Max stress Failure

(*10^4)

(*10^3)

(*10^4)

Sample 1

Sample 2

a).
Cohesion Pressure Min stress Max stress Failure

(*10^4)

(*10^3) (*10^4)

Sample 1

Sample 2

b)
Permeability Pressure Min stress Max stress Failure

(*10^4)

(*10^3)

(*10^4)

Sample 1

159
Sample 2

c).
Figure 32. Effect of Rock properties on simulated failure map (potential microseismic map). a). Two samples
of rock facies map along with their pressure and stress state. b). Two samples with similar facies map but
various constant rock property within each facies. c). Two samples with various within facies heterogeneity.

160

a). b).
Figure 33. Effect of uncertainty in rock property on a). Mohr-Coloumb failure line and b). Mohr circle.

161

Figure 34. Relationship between Mohr-Coulomb’s failure criterion and sigmoid function.

162

Figure 35. Simulated MEQ cloud using proposed stochastic model for simple homogeneous case.

163

a).

b).
Figure 36. Simulated MEQ cloud using both Deterministic and Stochastic approach. a). Rock with low rock
strength. b). Rock with high strength.

164

a).
Failure Probability MEQ Map (sample #1) MEQ Map (sample #2)

b).
Failure Probability MEQ Map (sample #1) MEQ Map (sample #2)

c).
Figure 37. a). Illustration of elipse around event location to account for error in recorded data. b). simulated
MEQ cloud using stochastic approach with consideration of error in recorded data for rock with low
strength. c). simulated MEQ cloud for rock with high strength.

165

Figure 38. Sensitivity of MEQ events to rock properties of geological model with one single facies type; C, E,
F, poi, K and pro stand for the cohesion, Young’s modulus, friction angle, Poisson’s ratio, permeability and
porosity, respectively.

166

Figure 39. Sensitivity of MEQ events to rock properties of geological map with two rock facies type (sand and
shale); C, E, F, poi, K and pro stand for the cohesion, Young’s modulus, friction angle, Poisson’s ratio,
permeability and porosity, respectively. Subscripts s refers to sand.

167

Figure 40. Proposed workflow for transformation and assimilation of discrete microseismic cloud.

168

Figure 41. Schematic representation for continuous transformation of discrete microseismic data in one
dimension. Blue crosses represent the location of MEQ events. Black symmetric curves display the Gaussian
kernels. The red curve shows spatial density of MEQ events as a continuous representation of MEQ cloud.

169

Figure 42. Two-dimensional Schematic representation for continuous transformation of discrete microseismic
data using kernel density function.

170

Figure 43. (top) Facies distribution of the fourth layer of the 3D model; (middle and bottom) the cohesion,
permeability, Young’s Modulus, friction angle and Poisson’s ratio before and after assimilation of discrete
microseismic data. The initial ensemble of parameters and updated ensemble are shown with grey and blue
histograms. The value of the reference parameter is shown with the red line. The initial and updated value of
two samples are shown with * and o markers. The initial and updated parameters for Sample 1 are shown in
black and cyan, while the initial and updated parameters of Sample 2 are shown in red and green,
respectively.

Freq.

Sand Cohesion (kpa) Sand Poisson’s ratio Sand Young’s Modulus (*10^6)
Freq.

Shale Cohesion (kpa) Shale Poisson’s ratio Shale Young’s Modulus (*10^5)
Freq.

Sand Friction Angle Sand Perm (md)
Freq.

Shale Friction Angle Shale Perm (md)

Freq.

Sand Cohesion (kpa) Sand Poisson’s ratio Sand Young’s Modulus (*10^6)
Freq.

Shale Cohesion (kpa) Shale Poisson’s ratio Shale Young’s Modulus (*10^5)
Freq.

Sand Friction Angle Sand Perm (md)
Freq.

Shale Friction Angle Shale Perm (md)

171

Figure 44. Stochastic simulated MEQ cloud before and after update for experiment 1. (top) Reference MEQ
distribution and its corresponding KDE transformation. a). MEQ events distribution, b). KDE
transformation (continuous domain).

Reference MEQ Map Reference KDE Map

Y(Km)

X(Km)
Sample 1 Sample 2 Mean Variance
Initial

Update

a).
Initial

Update

172

Figure 45. Proposed workflow for discrete rock facies model calibration by integration of discrete
microseismic data.

173

Figure 46. Evolution of the geological rock facies map by integrating discrete microseismic data.

Reference Map
Y(Km)

X(Km)

Sample 1 Sample 2 Mean Variance
Initial

Update

Reference Map
Y(Km)

X(Km)

Sample 1 Sample 2 Mean Variance
Initial

Update

174

Figure 47. Stochastically simulated microseismic clouds before and after update for Experiment 2. (top)
Reference microseismic distribution and its corresponding KDE transformation; (Rows 2 and 3)
microseismic events distribution, including two samples and mean and variance maps; (Rows 3 and 4) KDE
transformation results (continuous domain), including two samples, the mean and variance maps.

Reference MEQ Map Reference KDE Map

Y(Km)

X(Km)
Sample 1 Sample 2 Mean Variance
Initial

Update

a).
Initial

Update

175

Figure 48. Evolution of the geological rock facies map by integrating discrete microseismic data for
experiment 3.

Reference Map
Y(Km)

X(Km)

Sample 1 Sample 2 Mean Variance
Initial

Update

176

Figure 49. Experiment 3: (top) Facies distribution of the fourth layer of the 3D model; (middle and bottom)
the cohesion, permeability, Young’s Modulus, friction angle and Poisson’s ratio before and after assimilation
of discrete microseismic data. The initial ensemble of parameters and updated ensemble are shown with grey
and blue histograms. The value of the reference parameter is shown with the red line. The initial and updated
value of two samples are shown with * and o markers. The initial and updated parameters for Sample 1 are
shown in black and cyan, while the initial and updated parameters of Sample 2 are shown in red and green,
respectively.

Freq.

Sand Cohesion (kpa) Sand Poisson’s ratio Sand Young’s Modulus (*10^6)
Freq.

Shale Cohesion (kpa) Shale Poisson’s ratio Shale Young’s Modulus (*10^5)
Freq.

Sand Friction Angle Sand Perm (md)
Freq.

Shale Friction Angle Shale Perm (md)

Freq.

Sand Cohesion (kpa) Sand Poisson’s ratio Sand Young’s Modulus (*10^6)
Freq.

Shale Cohesion (kpa) Shale Poisson’s ratio Shale Young’s Modulus (*10^5)
Freq.

Sand Friction Angle Sand Perm (md)
Freq.

Shale Friction Angle Shale Perm (md)

177

Figure 50. Stochastically simulated microseismic clouds before and after update for Experiment 2. (top)
Reference microseismic distribution and its corresponding KDE transformation; (Rows 2 and 3)
microseismic events distribution, including two samples and mean and variance maps; (Rows 3 and 4) KDE
transformation results (continuous domain), including two samples, the mean and variance maps.

Reference MEQ Map Reference KDE Map

Y(Km)

X(Km)
Sample 1 Sample 2 Mean Variance
Initial

Update

a).
Initial

Update

178
(a) (b)

Figure B-1. Schematic representation of fracture and matrix grid blocks for transmissibility calculation in x
and y direction; a). Adjacent fracture grid blocks (left), b). Adjacent fracture and matrix grid blocks (right).

(a) Pressure distribution (b) Flow regimes in our simulation
L=dy
w=dx
dx2
dx1
dy

179
Y

RNP and RNP derivative

X True Time (hr)

Figure B-2. Different flow regimes in production from fractured reservoirs.

180
Cumulative Oil Pro (STB)

Time (day)

Figure B-3. Production sensitivity to fracture aperture cell size: Star (*), solid circle (o), and dashed line (--) plot
represent the production response for a fracture aperture equal to 3 ft, 1ft, and 0.1 ft.

181

Figure B-4. Comparison between the pressure (middle) and production (right) responses from our modeling approach
(top) with those from the MRST simulation package with embedded fracture models (bottom).

Oil production rate
Fracture model
Pressure map

Abstract (if available)

Abstract Construction of predictive subsurface flow models involves subjective interpretation and interpolation of spatially limited data, often using imperfect modeling assumptions. Hence, the process can introduce significant uncertainty and bias in predicting the flow and transport behavior of these systems. The uncertainty in the distribution of discrete fracture network and rock facies distribution in complex geologic environments along with their flow and mechanical properties can be consequential for forecasting the dynamic response of these systems to perturbations due to fluid injection and development activities. In this work, tracer test, performance and microseismic data are used for dynamic characterization of important parameters of dynamic reservoir models, including (1) matrix fluid and mechanical properties, (2) rock facies distribution (3) planar fracture half-length and hydraulic conductivity, (4) discrete fracture networks density and conductivity, and (5) fracture closing (conductivity decline) rate during production. Conventional model calibration techniques are designed to update continuous model parameters and have difficulty in estimating discrete parameters such as rock facies distribution (or distribution of discrete fracture network). Here, we present a distance transform approach for converting discrete parameters to continuous parameters that can be updated, using continuous model calibration methods, from flow and transport-related data. Distance-based transforms are widely used in discrete image processing, where the discrete values in each pixel are replaced with their distance (i.e., a continuous variable) to the nearest cell with a different value (i.e., nearest boundary cell). The resulting continuous distance map can then be back-transformed to retrieve the original discrete image. For facies model calibration, the facies maps are first transformed into continuous spatial maps of distances to facilitate model calibration with continuous inversion methods. For facies with large-scale connectivity, truncated singular value decomposition (SVD) parametrization can be used to represent the distance maps with low rank parameters. In this work, a variant of the ensemble smoother is used to update the final continuous parameters (either distance maps or their SVD coefficients if used). The resulting updated parameters are back-transformed to obtain calibrated models. The distance transform method addresses an important problem in facies model calibration where model updating can result in losing facies connectivity and discreteness. ❧ In general, high pressure fluid injection into geologic formations can trigger rock deformation and failure and lead to microseismic events. Injection induced seismicity has been proposed as a monitoring technique that can be used to constrain the description of rock flow and mechanical properties. To establish the correlations between rock properties and microseismic clouds, the developed framework combines coupled flow and geomechanics simulation outputs with Mohr-Coulomb’s failure criterion to describe the spatiotemporal distribution of seismicity potential in the formation. This model is then used with an ensemble smoother with multiple data assimilation, ES-MDA to estimate rock properties form observed microseismic data. Since the ensemble smoother is designed for assimilation of continuous data, kernel density estimation (KDE) is applied to convert discrete microseismic events to a continuous map of seismicity density. We present the developed formulation and discuss its application to inversion of microseismic data for characterization of reservoir flow and mechanical properties. Several examples are presented to evaluate the performance of the method for estimation of rock properties from microseismic data.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

$History matching for fractured reservoirs with Hough-transform-based parameterization$

PDF

History matching for fractured reservoirs with Hough-transform-based parameterization

PDF

Subsurface model calibration for complex facies models

PDF

Uncertainty quantification and data assimilation via transform process for strongly nonlinear problems

PDF

Assessing induced seismicity rate increase based on deterministic and probabilistic modeling

PDF

Integrated reservoir characterization for unconventional reservoirs using seismic, microseismic and well log data

PDF

Modeling and simulation of complex recovery processes

PDF

Real-time reservoir characterization and optimization during immiscible displacement processes

$Integrated workflow for characterizing and modeling fracture network in unconventional reservoirs using microseismic data$

PDF

Integrated workflow for characterizing and modeling fracture network in unconventional reservoirs using microseismic data

$An extended finite element method based modeling of hydraulic fracturing$

PDF

An extended finite element method based modeling of hydraulic fracturing

PDF

Deep learning for subsurface characterization and forecasting

PDF

Deep learning architectures for characterization and forecasting of fluid flow in subsurface systems

PDF

A method for characterizing reservoir heterogeneity using continuous measurement of tracer data

PDF

Stochastic data assimilation with application to multi-phase flow and health monitoring problems

PDF

Continuum modeling of reservoir permeability enhancement and rock degradation during pressurized injection

PDF

The role of counter-current flow in the modeling and simulation of multi-phase flow in porous media

PDF

Dynamics of CO₂ leakage through faults

$A coupled poromechanics-damage mechanics approach to model fracturing in multiphase flow$

PDF

A coupled poromechanics-damage mechanics approach to model fracturing in multiphase flow

PDF

Sparse feature learning for dynamic subsurface imaging

PDF

Reactivation of multiple faults in oilfields with injection and production

$Coarse-scale simulation of heterogeneous reservoirs and multi-fractured horizontal wells$

PDF

Coarse-scale simulation of heterogeneous reservoirs and multi-fractured horizontal wells

Asset Metadata

Creator Hakim Elahi, Siavash (author)

Core Title Integration of multi-physics data into dynamic reservoir model with ensemble Kalman filter

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Petroleum Engineering

Publication Date 04/10/2018

Defense Date 11/21/2017

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag data assimilation,dynamic reservoir model,ensemble Kalman filter,multi-physics data,OAI-PMH Harvest

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Jafarpour, Behnam (committee chair), Ershaghi, Iraj (committee member), Sammis, Charles (committee member)

Creator Email hakimela@usc.edu,siavashhakimds@gmail.com

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c89-11278

Unique identifier UC11670341

Identifier etd-HakimElahi-6186.pdf (filename),usctheses-c89-11278 (legacy record id)

Legacy Identifier etd-HakimElahi-6186.pdf

Dmrecord 11278

Document Type Dissertation

Rights Hakim Elahi, Siavash

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA