Simultaneous Inversion of Production Data and Seismic Attributes: Application to a Synthetic SAGD Produced Field Case.

Résumé — Inversion simultanée des données de production et des attributs sismiques : application à un champ synthétique produit par injection de vapeur — L’utilisation conjointe des données de production et des attributs de sismique répétée facilite la compréhension des mouvements de fluide dans les formations géologiques et aide à la construction de modèles numériques fiables représentant ces formations. Ceci a récemment motivé le développement de techniques d’inversion ou de calage dédiées à l’identification de modèles cohérents avec l’ensemble des données disponibles. La méthodologie décrite dans ce papier permet de déterminer les cartes de propriétés pétrophysiques comme les faciès, les porosités ou les perméabilités à partir des données de production et des attributs de sismique répétée. Elle est appliquée avec succès à un champ synthétique d’huile lourde produite par injection de vapeur. Ce cas d’étude, inspiré d’un cas réel, illustre comment définir une méthodologie d’inversion qui respecte l’ensemble des données disponibles. L’étude est centrée sur trois points particuliers de la méthodologie qui sont les clés du succès. Premièrement, le choix de la paramétrisation est essentiel. La paramétrisation doit permettre des variations locales avec peu de paramètres. Deuxièmement, cette étude met Abstract — Simultaneous Inversion of Production Data and Seismic Attributes: Application to a Synthetic SAGD Produced Field Case — The joint use of production data and time-lapse seismic attributes can help to understand fluid flows within geological formations and to build reliable numerical models for representing these formations. This concern recently motivated the development of dedicated inversion or matching techniques for identifying models consistent with all collected data. The methodology presented in this paper makes it possible to map petrophysical properties such as facies, porosity and permeability into reservoirs from production data and seismic attributes. It is successfully applied to a synthetic case describing a heavy oil field produced from steam assisted gravity drainage. This case study generated from a real case shows how to design the inversion methodology to match the entire set of available data. A few key points are highlighted since they drive the success of the proposed matching methodology. First, parameterization is essential. It must allow for locally varying petrophysical properties from a reduced number of parameters. Second, this study stresses the need for alternative formulations for quantifying the mismatch between reference seismic attributes and simulated seismic attributes. Last, two methods are compared for integrating reference seismic attributes either in time or in depth (after a time-to-depth conversion). In the case studied, it is shown that the two approaches are equivalent since time-to depth-conversion error is quite small.


INTRODUCTION
Fluid injection or production in underground geological formations induces changes in saturation, pressure and/or temperature, which can be detected from seismic data. This has been stirring the development of time-lapse (4D) seismic to monitor fluid displacements (Eastwood et al., 1994;Lumley, 2001;Arts et al., 2004;Roggero et al., 2007;Byerley et al., 2009). The recent advances in the quantitative processing of the so-obtained data favor the fusion of geophysics and reservoir engineering to build numerical models representative of reservoirs. A reservoir model is a grid populated by petrophysical properties such as facies, porosities, permeabilities, etc. When provided as an input to a flow simulator, it helps to understand how fluid flows. This modeling step is essential in several disciplines of geosciences: for instance, for optimizing oil production in reservoir engineering or for predicting the migration of pollutants in hydrology. A reservoir model is all the more relevant as it respects the data collected on the field: production data (bottom hole pressures, flow rates, tracer concentrations measured at wells) as well as 4D-seismic attributes (variations in impedances, velocities, time thicknesses). Identifying such a model is an inverse problem. It entails the iterative minimization of an objective function, that quantifies the data mismatch. Several techniques, referring to history-matching, were developed and shown to be efficient when dealing with production data (Carrera et al., 2005;Oliver and Chen, 2010). Recent works have been focusing on the joint integration of production data and 4D-seismic attributes into reservoir models (Vasco et al., 2004;Stephen et al., 2006;Roggero et al., 2007). However, formulating an efficient and tractable methodology is still a challenging issue.
In this paper, we introduce an inversion method to constrain reservoir models to both production data and 4D-seismic attributes. Then, we build a synthetic, but realistic case defined from a real bitumen field (Lerat et al., 2010) produced from Steam Assisted Gravity Drainage (SAGD). It is used in the last section to illustrate the applicability of the proposed joint inversion approach. The following questions have arisen. In what way should seismic attributes be used in the matching process? How should we quantify the mismatch between reference seismic attributes and the numerical seismic responses? Which parameterization is suitable for capturing geological heterogeneities?

