## Robust joint full-waveform inversion of time-lapse seismic data sets with total-variation regularization

Musa Maharramov (maharram@stanford.edu) Biondo Biondi (biondo@sep.stanford.edu)

## ABSTRACT

We present a technique for reconstructing subsurface velocity model changes from timelapse seismic survey data using full-waveform inversion (FWI). The technique is based on simultaneously inverting multiple survey vintages, with model difference regularization using the total variation (TV) seminorm. We compare the new TV-regularized time-lapse FWI with the L 2 -regularized joint inversion proposed in our earlier work, using synthetic data sets that exhibit survey repeatability issues. The results demonstrate clear advantages of the proposed TV-regularized joint inversion over alternatives methods for recovering production-induced model changes that are due to both fluid substitution and geomechanical effects.

## INTRODUCTION

Effective reservoir monitoring depends on successful tracking of production-induced fluid movement in the reservoir and overburden, using input from seismic imaging, geomechanics, geology and reservoir simulation (Biondi et al., 1996). To achieve this, most traditional methods rely on the conversion of picked time shifts and reflectivity differences between migrated images into reflector movement and impedance changes. Though effective in practical applications, this approach requires a significant amount of expert interpretation and relies on quality control in the conversion process. Wave-equation image-difference tomography has been proposed as a more automatic alternative method to recover velocity changes (Albertin et al., 2006); it allows localized target-oriented inversion of model perturbations (Maharramov and Albertin, 2007). An alternative approach is based on using the high-resolution power of the full-waveform inversion (Sirgue et al., 2010) to reconstruct production-induced changes from wide-offset seismic acquisitions, and is the subject of this paper.

Time-lapse full-waveform inversion (Watanabe et al., 2004; Denli and Huang, 2009; Routh et al., 2012) is a promising technique for time-lapse seismic imaging where production-induced subsurface model changes are within the resolution of FWI. However, as with alternative time-lapse techniques, time-lapse FWI is sensitive to repeatability issues (Asnaashari et al., 2012). Non-repeatable acquisition geometries (e.g., slightly shifted source and receiver positions), acquisition gaps (e.g., due to new obstacles), different source signatures and measurement noise-all contribute to differences in the data from different survey vintages. Differences in the input data sets due to repeatability issues may easily mask valuable production-induced changes. However, even with noise-free synthetic data without any acquisition repeatability issues, numerical artifacts may contaminate the inverted difference of monitor and baseline when practical limitations are imposed on solver

iteration count. Maharramov and Biondi (2013) devised a time-lapse FWI that minimized model differences outside of areas affected by production by jointly inverting for multiple models, and imposing a regularization condition on the model difference. The joint inversion can be performed simultaneously for multiple model vintages or using an empirical technique of 'cross-updating' (Maharramov and Biondi, 2013, 2014b). Maharramov and Biondi (2014b,a) applied these methods to noisy synthetic data and compared the results to alternative methods, demonstrating advantages of both target-oriented simultaneous inversion and cross-updating over alternative methods. Production-induced changes that are due to fluid substitution effects are spatially constrained to areas affected by fluid transport, while those due to geomechanical effects-e.g., stretching of the overburden-may result in smoothly varying velocity model differences. In either case, regularization terms that promote 'sparsity' of spatial model-difference gradients may be expected to improve the recovery of valuable production effects while suppressing undesirable oscillatory artifacts. In this work we present the results of using total-variation (TV) regularization of the model difference in the simultaneous time-lapse FWI and demonstrate the significant improvement of the inversion results compared to alternative methods.

## THE METHOD

Full-waveform inversion is defined as solving the following optimization problem (Tarantola, 1984; Virieux and Operto, 2009)

$$\| M u - \mathrm d \| _ { 2 } \, \to \, \min$$

where M d , are the measurement operator and data, u is the solution of a forward-modeling problem

$$D ( \mathbf m ) \mathbf u \, = \, \phi,$$

where D is the forward-modeling operator that depends on a model vector m as a parameter, and φ is a source. The minimization problem (1) is solved with respect to either both the model m and source φ or just the model. In the frequency-domain formulation of the acoustic waveform inversion, the forward-modeling equation (2) becomes

