A review and analysis of errors in post‐stack time‐shift interpretation

This study complements two earlier papers on the interpretation and estimation of post‐stack time‐shifts that detail the most popular measurement methods to date. In this current work, the focus is on describing the magnitude and identifying the origin of the uncertainty present in these time‐shift values. We also consider how the underlying assumptions behind conventional time‐shift estimation methods can contribute to the inadequate resolution of the subsurface time‐lapse changes. The various errors fall into three broad categories: (1) those related to intrinsic data limits that cannot be avoided; (2) those associated with seismic measurement methods that can be corrected with some effort; and finally (3) those that arise due to approximations to the wave physics made during the design or implementation of the methods. In the first category are limitations due to sampling rate and signal‐to‐noise ratio, and wavelet interferences resulting from the narrow band nature of the seismic data and the heterogeneous nature of the geology itself. The effects produce errors of typically a fraction of a millisecond, but exceptionally in 4D data with poor repeatability up to 1 ms. In the second category are acquisition and processing errors. These are usually of the order of a few milliseconds but can reach up to 10s of milliseconds and are, therefore, essential to correct. It is possible to correct the effects of tides and seawater variations in marine acquisition if these changes are monitored and measured. Navigation and timing are identified as issues to consider carefully. For land data, daily and seasonal near‐surface variations are still problematic, and source and sensors must be buried to deliver interpretable time‐shifts. The impact of choices made during processing can be significant but is specific to the workflow and dataset and thus cannot be generically assessed. However, the effects from residual multiples can be identified and treated. Of moderate importance is the third category of errors, which consists of four items. Of these, the effect of lateral image shifts and amplitude effects are judged to play a minor role. Deviation from the assumption of vertical ray‐paths during post‐stack analysis appears to be of concern only for reservoirs with noticeably dipping structures and strong contrasts. The effect of pre‐stack variations on the post‐stack measurements remains a topic to be more thoroughly examined, and the conclusion is unclear. Finally, it is apparent that uncertainties in all categories are strongly dependent on the field setting and geographical location.


INTRODUCTION
We focus on the estimation of post-stack time-lapse timeshifts between one or several repeated seismic surveys.These time-shifts provide an essential correction for the 4D seismic difference volumes to provide reliable amplitudes for interpretation.The time-shifts themselves have also proven to be a highly informative seismic attribute that can be directly related to the underlying velocity perturbations and/or subsurface displacements due to reservoir production and the development of stress or fluid migration in the overburden.In Paper 1 of our review series (MacBeth et al., 2019), we discussed the magnitude of time-lapse time-shifts, and how they may be used to interpret dynamic changes in the reservoir and overburden.In Paper 2 (MacBeth et al., 2020), we focused on describing and cataloguing various methods of measurement and discussed their relative performance.In the data examples used to illustrate these two reviews, together with those published in the literature and observed by us extensively in-house, time-shift measurements clearly contain some degree of error and are seldom the clean, trustworthy results we would wish.These errors originate not just from the mathematical approximations and assumptions underlying the chosen technique and its numerical implementation but are also due to the act of measurement on noisy bandlimited seismic data.They are also a product of the type of acquisition and processing workflow.The purpose of this new review, combined with some analyses, is to describe and quantify the factors that contribute to these errors in practice, with a view to guiding robust interpretation.Where possible, solutions will be offered to account for the error source.In other cases, the best that is possible is a measure of confidence in the results that can be propagated along the interpretation workflow.As before, our study is supported by not only several time-lapse seismic examples from producing fields but also some modelling examples.
In the sections below, we describe the major errors currently identifiable in the literature as being critical to the success of time-shift analysis, also when they might occur and their likely magnitude.Our list may not be fully exhaustive but hopes to broadly cover the spectrum of known errors to date.Time-shift errors are separated into three categories: (1) those associated with the fundamental limits of the performance of the methods themselves; (2) acquisition and processing-related errors; and (3) errors due to the inaccurate implementation or approximations of the underlying wave physics.First, we observe that the time-shift resolution must be limited by the signal-to-noise ratio and the bandwidth of the data.We also expect that the idealized assumptions made by the measurement methods (as detailed by the measurement model defined in MacBeth et al., 2020) may be challenged in geology, where wave physics is naturally more complex.Next, we consider errors in post-stack time-shifts caused by the seismic acquisition and induced during the pre-stack processing workflow.As time-shifts are typically only several milliseconds in magnitude, it is not unexpected that they may be easily swamped by phenomena, such as water velocity variations, tidal variations, clock drift and navigational errors.We would also anticipate that trim statics, de-multiple, de-noise/de-stripe, non-flat events pre-stack or wavelet variations to bias the character of time-shifts.This latter group of processing errors are entirely workflow dependent, will vary with company and project and cannot be discussed in detail, although some qualitative conclusions can be reached.Finally, we consider the impact of neglecting the full wave physics on the time-shifts.This includes lateral image shifts, non-vertical ray-paths, pre-stack variations and amplitude effects.As we shall see, in the majority of cases, these effects can be moderate to small.

INTRINSIC LIMITS ON ACCURACY
Within this category, we discuss the factors that govern the ultimate resolution of time-shift methods for noisy, bandlimited seismic data.Consideration is given to trace sampling and detection limits as defined by the signal processing literature.We also investigate the measurement and the accuracy of detecting changes in realistic thin-bedded geology.