INVERSE PROBLEM
The inverse problem is solved by comparing the actual data (production data and seismic attributes) to the corresponding numerical responses. It hinges on the minimization of an objective function usually defined in a least-square sense: (1) m is the set of uncertain parameters describing the reservoir model and d obs includes the data to be matched. Covariance matrix C D quantifies the experimental and theoretical uncertainties. It is often assumed to be diagonal although this assumption may be questionable especially for seismic attributes. g is the operator mapping the parameter space to the data space.

Forward Modeling
The g operator introduced above includes successive modeling steps, which yields the forward model. It basically involves geological modeling, upscaling, flow modeling and petroelastic modeling. These steps are briefly recapped hereafter.
Geological modeling consists of populating the grid with petrophysical properties using geostatistical simulation algorithms (Journel et al., 1998;Chilès and Delfiner, 1999). Facies proportions and the second order statistics of porosity and permeability per facies are derived from the analysis of well logs and the baseline seismic attributes. Then, the resulting geological model is upscaled to reduce the number of grid blocks and make flow simulation feasible in a reasonable amount of time (still of the order of a few hours for the case described in Sect. 2). The upscaling process involves the computation of equivalent petrophysical properties to be assigned to coarse grid blocks given the properties of the fine grid blocks. Then, we perform flow simulation. We use a thermal flow simulator, which yields the evolution with time of pressures and rates at wells in addition to the fluid pressure, saturation and temperature values over the whole reservoir grid at given times. The following step is petro-elastic modeling. It aims at determining the seismic attributes at time zero and their variations with time because of pressure, saturation and temperature changes. Elastic properties are predicted referring to Hertz-Mindlin relations (Mindlin, 1949) to account for pressure and to the generalized Gassmann equations (Ciz and Shapiro, 2007) for fluid substitution.
In a few words, Hertz-Mindlin formulas quantify the influence of effective stress changes on the drained bulk and shear moduli, denoted K M and G M : (2) σ eff (t) is the effective mean stress at time t and σ ini eff the initial effective mean stress. h K and h G are the Hertz exponents. Further, Gassmann equations are used to describe saturation effects on the undrained bulk and shear moduli, denoted by K(t) and G(t). The generalized formulation (Ciz and Shapiro, 2007) was preferred for dealing with the heavy oil application case presented hereafter as it accounts for the non negligible shear modulus of heavy oil: (3) ϕ is porosity. K Gr and G Gr are the bulk and shear moduli of the solid matrix whereas K F and G F are the bulk and shear moduli of the phase filling the pore space. Then, impedances are derived from elastic moduli and are transferred into seismic domain. We ensure the same frequency range as seismic signal using an appropriate band pass filter (Mezghani et al., 2004).

