Reservoir characterization based upon the onset of time-lapse amplitude changes

Time-lapse geophysical monitoring has potential as a tool for reservoir characterization, that is, for determining reservoir properties such as permeability. Onset times, the calendar times at which geophysical observations begin to deviate from their initial or background values, provide a useful basis for such characterization. We found that, in contrast to time-lapse amplitude changes, onset times were not sensitive to the exact method used to related changes in fluid saturation to changes in seismic velocities. As a consequence of this, we found that an inversion for effective permeability based upon onset times was robust with respect to variations in the rock-physics model. In particular, inversions of synthetic onset times calculated using Voigt and Reuss averaging techniques, but inverted using sensitivities from Hill ’ s averaging method, resulted in almost identical misfit reductions and similar permeability models. All solutions based on onset times recovered the large-scale, resolvable features of the reference model. Synthetic tests indicated that reliable onset times can be obtained from noisy seismic amplitudes. Testing also indicated that large-scale permeability variations can be recovered even if we used onset times from seismic surveys that were spaced as much as 300 days apart.


INTRODUCTION
Geophysics provides many methods that are sensitive to the effects of fluid flow in the subsurface.For example, there is a long history of time-lapse seismic techniques applied to detect fluid movement, in a variety of contexts from crosswell to surface seismic (Tura and Lumley, 1998;Landro, 2001;Hoversten et al., 2003;Calvert, 2005;Hatchell and Bourne, 2005;Hodgson et al., 2007).More recently, electromagnetic methods have been deployed for timelapse monitoring of fluid flow (Karaoulis et al., 2011).Geodetic methods, which typically involve repeated surveys over time, have also been used to detect fluid-induced ground motion (Rucci et al., 2010).Although such surveys can often detect changes in geophysical quantities, such as seismic velocities, electrical conductivity, or volume change, it has proven difficult to relate these changes to quantities that are of direct interest in understanding fluid flow.For example, it is extremely difficult to relate static geophysical properties to permeability.It is possible to relate changes in seismic velocity to changes in fluid saturation and pressure at depth.Such relationships still depend upon physical models of porous rocks and may rely on aspects that are not sufficiently well known, such as the nature of the fluid distribution in the pore space on the scale of centimeters to meters (Pride et al., 2004).
One difficulty is that we are trying to relate the magnitude or the size of changes in geophysical attributes to the degree of saturation or pressure change.Computing the amplitudes of changes can be problematic because they depend upon coupling terms that often vary spatially and are often poorly known, introducing another potential source of nonuniqueness.The elastic compressibility of a rock frame is an example of a coupling term, and its spatial variation is an issue when one is trying to relate bulk volume changes within a reservoir to fluid pressure changes (Rucci et al., 2010).Another difficulty is that geophysical surveys are often spaced so far apart in time that several important aspects of the fluid flow are aliased.For example, pressure, temperature, and saturation effects may be confounded due to inadequate temporal sampling.
In an effort to avoid some of these difficulties, we explore an approach to time-lapse reservoir monitoring and characterization that is not based upon the magnitude of changes in geophysical observations.Rather, our approach is rooted in the concept of an onset time, the calendar time at which an observed geophysical quantity begins to deviate from its initial or background value.The onset time should not be confused with a seismic traveltime nor with induced time shifts in seismic traces.The use of onset times requires sufficient temporal sampling to identify the initiation of a change in a measured quantity.For example, monthly seismic reflection surveys would provide a temporal resolution in which we could isolate the onset of significant changes in reflection amplitudes to within roughly 30 days.Even-finer resolution is provided by automated data collection systems such as the multilevel continuous activesource seismic monitoring (CASSM) system, which gathers an entire crosswell survey every three to four minutes (Ajo-Franklin et al., 2011).Given the difficulties associated with the aliasing of fluid flow effects, such improved time sampling is certainly desirable.As shown in this paper, an onset time is sensitive to aspects of the fluid flow and to flow properties and much less sensitive to the details of the rock-physics model.
Our approach, based upon onset times, follows a recent generalization of techniques developed for the interpretation of geodetic data as described in Vasco et al. (2014).In that paper, onset times were introduced in the context of seismic first-arrival times and crosswell imaging.The onset of changes in the traveltimes from a borehole source to receivers in an adjacent well were used to understand the migration of carbon dioxide within a reservoir.The movement of the carbon dioxide, in particular its breakthrough at various points within the reservoir, constrained the permeability variation among the boreholes.The technique was applied to CASSM data, gathered at the Frio pilot site near Houston, Texas.
Estimating reservoir properties such as permeability, based on onset times, contrasts with previous efforts that relied upon timelapse amplitude, traveltime, and waveform changes (Huang et al., 1998;Vasco et al., 2004;Vasco, 2004a;Dong and Oliver, 2005;MacBeth and Al-Maskeri, 2006;Dadashpour et al., 2008Dadashpour et al., , 2009;;Feng and Mannseth, 2010;Rey et al., 2012).Here, we demonstrate that onset times, extracted from seismic time-lapse reflection data recorded in a few surveys, can provide a reliable basis for reservoir characterization.The onset time approach complements existing time-lapse amplitude and time-lapse waveform inversion methods, for example, providing a robust initial permeability model.