$$- \, \omega ^ { 2 } u - v ^ { 2 } ( x ^ { 1 }, \dots, x ^ { n } ) \Delta u \, = \, \phi ( \omega, x ^ { 1 }, \dots, x ^ { n } )$$

where ω is a temporal frequency, n is the problem dimension, and v is the acoustic wave propagation velocity. Values of the slowness s = 1 /v at all the points of the modeling domain constitute the model parameter vector m . The direct problem (3) can be solved in the frequency domain, or in the time domain followed by a discrete Fourier transform in time (Virieux and Operto, 2009). The inverse problem (1) is typically solved using a multiscale approach, from low to high frequencies, supplying the output of each frequency inversion to the next step.

FWI applications in time-lapse problems seek to recover induced changes in the subsurface model using multiple data sets from different acquisition vintages. For two surveys sufficiently separated in time, we call such data sets (and the associated models) baseline and monitor .

Time-lapse FWI can be carried out by separately inverting the baseline and monitor models ( parallel difference ) or inverting them sequentially with, e.g., the baseline supplied

as a starting model for the monitor inversion ( sequential difference ). Another alternative is to apply the double-difference method, with a baseline model inversion followed by a monitor inversion that solves the following optimization problem

$$\| \left ( M _ { m } ^ { s } \mathbf u _ { m } - M _ { b } ^ { s } \mathbf u _ { b } \right ) - \left ( M _ { m } \mathbf d _ { m } - M _ { b } \mathbf d _ { b } \right ) \| _ { 2 } \, \to \, \min \, & & ( 4 )$$

by changing the monitor model (Watanabe et al., 2004; Denli and Huang, 2009; Zheng et al., 2011; Asnaashari et al., 2012; Raknes et al., 2013). The subscripts in equation (4) denote the baseline and monitor surveys, d denotes the field data, and the M 's are measurement operators that project the synthetic and field data onto a common grid. The superscript s indicates the measurement operators applied to the synthetic data.

In all of these techniques, optimization is carried out with respect to one model at a time, albeit of different vintages at different stages of the inversion. In our method we invert for the baseline and monitor models simultaneously by solving either one of the following two optimization problems:

$$\alpha \| \mathsf M _ { b } \mathsf u _ { b } - \mathsf d _ { b } \| _ { 2 } ^ { 2 } + \beta \| \mathsf M _ { m } \mathsf u _ { m } - \mathsf d _ { m } \| _ { 2 } ^ { 2 } +$$

$$\gamma \| \left ( M _ { m } ^ { s } u _ { m } - M _ { b } ^ { s } u _ { b } \right ) - \left ( M _ { m } d _ { m } - M _ { b } d _ { b } \right ) \| _ { 2 } ^ { 2 } +$$

$$\alpha _ { 1 } \| W _ { b } R _ { b } ( m _ { b } - m _ { \stackrel { b } { m - \dots - } } ^ { \text{PRIOR} } ) \| _ { 2 } ^ { 2 } +$$

$$\beta _ { 1 } \| \mathbf W _ { m } \mathbf R _ { m } ( \mathbf m _ { \stackrel { \dots } { \dots } \dots } - \mathbf m _ { m } ^ { \text{PRIOR} } ) \| _ { 2 } ^ { 2 } +$$

$$\delta \| \text{WR} ( m _ { m } - m _ { b } - \Delta m ^ { \text{PRIOR} } ) \| _ { 2 } ^ { 2 } \, \to \, \min,$$

$$\alpha \| \mathbf M _ { b } \mathbf u _ { b } - \mathbf d _ { b } \| _ { 2 } ^ { 2 } + \beta \| \mathbf M _ { m } \mathbf u _ { m } - \mathbf d _ { m } \| _ { 2 } ^ { 2 } + \mathbf ( 1 0 )$$

$$\gamma \| \left ( M _ { m } ^ { s } u _ { m } - M _ { b } ^ { s } u _ { b } \right ) - \left ( M _ { m } d _ { m } - M _ { b } d _ { b } \right ) \| _ { 2 } ^ { 2 } + \quad & ( 1 1 )$$