Computation of the Objective Function
Overall, the forward modeling workflow provides numerical production and seismic responses to be compared to the actual production data and seismic attributes within the objective function. Finding a model compatible with the observed data calls for the minimization of the objective function. In the case studied, the objective function is the sum of two main terms: one for the production data mismatch, denoted by J prod , and the other one for the seismic attributes mismatch, denoted by J seis . In practice, these terms are written as follows: where nprod is the number of production data series to be matched and ntime j the number of measurement times for data series j. Similarly, nseis is the number of seismic attribute series to be matched and nmesh j the number of measurement for data series j. Coefficients w are weights assigned to data series and σ is the standard deviation of data errors. Normalizing the misfit term by the number of data is a standard technique to balance the terms included in the objective function. Selecting appropriate values for weighting coefficients w is clearly case-dependent and driven by reservoir engineer experience. We selected a gradient-based optimizer (Tarantola, 1987) to drive the minimization process. More precisely, the opti- ϕ mization method used in this work is a Sequential Quadratic Programming and Augmented Lagrangian (SQPAL) solver (Sinoquet and Delbos, 2008). This is a local optimization algorithm where the gradients of the objective function with respect to parameters are estimated by finite differences. The model parameters are sequentially varied until the objective function is small enough or the number of iterations exceeds a maximum value. When J = J prod + J seis (Eq. 4) is used as objective function, the problem is usually ill-posed. A number of methods have been proposed to regularize similar inverse problems (Carrera et al., 2005). A possibility consists of narrowing the search space, which can be performed by selecting an appropriate parameterization technique (Oliver, 2010;Le Ravalec-Dupin, 2010). The choice of appropriate parameterization is decisive: it impacts the final matched model and the overall performance of the optimization algorithm (i.e., it helps to avoid local minima and can make it easier to decrease the objective function).
The subsequent sections discuss two problematical phases of the integration of seismic attributes into reservoir models.

Time-to-Depth Conversion
Flow simulation ends up with pressures, saturations and temperatures computed over a grid whose vertical axis is depth. Since these variables are given as inputs to the petro-elastic model, the resulting calculated seismic attributes are also in depth domain. However, the seismic attributes to be matched are usually known in time domain. Since the reference and the simulated seismic attributes are given in two different domains, the comparison between them is uncomfortable. Two alternatives can be envisioned. The first approach encompasses a time-to-depth conversion applied to the reference seismic attributes and a comparison in depth domain (Dong and Oliver, 2005;Stephen et al., 2006;Roggero et al., 2007). With this approach, the conversion is made once and for all before the history matching process. As a consequence, the velocity model used for conversion does not account for changes applied to parameters during optimization. A second possibility involves a depth-to-time conversion applied to the computed seismic attributes and a comparison in time domain (Fornel et al., 2007;Tillier et al., 2011). The depthto-time conversion is based upon the velocity model derived from the petro-elastic model at each iteration. To sum up, the petro-elastic model yields P-and S-wave velocities in every grid blocks. These values can be used to transfer the simulated attributes from the depth to time domain so that we end up with attributes in time domain at geological scale. This second approach accounts for the optimizer-driven changes previously applied to parameters, among which facies and porosities. Thus, it takes advantage of the persistent consistency of the depth-to-time conversion process with petrophysical properties. One of the issues addressed hereafter is about the suitability of an approach compared to the other one.
Shall we compare reference and simulated seismic attributes in time or depth domains?

