Geophysics-Based Fluid-Facies Predictions Using Ensemble Updating of Binary State Vectors

Fluid flow simulations are commonly used to predict the fluid displacement in subsurface reservoirs; however, model validation is challenging due to the lack of direct measurements. Geophysical data can be used to monitor the displacement of the fluid front. The updating of the fluid front location in two-phase flow problems based on time-lapse geophysical data can be formulated as an inverse problem, specifically a data assimilation problem, where the state is a vector of binary variables representing the fluid-facies and the observations are measurements of continuous geophysical properties, such as electrical or elastic properties. In geosciences, many data assimilation problems are solved using ensemble-based methods relying on the Kalman filter approach. However, for discrete variables, such approaches cannot be applied due to the Gaussian-linear assumption. An innovative approach for mixed discrete-continuous problems based on ensemble updating of binary state vectors is presented for fluid-facies prediction problems with time-lapse geophysical properties. The proposed inversion method is demonstrated in a synthetic two-dimensional simulation example where water is injected into a reservoir and hydrocarbon is produced. Resistivity values obtained from controlled-source electromagnetic data are assumed to be available at different times. According to the results, the proposed inversion method is to a large extent able to reproduce the true underlying binary field of fluid-facies.


Introduction
Monitoring the spatial distribution of fluids in the subsurface is a challenging modeling problem due to the uncertainty in rock and fluid properties. If the petrophysical properties, (e.g., porosity and permeability) and the subsurface conditions (e.g., temperature and pressure) are known, the fluid volumes and their spatial distribution can be predicted by solving partial differential equations for fluid flow in porous media (Bear 2013). However, in practical applications, the majority of the model parameters cannot be measured directly and can only be approximated from indirect geophysical measurements (Aki and Richards 1980;Doyen 2007). Hence, the accuracy in the model predictions depends on the quality of the data and the approximations in the physical models. For these reasons, it is necessary to update the fluid models every time new geophysical data are available.
One of the main challenges in monitoring fluid flow in the subsurface is obtaining accurate and precise predictions of the fluid front at different times. Geophysical surveys, acquired at the surface, provide measurements of physical properties whose changes in time reflect variations in rock and fluid properties. Hence, the location of the fluid front can be updated based on time-lapse geophysical data and fluid flow simulations. Examples of time-lapse geophysical surveys include seismic data that depend on changes in elastic properties and electromagnetic data that depend on changes in electrical properties.
Changes in fluid saturations modify the elastic and electrical response of reservoir rocks. Therefore, unknown fluid flow saturations can be predicted as the solution of an inverse problem where the data are geophysical observations and the governing equations are geophysical models (Doyen et al. 2000). However, the resolution of geophysical data measured at the surface is generally lower than the desired resolution of the fluid distribution model due to the bandlimited nature of the signal and the relatively low signal-to-noise ratio. Therefore, prediction of the fluid saturations is highly uncertain and saturation models can be inaccurate. In many applications, geophysical data can be used to interpret the fluid phase at a specific location; however, it is generally difficult to accurately predict fluid percentages. For this reason, rather than using continuous variables representing fluid volumes or fluid saturations, in this work, the inverse problem is formulated in terms of a discrete variable. For multiphase flow problems, the discrete variable represents fluid-facies (or fluid units), such as water-, air-, hydrocarbon-, or CO 2 -saturated rocks, depending on the application. In particular, for two-phase flow problems, fluid-facies can be represented by a binary variable. The proposed methodology is defined in a probabilistic setting; therefore, the model property of interest is a discrete random variable.
In this work, the focus is on fluid properties that cause changes in the electrical response measured by electromagnetic data, such as electrical resistivity tomography (ERT) and controlled-source electromagnetic (CSEM). The goal of this work is to predict fluid-facies in subsurface reservoirs based on estimated resistivity values from time-lapse geophysical data. From a mathematical point of view, this is a data assimilation problem where the unknown state is a discrete random variable (the fluid-facies) and the observed data are geophysical measurements (the resistivity values). The resistivity of fluid-saturated porous rocks is obtained from electromagnetic measurements.
Here, it is assumed that resistivity values have been computed from CSEM data as a result of a preliminary inversion process. In near surface geophysics applications, ERT data are commonly acquired to investigate the seasonal water displacement (Flinchum et al. 2018;Kotikian et al. 2019;Claes et al. 2020). In hydrocarbon exploration and carbon sequestration, CSEM data are often measured to study the fluid spatial distribution and monitor the fluid front location (Weitemeyer et al. 2006;Lien and Mannseth 2008;Orange et al. 2009;Constable 2010;MacGregor 2012;Lien et al. 2014). The prediction of resistivity from electromagnetic data is itself an inverse problem (Gasperikova and Hoversten 2006;Buland and Kolbjørnsen 2012;MacGregor 2012;Bhuyian et al. 2012). Probabilistic and deterministic approaches have been proposed with different model parameterizations, in terms of resistivity or saturation. Time-lapse electrical resistivity inversions have been proposed in several applications (Berre et al. 2011;Shahin et al. 2012;Lien et al. 2014;Tveit et al. 2015;Commer et al. 2016;Bergmann et al. 2017).
Classification methods for litho-facies and fluid-facies based on geophysical data have also been proposed. Several clustering analysis methods described in Hastie et al. (2009) have been used for facies classification in geophysics inverse problems (Doyen 2007). Classification methods include supervised and unsupervised techniques (Hastie et al. 2009;Martinez and Martinez 2015). Clustering and pattern recognition methods have been used to classify geophysical measurements. However, the majority of these applications focuses on static characterization problems, with a spatial correlation component (i.e. facies are spatially correlated to mimic the geological continuity) but without a temporal component (i.e. facies are predicted at a given time step, typically before dynamic processes start).
The focus of this work is time-dependent fluid-facies characterization problems where the spatial distribution of fluid-facies changes in time and is monitored using time-lapse geophysical properties. Therefore, in this study a dynamic fluid-facies classification is presented, and it is applicable to geological dynamic problems where one fluid (e.g., water) replaces another fluid (e.g., hydrocarbon) in rock formations. Reservoir modeling with time-lapse data is a data assimilation problem where the model variables are predicted and updated when new measurements become available. Several stochastic optimization methods have been proposed for data assimilation problems, and during the last decade, ensemble-based methods have become the most popular stochastic data assimilation method in geoscience applications.
Data assimilation can refer to a range of different inference procedures, of which the two most common are filtering and smoothing. In the present article, the focus is exclusively on the filtering problem. There are two main classes of ensemble-based filtering methods: particle filters (Doucet et al. 2001) and ensemble Kalman filters (EnKFs) (Evensen 2009). Hybrid versions of these filters have also been proposed. Particle filters have the advantage of being exact in the sense that as the ensemble size goes to infinity, the ensemble representations of the series of filtering distributions converge to the corresponding correct series of distributions. Particle filters are also very general as they do not rely on any specific assumptions about the distributions of the model variables. Hence, particle filter methods are, in principle, applicable to both discrete and continuous variables. In practice, however, particle filters are known to collapse when the dimension of the state vector is large. The EnKF is a filtering method which relies on a linear-Gaussian assumption about the underlying model. Despite the linear-Gaussian assumption, studies show that the EnKF provides good results even in non-linear, non-Gaussian problems, and unlike the particle filter it also scales well to problems with very high-dimensional state vectors. The EnKF has been applied to geophysical data assimilation and history matching problems using seismic and electromagnetic data (e.g. Tveit et al. 2015;Tveit et al. 2020). Recent publications focus on the integration of fluid flow simulation and geophysical data assimilation for the monitoring of the fluid front location (Trani et al. 2012;Leeuwenburgh and Arts 2014;Zhang and Leeuwenburgh 2017). However, the Gaussian approximations make the EnKF applicable only in situations with continuous variables. For problems with discrete variables, such as fluid-facies classification, the filter is not appropriate.
Ensemble filtering of discrete variables is a challenging problem which has received fairly little attention in the literature compared to filtering of continuous variables. Oliver et al. (2008) propose a strategy where the EnKF is used to update the discrete variables. Specifically, they propose a two-step strategy where, in the first step, the EnKF is used to update the discrete variables as if they were continuous, and in the second step, the updated continuous-valued variables are mapped back to the original discrete state space using the Viterbi algorithm (Viterbi 1967). Loe and Tjelmeland (2020) present an alternative updating method for binary vectors in one-dimensional space based on a generalized approach of the EnKF. Instead of using a linear-Gaussian model assumption in the ensemble update, as in the EnKF, they construct an update based on a hidden Markov model assumption. To capture as much information as possible from the forecast ensemble, including potential non-Markov properties, the expected number of components of the binary state vector that remain unchanged is maximized.
This paper presents an ensemble-based data assimilation method for a problem where the state vector at each time step is a vector of binary variables and the observations are continuous-valued estimated resistivity values. The binary variables of the state vector represent two different fluid-facies, for example water and hydrocarbon or CO 2 , and each binary variable is connected to a continuous-valued variable representing water saturation. High water saturation values indicate the presence of the water facies, while low saturation values indicate the presence of the other fluid-facies. The proposed ensemble filtering method alternates between a forecast step performed in the continuous state space of the saturation variable and an update step performed in the discrete state space of the fluid-facies variable, and between each step an appropriate mapping from one state space to the other is performed. The update step is performed according to the updating procedure for binary state vectors proposed in Loe and Tjelmeland (2020). The proposed inversion method is demonstrated in a synthetic two-dimensional example representing a two-phase flow problem with resistivity values available at different times. According to the results, the proposed procedure is to a large extent able to reproduce the true underlying binary field of fluid-facies. Larger ensemble sizes provide more accurate results, but the results obtained with smaller ensemble sizes are also satisfactory.
The remains of this paper takes the following outline. First, Sect. 2 formulates the inverse problem more formally and presents the proposed ensemble-based inversion method. Next, Sect. 3 presents numerical results based on a two-dimensional synthetic model for a two-phase fluid flow problem. Finally, a few closing remarks are given in Sect. 4.