METHODOLOGY
As noted in the Introduction, the onset time is that instant at which a time-lapse observation begins to deviate significantly from its background, or initial, value.In this section, we present a method for characterizing the flow properties of a reservoir using such onset times.To achieve this, we first relate the onset of a change in a time-lapse attribute to the flow properties of a porous medium.We then present a brief outline of an iterative algorithm for estimating permeability given a collection of onset times.The general approach is similar to the inversion technique described in Vasco et al. (2014).Therefore, we will abbreviate the initial discussion somewhat and focus on the application to the inversion of time-lapse amplitude changes.

Sensitivities: Relating onset times to flow properties
Our description is subdivided into three parts: First, we relate changes in the saturation within a producing reservoir to the flow properties that characterize it.Second, a perturbation in the onset time of a change in saturation is related to a perturbation in the effective permeability at a point in the reservoir.Third, we relate changes in the fluid saturation to seismic time-lapse amplitude changes.In particular, we relate the onset of time-lapse amplitude changes to the initiation of a change in reservoir saturation at that point.

Relating saturation changes to the properties of a porous medium
The governing equations for fluid flow lie at the heart of our approach.Although the methodology is general and applicable to any situation involving transient flow, even multicomponent flow, we will restrict our discussion to the case of two fluid phases.One component will represent an aqueous or water phase, denoted by a subscript w, whereas the other will consists of a nonaqueous phase, signified by a subscript n.We denote the saturation of the aqueous phase by S w ðx; tÞ, whereas the saturation of the nonaqueous phase is represented by S n ðx; tÞ.Because the pore space is completely filled by the two fluids, the saturations must satisfy the relationship S w þ S n ¼ 1.We begin with the equations for the conservation of mass for the nonaqueous phase (Peaceman [1977], p. 16): where q n is the source or sink term, ϕðxÞ is the porosity, and ρ n is the fluid density.An identical equation holds for the aqueous phase, and the phases are coupled by the constraint that the saturations sum to unity.The fluid flow velocity v n is given by two-phase generalizations of Darcy's law.Darcy's law states that the velocities for the aqueous and nonaqueous phases are driven by their respective pressure gradients and gravitational forces.For the nonaqueous phase, we have the relationship where KðxÞ is the absolute permeability, μ n is the fluid viscosity, p n is the fluid pressure in the nonaqueous phase, g is the gravitational acceleration, Z is a unit vector in the direction of the gravitational field, and k rn ðS w Þ is the relative permeabilities, given by the ratio of the effective permeability of each fluid to the absolute permeability K (Peaceman [1977], p. 15).An aqueous phase flow velocity vector v w is defined in a similar fashion.
It is possible to gain some physical insight by writing the system as an equation for the average pressure p avg ¼ ðp w þ p n Þ∕2 and an equation for the saturation of one of the phases (Peaceman [1977], p. 18).The equations are still coupled, but they do simplify under certain circumstances, for example, when the fluids are incompressible.An alternative decomposition into a system of equations for S w and p w is given by Vasco (2011).An equation for the phase of the coupled saturation and pressure front provided expressions for the velocities of the pressureand saturation-dominated disturbances.Typically, the pressure-dominated disturbance propagates much faster than the saturation front.This velocity difference has implications for time-lapse monitoring and reservoir characterization.With sufficient temporal resolution, it is possible to distinguish the pressure-and saturation-dominated disturbances.In most of today's time-lapse monitoring systems, the effects of saturation and pressure changes are aliased in time and cannot be distinguished.Relating perturbations in permeability to variations in the onset of saturation changes Now, let us consider the onset of changes within the reservoir due to production and/or injection.For a given reservoir model, well configuration, and production/injection source terms q n and q w , the governing equations can be solved using a numerical simulation code such as integrated finite differences (Pruess et al., 1999;Pruess, 2004) or streamline-based computation (Datta-Gupta and King, 2007).We shall concentrate on the onset of saturation changes due to the propagation of a two-phase saturation front from an injection well.That is, we will assume that the rapidly varying pressure changes have already migrated through the reservoir, or that the pressure effects can be neglected.
The numerical solution of the governing equations provides a pressure and saturation history for each grid block within the reservoir model.Here, we use a integral finite-difference-based simulator TOUGH2 (Pruess et al., 1999;Pruess, 2004) and the threelayer reservoir model plotted in Figure 1 to calculate the saturation and pressure changes due to the injection of carbon dioxide.Each layer in the model is 30-m thick for a total reservoir thickness of 90 m.In the uppermost layer (1.80-1.83km), higher permeabilities lie to the northwest and due east of the injection well, which is denoted by the black square near the center of the plot.The permeability varies by more than two orders of magnitude, and the lowest permeabilities are to the south of the injector and in the northeast corner of the layer.Each layer is subdivided into a 41 × 41 grid of cells, giving lateral dimensions of 62 × 62 m for each grid block.
The carbon dioxide is injected into the top two layers of the reservoir model.From the reservoir simulation histories, one can compute the onset of a saturation change for each grid block of the model.The onset time of the transition from a water-saturated state to a CO 2 -saturated state is plotted in Figure 2 for all three layers.
Due to buoyancy, the carbon dioxide flows preferentially into the top layer of the reservoir model.The reservoir simulator also provides the saturation phase velocities, for example, v n , the velocity of the nonaqueous phase.Using this information, we can construct flow paths from any given observation point x o that has experienced a saturation change back to the injection well.We denote this path by Xðx o Þ.The propagation time, which will be the onset time referenced to the start of injection, is then given by the line integral for the nonaqueous phase.One can draw out the explicit dependence upon the flow properties ϕðxÞ and KðxÞ, where we have defined the vector U n ¼ ∇p n − ρ n gZ, and we have made use of the definition of v n , given by equation 2. We are assuming that the fluid viscosity μ n is constant and that the relative permeability function k rn ðS n Þ is known for the reservoir under study.As noted by Vasco et al. (2004Vasco et al. ( , 2014)), the vector U n harbors an implicit dependence upon the flow properties.It is shown by Vasco et al. (2004) that, as long as the pressure fields are recalculated for each iteration of the inversion algorithm given below, the implicit dependence may be neglected.
For our iterative inversion algorithm, in which we update a given reservoir model to fit a collection of observations, we shall need to relate perturbations in the flow properties to perturbations in the propagation time τ.Assuming that the fluid properties, the relative permeability, and the capillary pressure properties are fixed, perturbations may occur in the porosity, the permeability, and, because of its implicit dependence on porosity and permeability, the nonaqueous phase velocity vector U n .As we recompute U n at each iteration of our inversion algorithm, we neglect the effect of its perturbation upon the propagation time.As shown by Vasco et al. (2014), under these restrictions, one may relate perturbations in τ to perturbations in the permeability via the integral Equation 5 provides a semianalytic relationship between a perturbation in the reservoir permeability along the path Xðx o Þ and a perturbation of the arrival time, or the onset time, of the saturation change at the observation point x o .The quantities in this expression are easily calculated using the results of a single numerical reservoir simulation.

