# Four-Dimensional History Matching Using ES-MDA and Flow-Based Distance-to-Front Measurement

Eduardo Barrela, Philippe Berthet, Mario Trani, Olivier Thual, Corentin Lapeyre

► To cite this version:

Eduardo Barrela, Philippe Berthet, Mario Trani, Olivier Thual, Corentin Lapeyre. Four-Dimensional History Matching Using ES-MDA and Flow-Based Distance-to-Front Measurement. Energies, 2023, 16 (24), pp.7984. 10.3390/en16247984. hal-04691117

HAL Id: hal-04691117

https://hal.science/hal-04691117v1

Submitted on 7 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Four-Dimensional History Matching Using ES-MDA and Flow-Based Distance-to-Front Measurement

Eduardo Barrela <fnref n="1" />*, Philippe Berthet <fnref n="1" />*, Mario Trani <fnref n="1" />, Olivier Thual <fnref n="2" /> and Coretin Lapeyre <fnref n="2" />*

**Abstract:** The use of 4D seismic data in history matching has been a topic of great interest in the hydrocarbon industry as it can provide important information regarding changes in subsurfaces caused by fluid substitution and other factors where well data is not available. However, the high dimensionality and uncertainty associated with seismic data make its integration into the history-matching process a challenging task. Methods for adequate data reduction have been proposed in the past, but most address 4D information mismatch from a purely mathematical or image distance-based standpoint. In this study, we propose a quantitative and flow-based approach for integrating 4D seismic data into the history-matching process. By introducing a novel distance parametrization technique for measuring front mismatch information using streamlines, we address the problem from a flow-based standpoint; at the same time, we maintain the amount of necessary front data at a reduced and manageable amount. The proposed method is tested, and its results are compared on a synthetic case against another traditional method based on the Hausdorff distance. The effectiveness of the method is also demonstrated on a semi-synthetic model based on a real-case scenario, where the standard Hausdorff methodology could not be applied due to high data dimensionality.

**Keywords:** four-dimensional seismic; history matching; ensemble smoother with multiple data assimilation; distance-to-front; streamlines

## 1. Introduction

In this paper, we aim to explore and demonstrate the effectiveness of front re-parametrization methods in the context of history matching for hydrocarbon reservoirs. Specifically, we focus on a novel distance-to-front method using streamlines obtained from full physics flow simulation aimed at enhancing the accuracy and efficiency of history-matching processes in reservoir management. History matching is an important stage in the development of every hydrocarbon field, playing a crucial role in both the modeling and simulation phases. During the process of history matching, key properties of the reservoir model are calibrated to match past production data. This not only allows inference of reservoir properties from production data but also, ultimately, allows obtaining predictive model(s) that reduce uncertainty in forecasting future production. The end goal is to use predictive reservoir models to aid the development of future production strategies and support management decisions. The reliability of a reservoir model can be estimated when it is able to reproduce all available data as accurately as possible. Currently, many oil and gas fields have long historical data that can help constrain the history matching procedure, as well as seismic data acquired during their production life.

The incorporation of time-lapse seismic data (4D seismic data) into the history-matching workflow has been a topic of great interest over the years. Four-dimensional seismic data can provide information in areas of the reservoir where no data are available; more specifically, the data can be used as a tool to monitor changes in the subsurface

<figure><img src="image_2.png" /></figure>

check for
updates

**Citation:** Barrela, E.; Berthet, P.; Trani, M.; Thual, O.; Lapeyre, C. Four-Dimensional History Matching Using ES-MDA and Flow-Based Distance-to-Front Measurement. *Energies* **2023**, *16*, 7984. https://doi.org/10.3390/en16247984

Academic Editors: Haifeng Zhao and Yang Xia

Received: 2 October 2023

Revised: 20 November 2023

Accepted: 25 November 2023

Published: 9 December 2023

<figure><img src="image_3.png" /></figure>

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

<footnote n="1">¹ TotalEnergies S.E.—Centre Scientifique & Technique Jean Féger, Av. Larribau, 64000 Pau, France</footnote>

<footnote n="2">² Centre Européen de Recherche et de Formation Avancée en Calcul Scientifique, 42 Av. Gaspard Coriolis, 31100 Toulouse, France</footnote>

<footnote>* Correspondence: eduardo-jose.airoso-barrela@totalenergies.com (E.B.); philippe.berthet@totalenergies.com (P.B.)</footnote>
originating from fluid substitution due to production [1–5] or changes in other subsurface
properties such as temperature or pressure [6]. These insights are vital for enhancing the
accuracy of history matching by aiding in more precisely calibrating reservoir models
to reflect actual conditions. Thus, the integration of 4D seismic data plays a pivotal role
in overcoming the limitations of conventional data, offering a more comprehensive un-
derstanding of the reservoir and ensuring more reliable and effective history-matching
outcomes.

From an industry standpoint, 4D seismic data have been primarily used as a qualitative constraint on the reservoir model [7–12] as the understanding of the reservoir evolves and revisions are made in a multidisciplinary framework, encompassing the domains of geology, seismic, and reservoir engineering [13]. However, while it is possible to identify the changes occurring on the subsurface, like saturation and pressure variations (e.g., [14,15]), using 4D data in a qualitative way, we have no information on the magnitude of the change or with what certainty.

On top of that, seismic interpretation can be subject to variations in its interpretation, even when facing easy-to-interpret seismic data. In a study by Rankey and Mitchell [16], the authors focus on the subjectivity of seismic data interpretation by presenting the same seismic data to different interpreters. The interesting conclusion was that although the data were unanimously considered easy to interpret, overconfidence in the interpreters led to variations in the interpretation, which in turn reflected considerable differences in volumetrics estimation. Apart from that, other uncertainties arising from non-repeatability effects [17], noise, and imaging are present in 4D seismic data. On the plus side, this can provide a great opportunity for the quantitative integration of 4D seismic data in history-matching workflows as it allows access to uncertainty quantification over model parameter estimates and fluid production forecasts analysis (e.g., [18–21]). Advancements over the past decades, leading to the proposal and development of a range of stochastic seismic inversion techniques, have provided ways of generating an ensemble of alternative heterogeneous impedance representations that agree with the 3D seismic volume, accounting for the non-uniqueness of the inversion process. Nonetheless, the full integration of quantitative 4D seismic data interpretation into the history-matching procedure is far from straightforward, remaining a challenge to be addressed as well as a topic of great interest. Several examples of quantitative approaches for integrating 4D seismic data exist in the literature (e.g., [10,22–31]). Problems identified with the adoption of such techniques on an industrial level are related to the practical feasibility and the inexistence of a fully integrated software solution that can easily handle the integration of both production and seismic data in a computer-aided history matching loop. However, the main concern is related to the computational feasibility of incorporating the large amount of data associated with seismic acquisition into existing workflows [32]. In particular, aspects related to the high nonlinearity of the problem at hand include the limitation on the number of degrees of freedom associated with the amount of data to be assimilated, questions on how to address the contribution in the assimilation procedure of the different types of data in consideration, how to elect relevant parameters for matching both seismic data and fluid flow production while staying within the boundaries of a plausible geological and physical model, over-conservative prior assumptions, and errors arising from the forward modelling of seismic data and attributes [33]. Therefore, due to the complexity and computational cost associated with the modeling of seismic attributes, presented workflows often fall into the categories of either being able to provide only unique solutions to the problem or requiring a significant reduction of the uncertainty space [34]. For this reason, alternative methods for adding seismic data information to the history-matching procedure have been explored. Recently, Rollmann et al. [35] presented a method using a convolutional neural network trained to fit observed seismic history. However, results were only shown in a synthetic case, and the overhead cost of gathering the necessary amount of training data (large amounts of models that need to be classified) as well the time spent in appropriate architecture development (which can be very case-specific) and computational costs associated with the
training of the network are still big disadvantages. Furthermore, the integration of such methods into the 4D history-matching workflow carries additional challenges (especially under real or realistic case scenarios), remaining a topic for future research.

