US011125905B2

( 12 ) United States Patent ( 10) Patent No.: US 11,125,905 B2

Cancelliere et al . (45) Date of Patent : Sep. 21 , 2021
( 54 ) METHODS FOR AUTOMATED HISTORY 7,531,791 B2 5/2009 Bryman
MATCHING UTILIZING MUON 8,288,721 B2 10/2012 Morris et al .
8,384,017 B2 2/2013 Botto
TOMOGRAPHY 8,536,527 B2 9/2013 Morris et al .
9,423,362 B2 8/2016 Sossong et al .
( 71 ) Applicant: SAUDI ARABIAN OIL COMPANY, 9,639,973
9,746,580
B2
B2
5/2017 Bai et al .
8/2017 Hayes et al .
Dhahran (SA ) 9,784,859 B2 10/2017 Blanpied et al .
9,915,626 B2 3/2018 Blanpied et al .
( 72 ) Inventors: Michel Cancelliere , Dhahran ( SA) ; 9,939,537 B2 4/2018 Sossong
Ozgur Kirlangic , Dhahran ( SA ) 10,042,079 B2 8/2018 Patnaik
(Continued )
( 73 ) Assignee : Saudi Arabian Oil Company , Dhahran
( SA) FOREIGN PATENT DOCUMENTS
(*) Notice: Subject to any disclaimer, the term of this JP
WO
2013002830 A
2013155075 Al
1/2013
10/2013
patent is extended or adjusted under 35 WO 2017194647 A1 11/2017
U.S.C. 154 ( b ) by 272 days .
( 21 ) Appl. No .: 16 /402,317 OTHER PUBLICATIONS
( 22 ) Filed : May 3 , 2019 Antonuccio , et al . , The Muon Portal Project: Development of an
Innovative Scanning Portal Based on Muon Tomography, 2013 3rd
( 65 ) Prior Publication Data International Conference on ANIMMA, 978-1-4799-1047 - Feb .
2013 , © 2013 IEEE .
US 2020/0348440 A1 Nov. 5 , 2020 (Continued )
( 51 ) Int . Cl . Primary Examiner Roy Y Yi
GOIV 5/00 ( 2006.01 ) (74 ) Attorney, Agent, or Firm Bracewell; Constance G.
=

GOIV 99/00 ( 2009.01 ) Rhebergen ; Christopher L. Drymalla

( 52 ) U.S. CI .
CPC GO1V 5/0075 (2013.01 ) ; GOIV 99/005 ( 57 ) ABSTRACT
(2013.01 ) Embodiments provide a method for surveying a hydrocar
( 58 ) Field of Classification Search bon reservoir utilizing a reservoir model . The method
None includes the step of establishing an ensemble of models
See application file for complete search history. reflecting attributes of the hydrocarbon reservoir based on
(56 ) References Cited the reservoir model in its present state . The method includes
the step of updating the reservoir model by utilizing a
U.S. PATENT DOCUMENTS volumetric density image of the hydrocarbon reservoir. The
volumetric density image can be constructed via muon
4,504,438 A * 3/1985 Levy GO1V 5/04 tomography
376/156
7,488,934 B2 2/2009 Bryman 23 Claims , 2 Drawing Sheets
1
140 1 1

140 110
s
1

130

150

120
 US 11,125,905 B2
Page 2

( 56 ) References Cited Jaenisch, et al . , Real Time Muon Tomography Imaging Simulation

and Fast Threat Target Identification , Proceedings of SPIE — The
U.S. PATENT DOCUMENTS International Society for Optical Engineering 7310 , May 2009 .
Jiang, et al . , Modelling and Monitoring of Geological Carbon
2008/0128604 A1 * 6/2008 Bryman GO1T 1/203 Storage , Applied Energy 112 ( 2013 ) 784-792 .
2008/0315091 A1 12/2008 Morris et al .
250/266 Kudryavtsev , et al ., Monitoring subsurface CO2 emplacement and
2010/0198570 A1 8/2010 Sarma et al .
security of storage using muon tomography, 11 International Journal
2011/0035151 A1 2/2011 Botto of Greenhouse Gas Control 21 ( 2012 ) .
2012/0215511 A1 * 8/2012 Sarma GO1V 1/308 Lesparre, et al . , 3D Density Imaging with Muon Flux Measurements
703/10 from Underground Galleries , Geophysical Journal International
2014/0091803 A1 4/2014 Dodds et al . Advance Access published Dec. 23 , 2016 , 36 pages.
2015/0247940 A1 9/2015 De Matos Ravanelli et al . Valestrand , et al . , The Effect of Including Tracer Data in the EnKF
2016/0018541 A1 1/2016 Stefani Approach , SPE - 113440 ( 2008 ) .
2016/0116630 A1 4/2016 Sossong Zagayevskiy , et al . , Assimilation of Time-Lapse Temperature Obser
2017/0329039 A1 11/2017 Kang et al . vations and 4D - Seismic Data With the EnKF in SAGD Petroleum
Reservoirs, SPE - 174547 ( 2015 ) .
OTHER PUBLICATIONS Zhong , et al . , A case study of using cosmic ray muons to monitor
supercritical CO2 migration in geological formations, 185 Applied
Energy 1450 ( 2017 ) .
Bonneville , Borehole Muon Detector for 4D Density Tomography Dabboor, O. , “ On the Detectability of Density Distributions of
of Subsurface Reservoirs, Presentation, Carbon Storage and Oil and Asgteroids and SAGD Reservoirs Using Gravimetry and Muon
Natural Gas Technologies Review Meeting ( 2016 ) . Tomography, ” Queen's University, a Thesis to the Department of
Bonneville , et al . , A novel muon detector for borehole density Geological Sciences and Geological Engineering, Oct. 31 , 2018 ,
tomography, 851 Nuclear Instruments and Methods in Physics 144 pages.
Research Section A 108 ( 2017 ) . International Search Report and Written Opinion for PCT Applica
Jaenisch , et al., Muon Imaging and Data Modeling , Proceedings of tion No. PCT /US2020 /031167 dated Sep. 18 , 2020 , 23 pages .
SPIE—The International Society for Optical Engineering , May
2007 . * cited by examiner
U.S. Patent Sep. 21 , 2021 Sheet 1 of 2 US 11,125,905 B2