Relating the onset of saturation changes to the onset of timelapse seismic amplitude changes
The actual data will consist of time-lapse geophysical observations that are indirectly related to the state of the reservoir.That is, we must rely upon a rock-physics model to connect saturation changes to changes in seismic properties.In the example below, the synthetic time-lapse amplitude changes are related to the fluid saturation changes through a velocity-saturation relationship.The reflection coefficients for the reservoir layer depend upon the seismic velocity V, which we write as a function of the spatial position, time, and the onset time of saturation change at the given location, given by equation 4. We include the onset time to underscore the connection between the arrival time of the fluid phase and the onset of a corresponding change in the geophysical properties.Consider the velocity of the compressional wave, given by Vðx; t; τÞ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi K s ðx; t; τÞ þ 4 3 G ρ s ðx; t; τÞ s (6) (Aki and Richards, 1980;Chapman, 2004), which is a function of the bulk modulus of the saturated rock K s ðx; t; τÞ, the shear modulus G, which we assume does not depend upon the fluid saturations, and the density of the fluid-saturated rock ρ s ðx; t; τÞ.The density of the saturated rock is given by a weighted linear average of the densities of the constituents where and ρ i is the density of the ith fluid under reservoir conditions.We will adopt Gassmann's (1951) relation for the bulk modulus of the fluid saturated rock: where K d is the dry frame bulk modulus, K g is the bulk modulus of the solid grains comprising the rock, and K f is the fluid bulk modulus.
We are now faced with the task of deciding which model to use for mapping the fluid saturations, S n and S w , into a composite fluid bulk modulus K f .Currently, there is no universally accepted method for such a mapping for an arbitrary distribution of fluids within the reservoir.In fact, there is reason to believe that the relationship between a variation in fluid saturation and the change in the effective bulk moduli depends upon the distribution of fluids within the porous medium on a scale smaller than a seismic wavelength (White, 1975;Dutta and Odé, 1979;Norris, 1993;Knight et al., 1998;Johnson, 2001;Pride et al., 2004).Thus, it may be difficult to determine an appropriate model to use if the medium to smallscale distribution of the fluids must be known a priori.
In the face of such ambiguity, we will consider upper and lower bounds on the moduli.Specifically, we will examine how varying between maximum and minimum values of the moduli influences estimates of flow properties obtained by the inversion of time-lapse amplitude changes and onset times.This approach will be discussed in more detail in the "Application" section below.For now, let us consider the particular bounds that will be used.The Voigt upper bound for the bulk moduli of a mixture of N fluids is given by    where S i ðx; t; τÞ is the saturation of the ith fluid and K i is the bulk modulus of the ith fluid (Mavko et al., 1998).The upper bound is obtained when the components are arranged in parallel to the direction of compression (propagation).The stiffest layer then controls the modulus.The Reuss lower bound for the effective fluid bulk modulus is given by This lower bound is obtained when the compression is in a direction perpendicular to layers of pure components.Then, the weakest layer determines the modulus.Such bounds are well established and used in a variety of fields, such as computing bounds on the effective thermal conductivity of a composite material (Wiener, 1910).
Other bounds, such as the Hashin-Shtrikman bounds (Hashin and Shtrikman, 1963), are possible, but the Voigt-Reuss bounds are the most conservative for a particular volumetric mixture of fluids.
In addition to upper and lower bounds, we consider the Hill estimate of the fluid bulk modulus (Hill, 1963).The Hill estimate is just the average of the Voigt and Reuss bounds: In Figure 3, we present the Voigt and Reuss upper and lower bounds and the Hill average as a function of the saturation of carbon dioxide.The bounds indicate that differences in velocity estimates due to the choice of the averaging method are larger than the maximum velocity change due to the fluid substitution.For example, at a carbon dioxide saturation of 10%, the difference between the Voigt and Reuss bounds is approximately 0.20 km∕s whereas the total change in the Hill estimate is around 0.15 km∕s.Given an elastic model for the overburden and a poroelastic model for the reservoir, such as the 11-layer model shown in Figure 4, we can compute the seismic response as a function of the changing state of the reservoir.For example, the impedance contrast at the top of the reservoir will be related to the velocity and density changes within the reservoir layer, given by equations 6 and 7, according to where ρ o and V o are the density and compressional velocity associated with the layer overlying the reservoir.One could use expression 13 to approximate the reflection amplitude changes associated with saturation changes within the reservoir.We will use the approach adopted by Vasco et al. (2004), which includes reverberations as well as reflections in the response.In particular, we use the approach of Kennett (1974Kennett ( , 1983) ) to compute the seismic response of a stack of elastic layers.Because we only consider near-offset reflection amplitudes and assume that the reservoir model consists of cells of significant lateral dimensions, we represent the reservoir and the surrounding layers as a collection of 1D columns.Each column may have a distinct set of fluid saturations within the reservoir interval, governed by the flow due to injection and production.
In Figure 5, we plot the calculated amplitude changes for reflections from the interface at the top of the reservoir model.The injected fluid flows preferentially into the high-permeability region of the uppermost layer, as indicated in Figure 2, leading to spatial variations in fluid saturations and corresponding variations in time-lapse amplitude changes.The amplitude changes are associated with saturation changes over the entire 900 days of injection.As indicated above, the velocity changes resulting from the saturation changes will vary, depending upon the small scale distribution of the fluid.For the calculation of amplitude changes plotted in Figure 5, we use the Hill average (equation 12) to compute the fluid bulk moduli.Then, Gassmann's relation and formula 7 for the density are used to compute the velocity changes according to equation 6.The velocities are used, along with the corresponding formulas for the shear velocities and the densities, in Kennett's (1974) method for computing the elastic response and the amplitude changes.The nonzero amplitude changes are restricted to regions in which there are saturation changes induced by the 900 days of injection, as is evident in comparing Figures 2 and 5. Now, consider how variations in the method used to map the fluid saturations into a composite fluid bulk modulus results in different time-lapse amplitude changes for reflections off the top of the reservoir.In Figure 6, we plot the amplitude changes after 900 days of injection, calculated using the Voigt upper bounds (equation 10) and the Reuss lower bounds (equation 11).Note the significant differences in the magnitude and the spatial variation of the amplitude changes.To gain some insight into the amplitude variations over time, let us consider the amplitude changes at two points, labeled A and B in Figure 1.In particular, in Figure 7, we plot the time series of the amplitude changes at points A and B for all surveys conducted over the 900-day interval.The surveys are spaced 30 days apart, simulating monthly time-lapse monitoring.There are clear differences in the magnitudes of the calculated amplitude changes depending upon which averaging technique, Voigt, Reuss, or Hill, is used.One notable feature in Figure 7 is the consistency of the onset times for all three of these approaches.That is, the calendar times at which the amplitudes start to deviate from their background values are similar for the Voigt-and Reuss-based amplitude changes.This agreement is seen for all reflection points (Figure 8), suggesting that the calculated onset times are not sensitive to the method used to estimate the fluid bulk moduli.We should note that there are several ways to define an onset time.For example, one may take the time at which the observable quantity deviates from its background value in a statistically significant manner.Alternatively, as in Vasco et al. (2014), one may define the onset time as the time at which the change in a quantity exceeds a preset percentage of the total change of the observed quantity.Here, we define the onset time as the calendar time at which the rate of change of the observed quantity is the greatest.
The robustness of the onset time with respect to variations in the model used to calculate fluid bulk moduli suggests that it might be useful for characterizing reservoir flow properties.In particular,