$$\alpha _ { 1 } \| W _ { b } \mathbf R _ { b } ( \mathbf m _ { b } - \mathbf m _ { \stackrel { \text{PRIOR} } { \text{-\dots } } } ) \| _ { 1 } + \text{ (12)}$$

$$\beta _ { 1 } \| \mathbf W _ { m } \mathbf R _ { m } ( \mathbf m _ { \stackrel { \dots } { \dots } \dots } \, \mathbf m _ { m } ^ { \text{PRIOR} } ) \| _ { 1 } +$$

$$\delta \| \text{WR} ( m _ { m } - m _ { b } - \Delta m ^ { \text{PRIOR} } ) \| _ { 1 } \, \to \, \min,$$

with respect to both the baseline and monitor models m b and m m . Problem (5-9) describes time-lapse FWI with L 2 regularization of the individual models (7,8) and model difference (9) (Maharramov and Biondi, 2014b). In this work we study the second formulation (10-14) that involves L 1 -regularization of the individual models and their difference. The terms (10) correspond to separate baseline and monitor inversions, the term (11) is the optional double difference term, the terms (12) and (13) are optional separate baseline and monitor inversion regularization terms (Aster et al., 2012), and the term (14) represents regularization of the model difference. In (12)-(14), R and W denote regularization and weighting operators respectively, with the subscript denoting the survey vintage where applicable. If R is the gradient magnitude operator

$$\mathbf R f ( x, y, z ) \, = \, \sqrt { f _ { x } ^ { 2 } + f _ { y } ^ { 2 } + f _ { z } ^ { 2 } },$$

then (12-14) become total-variation (TV) seminorms . The latter case is of particular interest in this work as the minimization of the L 1 norm of gradient may promote 'blockiness' of the model-difference, potentially reducing oscillatory artifacts (Aster et al., 2012).

A joint inversion approach has been applied earlier to the linearized waveform inversion (Ayeni and Biondi, 2012). In Maharramov and Biondi (2013, 2014b,a), a simultaneous or

full-waveform inversion problem (5,9) was studied with a single model difference L 2 regularization term (9).

An implementation of the proposed simultaneous inversion algorithm requires solving a nonlinear optimization problem with twice the data and model dimensions of problems (1) and (4). The model difference regularization weights W and, optionally, the prior ∆ m PRIOR may be obtained from prior geomechanical information. For example, a rough estimate of production-induced velocity changes can be obtained from time shifts (Hatchell and Bourne, 2005; Barkved and Kristiansen, 2005) and used to map subsurface regions of expected production-induced perturbation, and optionally provide a difference prior. However, successfully solving the L 1 -regularized problem (10-11) is less sensitive to choice of the weighting operator W . For example, we show below that the TV-regularization using (15) with W = 1 recovers non-oscillatory components of the model difference, while the L 2 approach would result in either smoothing or uniform reduction of the model difference.

Figure 1: Starting model used in the inversion.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0645v2-figure-8679901b-c35d-4632-be1e-6bedd7216569.jpg)

In addition to the fully simultaneous inversion, Maharramov and Biondi (2013, 2014b) proposed and tested a cross-updating technique that offers a simple but remarkably effective approximation to minimizing the objective function (5),(9), while obviating the difference regularization and weighting operators R and W for problem (5,9). This technique consists of one standard run of the sequential difference algorithm, followed by a second run with the inverted monitor model supplied as the starting model for the second baseline inversion

$$\begin{matrix} m _ { \text{INIT} } \rightarrow \text{baseline inversion} \rightarrow \text{monitor inversion} \rightarrow \\ & \text{baseline inversion} \rightarrow \text{monitor inversion}, \end{matrix}$$

and computing the difference of the latest inverted monitor and baseline models. Process (16) can be considered as an approximation to minimizing (5) and (9) because nonrepeatable footprints of both inversions are propagated to both models, canceling out in the difference. Both the simultaneous inversion and cross-updating minimize the model difference by tackling model artifacts that are in the null space of the Fr´chet derivative of e the forward modeling operator. The joint inversion minimizes the effect of such artifacts on the model difference by either minimizing the model difference term (9) in the simultaneous inversion, or by propagating these artifacts to both models in cross-updating (16). Note that this process is not guaranteed to improve the results of the baseline and monitor model inversions but was only proposed for improving the model difference. Maharramov