In an attempt to address the problem related to the varying nature of data within the 4D seismic history matching procedure, Tillier et al. [36] proposed a method based on the local modified Hausdorff distance for measuring the dissimilarity between observed and simulated seismic attributes. This idea was later expanded on by Abadpour et al. [37], where coupled with an Ensemble Kalman Filter (EnKF) [38] workflow, the same Hausdorff distance metric was used to compute the distance between observed and saturation fronts, synthesized by binary image-based data reduction. The method showed promising results in a synthetic case; however, unless some other precursory data reduction steps are considered, the direct application of this method, particularly within the scope of ensemble-based workflows, becomes difficult to achieve in large-scale models as it implies computing the inverse of a gain matrix that could be at least the square of the number of cells in the model. There are also other issues related to binary image-based approaches, such as the double-penalty effect (i.e., when a feature is predicted where it should not be and is not predicted where it should), which, together with the amount and extension of data, add complexity to the minimization procedure even if the mismatch can be easily quantifiable.

Other binary image-based methods have also been proposed over the past decade (e.g., [39–43]). Similarly, such methods mainly focus on the conversion of hardening and softening signals of 4D seismic data to a binary image. However, while providing an effective and straightforward way of integrating 4D seismic information by reducing the level of information in a continuous 4D signal into discrete states (0 and 1), the main drawback is still related to the amount of data resulting from the computation of the difference between observed and measured responses.

Observed seismic amplitude fronts can be compared to saturation fronts when under cases where pressure variation or compaction effects [1,34] or the effect of variation in porosity or net-to-gross ratio can be neglected [33]. These can capture the main information related to the drained area of a reservoir under production. Like binary image-based methods, they can be looked at as a solution for applying data reduction to the problem at hand. This type of front re-parametrization reduces the amount of seismic data to be used by representing swept regions through a saturation front, which can capture the most significant 4D seismic information. A reduction in the nonlinearity of the problem is also achieved since front positions are closely related to uncertain petrophysical parameters of the model (before or after the front location). Finally, full seismic inversion procedures are avoided, and the method remains an option even when facing low-quality seismic data sets. All of these qualities place such methods as good candidates for application on history-matching workflows with the potential for increasing the performance of history-matching workflows.

Kretz [44] proposed a history-matching workflow to match front positions based on streamlined simulation. In their proposal, model permeabilities along the streamlines were modified in order to match observed and simulated saturation front positions derived from 4D seismic data. The discrepancy between front positions was provided by the difference in time-of-flight measured from the streamlines. The method showed great promise and served as a starting point for research on other front parametrization-based methods used in 4D history-matching workflows (e.g., [34,39,45,46]). However, the main drawback with this method relates to the properties being updated only along cells intersected by the path of the streamlines. This leads to models that lose their geologic consistency, no longer honoring geostatistical assumptions, and could also lead to overfitting of the matched production data. Finally, the application of the proposed method to realistic 3D cases was not discussed, and results were presented only in simple 2D synthetic examples.

Trani et al. [34] proposed the re-parameterization of saturations extracted from 4D seismic data in terms of front arrival times. The main disadvantage of this method was the need to run the fluid flow simulations beyond the update time at which the seismic data
are available (i.e., until all cells of the observed front location have been flooded). Later, Trani [47,48] concluded that an ensemble of complete predictions with extended simulation times could be replaced by an approximation of late arrival time by an arbitrarily large value. In an attempt to address the problem with extended simulation times and based on the work by Tillier et al. [36], Leeuwenburgh and Arts [45] and Zhang and Leeuwenburgh [46,47] proposed re-parameterizing front time-of-arrivals into distance-to-fronts. In the work presented by Leeuwenburgh and Arts [45], the authors assumed that uniform velocities in a monotonically expanding 3D front could be converted to distances calculated using a fast-marching algorithm [49] for the solution of Eikonal equations to cartesian grids [50]. These computed distances were then used as innovations in an ensemble history-matching framework using the EnKF. Later, Zhang and Leeuwenburgh [46,47] proposed an improvement over the method presented by Leeuwenburgh and Arts [45], where the use of the fast-marching method was extended for applicability on corner point grids, thereby improving the overall accuracy of the method. They presented their method in a simple 2D synthetic case and applied it to the Norne field using the ensemble smoother with multiple data assimilation (ES-MDA) [51], Appendix A.

In this paper, we focus on the simplicity and applicability of front re-parametrization methods and present a flow-based alternative to previous distance-to-front methods. For this, based on the work presented by Kretz [44], we propose a distance-to-front method using streamlines obtained from full physics flow simulation.

In the following section, we introduce the concept of distance-to-front measurement (Section 2.1), then we introduce fluid flow streamlines as a post-processing of flow simulation and propose a method for calculating distances to fronts using streamline information obtained from full-physics flow simulation (Section 2.2). We then present a set of numerical experiments on a 2D synthetic case (Section 3.1) and a realistic 3D case based on a real case scenario (Section 3.2), where the proposed method is applied and the obtained results are discussed. Finally, the main conclusions of the work are presented in Section 4.

## 2. Methodology

### 2.1. Distance-to-Front Measurement

Measuring the distance to a front requires front extraction, normally from a seismic attribute that can capture spatial changes in saturation or pressure over the subsurface. These changes can often refer to a timelapse confirmation on features related to the displacement of oil by water and/or gas, dissolution effects, or significant pressure changes. By posing the problem in this way, we assume that the shape or boundary of these features is enough to capture the relevant phenomena as opposed to using the original individual cell amplitudes of 4D seismic data. Added to that, reducing the information to a relevant shape or boundary can often be more informative and reliable while being an advantage in terms of computational cost and efficiency.

Normally, within the feature-based 4D history-matching domain, grid-based geometrical distance measurements are used (e.g., Euclidean, Hamming, fast-marching methods, Chain-Vese, Hausdorff, etc.). However, one should address the choice of measuring subsurface changes from a dynamic perspective, as subsurface changes derive from dynamic mechanisms related to fluid production and subsurface geology.

### 2.2. Using Streamlines for Distance-to-Front Calculation

In this regard, streamlines can be a solution to link both subsurface geology and fluid production. Streamlines and streamline-based properties are valuable tools for understanding reservoir connectivity and fluid flow patterns for large, heterogeneous models and can be easily obtained through the post-processing of full-physic simulations. In order to define a fluid flow streamline, three key properties are necessary: flow rate, time of flight, and a cell ID pointing to a given cell in the reservoir grid. The grid cell ID is used to extract relevant information from the grid and map solution variables between the streamline and the global numerical grid.
Streamline geometries (used for visualization purposes only) are then obtained through the Pollock method [52,53], starting from a single cell of the model with the calculation of a flow rate calculated for each of the cell faces (assumed to be uniform along the faces). Along with pressure, the total flow rate in and out of each of the faces can then be calculated based on the total Darcy velocity:

$$ \vec{v}_t = -\lambda_t \vec{\nabla} P_0 + (\sum j \lambda_j \rho_j) \vec{g}, \qquad (1) $$

where $\vec{v}$ is the flow velocity, $\lambda$ is the phase mobility, $\vec{\nabla} P$ is the pressure gradient, $\rho$ is the mass density, and $\vec{g}$ is the acceleration due to gravity.