Amplitude inversion
Hill Reuss Voigt Figure 9.The squared average misfit to the amplitude changes as a function of the number of iterations of the inversion algorithm.The three curves delineate the different methods (Hill, Reuss, and Voigt) used to compute the fluid bulk moduli that are the basis for the calculation of the synthetic amplitude changes.All of the quantities used in the inversion algorithms follow the Hill (1963) approach of averaging the Reuss and Voigt bounds to compute the predicted amplitude changes and the model parameter sensitivities.

a) b)
Figure 8. Onset times associated with the amplitude changes for reflections from the top of the reservoir.The onset time is defined as the calendar time, in days, since the start of injection, at which the amplitude begins to deviate from its background value.The size and color of each box are proportional to the onset time.Panels (a and b) indicate onset times derived from amplitudes computed using Reuss and Voigt fluid bulk moduli, respectively.
Characterization based on onset times M7 Downloaded 10/22/17 to 166.87.199.141.Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/onset times may be primarily sensitive to saturation changes within the reservoir and not very sensitive to the rock-physics model used to calculate seismic properties.We will explore this prospect in greater detail in the "Applications" section below.First, we formulate an algorithm for estimating reservoir permeability based upon a collection of onset times.

An iterative algorithm for estimating reservoir permeability
In this subsection, we present an iterative scheme for updating a reservoir model to match a collection of onset times.The onset times may be obtained from time-lapse reflection data.For example, we may use the onset of a change in the amplitude of a reflection from the top of and/or from interfaces within a reservoir, if the reflection can be connected to a saturation change within the reservoir.The onset times are then related to flow within the reservoir via the path integral equation 3. We may then invoke relationship equation 5 between a perturbation in the onset time and a perturbation in the permeability along the path.A discrete version of equation 5 provides a linear constraint on the reservoir permeability perturbations for each onset time.Given a set of onset times, we have a corresponding collection of linear equations.We may write the system of equations in matrix-vector form: where the elements of A are given by the discrete form of the integrand in equation 5:  iteration, we seek updates that improve the match to the observed onset times or minimize the squared magnitude of the residual vector.Note that, due to the intrinsic nonuniqueness of the inverse problem, it is often desirable to include additional regularization terms, penalizing attributes of the model such as its spatial roughness or its vector norm (Menke, 1989).Thus, we append a regularization term that minimizes the squared magnitude of the perturbation to the squared misfit: R 2 ðδKÞ ¼ ðd − AδKÞ T ðd − AδKÞ þ W n δK T δK; (16) where W n is a scalar weighting that determines the relative importance of the norm minimization with respect to fitting the observations.The necessary conditions for a minimum determine a linear system of equations, for the vector δK (Noble and Daniel [1977], p. 60).We solve the linear system 17 using the iterative least-squares algorithm of Paige and Saunders (1982).The update for the mth iteration is given by For each iteration, we conduct a full reservoir simulation using the previous model K m−1 to compute the elements of the matrix A and the residuals.