and Biondi (2014b,a) demonstrated a significant improvement of model difference recovery by both the L 2 -regularized target-oriented simultaneous inversion and cross-updating compared to the parallel, sequential and double difference techniques. The simultaneous inversion and cross-updating yielded qualitatively similar results within the inversion target.

Here we compare joint simultaneous inversion with a TV-regularized model difference (10,14,15) to parallel difference and cross-updating.

Figure 2: Baseline model inverted from noise-free synthetic data.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0645v2-figure-dc4d1ea0-7b2c-47b0-afd6-5602c05f4a33.jpg)

## NUMERICAL EXAMPLES

The Marmousi velocity model is used as a baseline, over a 384 × 122 grid with a 24 m grid spacing. Production-induced velocity changes are modeled as a negative -150 m/s perturbation at about 4.5 km inline 800 m depth, and a positive 200 m/s perturbation at 6.5 km inline, 1 km depth. Additionally, a smoothly varying negative velocity change, peaking at -50 m/s, was included above the positive anomaly as shown in Figure 4. The whole Marmousi model is inverted, however, only model differences for the section between the approximate inline coordinates 4 km and 6.7 km to the depth of approximately 1.4 km are shown here. The inversion is carried out in the frequency domain for 3.0, 3.6, 4.3, 5.1, 6.2, 7.5, 9.0, 10.8, 12.8, and 15.5 Hz, where the frequencies are chosen based on the estimated offset to depth range of the data (Sirgue and Pratt, 2004). The baseline acquisition has 192 shots at a depth of 16 m with a 48 m spacing, and 381 receivers at a depth of 15 m with a 24 m spacing. The minimum offset is 48 m. The source function is a Ricker wavelet centered at 10.1 Hz. Absorbing boundary conditions are applied along the entire model boundary, including the surface (thus suppressing multiples). A smoothed true model (Figure 1) is used as a starting model for the initial baseline inversion (and for the initial monitor inversion in the parallel difference). The smoothing is performed using a triangular filter with a 20-sample half-window in both vertical and horizontal directions. The result of inverting the baseline model from the clean synthetic data is shown in Figure 2. Random Gaussian noise is added to the noise-free synthetic data to produce a noisy data set with 7 dB signal-to-noise ratio. The noisy monitor data set is generated for the model perturbation of Figure 4, using the same acquisition geometry and source wavelet. The results of baseline model inversion from the clean and 7 dB SNR synthetic data are shown in Figure 3a and Figure 3b, respectively. Results of model difference inversion from the clean

and 7 dB SNR synthetic data sets using various methods are shown in Figures 5a,5b,5c and Figures 6a,6b,6c, respectively. The simultaneous inversion objective function contains only terms (10) and (14) with no difference prior, i.e., ∆ m PRIOR = 0. The model-difference regularization weights W in (14) are set to 1 everywhere in the modeling domain. The two terms in (10) are of the same magnitude and therefore α and β are set to 1. Parameter δ is set to 10 -5 but can be varied for different acquisition source and geometry parameters. The result of the initial baseline inversion is supplied as a starting model for both m b and m m in the simultaneous inversion. In all the inversions, up to 10 iterations of the nonlinear conjugate gradients algorithm (Nocedal and Wright, 2006) are performed for each frequency. Neither regularization nor model priors are used in single-model inversions (i.e., in the crossupdating and parallel difference methods). Maharramov and Biondi (2014b) demonstrated significant improvement by cross-updating compared to sequential differencing, and rough qualitative equivalence of cross-updating and the L 2 -regularized simultaneous inversion. Therefore, in this work we compare the TV-regularized simultaneous inversion only against parallel differencing and cross-updating.

Figure 3: (a) Target area of the baseline model inverted from clean synthetic. (b) Target area of the baseline model inverted from a 7 dB SNR synthetic. In both cases the baseline model is reconstructed reasonably well, however, errors due to noise are comparable in magnitude to production-induced effects.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0645v2-figure-989d5e26-7ea6-4412-8a61-19ca837ec236.jpg)