Quantifying the Seismic Error Term
The parameters estimated from the inversion procedure are required to satisfy a weighted measure of both production data and seismic attributes (Eq. 4). The efficiency of the joint inversion of production data and seismic attributes strongly depends on the ability to properly capture the differences between the reference data and the corresponding simulated responses. The least-square formulation is known to be appropriate for production data, but production data and seismic attributes are very different by nature. First of all, the support of seismic attributes is evident across the whole grid while production data are collected at wells. Thus, we may have millions of seismic attributes to be balanced with hundreds of production data. In addition, seismic attributes are not raw data: they are derived from a preliminary inversion process with numerous uncertainties and no unique solution. The first inversion attempts, which will be recapped in the subsequent sections, stressed that the least-square formulation may not be relevant for seismic attributes. Keeping in mind that over-parameterization and over-fitting must be avoided when solving inverse problems, we suggest a different approach to measure the seismic mismatch. Instead of exactly matching seismic attributes in every grid blocks, we try to capture their main features. Thus, to measure the differences between two grids populated with seismic attributes (the grid with the reference seismic attributes and the grid with the simulated seismic attributes), we apply a new formulation identified hereafter as the Local Dissimilarity Map (LDM) formulation (Da Le Ravalec et al., 2011). We proceed following three steps. First, we apply filtering and clustering to the two reference and simulated seismic grids to detect the significant seismic features. The filter can be a simple moving average. Clustering involves the classification of seismic attributes into a reduced number of classes. This can be performed for instance with the wellknown K-means algorithm (MacQueen, 1967). We actually split attribute values into two classes so that the seismic attribute grid is converted into a binary black and white grid. The detected seismic features are assumed to be black. Second, we refer to the LDM approach introduced by Baudrier et al. (2008) to quantify the local dissimilarities between the two binary images. It relies on the computation of a local modified Hausdorff distance. The result is a grid with a dissimilarity value assigned to each grid block. For illustration purposes, let us consider a given grid block. Its local dissimilarity distance is zero when its two values in the reference and simulated binary grids are identical. Otherwise, it is set to the (non zero) distance to the nearest black grid block. The last step consists of calculating a global dissimilarity as the sum of the squared local dissimilarities over the whole grid.
The so-obtained value is incorporated into the objective function as the seismic mismatch term. A similar method was recently introduced by Jin et al. (2011) to account for 4D seismic attributes. Just as we did, they converted the seismic images into binary images. However, they only considered a pixel to pixel error while we refer to a sophisticated way to compute the distance between two binary images. Whatever the metric used to quantify the seismic mismatch, there is still a need to choose appropriate weights to balance the contribution of production and seismic data into the objective function. This choice depends on the case studied and is usually based upon engineer experience. The one made in this paper will be detailed in the following section.

Facies Parameterization
Parameterization is a key of success for history matching. In the test case described in this paper, we focus on facies proportions. One of our prerequisite was to apply a parameterization method, which makes it possible to globally or locally vary facies proportions. Keeping this goal in mind, we selected the ratio method developed by Ponsot-Jacquin et al. (2009). When using the truncated Gaussian method to generate facies realizations, we actually thresholds for truncation on the basis of the facies proportion matrix (Doligez et al., 2007). For history matching purposes, the question is then: how to adjust this proportion matrix with a small number of parameters which are modified by the optimizer? The idea suggested by Ponsot-Jacquin et al. (2009) is to define two groups of facies and to consider the ratio of proportions between these two groups. The ratio is then the parameter modified by the optimizer. When the optimizer changes the value of this ratio, all the proportions of the facies selected in groups are modified with respect to this value. Proportions are modified proportionally to their initial value as illustrated in Figure 1. For this simple stationary case, the ratio of interest is the ratio of facies "sand" proportion over the sum of facies "sand", "shale 1" and "shale 2" proportions. The initial value R 0 of this ratio is changed to a new value R 1 , and facies proportions are updated to respect this ratio. A straightforward generalization of this method to the non stationary case is given in Ponsot-Jacquin et al. (2009). This method is efficient since it makes it possible to adjust the whole facies map with a reduced number of parameters. In addition, it can be applied to a subdomain only. However, its main drawback is that proportions are locally modified without considering correlation lengths. Thus, it generates discontinuities with no geological meaning at the boundaries of the modified sub-domain. Then, Tillier et al. (2010) proposed to compute the new ratio value referring to co-kriging, thus preserving geological continuity. However, this involves a significant increase in computation time. This is why we prefer to use the original ratio method in the case described hereafter.