APPLICATIONS
In this section, we describe two sets of synthetic tests and we examine the influence of noisy amplitudes on estimates of onset times.In the first set of tests, we introduce variations in the relationship between the fluid saturation change and the seismic velocity variation.This constitutive law is a component of the rock-physics model and is critical in the interpretation of time-lapse seismic data in terms of reservoir processes.In the second set of tests, we vary the interval among time-lapse seismic surveys from every 15 days up to a maximum of 300 days.We invert the corresponding onset times for reservoir permeability variations in the vicinity of the well.

Variations in the rock-physics model
The robustness of onset times with respect to variations in the relationship of the seismic velocity to the fluid distribution, as indicated in Figures 7 and 8, suggests that they might be better observations on which to base an inversion for flow properties.We can test this conjecture by comparing amplitude inversions with onset time inversions for cases in which the relationship between the fluid saturations and the seismic moduli is incorrect.To do this, we consider the inversion of seismic amplitude changes, using the trajectory-based approach described in Vasco et al. (2004).In this iterative inversion technique, trajectories representing flow lines in the reservoir are used to compute sensitivities relating perturbations in flow properties to perturbations in saturation changes.Then, numerical differencing is used to relate the perturbations in saturations to perturbations in seismic observations, such as seismic amplitude changes.The advantage of this approach is its efficiency.This efficiency derives from the fact that the sensitivities are com- puted in a semianalytic fashion, as in equation 5.The approach also requires one reservoir simulation for each iteration of the algorithm.
To establish what can be achieved under optimal conditions, we first consider an inversion when the rock-physics model is correct.That is, the fluid averaging relationship used in the forward problem, in this case Hill's (1963) average (equation 12), is also used to compute sensitivities for the inverse problem.Our starting model was uniform: Each cell in the reservoir model had the same permeability.The resulting reduction in time-lapse amplitude misfit is shown in Figure 9, given by the curve labeled Hill, as a function of the number of iterations of the linearized inversion algorithm.The squared error is reduced by approximately an order of magnitude in 15 iterations.The final permeability estimates for the uppermost layer of the model are shown in Figure 10a.For the Hill inversion, the large-scale pattern of higher permeabilities, observed in Figure 1, to the northwest and to the east of the injector is recovered.As noted by Vasco et al. (2004), the only reliable estimates are for cells that have experienced saturation changes and a corresponding amplitude change.Because of the preferential flow in the uppermost layer, and the fact that the amplitudes for reflections off the top of the reservoir are not very sensitive to changes in the second and third layers of the model, the estimated permeability changes are essentially zero there.Now consider the more likely scenario in which the rock-physics model is, to some degree, incorrect.In this case, the Hill (1963) average is used for calculating residuals and sensitivities for the inverse problem, whereas the Reuss and Voigt approaches are used in the forward problem, to compute the synthetic amplitude changes.The misfit reduction as a function of the number of iterations for these two inversions are also shown in Figure 9, labeled as Voigt and Reuss.When the Voigt estimates are used to compute the synthetic reference amplitude changes, there is no overall reduction in the misfit.When the Reuss approach is used to compute the synthetic amplitude changes (Reuss), the misfit reduction is significantly less than the optimal case (Hill).The final models, shown in two lower panels of Figure 10, reveal that the inversion of the Reuss-based synthetic amplitudes results in a model (Figure 10b) that is very different from the reference permeability model (Figure 1).The model is dominated by a high-permeability region surrounding the injection well, and there is no indication of generally lower permeabilities to the south.The inversion of the Voigt-based synthetic amplitudes produces a model (Figure 10c) that contains anomalies that are of the opposite sign from the reference model.Now, consider an inversion of the onset times, the calendar time at which the reflection amplitudes start to change from their back-ground values.We plot the onset times, computed using the Reuss and Voigt algorithms, in Figure 8.As we just did for the amplitudes, we first consider an inversion in which the rock-physics model is correct.Specifically, we use Hill's (1963) average estimate (equation 12) for the forward and inverse calculations.Thus, at each step of the iterative algorithm, we solve the linear system: where A Hill signifies that the sensitivities were computed using the Hill estimate (equation 12), whereas d Hill indicates that the forward model used to compute the synthetic onset times was also based upon the Hill model.The resulting squared misfit as a function of the number of linearized iterations of the inversion algorithm is plotted in Figure 11.The squared misfit is reduced by greater than an order-of-magnitude in 10 iterations.
Of more interest is the performance of the iterative inversion algorithm, based upon onset time, in the presence of an incorrect constitutive model relating fluid saturations and seismic velocity.For this, we invert onset times computed using the Reuss and Voigt approaches, shown in Figure 8, using an algorithm with sensitivities based upon Hill's (1963) averaging.Instead of solving equation 19 at each step of the iterative updating algorithm, we solve the linear system in which the synthetic onset time residuals are based upon the Reuss, and Voigt, forward models.The decreases in squared onset time misfit as a function of the number of linearized iterations are also shown in Figure 11.The inversion of both sets of onset times (Reuss and Voigt) result in similar reductions in misfit.
In plotting the permeability estimates, we only show the uppermost layer of the three-layer model.Because the sensitivities in the lowest two layers are almost zero, their estimated permeability deviations are dominated by the norm minimization and are near zero.Furthermore, the only significant, resolvable permeability anomalies lie within the area in the layer that has experienced observable saturation changes during the 900 days of injection.This area is indicated by arrival times of 900 days or fewer in Figure 2a.The resulting permeability models for all three inversions (Figure 12) contain the high values to the northwest of the injection well.Figure 13.Time series for the amplitudes of reflections from the top of the reservoir at points A and B (Figure 1).The solid curve indicates the noise-free case, whereas the unfilled circles indicate amplitudes with 5% added Gaussian noise.
The filled circles denote amplitudes with Gaussian noise with a standard deviation that is 10% of the reflection amplitude.The Hill and the Reuss inversion contain the less pronounced higher permeability feature just to the east of the well.The inversion of the Voigt onset times contains smaller, isolated slightly higher permeability anomalies in the eastern area, separated by slightly lower values of K.All models contain significantly lower permeabilities directly to the south of the injection well and somewhat lower permeabilities to the northeast.