The results of applying cross-updating to the two data sets are shown in Figures 5b and 6b, respectively. The corresponding TV-regularized simultaneous inversion results are shown in Figures 5c and 6c. Since problem (1) is nonlinear, supplying the result of the highest frequency inversion back to the lowest frequency and repeating the whole inversion cycle for all frequencies may result in achieving a better data fit. In repeated cycles, lowerfrequency inversions usually terminate earlier but higher frequencies still deliver model updates. For an objective comparison of the joint inversion with the parallel difference method, the effects of insufficient iteration count are reduced by performing an extra cycle of baseline and monitor inversion (we call this approach 'iterated' parallel difference (Maharramov and Biondi, 2014b)). The results of applying the iterated parallel difference to the twodata sets are shown in Figures 5a and 6a. While cross-updating demonstrates certain robustness with regard to uncorrelated noise in the data and computational artifacts (note the significant quantitative improvement of reconstructed difference magnitudes in Figures 5b and 6b), the TV-regularized achieves further signifcant improvement by reducing oscilatory artifacts and honoring both smooth and blocky components of the model difference.

Figure 4: True velocity difference consists of a negative ( -150 m/s) perturbation at about 4.5 km inline 800 m depth, and a positive (200 m/s) perturbation at 6.5 km inline, 1 km depth.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0645v2-figure-24dad3d0-27ce-4c53-b11b-6e9c7c35855d.jpg)

## CONCLUSIONS

Our new TV-regularized simultaneous inversion technique is a more robust further development of our previous joint inversion method (Maharramov and Biondi, 2013, 2014b,a). Use of TV regularization in the simultaneous inversion allows recovery of production-induced changes without specifying variable weighting operator W , and penalizes unwanted model oscillations that may mask useful production-induced changes.

One potentially beneficial extension of this work is using TV regularization of individual models in (7,9). This approach may used for both time-lapse and standard FWI, and will be the subject of our next work.

## ACKNOWLEDGMENTS

The authors would like to thank Stewart Levin for a number of useful discussions, and the Stanford Center for Computational Earth and Environmental Science for providing computing resources.

## REFERENCES

Albertin, U., P. Sava, J. Etgen, and M. Maharramov, 2006, Adjoint wave-equation velocity analysis: 76th Annual International Meeting, SEG, Expanded Abstracts, 3345-3349. Asnaashari, A., R. Brossier, S. Garambois, F. Audebert, P. Thore, and J. Virieux, 2012, Time-lapse imaging using regularized FWI: A robustness study: 82nd Annual International Meeting, SEG, Expanded Abstracts, doi:10.1190/segam2012-0699.1, 1-5. Aster, R., B. Borders, and C. Thurber, 2012, Parameter estimation and inverse problems: Elsevier.

Ayeni, G., and B. Biondi, 2012, Time-lapse seismic imaging by linearized joint inversion-A Valhall Field case study: 82nd Annual International Meeting, SEG, Expanded Abstracts, doi:10.1190/segam2012-0903.1, 1-6.

Barkved, O., and T. Kristiansen, 2005, Seismic time-lapse effects and stress changes: Examples from a compacting reservoir: The Leading Edge, 24 , 1244-1248.

456

1176

1416

456

1176

1416

4056

4056

5256

5256

6456

6456

b

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0645v2-figure-67dc00ed-7e67-4f4b-98ba-8659a8b3f9c7.jpg)

1.35e+02

2.70e+02

1.01e+02

1.510+02

2.010+02

~9.14e+01

9.14e+01

1.22e+02

c

Figure 5: Model difference inverted from a clean synthetic with matching baseline and monitor acquisition geometries using (a) iterated parallel difference; (b) cross-updating; (c) TV-regularized simultaneous inversion.

a

456

1176

1416

456

1176

1416

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0645v2-figure-24cc12e1-15d6-4859-8fb8-d89c2dcf6e8e.jpg)

a

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0645v2-figure-918a4e46-2057-43f4-a8bf-d8198b197195.jpg)

b

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0645v2-figure-b5c0cab2-ef41-4899-96cd-72dfed7ea9b3.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0645v2-figure-c711a50b-850d-4ff0-b9c1-90f61d738b2d.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0645v2-figure-bc5bfafe-fe6b-45cb-9ea1-731b80bab53e.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0645v2-figure-619f0ee4-7b6e-4b67-91ee-78459c6e0e40.jpg)