Inverse Problem Setting
The problem addressed in this work is the prediction of fluid-facies from time-lapse resistivity values. Consider a time series vector of resistivity measurements recorded at time t i , and T ⊆ {1, . . . , N t }, the goal is, for each time step i = 1, . . . , N t , to assess the distribution of fluid-facies k i in the reservoir. Notice from the set T that an observation d i may be available at every time step i = 1, . . . , N t , or just a subset of them.
In this work, each component k j i of k i is assumed to be connected to a continuous variable m j i ∈ [s wi , 1] representing water saturation, where s wi is an irreducible water saturation value; that is, the fraction of water that a porous rock can retain due to non-connected porosity, low permeability and/or capillary forces. Here, the irreducible water saturation value s wi = 0.2 is assumed. Given {m i } N t i=1 , the resistivity data {d i } i∈T are assumed to be conditionally independent, so that the vector d i at time t i depends only on m i according to where f is a known, possibly non-linear function, and the variable e i ∈ R N e is an N e -dimensional vector of measurement random errors assumed to follow a known probability distribution. Similarly, the saturation m i+1 at time step i + 1, given all the saturation values m 1 , . . . m i up to time step i, depends only on m i according to a known forward model, for i = 1, . . . , N t , where g is the fluid flow simulation, generally given by a system of partial differential equations solved by finite difference methods (Aziz 1979). The goal of this work is, for each time step i = 1, . . . , N t , to assess the filtering distribution p(k i |d 1:i ), where d 1:i = d j ; j ∈ T ∩ { j ≤ i} , that is, the distribution of fluid-facies k i given all the resistivity data up to time t i . Only K = 2 fluid-facies are assumed in this work: facies 1 represents water and facies 0 represents another fluid-facies. The relationship between the fluid-facies k j i ∈ {0, 1} and the saturation value m j i ∈ [s wi , 1] is assumed as where r ∈ (s wi , 1) is some appropriate threshold. The value of the parameter r might vary from one application to another. A reasonable choice is to set r = 0.5 such that each fluid facies is named after the predominant fluid component. However, from a reservoir management perspective, the focus is generally on areas with a high concentration of hydrocarbon. Therefore one could choose a lower value to identify the regions that are economically valuable.