Onset times extracted from noisy seismic amplitude data
Actual time-lapse amplitude data are typically noisy due to the nonrepeatable nature of large 3D seismic surveys.For example, there are difficulties in reoccupying source and receiver sites, changes in near surface propagation due to variations in the water table or in the overlying water column, nonrepeatable sources, and so on.These variations lead to deviations in the amplitude data even when there are no changes within the reservoir.In an effort to examine the effects of amplitude perturbations, we added Gaussian noise to the amplitudes before estimating the onset times.The magnitude of the noise is specified with respect to the amplitude of the background reflection from the top of the reservoir.Thus, we considered noise with standard deviations of 5% and 10% of the reflection amplitudes, illustrated in Figure 13 for reflections from two points at the top of the reservoir.The locations of the two points are indicated by the open circles in Figure 1.The random amplitude deviations in Figure 13 are a significant fraction of the amplitude change due to substitution of carbon dioxide for water.To reduce the effect of the noise on the onset time estimates, we applied a smoothing operator over a five element moving window.In Figure 14, we plot the calculated onset times for cases in which 5% and 10% of Gaussian noise are added to the amplitudes.For comparison, we also include the calculated onset times for noise-free amplitude data.The overall pattern is similar for the three cases; however, there are a few outliers in the onset times when the amplitudes contain 10% added noise.This suggests that a robust estimator, such as an l 1 residual norm (Menke, 1989), might be a better choice for the inversion algorithm.Also, additional regularization such as a roughness penalty (Parker, 1994) might be advisable in the presence of significant noise.