c

Figure 6: Model difference inverted from a 7 dB SNR synthetic with matching baseline and monitor acquisition geometries using (a) iterated parallel difference; (b) cross-updating; (c) regularized simultaneous inversion.

Biondi, B., C. Deutsch, R. Gundes, D. Lumley, G. Mavko, T. Mukerji, J. Rickett, and M. Thiele, 1996, Reservoir monitoring: A multidisciplinary feasibility study: 66th Annual International Meeting, SEG, Expanded Abstracts, 1775-1778.

Denli, H., and L. Huang, 2009, Double-difference elastic waveform tomography in the time domain: 79th Annual International Meeting, SEG, Expanded Abstracts, 2302-2306.

Hatchell, P., and S. Bourne, 2005, Measuring reservoir compaction using time-lapse timeshifts: 75th Annual International Meeting, SEG, Expanded Abstracts, 2500-2503.

Maharramov, M., and U. Albertin, 2007, Localized image-difference wave equation tomography: 77th Annual International Meeting, SEG, Expanded Abstracts, 3009-3013.

Maharramov, M., and B. Biondi, 2013, Simultaneous time-lapse full waveform inversion: SEP Report, 150 , 63-70.

---, 2014a, Joint 4dfwi with model-difference regularization. seg-agu summer research workshop: SEG-AGU Summer Research Workshop. Advances in Active+Passive 'Full Wavefield' Seismic Imaging: From Reservoirs to Plate Tectonics.

- ---, 2014b, Joint full-waveform inversion of time-lapse seismic data sets: SEP Report, 152 , 19-28.

Nocedal, J., and S. J. Wright, 2006, Numerical optimization: Springer.

Raknes, E., W. Weibull, and B. Arntsen, 2013, Time-lapse full waveform inversion: Synthetic and real data examples: 83rd Annual International Meeting, SEG, Expanded Abstracts, 944-948.

Routh, P., G. Palacharla, I. Chikichev, and S. Lazaratos, 2012, Full wavefield inversion of time-lapse data for improved imaging and reservoir characterization: 82nd Annual International Meeting, SEG, Expanded Abstracts, doi:10.1190/segam2012-1043.1, 1-6.

Sirgue, L., O. Barkved, J. Dellinger, J. Etgen, U. Albertin, and J. Kommedndal, 2010, Full waveform inversion: the next leap forward in imaging at Valhall: First Break, 28 , 65-70. Sirgue, L., and R. Pratt, 2004, Efficient waveform inversion and imaging: A strategy for selecting temporal frequencies: Geophysics, 69 , 231-248.

Tarantola, A., 1984, Inversion of seismic reflection data in the acoustic approximation: Geophysics, 49 , 1259-1266.

Virieux, J., and S. Operto, 2009, An overview of full-waveform inversion in exploration geophysics: Geophysics, 74 , no. 6, WCC1-WCC26.

Watanabe, T., S. Shimizu, E. Asakawa, and T. Matsuoka, 2004, Differential waveform tomography for time-lapse crosswell seismic data with application to gas hydrate production monitoring: 74th Annual International Meeting, SEG, Expanded Abstracts, 2323-2326.

Zheng, Y., P. Barton, and S. Singh, 2011, Strategies for elastic full waveform inversion of timelapse ocean bottom cable (OBC) seismic data: 81st Annual International Meeting, SEG, Expanded Abstracts, 4195-4200.
Musa Maharramov Biondo Biondi
Robust joint full-waveform inversion of time-lapse seismic data sets with total-variation regularization
We present a technique for reconstructing subsurface velocity model changes
from time-lapse seismic survey data using full-waveform inversion (FWI). The
technique is based on simultaneously inverting multiple survey vintages, with
model difference regularization using the total variation (TV) seminorm. We
compare the new TV-regularized time-lapse FWI with the $L_2$-regularized joint
inversion proposed in our earlier work, using synthetic data sets that exhibit
survey repeatability issues. The results demonstrate clear advantages of the
proposed TV-regularized joint inversion over alternatives methods for
recovering production-induced model changes that are due to both fluid
substitution and geomechanical effects.