To conform with the orthogonal grid assumptions of the Pollock method (Figure 1), an isoparametric transformation is applied to all of the grid cells of the model onto a unit cube grid [54]. Consequently, for visualization purposes, this representation can be transformed back into the spatial coordinate system (e.g., corner point gridding format). This way, by being able to obtain a set of streamlines branching all cells of the model where flow is occurring from a given source to a sink, we can use streamlines as a flow-based real coordinate system medium to measure the distance between the true location of water-saturated fronts and simulated saturation fronts obtained from the flow simulation of candidate models.

<figure><img src="image_4.png" /><figcaption>Figure 1. Streamline tracing scheme where a single streamline enters the Y0 face of a single cell and exits the X1 face through the exit point. The curvature along the cell is guided by the calculated in/out flow velocities of all faces of the cell.</figcaption></figure>

As shown in Figure 2a, a simple 2D five-spot reservoir model was used to illustrate the method. Two fronts are presented, one originating from either the inversion or interpretation of 4D seismic data (observed front, in black) and a second one representing the same front on a given candidate model after flow simulation (simulated front, in grey). By computing the fluid-flow streamlines over the simulation period, we obtained a set of streamlines that connects both fronts. This way, it became possible to obtain the distance between both features using flow-based distance measurement supported by streamlines and to use the computed distance vector as an innovation in the history-matching workflow. The same concept is shown for cases where the resulting simulated flow patterns were obtained using an anisotropic permeability field (Figure 2b) and for the example of a more heterogeneous geological scenario representing a channelized structure (Figure 2c). From Figure 2c, it also becomes clear that by applying this method, we gained added advantage of capturing the distance between both fronts along the geological representation of our model (grid geometry and petrophysical properties driving fluid flow) as opposed to using a purely geometrical approach where such detail might be lost. A distance value of zero
was assigned for all positions where both observed and measured front data coincided; for the remaining unmatched locations, an added physical meaning was obtained by assigning positive distances for locations where the measured front was ahead of the observed front and vice versa.

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Streamline tracing with observed and predicted waterfronts for different model geologies at the same time step. The observed fronts are indicated in black, simulated fronts in grey and production streamlines in red. (a) streamlines obtained in an isotropic model, (b) streamlines obtained in an anisotropic model, (c) streamlines obtained in a model representing a channelized structure.</figcaption></figure>

Under cases where geology is uncertain (e.g., multiple different geological scenarios), the same distance metric can be used. For this case, we considered an uncertain geological scenario where we took two different channelized structures, honoring the well data but with different orientations.

Figure 3 illustrates the application of the proposed method, where we evaluated the 4D seismic distance between both model propositions on a given monitor date.

To obtain the waterfronts at a given time $t_f$ (Figure 3b,e), we calculated the difference between the water saturation at $t_f$ (Figure 3a,d) with the initial water saturation at $t_0$. A threshold was applied to binarize the information into drained/undrained regions, and finally, we extracted the contour of the drained region (a perimeter in 2D and a surface in 3D) (Figure 3c). Since different candidate models have different simulated front shapes (e.g., in 2D having smaller or larger perimeters over a different number of model cells), distances must be computed from the location of the observed front, which should be sourced from seismic 4D interpretation and expertise related to the mechanisms driving the production of a given reservoir. This way, the observed front location, comprising a set of cells to which the distances are to be computed, is input into the algorithm. The final innovation vector is calculated by gathering all the distance measurements obtained from the set of shortest distance paths provided by computed streamlines connecting the observed front to the simulated front. For illustration purposes, Figure 4 is a spatial representation of the obtained innovation vector on top of the observed front location.

The same calculation can be carried out for different models, with the result being a vector of innovations with the size of the number of cells on the observed front. By being able to obtain sets of equally sized innovation vectors for a variety of models, the proposed method can be easily integrated into any ensemble-based assimilation approach. Algorithm 1 summarizes the use of the innovation vector in the context of an ensemble-based history-matching workflow.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Example of front extraction (c) for a true model (a) and a candidate model (d), given the bi-narized $\Delta S_W$ obtained from both solutions after thresholding (b,e). The colormap for (a,d) represents $S_W$ on the monitor date (dimensionless).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Distance between observed and simulated fronts streamlines represented over all cells of the observed saturation front (i and j indexes are grid cells, and colorbar units are m).</figcaption></figure>
<table><caption>**Algorithm 1:** Four-dimensional distance-to-front using streamlines</caption><thead><tr><th></th><th></th></tr></thead><tbody><tr><td>FOR EACH ENSEMBLE UPDATE STEP</td><td></td></tr><tr><td>1.</td><td>Input array(s) of cells representing observed front(s) in the model grid.</td></tr><tr><td>2.</td><td>Full physics flow simulation of all candidate members of the ensemble.</td></tr><tr><td>3.</td><td>Calculate the difference between saturation data from the date of the seismic survey(s) to t<sub>0</sub>.</td></tr><tr><td>4.</td><td>Binarize the output of step 3 into flooded/non-flooded regions according to a threshold.</td></tr><tr><td>5.</td><td>Calculate the contour of the flooded region at the date of seismic survey(s).</td></tr><tr><td>6.</td><td>Matching locations for the output of step 5 and step 1 are assigned a distance of 0.</td></tr><tr><td>7.</td><td>Post-processing of streamlines coming from step 2.</td></tr><tr><td>8.</td><td>Extract the shortest distance given by the streamlines (step 7) connecting the observed front (step 1) with the simulated front (step 5).</td></tr><tr><td>9.</td><td>Merge array(s) of distances computed at step 6 and step 8.</td></tr><tr><td>CONTINUE TO THE NEXT UPDATE STEP</td><td></td></tr></tbody></table>

## 3. Results and Discussion

### 3.1. Synthetic 3D Case

In this section, we resort to a simple 3D channel reservoir to showcase the performance of the proposed method for 4D seismic history matching using streamlines for distance-to-front measurement. The reservoir comprises an anticlinal trap with three facies, containing a North-West to South-East trending leaking normal fault with a throw towards North-East. It has an extension of 5900 × 3800 m with a thickness of 52 m and is discretized into a 60 × 39 × 5 grid with 11,682 active cells in total. The used model attempted to represent the flow of two immiscible phases (oil and water) with a connate water saturation of 0.15, residual oil saturation of 0.15, and initial formation pressure of 400 bar. Out of the three facies in the model, two were permeable, and one acted as a horizontal barrier in the middle layer of the reservoir. The model had two injectors and two producing wells located in the best-quality facies. The true model was randomly sampled from an ensemble of realizations generated through truncated Gaussian simulation (TGS) [55,56] and conditioned to the information at well locations. The average permeability was 650 mD and 150 mD in the best and background facies, respectively. The average porosity in the best facies was 0.3, and in the background facies, it was 0.18. Field production spanned a total of 20 years, with production occurring through the intervention of two producing wells located on each side of the normal fault and an injector well located in the south.

The considered uncertain parameters in the reservoir were facies location, populated by the respective spatial distributions of porosity, the permeability and net-to-gross ratio (NTG) of each of the facies, and an ensemble of 100 realizations used to sufficiently avoid sampling errors and rank deficiency in the updated procedure.