Variations in survey intervals
One of the chief drawbacks associated with the use of onset times is the temporal sampling required to resolve them.The time and expense required for a large number of monitor surveys can be a barrier to their wider use.Therefore, one would like to know if it is possible to extract onset times from surveys that are widely spaced in time.In this subsection, we examine four sets of seismic surveys, in which each set is gathered using a different time interval among the surveys.In particular, we consider repeat times of 15, 60, 120, and 300 days and estimate the onset times, which are shown in Figure 15.The reservoir model and layered elastic overburden are identical to that used in the 30-day case presented above.The Hill model was used to relate changes in saturation produced by the reservoir simulator to changes in the elastic properties and, correspondingly, changes in the seismic reflections off the top of the reservoir.Overall, the onset times form a similar pattern, but one can see the truncation effects in the 120-and 300-day interval estimates in Figure 15.For example, there are obvious steps in the onset time   16.Generally, all estimates contain the largescale features of the model that are seen in the topmost layer in Figure 1: higher permeability to the northwest and lower permeability to the south.One can see the effects of the truncation as smearing of the anomalies along the trajectories between the injection well and the grid blocks of the model.Inclusion of a roughness penalty term in the misfit function 16 might improve these results.

DISCUSSION
The utility of onset times, as demonstrated here, has been noted in other contexts, for example, in the use of surface tilt to estimate flow properties in a fracture zone (Vasco, 2004b).The advantage associated with arrival or onset times over the inversion of amplitudes has also been noted in an inversion of transient surface deformation for reservoir flow properties (Rucci et al., 2010).Finally, the onset of changes in seismic traveltimes, due to the injection of carbon dioxide, was used to estimate lateral variations in reservoir permeability (Vasco et al., 2014).
The use of onset times provides an alternative to the direct inversion of time-lapse amplitude changes, akin to use of the arrival time of a fluid front.Therefore, the approach is particularly sensitive to flow properties such as reservoir porosity and permeability.The fundamental principle underlying the approach is a causal connection between changes in reservoir saturation and pressure and changes in the seismic observations.Therefore, the technique works when the geophysical property governing the observation, such as compressional velocity, has a nonzero sensitivity to the changes within the reservoir.For example, all three of the curves in Figure 3 have a nonzero slope for carbon dioxide saturations that are near zero.The approach will face difficulties if the curves flattened out for a significant range of saturation change.Note that in these circumstances, a direct inversion of time-lapse amplitude changes would also have trouble.We have seen how noisy amplitudes can introduce outliers into the onset time estimates.Thus, a robust estimator, such as an inversion based upon the minimization of the l 1 norm of the residuals, might be advisable.Also, long time intervals among surveys will degrade the onset time estimates and lead to poorer estimates of reservoir permeability.Such long time intervals will also increase the likelihood of aliasing changes due to saturation and changes due to variations in reservoir pressure.In Figure 17, we plot the time history of saturation and reduced pressure for the grid block beneath point B in Figure 1.In this figure, one can see an early variation in pressure due to the start of injection.This effect will be more pronounced closer to the injection well.After that, the pressure changes slowly, until the saturation front arrives and there is a change in reservoir pressure due to capillary effects and the saturation change.Without sufficient temporal sampling, it would not be possible to distinguish the two episodes of relatively rapid pressure change.In favorable circumstances, one could use the early onset of the pressure-related change as well as the onset of the saturation related changes.
The technology for frequent time-lapse monitoring is advancing or has advanced in particular areas of geophysics.New approaches allowing for more frequent seismic monitoring have appeared recently.For example, permanent, life-of-field, and continuous oceanbottom seismic arrays are increasingly common (van Gestel et al., 2008;Berg et al., 2012;Bertrand et al., 2014).These advancements should lead to improved reservoir monitoring and better signal-tonoise (Landro and Skopintseva, 2008).Also, permanently cemented shallow seismic sensors are being tested for frequent monitoring of the injection of carbon dioxide (Bakulin et al., 2012).Continuous active-source monitoring is being developed in a variety of settings, including crosswell and vertical seismic profiling configurations, for the monitoring of fluid movement within the subsurface (Ajo-Franklin et al., 2011;Daley et al., 2011).The monitoring of deformation induced by fluid injection and migration has a long history that includes many techniques with short sampling intervals such as tilt, the global positioning system, and interferometric synthetic aperture radar (InSAR).New technologies for such geodetic monitoring, such as high-precision ocean-bottom pressure sensors and X-band InSAR, are under development.