Forward Model
The prediction of the time-dependent electrical response of a reservoir model requires a rock-physics model to link the petrophysical properties, such as porosity and fluid saturations, to the resistivity of the saturated porous rocks and a fluid flow simulation model to compute the saturation at a given time step, given the saturation at the previous time step. In the proposed approach, porosity and permeability are assumed to be estimated from pre-injection geophysical measurements (e.g., seismic data). Alternatively, multiple geostatistical simulations of porosity and permeability can be generated to repeatedly apply the methodology to an ensemble of realizations; however, the computational cost would linearly increase with the number of realizations.
A rock-physics model is a relationship to predict the geophysical response of saturated porous rocks. Assuming that the porosity φ of the porous rock is known, the resistivity R (the measured data d in the inverse problem) of the porous rock saturated with water saturation s w can be predicted using Archie's law (Mavko et al. 2009 where R w is the resistivity of formation water, a is the cementation exponent, and b is the saturation exponent (Mavko et al. 2009). The parameters R w , a and b in Eq.
(4) are assumed to be constant in time. Archie's equation is valid for clean sandstone formations. For formations with a small to medium clay volume, Archie's equation can be modified to account for the conductivity of the clay mineral as in Simandoux and Poupon-Leveaux models (Mavko et al. 2009). The dynamic model that governs two-phase fluid flow in porous media is based on the constitutive equations of mass and momentum balance. The model is numerically solved using the black-oil framework to predict the saturation and pressure at each time step, given the initial rock and fluid parameters (Aziz 1979). In this work, the MATLAB Reservoir Simulation Toolbox (Lie 2019) is adopted.

