Adaptive optimal experimental designand inversion of a coupled uid owand geophysical imaging model forreservoir monitoringbyJennifer FohringB.Sc., The University of British Columbia, 2011A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinThe Faculty of Graduate and Postdoctoral Studies(Geophysics)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)May 2016c Jennifer Fohring 2016AbstractImaging and prediction of uid ow within the subsurface provides information crucialto decision making processes in elds such as groundwater management and enhancedoil recovery. The ow of a uid through a reservoir depends primarily on the perme-ability of the subsurface rock; a quantity that is often unknown throughout the entiredomain of the reservoir. One method for predicting ow is to estimate the permeabil-ity of the reservoir and simulate ow through a mathematical subsurface ow model.Given the model, ow data can be inverted to estimate the permeability. However,this inversion approach can lead to inaccurate results due to the sparse sampling ofow data, and thus inaccurate predictions.To acquire a higher sampling of data, geophysical survey techniques are appliedin order to eciently collect a higher density of data sampled at the surface. Thesedata are sensitive to changes to the geophysical properties of the reservoir due toow. Inversion of geophysical data then provides images of changes to the geophysicalproperties of the reservoir. In order to estimate the ow parameters using geophysicaldata, the two mathematical models require coupling.The thesis therefore proposes two approaches to improve the imaging and predic-tion of ow. First, a novel coupled inverse problem for estimating the uid velocityeld and the initial geophysical property model from geophysical data is developed.Second, a new method of optimally designing the geophysical survey for the coupledinverse problem is developed. The new adaptive design approach builds on traditionalA-Optimal design methods such that historic data are included in the design algo-iiAbstractrithm. This produces designs that adapt with ow in the subsurface and reduce thecollection of unnecessary data. Both the coupled inverse problem and adaptive sur-vey design method are demonstrated using a seismic tomography geophysical surveyand a tracer advection uid ow model. Numerical examples show that the coupledapproach yields an improved ow estimate as well as improved image quality, whilethe adaptive optimal designs provide sucient geophysical data.iiiPrefaceThis thesis contains original research conducted while studying at the University ofBritish Columbia, resulting in two publications and one expanded conference pro-ceeding.The idea of the coupled inverse problem presented in Chapter 4 came originallyfrom conversations with Dr. Eldad Haber. The subsequent derivations, code imple-mentation, numerical tests, and manuscript preparation were carried out in collabo-ration with Dr. Lars Ruthotto, a post doctoral researcher with Dr. Eldad Haber atthe time. The work was published in Fohring, J., Haber, E., and Ruthotto, L. (2014).Geophysical imaging of uid ow in porous media. SIAM Journal on Scientic Com-puting, 36(5):218{236), and various parts of were adapted from an SEG conferenceproceeding ( Fohring, J., Ruthotto, L., and Haber, E. (2013). Geophysical Imaging,Reservoir History Matching and Forecasting. In 2013 SEG Annual Meeting, Houston.Society of Exploration Geophysicists).The idea of adaptive experimental design also resulted from conversations with Dr.Eldad Haber. The development of this idea is presented in Chapter 5. The subsequentderivations, code implementation, numerical tests, and manuscript preparation foradaptive optimal design were carried out by myself with contributions, input, andadvice from Dr. Haber. The bulk of the text in Chapter 5 is adapted from thepaper (Fohring, J. and Haber, E. (2016). Adaptive A-optimal experimental design forlinear dynamical systems. SIAM Journal on Uncertainty Quantication, xx:1{19),which is still at the time of writing, going through the review process.ivTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Reservoir characterization . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Geophysics for reservoir characterization . . . . . . . . . . . . . . . . . 41.3 Optimal survey design . . . . . . . . . . . . . . . . . . . . . . . . . . . 101.4 Thesis objective and outline . . . . . . . . . . . . . . . . . . . . . . . . 122 A model reservoir monitoring experiment . . . . . . . . . . . . . . . . 152.1 Geophysics: seismic tomography . . . . . . . . . . . . . . . . . . . . . 172.2 Dynamics: tracer advection . . . . . . . . . . . . . . . . . . . . . . . . 182.3 Physical parameter relations . . . . . . . . . . . . . . . . . . . . . . . . 202.4 Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 222.4.1 Discretization of the tomography equation . . . . . . . . . . . . 22vTable of Contents2.4.2 Discretization of the mass balance . . . . . . . . . . . . . . . . 232.4.3 Discretization of the ux-balance conservation . . . . . . . . . 242.4.4 A Particle-In-Cell method for ow simulation . . . . . . . . . . 252.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 283 Error characterization and model estimation . . . . . . . . . . . . . . 303.1 Linear inversion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 303.2 Coupled inversion and error characterization . . . . . . . . . . . . . . . 343.2.1 Exact dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . 343.2.2 Inexact dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . 363.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 374 Coupled inversion for velocity and initial slowness . . . . . . . . . . 394.1 Flow constrained geophysical imaging . . . . . . . . . . . . . . . . . . 404.2 Discretization of regularization functionals . . . . . . . . . . . . . . . . 424.3 Numerical optimization . . . . . . . . . . . . . . . . . . . . . . . . . . 444.4 Initialization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 494.5 Experimental setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 504.6 Sequential reconstruction as initialization . . . . . . . . . . . . . . . . 504.7 Joint reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 544.8 Forecasting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 554.9 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 565 Dynamic optimal experimental design . . . . . . . . . . . . . . . . . . 575.1 Adaptive experimental design . . . . . . . . . . . . . . . . . . . . . . . 585.2 Design for dynamical systems . . . . . . . . . . . . . . . . . . . . . . . 645.2.1 Exact dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . 645.2.2 Inexact dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . 65viTable of Contents5.3 Numerical optimization . . . . . . . . . . . . . . . . . . . . . . . . . . 665.3.1 Exact dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . 665.3.2 Inexact dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . 695.4 Numerical examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . 705.4.1 Regularization . . . . . . . . . . . . . . . . . . . . . . . . . . . 725.4.2 Monitor function . . . . . . . . . . . . . . . . . . . . . . . . . . 745.5 Design results: exact dynamics . . . . . . . . . . . . . . . . . . . . . . 745.6 Results: inexact dynamics . . . . . . . . . . . . . . . . . . . . . . . . . 785.7 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 786 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89viiList of Figures1.1 Model of a reservoir with minimal wells. Estimation of the permeability(k(~x)) using only pressure (p) and hydraulic head (z) data of a uidmoving with velocity (~u), collected from sparsely distributed wells willgenerate highly inaccurate results. . . . . . . . . . . . . . . . . . . . . 31.2 Model of a reservoir with geophysical survey measurement locations de-picted in orange. Changes to the geophysical model parameters (m(~x))can be estimated through inversion of the densely sampled geophysicaldata (di). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.1 Tomography setting: sources (si) and receivers (rj) are placed in t-wo boreholes and the travel times of an acoustic wave traveling alongadjoining ray paths (i;j) are measured. . . . . . . . . . . . . . . . . . 172.2 Tracer advection: Tracer of concentration (c) is advected in a uidof density () with velocity eld (~u). The hydraulic conductivity (K)determines the velocity eld given a pressure gradient within the reser-voir. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 202.3 A grid cell with length h1and height h2, slowness value (mj) locatedat the cell center, and components of the uid velocity owing in (~u =[u+1; u+2]) and out (~u = [u1; u2]) on the cell edges. . . . . . . . . . . . 24viiiList of Figures2.4 Particle-In-Cell method for the advection and mass-conservation equa-tion. A particle located in the cell-centered point~xAon the regularmesh is pushed forward with velocity~u along the path~ut, to thenon-grid point~xB. The value (mass) associated with the particle isthen transferred to the cell centered grid points adjacent to~xBusinginterpolation weights. . . . . . . . . . . . . . . . . . . . . . . . . . . . 274.1 Visualizations of the ground true and reference hydraulic conductivitymodels and associated ow elds. The true and reference models arediscretized on a 100200 cell grid on a domain of 100m depth by 200macross. A static source, at 55m depth on the left border, and a sink,at 65m depth on the right border of the domain, are indicated by red.The reference model is a simplied model which interpolates boreholedata between the two wells. . . . . . . . . . . . . . . . . . . . . . . . . 514.2 Visualization of ground truth data for the domain with dimensions100m (y direction) depth and 200m (x direction) across (rst colum-n), individual reconstruction (second column) and results of the jointreconstruction approach. The rst row visualizes the initial slowness(m0). The second shows a quiver plot of the recovered velocity eld(~u). The components of the velocity eld [ux; uy] are pictured in thethird and fourth row and the last row shows the curl of the velocityeld (r~u). Note that there is almost no ow in the y direction, andthat at the boundaries of the reservoir there are non-zero values of thecurl of~u . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52ixList of Figures4.3 The evolution of the tracer is pictured for the ground-truth (left col-umn), individual reconstruction used for initialization (middle column),and the proposed joint reconstruction. The tracer evolution is simu-lated for 40 days with the respective initial slowness and ow eld.The front of the ground truth tracer is visualized by dashed lines forthe rst 3 rows. Note that the joint method results match the groundtruth tracer front to a greater degree than the predictions from theinitialization estimates. . . . . . . . . . . . . . . . . . . . . . . . . . . 535.1 Initial experiment setup. Tomography sources are pictured in greenon the left and receivers in blue on the right. The initial setup coversthe entire ow domain with 20 sources and 30 receivers. There is onesource of uid, pumped at a constant rate from the top green point,and a sink at the bottom (red) of the domain. Flow is in the downwarddirection. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 715.2 Plot of the relative error (kmtrue1bm1k=(kmtrue1k)) versus the regular-ization parameter for both the linear regularization and smoothedtotal variation. The TV regularization gives a better estimate of themodel for the case of a sharp target. . . . . . . . . . . . . . . . . . . . 735.3 Exact dynamics: plots of the adapted mean squared error amse vs.the number of nonzero weights (left column), and the weights used toconduct experiments 1, 3, and 5, at times t0; t2; t4. . . . . . . . . . . 755.4 Exact dynamics: relative error plots. Model error is compared for eachexperiment for both the linear and TV regularization . . . . . . . . . . 77xList of Figures5.5 Exact dynamics: optimal designs and recovered models for 9 experi-ments. The top row shows the rays and the recovered model, with thenumber of rays (#d = 297) given below, followed by the true modelsin the second row, the models recovered using total variation in thethird row, and nally the models recovered using the linear gradientregularization. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 795.6 Inexact dynamics: relative error plots. Model error is compared foreach experiment for both the linear and TV regularization . . . . . . . 805.7 Inexact dynamics: optimal designs and recovered models for 9 experi-ments. The top row shows the rays and the recovered model, with thenumber of rays given below. The second row pictures true models, thethird shows the models recovered using total variation, and nally thebottom shows models recovered using the linear gradient regularization. 82xiAcknowledgementsThis thesis would never have been possible without the guidance, assistance andsupport of a large group individuals.First and foremost, I would like to express my deepest gratitude to my supervi-sor Eldad Haber for all of the guidance, advice, and arguments he has provided. Ihave learned a great deal working with you over the past 6 years, not only aboutcomputational geophysics, but also about myself. Toda raba.During my PhD I had the privilege of working at IBM with Lior Horesh. ThankLior for your kindness, encouragement, and for providing me with so many dierentexplanations and perspectives on inverse problems and optimization. I am trulygrateful for the time spent with you and your lovely family in New York.And to Lars Ruthotto, the paper wouldn't have happened without you. I learnedso much working with you. Thank you.Many thanks to all of my oce mates and colleagues; particularly Klara and Luzfor your much appreciated reading and editing of my thesis, and for lending an ear tomany questions, concerns, and complaints! And also to Gudni, Sarah, Justin, Kris,Mike, Mike, Dikun, Rowan, Seogi,and Lindsey, for help, encouragement, editing, andfun!A very special thanks to my dear friends Marianne and JF: for campus lunches,hikes, dogs, skiing, beer, wine, rowing, science, and games.And nally, thank you to my best friend and partner Dave Marchant. For every-thing.xiiAcknowledgementsThe funding of this work was generously provided from a variety of sources.Thanks to University of British Columbia, the Department of Earth and Ocean Sci-ences, the Society of Exploration Geophysics, the Natural Sciences and EngineeringResearch Council of Canada, and MITACS for their support.xiiiTo my parents.xivChapter 1IntroductionSubsurface reservoirs are dened as regions of uid saturated porous rock and arecommonly depleted through uid extraction. Maintenance of groundwater aquifers(a reservoir of water) or oil producing reservoirs as they are depleted, often requiresecient methods of monitoring uid ow within the subsurface. The ability to trackand predict the ow is important for both environmental and economical reservoirmanagement. For example, in an oil recovery scenario, where a uid such as water isinjected into a depleted oil reservoir to enhance recovery, it is crucial to know wherethe oil and water are at a given time. Monitoring this process provides informationcrucial to future uid management decisions. For example, decisions such as fromwhich wells to extract, which wells to shut in, and from which wells to inject canaect production rates, uid loss, and overall eciency.One particular aspect of reservoir management programs involves predicting orforecasting uid ow within the subsurface. To do so, ow is simulated by building amathematical model of the subsurface, often referred to as reservoir characterization.1.1 Reservoir characterizationIn most cases, reservoirs consist of isolated regions of porous rock surrounded byregions of non-porous rock. See Figure 1.1 below for a diagram. Thus, the dynamicsof a reservoir are commonly governed by mathematical models describing uid owin a porous media, which come in many forms with varying degrees of complexity.11.1. Reservoir characterizationSome of the more common models include: single phase ow, two-phase ow, blackoil models, advection-diusion models, and reactive contaminant transport models.For example, the model describing single phase ow in a porous media as presentedin (Chen et al., 2006), given the spatial and temporal variables~x and t, is as follows@(ffi)@t= r (~u) + q; (1.1a)~u =k(rp grz): (1.1b)Equation 1.1a is derived by tracking mass conservation of a uid with viscosity owing through a rectangular three dimensional dierential volume, see (Chen et al.,2006) for the complete derivation. The porosity of the rock (the fraction of a rep-resentative elementary volume available for uid) is denoted by ffi, and the densityof the uid by (t;~x). The velocity of the uid~u(~x) is dened as the supercialDarcy velocity, and described by Equation 1.1b, the conservation of momentum e-quation also known as Darcy's law. The ow in this model is generated by a pressuregradient, rp(t;~x), resulting from the sources and sinks, q, and from a gravitationalgradient, grz, where rz is the dierence in uid heights, and g is the magnitude ofthe gravitational acceleration.If there is no compressibility in the rock then@(ffi)@t! 0, such that no uid is storedwithin the rock pore space. Additionally, the velocity depends on the permeabilitytensor k of the porous rock and the uid viscosity . The permeability is an averageproperty of the reservoir that measures the ability of a porous media to transmit uid.In most cases, k is a tensor, unless it is assumed that k(~x)11= k(~x)22= k(~x)33= k(~x).In this case the porous media is known as isotropic, where the ow is the same in alldirections. It is otherwise anisotropic.In order to simulate ow using a reservoir model and solve for the pressure p,the permeability k(~x), often referred to as a reservoir parameter function, must be21.1. Reservoir characterizationFigure 1.1: Model of a reservoir with minimal wells. Estimation of the permeability(k(~x)) using only pressure (p) and hydraulic head (z) data of a uid moving withvelocity (~u), collected from sparsely distributed wells will generate highly inaccurateresults.known. However, this function is often unknown and thus must be estimated fromborehole rock samples. In most cases boreholes are expensive to drill, and are thussparsely distributed throughout the reservoir. This leads to inaccurate estimates ofk(~x) for the entire region, as values must be extrapolated between boreholes overlarge distances (Roubinet et al., 2013).One approach to alleviating this problem is to use historic ow data, measurementsof the pressure p and hydraulic head z, to estimate the unknown parameter functionk(~x). This process is alternately known as parameter estimation, inversion, and inthe oil and gas industry is referred to as history matching (Oliver and Chen, 2010;Oliver et al., 2008). The basic idea of history matching is to estimate the reservoirparameters such that simulated ow data ts measured data.Technically, the inversion process is rather challenging as it involves the discretiza-tion of a system of partial dierential equations, the solution of the forward problem(ow simulation), the computation of the gradients of the simulator with respect tothe parameters, and the solution of an optimization problem. An excellent review of31.2. Geophysics for reservoir characterizationthe process can be found in (Oliver and Chen, 2010).Although it is possible to estimate the reservoir parameters through inversion andpredict ow data, these estimates can be a highly inaccurate representation of thereservoir. This inaccuracy stems from the large null space of the inverse problemassociated with the highly sparse spatial sampling of reservoir ow data.To overcome the null space issue subspace techniques have been applied, wherethe reservoir parameters are restricted to \live" in a small subspace spanned by a(relatively) small number of vectors (Abacioglu et al., 2001; Gerritsen and Durlof-sky, 2005; Oliver and Chen, 2010; Sarma et al., 2013). This approach decreases thevariance in the recovered estimates of k(~x) by increasing the bias of the estimatedreservoir parameters (Tenorio, 2001).A dierent approach to reduce the variance of the recovered reservoir parametermodels is to simply add additional ow data and further sample the reservoir. Whilethis is clearly one of the better ways to reduce the uncertainty in the estimates, itis not practical as it is unlikely that many more wells will be drilled just to improvethe simulation capability. However, assuming that the geophysical rock properties ofthe reservoir (density or electrical conductivity for example) will change with uidow, implies that geophysical survey techniques can potentially be used to map thesechanges. This makes geophysics for reservoir characterization a practical method toincrease data collection coverage without having to drill additional wells.1.2 Geophysics for reservoir characterizationGeophysical survey techniques measure variation in the geophysical properties of thesubsurface by systematically collecting geophysical data over a large spatial surfacegrid and within existing boreholes. Geophysical data are not only valuable on theirown for detecting variation in the subsurface, but can also be inverted to obtain an41.2. Geophysics for reservoir characterizationimage of the underlying geophysical property distribution. Thus, geophysical data cancomplement reservoir ow data by providing a denser grid of measurements than fromsubsurface ow data obtained from point locations at well sites alone (Hubbard andRubin, 2000). For example, Figure 1.2 provides an example of a geophysical surveywith data, di, collected at the surface (orange dots). Sensors could additionallybe placed in the existing boreholes to further increase data numbers and depth ofinvestigation.Figure 1.2: Model of a reservoir with geophysical survey measurement locations de-picted in orange. Changes to the geophysical model parameters (m(~x)) can be esti-mated through inversion of the densely sampled geophysical data (di).The spatial resolution and depth of investigation of geophysical data depends pre-dominantly on the choice of technique. Applied geophysical surveys are a widely usedmethod of imaging the subsurface in resource exploration. Therefore, there are manywell developed survey types available for imaging subsurface ow. In particular, timelapse seismic surveys, ground penetrating radar, time lapse gravity measurements,and electromagnetic imaging techniques have been used for oil reservoir and aquifercharacterization for some time (Alumbaugh and Morrison, 1995; Archie, 1942; Hub-bard and Rubin, 2000; Lumley, 2001a; Vasco et al., 2004; Vesnaver et al., 2003).51.2. Geophysics for reservoir characterizationAs an example, gravity methods, which are sensitive to shallow changes to thedensity of the subsurface rock, are a common choice in hydrology applications. Wateraquifers tend to be small and close to the surface and thus gravity methods candetect small variations in density as water is depleted, see for example (Alumbaughand Morrison, 1995; Slater et al., 2000; Wilt et al., 1995). In this case measuredgeophysical data consist of small changes to the local gravitational eld above theaquifer.Oil reservoirs in comparison, are often large, at, not necessarily close to the sur-face, and commonly discovered using seismic techniques. Thus seismic surveys arefrequently the rst choice for monitoring programs (Lumley, 2001b). Additionally,electromagnetic methods are attractive for monitoring injected uids because the elec-trical conductivity of the injected uid can be controlled. For example, the injecteduid can be selected such that it's conductivity maximizes the conductivity contrastbetween the existing media and the introduced uid. This then generates a greatersignal in the measured electric and magnetic potential data.Although geophysical surveys are an excellent method for imaging changes in thegeophysical properties of the subsurface, inferring reservoir parameter functions fromgeophysical data or recovered images of the geophysical property models is not sostraightforward. The use of geophysics for inferring reservoir parameters, sometimesknown as geophysical history matching, rests on the assumption that a mathemati-cal relationship (petrophysical relationship) between the geophysical properties andreservoir parameter function(s) exists and is known. Historically, much work has beendone to determine these relationships (Archie, 1942; Mavko et al., 2009). Althoughthese relations tend to be empirical and site specic, they have been successful inprogressing geophysical history matching work-ows.A second point of concern in the geophysical history matching process is that61.2. Geophysics for reservoir characterizationmost often the mathematical geophysical model and the reservoir uid model aredecoupled. Thus, the inversion process is decoupled. This decoupled inversion processhas a major shortcoming; since geophysical imaging is almost always ill-posed, theaccuracy of the imaging is typically low. Estimating the reservoir model parametersfrom the geophysical images alone, or in combination with low resolution modelsobtained from inverting uid ow data, can yield inaccurate and biased estimates.To address this shortcoming, (Hubbard and Rubin, 2000) propose several param-eter estimation approaches using statistical geophysical-hydrological techniques tocharacterize an aquifer with geophysical and hydrological data. One particular tech-nique takes a Bayesian iterative approach where the ow log-permeability (a com-monly measured reservoir parameter) is rst estimated through inversion of sparsehydrological data. This estimate is considered to be a random variable with a proba-bility distribution, and is then used as a prior distribution for the Bayesian inversionof geophysical data. Geophysical data along with the known parameter relationshipsare then inverted to recover an update to the log permeability while including theprior permeability estimate. Alternately, several geostatistical approaches have beenproposed to match structures within images of diering parameters to geologic datafrom core samples, or maximize correlations. These approaches suer the same pa-rameter function requirements and lower resolution biased estimates (Cassiani et al.,1998; Doyen, 1988; Pollock and Cirpka, 2012).In oil reservoir monitoring practice, seismic history matching is the predominan-t method (Abul, 2010; Emerick and Reynolds, 2012; Gosselin et al., 2003; Lumley,2001b; Mezghani et al., 2013; Nenna et al., 2011; Oliver and Chen, 2010; Sarma et al.,2013; Trani, 2012; Vasco et al., 2004). Similar to the geophysical-hydrological tech-niques, seismic history matching involves estimating geophysical properties throughinversion of geophysical data, mapping the estimated geophysical property model toa reservoir parameter model through a physical property relation, and nally simu-71.2. Geophysics for reservoir characterizationlating uid ow data and matching it with observed ow data collected. If the datamatch is poor, either the reservoir parameter model is adjusted, or new geophysicalinversions are carried out and the process is repeated. If the simulations using thereservoir model predict the historic ow data then forecasts can be generated withsome certainty. In general, the seismic history matching process suers as it requiresparameter mappings and the mathematical models are decoupled in the inversions.The concept of coupling the dynamic reservoir model with a geophysical modelin order to gain a better estimate of the unknown reservoir parameters is not new.For example, the Kalman lter was developed in the early 1960's to incorporate newobservation data into the estimation of the position of a moving object and becamea common method for applications in guidance, navigation and control of vehicles,particularly aircraft and spacecraft.Classical Kalman ltering is a Bayesian inversion method that iteratively uses atime-series of noisy observed measurements to estimate unknown model parameters.If the noise is assumed to be characterized by a Gaussian distribution then the lteryields the exact conditional probability following Bayes rule (recursive Bayesian esti-mation), which can also be interpreted as the solution to a linear regularized inverseproblem (Kalman, 1960). In the context of geophysical history matching, Kalmanltering is derived from an inverse problem that simultaneously minimizes the geo-physical measurement error and the error in the ow model simulation.Kalman ltering, and ensemble Kalman ltering (EnKF) (Evensen, 1994), haverecently been proposed as alternative approaches to seismic history matching as theycouple the reservoir model and geophysical model (Abul, 2010; Nvdal et al., 2002;Nenna et al., 2011; Vauhkonen et al., 1998). The Kalman lter (also known as theKalman gain matrix) couples the geophysics, the dynamics, and the prior error co-variance matrices. However, its use for geophysical imaging of ow requires a phys-ical parameter relationship. Additionally, the large dense covariance matrix which81.2. Geophysics for reservoir characterizationcharacterizes the prior distribution must be known, and is propagated in time. Forproblems such as seismic imaging, where the data set and the parameter space arelarge, Kalman ltering becomes computationally infeasible as a result of the storageand computation with the dense covariance matrices. To alleviate this, stochasticmethods such as the EnKF were developed to approximate the covariance matrices.The EnKF method is a Monte Carlo type sequential Bayesian inversion methodrst proposed by Evensen (1994) to reduce computations and memory. The EnKFmethod introduces ensembles of state (parameter models) and measurement vectors(geophysical data for example) sampled from Gaussian distributions such that theprior state covariance matrix can be approximated by computing the covariance ofthe ensemble set. This reduces the inverse problem to the dimension of the ensem-ble. Abul (2010) found that there is promise in EnKF methods for reservoir historymatching, however, it was limited by uncertainty in the initial condition and errorestimates, and by the computation required for integration of large amounts of da-ta and the non-linearity introduced through the relationship between reservoir andseismic parameters.In a similar avour to Kalman ltering, there has been some recent researchto couple non-linear geophysical imaging with ow models in the hydrogeophysicscommunity through regularization of the geophysical inverse problem with the owmodels, with promising results (Cockett and Haber, 2016; Hoversten et al., 2006;Steklova and Haber, 2015). These approaches are dierent in that they were developedfor specic ow models, require Archie's law, are not solved iteratively, and havespecically tailored regularization terms. Additionally, the inversions are carried outfor all times, therefore incorporating all historic data in one inversion.In light of current reservoir monitoring practice, there is room to improve up-on coupled geophysical imaging and subsurface ow inversion methods for reservoirmonitoring.91.3. Optimal survey designIn addition to coupling geophysical survey methods to estimate ow model param-eters, there is reason to consider optimizing the geophysical survey design to adaptto the ow dynamics, such that the best estimates of the reservoir parameters areobtained from a coupled inversion.1.3 Optimal survey designWhile the coupled ow and geophysical inverse problem can potentially provide betterestimates of reservoir parameter functions through inversion of large geophysical datasets, the geophysical data collection can also be designed in an optimal way. Forexample, if one is interested in monitoring the progression of a uid front as it movesthrough a reservoir, collecting geophysical data far downstream may not contributeany additional information to the inversion, and thus those measurements could beexcluded. In the context of reservoir monitoring, an optimal experimental designmight amount to adding or eliminating a well location as demonstrated by Gharamtiet al. (2015), or reducing the number of surface measurements from a maximumallowable set.In many realistic scenarios, the ow dynamics are suciently slow to allow forimproved designs after a rst experiment has been conducted. This is particularly truewhen monitoring ow in reservoirs where uid velocities can be so slow that changesoccur over months or even years (Vasco et al., 2004). Additionally, the availablehistoric experimental data and inversion results will contain valuable informationabout the subsurface. This information should thus be included when determiningan optimal survey design for a future survey.The mathematical method of determining an experimental design is known asoptimal experimental design (OED), and has existed for some time for both the well-posed and ill-posed parameter estimation cases (Aggarwal et al., 2015; Ajo-Franklin,101.3. Optimal survey design2009; Atkinson and Donev, 1992; Bardow, 2008; Chaloner and Verdinelli, 1995; Curtis,1999; Fedorov, 1972; Haber et al., 2008, 2010, 2011; Huan and Marzouk, 2013; Lucero,2013; Pukelsheim and Studden, 1993).Typically, optimal experimental designs are computed by minimizing the error inthe estimated parameter model according to some specied criteria. In this case themodel estimate is assumed to be the mean of a distribution with a covariance matrixcharacterizing the error. The type of design method depends on the optimality criteriaapplied to the covariance matrix. There are several standard optimality criteria; forexample, A-optimal design minimizes the trace of the covariance matrix or averagevariance, D-optimality minimizes the determinant of the covariance matrix, and E-optimality minimizes the eigenvalues of the inverse covariance matrix.There are some optimal design methods which consider dynamic systems, yetthere is little in the way of optimal design algorithms for ill-posed inverse problemswhere the model is governed by a dynamical system. Although there has been workon sequential design for the well-posed case, no work is known for the ill-posed case,(see the following for examples: (Cherno, 1972; Khodja et al., 2010; Papalambrosand Wilde, 2000; Wilkinson et al., 2015) and references within). In some cases,static design methods can be applied to time dependent partial dierential equations(PDEs) (Alexanderian et al., 2014; Haber et al., 2011), this approach results in ana-priori estimate of the designs for \all times" which depend only on informationprovided in the form of a regularization (or prior) to the ill-posed problem. Forexample, (Alexanderian et al., 2014) compute optimal designs for a tracer advection-diusion reservoir model. This example however, did not include any geophysicalimaging, nor was any historic data included in the experimental design computation.There has been some work applying optimal design methods to the Kalman l-ter (Gharamti et al., 2015; Sagnol and Harman, 2014); for example, Sagnol and Har-man (2014) propose a D-optimal design method for the Kalman lter, which results111.4. Thesis objective and outlinein designs dependent on the posterior state estimation covariance matrix. Alter-nately, Gharamti et al. (2015) present an optimal method for selecting aquifer welllocations to track the ow of a contaminant plume. In both cases, the prior covariancematrix to be minimized includes the transition matrix (a discretization of the dynam-ics), but the designs still do not depend on historic data, and thus these approachesare no dierent then the application of static methods.In fact, as far as can be discerned from the literature, apart from sequentialdesign approaches (Cherno, 1972; Khodja et al., 2010), which use the posterior modelestimate and its posterior covariance as the prior for the design at the next time step,current design techniques for linear inverse problems do not utilize information thatis collected at early times in order to design experiments at later times. This is aconsequence of optimizing over the posterior variance in the model estimates.For any linear ill-posed inverse problem, static or dynamic, the posterior covari-ance matrix will only depend on the physics of the problem, and the linear regulariza-tion operator (or prior distribution covariance matrix for the Bayesian perspective),regardless of how the inverse problem is formulated. Thus, in order to include historicdata to a standard optimal design problem, the posterior covariance matrix must bereconsidered.1.4 Thesis objective and outlineThe main objective of the research presented in this thesis is to improve reservoir mon-itoring and forecasting by coupling geophysical survey techniques with a dynamic uidow reservoir model and optimally designing the data collection in the geophysicalsurvey.Coupling the geophysics and ow model in one inverse problem allows for theestimation of uid parameter functions such as the permeability k(~x) from geophysical121.4. Thesis objective and outlinedata. To address the coupling problem, careful consideration of the choice of thereservoir uid ow model and its discretization is taken such that knowledge of thephysical parameter relationship is not necessary. Additionally, error in the reservoirmodel is assumed to be zero. This is contrary to the usual assumptions. In this way,the geophysical inversion problem will be reformulated as an ill-posed inverse problemwith appropriate regularization terms that is constrained by the uid ow model.Since the geophysical survey provides a much denser sampling of data, the uncertaintyin the recovered ow model parameter function is reduced and thus provides betteraccuracy when forecasting ow.Optimally designing the geophysical survey for the coupled problem, while in-cluding historic survey data, produces a survey that adapts with the dynamics of thereservoir while reducing the number of measurements required. This research buildson A-optimal sparsity constrained design algorithms presented in (Haber et al., 2011)for large-scale ill-posed inverse problems, by introducing a monitor function to thedesign optimization problem to include historic data.The main body of the thesis consists of ve chapters. Since the same geophysicalsurvey and ow model are used to demonstrate the ideas in the proposed research,the thesis rst outlines these models, and then presents the new ideas.To this end, Chapter 2 introduces the seismic tomography geophysical imagingsurvey and the tracer advection ow model chosen to represent a simplied reservoirmonitoring scenario, presents a mathematical argument for not requiring a physicalproperty relationship between models, and details the discretization of the two modelPDEs.Chapter 3 outlines the decoupled linear inversion process and discusses the issueswith coupling the two diering physical models. In particular, consideration mustbe taken when characterizing the error in each of the models, as this leads to very131.4. Thesis objective and outlinedierent formulations of the coupled inverse problems.Chapter 4 presents the seismic tomography tracer advection coupled inverse prob-lem formulated to solve for the uid velocity eld and initial tracer concentrationwithin the reservoir. Because of the ill-posed nature of the problem, special consid-eration of the regularization terms required to solve the problem are discussed. Inparticular this chapter highlights the construction of a unique regularization functionfor the uid velocity eld. Due to the non-linearity of the regularization functions,solutions to the inverse problem are not so straightforward. This chapter thereforedetails the discretization of the regularization and optimization of the inverse prob-lems, and nally presents results for the tomography tracer advection model. Theobjective of the example is to show that estimates of the initial tracer distributionand uid velocity eld are improved as new geophysical data are obtained, and thatforecasting improves as a consequence.Chapter 5 details the formulation of the adaptive A-optimal experimental designmethod for two cases of the error characterization in the dynamic uid ow model.The new formulation proposes the addition of a monitor function in the denitionof the posterior covariance matrix of the estimated model parameters. The non-linearity of the resulting optimization problem requires iterative solution techniques,which are outlined in Section 5.3. The adaptive design method is then demonstratedfor the coupled seismic tomography tracer advection reservoir monitoring example forboth cases of the error characterization. The results show that experimental designscomputed using this method adapt to the ow within the reservoir while still providingsucient data for optimal recovery of tracer distributions.Chapter 6 summarizes the key ndings of the research presented in this thesis,discusses the practical contributions, and comments on areas of future work.14Chapter 2A model reservoir monitoringexperimentTo demonstrate the research ideas presented in this thesis a simple model of a reservoirmonitoring scenario is used throughout. The model consists of a seismic tomographysurvey to image the reservoir, and a tracer advection ow model to characterize thesubsurface uid ow dynamics.Although there are several choices when considering geophysical surveys for moni-toring subsurface ow, such as gravity surveys or electrical methods, borehole seismictomography surveys provide a good starting point for monitoring a reservoir. Inparticular, the tomography survey utilizes already existing boreholes for source andreceiver locations, while receivers can be placed on the surface to expand the mea-surement area. It has also been established that seismic wave velocity is sensitive tochanges in rock uid properties (Alumbaugh and Morrison, 1995; Nolet, 1987; Oliverand Chen, 2010), and thus the data are sensitive to ow. An additional benet isthe linearity of the governing equations. The measured tomography data are a linearmapping of the subsurface rocks seismic properties, which can be inverted to recovera spatial image of the properties. Thus imaging the seismic wave velocity involvessolving a fairly straightforward linear inverse problem.There are several well established and researched models describing ow within areservoir, such as two-phase ow models or a black oil model, among many others.15Chapter 2. A model reservoir monitoring experimentSee (Chen et al., 2006) for an excellent outline of ow models. However, many of thesemodels are complicated and can be discretized with a limited number of techniqueswith time step and spatial discretization length restrictions. Since the goal is tocouple the ow model and the geophysical model to form an already computationallyexpensive ill-posed inverse problem, it is not benecial at this stage of proof of conceptto incorporate more physically accurate, and computationally dicult models.Although the tracer advection model might seem like an unlikely choice for sim-ulating ow within a reservoir, there is some justication; the tracer itself can bethought of as an analogy to the second phase (a diering uid than the alreadypresent uid) moving in a reservoir and mathematically resembles two-phase owmodels (Chen et al., 2006). In reality, the tracer advection model is a single phaseow model. Also, at very long time scales, where an injected uid front might takemonths to move even a few meters, it is not necessary for the purposes of this researchto incorporate a model that characterizes minute changes. Additionally, in hydrolo-gy, tracer experiments are common practice for estimating aquifer parameters such ashydraulic conductivity, and have been used in conjunction with seismic tomographysurveys in the past (Hyndman et al., 1994).Practically, from a computational perspective, the tracer advection models lin-earity lends itself to a variety of discretization techniques, including unconditionallystable Lagrangian particle in cell methods. A second benet becomes clear in Sec-tion 2.3 when discussing physical parameter relationships and coupling the ow modelwith the geophysical imaging.This chapter introduces the governing partial dierential equations (PDEs) for thetomography survey experiment and tracer advection model, discusses the relationshipbetween the dierent physical properties of the two PDEs, and nally presents adiscretization of the resulting governing system of equations.162.1. Geophysics: seismic tomography2.1 Geophysics: seismic tomographyClassical travel time tomography uses ray paths to model travel times and so assumesthat the data are of high frequency, meaning the wavelength is at least several timessmaller than the spatial variations of the velocity model (Yilmaz, 2013).A typical 2D seismic tomography experiment places nssources on one side of theregion to be imaged and nrreceivers on the other. Travel times of an acoustic wavetraveling from source to receiver constitute the measured data, where the ultimate goalis to recover an image of changes to the distribution of physical properties generatingchanges in these data. A schematic of the seismic experiment is illustrated in Fig. 2.1.Figure 2.1: Tomography setting: sources (si) and receivers (rj) are placed in twoboreholes and the travel times of an acoustic wave traveling along adjoining ray paths(i;j) are measured.Mathematically, the travel time data are calculated by integrating the inverse ofthe acoustic wave velocity, also know as the slowness, m(t;~x), over a ray path i;j.172.2. Dynamics: tracer advectionThe travel time along the ray path from the ith source to the jth receiver is given bydi;j(t) =Zi;jm(t;~x) d`: (2.1)where m 2 M is the slowness model belonging to the subspace, M. Assuming alayered earth with homogeneous isotropic porous rock layers, then the ray paths thewaves follow will be linear and independent of the slowness (this is true for smallperturbations in the slowness eld). Thus, the problem can be cast as a linear inverseproblem for m(t;~x) (Jones, 2010).Given noisy data vectors d(t) measured at times (t0; : : : ; tn) Equation (2.1) canbe written as,Fm(tk;~x) + k= d(tk); (2.2)where d is a ns nrvector of real numbers, F :M! Rnsnris the forward mappingoperator, and measurement error k. The receivers, or modern seismometers, measuretiny displacements in the earth. These instruments use electronic sensors, ampliers,and recording devices, which can cover a wide range of frequencies. In general, therst arrival time and initial displacement of earth is measured successively manytimes over a very short interval and stacked such that the measurement is given as acontinuous random variable, with error i;j ffi;j.2.2 Dynamics: tracer advectionThe tracer advection equations as presented in (Chen et al., 2006), describe the trans-port of a solute in a fully saturated uid phase with the assumption that uxes dueto dispersion and diusion are small relative to the advective transport of the so-lute. The motion of a tracer with conserved concentration c (volumetric fraction in182.2. Dynamics: tracer advectionthe uid phase), in a uid of density (t;~x), with spatial and temporal variables~x; trespectively, is described by the following relation,@(c)@t+r c~u = 0; (2.3a)s.t. c(0;~x) = c0(~x);where the ow eld~u satisesr ~u = q; (2.3b)~u = K(~x)rp; (2.3c)p(0;~x) = p0(~x): (2.3d)accompanied with either Dirichlet or Neumann boundary conditions for the pressure,p. The uid velocity eld~u is determined by Darcy ow in this approximation, whereonly the hydraulic conductivityK = k(~x)=, a function of the permeability k and uidviscosity , and pressure gradient describe the ow eld. Recalling Equations (1.1b),the gravitational term has been excluded on the assumption that it is very small incomparison to the pressure gradient. This implies that pumping rates are high enoughto generate signicant pressure within the reservoir. Note that Equation (2.3b) impliesthat what ows in ows out within the entire domain, such that the source pumpingrate in is equal to the pumping rate out, qin= qout). Thus the rock is assumed to beincompressible, such that no uid is stored. The tracer advection model is illustratedin Figure 2.2.In a typical reservoir characterization setting one would collect pressure and con-centration measurements and invert these data to recover K(~x). However, assumingthe presence of the tracer generates a change in the seismic wave velocity, then the ideais to rst estimate u, and later estimate K. To do this, the petrophysical parameter192.3. Physical parameter relationsFigure 2.2: Tracer advection: Tracer of concentration (c) is advected in a uid ofdensity () with velocity eld (~u). The hydraulic conductivity (K) determines thevelocity eld given a pressure gradient within the reservoir.relationship between c(~x) and the seismic slowness m(t;~x) must be considered.2.3 Physical parameter relationsIn order to use geophysical techniques to estimate uid ow parameters through in-version of geophysical data, there must be a relationship between the ow parametersand the geophysical rock properties, also known as a petrophysical relationship. Onefamous petrophysical function, known as Archie's law, relates the electrical conduc-tivity of a sedimentary rock to its porosity and brine saturation (Archie, 1942). Inthis way, changes in saturation will result in changes to the electrical conductivity, aproperty that governs the ow of electric current in the presence of an electric eld.However, Archie's law is an empirical relation that is site dependent and requiressignicant laboratory work which cannot be generated for the general case (Archie,202.3. Physical parameter relations1942; Gosselin et al., 2003). Many empirical relations between seismic velocities androck porosity have also been developed (Abul, 2010; Mavko et al., 2009). However,these are also rock and site specic, and are often highly nonlinear.For the seismic tomography experiment it may be that in the presence of thetracer, the seismic slowness changes relative to areas where there is tracer present.Thus, it is assumed that the tracer c^ = c and the seismic slowness m are related suchthat petrophysical relationship m(c^) exists, however, given the choice of ow modelit is not necessary to know m(c^).To see this consider the case where the ow is divergence free such that r~u = 0.The uid velocity model, Equation (2.3b), typically has a very sparse right hand sidewith delta functions for injection and extraction points. Therefore, the ow eld isdivergence free almost everywhere within the reservoir.Recalling Equation (2.3a), the equation describing the advection of the tracer is@(c^)@t+r c^~u = 0:Now, assuming that the only source of change to the subsurface rock comes from thetransport of the tracer, then the \transport" of the seismic slowness is governed bythe same conservation law as follows,@m(c^)@t+r (m(c^)~u);applying the chain rule to the rst term and expanding the divergence term yields,@m(c^)@c^@c@t+~u>rm(c^) + c^r ~u:212.4. DiscretizationExpanding the gradient term and recalling the assumption r ~u = 0 gives,@m(c^)@c^@c@t+~u>@m(c^)@c^rc^;and nally factoring out@m(c^)@c^while recalling that Equation (2.3a) is equal to zero,@m(c^)@c^@c^@t+~u>rc^=@m(c^)@c^(0) = 0Thus, it can be assumed from this result that the transport equation describingthe ow of the tracer can also be used for the \transport" of the slowness. Thisfeature of the tracer advection model is particularly advantageous as it eliminates thenecessity of knowing the petrophysical relationship between the tracer and seismicvelocity.Given the two models, the discretization is now outlined in the following sections.2.4 DiscretizationThe numerical discretization of the geophysical imaging and the ow dynamics ispresented in three parts: the discretization of the tomography experiment (Equa-tion (2.2)), the tracer advection (Equation (2.3a)), and the uid conservation law(Equation (2.3b)).For simplicity, the discretization is restricted to a two-dimensional setting wherethe computational domain = [0; L1] [0; L2] is rectangular, and divided into l1byl2rectangular cells of edge length h1= L1=l1and h2= L2=l2.2.4.1 Discretization of the tomography equationTo discretize the tomography experiment with discrete data dkmeasured at time k =1; : : : ; n the line integrals in (2.1) are approximated by summing over line segments222.4. Discretizationdlk. That is, the approximation of the integral in Equation 2.1 for the travel time di;jof a wave traveling along the ray path i;jbetween the ithsource and the jthreceiverisdi;j=Zi;jm(t;~x)dl pXk=1mkdlk;where p is the number of cells that i;jintersects in the mesh, mkis the value of theslowness at the cell center of the kthcell, and dlkare the line segments of i;jpassingthrough each cell respectively.Thus, the discrete tomography experiment can be written as a linear systemFmk+ k= dk; k = 1; : : : ; n; (2.4)where the entry of Fijis the length dlkof the intersection of i;jthrough the kth cell.The cell centered values of the grid are stored in the vector m.The operator F : RM! RN, maps from the model space to the data spaceN = ns nr (number of sources number of receivers), and thus dkis an N 1vector of discrete travel time measurements, and m is an M 1 vector of seismicslowness values with units s=m.2.4.2 Discretization of the mass balanceSince vector operators such as the curl and divergence will be used later, a Marker andCell (MAC) (Fletcher, 2012) grid is chosen for the discretization of the uid velocityvector eld~u and the initial seismic slowness m0. The MAC method is discretizedon an Eulerian staggered grid system, and is a nite dierence solution technique forinvestigating the dynamics of an incompressible viscous uid. One key feature is theuse of Lagrangian virtual particles, whose coordinates are stored, and which move232.4. Discretizationmju+1u1u+2u2Figure 2.3: A grid cell with length h1and height h2, slowness value (mj) located atthe cell center, and components of the uid velocity owing in (~u = [u+1; u+2]) and out(~u = [u1; u2]) on the cell edges.from a cell to the next according to the latest computed velocity eld (McKee et al.,2008).Figure 2.3 pictures a cell in the grid where~u = [u1; u2] is discretized on the cellfaces, and m at the cell centers using a cell-centered grid function m0. This yields anapproximation of the divergence in the cell-centers using short dierences.2.4.3 Discretization of the ux-balance conservationA nite volume discretization of the ux-balance equation r ~u = q on a staggeredgrid is obtained by using the divergence theorem;ZVcellr ~u dV =ZScell~u ~n dS;where, Vcellis the volume of the cell, and Scellis the cell surface area. The vector~nis the unit normal vector to the surface. See (Haber and Ascher, 2001) for furtherdetails.Referring to Figure 2.3, the divergence over a cell with dimensions of h1 h2in242.4. Discretizationthe mesh can be expressed as(r ~u)cellh2(u1 u+1) + h1(u2 u+2)h1h2;where the u+iis the magnitude the ux into the cell and uiis the magnitude ofthe ux out of the cell in ith coordinate direction. This discretization leads to thefollowing standard discretization of the divergenceDIV u = q; where DIV =Il2dl1l1+1dl2l2+1Il1; (2.5)whereis the matrix Kronecker product, Il2 Rlldenotes the identity matrix, l isthe number of cells in a respective direction, and where dll+12 Rll+1is a short nitedierence matrix; see (Modersitzki, 2009) for more details about the implementation.The discrete divergence operator maps from faces of the mesh to cell-centers.2.4.4 A Particle-In-Cell method for ow simulationThe tracer advection equation (Equation (2.3a)) can be discretized in a number ofways. One option is to use explicit Eulerian techniques such as up-winding, Lax-Wendro or Lax-Friedrichs (Ascher et al., 2006). However, explicit techniques suerfrom one main shortcoming; since the velocity of the equation is unknown, one has tomonitor the time step making sure stability is maintained. Monitoring and adjustingthe time step can add enormous complexity to the optimization technique used if at-tempting to estimate the velocity. Methods such as up-winding and other ux limitershave another shortcoming; they are not dierentiable with respect to the velocity andtherefore simple optimization techniques can run into diculties. Alternately, onehas the option of using implicit methods which require the solution of linear systemsat every time step, or using semi-Lagrangian techniques.252.4. DiscretizationSemi-Lagrangian techniques are particularly attractive for this application as theyare explicit, unconditionally stable, and can be easily made dierentiable. Semi-Lagrangian methods can suer from low accuracy, however for this application wherethe velocity eld is approximated and thus inaccurate, working with highly accuratediscretizations is usually computationally wasteful. Therefore a Particle-In-Cell (PIC)method, that can also be interpreted as a semi-Lagrangian technique, is applied forthe discretization.The idea of using PIC methods for the solution of conservation laws is not new andwas proposed in Evans and Harlow (1957). Recently, PIC methods have emerged incomputer graphics (Edwards and Bridson, 2012) and in ow in porous media (Roubi-net et al., 2013).To use PIC methods, a particle xjis associated with the midpoint of the jth gridcell. The particle is assigned a value mj= m0(~xj), and the advection of the particlexjis described by the following ordinary dierential equation (ODE)@xj(t)@t=~u(x): (2.6)Using a midpoint quadrature rule, the particles position after the rst time step isapproximated byxj(t1) = xj(0) + t(Acfu)j+O(t2); (2.7)where Acfdenotes an averaging operator from cell faces (where the velocity resides)to cell centers. Thus, the position of the particle at a later time is easily computedgiven a time step and velocity eld. Since the new position will in general not be themid-point of a cell, the value mjis distributed to the closest surrounding points asillustrated in Figure 2.4. To compute the distribution weights, bilinear basis functions(bilinear interpolation) were used, however, higher order interpolation techniques canand have been applied (Edwards and Bridson, 2012). Thus, the slowness eld mkat262.4. Discretization~xA~xB~utFigure 2.4: Particle-In-Cell method for the advection and mass-conservation equation.A particle located in the cell-centered point~xAon the regular mesh is pushed forwardwith velocity~u along the path~ut, to the non-grid point~xB. The value (mass)associated with the particle is then transferred to the cell centered grid points adjacentto~xBusing interpolation weights.272.5. Summarytime tk= tk1+t is given bymk= T(tu)mk1; (2.8)where T(tu) denotes the push-forward matrix containing the values of the bilinearhat functions associated with the particles at the grid points. Since the value ofeach particle is spread between neighboring cell-centers, this construction guaranteesexact mass preservation as long as no ux exits at the boundaries of the computationaldomain. Using the push forward matrix, the time-stepping process is thenmk= T(tu)mk1= Tkm0: (2.9)The slowness at the kthtime step is therefore a pushed-forward version of m0, wherethe path the particle takes is piece-wise linear.2.5 SummaryThis chapter presented the physical models selected to represent the geophysical imag-ing survey and tracer advection uid ow model for the reservoir monitoring experi-ment, discussed the relationship between the model physical properties, and detailedthe discretization.One novel and important result of the tracer advection model was discussed inSection 2.3, is that under certain conditions the \motion" of the seismic slownessis governed by the advection model. This result eliminates the need to know thepetrophysical relationship between the tracer and the slowness, and will be benecialwhen formulating the coupled inverse problems later in Chapters 3 and 4.In addition, the discretization of the physical models was detailed, and in particu-lar, an explicit unconditionally stable Lagrangian technique was applied to discretize282.5. Summarythe tracer advection equation. Because the method is explicit and unconditionallystable, it will not add additional complexity to the inverse problems to be proposed.29Chapter 3Error characterization andmodel estimationEstimating either the seismic slowness, or the uid ow parameter functions and lateroptimally designing the geophysical survey requires the solution of an inverse problem.Before formulating the coupled inverse problem to estimate ow parameters, anddeveloping an optimal experimental design method, it is important to understand thebasis of a linear inverse problem from both the Bayesian and frequentist perspectives,and additionally, how the characterization of the error in the dynamic ow modeleects the formulations. The main goal of this thesis is to give the reader a quickbackground in linear inverse theory in preparation for the inverse problem developedin Chapter 4, and the adaptive optimal experimental design method developed inChapter 5.This chapter rst presents a brief review of Bayesian linear inversion and itsrelationship to frequentist linear inversion, followed by two formulations of a coupledinverse problem that assumes knowledge of the velocity eld~u, for the case wherethere is no noise in the ow model, and the case where there is random noise.3.1 Linear inversionThe process of inferring model parameters from a set of noisy observations is com-monly known as inversion, or parameter estimation. In this context it is assumed303.1. Linear inversionthat given a physical model and model parameters, observation data can be simulat-ed by solving the forward problem. Thus, if one has some real observation data, thenan estimate of the model parameters can be obtained by minimizing the dierencebetween simulated data (predicted data) and the observations.Inverse problems can be thought of as statistical estimation problems that canbe studied from both Bayesian and frequentist perspectives (Biegler et al., 2011; S-tark and Tenorio, 2010). In either case, a stochastic model for the data is requiredand constraints on the unknown model parameters can be incorporated. The majordierence in the two perspectives is that Bayesian methods require that constraintsbe formulated as a prior probability distribution, whereas the frequentist method-ology does not. An important well known observation worth noting, is that for alinear Tikhonov regularized inverse problem, with some simple assumptions, the twomethodologies reduce to an equivalent problem. This result becomes important whenlater in Chapter 5 when the adaptive optimal experimental design of the ill-posedcoupled problem is examined.To begin a basic review of inversion, consider the dynamic process and measure-ment technique described by the following set of discrete linear equations describingthe tracer advection (Equation (3.1a)), and the seismic tomography (Equation (3.1b)),presented in Chapter 2,mk= Tmk1+ k; (3.1a)dk= Fmk+ k; (3.1b)where the noise vectors kand kare assumed to be Gaussian normal, k N(0;Qk)and k N(0;W1k). If the noise is uncorrelated, then the matrixW is a diagonalmatrix with entries wi= 1=ff2i, where ffiis the standard deviation of iaway from313.1. Linear inversionzero.Assuming for the time being that the goal is to recover only the geophysical modelproperties given geophysical data without consideration of the dynamics, then thegeophysical properties, or slowness m for the tomography experiment, given traveltime data set d, can be estimated by maximizing the Bayesian posterior likelihood ofm given d, which is proportional to the product of the probability of the data giventhe model, and the model likelihood:P (mjd) P (djm)P (m);where m is a random variable with Gaussian prior distribution m N(;). Thusthe probability of the model given the data is,P (mjd) / exp (12(Fm d)>W(Fm d)) exp (12(m )>1(m )):Taking the negative log yields,ln(P (mjd)) /12(Fm d)>W(Fm d) +12(m )>1(m );and minimizing with respect to m givesF>W(Fm+ d) +1(m ) = 0:Finally, solving for m, yields the maximum a-posterior estimate (MAP estimate)bm = (F>WF+1)1(F>Wd+1):Alternately, a common linear inverse problem for estimating the geophysical pa-323.1. Linear inversionrameters m from the frequentist perspective is given by the following regularizedTikhonov inverse problem (Hansen, 1998; Tikhonov and Arsenin, 1977),minm12kFm dk2+12 kL(mmp)k2: (3.2)The rst term in Equation (3.2) penalizes the dierence between the simulated da-ta and measured data, and the second term (the regularization term) is chosen topromote smoothness between the estimated model and a reference model m0. Theregularization term improves the convexity of the inverse problem as it is ill-posed inthe classical sense; without it solutions will t the measurement noise and are thusmeaningless.The choice of regularization is an interesting problem in itself, and is usuallyproblem specic. There is much research in the area of optimally selecting and con-structing regularization functionals, particularly in medical imaging (Hansen, 1998;Huang et al., 2012). Later, when developing the coupled inverse problem, speciccare will be taken in constructing an appropriate regularization functional for reser-voir parameter estimation.Minimizing Equation (3.2) yieldsbm = (F>F+ L>L)1(F>d+ L>Lmp):Note that if L>L = 1,W = 1=I, and mpis considered the prior mean, then thetwo problems are equivalent.This well known fact will allow for the application of Bayesian optimal experimen-tal design techniques to coupled inverse problems with linear regularization terms laterdiscussed in Chapter 5.However, in order to couple the Equations (3.1) to form a single inverse problem,333.2. Coupled inversion and error characterizationthe error kin the dynamic model must be considered.3.2 Coupled inversion and error characterizationBefore combining Equations (3.1a) and (3.1b) there are two instances of the noise inthe dynamics that must be addressed which result in dierent coupled inverse problemformulations.Assume rst that the dynamics are exact, that is k= 0; 8 k. Such a scenario isoften assumed in history matching (Oliver and Chen, 2010). In this case the dynamicsare determined by the initial model m0, and estimated using the measurements forall times. The second instance refers to the case that k6= 0, and therefore modelsm1; : : : ;mkfor all times are estimated from all available data.3.2.1 Exact dynamicsSetting k= 0 and assuming the velocity eld u is known, the model at time tkcanbe written as a linear mapping of the initial model m0, such thatmk= T(tu)mk1; or mk= T(tu)km0; (3.3)which can be written as the following linear system,0BBBBBBBBBB@IT(tu) IT(tu) I......T(tu) I1CCCCCCCCCCA| {z }=:A(u)0BBBBBBBBBB@m1m2m3...mn1CCCCCCCCCCA| {z }=:m+0BBBBBBBBBB@T(tu)00...01CCCCCCCCCCA| {z }=:B(u)m0= 0:(3.4)343.2. Coupled inversion and error characterizationIn compact form the system is given byA(u)m+B(u)m0= 0:Thus it is possible to write the time evolution of the initial model m0in terms of thefuture models and the dynamics, such thatm = A(u)1B(u)m0: (3.5)The data measurements at each time are then given bydk= Fkm0+ k; (3.6)where Fk= FTk, and for the large system for all times,d = (IF)A(u)1B(u)m0+; (3.7)whereis the matrix Kronecker product, d = [d1; ::dk]>and = [1; ::k]>.The coupled linear Tikhonov inverse problem to estimate m0, assuming that thedynamics are known, is thenminm012(IF)A(u)1B(u)m0 dobs2+12 kLm0k2: (3.8)The case where the parameters governing the dynamics are unknown, that is, thereservoir parameters and velocity eld are to be estimated, again results in a dierentformulation of the inverse problem. A new formulation for estimating these modelsis presented in Chapter 4.353.2. Coupled inversion and error characterization3.2.2 Inexact dynamicsTo address the case where the dynamics are inexact the geophysical imaging and thetracer advection of m0are written for all timesFm0 d0= 0;Fm1 d1= 1; m1Tm0= 1;......Fmk dk= k; mkTmk1= k:(3.9)Collecting all equations in matrix form results in the following system,0BBBB@F...F1CCCCA0BBBB@m0...mk1CCCCA0BBBB@d0...dk1CCCCA=0BBBB@0...k1CCCCA; (3.10)0BBBB@T I......T I1CCCCA0BBBB@m0...mk1CCCCA=0BBBB@0...k1CCCCA:Dening the transport matrix T for all times asT =0BBBB@T I......T I1CCCCA;mk= (m>0; : : : ;m>k)>, d = (d>0; : : : ;d>k)>, = (>0; : : : ; >k)>, and = (>0; : : : ;>k)>the compact system is363.3. Summaryd = (IF)m+ ; (3.11)Tm = : (3.12)Assuming again that the velocity eld u is known, the most recent model mkisestimated by solving the following regularized minimization problem which includesall historic data and model estimates m0; :::mk1,bmk= argmin12(IF)mk dk2Wk+12T mk2Q1k(3.13)+2k(IL)mkk22:A second option is to re-estimate all models given all current data. This is akin tominimizing over mk. Writing the problem in this way, where the inverse problemminimizing both error vectors, is also known as Kalman smoothing (Aravkin, 2010;Evensen, 1994).3.3 SummaryThis chapter provides a brief review of linear inverse theory from both the Bayesianand frequentist perspectives, highlights the relationship between the two methodolo-gies for linear problems, and discusses the error characterization in the dynamic owmodel.Given the exact and inexact noise considerations presented, the main questions ofthe thesis can now be addressed. Chapter 4 considers the exact case and presents theformulation of the inverse problem for estimating both the initial state of the dynamicmodel, m0and the velocity eld~u. Chapter 5 presents a new optimal design method373.3. Summarydeveloped from the Bayesian perspective for the frequentist coupled inverse problem.The design method is presented for both exact and inexact noise realizations in thedynamic ow model.38Chapter 4Coupled inversion for velocityand initial slownessIn many cases in reservoir characterization, it is the hydraulic conductivity K that issought in order forecast ow. For the tracer advection model, one might argue thatthis is not necessary, and that it may be a simpler problem to estimate the velocityeld u. In particular, once the velocity eld is known, and given an estimate of theinitial slowness, the ow can be predicted by marching the forward in time. This isthe motivation for the approach detailed here.Whereas the previous chapter outlined a method for estimating the initial slownessm0given seismic travel time data while the velocity eld is assumed to be known,this chapter highlights a novel approach for estimating both m0and the velocityeld u from seismic tomography data, through a coupled inverse problem. Thisformulation assumes that there is no error in the dynamic ow model. That is,mk= T(u)mk1+ k, where k= 0.The mathematical formulation developed in this chapter is similar to that of Sec-tion 3.2.1, with the exception that the velocity eld is unknown. The resulting param-eter estimation problem is posed as a ow constrained inverse problem with speci-cally tailored regularization functionals for both the initial modelm0and the velocityeld u. Following the discretization of the regularization functionals, the initializationand numerical optimization of the resulting non-linear inverse problem is discussed.394.1. Flow constrained geophysical imagingThe chapter concludes with a numerical example, where the initial slowness modeland velocity eld are recovered and used to forecast the ow of the tracer.4.1 Flow constrained geophysical imagingThe formulation presented in Section 3.2.1 can be thought of as a ow constrainedinverse problem for the case where the velocity eld u is known. In the case whereone does not know the ow eld, the following constrained problem can be solved toobtain estimates of both m0and u,minu;m012kXj=0kFmj djk2(4.1a)s.t. mk= T(tu)mk1; m(0;x) =m0; DIV u = q: (4.1b)Recall Equation (3.4) such that mk= T(tu)mk1can be written in matrix formfor all times in terms of the initial condition m0as m = A(u)1B(u)m0. Withthe assumption that the dynamics are exact, ie. the noise vector k= 0;8k, and thedata for all times are given by Equation (3.7), then Equation (4.1) reduces tominu;m012k(IF)A(u)1B(u)m0 dk2(4.2a)s.t. DIV u = q: (4.2b)Minimizing Equation (4.2) under the ow constraints yields a velocity eld and initialslowness distribution that ts the tomography data. However, since both the tomog-raphy and ow estimation problems are ill-posed, regularization is needed. HereRow(~u) is dened as a regularization functional on the velocity eld, and Rm(m0)as a regularization functional for the initial slowness.To regularize the ow eld estimation the following continuous regularization func-404.1. Flow constrained geophysical imagingtional for~u is dened:Row(~u) =Z1ffi(r~u) +22w(~x) j~u~urefj2d~x; (4.3)where the rst term penalizes the curl of the ow eld and the second term seeks tominimize the dierence between the recovered velocity,~u, and a reference velocity,~uref. The parameters 1; 2> 0 balance the contribution of both terms.The choice to control the curl of the ow in the regularization was made sincethe divergence and the curl complement each other. Given that the divergence ofthe ow eld is set by the constraint (4.1b) it makes sense to regularize only overits orthogonal complement. The penalty function ffi was chosen by noting that sharpcontrasts in hydraulic properties of dierent rock units are common; for example seethe visualization of the r~u in the left column in Figure 4.2. Taking the curl of thevelocity eld~u = Krp with r ~u = q givesr~u = r (Krp) = rprK;where kr ~uk = krpkkrKk sin():Thus, the curl of the ow eld has tangential discontinuities where K has jumpsbetween rock types. For the most part the pressure gradient is parallel to rK exceptalong the boundaries of rock types. For example, an aquifer might be surroundedby non porous rock with zero hydraulic conductivity, and thus the cross productalong that boundary will not necessarily vanish. Thus, assuming that K is piecewiseconstant, the curl of the velocity is expected to be sparse.The following convex approximation to the `1-norm is applied to promote the414.2. Discretization of regularization functionalssparsity of the curl of u,ffi(c) =pc2+ : (4.4)The reference model,~uref, is included in order to incorporate prior informationabout the subsurface. The weighting function w(~x) quanties the condence in~uref;see (Oldenburg and Pratt, 2002). One option to compute~urefis by solving (2.3b)and (2.3c) using a reference conductivity model K0constructed from a priori bore-hole data. In this case, w(~x) is large close to the boreholes and grows smaller as thedistance from the known drill site grows.It was mentioned in Section 2.2 that the diusion term is not included in theadvection-diusion model. Because of this assumption, it is unlikely that the initialtracer model will have diuse edges. Thus, smoothed Total Variation (TV) is usedas the regularization to promote sharp edges in the estimated models (Ascher et al.,2006),Rm(m0) = Zffi(jrm0j) d~x: (4.5)Finally, the regularization parameters 1, 2and should be selected such thatthe data mist is approximately equal to the norm of the noise; a quantity that isin general unknown, (Parker, 1994). Therefore, a cooling strategy starting with largeregularization parameters that are then decreased incrementally until a reasonablysmall data mist is achieved is applied, (see Section 4.6).4.2 Discretization of regularization functionalsThe discretization of the regularization functionals applies standard nite dierenceapproximations of the partial dierential operators on orthogonal staggered grids,following the discretization of the tomography experiment and the ow equations424.2. Discretization of regularization functionalsoutlined in Chapter 2.To discretize the curl of the ow eld~u, Stoke's theorem is used,ZSr~u n^dS =I~u d; (4.6)where, n^ is the unit normal vector to the surface dS, is the path around the surface,and d is the innitesimal path length.The discretization of the curl operator follows (Haber and Ascher, 2001; Moder-sitzki, 2009), such that the discrete curl operator is dened byCURL =dl2l2+1Il1Il2dl1l1+1: (4.7)Here dll+12 Rll+1is again a short nite dierence matrix. Note that the CURLoperates from the cell faces to the nodes. For a complete derivation of the matricessee (Haber and Ascher, 2001; Modersitzki, 2009).The integral in Equation (4.3) is discretized using a midpoint rule, resulting inthe discrete regularization,Row(u) = 1h2e>Acnffi(CURL u) +2h22(u uref)>W (u uref); (4.8)where the matrix Acnaverages from nodes to cell-centers, h is the cell length, ande 2 RMis a vector of ones.Next, the regularization functional for the initial slowness is discretized. Notingthat m0is discretized on a cell-centered grid, a standard discrete approximation forthe smoothed total variation regularization is applied (Ascher et al., 2006)Rm(m0) = h2e>qAcf((GRAD m0) (GRAD m0)) + ; (4.9)434.3. Numerical optimizationwhere is the Hadamard product, and GRAD is a standard 2-point discretizationof the gradient of a cell-centered variable which maps from cell-centers to faces, asdescribed in (Ascher et al., 2006; Haber and Ascher, 2001). Acfis an averaging matrixfrom cell-faces to cell-centers. Note that the notation is somewhat abused, and thathere the square root of a vector is the point-wise square root.Now that the regularization functionals have been established and discretized, theoptimization methods are discussed in the following section.4.3 Numerical optimizationIn this section the approach to solving the ow constrained discrete optimizationproblem is outlined. A variable projection method is chosen to solve for the initialslowness and velocity eld in turn. Because the TV regularization is used for m0theobjective function is non-linear with respect to m0, and thus a primal-dual Newtonmethod (Chan et al., 1999) is applied. To estimate u an approximate SequentialQuadratic Programming (SQP) method is utilized, (see Nocedal and Wright (2000)for further details).The discrete form of the coupled variational optimization problem (4.2) for u andm0with the discrete regularization functionals is given by,minu;m012k(IF)A(u)1B(u)m0 dk2+Row(u) +Rm(m0) (4.10a)s.t. DIV u = q: (4.10b)There are a number of options for the solution of such problems. One option is tosolve the problem directly with respect to u and m0. This approach has a num-ber of disadvantages. First, it requires solving a large coupled problem where the444.3. Numerical optimizationparameters may have dierent scales. Second, it does not take advantage of theexisting ability to solve the decoupled problems eciently. Finally, solving the cou-pled problem requires the simultaneous evaluation of two regularization parameters.An alternate, more attractive option is to use a variant of the variable projectionmethod (Golub and Pereyra, 2003). This method was applied successfully in (Chunget al., 2006) for solving the related super-resolution problem. Furthermore, it hasbeen shown in (Chung, 2009) that it is possible to choose regularization parametersfor the dierent variables in the algorithm separately, thus decoupling the problemof selecting regularization parameters. The variable projection method as applied tothe tomography-ow optimization problem is as follows.First, assuming u to be xed, the conditions for a minimum with respect to m0areG(u)>(G(u)m0 d) +rm0Rm(m0) = 0; (4.11)where G(u) = (IF)A(u)1B(u). In the standard variable projection method,where quadratic regularization is used, Equatin (4.11) is linear with respect to m0and therefore can be solved directly. In this application the equation is nonlineardue to the TV regularization applied to m0. Nonetheless, it is possible to solve thisproblem rather quickly using, for example, a primal-dual Newton method (Chan et al.,1999).Assuming thatm0solves Equation (4.11), thenm0=m0(u) and the optimizationproblem can be rewritten as a problem for u aloneminu12kG(u)m0(u) dk2+Row(u) +Rm(m(u)) (4.12a)s:t: DIVu = q: (4.12b)454.3. Numerical optimizationIntroducing the Lagrange multiplier , the Lagrangian L(u;) of the problem isL(u;) =12kG(u)m0(u) dk2+Row(u) +Rm(m0(u)) + >(DIVu q): (4.13)An important observation made in (Golub and Greif, 2003) was that if m0solvesthe system Equation (4.11) thenruL(u;;m0(u)) =@L@u+@L@m0@m0@u=@L@u;because Equation (4.11) implies that@L@m0= 0 and thus m0can be treated as aconstant in the objective function.This observation leads to the following conditions for a minimumJ(u)>(G(u)m0(u) d) +ruRow(u) + DIV> = 0 and DIVu = q; (4.14)where J(u) is the Jacobian (sensitivity) of the data mist term with respect to u. Tocompute the Jacobian the mist is rewritten as,D(u) =12r(u)>r(u);where r(u) := (IF)A(u)1B(u)m0 d and J(u) =@r@u:Using (3.4) to simplify notation, and recalling the denitionm = A(u)1B(u)m0,the Jacobian is then@r(u)@u= (IF)@m(u)@u: (4.15)In order to dierentiatemwith respect to u, Equation (3.4) is implicitly dierentiated,464.3. Numerical optimizationyielding@(A(u)m)@u+A(u)@m(u)@u=@(B(u)m0)@u: (4.16)For the computation of the derivatives of A(u)m and B(u)m0only the push forwardoperator times the ithmodel, T(u)mi, needs to be dierentiated for indexes i =1; : : : ; n. These derivatives depend on the employed interpolation basis functions.Dierentiating the product T(u)miwith respect to u has been done in (Chung et al.,2006). Furthermore, since the interpolation matrix is sparse, the derivative matrices@(A(u)m)@uand@(B(u)m0)@uare sparse.To summarize, Equation (4.16) was solved for the derivative of m with respectto u and substituted into Equation (4.15) to obtain the derivative of the residualfunction r(u) asJ(u) =@r(u)@u= (IF)A(u)1@(B(u)m0)@u@(A(u)m)@u: (4.17)An important observation needs to be made here. While the matrix J(u) is densein general, its product with a vector can be computed eciently using sparse matrixtechniques. To calculate J(u) times an arbitrary vector z one rst computes thematrix vector producty =@(B(u)m0)@u@(A(u)m)@uz:Since this matrix is sparse the computation can be done eciently. Next, the solutionx = A(u)1y is obtained by solving the linear systemA(u)x = y:474.3. Numerical optimizationSolving this linear system is equivalent to solving a single ow forward problem, thatis, advection in time. Finally, the matrix F is also sparse and thus computing thematrix vector product, (IF)x can also be done quickly and eciently.To complete the computation of the derivatives the dierentiation of the regular-ization term is now discussed. It is straight forward to verify that@Row@u= 1CURL>diag(1=ffi(CURLu))CURLu; (4.18)where the notation 1=ffi is a vector that is divided point wise.Given all the components described above the goal is to now solve the discreteoptimization problem. There are a number of options for the solution of the problem.Here an approximation to the Sequential Quadratic Programming approach (SQP) isapplied. By linearizing the system in (4.12) and using a Gauss-Newton approximationof the Hessian with respect to u the following linear system is obtained,0B@J>J+ 1CURL>diag(1=ffi(CURLu))CURL DIV>DIV 01CA0B@u1CA= 0B@LuL1CA: (4.19)This system is solved for u and and a soft (backtracking) line search is used forthe update of u and (Nocedal and Wright, 2000).To solve the linear system (4.19) note that the system has many similarities tothe Stokes problem with a positive-semi denite (1,1) block, as outlined in Goluband Greif (2003). Such systems can be solved by a combination of the augmentedLagrangian method and an approximation to the Schur complement (Benzi et al.,2005).484.4. Initialization4.4 InitializationThe optimization problem given by Equation (4.10) is a nonlinear non convex prob-lem that without an appropriate initialization will result in solutions of low quality.Therefore an initialization methodology is used that cheaply yields a \reasonable"starting point.To this end, consider the decoupling of the imaging of the slowness for all times.That is, consider the individual problemsFmk+ = dk; k = 1; : : : ; n:Inverting each data set to obtain n initial estimates of mkis accomplished by solvingn decoupled optimization problems;bmk= argminmk12kFmk dkk2+Rm(mk); k = 1; : : : ; n: (4.20)Given the estimates ofbmkan estimate of the initialization velocity is computed bysolving the optimization problembu = argminu12kA(u)mB(u)m0k2+Row(u): (4.21)This estimate is equivalent to obtaining a ow estimate without improving m, whichin this experiment leads to a good initial guess for both the slowness and velocity.In the following sections the approach to jointly estimate the initial slowness andthe velocity eld are presented for a small numerical example. The results show thatthe recovered ow eld and reconstructed initial slowness can be used to predict theow of an injected tracer within the subsurface.494.5. Experimental setup4.5 Experimental setupThe computational domain, = [0; 100][0; 200] meters, is divided into l = [200; 100]cells of width h = [1; 1]. The ground truth conductivity model consists ve layers,with conductivities ranging between 10 and 1000 m3day=kg. The ground truth modelis picture in the top panel of Figure 4.1. Two boreholes are used for the injectionand extraction of uids. A reference conductivity model was constructed by linearinterpolation of the borehole data. The reference conductivity model is a simplelayer of high conductivity (k = 1000m3 day=kg) surrounded by a background layer(k = 10m3day=kg). Both the true and reference conductivity models are pictured inFigure 4.1. The ground-truth initial tracer mgtis a piece-wise constant model withtwo regions of concentration 0.5 and 1; see top left plot in Figure 4.2. Fluid is injectedfrom the left well at 50 m depth and extracted on the right well at 60 m depth witha static pumping rate of 100 m3= day.Given the hydraulic conductivity model and the source term, the ground-truthow eld ugtand the reference ow eld urefare obtained by solving Equation (2.3);see quiver plot in Figure 4.1. Tomography travel-time data are simulated by solvingthe forward problem given by Equation (3.7) for mgtand urefand adding Gaussianwhite noise with a standard deviation ff = 0:5. The tomography experiment consistsof ns= nr= 35 transmitters and receivers equally spaced from 20m to 100m depthalong the left and right boundary of the domain, respectively. Data were simulatedfor 15 days with a time step t = 1 day. Snapshots of the advection of the tracer canbe seen in the rst column of Figure 4.3.4.6 Sequential reconstruction as initializationTo obtain a starting guess of the ow eld and the initial slowness m0, the proce-dure outlined in Section 4.4 is applied. First, the slowness evolution is individually504.6. Sequential reconstruction as initializationground truth hydraulic conductivityreference hydraulic conductivity 1002003004005006007008009001000Hydraulic conductivity (m3 day/kg)Figure 4.1: Visualizations of the ground true and reference hydraulic conductivitymodels and associated ow elds. The true and reference models are discretized ona 100 200 cell grid on a domain of 100m depth by 200m across. A static source,at 55m depth on the left border, and a sink, at 65m depth on the right border ofthe domain, are indicated by red. The reference model is a simplied model whichinterpolates borehole data between the two wells.514.6. Sequential reconstruction as initializationm0true model initialization coupled estimate uuyuxCURL uFigure 4.2: Visualization of ground truth data for the domain with dimensions 100m(y direction) depth and 200m (x direction) across (rst column), individual recon-struction (second column) and results of the joint reconstruction approach. The rstrow visualizes the initial slowness (m0). The second shows a quiver plot of the recov-ered velocity eld (~u). The components of the velocity eld [ux; uy] are pictured inthe third and fourth row and the last row shows the curl of the velocity eld (r~u).Note that there is almost no ow in the y direction, and that at the boundaries ofthe reservoir there are non-zero values of the curl of~u524.6. Sequential reconstruction as initializationday 1true model initialization coupled estimateday 10day 15day 25day 35Figure 4.3: The evolution of the tracer is pictured for the ground-truth (left column),individual reconstruction used for initialization (middle column), and the proposedjoint reconstruction. The tracer evolution is simulated for 40 days with the respectiveinitial slowness and ow eld. The front of the ground truth tracer is visualized bydashed lines for the rst 3 rows. Note that the joint method results match the groundtruth tracer front to a greater degree than the predictions from the initializationestimates.534.7. Joint reconstructionreconstructed for all 15 days solving Equation (4.20) using a standard Gauss-Newtonmethod with a regularization parameter = 0:4. Subsequently, a ow estimate bu wasobtained by solving Equation (4.21) using the described SQP method with regular-ization parameters 1= 103and 2= 4 103. The regularization parameters werefound using a cooling strategy which starts with large parameters that are decreasediteratively until the measurement noise became visually too prominent in both thereconstructed images and the estimated ow elds.A larger bias towards the reference ow eld, uref, is given close to the boreholesby using the weighting function w(x1; x2) = (x2 100)2=1002.Estimates of the initial slowness and ow eld are visualized in the second columnof Figure 4.2. The rst row shows the reconstruction of the initial slowness, estimatedby inverting the tomography data of the rst day. The ow eld estimate, visualizedin the second, third, forth, and fth rows; in a quiver plot, component wise and by itscurl, captures the main characteristics of the ground truth. However, the magnitudeof the velocity eld is small compared to the ground truth. This can be seen bylooking at the plots of uyand in the second column of Figure 4.3, where the frontof the high concentration plume is behind the ground truth. This demonstrates theimportance of using the ow to better guide the imaging which in turn yields moreaccurate ow predictions.4.7 Joint reconstructionTo improve the estimates ofm0and the velocity eld, the coupled tracer tomographyinverse problem (Equation (4.1a)) is solved. The regularization parameters wereadjusted to 1= 5 103; 2= 7 104and = 200. The result of the individualestimate of bu, was used as a starting guess for the optimization and ureffor thereference model.544.8. ForecastingSince errors in the ow eld accumulate in time, an outer loop is used to introducethe tomography data sequentially as the reconstruction of the velocity eld improves.To be precise, only the tomography data from the rst day is used in the estimationproblem for m0, and then ow eld is estimated by solving Equation (4.12) using theSQP method outlined earlier. Afterwards the estimate m0is updated by introducingtomography data from one additional day and updating the estimate of u. In total,ve outer iterations are performed using tomography data from the rst six days.Note that all tomography data is used in the ow estimation at all iterations. Itis important to note that computational time increases with the coupled recoverymethod since it requires solving the System (4.19) in addition to the imaging step.The initial slowness and velocity eld reconstruction results are visualized in thethird column of Figure 4.2. Note the signicant improvement in the reconstructionof m0, particularly the sharpening of the edges.4.8 ForecastingThe numerically estimated ow elds and initial slowness are used to forecast theevolution of the tracer beyond the rst 15 days by solving the forward ow model inEquation (2.9) for 35 days. The numerical solutions of the initialization and the jointhistory matching problems are compared to the ground-truth in Figure 4.3.In the considered application, the most interesting quantity is the prediction ofthe arrival time. The front of the tracer is thus visualized by a dashed line in thesubplots. It can be seen that the coupled estimate of u and m improves the predictionof the arrival time when compared with the estimates obtained in the initialization.554.9. Summary4.9 SummaryThis chapter presented a new method for the estimation and prediction of tracer owin porous media given seismic tomography geophysical data. The new coupled methodwas formulated as a ow constrained inverse problem to estimate both the initialseismic slowness and the velocity eld in which tracer is injected. As a consequenceof the coupled ill-posed inverse problem, a novel regularization for the velocity eldwas additionally constructed to promote discontinuities at rock property interfaces.The coupled method was demonstrated for the model borehole seismic tomographysurvey of tracer advection. The numerical results were compared with numericalresults obtained through a decoupled process and demonstrate that this new coupledapproach yields not only accurate reconstructions of the initial slowness and ow eld,but, especially, accurate predictions of the uid ow velocity and tracer evolution.In addition to coupling the ow and geophysics, the next approach to improv-ing reservoir monitoring is to consider what data are collected as the tracer movesthrough the subsurface. This optimal experiment design question is investigated inthe following chapter.56Chapter 5Dynamic optimal experimentaldesignThe previous chapter presented a new joint method for estimating the initial slownessm0, the velocity eld u, given all previous seismic tomography data sets for the reser-voir monitoring experiment outlined in Chapter 2. Although this method alleviatessome of the uncertainty and expense involved in reservoir monitoring programs bymaking use of geophysical data instead of sparse uid measurements, it is possiblethat the tomography data set could be reduced while maintaining an optimal estima-tion of the initial slowness. In this chapter, the idea is therefore to reduce the numberof necessary tomography data measurements in an optimal way. This is commonlyknown as optimal experimental design.The optimal design method will generate an experiment with a reduced numberof measurements while maintaining the best estimate of the initial slowness, incorpo-rate the historic tomography data sets, and couple the ow dynamics such that theexperimental design adapts with the moving tracer. For the formulation of the designmethod, the assumption is made that dynamics are well known, that is, it is assumedthat~u is known, or estimated to reasonable accuracy for forecasting.The new approach, dened as adaptive optimal experimental design, applies amethod similar to classic A-optimal design criteria, by minimizing the adapted meansquare error (amse) of a regularized model estimate, instead of the mean square575.1. Adaptive experimental designerror (Atkinson and Donev, 1992; Fedorov, 1972). The adapted mean square erroris dened by introducing a monitor function which scales the mean squared erroraccording to historic model estimates.This chapter introduces the adaptive design method from rst principals, thenumerical optimization of the resulting design objective function, and discusses theresults for the tomography tracer advection reservoir model example.5.1 Adaptive experimental designTo introduce the adaptive experimental design method it is necessary to understandA-optimal design criteria from rst principals. To begin, consider only two measure-ment vectors of the formF1m+ 1= d1and F2m+ 2= d2; (5.1)where F1is the forward tomography operator at time t1and F2is the forward operatorat time t2. As before, the measurement error vectors kare assumed to be normallydistributed uncorrelated noise, k N(0;W1k), with zero mean and a diagonalcovariance matrixW1k, for k = 1; 2.At this point consider only the static case, where the goal is to design the exper-iments at times t1and t2such that the \best" recovery of m is obtained by somecriteria, such as minimizing the Tikhonov optimization problem. This scenario tsthe problem where the dynamical system is \exact" (containing no noise as in Sec-tion 3.2.1) and m represents the dynamical system parameters, such as the initialcondition of the tracer in the reservoir modeling experiment.Two dierent designs can be considered. First, it is possible to perform a-prioridesign, that is, to design the experiments F1and F2prior to the data collection.585.1. Adaptive experimental designThis is the case when the time between t1and t2is much shorter than the timefor the processing of the data. A second approach is to use post-priori design. Theidea here is to use the results obtained from the rst experiment in order to designa \better" second experiment. The goal is to rst obtain an a-priori design for therst experiment and then, after the data is collected and processed, use the estimatorobtained at t1in order to design the data collection at time t2.Before discussing the design of the data acquisition for time t2, it is necessaryto review the a-priori design of the experiment at time t1. Consider the regularizedestimation of m given d1, accomplished by solving the following penalized-weightedleast squares (Tikhonov regularized) optimization problembm1= argmin12kF1m d1k2W1+2kL(mbm0)k22; (5.2)whereW1= diag(w1) is a matrix of inverse variances, that is, the inverse covarianceof the noise, L is a smoothing penalty matrix,bm0is the current estimate of the modelprior to having data, and is a regularization parameter. Recall from Chapter 3 thatfrom the Bayesian perspective,bm0can be considered the mean of a prior probabilitydistribution with covariance matrix (L>L)1. Minimizing Equation (5.2) yields theestimatorbm1= (F>1W1F1+ L>L)1(F>1W1d1+ L>Lbm0): (5.3)Dening the matrix C = (F>1W1F1+ L>L), and recalling that F1m+ 1= d1theerror in the recovery can be written asbm1m = C1F>1W1F1m+C1F>1W11+ C1L>Lbm0(5.4)+ (C1L>Lm C1L>Lm)m:595.1. Adaptive experimental designCollecting terms and using the denition of C gives,bm1m = C1F>1W11+ C1L>L(bm0m): (5.5)Squaring and taking the expectation over 1and recalling that 1 N(0;W11), themean square error is then given bymse(w1;m) = Ekbm1mk22= trace[F1C2F>1W1] + 2EkC1L>L(bm0m)k22:(5.6)Taking the Bayesian point of view assumes thatmbm0is Gaussian with a zero meanand a covariance matrix (L>L)1such that the mean squared error is equivalent tothe Bayesian risk.mse(w1) = ffi1(w1) = trace[C2F>1W1F1] + trace[C2L>L]: (5.7)Finally, using the linearity of the trace and the denition of C results inffi1(w1) = traceh(F>1W1F1+ L>L)1i: (5.8)That is, the Bayesian risk is the trace of the inverse precision matrix, or trace of theposterior covariance matrix; a well known result.The goal of the design problem seeks to obtain a better estimate form by assumingthat the (inverse) variances, w1, of the collected data can be controlled in some way.By controlled, it is assumed that a data measurement is accompanied by a s-tandard deviation determined by the instrumentation. If this is not the case, severalmeasurements at a particular location can be conducted to estimate the error. Similartreatment is given in (Alexanderian et al., 2014; Haber et al., 2008, 2010, 2011) where605.1. Adaptive experimental designthis point is further discussed. In this way, through optimization, the weights wiareestimated prior to the actual data collection, and a measurement dithat is assignedwith innite standard deviation (wi= 0), or zero variance, by the optimization isthus not collected.To estimate the variances, the Bayesian risk ffi1(w1) is minimized with respectto w. Additionally, since the goal is to obtain a sparse design (that is, collect onlya few measurements) an additional cost on w1is added. This cost is equivalent tothe 1-norm of w and promotes its sparsity (see (Alexanderian et al., 2014; Haberet al., 2008) for further discussion), resulting in the following optimization problemwhich balances the minimization of the mean squared error (mse) and the cost of theexperiment (Haber et al., 2008),ffi1(w1) = traceh(F>1W1F1+ L>L)1i+ e>w1; 0 w1;i; (5.9)where e>is a vector of ones.Consider now using the design problem of estimatingm given the estimated modelbm1. One could rewrite the recovery problem in a similar way, that isbm2= argmin12kF1m d1k2W1+12kF2m d2k2W2+2kL(mbm0)k22; (5.10)which leads to the estimatebm2(w2) = (F>1W1F1+ F>2W2F2+ L>L)1(F>1W1d1+ F>2W2d2+ L>Lbm0):(5.11)Note that w1is assumed to be xed and therefore the new estimate is a function ofw2alone. If the steps above are repeated then the A-optimal design for time t2yields615.1. Adaptive experimental designthe minimization of the functionffi2(w2) = traceh(F>1W1F1+ F>2W2F2+ L>L)1i+ e>w2; 0 w2;i: (5.12)At this point it is worth pausing for a moment and noting an important feature ofthe experimental design criteria. The design criteria do not depend on the data. Thisobservation is true for any of the design criteria (that is, C,D and E designs)(Fedorov,1972). Furthermore, it is easy to verify that any linear estimator of the data yields acovariance matrix that is independent of d1. This implies that using current designcriteria does not make use of the estimated model obtained at time t1to obtain abetter estimate at time t2.In order to use the information obtained at time t1for the design of the experimentat time t2the concept of adaptive design is presented.Assume that the model,bm1has some \interesting" features and some \boring"features. To be more specic, assume that the dierence1= jbm1bm0j (5.13)is small in some norm over a region and large in others. The goal then is to betterestimate the new features that appear in the model.To this end, the monitor function, fi () is introduced to measure the change inthe model. For example, to start, consider the functionfi = ; (5.14)that simply measures the change in the estimator compared to the previously knownestimator. A dierent function that was found to be useful in numerical experiments625.1. Adaptive experimental designisfi = (); (5.15)where is a smoothed characteristic function. In this case the optimization focusesonly on areas where is large and does not sample areas where is small.Given the monitor function, rather than minimizing the mean square error, theidea is to minimize the adaptive mean square error (amse), dened byamse(w2;m) = Ekbm2mk2fi1= E[(bm2m)>diag(fi1)(bm2m)]: (5.16)The idea behind the amse is to obtain a tighter bound on the model where theestimator exhibits large changes compared to the known a-priori estimator.Starting from Equation (5.4), modifying it to deal with the amse, and repeatingthe calculation above, the expectation over the amse isffi2(w2) = tracehdiag(fi1)(F>1W1F1+ F>2W2F2+ L>L)1i: (5.17)Similar to the non-adaptive case, the adaptive design is computed by minimizing the(penalized) function ffi2(w2).The above concept can easily be extended to any number of time steps. Considerk 1 experiments that are conducted at times t1; ::tk1using the forward operatorsF1; : : : ;Fk1and assume that the goal is to design an experiment for the problem attime tk. The design criteria in this case readsffik(wk) = trace24diag(fik1)kXj=1F>jWjFj+ L>L135+ e>wk; 0 wk;i;(5.18)635.2. Design for dynamical systemswhere fik1is the monitor function that measures the dierence between the modelestimated at time tk1and the model estimated at time tk2.5.2 Design for dynamical systemsThe previous section presented a formulation for adaptive A-optimal experimentaldesign that assumes the model, m is static. This section therefore addresses the casewhere the model is governed by a (potentially noisy) dynamical system.Recalling Chapter 3 and the system of Equations (3.1), the adaptive optimaldesign method is applied for both cases, where the dynamic tracer advection modelis assumed to be exact or inexact. Note that in these cases the velocity eld u isassumed to be known.5.2.1 Exact dynamicsConsidering the data vector given by Equation (3.6) in Chapter 3, and noting thatthis case is exactly the case described in the previous section, the amse criteria canby directly applied.Given the assumption that there is no error, the estimation of the initial modelm0given k data sets isbm0;k= argmin12kXj=1kFjm0 djk2Wj+2kL(m0bm0)k22: (5.19)and therefore the design function is identical to Equation (5.18)ffik(wk) = trace24diag(fik1)kXj=1F>jWjFj+ L>L135+ e>wk; 0 wk;i:(5.20)645.2. Design for dynamical systems5.2.2 Inexact dynamicsTo address the case where the dynamics are inexact recall Equations (3.11) and (3.12),which correspond to the noise vectors for the tomography imaging experiment andthe tracer dynamics model.In order to apply the amse method the problem is formulated as in Kalman s-moothing (Aravkin et al., 2013; Kalman, 1960). A similar approach that maximizesthe information gain was presented in (Gharamti et al., 2015) for the ensemble Kalmanlter.Recalling that the noise vector k N(0;Qk), and assuming that all Qkareknown, that is, that the dynamics are not changing in time, all modelsbmk=(bm>0; : : : ;bm>k)>are estimated by minimizing Equation (3.13) presented in Section 3.2.2,bmk= argmin12(IF)mk dk2Wk+12T mk2Q1k+2k(IL)mkk22;such that,bmk=(IF>)Wk(IF) +T>Q1kT+ IL>L1(IF>)Wk(IF)dk;(5.21)whereWk= diag(w0; : : : ;wk), and Q1k= blkdiag(Q11; : : : ;Q1k).It is assumed in this case that the rst experiment d0, will either be conductedsuch that all data are collected, or will be conducted in a naive way. In either case itis assumed that the initial experimental design,W0= diag(w0), is known or can bedesigned for as was suggested in previous work (Haber et al., 2011).In a straightforward extension of the previous section, the amse design optimiza-655.3. Numerical optimizationtion function for a design at time tkis thenk(wk) = tracediag(fik)(IF)>Wk(IF) +T>Q1kT+ I(L>L)1+ e>w;0 wk;i: (5.22)wherefik= (fi0(bm0); : : : ; fik(Tbmk1)).It is assumed here that all previous experiments 0; : : : ; k 1 have already beenconducted. Therefore, solutions are only for the current time point tk. At this point,bmkis not known. To predictbmkat time tkthe previous estimator is propagatedahead in time,bmk Tbmk1.The two formulations of the adaptive design method for a dynamical system yieldtwo dierent numerical optimization problems. However, in the limit that all Qk!0, the inexact optimization problem reduces to the exact formulation, see (Nocedaland Wright, 2000) for a discussion of penalty methods. The numerical optimizationtechniques for the solution of these problems are discussed in the following section.5.3 Numerical optimizationThis section presents the calculation of the gradients required to solve the designoptimization problems for both the exact and inexact tracer dynamics.5.3.1 Exact dynamicsFinding an optimal design for time step k requires the minimization of the objectivefunctional given by Equation (5.20), which involves computing the trace of a largedense matrix. Typically, when solving geophysical problems, the size of matrices canbe on the order of millions. Thus, forming and storing these matrices is memory andcomputationally intense.665.3. Numerical optimizationIn the non-linear optimization algorithm, at each iteration the computation of theobjective function and its gradient is required. In order to avoid numerous expen-sive calculations while storing and manipulating large dense matrices, a stochasticHutchinson trace estimator is applied to approximate the trace (Haber et al., 2011;Hutchinson, 1990). By doing so, the objective function reduces to matrix vectorproducts which can be carried out much more eciently.If the vectors v1; : : : ;vmare independent and each with independent entries takingthe values of 1 and 1 with equal probability, then the stochastic trace estimator ofthe trace of a matrix H is given bytrace(H) 1mmXi=1v>iHvi:Applying the trace estimator for m = 1, Equation (5.20) becomesffik(wk) = v>diag(fi )F>kWkFk+G1v + e>wk; (5.23)where e is a vector of ones, G =Pk1j=1F>jWjFj+ L>L, and v is a random vectorwith equally distributed values of 1 and 1.Note that computing the the gradient of ffik(wk) involves computing the gradientof an inverse matrix. This computation is carried out by dening z such thatz = (F>kWkFk+G)1v , (F>kWkFk+G)z = v; (5.24)recalling thatWk= diag(wk), and dierentiating implicitly to obtain thatF>kdiag(Fkz) + (F>kdiag(wk)Fk+G)J = 0: (5.25)675.3. Numerical optimizationSolving for the Jacobian J =dzdwkresults in,J = (F>kdiag(wk)Fk+G)1F>kdiag(Fkz): (5.26)Substituting J back intorffik, transposing, dening the matrixC = F>kWkFk+G,and the vectors, y and z by,Cz = v and Cy = fi v;yields the gradient of the design function,rffik(wk) = (Fkz) (Fky) + e; (5.27)where is the Hadamard product.Unfortunately due to the inclusion of the monitor function fi , the computation ofthe gradient is slightly more complicated compared with the classic A-optimal designcase. Here the optimization requires solutions to two linear systems, Cz = v andCy = fi v compared with a single linear system for the classic A design. However,since only matrix vector products are required, the conjugate gradient method is usedto solve these systems avoiding ever forming C, since it requires only matrix vectorproducts.Due to the non-linearity of the design problem, iterative solution methods arerequired. Here, a gradient descent method is used for the solution of the problemwith a backtracking line search to solve for wk.685.3. Numerical optimization5.3.2 Inexact dynamicsTo minimize the noisy design function k(wk), the Hutchinson trace estimator isagain used to approximate the trace operation in Equation (5.22),k(wk) = v>diag(fik)(IF)>Wk(IF) +T>Q1kT+ I(L>L)1v;(5.28)wherefi = (fi0; fi1; : : : ; fik)>and v is a (k + 1)N 1 vector of evenly distributedvalues of [1; 1]. As in the noiseless case, again note that the computation of thegradient of kinvolves the computation of the gradient of a large dense inversematrix. Therefore, the vector z = (z0; : : : ; zk)>is dened such that(IF)>Wk(IF) +T>Q1kT+ I(L>L)z = v: (5.29)Dierentiating with respect to wk, and solving for J =dzdvresults inJ = (IF)>Wk(IF) +T>Q1kT+ I(L>L)1B: (5.30)whereB =0BBBBBBB@00...F>diag(Fzk)1CCCCCCCA: (5.31)695.4. Numerical examplesSubstituting J back into the gradient of the design objective function and transposinggives,rk(wk) = B>y+ e; (5.32)wherey =(IF)>Wk(IF) +T>Q1kT+ I(L>L)>(fi v): (5.33)To solve for z and y, the preconditioned conjugate gradient method was used to avoidforming the large matrices. The steepest descent method was used for the non-lineardesign optimization.5.4 Numerical examplesThe adaptive design method is demonstrated for the reservoir monitoring simulationusing the tomography experiment and tracer advection. The survey was designed tofully demonstrate the ability of the method to track motion and areas of interest inthe model. The example consists of a circular volume of tracer moving vertically. Thedynamics governing the motion of the tracer are described by the tracer advectionmodel, and imaged using the borehole seismic tomography survey.The partial dierential equations governing both the seismic tomography experi-ment (Equation (2.1)) and the tracer dynamics, (Equation (2.3)) are discretized on a2D computational domain = [0; 400] [0; 100]m, divided into [100; 50] cells of widthh = [4; 2]m. The uid velocity eld u was computed by solving Equations (2.3b) and(2.3c) with a constant hydraulic conductivity for 10m3 day=kg, and pumping rate of150m3=s.705.4. Numerical examplesFigure 5.1: Initial experiment setup. Tomography sources are pictured in green onthe left and receivers in blue on the right. The initial setup covers the entire owdomain with 20 sources and 30 receivers. There is one source of uid, pumped at aconstant rate from the top green point, and a sink at the bottom (red) of the domain.Flow is in the downward direction.715.4. Numerical examples5.4.1 RegularizationTwo dierent inversions were carried out to recover models mkby applying twodierent regularization functions. First a linear inversion was performed where theregularization operator L = I was assigned to promote smallness in the recoveredmodels. Alternately, L = GRAD is a common choice to promote smooth edges in themodels. However, since the original model has sharp edges, using a smooth recoverymay lead to more inaccurate results.Second, a non-linear inversion was performed with smoothed total variation (TV)as the regularization (Ascher et al., 2006). The total variation regularization pro-motes discontinuous boundaries, or sharp edges. This is particularly valid for thisexperiment since the tracer advection model does not include a diusion term. Thetracer model is concentrated in the circle and zero everywhere else. In this ideal casethere should be no diusion of the tracer, and thus the boundaries should remainsharp as the tracer moves along in the uid velocity eld.In continuous space, smoothed total variation is described by the following rela-tion,R(m) =Zffi(jrmj) d~x; (5.34)where the convex function ffi(c) =pc2+ , is an approximation to the `1norm.A standard discrete approximation for the smoothed total variation regularizationfound in (Ascher et al., 2006) is,R(m) = h2e>qAcf((GRAD m) (GRAD m)) + ; (5.35)where the gradient operator, GRAD, maps from cell centers to cell faces, and theaveraging matrix Acfmaps from faces to cell centers.725.4. Numerical examplesAlthough the optimal design method was formulated based on a linear estimationof mk, it is not unreasonable to estimate mkusing the optimal design for the linearmodel estimation problem with a non-linear regularization, provided that this willcontribute to recovering the best estimates of mk. As can be seen in Figure 5.2below, the overall relative model error for estimates recovered with the total variationregularization is signicantly lower than that of the models estimated using the linearregularization. Thus the estimates obtained using the TV regularization are moreaccurate with respect to the true model for this particular problem. Since the monitorfunction has a very large impact on the future optimal designs, it is desirable to getthe best possible estimate of mk.Figure 5.2: Plot of the relative error (kmtrue1bm1k=(kmtrue1k)) versus the regular-ization parameter for both the linear regularization and smoothed total variation.The TV regularization gives a better estimate of the model for the case of a sharptarget.735.5. Design results: exact dynamics5.4.2 Monitor functionTo design only for the tracer and not the entire reservoir, the monitor function deltais dened by = jmkmbj, where mbis the background model (that is, the modelwithout the dynamic target). Because neither the estimated model nor the back-ground model are perfect, a threshold value was chosen to remove unwanted orerroneous information from where the background model mb, is the model used tocompute the velocity eld in Equation (2.3c).5.5 Design results: exact dynamicsFor the exact case, data were simulated by marching the tracer along in time for timesteps of 25 days and measured by conducting a tomography survey at each time with4% Gaussian noise added to each data set. In total 9 experiments were conducted fortimes t0; t1; :::; t8.The initial design for time point t0did not include the dynamics. This amounts tosetting fi0= 1 everywhere. The design optimization problem in this case is identicalto the static case and produces an optimal design based only on the physics andinformation from the linear regularization operator L = I. This is apparent in boththe plot of the weights for experiment 1 at time t0in Figure 5.3 and the image of therays in Figure 5.5, where the design species fewer rays that cover the entire domain.For each experiment the amse was plotted versus the number of nonzero weightsfor a set of penalty parameters , see Figure 5.3 for examples for times t0; t2; t4. Fromthese curves, the best set of weights were chosen such that the weights were sparse,but also so the amse was kept reasonably small. After a set of weights was chosen,both m0and the current model mkwere reconstructed from the reduced set of datadk, in addition to reconstruction obtained from the original A-optimal design from t0,and using all data for both the linear regularization and TV. The values of d which745.5. Design results: exact dynamicsFigure 5.3: Exact dynamics: plots of the adapted mean squared error amse vs. thenumber of nonzero weights (left column), and the weights used to conduct experiments1, 3, and 5, at times t0; t2; t4.755.5. Design results: exact dynamicscontribute to the reduced data set correspond to the non-zero weights of wk. Thetotal number of data required to image the models is signicantly reduced. In mostcases the number of data are on average 60 of a possible 600.The model error was then calculated for each case with respect to the groundtruth, true model for each experiment, and is plotted in Figure 5.4. It is clear fromthe gure that models estimated using all data have the lowest mse compared with thetrue model. However, those constructed with the reduced set, particularly when usingthe TV regularization, are not far o. This implies that one can gain a reasonableestimate of the model from a signicantly reduced set of data, without giving up toomuch on accuracy.In Figure 5.5 the rays corresponding to dkused to estimate the current modelsare pictured in blue over an image of the model. Note that the rays tend to follow thetarget as it moves through the domain and pass through the tracer. There are somespurious rays that do not pass through the tracer which are included in the optimaldesign. However, since the design depends on the monitor function fi (bm0), and thereconstructions ofbm0are never exact, these spurious rays are expected.It is also apparent that the number of spurious rays increases further in time eventhough estimates ofm0improve as more data are collected. In particular, at time t6itappears that the design algorithm has a harder time generating a design that capturesthe target well. There are many more rays which do not pass through the target. Thisis partially due to the increasing number of multiplications by the transport matrixT as time progresses, compounding errors, but also the loss of mass as the tracerexits the domain out of the sink. This is also apparent in the error calculations inFigure 5.4, where even while estimating models using all data, the error increases overtime. However, the designs still provide enough information to recover models whichindicate the location of the tracer. The reconstructions of mkwith smoothed totalvariation as the regularization are much closer to the true models, even when the765.5. Design results: exact dynamics0 5 10experiment (k)00.10.20.30.40.50.60.70.80.91relative errorLinear Regularization||mktrue-mkAll data||/||mktrue||||mktrue-mkadaptive linear||/||mktrue||||mktrue-mkA-Opt linear||/||mktrue||0 5 10experiment (k)00.10.20.30.40.50.60.70.80.91relative errorTotal Variation Regularization||mktrue-mkAll data TV||/||mktrue||||mktrue-mkadative TV||/||mktrue||||mktrue-mkA-Opt TV||/||mktrue||Figure 5.4: Exact dynamics: relative error plots. Model error is compared for eachexperiment for both the linear and TV regularization775.6. Results: inexact dynamicsoptimal design was computed for the linear regularization model estimation problem.Again, this is to be expected since the true model has discontinuous edges, and sincethere is no diusion in the ow model.5.6 Results: inexact dynamicsFor the rst experiment w0, the initial experiment computed in the noiseless examplewas used in the subsequent design estimations. The covariance matrices Qkwereassigned a scalar value for the variance for each time step to represent independentlyand identically distributed (iid) noise in the dynamics, such that Qk= ffkI. A con-stant value was assigned for the standard deviation, ffk= 0:2, for all experiments.Models were reconstructed using the non-linear smoothed total variation regulariza-tion, and also with the linear regularization. The model estimation error is plottedin Figure 5.6.Models were estimated using all data, the initial A-optimal design for each exper-iment, and nally using the adaptive method reduced data set. It is again clear thatinversion carried out using all data generated models with the lowest mse, however,the benet of reducing the number of measurements to approximately 60 from 600might be much greater than the loss in accuracy in the model estimates, depending onthe cost of data acquisition. The designs, and models estimates are pictured in Figure5.7, below. Notice that in this case there are fewer spurious rays. It is apparent thatadapted A-optimal designs were again able to track the motion of the tracer with asignicantly reduced set of measurements.5.7 SummaryIn this chapter, a new method for the design of experiments for dynamical systemswas presented. The method generates survey designs which adapt to the motion of785.7. Summaryexp:1 exp:2 exp:3 exp:4 exp:5 exp:6 exp:7 exp:8 exp:9mtrue0 mtrue1 mtrue2 mtrue3 mtrue4 mtrue5 mtrue6 mtrue7 mtrue8mTV0 mTV1 mTV2 mTV3 mTV4 mTV5 mTV6 mTV7 mTV8#d=297mlinear0#d=69mlinear1#d=73mlinear2#d=61mlinear3#d=77mlinear4#d=55mlinear5#d=76mlinear6#d=61mlinear7#d=35mlinear8Figure 5.5: Exact dynamics: optimal designs and recovered models for 9 experiments.The top row shows the rays and the recovered model, with the number of rays (#d =297) given below, followed by the true models in the second row, the models recoveredusing total variation in the third row, and nally the models recovered using the lineargradient regularization.795.7. Summary0 5 10experiment (k)00.10.20.30.40.50.60.70.80.91relative errorLinear Regularization||mktrue-mkAll data||/||mktrue||||mktrue-mkadaptive linear||/||mktrue||||mktrue-mkA-Opt linear||/||mktrue||0 5 10experiment (k)00.10.20.30.40.50.60.70.80.91relative errorTotal Variation Regularization||mktrue-mkAll data TV||/||mktrue||||mktrue-mkadative TV||/||mktrue||||mktrue-mkA-Opt TV||/||mktrue||Figure 5.6: Inexact dynamics: relative error plots. Model error is compared for eachexperiment for both the linear and TV regularizationspecic areas of the model while incorporating historic data. The motivation for notsimply applying static optimal design methods to the dynamic problem is that themean squared error of the current estimated model never contains information fromprevious data. To this end, the design optimization problem \knows" nothing aboutthe changes in the model due to the dynamics and thus designs that are based onclassical formulations only reect information from the physics and the regularization(or prior) of the estimation problem.The adaptive optimal experimental design approach is based on the introductionof a monitor function to scale the mean squared error of an estimator whose motion isgoverned by a dynamical system, according to historic data. To determine a design fora future experiment, the adaptive mse (amse) is minimized with the added constraint805.7. Summarythat the design is sparse. Two model estimation problems were considered, whoseconstruction depends on the characterization of the error in the dynamical system.This leads to experimental designs that track the changes in the model with a reducedset of measurements. The methodology was tested using seismic tomography to imagethe advection of a tracer in a reservoir.815.7. Summaryexp:1 exp:2 exp:3 exp:4 exp:5 exp:6 exp:7 exp:8 exp:9mtrue0 mtrue1 mtrue2 mtrue3 mtrue4 mtrue5 mtrue6 mtrue7 mtrue8mTV0 mTV1 mTV2 mTV3 mTV4 mTV5 mTV6 mTV7 mTV8#d=297mlinear0#d=79mlinear1#d=59mlinear2#d=67mlinear3#d=81mlinear4#d=62mlinear5#d=53mlinear6#d=51mlinear7#d=61mlinear8Figure 5.7: Inexact dynamics: optimal designs and recovered models for 9 experi-ments. The top row shows the rays and the recovered model, with the number ofrays given below. The second row pictures true models, the third shows the modelsrecovered using total variation, and nally the bottom shows models recovered usingthe linear gradient regularization.82Chapter 6ConclusionsThe main objective of the research presented in this thesis was to improve reservoirmonitoring and forecasting by coupling geophysical survey techniques with a dynamicuid ow model and optimally designing the geophysical survey.Forecasting of ow within a reservoir requires knowledge of uid ow parameters,such as hydraulic conductivity. Estimates of such parameters can be obtained by di-rect measurements from subsurface rock samples, or though the inversion of ow data.However, these estimates are often highly inaccurate due to sparse measurements ofboth the parameters and ow data.Therefore, one approach investigated in this thesis is to include a higher samplingof measurements without adding a large additional cost to the monitoring program,was to couple a geophysical survey with the dynamic uid ow model in a singleinverse problem, such that the unknown ow parameters could be estimated fromgeophysical data.Geophysical surveys tend to cover much larger spatial areas at lower cost sincesensors can be deployed at the surface as well as below ground. Collecting uid owdata requires expensive boreholes, and is thus sparsely sampled. In the context of anill-posed under-determined inverse problem, it is very dicult to estimate parametermodels for large subsurface domains with little data. Thus, the greater number ofdata and spatial coverage improves the estimates obtained by solving the coupledill-posed inverse problem.83Chapter 6. ConclusionsThe second approach to improve reservoir monitoring proposed in this thesis wasto compute optimal geophysical survey designs for data measurements for the coupledsubsurface ow inverse problem. In this case, it is inecient to be collecting data overregions of the subsurface where there is no change to the geophysical properties dueto ow. Ideally, one would want to be able to adapt the survey with the changesin the subsurface ow. However, traditional optimal design methods do not utilizehistoric estimates of the subsurface models in a way that includes them in the designoptimization problem. To this end, the idea was to formulate an optimal designestimation method for the coupled geophysics and ow model inverse problem thatincludes historic model estimates.To demonstrate the ideas proposed in the thesis, the physical models for thegeophysical survey and the uid ow model were rst presented in Chapter 2. Thechoice of models was particularly important for the coupling problem. The seismictomography experiment was chosen primarily due to its linearity, but also due the factthat data can be collected from existing boreholes as well as from sensors placed on thesurface. Borehole seismic tomography is also commonly used in practice, with a widevariety of applications, making it an ideal survey for the model reservoir monitoringproblem.The tracer advection model was used because of its linearity and because it is verysimilar mathematically to more complex multiple phase ow models. In addition, aswas shown in Section 2.3, the advection of the tracer can also be thought of as theadvection of the seismic slowness. This result eliminates the need for a petrophysi-cal function relating the geophysical properties (seismic slowness) to ow parameters(tracer or hydraulic conductivity). Technically, the tracer advection model lends itselfto ecient numerical optimization schemes, in particular a stable and dierentiablediscretization of the physical models was applied. For the discretization of the owa Particle-In-Cell method was used, whose stability is independent of the time-step84Chapter 6. Conclusionssize and therefore can be safely applied, even if the uid velocity eld is unknown.The uid velocity eld was discretized on staggered grids, resulting in a stable dis-cretization of the regularization operators and ow constraints, see Section 2.4 fordetails.Following the introduction of the physical models, Chapter 3 provides a briefintroduction to linear inversion and highlights the dierent formulations which resultfrom dierent characterizations of the noise in the dynamical uid ow model: zeronoise and non-zero noise. The choice of noise characterization is important not onlybecause it results in dierent formulations of the inverse problem, but also because itis often dicult to estimate a priori.In Chapter 4, a new methodology for estimating the initial tracer concentrationand uid velocity eld from geophysical data was presented for the case where thenoise in the dynamical ow model was assumed to be zero for all times. The newestimation problem was formulated as a constrained inverse problem with a speci-cally tailored regularization for the velocity eld that promotes discontinuities in thetangential components of the estimated ow eld. Estimates of the velocity were thenused for the estimation of the hydraulic conductivity.Numerical results for a simple layered earth model with a homogeneous isotropicreservoir demonstrated that this new approach yields not only accurate estimates ofthe initial change to the slowness but, especially, accurate predictions of the uid owvelocity tracer evolution. The work was published in the SIAM Journal on ScienticComputing (Fohring et al., 2014). However, a number of topics relating to the coupledgeophysics ow problem have yet to be studied.One important aspect of the inverse problem is the evaluation of uncertainty inthe ow. Although the error in the dynamic ow model was assumed to be zero, itis likely that measurement and numerical errors were propagated with the dynamicsthroughout the inverse problem. This was particularly apparent when the number85Chapter 6. Conclusionsof historic experiments, and thus time steps, increased. With increased times step-s, came an increase in the diculty in recovering an accurate approximation to thetrue velocity eld. This is counter to what one would expect; that more informationshould provide better results, and is likely due to the noise propagation. In partic-ular, the particle in cell (PIC) discretization of the tracer advection equation relieson a bilinear interpolation. Bilinear interpolation error is quadratic, and thus witheach time step, and therefore each multiplication of the tracer by the interpolationmatrix, the error grows. Thus experimenting further with higher order interpolationsmay help to reduce error propagation. However, higher order interpolation functionsmay simultaneously add complexity to the inverse problem because of the dicultyassociated with computing gradients.An additional topic of research involves applying the coupled problem with morecomplex ow models. However, the use of such models might require the knowledgeof petrophysical parameter functions and the application of diering discretizationtechniques that are not necessarily stable for all time step lengths. This would leadto added complexity in the inverse problem, which would require innovative solutiontechniques.Finally, to truly demonstrate the ability of the method to recover uid ow pa-rameters, the conduction of a eld or laboratory experiment to obtain a real data setis required to fully evaluate the contribution of the work to reservoir monitoring andgeophysical history matching.To address the second goal of the thesis, Chapter 5 presents a new adaptive optimalexperimental design method for the linear coupled reservoir monitoring problem. Thenew adapted mean square error (amse) method builds on classical A-Optimal designcriteria by altering the posterior mean squared error (mse) of a model estimate toinclude historic data. In the classic case the mse is minimized to recover an optimal set86Chapter 6. Conclusionsof measurements. Since the mse for the linear estimator includes only the physics andthe prior covariance matrix, designs do not adapt with the motion of the tracer. Theamse was therefore dened to include historic data through an introduced monitorfunction. Results presented in Section 5.4 demonstrate the amse method for thecoupled seismic tomography tracer advection. Applying the method produced designsthat tracked the motion of the tracer, while requiring signicantly fewer measurementsto recover images of the tracer compared with using all data available. The amseoptimal design method research presented in the thesis is currently in the nal reviewprocess for publication in SIAM Journal on Uncertainty Quantication (Fohring andHaber, 2016).The idea to include historic information in the amse through the introduction ofa monitor function originates from adaptive mesh methods. In this thesis only twooptions for the monitor function were chosen and tested. Thus there is signicantroom for further investigation of how diering monitor functions can aect optimaldesigns.The adapted design criteria presented in the thesis, only included reducing thenumber of measurements from a maximum set. However, in many cases it might bemore ecient to eliminate entire borehole locations, instead of just individual sensorlocations, and additionally impose a monitory cost on more expensive measurementlocations to quantitatively estimate cost reduction.Another avenue of investigation includes designing for dynamic coupled non-linearinverse problems. Current non-linear design methods are limited by the fact that thereis no closed form calculation of the posterior mean squared error. This results inexpensive bi-level design optimization problems that require training sets of possiblemodels. However, for the coupled dynamic problem there maybe room to includehistoric estimates through propagation of earlier time models to generate trainingsets.87Chapter 6. ConclusionsRecently (Alexanderian et al., 2015) presented a method for A-optimal design fornon-linear inverse problems where the posterior covariance matrix is approximated bya stochastic average over a set of data by an approximate Hessian. The method wasdemonstrated for an advection diusion ow problem. Building on this method mightbe an interesting approach to introducing the adaptive design method for non-linearproblems.One might also ask how sensitive the adaptive optimal design is to variations inthe uid velocity eld. This could be investigated by estimating and updating theow dynamics in time, and would amount to developing an adaptive optimal designmethod for the constrained variable projection problem, that is, for the non-linearow coupled inverse problem presented in Chapter 4.88BibliographyAbacioglu, Y., Oliver, D., and Reynolds, A. (2001). Ecient reservoir history match-ing using subspace vectors. Computational Geosciences, 5(2):151{172.Abul, F. (2010). 4D Seismic History Matching Using the Ensemble Kalman Filter(EnKF): Possibilities and Challenges. PhD thesis, University of Bergen.Aggarwal, R., Demkowicz, M., and Marzouk, Y. M. (2015). Information-DrivenExperimental Design in Materials Science. In Information science for materialsdiscovery and design, volume 225, chapter 2, pages 13{44. Springer InternationalPublishing.Ajo-Franklin, J. B. (2009). Optimal experiment design for time-lapse traveltime to-mography. Geophysics, 74(4):Q27{Q40.Alexanderian, A., Petra, N., Stadler, G., and Ghattas, O. (2014). A-Optimal designof experiments for innite-dimensional Bayesian linear inverse problems with regu-larized l0 - sparsication. SIAM Journal on Scientic Computing, 36(5):2122{2148.Alexanderian, A., Petra, N., Stadler, G., and Ghattas, O. (2015). A fast and scal-able method for a-optimal design of experiments for innite-dimensional bayesiannonlinear inverse problems. SIAM Journal on Scientic Computing, pages 1{29.Alumbaugh, D. L. and Morrison, H. F. (1995). Monitoring subsurface changesover time with cross-well electromagnetic tomography. Geophysical Prospecting,43(7):873{902.89BibliographyAravkin, A. (2010). Robust Methods for Kalman Filtering / Smoothing and BundleAdjustment. PhD thesis, University of Washington.Aravkin, A. A. Y., Bell, B. B. M. B., Burke, J. J. V., and Pillonetto, G. (2013).Kalman smoothing and block tridiagonal systems: new connections and numericalstability results. arXiv preprint arXiv:1303.5237.Archie, G. (1942). The Electrical Resistivity Log as an Aid in Determining SomeReservoir Characteristics. Petroleum Technology, 146(October):54{62.Ascher, U. M., Haber, E., and Huang, H. (2006). On ective methods for im-plicit piecewise smooth surface recovery. SIAM Journal on Scientic Computing,28(1):339{358.Atkinson, A. C. and Donev, A. N. (1992). Optimum Experimental Designs. ClarendonPress.Bardow, A. (2008). Optimal experimental design of ill-posed problems: The METERapproach. Computers & Chemical Engineering, 32(1-2):115{124.Benzi, M., Golub, G., and Liesen, J. (2005). Numerical solution of saddle pointproblems. Acta numerica.Biegler, L., Biros, G., Ghattas, O., Heinkenschloss, M., Keyes, D., Mallick, B., Mar-zouk, Y., Tenorio, L., Waanders, B. v. B., and Willcox, K., editors (2011). Large-Scale Inverse Problems and Quantication of Uncertainty. John Wiley & Sons.Cassiani, G., Bohm, G., Vesnaver, A., and Nicolich, R. (1998). A geostatisticalframework for incorporating seismic tomography auxiliary data into hydraulic con-ductivity estimation. Journal of Hydrology, 206(1-2):58{74.Chaloner, K. and Verdinelli, I. (1995). Bayesian Experimental Design : A Review.Statistical Science, 10(3):273{304.90BibliographyChan, T. F., Golub, G. H., and Mulet, P. (1999). A Nonlinear Primal-Dual Method forTotal Variation-Based Image Restoration. SIAM Journal on Scientic Computing,20(6):1964{1977.Chen, Z., Huan, G., and Ma, Y. (2006). Computational Methods for Multiphase Flowsin Porous Media. SIAM.Cherno, H. (1972). Sequential analysis and optimal design, volume 8. SIAM.Chung, J., Haber, E., and Nagy, J. (2006). Numerical methods for coupled super-resolution. Inverse Problems, 22(4):1261{1272.Chung, J. M.-L. (2009). Numerical Approaches for Large Scale Illposed Inverse Prob-lems. PhD thesis, Emory University.Cockett, R. and Haber, E. (2016). A Numerical Method for Large Scale Estimation ofDistributed Hydraulic Conductivity from Richards Equation. Journal of Hydrology,pages 1{15.Curtis, A. (1999). Optimal experiment design: cross-borehole tomographic examples.Geophysical Journal International, 136(3):637{650.Doyen, P. M. (1988). Porosity from seismic data: A geostatistical approach. Geo-physics, 53(10):1263{1275.Edwards, E. and Bridson, R. (2012). A high-order accurate particle-in-cell method.International Journal for Numerical Methods in Engineering, 90(9):1073{1088.Emerick, A. a. and Reynolds, A. C. (2012). History matching time-lapse seismic datausing the ensemble Kalman lter with multiple data assimilations. ComputationalGeosciences, 16(3):639{659.91BibliographyEvans, M. W. and Harlow, F. H. (1957). The Particle-in-Cell Method for Hydrody-namic Calculations. Science (New York, N.Y.), 178(4066):76.Evensen, G. (1994). Sequential data assimilation with a nonlinear quasi-geostrophicmodel using Monte Carlo methods to forecast error statistics. Journal of Geophys-ical Research, 99(C5):10143.Fedorov, V. (1972). Theory Of Optimal Experiments. Elsevier.Fletcher, C. (2012). Computational Techniques for Fluid Dynamics 2: Specic Tech-niques for Dierent Flow Categories. Springer Science & Business Media.Fohring, J. and Haber, E. (2016). Adaptive A-optimal experimental design for lineardynamical systems. SIAM Journal on Uncertainty Quantication, xx:1{19.Fohring, J., Haber, E., and Ruthotto, L. (2014). Geophysical imaging of uid ow inporous media. SIAM Journal on Scientic Computing, 36(5):218{236.Fohring, J., Ruthotto, L., and Haber, E. (2013). Geophysical Imaging, ReservoirHistory Matching and Forecasting. In 2013 SEG Annual Meeting, Houston. Societyof Exploration Geophysicists.Gerritsen, M. G. and Durlofsky, L. J. (2005). Modeling Fluid Flow in Oil Reservoirs.Annual Review of Fluid Mechanics, 37(1):211{238.Gharamti, M. E., Mazouk, Y. M., Huan, X., and Hoteit, I. (2015). A greedy approachfor placement of subsurface aquifer wells in an enseble ltering framework. InDynamic Data-Driven Environmental Systems Science, volume 8964, pages 301{309. Springer International Publishing, Switzerland.Golub, G. and Greif, C. (2003). On solving block-structured indenite linear systems.SIAM Journal on Scientic Computing.92BibliographyGolub, G. and Pereyra, V. (2003). Separable nonlinear least squares: the variableprojection method and its applications. Inverse Problems, 19(2):R1{R26.Gosselin, O., Total, E., Uk, P., Aanonsen, S. I., Aavatsmark, I., Hydro, N., andCominelli, A. (2003). History matching Using Time-lapse Seismic ( HUTS ). InSPE Annual Technical Conference and Exhibition, pages 1{15.Haber, E. and Ascher, U. (2001). Fast nite volume simulation of 3D electromag-netic problems with highly discontinuous coecients. SIAM Journal on ScienticComputing.Haber, E., Horesh, L., and Tenorio, L. (2008). Numerical methods for experi-mental design of large-scale linear ill-posed inverse problems. Inverse Problems,24(5):055012.Haber, E., Horesh, L., and Tenorio, L. (2010). Numerical methods for the de-sign of large-scale nonlinear discrete ill-posed inverse problems. Inverse Problems,26(2):025002.Haber, E., Magnant, Z., and Lucero, C. (2011). Numerical methods for A -optimaldesigns with a sparsity constraint. Computational Optimization and Applications,52(1):293{314.Hansen, P. (1998). Rank-decient and discrete ill-posed problems: numerical aspectsof linear inversion. SIAM.Hoversten, G. M., Cassassuce, F., Gasperikova, E., Newman, G. a., Chen, J., Rubin,Y., Hou, Z., and Vasco, D. (2006). Direct reservoir parameter estimation usingjoint inversion of marine seismic AVA and CSEM data. Geophysics, 71(3):C1.Huan, X. and Marzouk, Y. M. (2013). Simulation-based optimal Bayesian experimen-tal design for nonlinear systems. Journal of Computational Physics, 232(1):288{317.93BibliographyHuang, H., Haber, E., and Horesh, L. (2012). Optimal estimation of L1-regularizationprior from a regularized empirical bayesian risk standpoint. Inverse Problems andImaging, 6(0):447{464.Hubbard, S. S. and Rubin, Y. (2000). Hydrogeological parameter estimation usinggeophysical data: a review of selected techniques.Hutchinson, M. (1990). A stochastic estimator of the trace of the inuence matrixfor laplacian smoothing splines. Communications in Statistics - Simulation andComputation, 19(2)(April 2015):433{450.Hyndman, D. W., Harris, J. M., and Gorelick, S. M. (1994). Coupled seismic and trac-er test inversion for aquifer property characterization. Water Resources Research,30(7):1965.Jones, I. F. (2010). Tutorial: Velocity estimation via ray-based tomography. FirstBreak, 28(2):45{52.Kalman, R. E. (1960). A New Approach to Linear Filtering and Prediction Problems1. Transaction of the ASME - Journal of Basic Engineering, 82(Series D):35{45.Khodja, M. R., Prange, M. D., and Djikpesse, H. A. (2010). Guided Bayesian optimalexperimental design. Inverse Problems, 26(5):055008.Lucero, C. (2013). Bayes risk A-optimal experimental design for ill-posed inverseproblems. PhD thesis, Colorado School of Mines.Lumley, D. E. (2001a). The next wave in reservoir monitoring: The instrumented oileld. The Leading Edge, 20(6):640{648.Lumley, D. E. (2001b). Timelapse seismic reservoir monitoring. Geophysics, 66(1):50{53.94BibliographyMavko, G., Mukerji, T., and Dvorkin, J. (2009). The Rock Physics Handbook. Cam-bridge University Press.McKee, S., Tom, M. F., Ferreira, V. G., Cuminato, J. A., Castelo, A., Sousa, F. S.,and Mangiavacchi, N. (2008). The MAC method. Computers and Fluids, 37(8):907{930.Mezghani, M., Fornel, A., Langlais, V., Lucet, N., and Francais, I. (2013). HistoryMatching and Quantitative Use of 4D Seismic Data for an Improved ReservoirCharacterization. In SPE Annual Technical Conference and Exhibition. Society ofPetroleum Engineers.Modersitzki, J. (2009). FAIR: Flexible Algorithms for Image Registration. SIAM.Nvdal, G., Mannseth, T., and Vefring, E. H. (2002). SPE 75235 Near-Well ReservoirMonitoring Through Ensemble Kalman Filter. SPE International, 75235(1):1{9.Nenna, V., Pidlisecky, A., and Knight, R. (2011). Application of an extended Kalmanlter approach to inversion of time-lapse electrical resistivity imaging data for mon-itoring recharge. Water Resources Research, 47(10):1{13.Nocedal, J. and Wright, S. (2000). Numerical Optimization. Springer Science &Business Media.Nolet, G. (1987). Seismic Tomography: With Applications in Global Seismology andExploration Geophysics. Springer Science & Business Media.Oldenburg, D. and Pratt, D. (2002). Geophysical inversion for mineral exploration.Geophysical Inversion Facility.Oliver, D. S. and Chen, Y. (2010). Recent progress on reservoir history matching: areview. Computational Geosciences, 15(1):185{221.95BibliographyOliver, D. S., Reynolds, A. C., and Liu, N. (2008). Inverse Theory for PetroleumReservoir Characterization and History Matching. Cambridge University Press.Papalambros, P. Y. and Wilde, D. J. (2000). Principles of optimal design: modelingand computation. Cambridge university press.Parker, R. L. (1994). Geophysical Inverse Theory. Princeton University Press.Pollock, D. and Cirpka, O. a. (2012). Fully coupled hydrogeophysical inversion of alaboratory salt tracer experiment monitored by electrical resistivity tomography.Water Resources Research, 48(January):1{13.Pukelsheim, F. and Studden, W. J. (1993). E-Optimal Designs for Polynomial Re-gression. The Annals of Statistics, 21(1):402{415.Roubinet, D., de Dreuzy, J.-R., and Tartakovsky, D. M. (2013). Particle-trackingsimulations of anomalous transport in hierarchically fractured rocks. Computers &Geosciences, 50:52{58.Sagnol, G. and Harman, R. (2014). Optimal Designs for Steady-state Kalman lters.ZIB Report, 39(October).Sarma, P., Durlofsky, L. J., Aziz, K., and Chen, W. H. (2013). A New Approachto Automatic History Matching Using Kernel PCA. In SPE Reservoir SimulationSymposium. Society of Petroleum Engineers.Slater, L., a.M Binley, Daily, W., and Johnson, R. (2000). Cross-hole electrical imag-ing of a controlled saline tracer injection. Journal of Applied Geophysics, 44(2-3):85{102.Stark, P. B. and Tenorio, L. (2010). A Primer of Frequentist and Bayesian Inference inInverse Problems. Large-Scale Inverse Problems and Quantication of Uncertainty,pages 9{32.96BibliographySteklova, K. and Haber, E. (2015). Joint Hydrogeophysical Inversion : State Estima-tion for Seawater Intrusion Models in 3D. Journal of Hydrology.Tenorio, L. (2001). Statistical Regularization of Inverse Problems. SIAM Review,43(2):347{366.Tikhonov, A. and Arsenin, V. (1977). Solutions of Ill-Posed Problems. Winston, NewYork.Trani, M. (2012). From time-lapse seismic inversion to history matching of waterooded oil reservoirs. PhD thesis, Delft University.Vasco, D. W., DattaGupta, A., Behrens, R., Condon, P., and Rickett, J. (2004). Seis-mic imaging of reservoir ow properties: Timelapse amplitude changes. Geophysics,69(6):1425{1442.Vauhkonen, M., Karjalainen, P. a., and Kaipio, J. P. (1998). A Kalman lter ap-proach to track fast impedance changes in electrical impedance tomography. IEEEtransactions on bio-medical engineering, 45(4):486{93.Vesnaver, A. L., Accaino, F., Bohm, G., Madrussani, G., Pajchel, J., Rossi, G., andMoro, G. D. (2003). Timelapse tomography. Geophysics, 68(3):815{823.Wilkinson, P. B., Uhlemann, S., Meldrum, P. I., Chambers, J. E., Carriere, S., Oxby,L. S., and Loke, M. H. (2015). Adaptive time-lapse optimized survey design forelectrical resistivity tomography monitoring. Geophysical Journal International,203(1):755{766.Wilt, B. M., Morrison, H. F., Becker, A., and Torres-verdin, C. (1995). Crossholeelectromagnetic tomography : A new technology for oil eld characterization. TheLeading Edge, 14(3):173{177.97BibliographyYilmaz, O. (2013). Seismic data analysis Processing, In version, and Interpretationof Seismic Data. Society of Exploration Geophysicists.98
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Adaptive optimal experimental design and inversion...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Adaptive optimal experimental design and inversion of a coupled fluid flow and geophysical imaging model… Fohring, Jennifer 2016
pdf
Page Metadata
Item Metadata
Title | Adaptive optimal experimental design and inversion of a coupled fluid flow and geophysical imaging model for reservoir monitoring |
Creator |
Fohring, Jennifer |
Publisher | University of British Columbia |
Date Issued | 2016 |
Description | Imaging and prediction of fluid flow within the subsurface provides information crucial to decision making processes in fields such as groundwater management and enhanced oil recovery. The flow of a fluid through a reservoir depends primarily on the permeability of the subsurface rock; a quantity that is often unknown throughout the entire domain of the reservoir. One method for predicting flow is to estimate the permeability of the reservoir and simulate flow through a mathematical subsurface flow model. Given the model, flow data can be inverted to estimate the permeability. However, this inversion approach can lead to inaccurate results due to the sparse sampling of flow data, and thus inaccurate predictions. To acquire a higher sampling of data, geophysical survey techniques are applied in order to efficiently collect a higher density of data sampled at the surface. These data are sensitive to changes to the geophysical properties of the reservoir due to flow. Inversion of geophysical data then provides images of changes to the geophysical properties of the reservoir. In order to estimate the flow parameters using geophysical data, the two mathematical models require coupling. The thesis therefore proposes two approaches to improve the imaging and prediction of flow. First, a novel coupled inverse problem for estimating the fluid velocity field and the initial geophysical property model from geophysical data is developed. Second, a new method of optimally designing the geophysical survey for the coupled inverse problem is developed. The new adaptive design approach builds on traditional A-Optimal design methods such that historic data are included in the design algorithm. This produces designs that adapt with flow in the subsurface and reduce the collection of unnecessary data. Both the coupled inverse problem and adaptive survey design method are demonstrated using a seismic tomography geophysical survey and a tracer advection fluid flow model. Numerical examples show that the coupled approach yields an improved flow estimate as well as improved image quality, while the adaptive optimal designs provide sufficient geophysical data. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2016-05-11 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0300491 |
URI | http://hdl.handle.net/2429/58109 |
Degree |
Doctor of Philosophy - PhD |
Program |
Geophysics |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2016-09 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
Aggregated Source Repository | DSpace |
Download
- Media
- 24-ubc_2016_september_fohring_jennifer.pdf [ 4.28MB ]
- Metadata
- JSON: 24-1.0300491.json
- JSON-LD: 24-1.0300491-ld.json
- RDF/XML (Pretty): 24-1.0300491-rdf.xml
- RDF/JSON: 24-1.0300491-rdf.json
- Turtle: 24-1.0300491-turtle.txt
- N-Triples: 24-1.0300491-rdf-ntriples.txt
- Original Record: 24-1.0300491-source.json
- Full Text
- 24-1.0300491-fulltext.txt
- Citation
- 24-1.0300491.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0300491/manifest