Using onset times from frequent seismic surveys to understand fluid flow at the Peace River Field, Canada


 Our limited knowledge of the relationship between changes in the state of an aquifer or reservoir and the corresponding changes in the elastic moduli, that is the rock physics model, hampers the effective use of time-lapse seismic observations for estimating flow properties within the Earth. A central problem is the complicated dependence of the magnitude of time-lapse changes on the saturation, pressure, and temperature changes within an aquifer or reservoir. We describe an inversion methodology for reservoir characterization that uses onset times, the calendar time of the change in seismic attributes, rather than the magnitude of the changes. We find that onset times are much less sensitive than magnitudes to the rock physics model used to relate time-lapse observations to changes in saturation, temperature and fluid pressure. We apply the inversion scheme to observations from daily monitoring of enhanced oil recovery at the Peace River field in Canada. An array of 1492 buried hydrophones record seismic signals from 49 buried sources. Time-shifts for elastic waves traversing the reservoir are extracted from the daily time-lapse cubes. In our analysis 175 images of time-shifts are transformed into a single map of onset times, leading to a substantial reduction in the volume of data. These observations are used in conjunction with bottom hole pressure data to infer the initial conditions prior to the injection, and to update the reservoir permeability model. The combination of a global and local inversion scheme produces a collection of reservoir models that are best described by three clusters. The updated model leads to a nearly 70 percent reduction in seismic data misfit. The final set of solutions successfully predict the observed normalized pressure history during the soak and flow-back into the wells between 82 and 175 days into the cyclic steaming operation.


S U M M A R Y
Our limited knowledge of the relationship between changes in the state of an aquifer or reservoir and the corresponding changes in the elastic moduli, that is the rock physics model, hampers the effective use of time-lapse seismic observations for estimating flow properties within the Earth. A central problem is the complicated dependence of the magnitude of time-lapse changes on the saturation, pressure, and temperature changes within an aquifer or reservoir. We describe an inversion methodology for reservoir characterization that uses onset times, the calendar time of the change in seismic attributes, rather than the magnitude of the changes. We find that onset times are much less sensitive than magnitudes to the rock physics model used to relate time-lapse observations to changes in saturation, temperature and fluid pressure. We apply the inversion scheme to observations from daily monitoring of enhanced oil recovery at the Peace River field in Canada. An array of 1492 buried hydrophones record seismic signals from 49 buried sources. Time-shifts for elastic waves traversing the reservoir are extracted from the daily time-lapse cubes. In our analysis 175 images of time-shifts are transformed into a single map of onset times, leading to a substantial reduction in the volume of data. These observations are used in conjunction with bottom hole pressure data to infer the initial conditions prior to the injection, and to update the reservoir permeability model. The combination of a global and local inversion scheme produces a collection of reservoir models that are best described by three clusters. The updated model leads to a nearly 70 percent reduction in seismic data misfit. The final set of solutions successfully predict the observed normalized pressure history during the soak and flow-back into the wells between 82 and 175 days into the cyclic steaming operation.