Inversion Method
To solve the inverse problem presented above, an ensemble-based strategy where the forecast step is performed in the continuous domain of m i and the update step is performed in the discrete domain of k i is adopted. At each time t i , an ensemble of fluid-facies fields k represents the distribution of k i given the resistivity data up to time t i−1 , that is d 1:i−1 . Likewise an ensemble of saturation fields represents the distribution of m i given the same resistivity data. Correspondingly, the distributions of k i and m i given resistivity data up to time t i , that is d 1:i , are also represented by ensembles, which are denoted by , respectively. The main steps of the inversion procedure are summarized in Algorithm 1, while each step is studied in closer detail in the following sections.

The Update Step
As summarized in Algorithm 1, the update step of the proposed approach involves two parts. First, the ensemble using the assumed relation between k i and m i in Eq. (3). Second, is updated to take the new observation d i at time t i into account. In the following, the two parts of the update step are discussed in more detail.
The ensemble of saturation fields m by simply applying Eq.
(3) to each element in each of the ensemble members, that is, set for each location j = 1, . . . , N k for each ensemble member l = 1, . . . , M.
To update the ensemble k to take the new observation d i into account, the procedure proposed in Loe and Tjelmeland (2020) is adapted to the situation considered in the present article. In the present article, it is assumed that the fluid-facies k i at each time step i = 1, . . . , N t is defined on a two-dimensional lattice. However, the method in Loe and Tjelmeland (2020) is applicable only for vectors with a one-dimensional spatial arrangement. Therefore, in order to apply their procedure, the updating of each column in the lattice is done independently of the others. Of course, this is not an ideal approach since it means that some of the spatial correlation in the horizontal direction is lost; however, since the forecast step incorporates spatial correlation in both directions, one may still obtain satisfactory results. Let C denote the number of columns in the lattice and let k The distribution of the latent u j i should depend on the fluid-facies value k j i , and it is assumed that and where c 0 and c 1 are normalizing constants, and λ 0 and λ 1 are parameters specifying the level of noise in the resistivity measurements. Small values of λ 0 and λ 1 reflect noisy resistivity data, while higher of λ 0 and λ 1 reflect less noisy resistivity data. Essentially, the auxiliary variable u j i can be interpreted as a noisy realisation of the saturation value m j i . The logic behind the choice of distributions in Eqs. (8) and (9) is that it should be more likely to generate u The right plot in Fig. 1  Combining the estimated prior Markov chain for column c with the likelihood model specified above, the corresponding posterior distribution also becomes a non-stationary Markov chain. The properties of this posterior Markov chain are computationally easy to compute, and in particular the bivariate distributions for every two neighbor nodes in column c can be found. To update the prior ensemble members of column c, a distribution q k (l) i,c k (l) i,c which preserves these bivariate distributions is constructed. More specifically, under the assumption that the estimated prior Markov chain for column c is correct, k (l) i,c is updated by simulating from a conditional distribution q k (l) i,c k (l) i,c such that the bivariate distribution for every pair of neighbor nodes iñ k (l) i,c is equal to the corresponding bivariate distribution of the posterior Markov chain for column c.
The chosen distribution q k (l) i,c k (l) i,c for updating the prior ensemble members of column c can be expressed as where S is the number of rows in the lattice used to represent k i , and k i,c is a Markov chain with initial distribution specified by q 1 (·|·) and transition probabilities specified by q s (·|·, ·), s = 2, . . . , S. To specify the updating procedure completely it now remains to specify q 1 (·|·) and q s (·|·, ·), s = 2, . . . , S. These are specified to accomplish two goals. First, considering k (l) i,c as a sample from the estimated prior Markov chain, the marginal bivariate distributions for every pair k (l) i,(s−1,c) ,k (l) i,(s,c) , s = 2, . . . , S should be identical to the corresponding bivariate distribution in the posterior Markov chain discussed above. This requirement ensures that the updated fluid-facies values reflect the new resistivity data d i . However, with only this requirement many possible solutions exist for q 1 (·|·) and q s (·|·, ·), s = 2, . . . , S, so there is room for formulating another goal. Still considering k (l) i,c as a sample from the estimated prior Markov chain, the second goal for the updating of k (l) i,c is to maximise the expected number of elements in k (l) i,c that remain unchanged; that is, the goal is to maximize where I (A ) equals one if the event A is true, and zero otherwise, and the expectation is taken with respect to the joint distribution of k i,c . This requirement makes the updating robust with respect to the a priori Markov chain assumption made for k (l) i,c , l = 1, . . . , M. If the true distribution of k (l) i,c , l = 1, . . . , M is not a Markov chain, many of its non-Markov properties will prevail intok (l) i,c , l = 1, . . . , M since it is specified that as few changes as possible should be made to k (i) i,c in the generation ofk (l) i,c . Numerically, it turns out that in that the maximization of the expression in Eq. (13) under the constraints for the specified bivariate distributions for the pairs k (l) i, (s−1,c) ,k (l) i, (s,c) , s = 2, . . . , S can be efficiently computed using a combination of dynamic programming and linear programming. The details of the optimization algorithm are discussed in Loe and Tjelmeland (2020).