Sampling and signal-to-noise ratio
A fundamentally important question is: What is the smallest possible time-shift that can be measured from our data?We start by making a rough, general assessment of this limit by gathering statistics from the review of MacBeth et al. (2019).Table 1 shows the smallest measurable time-shift values extracted from the database of MacBeth et al. (2019), both for off-and onshore acquisitions, plotted with the normalized root mean square (NRMS) non-repeatability metric.The frailties of this metric as a measure of '4D noise' are well documented by Kragh and Christie (2002) and others, and this measure must therefore be taken as a general reference only.In offshore data, overburden slowdowns as small as 1.5 ms have been observed in the Curlew D field (Fehmers et al., 2007), 1 ms on the Skua field (Staples et al., 2007), 1 ms on the Cinguvu field (Reiser et al., 2018) and 0.8 ms on the Elgin field (Grandi et al., 2010).Reservoir or overburden speedups of 1.0 ms have been observed in the Kristin field (Hansen et al., 2009), 1.0 ms on Snorre (Røste & Ke, 2017), 0.4 ms on the Elgin field (Grandi et al., 2010), 0.3 ms on the Ekofisk Life of Field (LoFS) seismic data (Wong, 2017) and 0.5 ms on the Terra Nova field (Andersen et al., 2011).Brain (2017) and Brain et al. (2017) quoted time-shift detectability in the Southern Gas Basin at 0.5 ms for an NRMS of 18% in towed streamer data.It is also mentioned that the interpretation of sub-millisecond time-shifts is now routine (Brain, 2017).Wong (2017) carried out a noise floor study on Ekofisk towed streamer data (with an NRMS of 12%) and LoFS (with an NRMS of 3%-5%) data and found values of 0.3 and 0.08 ms, respectively.Collectively, these findings suggest a lower noise floor of, at most, several fractions of a millisecond for towed streamer data, whereas it appears that permanent seafloor sensors may return values with an accuracy of slightly lower than 0.1 ms.In onshore data, La Follet et al. (2015) reported slowdowns of 0.45 ms and speedups of 1.18 ms for monitoring steam injection into heavy oil using buried sources and receivers in Peace River with NRMS values of 9%.Burying the sensors in boreholes can yield 20% NRMS and detect time-shifts of 1 ms or lower from the injection of CO 2 into reservoirs in Saudi Arabia (Bakulin et al., 2018).Other small time-shift examples for land surveys are: Balol, in situ combustion (2 ms, Zadeh et al., 2010); Rainbow Lake, solvent injection (2 ms; Ng et al., 2005); and Cranfield, CO 2 injection (1 ms, Vasco et al., 2019).All these acquisitions have an NRMS greater than 50%.These suggest the lower limit on the noise floor for land data may be roughly 1 ms.Most land data have a high NRMS, and the non-repeatability errors usually do swamp the time-shifts.All seismic signals in commercial datasets are typically sampled every 2 or 4 ms (500 or 250 samples per second).This being above the Nyquist rate for the seismic data, the Shannon interpolation equation from signal processing literature tells us that time-delay estimates can achieve an accuracy approaching that of continuous data (Walker & Trahey, 1995).Laying aside digitization error, which can be considered negligible in today's applications, error at subsample accuracy is introduced because most time-shift algorithms provide estimates made by using finite-length windows (or trace segments) on data that have been corrupted by noise.This error places a fundamental limit on the ultimate performance of a time-shift estimation method.A theoretical lower bound exists that estimates the fundamental limit of all unbiased time-delay estimators.This limit is known as the Cramér-Rao lower bound (see, e.g., Carter, 1987;Silver et al., 2007;Walker & Trahey, 1995).Following the development by Walker and Trahey (1995), the lower bound (expressed as the standard deviation σ Δτ ), for the time-shift due to uncorrelated noise and a finite window length is (1) where ν = 2Βf 0 W is the number of degrees of freedom (Shannon number) of the data, and ( 2  2 0 + 12 2 0 ) is the second moment (up to a constant) of the uniform distribution with centre f 0 and width Bf 0 , which arises naturally from the assumption of a rectangular spectrum.f 0 is the centre frequency of the seismic signals (baseline and monitor traces), W is the trace window length, B is the fractional bandwidth (bandwidth divided by its central frequency), SNR 1 and SNR 2 are the signal-to-noise ratios (SNRs) of either baseline or monitor trace (i.e. the individual 3D noise levels defined as the RMS amplitude of the signal divided by the RMS amplitude of the noise) and ρ is the source correlation coefficient.The analytic expression in (1) is only possible by making several assumptions: the noise on each trace is uncorrelated but has identically shaped power spectra, signal and noise terms are uncorrelated and there are flat bandlimited power spectra for both noise and the signal.Importantly, the predicted error is the minimum achievable by any unbiased windowbased time-delay algorithm, including, but not limited to, cross-correlation.In situations where the signals are severely corrupted or there is correlated noise, this bound will clearly underestimate the error floor.In the case of seismic time-shift analysis, we assume SNR 1 = SNR 2 = SNR, and SNR ≫ 1 for both baseline and monitor traces, and additionally, that  ≈ 1 and  ≈ 1. Equation (1) can then be simplified to where  = 1∕ √  0  is a number dependent on the window length relative to central frequency f 0 (e.g. it is 0.5 when W = 200 ms and f 0 = 20 Hz).Thus, for our applications, the error is controlled mainly by the inverse of the SNR on the 3D traces and the central frequency.As an example, calculation from (1) shows that for a 200 ms time window, SNR = 10 and f 0 = 20 Hz that  Δ ≈ 0.4 ms.It is predicted from (1) and (2) that as the window length increases, σ Δτ decreases.In practice, smoothed estimates with long time windows also contain considerable bias.This, of course, assumes that the time-shift is stationary over the entire window and initiates the debate of long window information content versus accuracy (this is beyond the scope of the current discussion).However, the essence of this particular algorithm is that there is greater statistical sampling.We must therefore be careful not to read too much into the predictions of this formula but to use its guidance.1.
Figure 1 shows the error from (1) plotted against the NRMS (calculated from the 3D trace SNR -see the Appendix for details).This predicts the smallest time-shift that can be resolved from our acquisition and data sampling.In the figure, the lowest observable time-shifts from Table 1 are plotted for reference purposes.The figure indicates that our guess at the lower bound is roughly in agreement with the lowest values anticipated from Equation (1).For low NRMS values, the observed values are slightly larger, suggesting perhaps that dedicated 4D surveys do not report their ultimate resolution limit.The assumptions behind (1) should also be evaluated.For example, the equation assumes that the noise on the baseline and monitor traces is perfectly uncorrelated.We know that the assumption of uncorrelated noise is not perfectly encountered in practice (see the Appendix).The invalidity of the noise assumption will result in a reduction of the 4D noise relative to 3D noise and could, in turn, imply that our conclusions from (1) are, in fact, an overestimate of the possible time-shift floor.Clearly, although seismic does not quite satisfy the assumptions for this Cramér-Rao formula, as in Silver et al. (2007), we may use it as a rough guide for our seismic data.Finally, we should note that there will be a difference in response to noise depending on the time-shift method used, and we should contrast window-based methods, like cross-correlation and the correlated leakage method, with windowless methods such as the non-linear waveform inversions (WFM) or dynamic time warping (see MacBeth et al., 2020 for a description of these methods).Although Equation (1) applies to a window-based technique, an equation has been quoted by Rappin et al. (2019) originating from underwater acoustics, which contains no window length and may relate to nonlinear WFM: where H is an empirical scaling factor of 2/3, and B is now the bandwidth in Hz.These equations show not only the benefit of improving bandwidth and/or central frequency in our acquisition and processing but also the deleterious effects of reduced SNR.To deliver a result where this error can be neglected requires an optimal balance of bandwidth, central frequency and SNR.For current published examples from dedicated, high-fidelity 4D seismic datasets, we conclude that a low error has been achieved and can be neglected in these cases.This statement cannot, however, be considered a generality, and for each 4D project, the relations above must be used as a guide to deliver an optimal result from acquisition and processing.MacBeth et al. (2016) were the first to promote the awareness of the importance of geological variations in time-shift estimation.They showed, by modelling based on a field case study, that spurious (non-physical) time-shifts of several milliseconds can be present at, or slightly above, the top reservoir.This time-shift has the appearance of a thin signal of prominent time-shift with opposite polarity to that of the reservoir change (visually appearing as, e.g., a thin layer of slowdown over reservoir speedup due to a waterflood).This anomaly is due to waveform phase rotations and sidelobe interferences from reservoir changes, accompanied by tuning within the geological heterogeneity that cannot be interpreted as a simple time-shift.Our experience is that such spurious effects are also present in time-shifts evaluated through the reservoir section and arise from interferences generated by production-related fluid contact movements or velocity changes.MacBeth et al. (2016) concluded that this error varied across a range of different reservoir thicknesses and overburden/reservoir impedance contrasts and for a variety of different measurement methods.Accurate estimates of the time-shifts at the top reservoir require perfect knowledge of the seismic wavelet for the data, which must be input into a time-shift method that takes care of amplitude effects such as those identified by Williamson et al. (2007).Other possible fixes include recourse to the higher frequency content of the data, amplitude balancing or machine learning (Duan et al., 2020).The phenomenon may also be present at significant geological contrasts in the overburden and could therefore compromise the resolution of strains associated with depleting reservoirs.Another effect that can also output similar unexpected, non-physical time-shifts is the dimming of key reflectors on the monitor seismic below the noise floor (Brain, personal communication).In a further study of the underburden, Jaramillo (2020) generalized the spurious time-shift result by modelling the cumulative effect of internal multiples and wavelet interferences in realistic geology conditioned by wireline logs in interbedded sequences.The resultant spurious fluctuations prevent the time-shifts from being tracked accurately and magnify when time-strains are estimated, making true subsurface changes in velocity and compaction hard to quantify.We believe that although in most cases in the industry, we have been able to ignore the effect to date with success, it is present in most data to some degrees or other.It could be impossible to ignore when interpreting thick, stacked reservoir sequences.