OL
150
120

140 140 1B
.
FIG
1
TE

ra
}

130

LL L E
150

140
1A
.
FIG

110 130 120

U.S. Patent Sep. 21 , 2021 Sheet 2 of 2 US 11,125,905 B2

200

210 220
*** + ***********
230
CE *****

M M M MAN www

FIG . 2

360
00000000

370

310 320
- 330 350
M

340

FIG . 3
 US 11,125,905 B2
1 2
METHODS FOR AUTOMATED HISTORY minimized . The minimization of the cost function can be
MATCHING UTILIZING MUON obtained by applying an optimization algorithm . Although
TOMOGRAPHY techniques such as optimization and non - linear program
ming are not new in the art, the selection of the most
BACKGROUND 5
adequate optimization algorithm is not trivial, and the num
ber of independent variables involved in complex reservoir
Field of the Disclosure simulation does not make the solution of the optimization
problem a standard procedure.
Embodiments of the disclosure generally relate to hydro- 10 Muons originating from cosmic rays can generate time
carbon reservoir modelling. dependent images of reservoirs due to their ability to pen
Description of the Related Art etrate into the subsurface in the order of hundreds of meters .
These real time volumetric density images (or “ 4D data ” )
The development of a hydrocarbon reservoir may require matically from
obtained such muon sources can be utilized to auto
an understanding of the reservoir geology for appropriately 15 voir models . Suchorcosmic
update calibrate numerically calculated reser
ray muon tomography ( or muog
selecting means for hydrocarbon extraction and adequately raphy ) can be complemented with history matching
placing production wells , injection wells , and surface facili algorithms such as Ensemble Kalman Filtering (EnKF ) or
ties . Various reservoir modelling and visualization tools
have been developed to better understand the reservoir Ensemble Kalman Smoothing (EnKS).
geology in order to make decisions regarding reservoir 20 Advantageously, muography based history matching
development. Such modelling and visualization tools increases the accuracy of the calibration , provide multiple
include certain computational means to generate production calibrated models , and ease the history matching process ,
forecasts and optimize field performance . Typically, actual thus allowing more time and computing power for optimi
data are acquired at various observation locations of the zation versus not using such history matching techniques.
reservoir and are integrated into a reservoir simulation 25 Advantageously, an accurate reservoir model can be
model . This way, discrepancies between the reservoir simu- achieved using muography based history matching tech
lation model and the real reservoir can be reduced . niques even when the muographic images have low resolu
Conventional history matching techniques are used for tion , that is , the muographic images are pixelated . This is
identifying and minimizing such discrepancies, where a due to the muographic 4D data being continuously, or at
reservoir simulation model is iteratively calibrated until the 30 improve , fed into the history matching loops to
least periodically
the calibration process .
calibration results in an acceptable reproduction of reservoir
behavior prior to the adjustment. The accuracy of these Muon detectors combined with muon tomography allow
history matching techniques depends on the quality of the volumetric imaging of the subterranean density distribution
reservoir simulation model and the quality and quantity of and its evolution in real time in hydrocarbon reservoirs ,
the actual data . A history matched reservoir model can be 35 which may serve as an effective surveying method compared
used to predict future reservoir behavior with an enhanced to conventional methods such as 4D seismology. Seismic 4D
degree of confidence, especially in cases where the adjust- data are collected at specific points in time , where the
ments are constrained by actual geological properties of the sensitivity of the seismic data with respect to the varying
reservoir and actual fluid properties in the reservoir . density is lesser than that of muographic 4D data . In
40 addition , 4D seismology requires the generation of seismic
SUMMARY waves , which is not a requirement for 4D muography.
Embodiments of the disclosure provide a method for
Embodiments of the disclosure generally relate to hydro- surveying a hydrocarbon reservoir utilizing a reservoir
carbon reservoir modelling utilizing muon tomography. model. The method includes the step of establishing an
More specifically, embodiments of the disclosure relate to 45 ensemble of models reflecting attributes of the hydrocarbon
methods for modelling a hydrocarbon reservoir utilizing reservoir based on the reservoir model in its present state .
muon tomography by detecting muons and constructing The method includes the step of updating the reservoir
volumetric density images of a hydrocarbon reservoir to model by utilizing a volumetric density image of the hydro
enhance the hydrocarbon reservoir model . carbon reservoir.
Reservoir modelling involves inherent errors and approxi- 50 In some embodiments, the volumetric density image is
mations that result in uncertainties that cannot be fully constructed via muon tomography. In some embodiments,
eliminated . Because the actual data that can be acquired are the method further includes the step of forecasting a hydro
limited , the actual data cannot be solely relied on to visualize carbon production from the hydrocarbon reservoir with the
the entire spatial configuration of the reservoir model . reservoir model in the updating step . In some embodiments,
Conventional history matching techniques are used for 55 the method further includes the step of imaging the reservoir
identifying and minimizing such uncertainties. However, model in the updating step . In some embodiments, the
conventional history matching techniques require trial and updating step includes performing EnKF or EnKS . In some
error and can be consuming in terms of time and computing embodiments, the attributes of the hydrocarbon reservoir
power . Also , conventional history matching techniques may include static data or dynamic data . In some embodiments,
result in the calibration of the reservoir simulation model 60 the updating step includes providing production data to
deriving non - unique solutions . update the reservoir model . In some embodiments, the
Automated history matching techniques are introduced to method further includes the step of comparing the produc
alleviate the drawbacks associated with conventional history tion data and the reservoir model in the updating step . In
matching . Automated history matching techniques are used some embodiments , the comparing step includes correcting
to treat history matching as an optimization process, where 65 input parameters of the reservoir model when a difference
a cost function is defined representative of the discrepancy between the production data and the reservoir model devi
between actual and simulated data , and the cost function is ates from a predetermined value .
 US 11,125,905 B2