The Forecast Step
The forecast step of the proposed approach also involves two parts. First, the ensem- Second, the forward model in Eq.
(2) is used to generate m . . . ,m (M) i . In the following, the two parts of the forecast step are discussed in more detail.
To generate the saturation fieldm (l) i based on a given fluid-facies fieldk (l) i , the fluid-facies indicators ink (l) i are first used to define a lattice of distances, δ, from each node j to a node with the opposite value of node j. More precisely, the values in δ are  Fig. 2 To the left: True porosity and log permeability models, with the locations of the production and injection wells marked P and I, respectively. To the right: Assumed porosity and log permeability models defined sequentially as follows. First, δ j = 0 is set for all nodes j that has one or more neighbor node j so thatk . Thereafter, δ j = 1 is set for all nodes j for which δ j is still undefined and which has a neighbor node j with δ j = 0. Thereafter, δ j = 2 is set for all nodes j for which δ j is still undefined and which has a neighbor node j with δ j = 1. This process is continued until δ j is defined for all nodes j. The next step is to scale the δ j values into the [0, 1] interval. Letting denote the scaled field, the value for node j is defined as where δ max > 0 is a parameter controlling the size of the transition zone from s wi to 1. The larger the value of δ max , the larger the size of the transition zone. One should choose a value for δ max based on what one believes is a realistic transition for the application in consideration. In the numerical examples in Sect. 3, δ max = 8 is used. The field defines a trend for them (l) i values. To add noise to this trend a slightly modified version of the so-called smootherstep function is first used to transform the where Φ(·) is the cumulative distribution function of a standard normal distribution and Φ −1 (·) is the inverse of this function. The effect of Eq. (15) is that the value ν j is in the left tail of a normal distribution with mean Φ −1 (r ) and unit variance when j = 0, or in the right tail of the same distribution when j = 1, and with a smooth transition between these two extremes. Moreover, the last term in Eq. (15) ensures that ν j = Φ −1 (r ) when j = 0.5. A noisy version z of ν is then defined by setting Case 1 24 (every 6 months) Low λ 0 = 9.8, λ 1 = 5.0 Case 2 4 (every 3 years) Low λ 0 = 9.8, λ 1 = 5.0 Case 3 4 (every 3 years) High λ 0 = 7.8, λ 1 = 2.5 where ε is a Gaussian field with zero mean, unit variance and an exponential correlation function, and α > 0 is a parameter controlling the noise level. Finally, the saturation field is defined by transforming the z field back to the (s wi , 1) interval, The second part in the forecast step, to generate the ensemble m