Thin bed interferences
Accurate quantification of time-lapse velocity changes within, and hence time-shifts across, thin geological formations is hindered by the phenomenon of tuning.Although time-dependent phase variations in 3D data may be compensated to enhance geological resolution, these phase corrections are not commonly applied to 4D data, and the phase shifts remain present.To understand the thin layer effects on the seismic, consider at first a thin reservoir layer subject to time-lapse changes.Below tuning, when both layer time thickness and time-lapse time-shifts are small compared to the seismic period, the time-lapse difference signature Δd(t) is governed by the time derivative of the wavelet w(t) multiplied by the time-shift Δt across the layer For the case of 3D tuning at the same layer on the baseline data, the seismic trace d(t) is where T is the two-way time thickness of the thin bed.By inspection, it can be seen that the time-lapse signature is a scaled version of the baseline signature, that is Δd(t)/d(t) = Δt/T, and that for a thin bed, there is a very small time-shift between the baseline and monitor traces.Generalizing this result to a stack of layers and ignoring sidelobe interferences, we may assume that whenever a thin bed is encountered, the 4D time-shift signals drop to zero.For a bed with a thickness greater than the tuning thickness (expressed in time as T > 1/4f 0 where f 0 is the peak frequency of the seismic wavelet), the top and base of the geological layer are clearly separated, and estimation methods accurately output the time-shift, Δt.In practice, for realistic geology, the combination of signatures from many layers above and below tuning will cause these effects to overlap and lead to errors that are complex to interpret.Clearly, even the smallest error will magnify when time-shift is converted into time-strain or velocity change.Figure 2 shows the resultant error and the accuracy in tracking velocity changes based on a modelling exercise for discrete layers of variable thickness.Although the mean time-shift errors are 8%, relative velocity errors derived from time-strain can be over 100% and, therefore, cannot be ignored.The changes within individual thin layers are also not clearly identifiable.Figure 3 shows a possible example of time-shift drop-outs at thin beds.These are detected on the fast cross-correlation approach but are healed by the regularization of the NLI2D algorithm.The phenomena above may be better handled using a scheme that jointly inverts amplitudes and time-shifts.
To consider the above effect in a more realistic context, a high-frequency model is created by inverting baseline (preproduction) data from North Sea field data.The baseline impedance model (Figure 4a) is then perturbed by a velocity change profile estimated by 4D pre-stack tomography (Figure 4d; Izadian & MacBeth, 2021) to create the monitor model (Figure 4b).The seismic reflection response of both models is then computed using convolutional modelling (an example of the baseline is shown in Figure 4c), from which time-shifts are estimated with the NLI WFM method described in MacBeth et al. (2020).For the modelling, a 20 Hz zero-phase Ricker wavelet is chosen to replicate the frequency content of the original North Sea dataset.Time-shifts and their resultant time-strains (hence velocity changes) are now directly compared against the input model values (see Figure 5).The estimated time-shifts reveal an error of 21% just below the top reservoir, with these fluctuations amplifying considerably to over 200% when the time-shifts are converted to time-strain (relative velocity changes).Our North Sea example did not have large contrasts in impedance, but it is expected that in the presence of stratified gas layers, strong contrasting geology and internal multiples, interference effects will be more severe.Two of the examples from MacBeth et al. ( 2019) (in Figures 8a and 11a) show some evidence of spurious effects, and in this current paper, we see evidence in Figures 5, 6 and 15. Figure 6a shows one extreme example of CO 2 injection in the Sleipner field, where a spurious and un-interpretable pattern of large time-shifts is produced for the thick reservoir section (although mapped quantities make sense).This is a consequence of two effects: (a) The anticipated smooth time-shift behaviour is disrupted by wavelet interferences created by CO 2 accumulations under intra-reservoir shales; (b) the imprint of noise, as there is low repeatability (NRMS of 65%) for these towed streamer data.Figure 6b illustrates another class of unexpected anomaly in the North Sea towed streamer data with an NRMS of 11%, which may have a physical explanation as injectites in the near overburden (Marsden, personal communication).Unexpected speedups can be observed above injected gas, whereas in the underburden, the expected slowdown is apparent (Tian, personal communication).

ACQUISITION AND PROCESSING-RELATED ERRORS
A key point raised by MacBeth et al. (2019) is that the accuracy of both post-stack and pre-stack time-shifts is intimately connected to acquisition and processing repeatability, and that the numerous factors affecting repeatability also affect timeshifts (roughly, these are acquisition geometry differences, near-surface conditions, environment, noise and geologyfor more details on each item, see Johnston, 2013).Timeshifts were initially seen as part of the 'cross-equalization' workflow and not necessarily information in their own right (a small time-shift of 2 ms can induce a normalized root mean square (NRMS) as large as 25% in 20 Hz data).How-ever, as 4D seismic surveying has evolved over the past two decades, the quality and repeatability of our post-stack data product have improved considerably (see the examples of Campbell et al., 2011Campbell et al., , 2015;;Smith et al., 2021).One reason is that acquisition systems are now available with more repeatable geometries.Another reason is the careful measurement of the many controllable and avoidable adverse influences on the acquisition, which helps to deliver highquality co-processing solutions in the pre-stack domain that drive down the non-repeatability noise.As the data evolve through the now highly specialized, strongly time-shift-aware 4D processing workflow, repeatability metrics or their equivalent are tracked and successfully used to enhance 4D signal fidelity (Helgerud et al., 2011;Saint-Andres et al., 2012).Our desired time-shifts, created due to differences in subsurface velocity between baseline and monitor, can thus be propagated and isolated in the stacked time-domain data, to be measured independently of the choices made in the pre-stack processing workflow (of course, pre-stack time-shifts can also be measured, but it is not developed in depth in our current discussion).Thus, there is now considerable confidence that the time-shifts measured on the post-stack, migrated data are directly attributable to the underlying reservoir, over-or underburden time-lapse changes due to production.Unlike the previous section describing the possible intrinsic limits, acquisition and processing errors can be controlled by careful and judicious practice.However, inaccurate calibration data or sub-optimal parameter settings will create a residual contribution from the acquisition and processing nonrepeatability that propagates into the final product as an error.Of course, the post-stack 4D data that we previously used with confidence for 4D amplitude interpretation may not contain meaningful time-shift information, suggesting a quality control feedback loop is required to ensure that processing optimizes the stacked product for time-shift analysis (this is usually applied by specialized contractors).To illustrate what might be expected, Figure 7a shows an example of time-shifts estimates from full offset post-stack data generated by two slightly different pre-stack post-migration processing workflows.The reasoning behind perturbing these workflows is described in Hatab and MacBeth (2021).Workflows associated with realizations 1 and 2 differ in their treatment of noise, especially multiples.Subtle differences in the time-shift features are apparent, becoming much more pronounced when time-strains are evaluated (see Figure 7b).
Based on internal projects and a literature search, we identify several groups of effects which contain errors that are clearly present in the post-stack time domain.One wellknown category, that is now less of a problem, is related to source and receiver positioning.These errors may arise during acquisition geometry merger with navigation data, feathering, node planting, rough weather conditions or inaccurate understanding of source and receiver depths.In ocean bottom node (OBN), data source timing due to clock drift may be an issue.Large and apparent uncertainties are also present due to water velocity and tidal variations.Near-surface variability due to weather dependency and heterogeneity plays an important role on land.In the pre-stack processing workflow and post-stack conditioning steps for the seismic data, there are many processes that may potentially harm the time-shifts ultimately estimated on the full-offset, post-stack volume.Preservation of 4D consistency between baseline and monitor cannot always be maintained, due to, for example, trim statics (see Figure 8 for an example -although the intention in this case was not to generate interpretable time-shifts, but to use only the corrected amplitudes for interpretation), demultiple, de-noise components and de-striping to remove the acquisition footprint.Also included in this list are corrections applied during the early pre-processing of sail lines, treatment of offset classes or individual gathers or during migration and post-migration processing.Residual multiple energy may   In practice, service companies provide a robust and sophisticated 4D pre-stack processing that includes careful and considered treatment of time-shifts.It is not possible to comment on these treatments, but they will vary from company to company.Instead, in this section, we focus on describing the primary sources of error highlighted in the literature, which form a part of these solutions.