3 4
Embodiments of the disclosure also provide a method for density image . In some embodiments, the method further
surveying a hydrocarbon reservoir utilizing a reservoir includes the step of altering the detection period based on a
model . The method includes the step of positioning a muon resolution of the volumetric density image .
detector at a subterranean location proximate the hydrocar
bon reservoir. The method includes the step of collecting 5 BRIEF DESCRIPTION OF THE DRAWINGS
muon detection data from the muon detector over a detection
period. The method includes the step of constructing a So that the manner in which the previously -recited fea
volumetric density image of the hydrocarbon reservoir by tures , aspects , and advantages of the embodiments of this
processing the muon detection data via tomography. disclosure as well as others that will become apparent are
In some embodiments, the method further includes the 10 attained and can be understood in detail , a more particular
step of updating the reservoir model by utilizing the volu- description of the disclosure briefly summarized previously
metric density image . In some embodiments, the method may be had by reference to the embodiments that are
further includes the step of forecasting a hydrocarbon pro- illustrated in the drawings that form a of this specifi
duction from the hydrocarbon reservoir with the reservoir cation . However, it is to be noted that the appended drawings
model in the updating step . In some embodiments, the 15 illustrate only certain embodiments of the disclosure and are
method further includes the step of imaging the reservoir not to be considered limiting of the disclosure's scope as the
model in the updating step . In some embodiments , the disclosure may admit to other equally effective embodi
method further includes the step of collecting a subsequent ments .
muon detection data from the muon detector after a time FIGS . 1A - B are graphical representations showing muons
increment . The method further includes the step of con- 20 penetrating into the earth's surface, in accordance with an
structing a subsequent volumetric density image of the embodiment of the disclosure .
hydrocarbon reservoir by processing the subsequent muon FIG . 2 is a schematic diagram for history matching of
detection data via tomography. The method further includes production data and muographic 4D data, in accordance
the step of adjusting the reservoir model in the updating step with an embodiment of the disclosure .
by utilizing the subsequent volumetric density image . In 25 FIG . 3 is a schematic diagram showing a muography
some embodiments, the muon detector includes a scintillator assisted history matching reservoir modelling workflow , in
or a drift tube . In some embodiments , the method further accordance with an embodiment of the disclosure .
includes the step of altering the detection period based on a In the accompanying Figures, similar components or
resolution of the volumetric density image . features, or both, may have a similar reference label .
Embodiments of the disclosure also provide a method for 30
surveying a hydrocarbon reservoir utilizing a reservoir DETAILED DESCRIPTION
model. The method includes the step of establishing an
ensemble of models reflecting attributes of the hydrocarbon The disclosure refers to particular features, including
reservoir based on the reservoir model in its present state . process or method steps. Those of skill in the art understand
The method includes the step of positioning a muon detector 35 that the disclosure is not limited to or by the description of
at a subterranean location proximate the hydrocarbon res- embodiments given in the specification. The subject matter
ervoir . The method includes the step of collecting muon of the disclosure is not restricted except only in the spirit of
detection data from the muon detector over a detection the specification and appended claims .
period. The method includes the step of constructing a Those of skill in the art also understand that the termi
volumetric density image of the hydrocarbon reservoir by 40 nology used for describing particular embodiments does not
processing the muon detection data via tomography. The limit the scope or breadth of the embodiments of the
method includes the step of updating the reservoir model by disclosure . In interpreting the specification and appended
utilizing the volumetric density image. claims , all terms should be interpreted in the broadest
In some embodiments, the method further includes the possible manner consistent with the context of each term . All
step of forecasting a hydrocarbon production from the 45 technical and scientific terms used in the specification and
hydrocarbon reservoir with the reservoir model in the updat- appended claims have the same meaning as commonly
ing step . In some embodiments , the method further includes understood by one of ordinary skill in the art to which this
the step of imaging the reservoir model in the updating step . disclosure belongs unless defined otherwise .
In some embodiments, the updating step includes perform- Although the disclosure has been described with respect
ing EnKF or EnKS. In some embodiments , the attributes of 50 to certain features , it should be understood that the features
the hydrocarbon reservoir include static data or dynamic and embodiments of the features can be combined with other
data . In some embodiments, the updating step includes features and embodiments of those features.
providing production data to update the reservoir model . In Although the disclosure has been described in detail , it
some embodiments, the method further includes the step of should be understood that various changes, substitutions,
comparing the production data and the reservoir model in 55 and alternations can be made without departing from the
the updating step . In some embodiments, the comparing step principle and scope of the disclosure . Accordingly, the scope
includes correcting input parameters of the reservoir model of the present disclosure should be determined by the
when a difference between the production data and the following claims and their appropriate legal equivalents.
reservoir model deviates from a predetermined value . In As used throughout the disclosure, the singular forms “ a ,"
some embodiments, the method further includes the step of 60 “ an ,” and “ the ” include plural references unless the context
collecting a subsequent muon detection data from the muon clearly indicates otherwise.
detector after a time increment. The method further includes As used throughout the disclosure, the word " about”
the step of constructing a subsequent volumetric density includes +/- 5 % of the cited magnitude.
image of the hydrocarbon reservoir by processing the sub- As used throughout the disclosure , the words “ comprise,”
sequent muon detection data via tomography. The method 65 “ has, ” “ includes ,” and all other grammatical variations are
further includes the step of adjusting the reservoir model in each intended to have an open , non- limiting meaning that
the updating step by utilizing the subsequent volumetric does not exclude additional elements , components or steps .
 US 11,125,905 B2