Application
The proposed inversion method is tested using a synthetic reservoir model. The model consists of a two-dimensional reservoir, 25 m × 25 m, with constant thickness and four main channels with high-porosity rocks surrounded by low-porosity rocks; see Fig. 2. The fluid system includes two fluid phases: oil and water. Therefore, two fluid-facies are defined: oil-saturated rocks (corresponding to the value 0) and water-saturated rocks (corresponding to the value 1). The discretized reservoir is defined on a 128×128 grid, and the well configuration includes four injectors and six producers as shown in Fig. 2. The oil production mechanism is based on water injection simulated using the MATLAB Reservoir Simulation Toolbox (Lie 2019) for a time period of 12 years. The 12 year time period is discretized into 24 equidistant time points t 1 , . . . , t 24 such that each step of the simulation involves propagating the system six months forward in time. During the simulation, injection rates are kept constant at the injector locations, and bottom hole pressure is kept constant at the producer locations. Initially, the entire reservoir is filled with hydrocarbon with the irreducible water saturation value of  Figure 3a, b show the fluid-facies k i and saturation values m i , respectively, in the reservoir at the time steps i = 6, i = 12, i = 18 and i = 24; that is, after t 6 = 3, t 12 = 6, t 18 = 9 and t 24 = 12 years of the simulation. Figure 3c shows corresponding reference resistivity values (in log-scale); that is, the resistivity values one obtains by inserting the true water saturation values into Archie's law in Eq. (4). Pretending that the fluid-facies and saturation values used to generate the plots in Fig. 3 are unknown, the goal of the simulation experiment is to estimate the fluid-facies field at each time step based on noisy resistivity data. In this example, the resistivity data d i at time t i includes a two-dimensional map of resistivity measurements; specifically, the dimensionality N d of d i is equal to the dimensionality N k of k i (and m i ) so that an observation d j i is available for each variable k j i of k i . Since a 128 × 128 grid is considered, with a fluid-facies variable k j i in every cell j, the dimensions N k and N d are N k = N d = 128 · 128 = 16384.
The porosity and permeability models shown in Fig. 2 (left plots) are the true porosity and permeability models of the reservoir. These were the values used to generate the reference model shown in Fig. 3. Since porosity and permeability are generally not known, a reservoir model of assumed porosity and permeability models is built to mimic the resolution of a reservoir model estimated from pre-production seismic data. The assumed porosity and permeability models are shown in Fig. 2 (right plots).
Three case studies are presented, differing in the frequency with which the resistivity measurements are collected and the amount of noise in the measurements; see Table 1. The first case, referred to as case 1, represents an idealized situation where resistivity measurements are recorded frequently and the degree of noise in the data is small. Specifically, observations are assumed to be recorded every six months, or at every time step i = 1, . . . , N t . Hence the set T introduced in Sect. 2.1 is T = {1, 2, . . . , 24}. Figure 4a shows the simulated resistivity measurements d i (in log-scale) at the four time steps i = 6, 12, 18 and 24 for case 1. The resistivity data were generated with the likelihood model specified in Sect. 2, using the true fluid-facies shown in Fig. 3a and the assumed porosity model shown to the right in Fig. 2, and with the parameters λ 0 and λ 1 set to λ 0 = 9.8 and λ 1 = 5. These values for λ 0 and λ 1 represent optimistic noise conditions. In the second case, referred to as case 2, the same data as in case 1 are considered, but observations are assumed to be acquired only every three years of the simulation period; that is, an observation is recorded after 3, 6, 9, and 12 years, or at the time steps i = 6, 12, 18 and 24. Hence the set T is in this case T = {6, 12, 18, 24}, and the likelihood parameters λ 0 and λ 1 are the same as in case 1. In the third case, referred to as case 3, observations, as in case 2, are acquired only every 3 years, but a different set of data with a much higher level of noise is considered. Hence this case represents the most realistic of the three cases. Figure 4b shows the simulated resistivity measurements for case 3. Similarly to the resistivity data for cases 1 and 2, the resistivity measurements for case 3 were generated using the likelihood model specified in Sect. 2, but with the parameters λ 0 and λ 1 set to λ 0 = 7.8 and λ 1 = 2.5. These parameter values represent realistic noise conditions. For all three case studies the proposed inversion method is tested using three different ensemble sizes: M = 20, M = 100 and M = 500. The parameters δ max and α in Eqs. (14) and (16) are set to δ max = 8 and α = 0.2. The ensembles were initialised by first introducing an initial field of fluid facies k 0 for which it is assumed that k j 0 = 0 for every cell j in the reservoir, and thereafter generate each m = 1|d 1:i ) obtained in the ten runs of case B when using ensemble size M, the standard deviation of these ten estimates is computed. Results are shown in Figs. 8, 9 and 10 for cases 1, 2 and 3, respectively. Similarly to the other results presented above, the results obtained with the higher ensemble sizes M = 100 and M = 500 are overall smoother and less noisy than those obtained with the rather small ensemble size M = 20. A general trend, however, for all three ensemble sizes and all three cases, is that the standard deviations tend to be higher near the boundary of the fluid front, which is reasonable, since this is the most uncertain area where changes occur. Moreover, the results from case 3 in Fig. 10 are considerably noisier than the results from cases 2 and 3 in Figs. 8 and 9. This means that the case 3 results tend to vary more from one run to another. Again, this is reasonable, since the resistivity measurements in case 3 are more uncertain.

Conclusions
A novel method for monitoring and updating the evolution of fluid-facies from timelapse geophysical properties in a two-phase flow problem has been presented. The inversion method is based on an ensemble filtering method where the updating of the prior ensemble at each time step is performed using a particular updating method for binary vectors. The main novelty of the work is the extension of ensemble-based methods to mixed discrete-continuous problems to update the spatial distribution of fluid-facies. In the proposed application, the geophysical dataset includes time-lapse resistivity values that are assumed to have been estimated from CSEM data through a preliminary inversion process. The proposed method is tested in a synthetic example with a two-dimensional reservoir model. The results from this synthetic example are accurate and support the validation of the proposed methodology. In real data applications, the accuracy of the results depends on the quality of the data in terms of resolution and signal-to-noise ratio, and also on the accuracy of the fluid flow simulator. The main limitation of this work is that uncertainty in the estimation of porosity and permeability are not taken into account. Future research directions aim to extend the proposed method so that porosity and permeability are also treated as random variables and so that the geophysical dataset includes measured data such as electromagnetic amplitude and phase.
Funding Open Access funding provided by NTNU Norwegian University of Science and Technology (incl St. Olavs Hospital -Trondheim University Hospital).
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.