REFERENCE CASE AND DATA
An application case is presented in the following sections. Although synthetic, it is realistic as strongly inspired by a real Canadian field produced from SAGD. The production technique (Butler, 1991) entails two parallel horizontal wells at the bottom of a thick shaly sand reservoir: a production well and a steam injection well located 5 to 10 m above. Briefly, the injected steam forms a growing vapor phase chamber and delivers latent heat to the heavy oil, hence reducing its viscosity. Then, the heated oil drains under gravity within and along the edges of the steam chamber towards the production well.
The reservoir is about 260 m deep. The case studied focuses on a single well pair. Thus, the target sub-domain is 820 m by 98 m wide and about 50 m thick. It is discretized over a 100 × 41 × 98 grid. Near the well, one grid cell is approximately 0.5 m × 20 m × 1 m. The 98 layers were split into 3 units: the top, middle and bottom units with 30, 15 and 53 layers, respectively. The bottom unit consists of coarsegrained sands deposited by high-energy bradied channels. This is the main reservoir unit. The middle unit consists of medium to fine grained sandstones alternated with shale. It is interpreted as meandering channel facies association. The top unit consists of medium to fine grained heterolithic sandstones and mudstones from estuarine channels and tidal flats. The reservoir quality of this unit is poor. The hundreds of wells (SAGD well pairs and vertical observation wells) drilled into the real field provided many log data, from which five geological facies were identified (petrophysics and petro-elastic properties of those facies are given in Tab. 1): three sandy facies with good reservoir properties and two shaly facies with poor reservoir properties. Their proportions were estimated from logs and first assumed to be constant per layer. Log data were also used to describe the spatial variability of facies. A single variogram was inferred for all facies within each unit. Given this information, we applied the truncated Gaussian method (Le Loc'h and Galli, 1997) to generate the reference facies model. A cross-section perpendicular to wells was extracted from the reference facies realization: it is shown in Figure 8a in Section 4.2. For simplicity, each facies was attributed constant porosity and permeability values (Tab. 1). Following the work of Lerat et al. (2010), a thermal fluid flow simulation was performed at this fine scale, thus resulting in reference production responses. Note that the resulting simulated behavior was pretty close to the one observed for the real field. The flow simulation also yielded pressures, saturations and temperatures over the whole grid at successive times. The initial porosities being known, they were updated on the basis of a geomechanical simulation. Then, we applied the petro-elastic model to compute the corresponding P-wave impedance grids. This step was performed at the geological scale, thus resulting in seismic attributes at the geological scale. Therefore, we do not perform any seismic upscaling in this study, as the seismic and geological scales are the same.
Because of our concern for realism, we did not want the reference seismic attributes to be the straightforward results of the modeling workflow. Thus, we tried to go one step further and determined the seismic traces from the convolution of a seismic wavelet with the reflectivity coefficients (Lerat et al., 2010). At this point, we considered the previously estimated P-wave impedance grids as unknown and the seismic traces as real seismic data. We derived a prior velocity model from the velocities simulated at wells for the reference geological model and inverted the seismic traces following the Aki-Richards approximation (Aki and Richards, 2002). The inversion step (Delépine et al., 2010) yielded the P-wave impedances, which are considered as the reference seismic attributes. The seismic trace modeling/inversion step ensures that the reference impedances differ from the impedances provided at once by the petro-elastic model. Here, no supplementary noise was added to the data: uncertainty only comes from the inversion process.

Figure 1
Illustration of the ratio method in a stationary case.
Therefore, the reference case is characterized both by production data and seismic attributes. More precisely, production data include steam injection rate at the injector and pressure at the producer (Fig. 2, dark diamonds). The seismic data set is based upon three seismic vintages collected at time T 0 prior to any production and after one and three years of production (denoted by time T 1 and T 2 , respectively). The inversion of these seismic data resulted in P-wave impedances. The growth of the steam chamber being very restricted at time T 1 , the seismic attributes estimated for the same time provide only limited information about the spatial distribution of geological heterogeneities. This is the reason why we a) Steam injection rate and b) pressure at producer for the initial model (blue curve), the optimal model (red curve) and the reference case (dark diamonds). preferred to focus on the match of the P-impedances derived from the second monitor campaign (Fig. 3, middle column). The impedance changes along the well pair evidence the influence of shaly heterogeneities: steam chamber expansion is narrowed at the heel (Fig. 3b, slice 5, the steam chamber growth is at the same height as the well perforation), probably because of thin shale lenses. Such heterogeneities are usually so thin (< 3 m) that they cannot be detected from the baseline seismic alone.