Tide and seawater velocity variations
For marine towed streamers and permanent or retrievable seabed systems, temporal and spatial variations in the sea surface height and wave velocity in the seawater are of considerable concern in producing accurate and interpretable 4D (and 3D) seismic data.This is because of the time taken to acquire seismic surveys, which can be up to several months, coupled with the dynamic nature of the ocean environment.Apart from areal coverage, the acquisition time of seismic data depends on a broad range of practical operational and logistical considerations.These, in turn, influence the sequence and timing at which the surveys are shot (Omofoma, 2017).Tides vary across the day and change with the season but are fairly constant across the survey area.Water velocity can vary daily depending on the currents, geographical region, water depth and season.Velocity and tide heights will likely differ for sail lines within a particular survey and between subsequent monitor surveys.Tidal and water velocity variations are well known to create statics in towed streamer data (Mackay et al., 2003), permanent reservoir monitoring (PRM) systems (Bertrand & MacBeth, 2005) or OBN (Wang et al., 2013).Seawater velocity variations are particularly problematic in deep water, where a few metres per second fluctuation can give rise to many milliseconds of time-shift.Timing differences vary with reflector depth and recording offset and thus must be taken into account in the pre-migration data before completing the remainder of the processing workflow.The problem has mainly been resolved for primaries (our main concern for production-related time-shifts), but not for multiples that have additional ray-paths through the water layer.
In the data, the effect of tides on time-shifts is visible as obvious striping parallel to the acquisition of sail lines or 'jitter' in crossline sections (Barley, 1999;Lacombe et al., 2006).This arises as surveys are shot daily in batches of sail lines, which sample different environmental conditions.The extent of these variations can be appreciated by visualizing the specific time stamps of the seismic acquisition -see Omofoma (2017), Omofoma et al. (2019) and Figure 9 for illustrative examples.It should be mentioned that in the case of the towed streamer example of Figure 9b, different boat passes (i.e.stripes) are in different directions, which adds to the repeatability problems when acquiring such seismic data.With tidal variations of a few metres, time-shifts may reach 3 ms.The tidal static component can be largely removed, as tide height is readily measured to decimetre accuracy using differential GPS (see Barker et al., 2003) or other means.Recent research suggests that differential GPS remains popular in marine seismic acquisition and research disciplines beyond (Di et al., 2020).However, the seawater velocity change is more problematic -this is caused by well-understood temperature and salinity (and pressure) fluctuations throughout the water column (see, e.g., Medwin, 1975).The variations strongly depend on the weather, the mixing of currents and the time of year (Bertrand & MacBeth, 2003).Environmental conditions can generate significant depth gradients and also lateral variations in water velocity on the scale of the survey size.Temperature largely governs velocity variations and provides its range of variation in the open ocean.With general velocity fluctuations of a few metres per second to 10s of metres per second, time-shifts can be approximately 1 ms for shallow water applications and 14 ms for deeper water (Wang et al., 2013).Therefore, water velocity variations can produce larger time-shift errors than tidal variations (Hatchell et al., 2008).In post-stack data, neglect of water velocity variations gives rise to a long period imprint on the data.The synthetic example in Figure 10 models the expected impact of these water velocity changes on timeshifts from North Sea datasets acquired in summer and winter conditions.
Fortunately, the application of travel-time/velocity analysis of the seismic data itself can be used to provide a good direct estimate of the water velocity (Mackay et al., 2003).In a separate treatment, Lacombe et al. (2006) applied the method of Landrø and Stammeijer (2004) to the water layer by using travel-times of waves reflected from the water bottom to estimate the desired time-shift correction.This technique works for deep-to-moderate water depths of a few hundred metres.This approach is not accurate when targeting the shallow reflectors near the sea floor because of limited near offset data and a poorly recorded water bottom.Therefore, Ong et al. (2015) provided an extension to the method for the case of shallow water marine streamer datasets by including overburden reflectors in addition to the water bottom.The method will not work if there is a dip or 3D variations in the water layer.A more recent method of note is that of Teien (2017), who applied a time-shift inversion to the water layer arrivals to estimate tidal and seawater velocity variations.
Permanent installations or retrievable OBN acquisition in deep water applications experience more severe water column-related problems than towed streamer acquisition.For deep water settings, Wang et al. (2013) demonstrated that the desired corrections can be estimated by directly measuring both the water velocity and tidal height variations during seismic acquisition using specially manufactured seafloor devices.These devices are called Pressure Inverted Echo Sounders (PIES), and they can measure lateral and vertical water velocity profiles.To retrieve the necessary information, the device measures the two-way time of an acoustic signal reflected from the sea surface and combines it with depth data determined from pressure measurements.PIES have been deployed in a series of field trials in the deep waters of the Gulf of Mexico for Ocean Bottom Node surveys.Wang et al. (2013) noted that these devices can also be used for towed streamer data.

Navigation and timing differences
Despite best practices and high acquisition quality, in a statement made in 2013, Johnston concluded that it was not uncommon to find bulk shifts in positioning between seismic surveys.Ten years later, bulk shifts are now less of an issue and harder to detect from other sources of error.However, care must still be taken to prevent these from entering the workflow.If present, acquisition mispositioning in combination with structural dip will exaggerate the time-shift errors considerably.This effect was demonstrated in towed streamer data by Fehmers et al. (2007) for Curlew D, a North Sea gas condensate field at 3.3 km depth.They showed that despite NRMS of 27% for the data, small time-shifts of up to 1.5 ms could be observed from cross-correlation analysis.From their analysis, they conclude that the measured timeshifts are partly due to systematic errors at dipping horizons caused by a relative mispositioning error of only 6 m between the baseline and monitor surveys, exacerbated by steep dips.Once corrected, the time-shifts (now up to 1 ms in magnitude) show a more explicit correlation with production and highlight a depleted fault block.In further studies, Wang and Hatchell (2013) emphasized the importance of receiver repeatability for dedicated 4D OBN data in deep water Gulf of Mexico applications.Lecerf et al. (2011) described a method for ensuring reliable receiver positioning for 4D OBN in deep water in Brazil.Finally, Kiyaschchenko et al. ( 2020) described the importance of corrections for node timing errors, tidal variations and shot coordinate errors for deep water Brazilian 4D OBN applications.
Timing differences between surveys may be problematic when a baseline is not dedicated to 4D studies.Timing differences may also arise due to reference time variations from line to line or shot to shot.For OBN, additional factors of source/receiver positioning uncertainties and/or time drift of nodes collectively make the removal of water velocity effects in time-lapse seismic data processing difficult.According to Ceragioli et al. (2010) and Lecerf et al. (2011), these are the main factors that may impact measured time-shifts in OBN data.After acquisition, travel-times must be reconciled after node recovery through re-synchronization between the GPS times at the shot and internal OBN clock times.Although OBNs are positioned to a high degree of precision, some errors remain in those positions and are expressed by small relative time-shifts from one node to another (several milliseconds).Shot points are also subject to small errors, dependent on the state of the weather at the time of shooting.Thus, repositioning of each node based on first arrivals, residual clock drift correction and residual and fractional static corrections must be computed and applied.Lecerf et al. ( 2011) described a processing-based technique to take care of the many factors influencing time-shifts simultaneously.Kiyashchenko et al. (2020) presented two approaches for deep water 4D OBN applications in which seawater variations and bulk timing errors are decoupled using the travel-time information from the primaries and water bottom multiples (assuming the node depth and tidal factors are fixed by other non-seismic data).