True front positions for dates were obtained by running a flow simulation on a model outside of the ensemble until, respectively, 8 and 16 years after the initial production date. The resulting saturation maps were then processed in order to extract waterfront positions through the binarization of the saturation differences between monitor dates and initial condition, with a threshold (a threshold of 0.02 was used for this case) to obtain a swept region and interpretation of the front position over the binarized swept region (Figure 5). The Bayesian formalism on data assimilation problems normally requires the likelihood function to be responsible for assigning the weighting of data mismatch terms [33]. Regardless, there are examples in the literature of 4D history-matching attempts where observation errors were selected according to what the authors believed to be an "acceptable" result in terms of match quality [47,57–59]. In our case, we considered the interpreted front position to have an uncorrelated error with a standard deviation of 150 m, which consisted of approximately the size of three grid blocks along the XY direction. The resulting front positions for all seismic monitor dates, along with the historical production
data, were used in an assisted history-matching scheme using ES-MDA with five iterations and 100 candidate models per iteration to update the underlying uncertain parameters towards matching both production and seismic data. For all candidate models, we resorted to a full physic reservoir simulator to run fluid flow simulation and obtain saturation data on monitor dates as well as average bottom-hole-pressure (WBPC3), flow rates (oil and water), and water cut. A set of three experiments were run to show the capabilities of the proposed method. We considered a scenario where only production data were used for history matching (NO4D), a scenario where the Hausdorff distance [37] was used to measure waterfront position mismatches (4DHDF), and finally, our proposed method, using streamlines for distance-to-front measurement (4DSLN).

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Observed seismic information for the ground truth model (a) and initial ensemble grid block average of swept area (b), where the spread of initial candidate solutions can be observed. Colormap represents block average binarized $\Delta S_W$ measured over the initial ensemble (dimensionless).</figcaption></figure>

Figure 5a shows the observed seismic data in a binarized form (initial minus monitor date saturation) for one of the provided seismic monitors of the experiment and also the true front that we intended to match. Figure 5b shows the distribution of the simulated swept regions provided by the same binarized saturation information but averaged over all grid blocks for all members of the ensemble. We observed that there was a significant spread of the simulated swept regions (blue to red variation, Figure 5b), especially to the north of the reservoir where the model was less swept at this stage of the production schedule.

The synthesized results of the three ran experiments are presented in Figure 6, where we can observe the match of the swept regions at the end of the update procedure. From the obtained results, we can see that the absence of additional seismic information of the NO4D method provided a poor match concerning the observed seismic data (Figure 6b). The initial and final spreads of the swept regions also did not vary substantially, although we arrived at a fairly good match concerning the production data (Figure 7a). This reinforces the idea that seismic data integration is crucial for history matching as, while the solutions
provided might be considered moderately accurate in terms of matching well production data, they are, in fact, inaccurate in terms of matching real waterfront movement along the reservoir as the reservoir is being produced, greatly reducing any forecasting capabilities of the model. For the case of the experiment using 4DHDF, the final solution was a close match when compared to the observed seismic information (Figure 6c). This was also accompanied by a very good match in production data (Figure 7b). However, for cases where the grid cell count is higher, the applicability of the method is debatable. In fact, we were able to run the experiment on this simplistic and small 3D model, but the method became more unfeasible with increased model sizes at operational levels (often in the order of millions of cells). Under such conditions, calculating the ensemble update can rapidly become computationally intractable as it requires the inversion of a square matrix at least the size of the grid. Regarding the proposed method (4DSLN), we observed that the final match of the swept regions was nearly identical to the one obtained by the 4DHDF method (Figure 6d), with the same being observed on the match of production data (Figure 7c).

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Results obtained from running the three experiments in terms of the final swept area. Observed seismic information for the ground truth model (a) and final ensemble grid block average of swept areas for NO4D (b), 4DHDF (c), and 4DSLN (d). Colormap represents block average binarized $\Delta S_W$ measured over the final ensemble (dimensionless).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Evolution of production match quality from initial (blue) to final (red) ensembles, along with the available observed production data and associated uncertainty (black error bars). NO4D (a), 4DHDF (b), and 4DSLN (c).</figcaption></figure>
Figure 8 further explores the obtained results by showing the difference (grid block average) obtained between the swept regions of the final ensemble and the ground truth for all three ran experiments. We can again see that not using seismic information renders the final solution far from the truth (Figure 8a), whereas using seismic data provides a much closer match (Figure 8b,c). Moreover, we can see that using the proposed streamlined distance-to-front method provides very similar results when compared to 4DHDF, the advantage of using a substantially reduced amount of data (Figure 9).

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Grid block difference between the obtained swept areas of the final ensemble and the ground truth. NO4D (a), 4DHDF (b), and 4DSLN (c). Units are dimensionless.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Boxplot for RMSE of all production data mismatches (P1, P2 with WWCT, WWPR, and WBPC3). The blue bounding boxes represent the 25 and 75% quantiles; the whiskers are the extremes, the red line in the box is the median, the black dot is the mean, and the red plus signs are the outliers.</figcaption></figure>

Regarding the production match of the field-producing wells (two producers and one injector), Figure 8 shows the final obtained results in terms of well water cut (WWCT) for producer P2. Well P2 can be taken as a sufficient example as in all cases, the same match quality was observed for all present wells in the model. We can see that in all ran experiments, a reduction of the production ensemble spread was observed. As mentioned previously, and as can be observed in Figure 8, not using seismic data can potentially lead to a final ensemble where production data matching is improved; however, the inclusion of seismic information drastically increases the quality of the final match. In fact, both experiments using seismic data information for the update procedure provided very good matches with similar qualities when coupled with the 4DSLN method, obtaining residuals improved results when compared to 4DHDF.
The results for all ran experiments are further demonstrated in Figure 9, depicting the root mean square error (RMSE) defined as follows:

$$RMSE = \sqrt{\frac{\sum_{i=1}^{n} (\hat{y}_i - y_i)^2}{n}}, \quad (2)$$

where n is the number of observations, and ŷ and y are, respectively, the predicted and observed values. Figure 9 shows the minimization of the obtained production data mismatch for all run experiments. Immediately, we can observe that although the incorporation of production data considerably reduced production data mismatch relative to the prior ensemble (NO4D), the added value for the inclusion of seismic data was reflected in the ability to further improve predictions (4DHDF and 4DSLN). We also observed that the final production match quality was similar for both methods using seismic information.

In order to examine the quality of the updates to the facies model, we compared the obtained solutions to the elected ground truth model. In Figure 10, we show the probability map for good facies (shown in dark red) in layer 3 of the reservoir for the initial and final ensembles. Immediately, we can observe the potential for the inclusion of seismic data, as the experiment that used only well production data to constrain the update (NO4D) was unable to capture the spatial distribution patterns of the true model. Furthermore, the convergence towards the optimum solution is hardly visible, as almost no change occurred between the initial and final ensembles. This information, when considering the somewhat satisfactory results for production mismatch shown in Figures 6b and 8a for the same experiment, serves as a good example to stress the importance of using seismic information, illustrating a case of non-uniqueness observed in history matching. In fact, for the NO4D experiment, while many solutions may adequately fit the production data, the lack of seismic information to further constrain the update procedure will render the final ensemble of models less reliable for forecasting future production data. On the other hand, we can observe that 4DHDF and 4DSLN are able to arrive at a final ensemble that better resembles the true model. Using a full grid distance measurement, 4DHDF was able to accurately capture the spatial distribution of the facies locations in the final ensemble. Similarly, 4DSLN was able to arrive at a final ensemble that captured the most relevant spatial distribution patterns for fluid flow prediction, namely the connection of the North-East and South-East regions to the center of the model, with the advantage of using a much lower amount of seismic data for this purpose.

## 3.2. Realistic 3D Case

For the final example, we applied the proposed method to a realistic 3D case based on a real field, where the application of 4DHDF could not be achieved due to the large amount of data. The reservoir is a turbiditic depositional environment located in offshore Africa. The grid size was 194 × 203 × 48 with cell dimensions of 50 × 50 × 3 m on average over the i, j, and k directions. The field was produced by eight producing wells, and five injector wells provided pressure support. Twelve years of historical oil and water production, as well as bottom-hole pressure, were available. The model had 13 different flow units (Figure 11) populated by a total of three different facies types with different spatial continuity patterns (realizations are obtained using the TGS algorithm), with each facies type having specific petrophysical property distributions and dynamic parametrizations according to the quality of the rock they represented (poor to good sands). The reservoir was compartmentalized by a total of 25 faults with varying transmissibilities.

