SPE-125101-MS-P - SPE 125101 Data Assimilation Using The...

Info iconThis preview shows page 1. Sign up to view the full content.

View Full Document Right Arrow Icon
This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: SPE 125101 Data Assimilation Using The Constrained Ensemble Kalman Filter Hemant A. Phale, Dean S. Oliver, SPE, University of Oklahoma Copyright 2009, Society of Petroleum Engineers This paper was prepared for presentation at the 2009 SPE Annual Technical Conference and Exhibition held in New Orleans, LA, USA, 4–7 October 2009. This paper was selected for presentation by an SPE program committee following review of information contained in an abstract submitted by the author(s). Contents of the paper have not been reviewed by the Society of Petroleum Engineers and are subject to correction by the author(s). The material does not necessarily reflect any position of the Society of Petroleum Engineers, its officers, or members. Electronic reproduction, distribution, or storage of any part of this paper without the written consent of the Society of Petroleum Engineers is prohibited. Permission to reproduce in print is restricted to an abstract of not more than 300 words; illustrations may not be copied. The abstract must contain conspicuous acknowledgment of SPE copyright. Abstract When the ensemble Kalman filter (EnKF) is used for history matching, the resulting updates to reservoir properties sometimes exceed physical plausible bounds, especially when the problem is highly nonlinear. Problems of this type are often encountered during history matching compositional models using EnKF. In this paper, we illustrate the problem using an example in which the updated molar density of CO2 in some regions is observed to take negative values while molar densities of the remaining components are increased. Standard truncation schemes avoid negative values of molar densities, but do not address the problem of increased molar densities of other components. The results can include a spurious increase in reservoir pressure with a subsequent inability to maintain injection. In this paper, we present a method for constrained EnKF which takes into account the physical constraints on the plausible values of state variables during the data assimilation. The proposed method can be implemented in two different approaches, both of which convert inequality constraints to a small number of equality constraints. The first approach uses Lagrange multipliers to apply the active constraints. In the second approach, the constraints are used as virtual observations for calibrating the model parameters within plausible ranges. The constrained EnKF method is tested on a 2D compositional model and on a highly heterogeneous 3-phase flow reservoir model. The effect of the constraints on mass conservation are illustrated using a 1D Buckley-Leverett flow example. Results show that the constrained EnKF is able to enforce the nonnegativity constraints on molar densities and bound constraints on saturations (all phase saturations must be between 0 and 1), and achieve a better estimation of reservoir properties than is obtained using only truncation with EnKF. Introduction The reliability of reservoir models generally increases as more data are included in their construction. Geostatistical methods have traditionally been used to incorporate sparsely located static data (e.g. geological, geophysical, and well log or core data) into the construction of reservoir geological models (Deutsch and Journel, 1998). These methods, however, are typically limited by the requirement that the observations are linearly related to the model variables. On the other hand, dynamic production data (e.g. oil production rate, water production rate, gas oil ratio, and cumulative fluid production rates) are available in abundant quantity but are nonlinearly related to model variables. Methods for incorporating both static as well as dynamic data in geological reservoir model construction are often termed history matching. Assisted or automatic history matching is formulated as a minimization problem in which parameters of the model are adjusted to minimize the mismatch between actual observations and observations predicted from the model. In most cases a regularization term is included in the minimization function to penalize deviations from plausible values (Tarantola, 1987; Oliver et al., 2001; Li et al., 2003). If the data are assimilated sequentially in time, the process is called data assimilation. In recent years, the ensemble Kalman filter (EnKF) technique has established itself as a viable data assimilation method for models with large numbers of variables (Evensen, 1994). The EnKF relies on the Monte Carlo approach to forecast the error statistics and to compute an approximate Kalman gain matrix for updating model and state variables. In a recent review paper, Aanonsen et al. (2009) described a number of applications where the EnKF has been used for problems specific to the petroleum industry. The first application of the EnKF within the petroleum industry was presented by Lorentzen et al. (2001) and the first application to reservoir monitoring was presented by 2 SPE 125101 Nævdal et al. (2002). Since its introduction to petroleum engineering field, EnKF technique has gained increasing attention for history matching and continuous reservoir model updating (Nævdal et al., 2005; Haugen et al., 2006; Bianco et al., 2007; Evensen et al., 2007; Skjervheim et al., 2007; Zhang and Oliver, 2009). Despite the general success, one problem that can occur when using the EnKF to assimilate reservoir data is that the updated saturations may take values which are negative or are greater than unity. Gu and Oliver (2006) showed examples in which this problem occurred and proposed three different methods for preventing nonphysical water saturations. Using normal score transformed values of the saturations as state variables ensures that saturation values can never exceed plausible limits but it did not prevent some nonphysical oscillations in the values. A normal score transformation of one set of state variables (termed Gaussian anamorphosis) was reported to give good results in an oceanographic application (Simon and Bertino, 2009). A much different type of reparameterization was applied by Chen et al. (2009) in an attempt to limit problems with the updates of saturation. Although the method was successful in improving the quality of the updates, it required substantial additional computation. The most generally successful method are iterative forms of the ensemble Kalman filter (Gu and Oliver, 2007; Li and Reynolds, 2009) in which the Kalman gain is used only to update model variables. Because state variables are updated using the simulator they naturally satisfy the inequality constraints. The cost of iteration can be substantial, however, so the goal should be to avoid recourse to this solution. One approach that has sometimes proved successful with the Kalman filter is explicit incorporation of the constraints into the analysis step. Equality constraints have been addressed by a number of authors in the context of the Kalman filter (Qureshi et al., 1989; De Geeter et al., 1997; Simon and Chia, 2002; Ko and Bitmead, 2007) but while these approaches could generally be easily extended to the ensemble Kalman filter, the problem of inequality constraints is generally more important for data assimilation and certainly more challenging. The problem of applying inequality constraints to estimates of state variables has also been addressed (Massicotte et al., 1995; Rao et al., 2003; Gupta and Hauser, 2007; Rotea and Lana, 2008). Ungarala et al. (2007) presented a formulation for incorporating inequality constraints on state variables in an extended Kalman filter technique. In their approach, active constraints were identified by an exhaustive search of all combinations of state variables. Their approach was successfully applied to a non-linear batch reactor with small number of state variables and limited number of constraints, but direct application to reservoir problems with few hundred thousands of state variables and constraints is not feasible. Thacker (2007) discussed two different approaches for incorporating inequality constraints into a Kalman filter framework. In his investigation, it is determined that when the number of state variables is relatively large and correlated, it may not always be important to determine the active set of constraints precisely. As a result, the extension of the approach to the ensemble Kalman filter is relatively straightforward. In this paper, we describe a new method based on the work by Thacker (2007) for updating reservoir models using the ensemble Kalman filter while avoiding the problem of non-physical values for saturations and densities. This paper includes three test problems designed to validate the performance of the proposed method against the standard implementation of the ensemble Kalman filter. The first test problem shows the benefits of the constrained solution in a compositional model for which truncation of negative densities results in significant problems with the mass of fluid in the model. The second test problem is a three-phase black oil example for which the state variables have both lower and upper bounds. The final example is a simple 1D Buckley-Leverett flow example that investigates the effects of constraints on mass conservation during the update step. Methodology The ensemble Kalman filter The EnKF methodology is a two step process. The first step is a forecast step in which the system is evolved forward in time. The second step is a data assimilation step in which the model and state variables that describe the system are updated to reduce the observation data mismatch. At a data assimilation time tk , each state vector of the ensemble is updated using the ensemble approximation to the Kalman gain in the following manner: f f u yk,j = yk,j + Kk (dobs,k,j − Hk yk,j ), for j = 1, 2, . . . , Ne (1) where, the superscripts f and u stand for forecast and update, respectively; the subscript k stands for the time index; and H is the measurement operator that maps the the state vector yf to the simulated data. The forecast state is obtained by running the simulator forward in time up to the time at which observations are available. In Eq. 1, Kk is the ensemble approximation to the Kalman gain and is computed as: Kk = Cf HT (Hk Cf HT + CD,k )−1 Y,k k Y,k k (2) The collection of all ensemble state vectors is denoted by Y and is a matrix of dimensions Ny × Ne . In Eq. 2, SPE 125101 3 Cf denotes the prior covariance matrix for the state variables which can be estimated from the ensemble using Y,k the following expression: 1 f f f Cf = (Yf − Yk )(Yk − Yk )T (3) Y,k Ne − 1 k f where, Yk is the mean of state variables calculated across the ensemble members and is a vector of dimensions Ny . The state covariance matrix can be very large, but it is never necessary to compute the entire matrix. Instead, f f we compute products of the form (Yk − Yk )T HT which are typically much smaller in dimension than the state k covariance matrix. The Constrained ensemble Kalman filter The kinds of constraints that are of particular interest in this paper are linear inequality constraints of the form ymin ≤ yi ≤ ymax which is of the type satisfied by phase saturations and molar fractions, or of the form yi ≥ ymin which is of the type satisfied by pressures and densities. It is not easy to determine which of the inequality constraints should be converted to equality constraints and which might be ignored when the problem has many variables. A relatively easy approach would be to first assimilate the data without any constraints (data assimilation using the standard EnKF) and then check for constraint violations. In this paper, we apply an iterative procedure for identifying the constraints that are active in a particular data assimilation step. Once the active set is identified, equality constraints are applied to the problem in the update step. If constraints are still violated, the procedure is repeated. When it is determined that a state variable is on the boundary of the feasible region, then the inequality constraint can be replaced by equality constraints. The system of active constraints can be written in the following form, Py = c (4) where, the vector c contains the constraint values. Thacker (2007) discussed two different approaches for incorporating equality constraints into a Kalman filter framework. The first approach uses the undetermined Lagrange multipliers and it is shown that, the resulting formulation has a form very close to the Kalman filter updating equation. The second approach incorporates the constraints along with data in an expanded data vector and subsequently, the entire system of equations is expanded to reflect the constraints. Following is a brief description of the two approaches mentioned above. Lagrange multipliers This approach uses the undetermined Lagrange multipliers (Thacker and Long, 1988) for incorporating equality constraints into the data assimilation formalism. While the standard Kalman filter computes variables that minimize J(y) = 1 1 (y − yf )T (Cf )−1 (y − yf ) + (Hy − d)T C−1 (Hy − d) Y D 2 2 (5) the constrained Kalman filter will compute variables that minimize L(x, y) = J(y) − (Py − c)T x. (6) The vector x is a vector of the Lagrange multipliers that enforce the constraints. The values of the Lagrange multipliers must be determined as a part of the solution. The constrained minimum solution can be written as a correction to the unconstrained updated minimum solution yu (Thacker, 2007): u,con u u yj = yj + Cu PT (Pj Cu PT )−1 (cj − Pj yj ) (7) Y j Y j u,con u where, yj and yj are the constrained and unconstrained updated solutions for the jth ensemble member, u respectively; CY is the posterior covariance matrix; Pj is the coefficient matrix for the constraint equations for the jth ensemble member; and cj is the vector of constraint values for jth ensemble member. The constrained solution (Eq. 7) is clearly similar in form to the unconstrained solution from the EnKF (Eq. 1). The largest difference is that, while in the standard EnKF the same Kalman gain is used to update each ensemble member, here each ensemble member may have a different set of active constraints, so terms such as Pj , Cu PT , Y j u Pj Cu PT , cj , and Pj yj must be computed for each ensemble member at each data assimilation time step in order Y j to compute the constrained solution. 4 SPE 125101 Augmented data This approach also starts by computing the unconstrained solution. Once the unconstrained solution is obtained, the next step includes checking for any constraint violations that need to be enforced. The constraint data cj are then appended to the data vector dobs,j . The data-error covariance matrix CD should be enlarged with zero values on the diagonal corresponding to the constraint errors, and H should be enlarged to reflect the constraints. With these modifications, an augmented Kalman gain that includes columns corresponding to data and active constraints can be constructed for computing the constrained updates to the original forecast f yj . The modified form of Eq. 1 for implementing the constrained EnKF method can be given as: f u,con ˜ ˜ ˜ ˜ ˜ ˜ f = yj + Cf HT (HCf HT + CD )−1 (dobs,j − Hyj ) yj Y Y (8) ˜ where, HT is the modified form of the data extraction matrix which extracts the state variables that are subject ˜ to constraint enforcement along with the simulated data; CD is the modified form of the data-error covariance ˜ obs,j is the modified vector of the observations appended with the matrix which also reflects the constraints; and d constraint values. As we can notice from Eq. 7 and Eq. 8, the primary difference between the two approaches is in the sequence of incorporating the equality constraints. In the Lagrangian approach, the constraints are applied after assimilating data while in the augmented data approach, the constraints are assimilated simultaneously with data. For linear observation operators and Gaussian probability densities, the two approaches are equivalent. In the examples that follow, we use the augmented data approach. Algorithm description for the constrained EnKF Following is a brief generalized algorithm for implementing the augmented data approach of the constrained EnKF method. This algorithm describes various steps which need to be carried out in order to enforce the constraints on any type of state variable. 1. At each data assimilation time step, compute the unconstrained solution for each ensemble member using the standard EnKF methodology as given by Eq. 1. 2. Identify the state variable with the largest degree of constraint violation. 3. Append the constraint to the vector of actual observations and modify the observation operator matrix to include the constraint relation. Also expand the data-error covariance matrix, CD for the constraint. ˜ ˜ 4. With the newly expanded matrices H and CD , compute the modified Kalman gain for the particular ensemble member under consideration. 5. Update the state vector for that ensemble member using Eq. 8. 6. If the state variables of the updated realization still violate the inequality constraints by a significant amount, repeat steps 2 through 5. During each additional iteration, expand all the matrices and vectors in order to reflect the new constraints while keeping the old constraints from the previous iterations. 7. Repeat steps 2 through 6 for all ensemble members. 8. Repeat steps 1 through 7 for all the data assimilation time steps in the history matching period. Application to a Compositional Model In this section, we describe the application of constrained EnKF to a data assimilation problem of CO2 injection into a compositional reservoir model with 12 components, 11 in the hydrocarbon phase. The number of components is set to be large in order to test the robustness of the method for difficult problems. Results obtained from the implementation of the standard EnKF will be discussed first. These results illustrate problems that can occur with the application of the standard EnKF when the state variables must satisfy inequality constraints. Results from the implementation of the constrained EnKF method are presented subsequently. Reservoir Model Description The synthetic 2-dimensional (x-z cross-section) compositional reservoir model has immobile water saturation with no initial gas phase. The 11 components in the oil phase range from C1 to C7+ along with CO2 and N2 . There is one producer and one injector in the field. The producer is located on the east boundary of the field and is operated with a bottom hole pressure (BHP) constraint of 1600 psia. The CO2 injector is located on the west boundary of the field with a BHP constraint of 2200 psia. The reservoir model is divided into a uniform grid of dimensions 50 × 1 × 30. Grid block dimensions are 50 ft × 30 ft × 20 ft. The true porosity and log permeability are realizations of a correlated anisotropic Gaussian random field. The mean porosity is 0.32 SPE 125101 5 with standard deviation of 0.06. The mean log permeability is 4.0 with standard deviation of 1.4. Figs. 1(a) and 1(b) show the porosity and log permeability fields for the reference case, respectively. Completion locations are also shown in Fig. 1(b). The compositional suite from a commercial reservoir simulator ECLIPSE (Schlumberger, 2007) is used for the forward modeling of flow and transport. Injector -5 0 Porosity 5 Z-direction Producer 0.48 0.44 0.40 0.36 0.34 0.32 0.28 0.24 0.20 0.16 10 15 20 25 30 35 -10 0 10 20 30 40 50 60 X-direction (a) Reference: Porosity (b) Reference: Log Perm Figure 1: The porosity and log permeability fields for the compositional model (reference). Ensemble Kalman filter Forty initial realizations of porosity and log permeability are generated using the same geostatistical parameters as the reference model. None of the realizations are conditioned to static well data. Measurements used in the data assimilation include oil production rate, producing gas-oil ratio, CO2 injection rate, gas-phase CO2 mole fraction at the producer, and liquid-phase CO2 mole fraction at the producer. An additional observation of CO2 molar density at grid block (16,22) is assimilated at the first data assimilation time. The state vector for model updating includes the molar densities of 12 components from the reservoir fluid system, pressure, porosity and log permeability. Thus, there are a total of 15 variables in each grid block. The total simulated production history time is 1106 days and production data are assimilated approximately every 100 days. The noise associated with the production data is sampled from a normal distribution with a mean of 0 and standard deviation which varies by type of observation. Following is a list of the standard deviation values used for different types of production data: • oil production rate: 10.0 STB/day • CO2 injection rate: 20.0 lb moles/day • producing gas-oil ratio: 0.50 MSCF/STB • producing gas-phase CO2 mole fraction: 0.0005 • producing liquid-phase CO2 mole fraction: 0.00005 Figure 2 shows the distribution of CO2 throughout the reservoir at day 102 before data assimilation. Because of the relatively large degree of heterogeneity in permeability, the molar densities are highly variable. Results from standard EnKF At the first data assimilation time (day 102), in addition to the observations of production data, we assume that the CO2 molar density is measured at one location in the reservoir at that time. Fig. 3 shows the CO2 molar densities for the 22nd layer of the reservoir model before and after data assimilation. The updated CO2 molar density from the analysis step in EnKF (Fig. 3(b)) is negative in some locations for most of the realizations. The standard solution in the EnKF is to truncate state variables that violate constraints to the value at the boundary of the feasible region. Molar density cannot be negative so all CO2 densities in the model that are less than 0 would be set equal to 0 at the end of the analysis step. When that approach is taken in this case, the consequences on reservoir behavior are striking and are most easily understood through investigation of the behavior of one particular realization. We will examine the effect of the use of the standard EnKF update with truncation on the reservoir behavior of realization #31. At the end of the first data assimilation step, the values of all state variables are updated and returned to the reservoir simulator. Fig. 4 shows the updated pressure and molar densities of two of the components. The 6 SPE 125101 (a) Realization 6 (b) Realization 18 (c) Realization 32 Figure 2: Forecast CO2 molar density distributions at t=102 days for various ensemble members. 3 Molar density Molar density 2.0 1.5 1.0 0.5 0.0 10 20 30 40 50 2 1 0 1 Grid block index 10 20 30 40 50 Grid block index (a) Forecast CO2 molar density (b) Updated CO2 molar density (without truncation) Figure 3: The prior and posterior ensemble of molar density of CO2 for the 22nd layer. pressure after updating is quite reasonable with high pressure on the left and low pressure on the right near the producing completions. The molar density of NC4 and C7 are somewhat odd, however, as they show anomalously high densities to the right of the measurement location. These densities are higher than at any other location in the reservoir. (a) Pressure: Day 102 (at restart) (b) NC4 molar density: Posterior (c) C7 molar density: Posterior Figure 4: Realization #31 - Pressure and molar density distributions in the reservoir after the first data assimilation step. The model is restarted using all the static and dynamic variables updated at the first data assimilation time step (day 102). At day 102.5, a region of high pressure can be seen between the injector and the producer. The magnitude of the maximum pressure in this region is about 2700 psia which is substantially higher than the maximum allowable bottom hole pressure for the injector (2200 psia). The anomalously high pressure due to the excess mass is significant enough that it results in an inability to inject fluid into the reservoir (Fig. 5(b)). The SPE 125101 7 CO2 injection rate lb moles day regions with higher molar densities of components NC4 and C7 in Figs. 4(b) and 4(c) overlap quite well with the region of pressure build-up shown in Fig. 5(a). Similar observation of increased molar density is made for all of the other components from the oil phase with the exception of CO2 . 5000 0 5000 10 000 15 000 20 000 25 000 30 000 0 200 400 600 800 1000 1200 Time Days (a) Day 102.5251 (immediately after restart) (b) CO2 Injection Rate (c) Legend Figure 5: Pressure in Realization #31 at the end of the first forecast time step after data assimilation (left) and the ensemble of injection rates after data assimilation. All of the components in the oil phase for realization 31, except CO2 , exhibited an anomalous increase in the molar density after the update step. The increase in mass resulted in a subsequent increase in the pressure which then resulted in other problems. Truncation of the negative CO2 molar densities contributed to the mass excess problem by reducing only the amount of negative mass change without affecting the positive changes. In the next section, we show that the constrained EnKF solves the problem of state variables that violate constraints, without introducing the problem of excess mass. Results from constrained EnKF The constrained EnKF uses a different Kalman gain matrix and a different observation vector for each ensemble member. Here we describe in detail the application of the constrained EnKF to the problem of data assimilation for the same realization studied in the previous section. For the first data assimilation time step, the unconstrained solution for realization 31 showed a value of -0.5632 for the updated CO2 molar density at the grid block (27,23) as the largest violation of the positivity constraint on CO2 molar density value. Fig. 6 shows the results for the forecast, unconstrained and the constrained solutions for the components CO2 , NC4 , and C7 obtained from iteration 1 for the 23rd layer in the model. The green vertical line indicates the location of the grid block (27,23) at which the constraint is enforced. At subsequent iterations, the largest constraint violation occurs in other layers, but there is still a small effect from the application of constraints on CO2 molar density in other layers on the density in layer 23 as can be seen in Fig. 6. In the regions where the CO2 molar density takes negative values, the unconstrained EnKF analysis step gives anomalously high molar densities for all other components in the hydrocarbon phase. Fig. 6 shows that the consequence of applying the constraints on the molar density of CO2 is to reduce the increase in molar densities of components NC4 and C7 . The constrained EnKF technique helps to ensure that the magnitudes of the increases and decreases in molar densities are applied in a consistent manner. This consistency is missing in case of the simple ad hoc truncation scheme used during the implementation of the standard EnKF. The results showed that the constrained EnKF method is able to solve the problem of increased molar density (and thereby, additional mass) of components in a particular region of the reservoir field which was responsible for the high pressure build up in that region ultimately shutting the injector off. Fig. 7 shows the CO2 molar density values with maximum degree of constraint violation within the entire reservoir model for realization 31 at different iterations of the constrained EnKF method for some of the data assimilation time steps. With the application of the constrained EnKF method in an iterative manner, the magnitude of the CO2 molar density with the maximum degree of constraint violation moves in the direction of the enforced constraint value. SPE 125101 2.5 2.0 1.5 1.0 0.5 0.0 0.5 1.0 0.4 Molar density 0.4 Molar density Molar density 8 0.3 0.2 0.1 0.0 0 10 20 30 40 50 0.3 0.2 0.1 0.0 0 10 Grid block index 20 30 40 50 0 10 20 30 40 Grid block index (b) NC4 (a) CO2 50 Grid block index (c) C7 Figure 6: Layer 23: The forecast, unconstrained (standard EnKF), and constrained (CEnKF) solutions obtained at the 1st data assimilation time step using iteration 1 for the 31st ensemble member. The black curve stands for the forecast solution, red curve denotes the unconstrained EnKF update, and the blue curve represents the constrained EnKF update. 0.00 CO2 molar density CO2 molar density 0.1 0.2 0.3 0.4 0.5 0.6 0 2 4 6 iteration index (a) assimilation time step 1 8 10 0.05 0.10 0.15 0.20 0 2 4 6 8 10 iteration index (b) assimilation time step 4 Figure 7: Realization #31: The minimum value of CO2 molar density at different iterations of the constrained EnKF method for some of the data assimilation time steps. Matching Production Data The final updated ensembles of porosity and permeability fields obtained from the standard EnKF and the constrained EnKF method are run from the initial time (day 0) up to the end of the prediction period (day 1735). Fig. 9 and Fig. 10 show the comparison of the history match of different types of production data from the two methods. As the standard EnKF method suffered from the issue where most of the ensemble members had problem with continuous fluid injection during the data assimilation process, it had adverse effect on updating the static model parameters. This is reflected in the history matching results obtained from the standard EnKF which are shown in Fig. 9. The results for CO2 injection rate obtained from the standard EnKF method are biased and do not match the historical injection rate data. Similarly, the results for the oil production rate obtained from the standard EnKF (Fig. 9(a)) show bias in the early production period where most of the ensemble members produce significantly more amount of oil than the observations. The constrained EnKF method gave significantly better history matching results compared to the standard EnKF. The results for the history match of the gas-phase CO2 mole fraction (WYMF) data are shown in Fig. 10(c) and Fig. 10(d) for the two methods. The results indicate that the constrained EnKF resulted in a better match to the CO2 breakthrough time, compared to the standard EnKF. The root mean squared error (RMSE) values for the production data from the two methods are given in Table 1. For all types of data, results from the constrained EnKF are better than the results from the standard EnKF. Application to a Black Oil Model In this section, the results from the implementation of the ensemble Kalman filter in its standard form and the implementation of the constrained EnKF method (augmented data approach) to a data assimilation problem of black oil model are discussed. The reservoir model description is given below. Reservoir Model Description The synthetic 3-phase, 2-dimensional (x-y areal) black oil reservoir model has connate water saturation of 0.15 with no initial gas saturation. The initial reservoir pressure is 3600 psia, which is 5000 0 5000 10 000 15 000 20 000 25 000 30 000 0 200 400 600 800 1000 1200 CO2 injection rate lb moles day 9 CO2 injection rate lb moles day SPE 125101 0 5000 10 000 15 000 20 000 0 200 400 Time Days 600 800 1000 1200 Time Days (a) Standard EnKF (b) Constrained EnKF Figure 8: The CO2 injection rate profiles for all the ensemble members during data assimilation obtained from the standard EnKF and the constrained EnKF 1000 Oil Prod Rate STB day Oil Prod Rate STB day 1000 800 600 400 200 0 0 500 1000 800 600 400 200 0 0 1500 500 Time Days 0 2000 4000 6000 8000 10 000 12 000 14 000 500 1000 1500 (b) Constrained EnKF: OPR 1500 CO2 Injection rate lb moles day CO2 Injection rate lb moles day (a) Standard EnKF: OPR 0 1000 Time Days Time Days 0 2000 4000 6000 8000 10 000 12 000 14 000 0 500 1000 1500 Time Days (c) Standard EnKF: WCMIR (inj) (d) Constrained EnKF: WCMIR (inj) Figure 9: Comparison of the oil production rate (OPR) and CO2 injector rate (WCMIR) data obtained from rerunning the final updated ensembles obtained from the standard EnKF and the constrained EnKF. The black lines represent ensemble members, the red dots represent the observations from the reference model, and the vertical line indicates the end of the history matching period (day 1106). Table 1: Compositional model - RMSE values for different types of observations Observation WOPR WCMIR WGOR WYMF WXMF Standard EnKF 79.3226 1898.52 8.4097 0.1089 0.00185 Constrained EnKF 49.6641 667.61 3.2724 0.0625 0.00106 10 SPE 125101 Gas Oil Ratio MSCF STB Gas Oil Ratio MSCF STB 140 350 300 250 200 150 100 50 120 100 80 60 40 20 0 0 0 500 1000 0 1500 500 1.0 0.8 0.6 0.4 0.2 500 1000 1500 (b) Constrained EnKF: GOR Gas phase CO2 mole fraction Gas phase CO2 mole fraction (a) Standard EnKF: GOR 0.0 0 1000 Time Days Time Days 1500 Time Days (c) Standard EnKF: WYMF 1.0 0.8 0.6 0.4 0.2 0.0 0 500 1000 1500 Time Days (d) Constrained EnKF: WYMF Figure 10: Comparison of the gas oil ratio (GOR) and gas-phase CO2 mole fraction (WYMF) data obtained from rerunning the final updated ensembles obtained from the standard EnKF and the constrained EnKF. The black lines represent ensemble members, the red dots represent the observations from the reference model, and the vertical line indicates the end of the history matching period (day 1106). also the bubble point pressure. The reservoir is divided into a uniform cartesian grid of dimensions 50 × 50 × 1. The grid block dimensions are 30 ft × 30 ft × 20 ft. There are 4 producers and 1 injector in the field arranged in a five-spot well pattern. The porosity field has a mean of 0.25 and a standard deviation of 0.03. The log permeability field is generated with a mean and a standard deviation of 4.8 and 1.4, respectively. An anisotropic Gaussian variogram model with practical ranges of 30 and 5 gridblocks in the two principal directions was used to generate the fields. Fig. 11(a) and 11(b) show the porosity and log permeability fields for the reference case, respectively. All the producers are operated with a maximum oil production rate of 1500 STB/day as the primary constraint and a minimum bottom-hole pressure of 1000 psia as the secondary constraint. The injector, located in the center of the field, is operated with a maximum water injection rate of 8000 STB/day as the primary constraint and a maximum bottom hole pressure of 6000 psia as the secondary constraint. Data Assimilation The porosity and log permeability fields for the reference model and the initial ensemble are generated using the sequential Gaussian simulation algorithm with a correlation coefficient of 0.75. An ensemble of 100 realizations is conditioned to values of permeability and porosity at well locations from the reference model. The state vector includes static model parameters such as porosity and log permeability along with dynamic state variables such as pressure, water saturation, and gas saturation. The simulated production data, including oil production rate, producing gas-oil ratio, water production rate, and water injection rate, are used in the data assimilation. The total simulated production history is 280 days. The production data are assimilated at day 10, 20, 40, and thereafter at every 40 days. The measurement noise associated with the production data is sampled from a normal distribution with a mean of 0 and a standard deviation which varies by type of observation. Following is a list of the standard deviation values used for different types of production data: • oil production rate: 4.0 STB/day • producing gas-oil ratio: 0.003 MSCF/STB • water production rate: 0.0001 STB/day when the true rate is less than 1.0 STB/day; otherwise 4.0 STB/day • water injection rate: 5.0 STB/day SPE 125101 11 0 y-direction 10 0.34 0.32 0.30 0.28 0.26 0.24 0.22 0.20 0.18 PROD3 20 INJ 30 40 PROD1 0 PROD4 10 y-direction PROD4 20 INJ 30 40 PROD2 50 PROD3 PROD1 8.00 7.00 6.00 5.00 4.00 3.00 2.00 1.00 0.00 -1.00 PROD2 50 0 10 20 30 40 50 0 10 20 30 40 50 x-direction x-direction (a) Reference: Porosity (b) Reference: Log Perm Figure 11: The porosity and log permeability fields for the black oil model (reference). In the implementation of the constrained EnKF, the updated water saturation is constrained to the connate water saturation value of 0.15 whenever the updated value is smaller than that. The gas saturation is constrained to a value of zero upon taking a negative value after data assimilation. The water and gas saturations from the updated state vector (unconstrained solution) are analyzed simultaneously for any constraint violation and the constraints are enforced together. Results from Data Assimilation Fig. 12 shows the prior and posterior water saturations at the 3rd data assimilation time step for realization #68 using the standard EnKF. Fig. 12(b) indicates that the posterior water saturation from the data assimilation step is less than the connate water saturation (0.15087) in some part of the reservoir. As the water saturation cannot go lower than the connate water saturation, the truncation approach used in the implementation of the standard EnKF method would set all such water saturations in the model to the connate water saturation. Truncation or constraints had also been used at previous assimilations steps. Although the standard truncation scheme handles the issue of constraint violation, it does so by addressing the state variables (water saturation) only in the affected region and does not guarantee the consistency between the state variables in the remaining part of the model. 50 40 50 PROD1 PROD2 40 PROD1 PROD2 0.2 0.2 0.5 0.3 0.5 0.6 0.4 0.6 30 0.2 20 30 0.3 0.8 INJ 0.7 0.7 0.8 INJ 0.1 0.1 0.1 20 0.6 0.1 0.4 10 10 PROD4 0 PROD3 10 20 30 (a) Prior 40 PROD4 50 0 PROD3 10 20 30 40 50 (b) Posterior Figure 12: The prior and posterior water saturations at the 3rd data assimilation time step for realization #68 using the standard EnKF The application of the constrained EnKF method ensures that the inequality constraints on the water and gas saturations are honored at each data assimilation time step. Fig. 13 shows the prior and posterior water saturations for realization #68 using the constrained EnKF. The posterior water saturations from the standard EnKF and the constrained EnKF are shown in Fig. 12(b) and Fig. 13(b), respectively, and they both match the production data from the reference model. The posterior water saturations from the constrained EnKF, however, also honor the 12 SPE 125101 inequality constraints. 50 40 50 PROD1 PROD2 40 PROD1 PROD2 0.3 0.5 30 20 30 0.4 0.7 0.8 INJ 0.7 0.8 INJ 0.2 0.2 0.6 20 0.2 0.6 0.3 0.5 0.4 10 10 PROD4 0 PROD3 10 20 30 (a) Prior 40 PROD4 50 0 PROD3 10 20 30 40 50 (b) Posterior Figure 13: The prior and posterior water saturations at the 3rd data assimilation time step for realization #68 using the constrained EnKF Matching Production Data The production data are continually assimilated using the standard EnKF and the constrained EnKF. Fig. 14 shows the production data match for the oil production rate for producers PROD1, PROD2, and PROD3 from the application of the two methods. The constrained EnKF performed better in terms of matching the production data compared to the standard EnKF. Even though the ensemble obtained from the standard EnKF appears to cover the observations from the reference model, the spread of the ensemble members around the measurements is greater in case of the standard EnKF than the constrained EnKF. The ensemble obtained from the constrained EnKF appears to have enough variability in the prediction phase which is important for assimilating additional production data in the future. The results of the production data match for the producing gas-oil ratio for some of the producers from the two methods are shown in Fig. 15. The constrained EnKF method generally performed better than the standard EnKF in terms of matching the observations of gas-oil ratio for producers PROD2 and PROD3. Similar results were observed for water production and injection rates. Mass conservation The most significant problem that was observed with the application of the standard EnKF to the compositional data assimilation problem was that excess mass was added during the update step at the first data assimilation time. In general, mass should not be expected to be conserved in data assimilation or history matching, as the mass of a component is typically an unknown quantity that must be estimated. On the other hand, when a substance is injected into the reservoir at a known rate, it would be desirable for that quantity to be conserved in the updating process. In this section, we will investigate experimentally the consequences of the analysis step on mass conservation in a simple model. Consider one-dimensional, constant-rate injection of water into a reservoir of length 50 that is initially at irreducible water saturation (0.2) and assume that porosity is a correlated random field with mean 0.25 and standard deviation 0.07. The covariance for porosity is nearly Gaussian with a practical range of 25. At some particular time (560 in dimensionless units), we measure the saturation at the center of the reservoir. Fig. 16(a) shows the distribution of forecast water saturation from a large number of reservoir realizations at the measurement time. If the true water saturation at x = 25 is 0.2, then the result of a standard ensemble Kalman filter analysis step is as shown in Fig. 16(b). Approximately 10% of the updated realizations contain negative saturations after updating. Fig. 16(c) shows the same saturation profiles after updating using the constrained EnKF. Constraining the saturation values to lie in the feasible region is not sufficient to ensure that the profiles are monotonic and the advantage over truncation is not obvious because in the two-phase incompressible model, the problems with volume conservation are not so significant as in the compositional model. Truncation of negative water saturations results in a corresponding reduction in the updated oil saturation, which does not occur in the compositional model. While the total volume will not be effected by truncation in this case, the volume of the individual phases will be. Fig. 17 compares the volume of moveable water at the end of the forecast with the volume of moveable water after EnKF with truncation and with the constrained EnKF. The small variability in the volume at the end of the 1000 800 600 400 200 100 200 300 400 500 1000 500 0 0 100 Time Days 600 400 200 300 400 500 Oil production rate STB day Oil production rate STB day 800 200 400 500 100 1000 500 100 Time Days 200 300 400 500 1400 1200 1000 800 600 400 200 0 0 100 Time Days (d) Constrained EnKF: PROD1 300 400 500 (c) Standard EnKF: PROD3 1500 0 0 200 Time Days (b) Standard EnKF: PROD2 1000 100 300 1400 1200 1000 800 600 400 200 0 0 Time Days (a) Standard EnKF: PROD1 0 0 200 Oil production rate STB day 0 0 1500 Oil production rate STB day 13 Oil production rate STB day Oil production rate STB day SPE 125101 200 300 400 500 Time Days (e) Constrained EnKF: PROD2 (f) Constrained EnKF: PROD3 2.4 2.2 2.0 1.8 1.6 1.4 1.2 1.0 0 100 200 300 400 500 Gas oil ratio MSCF STB 4 Gas oil ratio MSCF STB Gas oil ratio MSCF STB Figure 14: Comparison of the oil production rate (OPR) for producers PROD1, PROD2, and PROD3 obtained from rerunning the final updated ensemble from the standard EnKF and the constrained EnKF 3 2 1 0 0 100 Time Days 400 500 2.4 1.8 1.6 1.4 1.2 300 0.5 0 100 400 Time Days (d) Constrained EnKF: PROD1 500 3 2 1 0 0 100 200 300 200 300 400 500 (c) Standard EnKF: PROD3 Gas oil ratio MSCF STB 2.0 200 1.0 Time Days 4 2.2 100 1.5 (b) Standard EnKF: PROD2 Gas oil ratio MSCF STB Gas oil ratio MSCF STB 300 2.0 Time Days (a) Standard EnKF: PROD1 1.0 0 200 2.5 400 Time Days (e) Constrained EnKF: PROD2 500 2.5 2.0 1.5 1.0 0.5 0 100 200 300 400 500 Time Days (f) Constrained EnKF: PROD3 Figure 15: Comparison of the gas oil ratio (GOR) for producers PROD1, PROD2, and PROD3 obtained from rerunning the final updated ensemble from the standard EnKF and the constrained EnKF 14 SPE 125101 0.8 0.8 0.8 0.7 0.7 0.6 0.6 0.6 0.4 0.5 0.4 0.5 0.4 0.2 0.3 0.3 10 10 20 30 40 20 30 40 50 50 10 (a) Forecast (b) Standard EnKF 20 30 40 50 (c) Constrained EnKF Figure 16: Water saturation profiles forecast is a result of discretization of the reservoir and the use of grid block center saturation values. The standard EnKF with truncation shows a bias towards higher water volume because the negative values are truncated. The constrained EnKF does not remove the bias, but does reduce the problem substantially. 12.5 12.0 11.5 11.0 Forecast EnKF constrained EnKF Figure 17: Distributions of volume of moveable water in the reservoir at the end of several operations: the forecast step, the standard EnKF update with truncation of saturations, and the constrained EnKF. If mass is explicitly included in the state vector, and all forecast realizations have the same mass, then the mass in all analyzed realizations will also be identical. The mass of a component in each gridblock must of course be non-negative, but mass will not necessarily be non-negative after an EnKF update. The consequence of truncating negative values of mass is that mass will not be conserved in the analysis step. Application of constraints to the mass variable would necessarily conserve the mass of each component if all realizations initially contained the same mass. In general, however, the variables that are included in the state vectors are saturations, porosities, molar densities, or mass fractions. When saturations or molar densities are used in the state vector, mass is not necessarily conserved even when truncation is not used. When truncation is used, the results are typically much worse. Fig. 18 shows results from the same Buckley-Leverett displacement experiment whose results were shown in Fig. 16. In this case, however, the molar densities and porosities are included in the state variable instead of saturations and porosities. One consequence of this choice of parameterization is that unlike saturation, there is no simple upper limit on the molar density. Note that truncation of the low values of molar density of water (solid curve in Fig. 18(a)), does not effect the molar density of oil (solid curve in Fig. 18(b)). Also, truncation of the molar density of water has no effect on the porosity field (solid curve in Fig. 18(c)). As a result, there is 9% more total moles of both components in the realization after truncation than in the realization after the analysis step. If the constrained EnKF is used, the molar density of water is quite similar to the molar density after truncation (dashed curve in Fig. 18(a)) while the molar density of oil and the porosity field are both reduced substantially (dashed curves in Fig. 18(b) and (c)). The net effect is that the total number of moles is only 2% less than the forecast, so the mass is much more nearly conserved, even though the mass (or number of moles) was not used directly in the state vector. SPE 125101 15 45 10 0.40 8 0.35 6 0.30 4 0.25 40 35 30 25 20 10 20 30 40 50 10 (a) Molar density of water 20 30 (b) Molar density of oil 40 50 10 20 30 40 50 (c) Porosity Figure 18: Comparison of molar densities and porosity for Realization 1, with simple truncation (solid) and with constrained EnKF (dashed). Conclusions In this paper we describe a new method for sequential data assimilation that ensures that the updated state variables satisfy inequality constraints. The method is a constrained version of the ensemble Kalman filter that is efficient in that it does not require additional model evaluations for constraint application. The constrained EnKF method is demonstrated on a compositional model for which state variables have lower limits, and on a threephase black oil model for which saturations are constrained with both upper and lower limits. When the standard EnKF approach is used for history matching with a compositional reservoir model having the molar densities of the components as a part of the state vector, it resulted in problems with the updated molar densities. The use of simple truncation in this case solved the issue of non-physical (negative) molar density values for the reservoir simulator, but it did not address the issue of the increased molar densities of some of the other components. The ultimate result was a substantial pressure build-up in the reservoir and subsequent loss of injection. The application of the constrained EnKF with a constraint on the CO2 molar density reduced the molar densities of other components in a systematic manner and solved the issues of excess mass and pressure. We also showed that in cases where mass conservation of a component is important, the constrained EnKF gives superior results to the standard truncation scheme. The proposed constrained EnKF method can be easily implemented with some minor modifications to the pre-existing codes for the standard EnKF. The additional computational cost and time for using the constrained EnKF is significantly lower compared to the time required for running the reservoir simulator. Acknowledgment Support for Hemant Phale was provided by the OU Consortium on Ensemble Methods (OUCEM). The authors would like to acknowledge the donation of multiple licenses of ECLIPSE by Schlumberger and the computational resources provided by the OU Supercomputing Center for Education and Research (OSCER). References Aanonsen, S. I., G. Nævdal, D. S. Oliver, A. C. Reynolds, and B. Vall`s, Ensemble Kalman filter in reservoir e engineering — a review, SPE Journal, approved for publication, 2009. Bianco, A., A. Cominelli, L. Dovera, G. Nævdal, and B. Vall`s, History matching and production forecast e uncertainty by means of the ensemble Kalman filter: A real field application (SPE-107161), in EUROPEC/EAGE Conference and Exhibition, 2007. Chen, Y., D. S. Oliver, and D. Zhang, Data assimilation for nonlinear problems by ensemble Kalman filter with reparameterization, Journal of Petroleum Science and Engineering, 66, 1–14, 2009. De Geeter, J., H. Van Brussel, and J. De Schutter, A smoothly constrained Kalman filter, IEEE Transactions on Pattern Analysis and Machine Intelligence, 19, 1171–1177, 1997. Deutsch, C. V. and A. G. Journel, GSLIB: Geostatistical Software Library and User’s Guide, second edn., Oxford University Press, New York, 1998. Evensen, G., Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics, Journal of Geophysical Research, 99(C5), 10,143–10,162, 1994. 16 SPE 125101 Evensen, G., J. Hove, H. C. Meisingset, E. Reiso, K. S. Seim, and Ø. Espelid, Using the EnKF for assisted history matching of a North Sea reservoir model (SPE 106184), in Proceedings of the 2007 SPE Reservoir Simulation Symposium,, 2007. Gu, Y. and D. S. Oliver, The ensemble Kalman filter for continuous updating of reservoir simulation models, Journal of Energy Resources Technology, 128(1), 79–87, 2006. Gu, Y. and D. S. Oliver, An iterative ensemble Kalman filter for multiphase fluid flow data assimilation, SPE Journal, 12(4), 438–446, 2007. Gupta, N. and R. Hauser, Kalman filtering with equality and inequality state constraints, Tech. Rep. Report No. 07/18, Numerical Analysis Group, Oxford Univ. Computing Laboratory, Univ. Oxford, Oxford, UK, 2007. Haugen, V., L.-J. Natvik, G. Evensen, A. Berg, K. Flornes, and G. Nævdal, History matching using the ensemble Kalman filter on a North Sea field case (SPE-102430), in SPE Annual Technical Conference and Exhibition, 2006. Ko, S. and R. R. Bitmead, State estimation for linear systems with state equality constraints, Automatica, 43(8), 1363–1368, 2007. Li, G. and A. C. Reynolds, Iterative ensemble Kalman filters for data assimilation, SPE Journal, approved for publication, 2009. Li, R., A. C. Reynolds, and D. S. Oliver, History matching of three-phase flow production data, SPE Journal, 8(4), 328–340, 2003. Lorentzen, R. J., K. K. Fjelde, J. Frøyen, A. C. V. M. Lage, G. Nævdal, and E. H. Vefring, Underbalanced drilling: Real time data interpretation and decision support, in SPE/IADC Drilling Conference, 2001. Massicotte, D., R. Z. Morawski, and A. Barwicz, Incorporation of a positivity constraint into a Kalman filterbased algorithm for correction of spectrometric data, IEEE Transactions on Instrumentation and Measurement, 44(1), 2–7, 1995. Nævdal, G., L. M. Johnsen, S. I. Aanonsen, and E. H. Vefring, Reservoir monitoring and continuous model updating using ensemble Kalman filter, SPE Journal, 10(1), 66–74, 2005. Nævdal, G., T. Mannseth, and E. H. Vefring, Near-well reservoir monitoring through ensemble Kalman filter: SPE 75235, in Proceedings of SPE/DOE Improved Oil Recovery Symposium, 2002. Oliver, D. S., A. C. Reynolds, Z. Bi, and Y. Abacioglu, Integration of production data into reservoir models, Petroleum Geoscience, 7(SUPP), 65–73, 2001. Qureshi, A. G., M. P. Res, and B. C. Burnaby, Constrained Kalman filtering for image restoration, in Acoustics, Speech, and Signal Processing, 1989. ICASSP-89., 1989 International Conference on, pp. 1405–1408, 1989. Rao, C. V., J. B. Rawlings, and D. Q. Mayne, Constrained state estimation for nonlinear discrete-time systems: Stability and moving horizon approximations, IEEE Transactions on Automatic Control, 48(2), 246–258, 2003. Rotea, M. and C. Lana, State estimation with probability constraints, International Journal of Control, 81(6), 920–930, 2008. Schlumberger, Eclipse Technical and Reference Manual, 2007. Simon, D. and T. L. Chia, Kalman filtering with state equality constraints, IEEE Transactions on Aerospace and Electronic Systems, 38(1), 128–136, 2002. Simon, E. and L. Bertino, Application of the Gaussian anamorphosis to assimilation in a 3-D coupled physicalecosystem model of the North Atlantic with the EnKF: a twin experiment, Ocean Science Discussions, 6(1), 617–652, 2009. Skjervheim, J.-A., G. Evensen, S. I. Aanonsen, B. O. Ruud, and T. A. Johansen, Incorporating 4D seismic data in reservoir simulation models using ensemble Kalman filter, SPE Journal, 12(3), 282–292, 2007. Tarantola, A., Inverse Problem Theory: Methods for Data Fitting and Model Parameter Estimation, Elsevier, Amsterdam, The Netherlands, 1987. SPE 125101 17 Thacker, W. C., Data assimilation with inequality constraints, Ocean Modelling, 16(3-4), 264–276, 2007. Thacker, W. C. and R. B. Long, Fitting dynamics to data, Journal of Geophysical Research, 93(1), 1227–1240, 1988. Ungarala, S., E. Dolence, and K. Li, Constrained extended Kalman filter for nonlinear state estimation, in Proccedings of the 8th International IFAC Symposium on Dynamics and Control of Process Systems, 2007. Zhang, Y. and D. S. Oliver, History matching using a hierarchical stochastic model with the ensemble Kalman filter: A field case study, SPE-118879, in Proceedings of the 2009 SPE Reservoir Simulation Symposium, The Woodlands, February 2–4, 2009. ...
View Full Document

Ask a homework question - tutors are online