5 6
Embodiments of the disclosure may suitably “ comprise,” As used throughout the disclosure , the terms “ muon
" consist , " or " consist essentially of ” the limiting features tomography , " " cosmic ray muon tomography,” and “muog
disclosed, and may be practiced in the absence of a limiting raphy " are used interchangeably.
feature not disclosed . For example, it can be recognized by As used throughout the disclosure, the term “ resolution”
those skilled in the art that certain steps can be combined 5 refers to the degree of spatial detail an image holds . Accord
into a single step . ingly , a high resolution image includes more spatial detail
As used throughout the disclosure, the words " optional” than a low resolution image . The degree of spatial detail
or “ optionally ” means that the subsequently described event generally corresponds to a pixel count of the image . For
or circumstances can or may not occur . The description 10 example, one skilled in the art may recognize that images
includes instances where the event or circumstance occurs having less than about 75 pixels per inch (ppi ), alternately
and instances where it does not occur. less than about 150 ppi , or alternately less than about 300 ppi
Where a range of values is provided in the specification or are considered low resolution images. Respectively, one
in the appended claims , it is understood that the interval skilled in the art may recognize that images having greater
encompasses each intervening value between the upper limit 15 than about 75 ppi , alternately greater than about 150 ppi , or
and the lower limit as well as the upper limit and the lower alternately greater than about 300 ppi are considered high
limit . The disclosure encompasses and bounds smaller resolution images . As used throughout the disclosure , the
ranges of the interval subject to any specific exclusion term “ pixelated image ” refers to a low resolution image
provided . “ Substantial” means equal to or greater than 1 % having less than about 75 ppi , alternately less than about 150
by the indicated unit of measure . “ Significant ” means equal 20 ppi , or alternately less than about 300 ppi .
to or greater than 0.1 % of the indicated unit of measure . Muon Tomography
Where reference is made in the specification and Embodiments of the disclosure relate to using muons
appended claims to a method comprising two or more originating from cosmic rays to generate time-dependent
defined steps , the defined steps can be carried out in any volumetric images to automatically update or calibrate
order or simultaneously except where the context excludes 25 numerical reservoir models . Muons generated in the earth's
that possibility. atmosphere can penetrate several hundreds of meters into
As used throughout the disclosure, terms such as “ first ” the subsurface . The flux of muons generally decreases with
and " second " are arbitrarily assigned and are merely depth . Muons can be detected using scintillator arrays
intended to differentiate between two or more components placed in horizontal or vertical wells where each detected
of an apparatus. It is to be understood that the words “ first” 30 muon provides information about the average density of the
formation path it traveled . Each muon travels a different
and “ second” serve no other purpose and are not part of the path , which can be inverted using image tomography tech
name or description of the component, nor do they neces
sarily define a relative location or position of the component. niques to obtain a volumetric description of subsurface,
including geological structures and fluid movements . The
Furthermore , it is to be understood that that the mere use of 35 resolution of the obtained image depends on the detected
the term “ first ” and “ second ” does not require that there be muon flux , which generally decreases with depth . These
any “third” component, although that possibility is contem images can be continuously , or at least periodically, used to
plated under the scope of the present disclosure. track movement of fluids in the hydrocarbon reservoir.
As used throughout the disclosure, spatial terms described FIGS . 1A - B show muons penetrating into the earth's
the relative position of an object or a group of objects 40 surface 110 , in accordance with an embodiment of the
relative to another object or group of objects. The spatial disclosure . The muons may penetrate to various depths
relationships apply along vertical and horizontal axes . depending on their energy and the amount of material
As used throughout the disclosure , the term “ state data ” encountered along their paths. This results in a reduction of
include static data , dynamic data , and production data . muon flux . Surviving muons may reach a muon detector 120
As used throughout the disclosure , the term " static data ” 45 or an array of muon detectors 120 that are placed in a
refer to data such as permeability and porosity field data that borehole 130 along the trajectories 140 of the muons as
were conventionally considered not to vary with time . It shown as dotted lines . The muon trajectories 140 may
should be noted that static data can be updated with time in extend through a hydrocarbon reservoir 150 which has a
ensemble based Bayesian filtering methods. Static data may different density than that of the surrounding geology. FIG .
also include well log data , seismic data , core data , geophysi- 50 1A shows a borehole 130 that is substantially vertical. FIG .
cal data, and petrophysical data. 1B shows a borehole 130 that is substantially horizontal. In
As used throughout the disclosure , the term " dynamic some embodiments, more than one boreholes 130 may exist
data” refer to data such as pressure and phase saturation data to place multiple muon detectors 120 proximate the hydro
of the entire model that correspond to solutions of the flow carbon reservoir 150 for enhanced precision . For example ,
equations. Dynamic data may also include fluid flow prop- 55 two horizontal boreholes 130 can be drilled , one above the
erties such as viscosity and density. Dynamic data may also hydrocarbon reservoir 150 and one below the hydrocarbon
include pressure, volume , and temperature dependencies for reservoir 150. Muon detectors 120 can be placed in each
the fluid flow properties. borehole 130 such that muon detection data of multiple
As used throughout the disclosure, the term “ production muon detectors 120 can be compared to reconstruct a
data” refer to data related to hydrocarbon production that 60 volumetric density image of the hydrocarbon reservoir 150 .
can be measured at wells such as historical data including One skilled in the art would recognize that the degree of
well production rates, bottom - hole pressure, phase produc- proximity of placing the muon detector 120 with respect to
tion, injection rate , and water cut . Production data may also the hydrocarbon reservoir may vary depending on factors
include “ i - field data . ” Such “ i - field data” includes well data, such as the sensitivity and the depth of the muon detector
group data , and separator data , all of which can be delivered 65 120 , the density of hydrocarbons existing in the hydrocarbon
in real time to a centralized information technology ( IT ) reservoir, and the density of geological structures surround
system and all of which are accessible to the modeler. ing the hydrocarbon reservoir.
 US 11,125,905 B2