For our experiment, we considered both geological and engineering uncertain parameters. Facies locations over the different flow units were considered to be uncertain geologic parameters, and a set of connectivity regulating parameters (fault transmissibilities, sedimentological and aquifer connection, and region transmissibilities), pore volume and productivity multipliers, were assumed to be uncertain engineering parameters, 28 in
total. Table 1 presents a summary of the uncertain parameters as well as their assumed uncertainty bounds.

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Comparison of match estimates for facies location for all experiments.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 11.** Top view (**left**) and 3D view (**right**) of the flow units for the realistic 3D case. Different colors represent different flow units with different rock properties and spatial continuity patterns, 13 in total.</figcaption></figure>
<table><caption>**Table 1.** Uncertain parameters and their associated uncertainty bounds.</caption><thead><tr><th>Parameter Type</th><th>Count</th><th>Minimum</th><th>Maximum</th></tr></thead><tbody><tr><td>Pore Volume Multipliers</td><td>13</td><td>0.85</td><td>1</td></tr><tr><td>Fault Transmissibility Multipliers</td><td>10</td><td>1 × 10<sup>-6</sup></td><td>1</td></tr><tr><td>Region Transmissibility Multipliers</td><td>3</td><td>1 × 10<sup>-6</sup></td><td>1</td></tr><tr><td>Productivity Index Multipliers</td><td>2</td><td>0.001</td><td>1</td></tr><tr><td>Total</td><td>28 Parameters</td><td></td><td></td></tr></tbody></table>

We generated an initial ensemble of 100 equiprobable candidate models using TGS to model facies locations based on variograms estimated from well log data. To sample over the uncertain engineering parameters, we resorted to Latin hypercube sampling based on the available prior knowledge. An extra model, elected as the ground truth, was also generated based on the same geological and engineering prior assumption. The location of the true seismic fronts at both monitor dates was assumed to be obtained from a standard seismic inversion procedure followed by additional processing and interpretation and corresponding to the information obtained after 6.2 and 9.7 years of fluid flow simulation on the ground truth model. The control parameters for the fluid flow simulation were set to operate the producers at historic reservoir volume rates and the injectors at historic fluid rates. The interpreted front locations were assumed to have uncorrelated errors with a standard deviation of 200 m (around four grid blocks).

Figure 12 shows the initial water saturation on a randomly selected model from the initial ensemble and the corresponding streamlines obtained on both monitor dates. The interpreted waterfronts for the ground truth model are also superimposed on the image.

<figure><img src="image_15.png" /><figcaption>**Figure 12.** Top view of the initial $S_{w}$ (blue—0 to red—1) for a realization of the realistic 3D. Well locations (circles with black borders) and two arbitrarily interpreted fronts (dashed black lines) extracted at two different seismic monitoring dates (a). Post-processed fluid flow streamlines (colors represent the time of flight) obtained at the dates of the first and second seismic monitors (b,c).</figcaption></figure>

In the same randomly selected model from the initial ensemble, using a threshold to binarize the change in saturation at the seismic monitor dates, we obtained the swept regions on each monitor date. Figure 13 shows the swept regions and front locations of the same randomly selected model on both monitor dates against the ones obtained from the ground truth. The mismatch between the observed and simulated fronts was observed, resulting from excessive water sweeping over the model.
<figure><img src="image_16.png" /><figcaption>**Figure 13.** Top view of the front mismatch at monitor date one (a) and two (b). The property shown (dark gray) is the simulated response in terms of binarized $\Delta S_w$ along with the associated simulated front (red) obtained for this realization.</figcaption></figure>

We set up the experiment using an ES-MDA update scheme of the 100 members of the initial ensemble over five iterations. Figure 14 compares, at the first monitor date and on layer 12 of the model, the true front location with a random unmatched realization from the initial ensemble (center) and the best solution found on the final ensemble. We can observe that the final best-matching solution is a clear improvement when compared to the initial guess example. Further, we superimposed the water cut ratio of a producing well (PROD 4) associated with the front arrival at the date of the seismic acquisition. The best matching realization (shown in red) was able to match the water cut ratio perfectly. We can also see that the front arrival time, closely related to the date of the seismic monitor (red dashed line), matches the water breakthrough of the well, meaning we successfully matched the correct time of waterfront arrival and the correct water volume production.

<figure><img src="image_17.png" /><figcaption>**Figure 14.** (a) Ground truth $S_w$ along with the associated true front (black) for the first monitor date. (b) Example of simulated $S_w$ of a randomly selected model from the initial ensemble, highlighting an excessive water sweep going over the front (red arrows). (c) Final best-matched model of $S_w$ along with WWCT production curves for PROD4 (c-top). Injectors producers are represented respectively in blue and red circles.</figcaption></figure>

The same match quality was observed for the second monitor date (Figure 15) on a different layer of the reservoir (layer 14).

By comparing the swept regions and front location of the ground truth with the best-obtained solution, we observed a very good match concerning the front location and the production data of an associated well at the true front location (PROD 6).
<figure><img src="image_18.png" /><figcaption>**Figure 15.** (a) Ground truth $S_w$ along with the associated true front (black) for the second monitor date. (b) Final best-matched model of $S_w$ along with WWCT production curves for PROD6 (b-top). Injectors producers are represented respectively in blue and red circles.</figcaption></figure>

A summary of the run, in terms of front location matches, from initial to final ensemble is illustrated in Figure 16. We can see the average and standard deviations of the distances measured to the front at every grid block point where the true front is located (up). We can also observe that before and after history matching, the final obtained front locations much more closely matched the truth (as they are closer to 0) and that, at the same time, the uncertainty over the final ensemble of solutions was substantially reduced.

<figure><img src="image_19.png" /><figcaption>**Figure 16.** Summary of the front matches for a given layer of the model at the first (a) and second (b) monitor dates with the grid block mean and standard deviations of distance-to-front obtained in the initial and final ensembles (top). A plot of the most representative initial and final waterfront locations</figcaption></figure>
in regard to the mean distances measured over the ensemble is represented in gray shade. In the bottom, the initial and final grid block average (red and black dots) and standard deviation (red and grey area) of distance-to-front obtained for all layers of the model where a front was interpreted.

Regarding the convergence on uncertain parameters over the run, Figure 17 shows the probability of facies locations from the initial to final ensemble, compared to the ground truth. We can immediately observe that, for this particular model, the most predominant facies type is facies 1 (in blue), taking the bulk of the task of governing fluid flow production and reservoir connectivity. The final ensemble shows the convergence towards the true locations of the three facies types, but most importantly, the correct spatial continuity and connectivity patterns that can be observed in the true model.

<figure><img src="image_20.png" /><figcaption>**Figure 17.** Summary of the probability of rock type obtained from the run experiment. The initial (top) and final (middle) matches of facies type (blue, green, red) match are compared against the ground truth (bottom).</figcaption></figure>

An example of the convergence that can be observed on selected uncertain engineering parameters for the experiment is shown in Figure 18. During the course of the run, we observed the convergence of the parameter distributions towards true values at the same time that the misfit was also being reduced.