I N T RO D U C T I O N
Time-lapse geophysical data, observations gathered from repeated geophysical surveys, are well suited for the monitoring of fluid flow within the Earth (Calvert 2005). As a result, time-lapse seismic data have been used to monitor the injection of carbon dioxide for underground storage (Arts et al. 2000), geothermal energy production, as well as to image fluid saturation and pressure changes in due to oil and gas production (e.g. Eastwood et al. 1994;Johnson et al. 1998;Tura & Lumley 1999;Landro et al. 2001;Behrens et al. 2002). The dynamic nature of time-lapse data, the fact that they are often related to saturation and fluid pressure changes in the reservoir, suggests that they could be used for aquifer and reser-voir characterization, as noted in Landa & Horne (1997), Vasco et al. (2004), and Dadashpour et al. (2008Dadashpour et al. ( , 2009Dadashpour et al. ( , 2010. The major impediment to successful characterization is the indirect relationship between the observations and the state of an aquifer or reservoir. To address this issue, a rock physics model is invoked to map the current state of the reservoir into seismic properties or attributes. This act introduces additional parameters that are necessary to characterize the poroelastic properties of the in situ rock. These parameters are usually not well constrained, determined from a few cores or laboratory measurements. Furthermore, the properties almost always vary spatially, particularly between formations. Thus, the introduction of rock physics parameters presents yet another level of non-uniqueness (Chen & Dickens 2009). This is a 1610 Published by Oxford University Press on behalf of The Royal Astronomical Society 2020. This work is written by (a) US Government employee(s) and is in the public domain in the US. barrier to aquifer and reservoir characterization that can be difficult to overcome.
As pointed out in Vasco et al. (2014Vasco et al. ( , 2015 with sufficient temporal sampling it is possible to adopt an approach that mitigates some of the issues associated with the intervening rock physics model. In particular, it is possible to define onset times, the calendar time at which a geophysical quantity changes from its background value. Given a weak causality requirement, the onset time can often be related to the time at which saturation, fluid pressure and temperature change within an aquifer or reservoir. Thus, the onset time is typically related to the arrival time of a fluid pressure and/or saturation front, and hence to the propagation time of the fluid front, rather than the magnitude of such changes. As a result, onset times are sensitive to flow-related properties and relatively insensitive to the parameters of the rock physics model, as demonstrated in Vasco et al. (2014Vasco et al. ( , 2015. In this paper, we illustrate the advantages of onset times for reservoir characterization by examining their use at a cyclic steam stimulation operation at one well from Pad 31, within the Peace River field (Fig. 1) in Alberta, Canada (Lopez et al. 2015;Przybysz-Jarnut et al. 2016). This is a very complicated setting, with heterogeneity, prior production, and documented changes in pressure, temperature, and saturation. In fact, there were four periods of enhanced oil recovery at Pad 31, which covers the area that we will examine: a pad-wide cyclic steam injection from 2001 to 2011, a horizontal steam drive from 2012 until 2013, a pad-wide top-down steam stimulation starting in 2014 and extending beyond 2016, and a localized cyclic steam injection at just one well pattern (31-08) from August 2015 until February 2016. Fortunately, we have a rich set of seismic monitoring data from a dense surface array to aid in our analysis (Fig. 2). The seismic array gathered daily surveys to monitor the fluid-induced changes. The techniques that we have developed allow for the compression of this multitude of seismic surveys into a single map of onset times, which is used to image the heterogeneity within the reservoir. Pad 31 horizontal production wells (red) and injection wells (green). The area covered by the production wells is 1.5 km by 1.5 km in the north-south and east-west directions, respectively. Also shown are seismic sources (red spheres) and receivers (blue dots) of the seismic monitoring system. Map modified from Shell Canada (2016).

M E T H O D O L O G Y
In this section, we describe our approach for using repeat time-lapse geophysical observations, recorded by a permanently buried seismic system, to monitor fluid flow and to characterize the reservoir. The fact that the region has undergone previous production, necessitates a two-stage approach. First, we estimate the values of a set of global parameters, primarily describing the initial pressure, temperature, and saturation of the reservoir and its large-scale permeability structure, as described by low order basis functions. Secondly, the finer-scale permeability variations are determined from both reservoir production data and time-lapse observations of the traveltimes of seismic waves that propagate through the reservoir.

Governing equations
Here we outline the equations governing the conditions within the reservoir and the changes due to fluid injection and production. Such changes lead to temporal and spatial variations in the seismic properties and we briefly describe Gassmann's (1951) approach for estimating elastic moduli in fluid saturated rock. Difficulties associated with estimating the appropriate effective fluid moduli lead us to the concept of the onset time of a change in a geophysical observable. We end this section with a brief description of the global and local updating schemes that will be used to estimate the reservoir properties.

Multiphase flow and thermal stimulation
The reservoir operations at the Peace River field, involving cyclic steam stimulation to extract very viscous bitumen, are described and modelled using the equations of non-isothermal multicomponent flow (Lopez et al. 2015;Przybysz-Jarnut et al. 2016). The multicomponent mass and energy balance equations may be written succinctly using index notation (Pruess et al. 2011), where the index κ indicates one of the N k fluid components and the (N k + 1)th component signifies the heat that flows within the reservoir. The mass and energy balances in the reservoir are given by where the quantities q κ in eq. (1) represent source or sink terms, often associated with injection or production wells. The dependent variable M κ is a mass accumulation term for the chemical component κ. This term is written in the form of a sum, given in terms of the porosity φ, the saturation S β , the density ρ β and the mass fraction X κ β of the fluid phase β. In the equation describing the energy balance, the heat accumulation term for a multiphase system is given by where ρ r is the grain density of the rock, C r is its specific heat, T is the temperature and u β is the specific internal energy in phase β. The advective flux vector for component κ, F κ , is a sum over all of the fluid phases, given by a multiphase form of Darcy's law, for a fluid phase travelling with the Darcy velocity w β . The absolute permeability k is a particularly important quantity, one of the main factors controlling fluid flow within the reservoir. The relative permeability k rβ is usually determined from laboratory experiments on cores from the main formations of the reservoir. The fluid pressure for phase β, P β , is one of the dependent variables along with the fluid saturation of the phase, S β , and the temperature T . The fluid viscosity μ β is determined from laboratory experiments on a given fluid at the appropriate temperatures and pressures of interest. Finally, g = gz is the gravitational force vector that alters the flow in the presence of fluid density variations. The vector for the heat flux is given by where λ is the thermal conductivity of the formation and h β is the specific enthalpy in phase β.
In the forward problem we are given an aquifer or reservoir model and we solve the governing equations and accompanying equations-of-state, initial and boundary conditions for the evolution of the saturation and pressure. The solution is usually constructed using a numerical fluid flow simulator (Peaceman 1977;Datta-Gupta & King 2007). For a realistic model solving the forward problem requires significant effort and is often computationally intensive because the governing equations are non-linear partial differential equations. Here, we shall tackle the inverse problem, in which we are given observations, both flow-related measurements and geophysical data, and tasked with estimating the characteristics of the aquifer or reservoir. This is typically a much greater challenge than the forward problem, requiring at least an order of magnitude more computation. Next, we develop a relationship between timevarying fluid saturations, pressures and temperatures within the Earth, and changes in the seismic properties at depth.

Relating velocities and elastic moduli to saturation, temperature and pressure changes
It is well known (Tura & Lumley 1999;Landro et al. 2001;Lumley 2001;Calvert 2005) that fluid saturation, pressure and temperature changes within and around an aquifer or reservoir will lead to changes in the elastic moduli of the fluid-filled porous medium and thus change its seismic characteristics. For example, the speed of a compressional wave transiting a saturated porous material, V p , depends upon the saturated bulk modulus, K sat , shear modulus, G f r and the density ρ sat of the fluid-filled rock according to (Mavko et al. 2009) We adopt Gassmann's equations (1951) to model the changes in elastic properties due to variations in fluid saturations, as they are generally accepted and widely used and found to agree with observations at seismic frequencies (Landro et al. 2001;Lumley 2001;Calvert 2005;Foster 2007). In Gassmann's approach the Velocity of a compressional wave as a function of temperature for different gas saturations while the pressure and water saturation are fixed. At temperatures less than 48 • del depends upon linear correlations, as discussed in Barker & Xue (2016). The vertical black line denotes a transition from these linear correlations to an equation-of-state. (c) Velocity of a compressional wave as a function of the water saturation. All velocity estimates are computed using Gassmann's approach but with four different combinations of the Reuss and Voigt models for computing the fluid bulk modulus. The labels indicate the weighted fractions of each of these models.
shear modulus is not influenced by the presence or absence of the fluid. Furthermore, the density of the fluid infiltrated rock is simply the weighted average of the component densities and for a multicomponent fluid the composite fluid density is given by the weighted sum The bulk modulus of the fluid saturated rock has a more complicated dependence on the component properties, given by the function where K f r is the bulk modulus of the porous rock frame, φ is the effective porosity of the medium and the bulk modulus of the mineral K g , which in the simplest case of a consolidated sandstone can be taken to be the bulk modulus of quartz. The parameter K f is the bulk modulus of the pore-filling fluids. The original formulation of Gassmann only considered a single fluid saturating a porous rock. In order to generalize the approach, the fluid modulus, K f , has been extended to cover the case in which the fluid is a mixture of several liquids and possibly gases. This leads to additional complications, with the essential difficulty that the composite modulus can depend upon how the fluids are distributed within the porous medium at length scales that are less than a seismic wavelength. By considering two extreme distributions one can derive upper and lower bounds on the effective fluid bulk modulus for a given fluid saturation, known as the Voigt and Reuss bounds, respectively (Mavko et al. 2009). In Fig. 3, we plot the velocity variation based upon the Voigt and Reuss composite fluid moduli as a function of the water saturation, S w . In a complex geological setting, including oriented fracture systems, it can be difficult to determine which modulus is most representative. One compromise estimate involves taking the average of the two moduli, the so-called Hill average (Fig. 3). The differences in the calculated values of the three models shown in Fig. 3 are almost as large as the entire variation due to the saturation change. We will use these limiting moduli to illustrate variations in rock physics models and how they can impact calculated changes in seismic properties associated with changes in fluid saturations.
In addition, the cyclic steam stimulation process used at the Peace River field also produces coupled changes in pressure and temperature within the reservoir, leading to complicated rock physics models (Das & Batzle 2010;Kato et al. 2010). The model of Barker & Xue (2016) was used to map the saturation, temperature and pressure changes into corresponding variations in elastic properties. The sensitivity of the seismic velocity variations as functions of the gas and water saturations, pressure, and temperature are presented in Fig. 4, showing the difficulty in interpreting velocity changes and hence seismic traveltime and amplitude changes in terms of unique variations in saturation, pressure and temperature within the reservoir. This difficulty is compounded by the fact that this area of the Peace River field has undergone earlier production, including a previous pad-wide cyclic steam stimulation that started in 2001 and lasted until the end of 2011. The cyclic steam injection was followed by a brief implementation of a horizontal steam drive operation from 2012 to the end of 2013. These earlier production efforts resulted in spatially varying temperatures, pressures, and saturations prior to the initiation of the top-down stream drive recovery process that ran from 2014 to the end of 2016, and the subsequent follow-up cyclic stream simulation on a single well pattern (31-08) that we shall analyse. The extreme heterogeneity in the initial conditions of the reservoir is indicated by the seismic amplitude variations in a regional time-lapse survey (Fig. 5) used to diagnose well problems during the cyclic stream stimulation, conducted in March 2009. This legacy seismic reflection survey was conducted prior to the daily seismic monitoring that is the focus of our work. In Fig. 5 one can observe large amplitude anomalies, associated with the appearance of gas that was expelled from the volatized oil as the pressure was reduced around the production wells, denoted by the black lines in     the figure. Thus, among other complexities, there are initial variations in gas saturation, temperature, and fluid pressures to contend with. This fact necessitates making these initial conditions a part of the inverse problem. That is, our workflow will use a global inversion approach as an initial step, in order to estimate the initial reservoir conditions and global properties.

The onset of a time-lapse change and its relationship to reservoir dynamics
The magnitudes of seismic velocity changes are influenced by the nature of the fluid distribution within the reservoir at length-scales that are less than the typical seismic wavelengths. Therefore, it can be difficult to relate changes in the magnitude of seismic velocities to changes in fluid saturation, pressure and temperature in a quantitative sense. To overcome these problems, we use an onset time methodology to relate the time-lapse seismic data to the propagating fluid fronts. We will describe this approach using the Peace River reservoir monitoring program as an illustration (Lopez et al. 2015;Przybysz-Jarnut et al. 2015, 2016. A permanent seismic reservoir monitoring system was installed at the field consisting of 49 buried sources, in a rough grid with 200-220 m spacing, at a depth of 25 m (Fig. 1). The 1492 receivers (hydrophones) are situated in a denser grid with 40 m spacing, in 20 m deep boreholes and packed in bentonite. The sources consisted of a set of 37.6 s long single frequency sweeps from 0.4 to 216 Hz. The entire set of 540 sweeps took 6 hr to complete for a single survey.
A time-lapse monitoring program was applied to a top-down steam drive oil recovery process that began in 2014, in which six new horizontal steam injection wells were drilled and operated above existing production wells drilled for an earlier cyclic stream stimulation (Lopez et al. 2015;Przybysz-Jarnut et al. 2016). The data are acquired in a continuous fashion and automatically processed to generate four complete time-lapse data cubes every 24 hr (Lopez et al. 2015;Przybysz-Jarnut et al. 2016). These four cubes are stacked to produce a single daily estimate. A vertical time Downloaded from https://academic.oup.com/gji/article/223/3/1610/5896957 by Lawrence Berkeley National Laboratory user on 24 November 2020 Figure 9. Interpretation of time-shift anomalies due to enhanced oil production at the Peace River field. The area covers the same region as that shown in Fig. 8. The black triangle denotes the location of the set of wells (31-08) analysed in this study. This map shows the cumulative time-shifts in milliseconds since the start of top-down injection until August 2015. The north-south trends of positive time-shift, indicating slow down, are due to the increase in pressure and temperature associated with the overlying stream injectors. The blue anomaly indicating speed up is thought to be due to water breakthrough and filling the area containing gas, seen in Fig. 5. Figure 10. Normalized bottom hole pressure variation during a steam injection and production cycle that was conducted at the isolated well pattern 31-08 at the southern edge of pad 31. The production in this figure extends from August 2015 until mid-January 2016. The peak pressure attained in the area was 7.5 MPa and the minimum pressure was slightly below 2.5 MPa.
section through one such data cube is shown in Fig. 6, along with a density log from a well that is intersected by the cross-section. The pink curve is the top of the reservoir (Bluesky formation) while the blue line indicates the base (Debolt formation). Small but visible traveltime-shifts, for reflections from layers at the bottom of the reservoir, are evident in the two snapshots plotted in Fig. 7. The reservoir appears to be thicker than the dominant wavelength of the seismic traces so tuning effects (Ghaderi & Landro 2009;Zhang & Castagna 2011), which occur when the top and bottom reflections from the reservoir interfere, are probably not an issue.
The daily monitoring allowed for the systematic extraction of small traveltime and amplitude changes for reservoir monitoring. Traveltime-shifts were extracted from the migrated time-lapse cubes using a cross-correlation technique over a 120 ms window that extends beyond the bottom of the reservoir. A triangular-weighting filter was applied to remove edge effects in the cross-correlation estimates. An example of the time-shifts generated by the top-down stream drive, gathered between 14 April 2014 and 30 March 2015, are shown in Fig. 8. There are clear coherent anomalies within the area of interest, generally positive time-shifts are colocated with the overlying steam injection wells, and a large negative time-shift anomaly is associated with the production wells just south of the centre of the well pad. In addition to the time-shifts, we also plot the time-lapse amplitude changes during this time interval. Note the small amplitude decreases associated with the injection wells. However, there are much larger amplitude increases that correlate with the large negative time-shifts. These amplitude increases, and the negative time-shifts, are thought to be due to water from the condensed steam replacing gas that had been generated during the earlier cyclic steam injection. The areas containing this gas are indicated by amplitude anomalies in the legacy seismic survey from 2009 that is plotted in Fig. 5. The time-shifts associated with water encroaching on this region of accumulated gas are plotted in Fig. 9, where we note the major contributing factors to the traveltime-shifts: fluid substitution, temperature and fluid pressure changes.
The traveltime-shifts are sensitive to velocity changes and possibly deformation within the reservoir itself. In the manner of seismic tomography, the time-shifts of waves propagating through the reservoir are a sum of the changes within each grid block of the reservoir model. Thus, if we consider a restricted segment of a seismic wave, propagating from a reference point just above the reservoir to the base of the reservoir, and then reflecting from the base of the reservoir and returning to the reference point, the total traveltime-shift is given by where B(x, y) denotes the indices of the grid blocks that are traversed by the seismic wave that is observed at location x, y of the seismic array, L n is the propagation length within the specified grid block and V p (x, y, τ, n) is the seismic velocity within the grid block, V o (x, y, n) is the baseline velocity. For two surveys that are closely spaced in time we are assuming that the ray paths do not change significantly in eq. (10). As noted above, the velocity is time-dependent due to the changes in fluid saturation, pressure and temperature induced by the injection and production. For our analysis of onset times we shall focus on an even later redevelopment on the southern-most portion of the production pad, where a single well set (31-08) underwent an additional cyclic steam stimulation (CSS). In this process steam was injected for 82 days, allowed to soak in and heat up the viscous oil, and then pumped out along with the mobilized oil. The normalized pressure response in the well, associated with one complete cycle which lasted for 175 d, is shown in Fig. 10. The seismic data is translated into transit timeshift maps, expressing the traveltime changes for the seismic waves that propagate across the reservoir between a chosen baseline survey (e.g. the start of the cycle) and subsequent monitor surveys. Over the stimulation cycle shown in Fig. 10, a total of 175 time lapse seismic surveys were available for integration (Fig. 11). Note the temporal and spatial complexity of the time-shifts around the two wells, as shown in Fig. 11. The interpretation of the time-shifts is based upon a rock physics model, using expressions (6), (9) and (10) for T (x, y, τ ) given above, coupled with the relationship between the seismic velocity V p (x, y, τ, n) and the changes in fluid saturation, pressure and temperature provided by Gassmann's equation and the other rock physics expressions from Barker & Xue (2016) noted above. As an illustration of the complex nature of the time-shifts, consider the temporal variations in the size of the traveltime-shift at a location near the injection well 31-08 (Fig. 12). The complicated spatial variation in the pattern of time-shifts is indicative of the various physical mechanisms at play in the field. For example, pressure induced velocity changes can arrive at a location much faster than saturation changes (Vasco 2011). Similar considerations apply to thermal fronts which can take even longer to move through a porous medium (Vasco 2010). Such transient fluid fronts are usually aliased by conventional time-lapse surveys, that are most often taken years apart, but they can be reliably imaged by a daily monitoring program. The fact that the magnitudes of the recorded time-shift data combine several processes involving pressure change, thermal effects and saturation variations, makes it extremely challenging to incorporate them directly into a history matching procedure. For this reason, we utilize the onset time idea to integrate the seismic time-shifts into a reservoir characterization scheme.
The onset time is defined as the calendar time at which the traveltime-shift exceeds a chosen threshold value. The first step is to define a threshold value that results in a meaningful definition of an onset time, as illustrated in Fig. 12. This pre-defined threshold has two main roles: (1) to ensure that the magnitude of the seismic observation is above the background noise level (2) to define the physical process that is being tracked, which often decides the sign of the threshold value. Time-lapse seismic data are typically noisy due to non-repeatable environmental noise, source and sensor issues, and changes in near surface propagation due to variations in the water table or in the overlying water column. These variations lead to changes in the seismic characteristics even when there are no dynamic changes within the reservoir, thus we need a threshold value that distinguishes between the noise and a meaningful signal. Based upon the calculated signal-to-noise ratio for data from the array at Peace River, the threshold was defined as a time-shift decrease of 0.1 ms (Fig. 12). To cross-validate the threshold value, we compared the signal with those from locations that are far from the well and where no changes are expected within the reservoir, as shown in Fig. 13.
The use of onset times not only leads to a significant data reduction, collapsing the 175 daily time-shift maps into a single spatial distribution of onset times, but has also been found to be less sensitive to the rock physics model used to interpret the seismic data (Vasco et al. 2014(Vasco et al. , 2015. As a demonstration of this, four different rock physics models were generated by linearly averaging the Reuss and Voigt estimates of the fluid modulus (Figs 3 and 4) to calculate the P-wave velocity. The variations of four such models with water saturation (S w ) are shown in Fig. 4(c), where we observe that the P-wave velocity is very sensitive to the method used to average the fluid moduli. In Fig. 14, we plot the size of the time-shift changes over the injection period (e.g. the first 82 surveys), calculated using the four models. In Fig. 15, we generate the corresponding onset time maps for the four rock physics models. There are no notable differences between the calculated onset times (Fig. 15) for the different models and all models display areal propagation of changes, related in this particular case to steam/fluid propagation. The similarity of the onset times stands in sharp contrast to the patterns of the magnitude of the traveltime-shifts which are strongly influenced by the particular rock physics model used for the calculations (Fig. 14).

Inversion strategy
Due to the coupled nature of our inverse problem, involving both fluid flow and seismic wave propagation, and the complicated processes and initial conditions, we adopt a two-stage inversion procedure. In the first step we conduct a sensitivity analysis in order to find the most important factor influencing our observations. Secondly, we implement an efficient parametrization for both the initial  conditions and properties that is based upon an eigenvalue decomposition of the grid Laplacian matrix. Third, we determine both large-scale properties and initial conditions that are necessary for the fluid flow simulation and the calculation of the seismic timeshifts using an evolution algorithm followed by a cluster analysis of the final population. In the last step we adjust the individual grid-block permeabilities in the reservoir model using an efficient tomographic-like approach to match the onset times. Because the focus of this paper is on matching the onset times, we discuss this step in some detail. Our description of the first steps of the inversion procedure is rather brief, with more detail provided in Appendices A and B. Furthermore, an in-depth discussion of the inversion approach is also given in Hetz (2017) and in Hetz et al. (2017b).

Initial determination of the global parameters
The coupled flow model contains a large number of parameters that need to be specified in order to conduct a numerical simulation. Some properties will be more important than others in controlling the simulation results. In order to discern those parameters that are to be included in the initial global inversion we conducted a sensitivity analysis as described in Hetz (2017). For the sensitivity study, the objective function was defined as the summation of misfits in the onset time seismic response, OT i , and the bottom hole pressure (BHP): By perturbing each parameter and examining the changes in the misfit we constructed a tornado diagram that indicates the relative importance of each major class of parameters in the misfit functional. Based on the sensitivity analysis we found that all of parameters have some influence on the objective function, but the completion interval is the most important parameter indicating the need to adjust the size of stimulated zone. Other important parameters include the permeability and the initial gas saturation.
In order to successfully simulate fluid flow and seismic wave propagation in the reservoir, we need to specify the initial state of the reservoir, including the pressure, temperature and saturation fields, and the large-scale properties of the model. A key element of this first step is a judicious representation of the fields in the initial model in order to maintain flexibility and prevent a proliferation of model parameters. In Appendix A, we discuss a representation in terms of the eigenvectors of the Laplacian of the simulation grid (Bhark et al. 2011). That is, we represent the model x as a linear combination of M Laplacian eigenvectors where φ i are the weighting factors that are to be found in the inversion. This parametrization has the advantage that it is tied to the fluid flow simulation grid, which may be quite irregular in order to represent a complicated geological model. Furthermore, the representation provides a flexible parametrization that can describe a uniform model, a layered model, and a fully 3-D model, and all models in between these end-members. The lowest order basis functions are constants for each layer, while the second set of functions are composed of linear variations within the given layer. The higher order basis functions contain increasingly rapid spatial variations in properties. The weighted summation of the first ten basis functions gives the fields of initial properties. The updating scheme for the global parameters is described in Appendix B. It is based upon the general notion of a set of Pareto optimal solutions (Lobato & Steffen 2017). We adopt this approach in order to treat the inverse problem as a multi-objective optimization task. That is, we are given two primary classes of observations, namely time-lapse seismic data and bottom hole pressure measurements, leading to two distinct misfit functions, given by (B1) and (B2) in Appendix B. We wish to determine models that minimize the misfit to the N s observed onset times (OT) and N b bottom hole pressures (BP) given by respectively. To some degree the set of Pareto optimal solutions generalizes the notion of a trade-off curve in geophysical linear inverse theory (Menke 2018). In particular, Pareto optimal solutions cannot be improved with respect to a given objective function, such as the fit to the seismic onset times, without increasing the value of at least one of the other objective functions. As noted in Appendix B Pareto optimal solutions lie on the boundary of the set of feasible solutions, the Pareto front. We generate the set of feasible solutions using a stochastic evolutionary technique, the genetic algorithm (Park et al. 2015). In this approach we represent a model in terms of binary strings. A randomly generated collection of models progressively evolves from one generation to another by mutation (random changes) and recombination (joining of portions of the models). The misfit functions M s (x) and M b (x) contribute to the definition of a fitness function that is used to select the models that are retained in the succeeding generation. The Pareto optimal solutions are defined with respect to the population of models in a generation. This set of solutions is further subdivided into groups of solutions using a clustering algorithm (Appendix B).

Local updates of the solution clusters
The second major step in the inversion algorithm involves adjusting the clusters of solutions through an iterative updating scheme. The entire process takes place on a fine-scale reservoir model that may consist of tens of thousands to millions of grid blocks. Therefore, efficiency is a paramount consideration. To this end we adopt a semi-analytic, streamline-based technique for calculating model parameter sensitivities, first presented in Vasco et al. (2004). The general idea, as it relates to the onset of changes in the time-shifts, is that the injected fluids or transient pressure fronts propagate outward from the source well to various points within the reservoir. For the cyclic steam stimulation associated with the wells in the pattern 31-08 in the Peace River field, we will be concerned with injected steam that may condense into water, and associated pressure and temperature changes. The changes in the elastic moduli resulting from the arriving fluid fronts lead to changes in the seismic waves propagating through the reservoir and alter the traveltimes of these waves. In the absence of significant deformation, the onset of a change in the seismic traveltime is directly related to the arrival time of the fluid front. For a 3-D model we can compute trajectories from each grid block where the saturation has changed to a point on the injection well by streamline simulation (Datta-Gupta & King 2007). We can use time-of-flight or streamline simulation methods to relate the arrival or onset time to reservoir properties along the flow paths or streamlines (Vasco et al. 1999(Vasco et al. , 2004Rey et al. 2012;Vasco & Datta-Gupta 2016;Watanabe et al. 2017). As an illustration, consider the movement of a thermal front due to the injection of steam or hot water along the streamline trajectories shown in Fig. 16. The traveltime for the injected steam, after condensing to hot water, τ (r o , t), from a point on the injector, r i to a location in the reservoir where we observe a change in a geophysical observation, r o , is given by an integral along the flow path where q w (r, t) is the velocity vector of the hot water at the leading edge of the coupled fluid front. This vector follows from the form of the flux vector in eq. (4) and is given by (Vasco & Datta-Gupta 2016) where and z is a unit vector pointing in the downward direction. The quantity κ(S w , S g , S o , P, T ) is the total fluid mobility given by that is usually considered to vary by formation but is most often taken as constant in a given formation. Note that the total fluid mobility will depend upon the reservoir conditions, through its dependence on the formation relative permeability curves and the   (Peaceman 1977;Vasco & Datta-Gupta 2016) with respect to the water saturation S w , which is also a function of the reservoir conditions and usually specified for each formation or lithology. The velocity of the thermal front is also a function of the porosity φ(r) and the absolute permeability k(r). Finally, the front propagation is controlled by the pressure field that is established during the injection. As mentioned above, the transient behaviour of the pressure field can be rapid in comparison to the propagation time of the saturation or thermal front. Therefore, we shall assume that, after the pressure transients have decayed, the average fluid pressure is primarily a function of spatial position r and will calculate it using a numerical reservoir simulator for a given initial or background reservoir model. In each iterative step we seek local, or grid-block, updates to the permeability model that further refine the fits to the seismic onset times and the bottom hole pressure data. In computing model parameter sensitivities, we fix the relative permeability functions and the capillary pressure curves for the formation, using values obtained from the initial geological data and the global update, as well as the initial saturation, pressure, and temperature conditions of the reservoir and the large-scale porosity and permeability variations. The sensitivities, relating a perturbation in the permeability at a location in the reservoir to a deviation in the onset time are obtained from the path integral for saturation front traveltime, τ (r o , t). That is, substituting the perturbed absolute permeability k = k o + δk into the expression for q w (r, t) and then into the integral gives where q o signifies the fluid velocity in the background or current reservoir model. The semi-analytic expression for the sensitivity of the onset time is given by and provides the basis for an efficient, tomographic approach to refining the local permeability model using onset times (Vasco et al. 2014(Vasco et al. , 2015. For each column of cells, the trajectory that represents the first arrival of the injected steam to a grid block in the column, is the path that determines the onset time for that location. Fig. 16 shows the correlation between the time-shift onset time and the saturation and the time-of-flight in days for a neutral tracer injected with the water. This traveltime is proportional to the propagation time of the injected water. In order to map the time-of-flight of a neutral tracer to the traveltime of the injected water we must multiply by the derivative of the fractional flow curve and the total mobility, as indicate above. The main purpose of Fig. 16 is to illustrate some of the trajectories that are the basis for the semi-analytic sensitivities, given by eqs (17) and (18), that for the basis for an efficient local inversion algorithm. Using a reservoir model, we may discretize the integral for the perturbed onset time associated with trajectory of the lth streamline, τ l , into a sum over the segment in each grid block of the reservoir model traversed by the path: The set of paths to each point in the model where we have estimated an onset time leads to a system of equations, δτ = Mδk, that may be solved in a least squares sense. That is, we minimize the sum of the squares of the residuals The conditions for an extremum of R 2 , the vanishing of the gradient with respect to the model parameters leads to the system of equations that may be solved for δk. The system of equations could be illposed if there are effectively fewer equations than unknowns. The usual remedy is to introduce additional regularization requirements, such as specifying that the magnitude of the model updates remains small if it is not constrained by the data, and, because the data cannot resolve small features, the spatial variations of the updates are often assumed to be smooth (Menke 2018). Such considerations, encapsulated in quadratic penalty terms lead to an augmented system of equations, as discussed in Vasco & Datta-Gupta (2016, p. 212). This matrix M is sparse because individual trajectories only intersects a small percentage of the grid blocks in the model. Therefore, the system of eq. (21) is solved using a least squares QR algorithm (LSQR) designed for large, sparse linear systems (Paige & Saunders 1982). In order to solve the nonlinear inverse problem, we iteratively update the model, adding perturbations and then recompute the residuals and quantities used in the linearized inversion, such as the saturation, pressure, and temperature fields. After a sufficient number of iterative updates, the misfit tends to level off, and the algorithm is terminated. In the section below, we illustrate the application of this approach to the data from the Peace River field.

A P P L I C AT I O N T O T I M E -L A P S E M O N I T O R I N G AT P E A C E R I V E R
The methodology was applied to the southernmost region of Pad-31 in the Peace River field, focusing on well set 31-08, three wells (31-8E1, 31-8E2 and 31-8W1) forming a 'tuning fork pattern' shown at the bottom of Fig. 4. One positive feature of the area around Pad-31 was that it lacked some of the vertical heterogeneity seen in other parts of the field. In particular, it did not contain shale baffles that had complicated the vertical flow in many other areas of the Peace River field. Our analysis of the monitoring data from the Peace River field begins with the initial global history match in which we determine the initial, temperatures, pressures and saturations as well as the large-scale variations in permeability and porosity. The inversion is based upon the initial portion of the cyclic steam stimulation involving the injection of hot steam, the first 82 days of the cycle. We use observations from the final soak and flow back to the well for validation purposes, attempting to predict the bottom hole pressure during this process using the history matched models. The initial water and gas saturations, porosity and permeability were taken from a geological model provided by the operator, and the initial temperatures were obtained by interpolating the observed tubing head temperatures at the beginning of the cycle. The reservoir simulation model consisted of an irregular grid with 21 layers with variable boundaries.
The model representation of the global properties is in terms of the eigenvectors of the grid Laplacian matrix, the adjacency-based parametrization described above. A total of ten basis functions, eigenvectors of the Laplacian, were used in the representation of the porosity, permeability, initial water and gas saturations, and the initial temperature. The genetic algorithm used to approximate the Pareto front and to determine the initial set of global parameters for the first step of our inversion scheme ran over 30 generations with population of 150 members per generation. The initial 150 models were generated stochastically by uniformly sampling from expected intervals of parameter values. The values of the model misfit functions associated with the seismic onset times, M s (x), and the reservoir bottom hole pressure data, M b (x), are plotted in Fig. 17(a). The initial scatter in the models, due to sampling randomly from the expected ranges of the parameters provides an indication of the variation in the two misfits expected in the model space for the range of all possible models. After 30 generations the genetic algorithm has reduced the misfit to both the seismic onset time observations and the bottom hole pressure data significantly in comparison to the prior cloud of solutions. The resulting suite of 150 models appear to define a trade-off curve between the two misfit functions, the Pareto front (Fig. 17b). An application of the K-means cluster analysis algorithm generates three clusters that are colour-coded in Fig. 17(b).
By applying a cluster analysis we further investigate the objective space. In particular, Fig. 18 shows the updated onset time maps of selected models in cluster 1, cluster 2 and cluster 3, respectively. For all clusters, we observed some improvement from the initial onset time map calculated using the prior model. The improvements in the match to the bottom hole pressure data are shown in Fig. 19, where we plot the calculated values for 40 models. One notable feature is the consistent pressure match during the soak validation interval, where we used the history matched models to predict the pressure behaviour, indicating that the models are able to adequately represent the saturation changes within the reservoir.
By looking at the parameter changes after the global update, as shown in Fig. 20, we can gain some insight on the different physical mechanisms that are associated with the clusters. For example, we observe that cluster 1 contains the greatest permeability decrease around the well. Also, the water saturation at the base of the reservoir increase more in cluster 1 as compared to clusters 2 and 3. This may explain the overestimation of the well pressure associated with cluster 1 (Fig. 19). Furthermore, the change in the temperature and gas saturation around the well in clusters 2 and 3 indicates different spatial flow patterns for these two models. These differences are reflected in the onset time maps as an underestimation of the propagation time.
The next stage of the inversion workflow involves adjustments to the reservoir permeabilities on the fine-scale grid in order to match the onset time observations for the first 82 d of steam injection and the bottom hole pressure. We apply the iterative linearized inversion algorithm to three candidate models which were selected based upon the cluster analysis. In Fig. 21 we plot the normalized misfit as a function of the number of iterations of the algorithm. The misfit is reduced to almost 30 per cent of its original value. Our iterative linearized algorithm is rather simple and uses a fixed step length for each iteration. The convergence is influenced by the weighting of the regularization and the characteristics of the linear solver that is applied at each step of the iteration. The updated onset time responses from the local step significantly improves the results due to the individual grid-block adjustments of reservoir flow properties. The changes made to permeability field, shown in Fig. 22, reveal that models from both clusters share common characteristics with similar large-scale increases and decreases. These updates imply that the stimulated zones are located mostly around the vertical part of the well. Fig. 23 displays the improvement in the pressure match and prediction as a result of the local updates for clusters 1 and 2 in Figs 23(a) and (b), respectively. Most notably, after the local update the excess pressures associated with first cluster from the global   update are reduced to values much closer to the observed pressures (Fig. 23a).
The final reservoir model produced by the inversion methodology is not only a useful tool for better matching the observations, but also gives additional insight into the state of the reservoir during the cyclic steaming operation. Fig. 24, a plot of the water saturation changes over the injection period, shows that the distribution of water is much less dispersed in the final clusters than it is in the initial model. The final models also help us to identify steam override during production, a common phenomenon in steam injection processes. The reason for this phenomenon is that mobility of displaced fluid is much lower than that of the displacing fluid (steam). Due to the differences in density between steam and the oil and water, steam override occurs. Fig. 25 shows the water saturation along the streamlines over the injection period. At the beginning of the cycle (Fig. 25a) the steam starts moving upward as soon as it is injected inside the model. This movement is captured by the onset time map. The gravity override phenomenon becomes less severe over time, as the fluid starts to move downward at later times (Fig. 25b). Overall our hierarchical history matching approach significantly reduces the misfit associated with the time-varying seismic and pressure data, and provides an improved representation of reservoir sweep through the identification of limits on the distribution of water and the detection of steam override.

D I S C U S S I O N
The use of onset times should be viewed as the first step in the construction of a detailed reservoir model, whereby flow properties are obtained from geophysical time-lapse data. Because onset times are chiefly sensitive to the flow properties of a reservoir or aquifer, and much less sensitive to the parameters of the rock physics model, they are well suited for estimating hydraulic conductivity or permeability. Furthermore, the onset times are related to the traveltimes of fluid fronts, which have a quasi-linear relationship to properties such as hydraulic conductivity (He et al. 2006). As a result, inversions of onset times for properties such as permeability are much less sensitive to the initial or starting model, and inversion algorithms based upon them are much less prone to becoming trapped in a local minimum, similar to seismic traveltime tomographic imaging. The next step would be to use the magnitudes of the time-shifts and reflection amplitudes to further refine the model and to estimate the poroelastic properties of the rock physics model. The final step would be to combine all of the data to construct the final reservoir model.
Like most surface seismic monitoring efforts our study was hampered by issues related to vertical resolution, due to the averaging of seismic waves and their dominantly vertical propagation. It may be possible to improve the resolution by including broad-band data, larger offsets and utilizing the full seismic waveform. Another option would be to use the pre-stack data directly for a tomographic  estimation of the time-shifts or the velocity changes. These enhancements should be topics for future research, as should be development of automated systems for seismic monitoring such as the continuous active seismic source monitoring system (Ajo-Franklin et al. 2011;Vasco et al. 2014). Such systems augment existing permanent arrays for monitoring reservoirs that exist in various fields around the world. While daily monitoring was possible at the Peace River field, the onset time approach is applicable to surveys that are separated by much longer intervals, such as yearly repeats (Vasco et al. 2015).

C O N C L U S I O N
This study demonstrates the advantages of onset times, the recorded times at which a set of time-lapse geophysical data begin to deviate from their initial or background values, for high resolution reservoir characterization. A synthetic test shows that, in comparison to seismic time-shift magnitudes, the onset times are insensitive to the details of the rock physics model used to relate the state of the reservoir to the seismic moduli. The methodology allows for the compression of multiple seismic surveys into a single map of onset times, that are directly related to fluid front propagation times. The compression of the frequent seismic surveys into a single set of onsets assists in the development of an efficient globally convergent stochastic inversion technique, in this case the genetic algorithm.
The Peace River field case treated here displays all of the complexity that one can encounter in enhanced oil recovery, including temperature and pressure variations, saturation changes and complicated reservoir initial conditions. Using a hierarchical workflow, we were able to construct a set of initial models satisfying both the onset times and the well pressure data. The Pareto surface defines a set of feasible solutions, generalizing the concept of a trade-off curve used in linear inverse problems. Using local model updates, where the flow properties were adjusted on a cell-by-cell basis, the algorithm was able to improve upon the global stochastic solution. The final set of reservoir models not only match the data used in the inversion, they also successfully predict well pressure data left aside for validation. Finally, the reservoir models provide insight into the processes operating in the reservoir during the cyclic steaming operation. In particular, the models predict a much sharper water/steam front and reveal steam override due to the influence of gravity. Finally, the estimates of the initial conditions and local permeabilities, allow us to construct an improved injectivity profile along the horizontal well, which is crucial for further development considerations.

A C K N O W L E D G E M E N T S
The work of G. Hetz and A. Datta-Gupta was supported by the Department of Energy under Award Number DE-FE0031625 and Downloaded from https://academic.oup.com/gji/article/223/3/1610/5896957 by Lawrence Berkeley National Laboratory user on 24 November 2020

A P P E N D I X A : M O D E L R E P R E S E N TAT I O N
Because the inversion approach contains what we are calling global and local model updates, essentially large-scale and fine-scale spatial variations of in the properties of the model, we need a flexible model representation that allows for a seamless transition between spatial scales. In this Appendix we briefly describe one such parametrization, as we will incorporate it into our two stage inversion scheme. We shall define the grid by its set of vertices V and edges E, characterizing it by a graph G = (V,E). The N vertices, V = {1,2,. . . ,N}, represent the centre of the grid cells at which the reservoir properties are defined. The edges E represent connections between vertices and one can specify the set of edges using an adjacency matrix a i j , where the non-zero entries denote a connection between vertices v i and v j . Specifically, the entries of the N × N grid adjacency matrix A are given by because we are only considering unweighted connections between vertices. Jafarpour & McLaughlin (2009) showed that a low dimensional approximation may be given by the lowest frequency Fourier components. In order to extend this approach to an irregular mesh, we make use of the association, first noted by Taubin (1995), between the discrete Fourier transform of a function and the decomposition of the function into a linear combination of the eigenvectors of the Laplacian of the grid. The grid Laplacian is a discrete second-order differencing operator given by where d i is the degree of the ith vertex a measure of the number of edges connected to the vertex. The Laplacian provides a measure of the connectivity of the grid and for many commonly encountered boundary conditions the discrete operator is a positive semi-definite, symmetric matrix (Bhark et al. 2011). Given these properties we may use the spectral theorem to construct an eigen-decomposition of the Laplacian matrix where the vectors v i are pairwise orthogonal unit eigenvectors. The eigenvalues λ i are the modal frequencies associated with the Laplacian eigenvectors, a direct consequence of the equivalence between the Laplacian eigenvectors and the basis set of the discrete Fourier transform (Taubin 1995). Here, we will represent the model x as a linear combination of basis vectors that consist of the Laplacian eigenvectors where M is small for a large-scale global specification of properties and equal to N for a full-scale representation of the model. This representation is referred to as an adjacency-based transformation or parametrization. The low frequencies, or small values of M, can be used to represent the global properties of the model, such as a uniform layer velocity, while the highest frequencies account for much more rapid local variations in properties. By changing the value of M used in our representation we can switch between inversions for local and global properties.

A P P E N D I X B : D E T E R M I N AT I O N O F T H E G L O B A L PA R A M E T E R S
In this Appendix, we describe our approach for determining the global properties of the model, including the initial saturations, pressure and temperature of the reservoir, as well as the large-scale porosity and permeability values at the beginning of the stimulation cycle. As a first step, a sensitivity analysis is conducted in order to identify the parameters to be considered in the global updating scheme. A description of that effort is presented in Hetz et al. (2017a,b) and Hetz (2017) and will not be repeated here. We consider the calibration or inversion procedure to be a multi-objective optimization problem. There are two classes of observations, geophysical measurements and hydrological or reservoir engineering where B P o i are the observed bottom hole pressure and B P c i are the calculated bottom hole pressures for the ith observation point. We can linearly combine the misfit functionals to produce a composite measure of sum of the squared residuals. However, it can be a challenge to correctly weight the two classes of data in order to produce a meaningful model. The conventional approach in geophysics is to construct a trade-off curve through a series of inversions, and to pick a point that balances the fits to each data class (Menke 2018). While this technique is useful linear inverse problems, it can encounter difficulties for nonlinear inverse problems, such as our inversion for reservoir properties.
One alternative to the minimization of a composite misfit is to consider multi-objective optimization techniques characterizing the trade-off between different objective functions. A general approach is provided by the notion of Pareto optimal solutions (Lobato & Steffen 2017). These are solutions that cannot be improved with respect to any particular objective function without degrading at least one of the other objective functions. To describe such solutions, consider a multi-objective optimization problem formulated as the minimization of a vector of m objective functions where X is the set of feasible solutions. One may also characterize Pareto optimal models using the notion of solution dominance. A feasible solution x 1 is said to Pareto dominate another feasible solution x 2 if for all indices i ∈ {1, 2, . . . , m} and for at least one index j ∈ {1, 2, . . . , m}. A solution is called Pareto optimal if there does not exist another solution that dominates it. The set of optimal solutions constitutes the Pareto front or boundary and characterize the trade-off between the various objective functions. A class of stochastically driven techniques, known as evolutionary algorithms, provide an a means of generating a Pareto front (Deb 2001). The genetic algorithm, perhaps the most widely used of these techniques, was motivated by an analogy with biological evolution. In particular, an initial set of models is constructed using a random number generator. The parameters describing each model are converted to binary strings, the full description of each model is referred to as a genome or chromosome. The family of models is successively updated by recombination and mutation. Recombination involves taking selected pairs of individuals and forming new members by randomly combining various segments from the two models. The new model will thus be a hybrid model with characteristics of both parent models. In addition, the process of mutation introduces random changes into the genomes of some subset of the new models. The evolution of the population of models is governed by a fitness function of the form exp[− f (x i )] where f (x i ) is the objective function. That is, the probability of selecting a particular model to take part in the construction of the next generation is given by a function of the general form once a new generation of models is produced it is used in the next iteration of the algorithm. The process is repeated until the overall fitness of the population reaches a satisfactory level of some maximum number of generations has been produced. One issue associated with this approach is that it can fail to adequately define non-convex Pareto fronts such as those associated with non-linear inverse problems. We use a stochastically driven technique to address this problem by: (i) Assigning fitness to population members based on nondominated sorting and ranking.
(ii) Preserving diversity among solutions on the same front by examining the distance between solutions (Deb et al. 2002). The models are first sorted according to their dominance rank (Fig. 10).
That is, solutions that are not dominated by any other models with respect to the given objective functions, those lying on or closest to the Pareto front are considered Rank 1 or belonging to Front 1. A model of Rank 2, or lying in Front 2, is only dominated by those of Rank 1 and no others. Generally, a model in Front k + 1: (1) Should be dominated by at least one model in Front k; (2) May or may not dominate solutions in Front k + 2 (Park et al. 2013).
Rather than use the expression P(x n ) given above, the fitness is equal to the rank of the model. If two models have equal rank than the model with the larger crowding distance [cdist i in Fig. B1] is selected to take part in the construction of the next generation through cross-over and mutation.
Finally, the K-means clustering algorithm (James et al. 2017, p. 386) is used to define clusters of solutions that share similar characteristics. We start by assuming that the solutions may be grouped in some number, say K, of clusters. This number may change as we try to find the optimal set of clusters. The main goal is to partition the data set into internally homogeneous and externally distinct groups. The idea is to minimize the within-cluster different between solutions, W (C k ), usually defined by the sum of the square of the distance between each solution in cluster k.
where x k is the cluster centroid and |C k | denotes the number of solutions in cluster k. The approach is initialized by randomly assigning solutions to one of the clusters and computing the centroid of the K clusters. Then the following steps are repeated until the within-cluster distances and cluster assignments stop changing: (1) reassign the observations to the centroid which lies closest to that solution (2) after the reassignment, recompute the cluster centroids. This approach is guaranteed to decrease the measure of total withincluster distances as explained in James et al. (2017, p. 402). The clusters provide an initial set of solutions which we can update in order to match the seismic onset times and bottom hole pressure observations, as described next.