Residual multiples
Multiples are highly non-repeatable due to the water velocity and tidal variations, and/or subsurface variations, but during 4D cross-equalization, only the primaries are aligned (Calvert, 2005a;Hatchell et al., 2008).The seismic data, therefore, always contain a degree of residual multiple energy, and this coherent noise may obfuscate our time-lapse interpretation of the primary events.The problem is more pronounced in some cases as the residual multiple energy can possess significant energy and overlap entirely with the reservoir of interest (van Gestel, 2021).Calvert (2005aCalvert ( , 2005b) developed a method for removing multiple differences from time-lapse seismic data by shooting one (or both) of the baseline and monitor surveys twice at different tidal or two-way water time states and deriving an operator that subsequently removes the multiples on the time-lapse difference data.The Calvert technique is a solution to the multiple removal problem but adds additional acquisition costs because the baseline or monitor must be acquired twice.Hatchell et al. (2008) examined the impact of long-period free-surface-and water column-related multiples on data from offshore PRMs shot over the compacting chalk reservoirs of the Valhall field.They demonstrated via a modelling exercise for post-stack data that multiples create a linear bias on the time-shift estimates and that this correction factor scales by an amount directly proportional to the ratio of the residual multiple and the primary energy.Multiple energy leads to anomalously high time-shifts in the overand underburden and large time-shifts near the top reservoir.Access to multiple surveys as on Valhall, in addition to knowledge of water velocity variations, permitted this relation to be exploited and the desired correction term to be applied to the data.In practice, the time-shift correction for multiples can be of the order of 2 ms (in a total measurement of 2 ms) and is thus a very significant part of the time-shift signal for the monitoring of short production periods.In a separate study for marine towed streamer data, Fehmers (2010) observed that in the underburden of the Tyra Field in the Dan-ish North Sea, a compacting chalk reservoir, there are large, localized time-shifts related to the misalignment of multiples.The misalignment, in this case, is caused by velocity changes in the overburden and reservoir generated by geomechanical changes.Time-shift vanishes below the reservoir and picks up again to the original 7 ms in the deeper underburden.The study focuses on free-surface-related multiples, but the author argues that the phenomenon described is present to a greater or lesser degree for all multiples formed at reflectors between the seafloor and top reservoir.Figure 11 recreates this effect in a synthetic by modelling the 4D seismic traces from the finescale model shown in Figure 4. Residual multiples are added to the 4D signal according to the recipe of Fehmers ( 2010), and time-shifts then re-estimated in the presence of residual multiple amplitudes of 10%, 20% and 30%.The resultant error manifests as a low-frequency vertical modulation overlying the estimated time-shifts.The overburden time-shifts (and hence time-strains) are impacted strongly, and the strength of discrete 4D signals in the reservoir interval and the underburden time-shifts are also affected.In our example, the errors are roughly 0.1-0.5 ms in an overall time-shift profile with a maximum of 3 ms.Residual multiples remain an important and somewhat elusive source of error in time-shifts estimated from post-stack seismic data.
As a footnote, it is now well known that multiples can provide a wealth of information on the subsurface and can be used to improve imaging (Verschuur, 2013) and assist with time-lapse studies (Qu & Verschuur, 2020).Indeed, for deep water 4D OBN, data they may form part of the imaging solution.Thus, for quantitative interpretation, we anticipate that time-shifts derived from multiples will be used as a source of information in the future, but for now, we must treat them as a source of significant error in our time-shift estimates.

Near-surface variability (land data)
On land, it is known historically that there is more difficulty in applying 4D technology with success.Repeatability issues relate to scattering in the near-surface combined with daily and seasonal variations in near-surface rock properties (Hornman et al. 2011) and receiver coupling variations.The near-surface changes also impact the repeatability of ghosts and multiples.In desert conditions, additional considerations arise due to dune movement (shifts of up to 1 m) and a very complex near-surface of varying sand thickness, highly variable near-surface velocities and underlying limestone karsts (see, e.g., Bakulin et al., 2012;Jervis et al., 2012).Such effects spread across the pre-stack volume and contaminate the measurement of subsurface time-shifts estimated poststack unpredictably.It is observed that surface sensors are more sensitive to daily changes in the weathered layer than surface sources, and that by burying sensors in boreholes, the near-surface effects can be reduced.If possible, the ultimate solution on land is to bury both receivers and sources (Berron et al., 2015;Hornmann & Wills, 2011).Even with such ideal acquisitions, near-surface time-shifts of the order of a millisecond have been observed to fluctuate daily in such desert environments (Jervis et al., 2018)

APPROXIMATIONS TO WAVE PHYSICS
A final category of error involves the underlying assumptions of the time-shift algorithms which invariably approximate the true seismic wave propagation to simplify the method for data application.For many datasets and field interpretations, this may be sufficient as an '80% solution', but the user should be aware of the errors involved and how they manifest, to prevent the over-interpretation of potential artefacts.Here we will address four categories of problems that limit the potential of post-stack time-shift methods.

Pre-stack variations
Failure to properly account for TVO (time-shift variation with offset, or more accurately angle: TVA) can introduce bias in time-shifts estimated from post-stack data, preventing true recovery of the subsurface velocity changes via the procedure of vertical differentiation to time-strain.For an overburden with 4D changes that are evenly spread and homogeneous in nature, a natural variation of time-shift with offset is predicted to occur at every subsurface location due simply to increasing ray-path length (Landrø & Stammeijer, 2004).Thus, for example, a 2 km thick overburden with an evenly distributed strain change of 10 −4 is calculated to have a timeshift of 1.20 ms at a ray angle of 0˚, and 1.53 ms at 30˚. (We assume in this calculation an overburden R-factor of 5 and background velocity of 2 km/s.)This is an idealization; however, as in practical operational settings, subsurface time-lapse changes are likely to be unevenly distributed in a background of no change.Long wavelength time-lapse velocity changes close to the reservoir can arise due to the stress arches and pillars induced by reservoir depletion or injection (Hatchell & Bourne, 2005;Røste et al., 2005).Shorter period velocity variations in the overburden can also occur due to, for example, cuttings injection (Narongsirikul, 2019), gas injection (Falahat, 2012), gas blowout (Zadeh & Landrø, 2011) or fault reactivations (Røste et al., 2007).Additionally, time-lapse variations in the reservoir due to pressure and saturation are somewhat mixed in terms of the scale of spatial variation.Røste et al. (2007) described a ray-path sampling of a small-scale fault reactivation that can create significant TVO.The studies of Kudarova et al. (2016) and Dvorak et al. (2018) suggest that offset (and also azimuthal) variations can range from negligible to significant (MacBeth et al., 2018).When the data are stacked, pre-stack time-shift variations due to ray-path sampling of an uneven distribution of time-lapse subsurface changes as described above are also stacked.Therefore, the post-stack time-shift estimates will be a stacked average of the time-shifts present at each offset.The TVO gradient can be used as a proxy metric for post-stack confidence (Hatab, personal communication).Offset-dependent time-shifts also affect the flatness of the monitor gathers, which have been aligned using the baseline velocity model, and this variation is detrimental to the final post-stack sections.The extent of this problem will depend on the noise levels present at each offset and the exact mag-nitude of the TVO at each reflector.This general problem is exaggerated if the baseline traces have not been adequately flattened or if amplitudes vary substantially from baseline to monitor.
The pre-stack variations described above impact the reservoir and overburden time-shifts estimated on post-stack data and hence their conventional (most popular) interpretation in terms of time-strains and velocity changes (as detailed in Mac-Beth et al., 2019).To illustrate this point, Figure 12a,b shows the expected ray-path coverage for an idealized, small-scale rectangular velocity perturbation embedded in a subsurface with no other time-lapse effects.This situation may represent a change in the overburden or a localized pressure or saturation change within the reservoir.It is observed that there is full offset sampling for reflection points inside the anomaly itself, but for other subsurface reflection locations, the ray-paths sparsely sample the anomaly.The exact offset dependence depends on the acquisition geometry, the overburden changes, the common mid point position relative to the anomaly of change and anomaly depths and dimensions.Figure 12c,d shows the expected time-shifts when estimates are made on full offset, stacked data, for a compact and extended reservoir change.In practice, the time-shifts will also depend on the noise levels at each offset, the strength of the TVO at each reflector and which offsets are most representative of the final stack.In contrast to a laterally continuous spread of time-lapse effects (most models of reservoir behaviour in the literature), a more laterally compact velocity change gives rise to significant post-stack time-shifts only inside the anomaly itself.Outside the anomaly boundary, post-stack time-shifts are predicted to decay in magnitude with depth and laterally as the number of time-shifted traces in the stacked gathers diminish against the background noise (Dvorak & MacBeth, 2019).This is because the ray-paths for the medium-to-large offset unevenly sample the anomaly.Compact time-lapse anomalies will not show the vertical tail of constant time-shift that might be anticipated from vertical ray-paths only (MacBeth et al., 2019), but a localized 'blob' of time-shift with a tail that decays with depth.With depth, the near-offsets sample the reservoir change but the far offsets start to miss, hence on stacking the traces in the gather, events with no time-shift are added to events with time-shift -hence, the time-shift of the composite post-stack trace is diluted.The result may be considered a consequence of varying TVO responses combined with noise levels.Conversion to time-strain via vertical differentiation will therefore produce strong artefacts at the blob edges that can be misinterpreted as property changes.Figure 13a-d shows several examples of post-stack time-shifts estimated in field data from producing reservoirs in which a laterally compact region of concentrated change occurs.In each case, the tail of the time-shifts decays with depth in accordance with our understanding.Occasionally, low repeatability combines with ray-path coverage which creates blank patches in the time-shift as in Figure 13d.This pre-stack 'bias' will also affect attributes formed by evaluating the lateral gradient.For data with significant noise on the near offsets, the central core of the blob may also be misplaced or missing, as the nears fail to register the time-shift signal, and medium-to-far offsets cannot focus on the result.
The ideal solution to deal with pre-stack time-shifts is to pose the velocity change estimation problem in the prestack domain.This can be achieved, for example, using pre-stack time-shift tomography (Edgar & Blanchard, 2015;Izadian & MacBeth, 2021), which is best implemented postmigration and post-residual moveout correction (Izadian & MacBeth, 2021).This technique must initially tackle the challenge of estimating time-shifts from pre-stack traces with low signal-to-noise, which is not without complications and some pre-conditioning (Izadian & MacBeth, 2022a, 2022b).Encouraging studies have shown that fractional velocity change can be recovered by time-shift tomography if care is taken to estimate the pre-stack time-shifts, and that uplift can be achieved against post-stack estimates (Izadian & MacBeth, 2021).Figure 14 shows one such result, in which fractional velocity changes calculated from both post-stack time-shifts and the pre-stack tomographic approach are compared.The result indicates that artefacts associated with calculation in the post-stack domain are reduced by capturing and utilizing the pre-stack variations.Another potential solution is to evaluate and apply time-shift pre-stack.This requires careful conditioning in the pre-stack domain, but if implemented effectively can provide a noticeable uplift to the final stacked amplitudes (Hatab & MacBeth, 2023).Finally, it should be mentioned that pre-stack time-shift analysis and its advantages are active and promising research areas that will be discussed in more detail in future publications.

Non-vertical ray-paths
This point is inextricably linked with the previous section on pre-stack variations.Literature on this topic focuses on the impact of non-flat reflectors assuming that post-stack data are interpreted as zero-offset normal rays to the reflector.Most post-stack estimation algorithms to date assume as a conceptual model that the time-shifts are measured along vertical ray-paths from the surface to each subsurface location (Mac-Beth et al., 2020).The relative velocity changes may therefore be recovered by vertical differentiation or a layer stripping inversion of the estimated time-shift volume.This practice is strictly valid only for sub-horizontal dips and weak velocity gradients in the subsurface and can lead to large errors in the positioning (in two-way time and horizontal distance) of the estimated velocity changes/time-strains in the case of larger dips.This is because zero-offset ray-paths normal to reflectors may not be vertical in the presence of a structurally complex overburden or dipping reservoir units, and for strong ray bending effects.Although pre-stack 3D processing corrects for reflector position in the final images, time-shifts created by a true subsurface sampling of the time-lapse effects along ray-paths normal to the reflector are not correctively aligned in this process.Thus, inversion to relative velocity change by applying algorithms that operate only vertically is erroneous.This point was made by Thore et al. (2010Thore et al. ( , 2012)), who observed that the 'wake' created by a time-lapse velocity perturbation in an overburden extension scenario for a dipping compacting reservoir did not occur along the vertical direction, but rather in a direction normal to the reservoir formation.To solve for the true velocity perturbations, they introduce a transformation that maps the image data to a pseudo-1D space where 'conventional' time-shift measurement techniques (as defined by the algorithms described in MacBeth et al., 2020) are valid.The transformation of the data consists of trace extractions along ray-paths starting normally to the reflector, which moves forwards and backwards to the surface.This approach does not solve all imaging problems, especially when sequences of bedding dips are non-parallel, and it is limited to a short range of offsets.Audebert and Agut (2014) later implemented the idea of Thore et al. (2010) in a slightly different way by working in the depth image domain and considering many local reflectors with small spatial extensions in parallel rather than a large discrete reflector as defined in the original method.The revised method reduces the computational cost and creates results that are localized to the reservoirs of interest.Following this, Karimi et al. (2016) introduced an alternate method based on time-lapse image registration in a stratigraphic coordinate system.They use a local similarity attribute to estimate time-shifts after flattening the time-lapse seismic images using the stratigraphic coordinate transformation of Karimi and Fomel (2015).The development above is followed by Khalil and Hoeber (2016) who consider a wave-equation-based image warping.They attempt to drive the warping/time-shift measurement more directly by the physics of seismic imaging to better handle dipping reflectors.Baseline and monitor images in depth are visualized as snapshots of a wavefield at different time steps.A series of seismic image waves respecting local structures can be propagated from this understanding.Spatially varying time-shifts are then computed by applying conventional 1D algorithms along this newly formed 'orthogonal' time-axis.1D time-shift measurement methods can now be employed.The advantage of the approach is that time-shifts are calculated in the direction normal to the wavefield, as recommended by Thore et al. (2012).This image-consistent technique is applied by Hoeber et al. (2016) to depth-migrated data from the Baobab field with a moderate degree of uplift over the conventional 1D approach applied vertically and post-stack.
As noted at the beginning of this section, the pre-stack timeshift tomography recommended in the previous section is a way of capturing the desired variations, including offset variations and dipping structure, and should be favoured against the recipe described above.This is because the corrections discussed in the above section are approximations, whereas the 4D tomography fully captures the offset-dependent behaviour.

Lateral shifts
Lateral shifts are less commonly reported than vertical timeshifts but have been considered and discussed in the past.These are observed in practice across several 4D datasets (Cox & Hatchell, 2008;Druzhinin & MacBeth, 2001;Fehmers et al., 2007;Hall et al., 2002Hall et al., ,2006;;Ji, 2017;Rickett & Lumley, 2001) and are typically quoted as up to 10 m in magnitude but have been measured as ±2 bins (Hall et al., 2006).Such shifts are not unsurprising, as well-known imaging considerations tell us that in the presence of a velocity change (particularly in the overburden), lateral shifts should always accompany vertical time-shifts in our data (e.g.Grubb et al., 2001).Indeed, the application of lateral and vertical shifts to the monitor data ('3D warping'), as in Hall et al. (2005), is in essence an inexpensive residual migration.When producing reservoirs experience significant geomechanical effects, additional shifts from subsurface displacement are superimposed on these imaging shifts due to velocity change.It is likely that the origin of velocity changes in the overburden is geomechanical in nature, so these shifts must always be superimposed.The image-based and geomechanics-based shifts appear to act in opposite directions (Cox & Hatchell, 2008).Methods for measurement of the full 3D (vertical, inline and crossline) warping have been discussed in MacBeth et al. (2020).For the application of such algorithms to data with horizontally lying geology, it is possible to calculate the vertical time-shifts independently from the horizontal shifts.These vertical shifts will be accurate, but their horizontal locations are inaccurate due to the lateral shifts.However, for the dipping structures, there is also a systematic error introduced into the measured vertical shift arising from the horizontal shift component which cannot be ignored (Leggott et al., 2000).In the presence of a dip, a degree of coupling between vertical and horizontal shifts is expected, and therefore, the magnitude of the vertical shift error is directly proportional to the dip and the magnitude of the horizontal shift component.Fehmers et al. (2007) showed that for small time-shifts, dipping structures can further obfuscate the true vertical timeshifts in the presence of navigational error.Cox and Hatchell (2008) modelled this error for geomechanical effects on the Shearwater field, UKCS.Errors at shallow and deeper target reflectors can range from a fraction of a millisecond to 1 ms.As more complex reservoirs with steeper dips and more subtle time-lapse signals of a few milliseconds are accessed with 4D seismic technology, we would expect 3D warping (or an equivalent image-based solution) to become a necessary consideration.To achieve reliable warping results, we will require algorithms that are sensitive to amplitude (or attribute) variation across the events of interest and perhaps some physical constraints on our solutions.

Amplitude effects
In practice, time-shifts estimated from the algorithms described in MacBeth et al. ( 2020) are most commonly applied as corrective dynamic alignments to the monitor volume, from which it is hoped to output reliable amplitude differences.The time-shifts and corrected amplitude differences can then be separately visualized as reservoir attribute volumes.The amplitude difference volumes may be further inverted into impedance changes using 4D inversion methods (early examples of this are Buland and El Quair (2006) and Sarkar et al. (2003), or more recently Zhan et al. (2017)).The basic premise behind this inversion is that the first stage of the workflow measures time-shifts, assuming amplitudes can be neglected, whereas the second stage inverts amplitudes to impedance changes assuming time-shifts can be neglected.
Although this workflow is perhaps the most common to date, Williamson et al. (2007) recognized that both time-shifts and amplitudes are physically coupled and proposed a joint inversion to 4D velocity changes.This coupling implies that some degrees of amplitude leakage must appear on the time-shifts, and vice-versa.In our study on tuning interferences (see section above), it was noted that most techniques performed well when applied iteratively as several cycles of estimate and correction, and that the effect due to amplitudes was minimal.Our previous comparison of methods with and without this time-shift-amplitude coupling suggested marginal improvement in the final results when amplitude was included as part of the time-shift solution (MacBeth et al., 2020).To check this further, we calculate the 4D seismic response from our fine-scale model in Figure 4 with and without changes in the reflectivity coefficients, that is the modelling includes only kinematic changes in the position of the reflection spikes.Time-shifts are then estimated from both synthetics.The results of time-shift estimation using NLI2D are shown in Figure 15a,c, whereas the differences from the true (model) values (shown in Figure 5a) are displayed in Figure 15b,d.The corresponding time-strains are given in Figure 15e,f.We note that errors in tracking the time-shifts inside the reservoir zone are amplified when the reflectivity changes are included, although estimates with and without these changes do look visually very similar, especially when a smooth gradational colour bar is employed.The differences become much more noticeable when comparing time-strains inside the reservoir and, to a lesser extent, in the overburden.This is because inside the reservoir interval itself, there are strong changes in reflection amplitude (sharp increase or decrease) and also the possibility of phase changes as waveforms constructively or destructively interfere due to 4D tuning between the baseline and the monitor.It has been shown above those interpretation difficulties (of a spurious time-shift) also arise around the tuning thickness and when there are wavelet interferences.An essential ingredient of approaches to negate this situation is a known seismic wavelet.The above study suggests that a clean separation of amplitude from time-shift effects becomes harder to identify near and inside the reservoir interval.4D inversion procedures that perform a cascaded inversion to 4D velocity changes using time-shifts to establish an initial model and then amplitudes in an iterative joint workflow (e.g.Li et al., 2017;Micksch et al., 2018;Zhan et al., 2017) or in a directly coupled scheme (Baek & Keho, 2015;Williamson et al., 2007) are the natural progression for this topic.Coupled approaches work well provided both kinematic and dynamic information has been treated equitably throughout the processing workflow, and one can rely on the wavelet.However, within the processing workflow, a different emphasis is placed on kinematics and amplitudes, with the kinematics tending to be more reliable.This explains why the time-shifts are often very good indicators of the underlying changes.To pro-  ceed further to eliminate the artefacts we observe, we need true amplitude processing, better velocity models, better event picking and time-shift corrections and better treatment of multiples.

DISCUSSION AND CONCLUSIONS
Errors are present in time-shifts (and hence also time-strains) estimated from post-stack seismic data.These arise from the intrinsic limitations of our seismic measurements, inaccuracies in our knowledge of the seismic acquisition or our inability to adequately capture the true physics in our methods.Table 2 summarizes our current understanding of the sources of this error that have been described in this study and their relative importance.It is understood that errors of many milliseconds are due to our seismic acquisition (which include navigation errors, timing errors, seawater velocity changes and tidal corrections).These can be measured and corrected by the careful and highly developed algorithms that currently exist within the industry.Although mistakes may inevitably be made in this process, they can be recognized and corrected due to judicious monitoring of 4D metrics during the pre-stack processing workflow.However, near-surface variability in land 4D seismic remains a problem despite significant advances in the past 5 years.The best course of action is to bury both receivers and sources, but this escalates project cost and may not be necessary if the 4D reservoir effects are large and time-shift errors can be tolerated.By comparison, intrinsic errors are unavoidable due to the nature of the seismic method itself and the subsurface geological heterogeneity.Small-scale fluctuations may arise as a result, which can be smoothed by regulariza-tion from time-shift methods.High resolution of time-lapse changes in thin beds may be obscured and needs specialized measurement.Fortunately, the impact of sampling error and non-repeatability noise for dedicated 4D seismic surveys is likely to be small (fractions of a millisecond) and low enough not to interfere with reservoir interpretation (normalized root mean square [NRMS] less than 20%).This may not be the case for land data in regions where the lack of repeatability may still be a major issue (NRMS of up to 60%).Thin bed interferences appear not to be problematic as they scale with the magnitude of the time-shifts and leave characteristic and identifiable artefacts that track along major events.These may be smoothed over by techniques that deploy spatial regularization, at the cost of resolution.Finally, most time-shift algorithms are based on approximations to the wave physics, which may introduce a modest degree of error (up to 1 ms) due to non-vertical normal rays at dipping reflectors, prestack variations due to overburden changes, lateral shifts and amplitude effects.Pre-stack variations due to ray-path sampling create the often-localized character of the time-shift measurements we typically observe.These items are known and can be avoided or incorporated with some degrees of cost to the user by extending the current post-stack estimation methods.
It is observed that despite the catalogue of errors detailed in this study, post-stack time-shifts still remain a reliable and cost-efficient indicator of subsurface changes due to production.These time-shifts remain robust often despite very complex overburden geology and are comprehensively interpreted across a range of fields (MacBeth et al., 2019).However, as 4D studies push the limits to deeper, more complex fields, with thinner, stacked reservoir sequences or salt body geometries with small, subtle time-shift signatures (van Gestel et al., 2019), it is probable that the inherent error sources we describe in this paper will become increasingly more important.These challenges will force us to improve the wave physics input into our techniques and shift our focus towards pre-stack analysis techniques such as 4D imagebased time-shift tomography (Izadian & MacBeth, 2021), and 4D full waveform inversion (Hicks et al., 2016).The benefit of such methods is beginning to reveal promise, but time will tell whether we encounter subsurface situations in which post-stack time-shifts are challenged as a solid source of subsurface information.Finally, we note that although the focus in this study has been on errors per se, there are still many time-shift signals that cannot be fully reconciled relative to the known history of reservoir production.This may be partly due to our lack of knowledge of the dynamic reservoir description, or uncertainty in interpreting these signals using existing fluid flow and geomechanical simulation models.These un-interpretable time-shifts require cataloguing and will be detailed in subsequent papers.

A C K N O W L E D G E M E N T S
We thank the sponsors of the Edinburgh Time-Lapse Project (ETLP), Phase VII and VIII (ADNOC, Aker BP, BHP, BP, CGG, Chevron, ConocoPhillips, Vår Energi AS, Exxon-Mobil, Halliburton, Harbour Energy, CNOOC, Neptune Energy, NORSAR, Petoro, Petrobras, Sharp Reflections, Shell, Equinor, TAQA, Tullow, Woodside) for supporting this research.We thank the sponsors for their support and Harbour Energy, ConocoPhillips and Equinor for the provision of data for the examples shown.We also thank Henning Hoeber and Celine Lacombe for helpful discussions on various aspects of this article, and an anonymous reviewer for offering improvements to Equations (1) and (2).Thanks also to Boshara Arshin, Mohamed Hatab and Sean Tian from ETLP for providing the data examples in Figures 6 and 7.The data examples are created with Schlumberger's Petrel software and ETLP's in-house software.

D A T A AVA I L A B I L I T Y S T A T E M E N T
The data that support the findings of this study are available from Harbour Energy, ConocoPhillips and Equinor.Confidentiality restrictions apply to the availability of these data, which were used under a contractual agreement for this study.Data are available from these companies by a written request for permission from the relevant company.

R E F E R E N C E S
not encountered in practice.Both uncorrelated and correlated components are present in the time-lapse data, leading to some noise cancellations when differencing (Hubans, 2016;Leguijt & Jenkins, 2013) and a 4D noise lower than √ 2N.In other words, the 4D effects appear better than that predicted from the 3D data alone.The invalidity of the noise assumption will result in a reduction of the 4D noise relative to 3D noise and could in turn imply that our conclusions are an overestimate of the true time-shift floor.
Sidenote: Although the discussion in the main text is focused on the error in the measured time-shift per se, attention must also be given to the noise levels in the 4D difference data generated directly by the time-shifts we are trying to measure.It is already known that time-shifts between the baseline and monitor generate significant 4D noise and errors on the difference amplitudes (Ross & Altan, 1997).We know that this creates geological leakage and magnification of the amplitudes (Dyce et al., 2004).The contribution to the NRMS is a well-documented function of bandwidth, SNR and energy ratio of the data (Burren & Lecerf, 2015;Lecerf et al., 2015).In published literature, observed NRMS values may be reported prior to time-shift correction, thus give a false impression of the true noise conditions.This will also be true of our time-shift algorithms until sufficient iterations have taken place to reach the noise floor.

T
A B L E 1 A selection of the smallest measurable time-shifts available from off-and onshore data, together with their corresponding normalized root mean square (NRMS) values.Also shown is the calculated Cramér-Rao (C-R) lower bound prediction for these data based on a seismic bandwidth of 30 Hz, central frequency of 20 Hz and an estimation window of 200 ms.An asterisk (*) indicates that the NRMS has been guessed in the absence of published values.

F
Predicted lower limit on the time-shift error calculated by the Cramér-Rao equation versus the normalized root mean square (NRMS) metric (from which the 3D signal-to-noise ratio (SNR) has been calculated).Lines correspond to: (1) -central frequency 20 Hz, window length 100 ms; (2) -central frequency 20 Hz, window length 200 ms; (3) -central frequency 35 Hz, window length 100 ms; (4) -central frequency 35 Hz, window length 200 ms.Solid circles correspond to the smallest values reported from field data shown in Table

F
I G U R E 2 (a) Time-shift estimates and (b) corresponding ΔV/V estimates for (c) a stack of layers of varying thickness and velocity with time-lapse perturbations applied.For (a)-(c), true time-shifts and ΔV/V are in red compared against the true model values in blue; (d) the baseline and monitor traces from which the time-shifts are estimated using the non-linear inversion (NLI) method.

F
I G U R E 3 Time-shift estimates for a North Sea dataset: (a) estimates using fast cross-correlation (FCC) method; (b) estimates using the non-linear inversion (NLI) method with spatial regularization.Details of these methods can be found in MacBeth et al. (2020).

F
I G U R E 4 (a) Inverted impedance from a North Sea dataset (baseline) together with the corresponding modelled seismic (c).The monitor impedance (b) has been adjusted to include the relative impedance changes originally derived 4D time-shift tomography of these data (d).

F
I G U R E 5 (a) Time-shift calculated from the model in Figure 4 directly by numerical integration; (b) time-shift estimated from the seismic data using NLI2D; (c) the difference between (a) and (b); (d) relative time-shift Δt/t derived from the estimates of (b), which should be compared to the model input in Figure 4d.

F
I G U R E 6 Two examples of spurious time-shifts: (a) an extreme example from the Sleipner field due to CO 2 injection.Time-shifts are estimated between towed streamer surveys shot in 2001 and 1994 with an normalized root mean square (NRMS) of 65%; (b) time-shifts estimated in a recent North Sea towed streamer 4D dataset (NRMS is 11%, and signal is after 4 years of production).Source: (a) Arshin, personal communication..cause unreliable and unusual time-shift signatures, and how this is treated pre-stack may influence the post-stack result.Finally, depth or lateral variations in the seismic wavelet must be minimized, as small phase rotations can be misconstrued by the algorithm as time-shifts.The above effects are magnified in the presence of complex, dipping overburden and reservoir structure, and deep water.Scattering and attenuative near-surface layers are also factors for land data.

F
I G U R E 7 Two time-shift estimates (a); and their corresponding time-strains (b); from the same North Sea data with different, but equally plausible, pre-stack processing workflows generating post-stack data realizations.We select realizations 1 and 2 for the purposes of this example.SeeHatab and MacBeth (2021) for details of the workflow perturbations.Top and base of main production zone are shown as black horizons.FI G U R E 8 Time-shifts estimated from a thin UKCS reservoir, in which the imprint of trim statics has remained.

F
Examples of typical survey time stamps for two North Sea acquisitions: (a) permanent reservoir monitoring (PRM) and (b) towed streamer.Both show the characteristic striping that inputs into the tidal (and to a lesser extent the water velocity-related) time-shift analysis, colour-coded by the time sequence of shooting (in days) of each shot line.Source: Adapted after Omofoma (2017).

F
Results of a modelling exercise similar to that of Figure4, in which an 80 m water layer has been added to the model: (a) water velocity profiles used, corresponding to winter (w) and summer (s) salinity and temperature conditions in the North Sea(Bertrand, 2005); (b) time-shifts measured for an unchanging water layer (in this case the winter-to-winter velocity profile); (c) time-shifts estimated when the seismic is modelled with summer (baseline) and then winter (monitor) profiles.

F
Time-shifts estimated from synthetic traces calculated for the fine-scale synthetic model of Figure 4b, compared against the true model values.The synthetics have been calculated with a residual component of free-surface multiple: (a) 10% residual amplitude; (b) 20% residual multiple; and (c) 30% residual multiple.
. Despite the accepted poor repeatability for surface-based acquisitions, there are also examples where robust processing flows have ensured confidence in the measured time-shifts for predominantly shallow, but also deep, reservoirs with fluid-related or geomechanical reservoir signals (see the examples shown by MacBeth et al., 2019).It is unfortunate that, in general, most land data acquired with surface sources and receivers still do have a high NRMS (40%-60%), and non-repeatability errors dominate time-shifts.

F
I G U R E 1 2 (a) Schematic illustrating offset-dependent ray-path coverage of a 4D anomaly.Parts (a) and (b) are post-stack time-shifts calculated by ray-path integration and summation over an offset range of 100-2350 m, for 4D anomalies embedded at 2000 m depth.In this case, the baseline velocity model is a constant 2000 m/s.Reservoir velocity changes, a speedup of 100 m/s, are confined to (b) a small rectangular region of 35 m radius and (c) a thin rectangular region of 525 m width and 100 m thickness.

F
Post-stack time-shifts estimated from four different field datasets showing the impact of pre-stack variations as evidenced by the decay of the time-shift tail with depth: (a) steam chamber development in Long Lake, Canada; (b) hydrocarbon gas injection into a saline aquifer in the UKCS; (c) pressure increase due to seawater injection into a UKCS field (Amini & Marsden, personal communication); (d) pressure increase from seawater injection on the Norne Field.Source: Examples (a), (b) and (d) are from MacBeth et al. (2019).F I G U R E 1 4 Relative velocity changes estimated: (a) post-stack on a field dataset using non-linear inversion (NLI); (b) from pre-stack time-shift tomography.Source: After Izadian and MacBeth (2021).

F
I G U R E 1 5 (a) Time-shift estimate from synthetics calculated using the fine-scale model of Figure 4b.The effect of changing reflection coefficient is absent; (b) difference of this estimate from true (model) value; (c) time-shift estimate with amplitude effects now included in the modelling; (d) corresponding difference between results of (c) and true (model) values; (e) relative velocity changes determined from (a); (f) relative velocity changes determined from (c).T A B L E 2 Estimated upper values for the various categories of time-shift error present in 4D seismic data.