The same can be seen for the production plots of water cut and bottom-hole pressure for all producing wells of the model (Figure 19). We can clearly observe the gradual reduction of the mismatch from the initial to the final iteration, with the final ensemble providing good match quality over all wells.
<figure><img src="image_21.png" /><figcaption>Figure 18. Scatter plot of the parameter value (x-axis) vs. misfit (y-axis) for the convergence of a selected set of parameters (true value in green circle) from the initial (blue) to the final iteration (red).</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Figure 19. Run evolution in terms of production match for WWCT (first and second rows) and WBPC3 (third and fourth rows) from the initial (blue) to the final iteration (red), along with the observations (black error bars) and acquisition date of used seismic monitors (dashed green lines).</figcaption></figure>
## 4. Conclusions

A novel parametrization scheme was proposed to quantify the differences between simulated and measured seismic data in terms of distances between fluid front positions measured by streamlines. The methodology was presented with the help of a simplistic 2D case and applied both on a synthetic and realistic 3D case with seismic and production data, highlighting the advantages and good performance of the proposed scheme when compared to other commonly used schemes. Although the presented cases were of synthetic and semi-synthetic nature, the application of the proposed method to real-case scenarios should not inherently present any disadvantages. Naturally, the success in achieving positive outcomes, as with all methods, hinges on the quality of the models and the data at hand.

The obtained results show that the proposed parametrization can achieve similar results when compared to other methods resorting to a reduced amount of data. The streamlined information used for distance measurement can be easily obtained as post-processing of standard full-physics simulation outputs.

Despite the advantages of using streamlines for distance measurement, there could be some potential disadvantages when processing a large number of streamlines, as could be the case with very large models. However, these challenges can be mitigated through parallelization methodologies or, for example, by refining the way streamlines are post-processed. For instance, focusing the processing on the cells where fluid fronts are located can significantly speed up the process. Some limitations in terms of precision may potentially arise when the models, and consequently the simulated fronts, are too far from the ground truth. This is due to the inherent requirement of obtaining streamlines that intersect both fronts in order to accurately compute distance data through them. Despite this, any missing information can be easily complemented by other distance metrics.

The capabilities of the method were not only showcased on a simplistic 3D model with simple reservoir conditions and production schedules but also on a more complex and realistic 3D model, where the Hausdorff distance method becomes computationally intractable, having obtained encouraging results both on fluid front and production data match. The method is expected to be effectively applicable to a range of reservoir types beyond those demonstrated in this study. Furthermore, the parametrization scheme is able to be adaptable and functional across diverse geological settings.

While the application of the innovative parametrization scheme was showcased within the framework of assisted history matching using an ensemble history matching methodology (ES-MDA), its application can be easily extended to any other history-matching workflow or even to different domains of application other than hydrocarbon exploration, e.g., CO₂ monitoring.

Future areas of research may involve advancements in how streamlines are calculated to be fit-for-purpose for the type of methodology presented here. This may include methods relying on the focusing of information on regions pertaining to both the simulated and real fronts. Another avenue of investigation that could prove interesting is implementing localization techniques using information derived from streamlines. This approach could potentially enhance the method's performance by refining the accuracy and efficiency of the parameterization scheme, especially in complex geological settings where the alignment of simulated and real fluid fronts is critical. These developments could lead to more nuanced and effective history matching in reservoir simulation, thus broadening the scope and applicability of the proposed method.

## 5. Patents

Berthet, P., Trani, M. [60]. A method for obtaining at least one physical property of a subsurface volume of a hydrocarbon reservoir over time (European priority Application, filing date 4 December 2020, publication number EP4009086).
**Author Contributions:** Conceptualization, E.B., P.B. and M.T.; methodology, E.B., P.B. and M.T.; validation, E.B., P.B. and O.T.; formal analysis, E.B.; investigation, E.B.; resources, P.B.; data curation, E.B.; writing—original draft preparation, E.B.; writing—review and editing, E.B., P.B., M.T., O.T. and C.L.; visualization, E.B.; supervision, P.B., O.T. and C.L.; project administration, P.B.; funding acquisition, P.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by TotalEnergies S.E. and the ANRT—Association Nationale de la Recherche et de la Technologie, grant number 2020/0558.

**Data Availability Statement:** Data are contained within the article.

**Acknowledgments:** The authors acknowledge the support from the “Association Nationale de la Recherche et de la Technologie” for funding this project and colleagues from the TotalEnergies S.E. affiliate for the valuable discussions and for providing the realistic 3D dataset used in this work.

**Conflicts of Interest:** The authors declare no conflict of interest.

## Appendix A

First introduced by Emerick and Reynolds [51,52] as an extension to the standard ES method, the multiple data assimilation (MDA) method was developed with the purpose of enhancing the performance of the ensemble Kalman filter [38] and ensemblesmoother [58] under nonlinear conditions. In this iterative approach, all data are assimilated multiple times by applying an inflation parameter $\alpha_k$ to the covariance matrix of measurement errors. This was proven to be equivalent to the single data assimilation case for linear-Gaussian systems, given that the measurement error covariance matrix is appropriately scaled [51,52].

Therefore, considering the standard ES update equation, the inflation parameter $\alpha_k$ is included as follows:

$$m_j^{k+1} = m_j^k + \hat{C}_{md}^k (\hat{C}_{dd}^k + \alpha_k C_d)^{-1} (d_{obs} - g(m_j^k) + e_j), \quad (A1)$$

where $k + 1$ and $k$ are the indexes of the iterative procedure, where $k = 1, 2, \dots, N_a$ with $N_a$ being the total number of assimilations; $\hat{C}_{md}^f$ is the cross-covariance matrix between the a priori vector of model parameters, $m^f$, and the predicted data vector, $g(m^f)$; $\hat{C}_{dd}^f$ is the covariance matrix of the predicted data of size $N_d \times N_d$; $d_{obs}$ is the vector of observed data; $e_j$ is the perturbation vector added to the observed data; $C_d$ is the covariance matrix of the observed data error of size $N_d \times N_d$.

To implement ES-MDA, the values of the inflation factor $\alpha_k$ in each iteration need to be defined. The necessary condition for choosing the inflation factor is the following:

$$\sum_{k=1}^{N_a} \frac{1}{\alpha_k} = 1. \qquad (A2)$$

A lack of consensus exists on how to choose the values of $\alpha_k$, with $\alpha_k = N_a$ being commonly used. On the same note, the number of iterations $N_a$ must be set beforehand, with generally accepted values ranging from 4 to 10 iterations. Improved versions of the method have been proposed where $\alpha_k$ and $N_a$ are chosen automatically [59]. For the present work, the authors opted to use $\alpha_k = N_a = 5$.

## References

1. Landrø, M. Discrimination between pressure and fluid saturation changes from time-lapse seismic data. *Geophysics* **2001**, 66, 836–844. [CrossRef]

2. Lumley, D.E. Time-lapse seismic reservoir monitoring. *Geophysics* **2001**, 66, 50–53. [CrossRef]

3. Thore, P.; Hubans, C. 4D seismic-to-well tying, a key step towards 4D inversion. *Geophysics* **2012**, 77, R227–R238. [CrossRef]

4. Thore, P.; Blanchard, T.D. 4D propagated layer-based inversion. *Geophysics* **2015**, 80, R15–R29. [CrossRef]

5. Grana, D.; Mukerji, T. Bayesian inversion of time-lapse seismic data for the estimation of static reservoir properties and dynamic property changes. *Geophys. Prospect.* **2014**, 63, 637–655. [CrossRef]
6. Maharramov, M.; Biondi, B.L.; Meadows, M.A. Time-lapse inverse theory with applications. *Geophysics* **2016**, *81*, R485–R501.
[CrossRef]

7. Kazemi, A.; Stephen, K.D.; Shams, A. Seismic History Matching of Nelson Using Time-Lapse Seismic Data: An Investigation of 4D Signature Normalization. *SPE Reserv. Evaluation Eng.* **2011**, *14*, 621–633. [CrossRef]