7 8
Image reconstruction of the hydrocarbon reservoir 150 matching tool incorporating EnKF or EnKS is provided so
can be performed by determining the incoming and outgoing that accurate estimates of permeability and porosity field
muon trajectories 140 as the muons pass through the data are available for the history matching process . The
scanned volume ( corresponding to the hydrocarbon reser- history matching process involves several geological real
voir 150 ) and the muon detector 120. The incoming and 5 izations accounting for the uncertainties in the initial model
outgoing muon trajectories 140 are determined using known through a probabilistic approach .
locations of the muon detectors 120 at which the muons The base survey is initially acquired before production is
were incident. Muon detection data acquired by such means begun in the reservoir. Subsequent surveys can be performed
can be inverted to reconstruct a volumetric density image later during production to quantify the fluid displacement in
which shows the spatial location of the hydrocarbon reser- 10 the reservoir . The time -dependent difference ( the difference
voir 150. Without being bound by any theory, the resolution between
of the volumetric density image may depend on the prox differencethebetween monitor survey and the base survey or the
the monitor survey and the subsequent
imity of the hydrocarbon reservoir 150 relative to the earth's survey ) is used for history matching to improve the quality
surface due to the quantity of muons that can be detected .
Also , the resolution of the volumetric density image may 15 ofmodel
the, match between actual productionin and
thus reducing the uncertainty the the reservoir
porosity and
depend on the proximity of the muon detector 120 relative permeability field data .
to the earth's surface due to the quantity of muons that can
be detected . In some embodiments, one may prolong the Ensemble Kalman Filtering
predetermined muon detection period. A time -averaged In some embodiments, to provide accurate estimates of
volumetric density image can be obtained . 20 state variables the initial ensemble is conditioned with
In some embodiments, the muon detector 120 may production data and model dynamics iteratively over time .
include a scintillator. Photons are generated when ionizing This can be achieved by using methodologies such as EnKF .
radiation is deposited in the scintillator. The photons propa- EnKF is an inverse -modelling local optimization tech
gate in the scintillator material while undergoing multiple nique originated from the Kalman Filter ( KF ) , which has
reflections when scattering off inner walls of the scintillator 25 been designed originally for electrical -signal processing.
material. The scintillated photons reach the entrance of an EnKF provides sequential assimilation of both static and
optical detector which converts the photons into an electrical dynamic data into a model . Static and dynamic data are
signal. This process occurs over a measurable time scale of assimilated sequentially from only the current timestep .
a few nanoseconds per meter of an optical path - length . The EnKF is based on a Bayesian framework and Monte Carlo
scintillator may detect muons from all directions , including 30 simulation that stochastically generates reservoir models
muon trajectories 140 where the hydrocarbon reservoir 150 that are integrated over time to estimate probability density
density does not change over time . functions ( or PDFs ) . A previous PDF is updated to a recent
In other mbodiments , the muon detector 120 may PDF incorporating recent data . This updating can be imple
include a drift tube. The drift tube includes a drift region mented independent of a reservoir simulator. The model
filled with gas and a high voltage plane. A muon traverses 35 estimate of EnKF is derived by maximizing analytically the
the drift region and leaves a track of ionization, where posterior probability without any numerical optimization
ionization electrons drift towards one or both ends of the algorithm . Because EnKF is devised accounting for the
drift region due to the high voltage bias provided by the high stochastic nature of the modelling system , the system is
voltage plane. The drifting electrons reach the end to pro- represented by an ensemble of equally and likely to be
duce an electrical signal . Temporal information of the elec- 40 drawn realizations of model estimates . The model covari
trical signals can be used to derive positional information of ance matrix is replaced by a sample covariance matrix ,
the incident muon . which is computed from the ensemble members. EnKF
History Matching requires storing a portion of the covariance matrix that
FIG . 2 shows an example workflow 200 for history describes the model- to - data and data - to -data relationship .
matching of production data and muographic 4D data in 45 Instead of deriving and solving analytical equations of
accordance with an embodiment of the disclosure . The hydrocarbon reservoir behavior, any suitable flow simula
workflow 200 is an iterative process as shown as the dotted tion methodology can be coupled with EnKF to constrain the
arrow that updates over time as new production data , new reservoir models and to predict future reservoir perfor
muographic 4D data , or both , are obtained over time . mance . Gradient computation is not necessary.
The first step in this workflow 200 as indicated by block 50 In some embodiments, EnKF provides a set of geological
210 is the utilization of the previous parameters and certain realizations which are generated to take in account the
state data estimates ( such as static data and dynamic data ) to uncertainties in the geological model . This set of realizations
perform the reservoir simulations to establish a set of is forwarded in time and corrected dynamically as new
reservoir models . The simulations are forwarded in time to production data ( such as true reservoir performance ), new
predict future reservoir performance as indicated by block 55 time-dependent data ( such as muographic 4D data ), or both
220. À monitor survey is computed based on the recently are obtained over time . For each realization , a prediction of
updated reservoir state data at the end of the predicted period reservoir performance and a time -dependent response are
and according to the certain reservoir parameters. The generated. The prediction of reservoir performance and
time-dependent difference (that is , the difference between time- dependent response are used during the EnKF process
the monitor survey and a base survey) is then computed and 60 ing to estimate the covariance matrix. The production data
used in conjunction with the predicted reservoir production and time-dependent data obtained from the field are used to
data in block 230 to update the previous reservoir model correct states of the realizations as they are formed . Conse
estimates , as well as to serve as a starting point for the next quently, forecast values generated by these realizations
iteration of history matching . approximate the observed values . By generating an
Embodiments of the disclosure provide quantitative 65 ensemble of realizations and updating the ensemble over
incorporation of the production data and muographic 4D time , an error estimate of the forecast can be obtained
data in the history matching process . An automated history providing additional information to the decision maker.
 US 11,125,905 B2
9 10
The EnKF uses a sample or ensemble of state vectors as Parameters including porosity, permeability, saturation,
provided in Equation ( 1 ) : and pressure are included in the state vector during process
{ x',i= 1,2 ,... } ing according to embodiments of the disclosure. Permeabil
Equation ( 1 )
ity and porosity correspond to static variables. Pressure and
where N, denotes the number of ensemble members. 5 phase saturation correspond to dynamic variables, that is ,
Assume that a set of an analysis ensemble {xx-19,1 , i= 1 , variables produced by the reservoir simulator in runtime.
2 , ... , Ne } is available at a given time tk -1 ( superscript a The assimilated data include production data obtained from
referring to analysis ). At every forecast step , all ensemble the original model perturbed with Gaussian errors .
members are integrated forward in time with the reservoir Ensemble Kalman Smoothing
model of the next observation (that is , the reservoir model at 10 In some embodiments, to provide accurate estimates of
time t) compute a forecast ensemble { x 1, i = 1 , 2 , ... , state variables the initial ensemble is conditioned with
Ne } ( superscript f referring to forecast) represented by historical data and model dynamics iteratively over time .
Equation (2 ) : This can be achieved by using methodologies such as EnKS .
xx{ i= M k-126 * +1%-A)+nx'+ Equation ( 2 ) EnKS is a modification of EnKF, where static and
where M k - 1 , k is the dynamic model forward operator cor 15 dynamic updatato the
are integrated sequentially and the model is
responding to times tk -1 and tx , and nk , is the ith column updated current timestep conditional to static and
dynamic data from current and previous timesteps. In terms
noise at time tk. The state estimate ( taken as a sample mean of reservoir modelling, EnKS can provide realistic forecasts
of the forecast ensemble ) and the associated covariance
matrix are respectively estimated using Equation (3 ) and 20 because the state vector generally involves static properties
such as permeability and porosity, which seldom change
Equation ( 4 ) :
over time . EnKF would require the simulator to restart at
every single assimilation step . EnKF would also require a
Equation ( 3 ) fall back mechanism for dynamic variables such as pressure
? =?? tim and phase saturation deviating from a predetermined range
25 or value . On the other hand, EnKS makes all changes at the
Equation (4 ) end of the simulation cycle by integrating over the entire
P ( simulation time and undergoing reiteration if necessary . The
i= 1 simulations are submitted in batch and the results are
analyzed after the full simulation batch has concluded .
30 Muography Assisted History Matching
where X is the forecast ensemble average at time to Pris History matching algorithms such as EnKF or EnKS can
the associated covariance matrix of the forecast ensemble at assimilate muographic 4D data even in cases where the
time tz , and superscript T is the transpose of a given matrix . generated volumetric image results in a pixelated image . The
The associated covariance matrix can be written as in muographic 4D data can be used to reduce the uncertainty of
Equation ( 5 ) : 35 the ensemble as the amount of muographic 4D data accu
PA= X / X Equation ( 5 ) mulates over time ( for example, from the muon detectors ) .
where the ith column of X is The muographic 4D data can be continuously, or at least
periodically, fed into the history matching loops of EnKF or
EnKS such that the reservoir model can be updated. Accord
40 ingly , the updated reservoir model provides an accurate
(Ne - 1 )2 ( x* - ). forecast to the decision maker. The accuracy of history
matching techniques relies on the quality of the reservoir
The analysis step is performed to every member of the model and the quality and quantity of the actual data .
FIG . 3 shows a process 300 illustrating the muography
forecast ensemble using the linear KF update Equation (6) : 45 assisted history matching reservoir modelling workflow in
x24.1 = x + * + KxDX -HX ) Equation ( 6 ) accordance with an embodiment of the disclosure . FIG . 3
where K is the ensemble -based approximate Kalman gain at also during
illustrates the computational methodology taking place
history matching of muographic 4D data using EnKF
time t , provided in Equation ( 7 ) :
or EnKS in accordance with an embodiment of the disclo
Rx = Xx (H XDT[(H XO (H2XD7 +Rx]-1 Equation (7) 50 sure .
ya is the observation perturbed with noise sampled from the utilizing In block 310 , a static geological model is constructed by
distribution of the observational error at time tk, Hk is the geophysicalstaticanddata . Actual static data including certain
petrophysical data can be acquired at vari
linear measurement operator at time tze, and Rx is a covari ous observation locations of the reservoir. In some embodi
ance matrix corresponding to the observational error at time 55 ments , the data include data sets of well logs , seismic , core,
tk.
The analysis state and its covariance matrix are then and petrophysical information of the hydrocarbon reservoir.
expressed by Equation ( 8 ) and Equation ( 9 ) , respectively : An initial estimate of several possible realizations of a
reservoir model is constructed . In some embodiments, the
static model can include an ensemble of models to reduce
Ne Equation ( 8 ) 60 uncertainty in the actual data. The uncertainty associated
=
have to with the actual static data can be addressed by creating a set
of equiprobable scenarios. For example , a set of equiprob
Equation ( 9 ) able scenarios can be created by implementing two scenarios
P = Ne (x - 1 )(x - X having different permeability distributions generated by a
65 sequential Gaussian simulation where both distributions
reflect the well data and variograms specified by the geo
modeler. The variability ( corresponding to uncertainty )
 US 11,125,905 B2
11 12
between the two scenarios may increase when the perme- lation time and computing power can be reduced by updat
ability distributions significantly deviate from data points ing the reservoir model when significant changes in density
provided by quantitative data . In other embodiments , the are observed in the reservoir.
static model can include a single model . In block 350 , the updated ensemble run by an iteration in
In block 320 , a dynamic model is constructed by utilizing 5 block 330 undergoes a consistency check . Production and
the single static model or the ensemble of static models pressure responses generated as simulated outputs in block
established in block 310. In addition, the dynamic model is 330 are compared with available historical data and i - field
constructed by incorporating dynamic data such as fluid flow data provided in block 340 to determine a time -dependent
difference . If the time-dependent difference does not deviate
properties. Fluid flow properties may include viscosity and 10 more
density. The dynamic model can include pressure , volume , tion and than a predetermined range or value between produc
and temperature dependencies for the fluid flow properties . the timepressure responses, historical data , and i - field data ,
- dependent difference and the predicted reservoir
The dynamic model can include flow equations and solu model simulated for each ensemble member is incorporated
tions of the flow equations. In some embodiments, the static in the history matching loop by EnKF or EnKS . The
model can be upscaled to reduce the number of cells used in 15 observed data ( reservoir performance and time -dependent
the simulation . As used throughout the disclosure, the term data ) is perturbed with Gaussian errors by a user specified
" cell ” refers to an elementary spatial discretization within a standard deviation . The square of these same values are used
partial differential equation where unknown multivariable as given variances of the observational errors in the com
functions are considered constant. The number of cells (or putation of the Kalman gain matrix in Equation (7 ) , which
grid blocks ) correspond to the complexity of the reservoir 20 are assumed uncorrelated ( that is , covariance matrix R is
model . The number of cells correspond to the complexity of diagonal). The Kalman gain matrix is computed taking as
the linear system to be solved in a simulation model . As used input the predicted reservoir performance, the time-depen
throughout the disclosure, the term “ upscaling” refers to a dent difference for each ensemble member, and the covari
process of averaging a series of cells into a reduced number ance matrix as expressed by Equation (7 ) . The Kalman gain
of cells . The objective is to reduce the overall number of 25 is used to update the state vector by means of Equation ( 6 ) .
cells while preserving the underlying flow characteristics . Still in block 350 , if the time-dependent difference devi
Still in block 320 , the dynamic model can be run reflecting ates more than a predetermined range or value between
a subset of historical production constraints such as produc production and pressure responses, historical data , and
tion data including historical oil or liquid production rates. i - field data ( that is , the predicted reservoir model does not
In block 330, the ensemble corresponding to the dynamic 30 pass
dynamicthe model
consistency
can be check ), input
manually parameters
corrected such thatof the
the
reservoir model realizations constructed in block 320 is simulation can be restarted .
iteratively updated ( or calibrated ) by inputs provided in
blocks 320 and 370. In some embodiments, the updating byInobtaining blocks 360 and 370 , volumetric density data provided
saturation fields generated by muons can be
process ( or history matching) can be executed manually . In 35 incorporated into
other embodiments, history matching can be executed in an improve the qualitytheof the assisted history matching loop to
calibration process .
automatic fashion ( automated history matching ) by incor In block 360 , each muon detection event provides infor
porating certain mathematical algorithms and specialized mation regarding the muon trajectory and the average den
software packages. The ensemble can be calibrated by sity along the trajectory path as shown for example in FIG .
history matching techniques such as an EnKF or EnKS . 40 1. In block 370 , muon detection data obtained in block 360
These techniques can assimilate data while accounting for is used for tomographic inversion to create a time-dependent
measurement and model errors . In some embodiments, volumetric density image of the hydrocarbon reservoir.
block 330 includes forecasting hydrocarbon production Optionally, muon detection period can be adjusted based on
based on the updated ensemble. In some embodiments, the the uncertainty or resolution of the inverted tomographic
updated reservoir model can be imaged as an output. 45 image . The adjustment of the muon detection period may
Still in block 330 , as part of the history matching utilizing alter the waiting time for each ensemble update . For
EnKF or EnKS processing , the dynamic reservoir model example, the detection period can be extended to obtain
realizations initially constructed in block 320 are iterated more muon detections in block 360 before the ensemble
sequentially based on the production data gathered since the updating process in block 330 is initiated in the next
last iteration in block 340. In block 340 , the production data 50 iteration . In the next iteration, a time -averaged volumetric
includes historical data and i - field data . The dynamic reser- density image can be obtained in block 370. In some
voir model realizations initially formed in block 320 are embodiments, such time -averaged volumetric density image
iterated sequentially based on the muographic 4D data can be a pixelated image that still can be assimilated by
gathered since the last iteration in block 370. Muography EnKF or EnKS .
assisted history matching involves a constant input flow of 55 Further modifications and alternative embodiments of
muographic 4D data , which can be averaged over time . The various aspects of the disclosure will be apparent to those
muographic 4D data input provided in block 370 includes a skilled in the art in view of this description . Accordingly, this
volumetric image of fluid saturations with a varying range of description is to be construed as illustrative only and is for
uncertainty based on the number of muon detection events the purpose of teaching those skilled in the art the general
over the predetermined timeframe. Both the muographic 4D 60 manner of carrying out the embodiments described in the
data provided in block 370 and historical data and i - field disclosure . It is to be understood that the forms shown and
data provided in block 340 are available in real time or at described in the disclosure are to be taken as examples of
least periodically. This allows the ensemble -based reservoir embodiments. Elements and materials may be substituted
model to continuously or at least periodically update itself as for those illustrated and described in the disclosure , parts
new inputs of muographic 4D data , historical data , and 65 and processes may be reversed or omitted , and certain
i - field data are introduced over time such that the reservoir features may be utilized independently, all as would be
model can assimilate the data . In some embodiments, simu- apparent to one skilled in the art after having the benefit of
 US 11,125,905 B2
9

13 14
this description. Changes may be made in the elements constructing a subsequent volumetric density image of the
described in the disclosure without departing from the spirit hydrocarbon reservoir by processing the subsequent
and scope of the disclosure as described in the following muon detection data via tomography ; and
claims . Headings used described in the disclosure are for updating the reservoir model utilizing the volumetric
organizational purposes only and are not meant to be used to 5 density image , the updating comprising adjusting the
limit the scope of the description . reservoir model utilizing the subsequent volumetric
density image .
What is claimed is : 11. The method of claim 10 , further comprising the step
1. A method for surveying a hydrocarbon reservoir uti- of:
lizing a reservoir model, the method comprising the steps of: 10 forecasting a hydrocarbon production from the hydrocar
establishing an ensemble of models reflecting attributes of bon reservoir with the reservoir model in the updating
the hydrocarbon reservoir based on the reservoir model step .
in its present state; 12. The method of claim 10 , further comprising the step
obtaining muon detection data collected over a detection of:
period from a muon detector at a subterranean location 15 imaging the reservoir model in the updating step .
proximate the hydrocarbon reservoir ; 13. The method of claim 10 , where the muon detector
constructing, via tomography processing of the muon includes one selected from the group consisting of: a scin
detection data, a volumetric density image of the tillator, a drift tube, and combinations of the same .
hydrocarbon reservoir ; 14. The method of claim 10 , further comprising the step
obtaining subsequent muon detection data collected after 20 of:
a time increment from a muon detector at a subterra- altering the detection period based on a resolution of the
nean location proximate the hydrocarbon reservoir; volumetric density image.
constructing, via tomography processing of the subse- 15. A method for surveying a hydrocarbon reservoir
quent muon detection data , a subsequent volumetric utilizing a reservoir model , the method comprising the steps
density image of the hydrocarbon reservoir ; and 25 of:
updating the reservoir model by utilizing a volumetric establishing an ensemble of models reflecting attributes of
density image of the hydrocarbon reservoir, the updat- the hydrocarbon reservoir based on the reservoir model
ing comprising adjusting the reservoir model utilizing in its present state ;
the subsequent volumetric density image . positioning a muon detector at a subterranean location
2. The method of claim 1 , where the volumetric density 30 proximate the hydrocarbon reservoir;
image is constructed via muon tomography. collecting muon detection data from the muon detector
3. The method of claim 1 , further comprising the step of: over a detection period ;
forecasting a hydrocarbon production from the hydrocar- constructing a volumetric density image of the hydrocar
bon reservoir with the reservoir model in the updating bon reservoir by processing the muon detection data via
step . 35 tomography;
4. The method of claim 1 , further comprising the step of: collecting a subsequent muon detection data from the
imaging the reservoir model in the updating step . muon detector after a time increment;
5. The method of claim 1 , where the updating step constructing a subsequent volumetric density image of the
includes performing one selected from the group consisting hydrocarbon reservoir by processing the subsequent
of: Ensemble Kalman Filtering , Ensemble Kalman Smooth- 40 muon detection data via tomography; and
ing , and combinations of the same. updating the reservoir model by utilizing the volumetric
6. The method of claim 1 , where the attributes of the density image, the updating comprising adjusting the
hydrocarbon reservoir include one selected from the group reservoir model utilizing the subsequent volumetric
consisting of: static data , dynamic data, and combinations of density image.
the same. 45 16. The method of claim 15 , further comprising the step
7. The method of claim 1 , where the updating step of:
includes providing production data to update the reservoir forecasting a hydrocarbon production from the hydrocar
model . bon reservoir with the reservoir model in the updating
8. The method of claim 7 , further comprising the step of: step .
comparing the production data and the reservoir model in 50 17. The method of claim 15 , further comprising the step
the updating step . of:
9. The method of claim 8 , where the comparing step imaging the reservoir model in the updating step .
includes correcting input parameters of the reservoir model 18. The method of claim 15 , where the updating step
when a difference between the production data and the includes performing one selected from the group consisting
reservoir model deviates from a predetermined value . 55 of: Ensemble Kalman Filtering, Ensemble Kalman Smooth
10. A method for surveying a hydrocarbon reservoir ing , and combinations of the same .
utilizing a reservoir model , the method comprising the steps 19. The method of claim 15 , where the attributes of the
of: hydrocarbon reservoir include one selected from the group
positioning a muon detector at a subterranean location consisting of: static data , dynamic data , and combinations of
proximate the hydrocarbon reservoir; 60 the same.
collecting muon detection data from the muon detector 20. The method of claim 15 , where the updating step
over a detection period ; and includes providing production data to update the reservoir
constructing a volumetric density image of the hydrocar- model.
bon reservoir by processing the muon detection data via 21. The method of claim 20 , further comprising the step
tomography, 65 of:
collecting subsequent muon detection data from the muon comparing the production data and the reservoir model in
detector after a time increment; the updating step .
 US 11,125,905 B2
15 16
22. The method of claim 21 , where the comparing step
includes correcting input parameters of the reservoir model
when a difference between the production data and the
reservoir model deviates from a predetermined value .
23. The method of claim 15 , further comprising the step 5
of:
altering the detection period based on a resolution of the
volumetric density image.
*