MATCHING METHODOLOGY
The purpose of the matching process is to eventually identify a reservoir model respecting all available data: production data and seismic attributes. A starting model guess is sequentially varied so as to decrease the data mismatch. At each iteration, the uncertain parameters are adjusted depending on the gradient values estimated from the optimizer and the sequence of simulation components, which create the forward modeling workflow (see Sect. 1.1), is run. It is worthwhile emphasizing that this modeling workflow departs from the one applied to define the reference case. There are four main differences. First, the number of facies is reduced. As stated above, the reference case encompasses 3 sands and 2 shales. As thermal fluid flow simulation is CPU-time intensive, we reduce the number of unknowns, hence the number of iterations required for minimizing the data mismatch. The three sands, whose petrophysical and elastic properties are similar, are grouped into a single mean sand facies. Therefore, the successive geological models built during the matching phase are populated by three facies: one sand and two shales. Second, the forward modeling workflow repeated at each iteration stops right after the petro-elastic model: the seismic traces are never considered. Third, the geo-mechanical effects are neglected since the Variations of P-impedances at time T 2 for the initial model (hand-left), the reference case (middle) and the optimal model (hand-right). Cross-sections 5, 17 and 27 are selected perpendicularly to the well pair, at the heel, in the middle and at the toe. required simulation CPU-time is prohibitive. Fourth, still to reduce the flow simulation CPU-time, the fine 100 × 41 × 98 geological model is upscaled and flow simulation is performed over a 21 × 41 × 42 grid. The computation of the objective function thus boils down to the comparison of reference production data generated for a 100 × 41 × 98 grid with production responses simulated for the 21 × 41 × 42 grid. These differences make it more difficult to determine a matched model, but emphasize the realism of the case studied. Overall, we aim at determining a matched model and to investigate how accurately it reproduces the reference data. A particular attention is paid to seismic attributes as they yield information about the shape of the steam chamber. The chamber growth for the reference case was shown to be poor at the heel (slice 5) and good at the toe (slice 27) of the well pair (Fig. 3, middle column). Is the matched model able to capture this non uniform development?

MATCHING RESULT
The matching approach developed hereafter includes two successive steps: first, the match of production data starting from the initial model (M 0 ) and resulting in an intermediate model (M 1 ), and second, the match of both production data and seismic attributes, starting from the intermediate model (M 1 ) and resulting in the optimal model (M 2 ). The reference data to be matched are the synthetic data built as explained above.

Match of Production Data
For this preliminary matching phase, the parameters to be adjusted are the horizontal permeability, the vertical to horizontal permeability ratio, the vertical correlation length characterizing the spatial facies variability in the bottom unit and the sandstone proportions in several regions. Initially, these proportions are assumed to be identical. Then, they are allowed to vary independently from each other. This feature adds degrees of freedom to the matching process and makes it more efficient. Facies proportions are varied using the ratio method recapped in Section 1.5. The ratio defined as an optimization parameter is the ratio of sand proportion over the sum of sand and both shales proportions. Then, we apply the truncated Gaussian method with the proportion matrix defined according to the ratio. Proportions are independently changed in each of the selected sub-domain.
The data to be matched in this first phase are the production data only: they are not informative enough to drive the definition of the sub-domains. Therefore, we just split the bottom unit of the reservoir model into 6 regions with roughly the same volume along the well pair (see Fig. 4). One parameter for varying the sandstone proportion is defined for each region, thus, at last, this history matching process is driven by nine parameters.
The reference production data are the steam injection rate and the pressure in the producer (see Sect. 2). The data mismatch is expressed following the usual least-square formulation (Eq. 2). The weights were chosen to properly balance the initial contribution of the injection rate and pressure terms. The history matching process can be stopped according to several criteria. Either the convergence is obtained, i.e., the objective function stops decreasing, or the maximal number of objective function evaluations is reached. Here, the optimization stopped due to the first criterion after 30 flow simulations: the objective function was reduced by 87%. Figure 2 shows that the production responses simulated for the resulting matched model (red curve) are in good agreement with the reference data (dark diamonds).
The changes applied to obtain the matched model (M 1 ) are an increase of the horizontal permeability of sand (denoted by Kh), a decrease in the ratio of vertical permeability over horizontal permeability for sand (denoted by Kv/Kh) and an increase in the vertical correlation length (denoted by Lv). Finally, the ratio of the sand proportion (denoted by p_Z1 to p_Z6 for regions 1 to 6) was increased in five of the six selected areas. The variations in uncertain parameters are reported in Table 2. Although not included into the objective function, the seismic attributes were also simulated for the matched model (M 1 ) (Fig. 3, left-hand column). Clearly, they differ from the reference seismic attributes (Fig. 3, middle  column). This illustrates that production data match does not ensure the seismic attribute match. An additional matching step is required.