8. Maleki, M.; Davolio, A.; Schiozer, D.J. Qualitative time-lapse seismic interpretation of Norne Field to assess challenges of 4D seismic attributes. *Lead. Edge.* **2018**, *37*, 754–762. [CrossRef]

9. Lygren, M.; Husby, O.; Osdal, B.; El Ouair, Y.; Springer, M. History matching using 4D seismic and pressure data on the Norne field. In Proceedings of the 67th EAGE Conference & Exhibition, Madrid, Spain, 13–16 June 2005; European Association of Geoscientists & Engineers: Utrecht, The Netherlands, 2005; p. cp-1. [CrossRef]

10. Roggero, F.; Ding, D.Y.; Berthet, P.; Lerat, O.; Cap, J.; Schreiber, P.-E. Matching of Production History and 4D Seismic Data--Application to the Girassol Field, Offshore Angola. In Proceedings of the SPE Annual Technical Conference and Exhibition, Anaheim, CA, USA, 11–14 November 2007; OnePetro: Richardson, TX, USA, 2007. [CrossRef]

11. Castro, S.A.; Caers, J.; Otterlei, C.; Meisingset, H.; Hoye, T.; Gomel, P.; Zachariassen, E. Incorporating 4D seismic data into reservoir models while honoring production and geologic data: A case study. *Lead. Edge.* **2009**, *28*, 1498–1505. [CrossRef]

12. Le Ravalec, M.; Tillier, E.; Da Veiga, S.; Enchery, G.; Gervais, V. Advanced integrated workflows for incorporating both production and 4d seismic-related data into reservoir models. *Oil Gas Sci. Technol. Rev. d'IFP Energies Nouv.* **2012**, *67*, 207–220. [CrossRef]

13. Roggero, F.; Lerat, O.; Ding, D.; Berthet, P.; Bordenave, C.; Lefeuvre, F.; Perfetti, P. History matching of production and 4D seismic data: Application to the Girassol field, offshore Angola. *Oil Gas Sci. Technol. Rev. d'IFP Energies Nouv.* **2012**, *67*, 237–262.
[CrossRef]

14. Byerley, G.; Singer, L.; Rose, P. Resaturated pay: A new infill target type identified through the application and continuous improvement of 4D seismic at the Forties Field. *Lead. Edge.* **2016**, *35*, 831–838. [CrossRef]

15. Calvert, M.A.; Hoover, A.R.; Vagg, L.D.; Ooi, K.C.; Hirsch, K.K. Halfdan 4D workflow and results leading to increased recovery. *Lead. Edge.* **2016**, *35*, 840–848. [CrossRef]

16. Rankey, E.C.; Mitchell, J.C. That's why it's called interpretation: Impact of horizon uncertainty on seismic attribute analysis. *Lead. Edge.* **2003**, *22*, 820–828. [CrossRef]

17. Zhou, W.; Lumley, D. Nonrepeatability effects on time-lapse 4D seismic full-waveform inversion for ocean-bottom node data. *Geophysics* **2021**, *86*, R547–R561. [CrossRef]

18. Sarkar, S.; Gouveia, W.P.; Johnston, D.H. On the inversion of time-lapse seismic data. In Proceedings of the 2003 SEG Annual Meeting, Dallas, TX, USA, 26–31 October 2003; OnePetro: Richardson, TX, USA, 2003. [CrossRef]

19. Buland, A.; El Ouair, Y. Bayesian time-lapse inversion. *Geophysics* **2006**, *71*, R43–R48. [CrossRef]

20. Suman, A.; Fernández-Martínez, J.L.; Mukerji, T. *Joint Inversion of Production and Time-Lapse Seismic Data: Application to Norne Field*; Stanford University: Stanford, CA, USA, 2013. [CrossRef]

21. Alvarez, E.; MacBeth, C.; Brain, J. Quantifying remaining oil saturation using time-lapse seismic amplitude changes at fluid contacts. *Pet. Geosci.* **2016**, *23*, 238–250. [CrossRef]

22. Arenas, E.; van Kruijsdijk, C.; Oldenziel, T. Semi-automatic history matching using the pilot point method including time-lapse seismic data. In Proceedings of the SPE Annual Technical Conference and Exhibition, New Orleans, LA, USA, 30 September–3 October 2001; OnePetro: Richardson, TX, USA, 2001. [CrossRef]

23. Fagervik, K.; Lygren, M.; Valen, T.S.; Hetlelid, A.; Berge, G.; Dahl, G.V.; Sønneland, L.; Lie, H.E.; Magnus, I. A method for performing history matching of reservoir flow models using 4d seismic. In Proceedings of the 2001 SEG Annual Meeting, San Antonio, TX, USA, 9–14 September 2001; OnePetro: Richardson, TX, USA, 2001. [CrossRef]

24. Gosselin, O.; Berg, S.v.D.; Cominelli, A. Integrated history-matching of production and 4D seismic data. In Proceedings of the SPE Annual Technical Conference and Exhibition, New Orleans, LA, USA, 30 September–3 October 2001; OnePetro: Richardson, TX, USA, 2001. [CrossRef]

25. Gosselin, O.; Aanonsen, S.I.; Aavatsmark, I.; Cominelli, A.; Gonard, R.; Kolasinski, M.; Ferdinandi, F.; Kovacic, L.; Neylon, K. History matching using time-lapse seismic (HUTS). In Proceedings of the SPE Annual Technical Conference and Exhibition, Denver, CO, USA, 5 October 2003; OnePetro: Richardson, TX, USA, 2003. [CrossRef]

26. van Ditzhuijzen, R.; Oldenziel, T.; van Kruijsdijk, C. Geological parameterization of a reservoir model for history matching incorporating time-lapse seismic based on a case study of the Statfjord field. In Proceedings of the SPE Annual Technical Conference and Exhibition, New Orleans, LA, USA, 30 September–3 October 2001; OnePetro: Richardson, TX, USA, 2001. [CrossRef]

27. Dong, Y.; Oliver, D.S. Quantitative Use of 4D Seismic Data for Reservoir Description. *SPE J.* **2005**, *10*, 91–99. [CrossRef]

28. Haverl, M.; Aga, M.; Reiso, E. Integrated Workflow for Quantitative Use of Time-Lapse Seismic Data in History Matching–A North Sea Field Case (SPE94453). In Proceedings of the 67th EAGE Conference & Exhibition, Madrid, Spain, 13–16 June 2005; European Association of Geoscientists & Engineers: Utrecht, The Netherlands, 2005; p. cp-1. [CrossRef]

29. Portella, R.C.M.; Emerick, A.A. Use of Quantitative 4D-Seismic Data in Automatic History Match. In Proceedings of the SPE Latin American and Caribbean Petroleum Engineering Conference, Rio de Janeiro, Brazil, 20 June 2005; OnePetro: Richardson, TX, USA, 2005. [CrossRef]

30. Stephen, K.D.; Soldo, J.; Macbeth, C.; Christie, M.A. Multiple model seismic and production history matching: A case study. *SPE J.* **2006**, *11*, 418–430. [CrossRef]
31. Dadashpour, M.; Kleppe, J.; Landro, M. Porosity and permeability estimation by gradient based history matching using time-lapse seismic data. In Proceedings of the SPE Middle East Oil and Gas Show and Conference, Red Hook, NY, USA, 11–14 March 2007; OnePetro: Richardson, TX, USA, 2007. [CrossRef]

32. Luo, X.; Bhakta, T.; Jakobsen, M.; Nævdal, G. Efficient big data assimilation through sparse representation: A 3D benchmark case study in petroleum engineering. *PLoS ONE* **2018**, *13*, e0198586. [CrossRef]