CONCLUSIONS
The results of this study suggest that the onset time, the time at which the seismic observations begin to deviate from their initial or background values, is a useful attribute for reservoir characterization.Calculated onset times are found to be more stable than predicted time-lapse amplitudes with respect to variations in the rock-physics model used for their computation.In particular, onset times appear to be sensitive to fluid break through times but insensitive to the exact relationship between fluid saturations and seismic velocities.The use of onset times is particularly well suited for data from frequent monitoring surveys and a permanent array of sensors.Furthermore, onset times provide a means of data reduction, reducing a multitude of time-lapse surveys into a single set of onsets.
Figure 2. Distribution of fluid front arrival times, associated with the transition from water-saturated to CO 2 -saturated grid blocks, within the three layers of the reservoir model.

Figure 1 .
Figure 1.Reference permeability model used to generate synthetic amplitude changes due to the injection of carbon dioxide into a central well, denoted by the filled central square.The open circles indicate the locations at which we will extract time series of the reflection amplitudes from each survey.The variations in tone denote changes in the logarithm of permeability.

Figure 3 .
Figure3.The velocity of a compressional wave as a function of the saturation of carbon dioxide.All velocity estimates are computed by Gassmann's approach but using different methods for calculating the composite fluid bulk modulus.In particular, the Voigt upper bounds follow from equation 10, the Reuss lower bounds are computed using equation 11, and the Hill estimates (equation 12) are the average of the two bounds.

Figure 4 .
Figure 4. Elastic model of the reservoir and the surrounding layers.The 90-m-thick reservoir interval is indicated by the parallel horizontal lines.Two models are shown, the baseline velocities (indicated by the open circles) and the monitor survey 900 days later (indicated by the black squares).

Figure 5 .
Figure 5. Amplitude changes, calculated using the Hill model, after 900 days of injection into the central well, indicated by the filled square.The dimensions of each box are proportional to the amplitude changes for a reflection at that point.The color scale indicates the fractional amplitude change.

Figure 7 .Figure 6 .
Figure 7. Time variation of the relative amplitude changes calculated for reflections from two points (points A and B in Figure 1) at the top of the reservoir.The time series denote amplitude changes based upon the Reuss, Voigt, and Hill composite fluid bulk modulus estimates.The onset of the time-lapse amplitude changes is indicated in this figure.

)
Figure11.The sum of the squared onset time residuals plotted as a function of the number of linearized iterations.The labels denote inversions in whichHill's (1963), Reuss's, and Voigt's methods were used to compute the synthetic onset times.In each inversion, Hill's averaging approach was used to calculate the predicted onset times and the sensitivities.

Figure 10 .
Figure 10.Permeability variations in the uppermost layer of the reservoir model, resulting from an inversion of the Hill-based amplitude changes (a), Reuss-based synthetic amplitudes (b), and Voigt-based synthetic amplitudes (c).These permeability estimates should be compared with Figure 1a.

Figure 12 .
Figure 12.Final permeability models from the inversions of the Hill-, Reuss-, and Voigt-based onset times using an algorithm that assumes Hill's (1963) approach is the correct way to related fluid saturations and the composite fluid bulk moduli.(a) Result of inverting Hillbased onset times.(b) Result of inverting Reuss-based onset times.(c) Result of inverting Voigt-based onset times.

Figure 14 .Figure 16 .
Figure 14.Onset times estimated from noise-free amplitude data (a), amplitudes with 5% added Gaussian noise (b), and amplitudes with 10% added Gaussian noise (c).The size and color of the filled boxes indicate the magnitude of the onset time.

Figure 17
Figure 17.Saturation (dashed curve with filled squares) and reduced pressure (solid line with open circles) as a function of time from the start of injection.The pressure is reduced by subtracting the background pressure, and it has units of 0.1 MPa.These two histories are associated with the reservoir grid block beneath point B in Figure 1.
Figure 15.Onset times estimated using timelapse data from surveys with different spacing in time.The size and color of the filled boxes indicate the magnitude of the onset time.estimates based upon surveys that are 300 days apart.Using these test data, we conducted four inversions for reservoir permeability variations in the vicinity of the injection well.The inversions are based upon the iterative algorithm with model norm regularization, with each update determined by solving equation 19.The resulting permeability estimates in the topmost, resolvable layer of the model are shown in Figure M12 Vasco et al.Downloaded 10/22/17 to 166.87.199.141.Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/