Match of Both Production Data and Seismic Attributes
The second matching phase, which involves the seismic attributes on top of the production data, starts from the  Parameterization is a key issue for any inverse problem. This is even more decisive when dealing with seismic data. Therefore, before starting any optimization process, we investigated the values of the initial seismic attribute mismatch over the grid, that is the difference between the seismic attributes characterizing the reference case and the ones simulated for the intermediate model (M 1 ). This allowed us to identify cross-sections with poor seismic match and for customizing the parameterization to focus on these crosssections. The reference seismic attributes (Fig. 3, middle  column) show that the development of the steam chamber is very restricted at the heel of the well (slice 5) probably because of intercalated shale baffles. This behavior is not reproduced by the intermediate model (M 1 ) (Fig. 3, handleft column). Therefore, in the bottom unit, the first five vertical slices at the heel are split into 10 sub-domains in the bottom unit (Fig. 5). The sand proportion in each subdomain is considered as a parameter to be adjusted. In addition, the steam chamber conformance is good at the toe for the reference case, but not for the intermediate model (M 1 ) (Fig. 3). Thus, a last sub-domain is identified from cross-sections 27 to 29 (Fig. 5) with sand proportion selected as an uncertain parameter. All around, facies proportions are unchanged even if the seismic match is not perfect. Again, facies proportions are varied using the ratio method with a ratio defined as proportion of sand over proportion of all other facies. Finally, the second matching process is driven by 11 parameters.
As explained in Section 2, the objective function includes the production data and the seismic attributes (variations in P-impedances between times T 0 and T 2 ). At this stage, the reference P-impedances are expressed in time domain. The predicted data are converted in the time domain before evaluating the data misfit as explained in Section 1.3. The production data term is kept to prevent the final model (M 2 ) from departing from the already matched reference data. The E Tillier et al. / Simultaneous Inversion of Production Data and Seismic Attributes: Application to a Synthetic SAGD Produced Field Case 297  Figure 5 Sub-domains (one different color for each sub-domain considered) where sand proportions are varied for matching seismic attributes. The well pair trajectory is represented by the blue line. The rectangle represents the bounding box of the reservoir grid.
A preliminary matching attempt was carried out with the usual least-square formulation for quantifying the data misfit for both production data and seismic attributes. It induced an objective function decrease of 3% only. We could check that the seismic impedances simulated for the final matched model were the same as the ones for the intermediate model (M 1 ). The optimal values of uncertain parameters are quite the same as the initial ones (Tab. 3). Many reasons can explain why this history matching failed. First, the drastic upscaling between the reservoir grid used for the history matching and the reference grid used for building the synthetic data implies that it is impossible to obtain a zero misfit. Nevertheless, we were expecting to retrieve the appropriate behavior. Second, specification of the weighting coefficients in the objective function heavily influences the final matched model and the performance of the optimization algorithm. We conducted an extensive study where several values for the weights were investigated, but we have never been able to improve the matching process. In this case, we concluded that the history matching mainly failed due to the way of measuring the misfit between reference and simulated seismic attributes. This stresses that an alternative formulation must be used to measure the seismic error. We thus moved to weight attributed to production data is small and the seismic attributes mismatch is the major contribution to the objective function. The production data mismatch being initially very low, it strongly impacts the objective function as soon as it departs from its low value, whatever the weight it is assigned. the LDM formulation introduced in Section 1.4 and repeated the same matching process, thus determining the optimal model (M 2 ). The way of measuring misfit is the only change made into the history matching process. The parameterization is unchanged as well as the way of obtaining the predicted data. The history matching results obtained are then much better. The objective function now decreases by 30%, which is very significant when accounting for seismic attributes. The production data match is still good while the seismic attributes match is clearly improved. The seismic attributes simulated for the so-obtained optimal matched model (M 2 ) reproduce the reference steam chamber conformance along the well pair (Fig. 3, middle and hand-right columns). Pressure, gas saturation and temperature maps are displayed for three slices of the intermediate model (M 1 ) (Fig. 6) and of the optimal model (M 2 ) (Fig. 7). These results could be achieved thanks to the flexibility of the parameterization used for describing the spatial distribution of facies. Figure 8 shows that the matching process resulted in more shale at the heel (see red circle in Fig. 8c) and less at the toe (see green circle in Fig. 8c). The initial and optimal values for the ratio of facies proportions are given in Table 3. Facies did not change all around the target sub-domains since the selected parameterization makes it possible to adjust proportions in given sub-domains only (Fig. 5). We refer here to the sub-domains defined during the parameterization step. However, it generates discontinuities between the target subdomains and the embedding medium. This feature could be avoided using a similar, but more time-consuming parameterization method (Tillier et al., 2010) as explained in Section 1.5.
Last, we developed the same matching methodology, but with the reference impedances given in depth domain instead of time domain. In other words, a preliminary timeto-depth conversion was carried out before proceeding to any matching process so that the reference impedances were shifted over the reservoir grid. The matching process performed with the seismic attributes in depth domain resulted in a 37% objective function decrease. It was shown that the impedances computed for the so-achieved matched model were almost the same as the ones identified from the impedances in time domain. Such a result can be explained. For the case studied, the error due to time-to-depth conversion is approximately 0.6 m, which is much less than the total reservoir thickness (50 m). This error is surely negligible compared to others (low seismic resolution, inversion of seismic traces, etc.). In the context of this study, the time-to-depth conversion is almost perfect and thus performing the history matching in one or the other domain as no influence. As a result, the match obtained with this second approach is neither better nor worst than the one obtained with the first approach.  Temperature, pressure and gas saturation maps for the optimal model (M 2 ) given on three vertical slices perpendicular to the well.

CONCLUSIONS
In this paper, we presented a matching methodology for identifying reservoir models respecting data with very distinct natures: production data at wells and seismic attributes over the whole reservoir grid. Then, it was successfully applied to a synthetic, but realistic SAGDproduced case. This study evidenced the need for an alternative formulation to measure the mismatch between the simulated and reference seismic attributes. The usual least-square formulation was shown to be unsuitable. Therefore, we developed a formulation relying on the local modified Hausdorff distance, which turned out to be efficient.
The parameterization strategy was shown to be a key to success. A good match could be achieved because the selected parameterization method allowed us for varying facies within the reservoir model. The variations were localized in target sub-domains and driven from a reduced number of parameters.
Last, we addressed issues related to time-to-depth conversions. The results obtained for the case studied showed that the time-to-depth conversion error was quite small with respect to others. Therefore, matching seismic attributes in depth domain or time domain was more or less equivalent in terms of results. However, considering the impedances in time domain as the data to be matched makes the matching methodology more consistent.