33. Oliver, D.S.; Fossum, K.; Bhakta, T.; Sandø, I.; Nævdal, G.; Lorentzen, R.J. 4D seismic history matching. *J. Pet. Sci. Eng.* **2021**, *207*, 109119. [CrossRef]

34. Trani, M.; Arts, R.; Leeuwenburgh, O. Seismic history matching of fluid fronts using the ensemble Kalman filter. *SPE J.* **2012**, *18*, 159–171. [CrossRef]

35. Rollmann, K.; Soriano-Vargas, A.; Almeida, F.; Davolio, A.; Schiozer, D.J.; Rocha, A. Convolutional Neural Network Formulation to Compare 4-D Seismic and Reservoir Simulation Models. *IEEE Trans. Syst. Man, Cybern. Syst.* **2021**, *52*, 3052–3065. [CrossRef]

36. Tillier, E.; Le Ravalec, M.; Da Veiga, S. simultaneous inversion of production data and seismic attributes: Application to a synthetic sagd produced field case. *Oil Gas Sci. Technol. Rev. d'IFP Energies Nouv.* **2012**, *67*, 289–301. [CrossRef]

37. Abadpour, A.; Bergey, P.; Piasecki, R. 4D seismic history matching with ensemble Kalman filter-assimilation on Hausdorff distance to saturation front. In Proceedings of the SPE Reservoir Simulation Symposium, Woodlands, TX, USA, 18–20 February 2013; OnePetro: Richardson, TX, USA, 2013. [CrossRef]

38. Evensen, G. The Ensemble Kalman Filter: Theoretical formulation and practical implementation. *Ocean Dyn.* **2003**, *53*, 343–367. [CrossRef]

39. Jin, L.; Alpak, F.O.; Hoek, P.J.v.D.; Pirmez, C.; Fehintola, T.; Tendo, F.; Olaniyan, E.E. A comparison of stochastic data-integration algorithms for the joint history matching of production and time-lapse-seismic data. *SPE Reserv. Eval. Eng.* **2012**, *15*, 498–512. [CrossRef]

40. Jin, L.; Weber, D.; Hoek, P.v.D.; Alpak, F.; Pirmez, C. 4D Seismic history matching using information from the flooded zone. *First Break.* **2012**, *30*, 11. [CrossRef]

41. Obidegwu, D.; Chassagne, R.; MacBeth, C. Seismic assisted history matching using binary maps. *J. Nat. Gas Sci. Eng.* **2017**, *42*, 69–84. [CrossRef]

42. Tillier, E.; Da Veiga, S.; Derfoul, R. Appropriate formulation of the objective function for the history matching of seismic attributes. *Comput. Geosci.* **2012**, *51*, 64–73. [CrossRef]

43. Davolio, A.; Schiozer, D.J. Probabilistic seismic history matching using binary images. *J. Geophys. Eng.* **2017**, *15*, 261–274. [CrossRef]

44. Kretz, V.; Vallès, B.; Sonneland, L. Fluid front history matching using 4D seismic and streamline simulation. In Proceedings of the SPE Annual Technical Conference and Exhibition, Houston, TX, USA, 26–29 September 2004; OnePetro: Richardson, TX, USA, 2004. [CrossRef]

45. Leeuwenburgh, O.; Arts, R. Distance parameterization for efficient seismic history matching with the ensemble Kalman Filter. *Comput. Geosci.* **2014**, *18*, 535–548. [CrossRef]

46. Zhang, Y.; Leeuwenburgh, O. Ensemble-based seismic history matching with distance parameterization for complex grids. In Proceedings of the ECMOR XV-15th European Conference on the Mathematics of Oil Recovery, Amsterdam, The Netherlands, 29 August–1 September 2016; European Association of Geoscientists & Engineers: Utrecht, The Netherlands, 2016; p. cp-494. [CrossRef]

47. Zhang, Y.; Leeuwenburgh, O. Image-oriented distance parameterization for ensemble-based seismic history matching. *Comput. Geosci.* **2017**, *21*, 713–731. [CrossRef]

48. Trani, M.; Moncorgé, A.; Bergey, P.; Chen, Y. Fluid Front History Matching Using an Iterative Ensemble Smoother. In Proceedings of the 77th EAGE Conference and Exhibition 2015, Madrid, Spain, 1–4 June 2015; European Association of Geoscientists & Engineers: Utrecht, The Netherlands, 2015; Volume 2015, pp. 1–5. [CrossRef]

49. Sethian, J.A. A fast marching level set method for monotonically advancing fronts. *Proc. Natl. Acad. Sci. USA* **1996**, *93*, 1591–1595. [CrossRef] [PubMed]

50. Hassouna, M.S.; Farag, A.A. MultiStencils fast marching methods: A highly accurate solution to the eikonal equation on cartesian domains. *IEEE Trans. Pattern Anal. Mach. Intell.* **2007**, *29*, 1563–1574. [CrossRef]

51. Emerick, A.A.; Reynolds, A.C. History matching time-lapse seismic data using the ensemble Kalman filter with multiple data assimilations. *Comput. Geosci.* **2012**, *16*, 639–659. [CrossRef]

52. Emerick, A.A.; Reynolds, A.C. Ensemble smoother with multiple data assimilation. *Comput. Geosci.* **2013**, *55*, 3–15. [CrossRef]

53. Pollock, D.W. Semianalytical computation of path lines for finite-difference models. *Groundwater* **1988**, *26*, 743–750. [CrossRef]

54. Sovold, K.; Rian, D.T.; Sandvik, A. Front Tracking Applied to the Simulation of Water Flooding in a Braided River System. In Proceedings of the SPE Latin America Petroleum Engineering Conference, Rio de Janeiro, Brazil, 14–19 October 1990; OnePetro: Richardson, TX, USA, 1990. [CrossRef]

55. Journel, A.G.; Isaaks, E.H. Conditional indicator simulation: Application to a Saskatchewan uranium deposit. *J. Int. Assoc. Math. Geol.* **1984**, *16*, 685–718. [CrossRef]

56. Matheron, G.; Beucher, H.; de Fouquet, C.; Galli, A.; Guerillot, D.; Ravenne, C. Conditional simulation of the geometry of fluvio-deltaic reservoirs. In Proceedings of the SPE Annual Technical Conference and Exhibition, Dallas, TX, USA, 27–30 September 1987; OnePetro: Richardson, TX, USA, 1987. [CrossRef]
57. Avansi, G.D.; Maschio, C.; Schiozer, D.J. Simultaneous history-matching approach by use of reservoir-characterization and reservoir-simulation studies. *SPE Reserv. Eval. Eng.* **2016**, *19*, 694–712. [CrossRef]

58. Skjervheim, J.-A.; Evensen, G.; Hove, J.; Vabø, J.G. An ensemble smoother for assisted history matching. In Proceedings of the SPE Reservoir Simulation Symposium, The Woodlands, TX, USA, 21–23 February 2011; OnePetro: Richardson, TX, USA, 2011. [CrossRef]

59. Emerick, A.A. Analysis of the performance of ensemble-based assimilation of production and seismic data. *J. Pet. Sci. Eng.* **2016**, *139*, 219–239. [CrossRef]

60. Berthet, P.; Trani, M. A Method for Obtaining at Least One Physical Property of a Subsurface Volume of a Hydrocarbon Reservoir Over Time (European Priority Application, Filing Date 4 December 2020, Publication Number EP4009086). Available online: https://www.sumobrain.com/patents/wipo/Method-obtaining-at-least-one/WO2022117735A1.html (accessed on 1 November 2023).

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.