- Email: [email protected]

Contents lists available at ScienceDirect

Advances in Water Resources journal homepage: www.elsevier.com/locate/advwatres

Flow relevant covariance localization during dynamic data assimilation using EnKF Deepak Devegowda a,*, Elkin Arroyo-Negrete b, Akhil Datta-Gupta c,1 a

University of Oklahoma, Mewbourne School of Petroleum and Geological Engineering, T-301, 100 E. Boyd Street, Norman, 73019 OK, USA Oxy Oil and Gas Corp., Houston, TX, USA c Texas A&M University, Department of Petroleum Engineering, College Station, TX 77843-3116, USA b

a r t i c l e

i n f o

Article history: Received 20 April 2009 Received in revised form 29 September 2009 Accepted 2 October 2009 Available online 9 October 2009 Keywords: Sequential data assimilation Ensemble Kalman ﬁlter Covariance localization Subsurface characterization Streamlines

a b s t r a c t Multiphase dynamic data integration into high resolution subsurface models is an integral aspect of reservoir and groundwater management strategies and uncertainty assessment. Over the past two decades, advances in computing and the development and implementation of robust algorithms for automatic history matching have considerably reduced the time and effort associated with subsurface characterization and reduced the subjectivity associated with manual model calibration. However, reliable and accurate subsurface characterization continues to be challenging due to the large number of model unknowns to be estimated using a relatively smaller set of measurements. For ensemble-based methods in particular, the difﬁculties are compounded by the need for a large number of model replicates to estimate sample-based statistical measures, speciﬁcally the covariances and cross-covariances that directly impact the spread of information from the measurement locations to the model parameters. Statistical noise resulting from modest ensemble sizes can overwhelm and degrade the model updates leading to geologically inconsistent subsurface models. In this work we propose to address the difﬁculties in the implementation of the ensemble Kalman ﬁlter (EnKF) for operational data integration problems. The methods described here use streamline-derived information to identify regions within the reservoir that will have a maximum impact on the dynamic response. This is achieved through spatial localization of the sample-based cross-covariance estimates between the measurements and the model unknowns using streamline trajectories. We illustrate the approach with a synthetic example and a large ﬁeld-study that demonstrate the difﬁculties with the traditional EnKF implementation. In both the numerical experiments, it is shown that these challenges are addressed using ﬂow relevant conditioning of the cross-covariance matrix. By mitigating sampling error in the cross-covariance estimates, the proposed approach provides signiﬁcant computational savings through the use of modest ensemble sizes, and consequently offers the opportunity for use with large ﬁeld-scale groundwater and reservoir characterization studies. Ó 2009 Elsevier Ltd. All rights reserved.

1. Introduction The adverse impact of inaccurate forecasts in groundwater hydrology and in petroleum reservoir applications can be substantial. Thus, proper characterization of the subsurface heterogeneity and quantiﬁcation of uncertainty are critical aspects of subsurface resource management strategies. This is achieved through the use of model calibration algorithms that integrate observations to improve the estimates of the model unknowns, thereby improving the predictive capability of the resulting model realizations. The recent trend in acquiring real-time performance data using down-hole well monitors and permanent sensors has driven the need for recursive and continuous model updating. This has led

* Corresponding author. Tel.: +1 405 325 3081; fax: +1 405 325 7477. E-mail addresses: [email protected] (D. Devegowda), Elkin_ArroyoNe [email protected] (E. Arroyo-Negrete), [email protected] (A. Datta-Gupta). 1 Tel.: +1 979 847 9030. 0309-1708/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.advwatres.2009.10.001

to the increased interest in the ensemble Kalman ﬁlter (EnKF) as a method for data assimilation [15,22]. In the past few years, the EnKF has been implemented for dynamic data assimilation and parameter estimation in a diverse set of disciplines [2,4,15,19,22,26]. For petroleum reservoir engineering, the primary focus has been to update geologic models via history matching or inverse modeling of production response. Starting with an ensemble of model realizations (typically permeability distributions) conditioned to all previously obtained data, the updating process utilizes the ensemble-derived statistics to calibrate the model parameters whenever measurements become available. Consequently, it is possible to maintain an ensemble of model realizations that are ‘current’ with ﬁeld observations and can be used for reservoir performance forecasting. The ease of implementation, the ability for uncertainty quantiﬁcation and the versatility in handling diverse data types make the EnKF a very appealing approach for data assimilation problems.

130

D. Devegowda et al. / Advances in Water Resources 33 (2010) 129–145

However, for ﬁeld-scale subsurface characterization problems, the EnKF can be computationally challenging because it requires ﬂow simulation for every ensemble member. In particular, for multiphase ﬂow problems and high resolution geologic models, these ﬂow simulations can be computationally expensive. A crucial aspect of the EnKF update is the ensemble-based estimation of various covariances and cross-covariances that relate the data to model parameters. These estimates can be severely compromised when using small ensemble sizes because of sampling errors [4,17,27]. The difﬁculties are exacerbated when we need to estimate a large number of model unknowns. The underlying non-linearity of the problem, for e.g. in multiphase ﬂow applications, can further degrade the EnKF update and lead to geologically inconsistent changes to the models [5,12,19]. A common observation with the use of the EnKF is the lack of geological and physical consistency between the pre- and post-assimilation model realizations. In addition to these challenges, the loss of ensemble variability can severely underestimate the underlying uncertainties [12,13]. The loss of variability causes the EnKF to ignore newer observations because the uncertainty in the prior models as reﬂected by the sample-derived covariances becomes negligible [13,20]. These difﬁculties can be somewhat mitigated by using larger ensemble sizes and hence, more reliable representations of the ensemble-derived statistics. However, this choice can be computationally infeasible for large problems. The use of improved variants of the EnKF for practical data assimilation is common throughout the published literature [2,6,23,25,36]. These methods rely on addressing either the issue of the loss of variance or the sampling related difﬁculties or both. Covariance inﬂation schemes [2,3] target the variance deﬁciency of the ensemble observed following a sequence of updates. On the other hand, covariance localization schemes [20] utilize prior or ad-hoc knowledge about the nature of the parameter and data correlations to obtain improved estimates of the various covariance matrices. Distance-based covariance localization assumes that correlations between variables at different locations become insigniﬁcant beyond a certain cut-off distance [17,20]. Within the selected zone centered at the observation location and extending out to the cut-off radius, the ensemble-derived correlations are assumed to decrease with distance. The application of distance-based localization, however, may not be effective for 3-dimensional and highly heterogeneous reservoir and groundwater models [5]. To adequately capture the relationship between the data and the model unknowns, the localizing function should be closely related to the underlying ﬂow ﬁeld. In this paper, we present two novel covariance localization schemes using streamline-derived information. The streamlines are readily obtained from the underlying ﬂow ﬁeld by integrating the ﬂuid ﬂuxes [9,10,24,30]. The ﬁrst method simply utilizes streamline trajectories to retain relevant spatial correlations and eliminate spurious terms in the cross-covariance estimates. The second approach utilizes the same underlying principles, but goes a step further by weighting the streamline trajectory based estimates using a localizing function based on the parameter sensitivities. The streamline trajectories and sensitivities can be computed easily and efﬁciently from the underlying velocity ﬁeld and adds very little computational overhead [10,24,32,33]. Because the streamline-based localizing schemes are tied to the underlying geology and the accompanying ﬂow ﬁeld, the primary beneﬁts of our approach are improved geologic realism and reduced loss of variance during assimilation. It also greatly enhances the computational efﬁciency because the streamline-based localization allows for the use of a smaller ensemble size compared to the traditional ensemble Kalman ﬁlter. In the next section, we provide a brief introduction to the EnKF for groundwater and petroleum engineering applications. We fol-

low the discussion with a section devoted to data assimilation for a synthetic test example and highlight the difﬁculties with the EnKF for non-linear multiphase problems and non-Gaussian prior models. Next, we examine the distance-based localization and address its limitations for subsurface characterization. Finally, we discuss the two proposed ﬂow relevant localization approaches and illustrate their beneﬁts using a synthetic example and a ﬁeldstudy. 2. EnKF: background There have been a variety of applications of the EnKF for a diverse set of problems since the original implementation by Evensen in 2003. Consequently, this section will provide only the mathematical details of the EnKF in relation to subsurface modeling. The basic derivation remains the same and the optimal estimates of the model unknowns can be derived from the Bayesian formalism or from the standpoint of minimizing the posterior uncertainty of the updated model variables [1,15,16,19]. It should be noted that the EnKF is optimal only for Gaussian model statistics and linear model dynamics. In groundwater and petroleum related applications, the augmented EnKF state vector, y at time, k, is deﬁned by a vector of the model unknowns (state and parameters) and the observed data. The model unknowns, comprising, for example, the grid block permeabilities, porosities, phase saturations and pressures, are denoted by the vector, m of length Nm and the data is denoted by the vector, d of size Nd. Each model realization within the ensemble of size, Ne, is represented by the associated state vector. The goal of the EnKF is to modify the model state and parameters so that the individual ensemble members honor the observations, dobs. Additionally, the observed data is assumed to have a Gaussian noise with zero mean and variance, R. Rather than assimilating the data, we now assimilate ‘realizations’ of the data with the associated noise. This implies that any sensor bias is removed from the observations in a pre-processing step before they are available for assimilation

dk;j ¼ dobs;k þ ej ¼ Hytrue;k;j þ ej ¼ H

mk;j dobs;k

þ ej ;

j ¼ 1; N e

ð1Þ

The trivial matrix H that relates the true, but unknown, state vector, ytrue, to the observed data is expressed as

H ¼ ½0 I

ð2Þ

I is the identity matrix of dimensions Nd Nd and 0 has the dimensions Nm Nd, so the action of H on the state vector, y is to simply select the data component, d. The initial suite of model realizations (permeability distributions in our case) that constitutes the ensemble is constructed using geostatistical interpolation methods [11] that utilize core measurements, well log-derived rock properties and even seismic data in conjunction with apriori geologic information. These data types are classiﬁed as static data to reﬂect the time-invariant nature of the information. The spatial relationships between the model parameters are assumed to be known through a covariance model-derived from the statistics of the data. In a probabilistic sense, the initial samples can be said to be drawn from a prior probability distribution that reﬂects our uncertainty in the initial estimates of the spatially varying properties. Starting with an ensemble of subsurface models that already incorporates prior information, the EnKF recursively integrates additional data when it becomes available. At assimilation step, k, if a new measurement is obtained, the update equation in the EnKF for each realization within the ensemble can be expressed as

D. Devegowda et al. / Advances in Water Resources 33 (2010) 129–145 apriori ypost þ Kðdk;j Hyapriori Þ; k;j ¼ yk;j k;j

j ¼ 1; N e

ð3Þ

where K is known as the Kalman gain and directly relates the data mismatch between model predictions and observations to the required changes in the model unknowns. The updated models now represent our best knowledge of the state of the system at assimilation step, k, given all previously obtained data. A numerical ﬂow and transport simulation then provides the dynamic responses and the forecast state at the next assimilation step, k + 1. Successive steps of prediction and model corrections consequently lead to the ensemble of models that are consistent with observed dynamic data and prior static information. The best estimate of the current state should be optimal in some pre-deﬁned sense, either as a minimum variance estimate or as the estimate that maximizes the posterior probability [1,16]. Therefore, the Kalman gain can be expressed as,

K ¼ Cy HT ðHCy HT þ RÞ1

ð4Þ

In the above equation, Cy is the state covariance matrix and by virtue of our deﬁnition of the augmented state vector, y, we have

Cy ¼

CM

CM;d

Cd;M

Cd;d

ð5Þ

In Eq. (4), R is a covariance matrix that reﬂects our uncertainty in the measurements. The model covariance, CM provides the spatial correlation between the model parameters. The physical relationships between the data and the model parameters are embedded in the cross-covariance estimate, CM,d and the ensemble-derived product CyHT in the expression for the Kalman gain can be seen to be dependent primarily on the cross-covariance. The model updates via the EnKF, in effect, are therefore strongly dependent on the accuracy of the ensemble-derived cross-covariance calculations. Sampling related error arising from small ensemble sizes can signiﬁcantly degrade EnKF performance. In addition, model variability, represented by Cy, is a key component of the update expression. In particular, when there is a signiﬁcant loss in variance during data assimilation (Cy 0), the Kalman gain becomes negligible and therefore, the EnKF completely disregards future observations. Both the loss of variance and the degraded ﬁlter performance are closely tied to the accuracy of the ensemblebased cross-covariance calculations. In the next section, we illustrate these effects using a test case and examine the role of covariance localization in mitigating these effects.

131

2.1. An illustration: parameter estimation in multiphase ﬂow In this example we illustrate the application of the EnKF for integrating dynamic data during multiphase ﬂow in an oil reservoir undergoing water injection. We also highlight some of the problems encountered by EnKF for non-linear problems and nonGaussian priors. The reservoir is deﬁned on a 50 50 grid with four producers located at the corners and a water injection well in the center. The observations are generated from a reference model shown in Fig. 1 and comprise of water-cut measurements at each of the four producing wells. The water-cut is simply the ratio of the water ﬂow rate to the total (oil and water) ﬂow rate. The wells are constrained to produce a certain volume of total ﬂuid per day and the injector is constrained to inject a pre-speciﬁed volume of water. An initial ensemble of model realizations (permeabilities) was generated using the Sequential Gaussian Simulation [11] with an assumed prior covariance model. Also, the individual realizations are chosen to have a bi-modal log-permeability histogram in order to highlight some of the difﬁculties encountered by the EnKF with non-Gaussian priors. We chose to assimilate data for 3000 days and forecast for an additional 1000 days. The plot in Fig. 2 shows the initial ensemble predictions for water-cut and the corresponding matches to the observed history after assimilation for an ensemble size of 50. As expected, a signiﬁcant reduction in the spread of the forecast is seen because of the water-cut data assimilation. Next, we performed the same experiment using different ensemble sizes, starting from an ensemble size of 25, then 50 and ﬁnally with a 100-member ensemble. In all the cases studied, the EnKF was able to satisfactorily reproduce the production history. However, a closer examination of the updated permeability ﬁelds reveals some of the shortcomings of the EnKF for this nonlinear and non-Gaussian example. In Fig. 3, we have shown four arbitrarily chosen realizations from the ensemble of 25 models for the pre- and post-assimilation steps. The apparent resemblance of each of the updated realizations is a consequence of the dramatic loss of variability during assimilation, partly caused by the small ensemble size. Although our conﬁdence in the posterior estimates would improve through data fusion, it is unrealistic to expect complete resolution of all the subsurface unknowns; ﬁrstly, because of the inherent low resolution of the production data and secondly, because the number of model unknowns is signiﬁcantly larger than the amount of available data. Additionally, artifacts such as localized regions with extreme values of permeability in the updated models indicate a signiﬁcant departure from prior spatial distribution resulting in a

Fig. 1. The reference permeability for the synthetic example.

132

D. Devegowda et al. / Advances in Water Resources 33 (2010) 129–145

Fig. 2. The matches to the water-cut measurements (in red) at each of the three producing wells in the bottom row. The initial model predictions are in the top row. (For interpretation of the references to colour in this ﬁgure legend, the reader is referred to the web version of this article.)

Fig. 3. For an ensemble size of 25, this ﬁgure shows a comparison of four randomly selected realizations updated with water-cut measurements via the EnKF with the corresponding initial permeability ﬁelds.

loss of geologic realism. The log-permeability histogram from one randomly selected realization is shown in Fig. 4. Clearly it shows a signiﬁcant departure from the prior histogram and the relative proportion of rock facies speciﬁed in the initial model. The lack of geologic consistency can adversely impact the predictive capability of the resulting model realizations. This is especially true for a variety of applications, e.g. when trying to identify suitable locations for additional wells, for the design of enhanced oil recovery processes or for inferring the predominant pathways for tracer and contaminant transport. The sensitivity to sample size is best illustrated with Figs. 5 and 6. It is evident that some of the difﬁculties with the EnKF are some-

what mitigated as the ensemble size is increased [12,13,20]. As expected, the estimates based on sample statistics become more reliable as the sample size is increased. The inaccuracy of the cross-covariance estimates adversely affects the performance of the EnKF when using a small number of model replicates. This conclusion is reinforced by Fig. 7 which plots the total variance contained in the ensemble as a function of assimilation time. However, our ability to increase the ensemble size is limited by the computational requirements which can become prohibitive for large-scale ﬁeld studies. For practical applications, it is important to keep the ensemble size relatively small and at the same time improve the cross-covariance estimates. This has been the

D. Devegowda et al. / Advances in Water Resources 33 (2010) 129–145

133

Fig. 4. (a) The log-permeability histogram for one arbitrarily chosen realization. (b) The Gaussian characteristics of the corresponding histogram post-assimilation.

Fig. 5. For an ensemble size of 100, this ﬁgure shows a comparison of four randomly selected realizations updated with water-cut measurements via the EnKF with the corresponding initial permeability ﬁelds.

motivation behind the concept of spatially localizing the effects of measurements, otherwise known as covariance localization [5,12,13,20,27]. 3. Covariance localization for ensemble-based data assimilation The goal of covariance localization for ensemble-based estimation problems is to eliminate or dampen the effects of spurious correlations between the model parameter, state and the measurements, typically arising from the limited ensemble size. Mathematically, this can be accomplished using the Schur product (element-by-element multiplication) of a localizing function, q, with the sample-based Kalman gain estimate [20,34],

ypost ¼ yapriori þ q Kðdobs Hyapriori Þ k k k

ð6Þ

The resulting expression for the Kalman gain can be written as,

KLocalized ¼ q K ¼ q Cy HT ðHCy HT þ RÞ1

ð7Þ

The justiﬁcation for a particular form of the localizing function will be dictated by the type of application and to some extent by the ease of implementation. The following section discusses the effectiveness of covariance localization for data assimilation using ensemble methods. Speciﬁcally, we compare the commonly used distance-based localization scheme with two proposed ﬂow relevant localizations schemes based on the physics of multiphase ﬂow. 3.1. The need for covariance localization: an illustration To further motivate the need for covariance localization, let us consider an ensemble of two-dimensional reservoir models with spatially uncorrelated permeability distribution. Essentially, these represent samples from a multi-Gaussian distribution with a ﬁxed mean and negligible spatial correlation length. As before, we have four producing wells at each of the corners and one injection well at the center of the ﬁeld. At an arbitrarily chosen assimilation step, the cross-covariance between the permeability values and the

134

D. Devegowda et al. / Advances in Water Resources 33 (2010) 129–145

Fig. 6. Sensitivity to ensemble size for the synthetic example viewed in terms of the log-permeability histograms. In spite of increasing ensemble size to a 100 members, we do not maintain prior model statistics.

Fig. 7. Total variance versus time during assimilation. Increasing the ensemble size reduces the loss of variance.

based on a 50-member ensemble are severely degraded by the small sample size. Spurious correlations at distant grid locations dominate the structure of the estimates. Additionally, the covariance structure compared to the 1000-member ensemble is signiﬁcantly corrupted by noise. The inaccuracy in the cross-covariance computations will adversely impact parameter estimates leading to geologically inconsistent posterior models as seen in our previous examples. Previous studies [19,16] have indicated that for bounded error growth in the model covariance estimation, the ensemble size should vary as the square of the number of variables to be estimated. This is clearly an impractical proposition for largescale systems. This underscores the need for a robust and efﬁcient covariance localization method that minimizes spurious correlations arising from modest ensemble sizes dictated by practical applications. In the sections below, we ﬁrst discuss a previously proposed distance-based covariance localization scheme [17,20,34] and some of its difﬁculties and limitations. We then propose two streamline-based localization schemes and demonstrate their effectiveness using synthetic and ﬁeld examples. 3.2. Distance-dependent covariance localization

water-cut response can be estimated from the ensemble and the results are shown in Fig. 8 for one of the producing wells (top left corner). We compare the estimates from a 50-member ensemble and a much larger 1000-member ensemble. This particular example with no spatial correlation was deliberately chosen to highlight the effects of ensemble size on cross-covariance estimates. In fact, the ensemble-based estimate is expected to be similar to that of a homogeneous permeability ﬁeld and consequently, will have prominent highs along the diagonal leading from the injection well to the selected producer with reduced values elsewhere. The cross-covariance estimates from the 1000-member ensemble which can be considered closer to the truth, shows the expected behavior as shown in Fig. 8. However, the calculations

Intuitively, the inﬂuence of a measurement on the state and parameter estimation is expected to decrease with increasing distance. This assumption forms the basis for the distance-dependent covariance localization. The impact of a measurement is assumed to be negligible beyond a pre-speciﬁed radius of inﬂuence [20]. Within this region of inﬂuence, typically a monotonically decreasing correlation function centered at the observation is used for covariance localization [18]. An example of a localizing function based on this concept is shown in Fig. 9. In this section, we examine the application of the distancebased localization to the waterﬂood example introduced in Section 2.1 with the reference permeability shown in Fig. 1. An ensemble

D. Devegowda et al. / Advances in Water Resources 33 (2010) 129–145

135

Fig. 8. Estimates of cross-covariance between water-cut data and permeability at a certain assimilation step for differing ensemble sizes.

localization is additionally, a function of the ensemble size. For smaller sample sizes, even short-range correlations can be signiﬁcantly degraded by spurious and random noise and the correlations may need to be damped at a shorter distance. Furthermore, the simple notion of a distance-dependent localization may not be consistent with the subsurface heterogeneity that can be anisotropic and might include high conductivity channels that contribute to long-range correlations in speciﬁc directions [5]. For 3dimensional applications, these problems can be compounded by the requirements of different localizing distance in the horizontal and vertical directions because of the differences in the correlation lengths. In the light of these drawbacks, distance-based localization generally have limited applicability, particularly for subsurface characterization using EnKF. There is a need for more reliable localization techniques that are tied to the underlying geology and represent the ﬂow ﬁeld more realistically. 3.3. Streamline trajectory based covariance localization

Fig. 9. The distance-dependent localizing function for the well (black dot) showing decreasing levels of inﬂuence with increasing distance from the well.

size of 50 is chosen for this example. An initial cut-off radius that bounds the region of inﬂuence for each well is assumed to be 200 feet. For comparison, the ﬁeld dimensions are 500 feet on each side of the square-shaped reservoir. The measurements are the water-cut information at each of the three production wells. The water-cut data for the updated models are shown along with the initial ensemble predictions in Fig. 10. Although, there is a considerable reduction in the spread of the prediction envelope, there are still signiﬁcant differences between the observations and modelderived water-cut performance data. The reason for a poor match is related to the inappropriate and arbitrary choice of a cut-off radius. This underscores the need for a proper choice of a radius of inﬂuence for localizing the effects of measurements. For this example case, a cut-off distance of 400 feet seems to be more appropriate as seen in Fig. 11. In general, the choice of a suitable characteristic length-scale associated with the observations is subjective and requires some trial-and-error. Although there are some guidelines in the literatures [3], the optimal radius of inﬂuence is likely to be application dependent and may be difﬁcult to determine a priori. Hamill and Whitaker [20] observed that the length-scale chosen to achieve

In this section we introduce the streamline trajectory based covariance localization. The primary goal here is to eliminate erroneous ﬂuctuations in the cross-covariance matrix because of sampling error, thereby allowing for parameter updates that are consistent with the underlying geology and the ﬂow ﬁeld. The use of streamline-based methods is physically intuitive and the streamline time of ﬂight quantitatively deﬁnes the region of inﬂuence for measurements accounting for the interactions between the heterogeneity and the ﬂow ﬁeld. Adjusting the cross-covariance by the use of streamlines moreover keeps the changes targeted to areas that have the greatest impact on the ﬂuid ﬂow and transport. Conceptually, this is analogous to streamline-assisted history matching which has been successfully used to identify regions within the reservoir that have signiﬁcant impact on the observations [14]. The streamline-based approach does not make any a priori assumptions about the strength of correlations at different locations and thus, naturally retains relevant short and long-range correlations. The procedure is simple and computationally inexpensive to implement. Fundamentally, the Kalman update expression remains the same as shown in Eq. (1); however, the localizing function now becomes a function of the underlying ﬂow ﬁeld. At any given assimilation step, to effectively localize regions corresponding to each observation, we establish whether or not a grid cell within the region is intersected by a streamline originating from the source of the observation. If the grid cell is intersected

136

D. Devegowda et al. / Advances in Water Resources 33 (2010) 129–145

Fig. 10. Using a cut-off radius of 200 feet, the matches to the observations are not satisfactory. This results from the use of a very short radius of inﬂuence which eliminates relevant correlations.

Fig. 11. By increasing the cut-off radius to 400 feet, the matches to the observed water-cut data is signiﬁcantly better.

by the relevant streamline, it is assigned a value of 1, otherwise 0. These results form the column corresponding to the observation in the matrix, q, where a value of one refers to a location relevant to the observation. By successively accounting for each observation, we populate the localizing function, q. Within this framework, noisy terms arising due to sampling errors are eliminated if they

are not compatible with the underlying ﬂow ﬁeld. This is repeated for each model replicate within the ensemble. Consequently, this results in an ensemble of regions of inﬂuence that represents potential target zones for each realization. To implement this form of covariance localization to the entire ensemble, an aggregate zone common to the ensemble is selected

D. Devegowda et al. / Advances in Water Resources 33 (2010) 129–145

137

Fig. 12. The spatial localization using streamline trajectories shows a complicated but ﬂow relevant pattern that cannot be replicated by a notion of distance. On the left (a and b) are the streamline patterns at a speciﬁed assimilation step and the corresponding grid blocks for two randomly chosen realizations. On the right in ‘c’ is the region selected by the streamline trajectories over all the realizations which is referred to as the ‘region of inﬂuence’ in the text.

by stacking the regions of inﬂuence of each realization. Potentially, additional constraints may be applied to restrict the targeted zone of inﬂuence. For example in multiphase ﬂow, the region may be restricted to areas behind a ﬂood front. For recursive data assimilation, this implies that the target area is redrawn at every assimilation step to account for changing subsurface ﬂow conditions, additional wells or simply to reﬂect modiﬁed well conﬁgurations. A visual representation of the spatial localization region for one randomly selected well in the bottom right of a synthetic reservoir is shown in Fig. 12 at a certain assimilation step. The grid blocks selected by the streamlines highlight the complicated ﬂow pattern that is a function of the underlying heterogeneities and the well locations. Clearly, the information provided by the streamline trajectories cannot be replicated by a Naïve notion of distance. Some segments of the reservoir adjacent to the wells are not covered by streamlines and on the other hand, distant locations are included in the area of inﬂuence as dictated by the ﬂow pattern. This implies

that the EnKF-predicted changes will be targeted and restricted to areas that are strongly correlated with the observations. Mathematically, the localizing function q now consists of a matrix with columns of 0’s and 1’s designed to select the appropriate grid cells relevant to the EnKF update. To examine the beneﬁts of this approach, streamline trajectory based localization is implemented for this synthetic example. Fig. 13 shows four randomly selected permeability realizations that are calibrated with water-cut data using EnKF and streamline-based localization. Although an ensemble size of 50 was used here, the results seem to outperform the traditional EnKF with 100 members as shown in Figs. 5 and 6. The updated permeability realizations show no overshooting or undershooting and appear to preserve all the important features of the prior geologic model including the bi-modality of the permeability histogram (Fig. 14). This is a consequence of restricting the model adjustments to areas relevant to the ﬂow through the use of streamline trajectories. Another important beneﬁt of ﬂow relevant covariance localization is

Fig. 13. Four selected realizations resulting from the use of streamline trajectory based covariance localization with the EnKF in comparison with the prior models.

138

D. Devegowda et al. / Advances in Water Resources 33 (2010) 129–145

Fig. 14. The prior model statistics in (a) are better preserved with the use of streamline trajectory based covariance localization. (b) Plots the corresponding log-permeability histogram for one selected realization.

to prevent the loss of ensemble variability and thus, avoiding ﬁlter divergence. This is shown in Fig. 15 and compared to the standard EnKF, a signiﬁcantly improved result is obtained with only half the number of ensemble members. Implementation of this form of covariance localization using streamline trajectories entails very little additional computational burden. The streamlines are traced using the ﬂuid ﬂuxes that are readily available from the solution of the ﬂow problem numerically [10,24,34]. Consequently, streamline trajectory based spatial localization forms the basis of an efﬁcient means to enhance state and parameter estimation using the EnKF.

Sensitivities deﬁne the relationship between the model parameters and the response and are routinely employed for dynamic data integration using inverse modeling algorithms to update the spatial distribution of porosity and permeability based on the ﬁeld response [7,8,10,28,29,31,32,35]. In this section, we ﬁrst review the basic methodology behind the streamline-derived sensitivities for covariance localization in the EnKF. We illustrate the technique with a simple example to demonstrate improved estimates of the cross-covariance. We then provide the mathematical details behind the sensitivity computations and the application of the methodology in conjunction with the EnKF for the purposes of covariance localization.

3.4. Streamline-derived sensitivity-based covariance localization Distance-dependent localization attempts to taper the crosscovariance on the assumption that the correlations are a function of the spatial distance only. With streamline trajectories, the cross-covariance calculations become more accurate because spurious terms unrelated to the physical solution are eliminated. In using streamline-derived sensitivity-based covariance localization, we attempt to achieve the goal of tapering the cross-covariance calculations within the area of inﬂuence deﬁned by streamline trajectories by the use of a localizing function derived from the parameter sensitivities [12,13]. The sensitivities are computed efﬁciently from the streamlines.

Fig. 15. Plot of the total variance versus time during assimilation. For the same performance, the standard EnKF implementation requires a much larger ensemble in comparison with trajectory based localization of the cross-covariance.

3.4.1. Sensitivity-based localization: an illustration To illustrate the effects of ﬁnite sample size on the cross-covariance estimation and the beneﬁts of streamline-derived sensitivitybased covariance localization, we consider the example introduced in Section 3.1. For two arbitrarily chosen permeability realizations, the streamline patterns for two wells at a selected step during assimilation are shown in Fig. 16. Barring minor deviations because of the random nature of the underlying permeability values, the streamline patterns are similar to that of a homogeneous permeability ﬁeld. This is, of course, the consequence of the uncorrelated nature of the permeability distribution in this case. The cross-covariance calculations between water-cut data at one of the wells and the permeabilities for 50 and 1000 members discussed in Section 3.1 are reproduced in Fig. 17. Also shown are the localized cross-covariance calculations using streamline-based sensitivities. The 1000-member ensemble has a sample size roughly on the order of the model parameter dimensions (= 2500) and consequently, will be closely related to the true value of the cross-covariance. Additionally, the relationship between the well data and the permeabilities is seen to have the expected structure with the highs located along the diagonal and the lows further away from the diagonal. This is physically intuitive but is not apparent for the 50-member ensemble. However, spatially localizing the cross-covariance based on streamline-derived sensitivities described in the next section for the 50-member ensemble results in the cross-covariance structure very similar to the estimate from the much larger 1000 member ensemble. We now discuss the mathematical background behind the sensitivity-based localization. 3.4.2. Mathematical background and formulation Our approach relies on inexpensive streamline-based sensitivity calculations that can be used to modify the ensemble-derived cross-covariance estimates. If we denote the sensitivity matrix by

D. Devegowda et al. / Advances in Water Resources 33 (2010) 129–145

139

Fig. 16. Two randomly selected realizations with the associated streamline trajectories at two wells clearly showing the resemblance to the patterns obtained from homogenous parameter ﬁelds.

Fig. 17. Sensitivity-based localization offers an order reduction in computational effort. With just 50 members, the cross-covariance estimates resemble the unconditioned estimate from a 1000-member ensemble. Moreover, the unconditioned estimate is largely noise with smaller ensemble sizes.

G and the corresponding model covariance matrix by CM, the crosscovariance can be represented as the product

CM;d ¼ CM GT

ð8Þ

From Fig. 17, we can see that for the case with ensemble size of 50, the cross-covariance estimate is dominated by spurious correlations and a simple screening of the sample-derived correlations may not be sufﬁcient. In fact, it might be necessary to include a weighting factor to account for the relative signiﬁcance of various

spatial locations to the production response. This can be efﬁciently achieved with the use of parameter sensitivities. To impose the structure of the sensitivity matrix G on the crosscovariance, we simply multiply the sample-based correlations with the normalized sensitivities. This is illustrated in Fig. 18 which shows the localizing function for the grid cells selected by streamlines. The localizing function, q, modiﬁes the ensemble-derived cross-covariance using a Schur product in the following manner,

140

D. Devegowda et al. / Advances in Water Resources 33 (2010) 129–145

Fig. 18. The locations selected by the streamlines in (a) enter into the cross-covariance estimate and the weighting function shown in (b) tapers the resulting crosscovariances. This is the basis of streamline-derived sensitivity-based covariance localization.

Fig. 19. The changes to the models are kept targeted and minimal and the matches to the observations are reasonable.

CM;d ¼ q CM;d;Ensemble ¼ GTN CM;d;Ensemble ¼ CM;d;Ensemble GTN

ð9Þ

The last equality is a consequence of the interchangeability of the Schur product operator which simply is an element-by-element multiplication and GN refers to the normalized streamlinederived sensitivities. The elements of the localization matrix q are obtained by summation of the sensitivity coefﬁcient at each

grid location obtained from each realization within the ensemble and arranging the resulting values in column vectors corresponding to each measurement. The values within each column vector are now normalized on a scale of 0–1 using the maximum absolute sensitivity within each column. This procedure essentially forms a matrix of localizing coefﬁcients that quantify the relative importance of each grid parameter to the individual production

D. Devegowda et al. / Advances in Water Resources 33 (2010) 129–145

141

Fig. 20. The updated model realizations shown on the right for four randomly selected members clearly show a high correlation with the prior geological models.

Fig. 21. The log-permeability histogram (b) for one selected realization in comparison with the prior in (a).

Fig. 22. Streamline-derived sensitivity-based covariance localization restricts model updates to ﬂow relevant regions and the loss of variance in the updated members is very minimal.

Fig. 23. Without localization, the posterior ensemble-derived covariance matrix shows a loss of variance in the trailing eigen-directions. However, by using ﬂow relevant conditioning, the eigen-spectrum shows minimal loss of variance.

142

D. Devegowda et al. / Advances in Water Resources 33 (2010) 129–145

Fig. 24. The Goldsmith pilot study area showing the location of the producers in yellow and the injectors in blue. (For interpretation of the references to colour in this ﬁgure legend, the reader is referred to the web version of this article.)

Fig. 25. (a) The initial log-permeability distribution and the corresponding histogram for one randomly selected realization. (b) The posterior estimate following data assimilation using EnKF while and (c) The posterior estimate using EnKF with streamline-based covariance localization.

responses. In all cases studied here, the sensitivity to permeability variation is used to update the parameter ﬁeld (the permeability). The computation of the matrix q is performed at each update step to reﬂect changing ﬁeld conditions or updated well conﬁgurations.

While there are several approaches to compute sensitivity coefﬁcients, streamline models are particularly well suited for this purpose [7,8,10,28,29,35]. By taking advantage of the analytic relationship between streamline characteristics such as the time-

D. Devegowda et al. / Advances in Water Resources 33 (2010) 129–145

143

of-ﬂight and reservoir parameters like permeability and porosity, sensitivity estimation can be obtained simultaneously with the streamline simulation with very little additional computational burden [8,10]. The efﬁciency of the computation is a direct result of the fact that the parameter sensitivities can be expressed as 1D integrals along the streamline trajectories [7,8,10,32,35].

based localization appears to be promising and offers the possibility of dramatically reducing the associated computational burden by allowing the use of modest ensemble size in practical applications.

3.4.3. Application: synthetic example We now discuss the performance of the streamline-sensitivitybased covariance localization scheme using the synthetic example previously introduced in Section 2.1. We use 50 ensemble members here. The matches to the water-cut performance data at each well are shown in Fig. 19. Fig. 20 shows a randomly selected set of four updated permeability ﬁelds. For comparison purposes, we have also shown the corresponding initial permeability. In contrast to the traditional EnKF, the changes here are more localized and much of the initial geologic features are preserved. Also, parameter overshooting is largely eliminated. There are no localized zones of high and low permeability in the updated realizations of Fig. 20. Additionally, without covariance localization, the prior multimodal parameter histograms tend to become Gaussian. However, in this case, prior model statistics are preserved and this is reﬂected in the respective histograms in Fig. 21. It is important to note that all of this was achieved using a modest ensemble size of 50 members. The tendency of the EnKF to create a variance-deﬁcient ensemble is also mitigated to a large extent. This is indicated in Fig. 22, which plots the total variance versus assimilation time with and without localization. The loss of variability is reﬂected by the eigen-spectrum of the model parameter covariance matrices which is an indicator of the energy contained in different eigen-directions. Dramatic decay rates of the eigen-spectrum can be viewed as insufﬁcient projection of the variance along certain axes of the hyper-ellipse that represents the parameter covariance structure. To illustrate this, the eigen-value decomposition of the ensemble-derived covariance matrix is computed at each assimilation step with and without localization. This is shown in Fig. 23. The degraded covariance structure associated with the lower energy trailing eigen-values implies a smaller effective ensemble size. In conclusion, for efﬁcient data fusion for practical hydrologic and reservoir characterization problems, the use of sensitivity-

In this section, we discuss the application of streamline-derived sensitivities for covariance localization to a ﬁeld-scale example. The study area is the Goldsmith CO2 pilot project, a predominantly dolomitic reservoir in West Texas. The pilot area comprises of nine inverted 5-spot patterns covering around 320 acres and the average formation thickness is 100 feet. Signiﬁcant water-cut data for 20 years of production prior to the start of the CO2 injection is available at nine production wells and these were used to condition the permeability ﬁelds. The detailed production rates, well scheduling, well conversions and shut-in can be found elsewhere [21]. An ensemble of 50 permeability realizations was constructed using Sequential Gaussian co-simulation with well and seismic data [10]. Fig. 24 is an areal plot of the location of the Goldsmith CO2 pilot area within an extended study area showing the location of the surrounding producers and injectors. In Fig. 25a, we plot the spatial distribution of the log-permeability ﬁeld and the corresponding histogram for one randomly selected realization prior to data assimilation. The multimodal histogram of the initial log-permeability ﬁeld reﬂects the presence of multiple geologic facies with a signiﬁcant peak at lower permeability values indicating a fairly large proportion of poor reservoir quality rock. Using the standard EnKF, the ensemble of model realizations are conditioned to water-cut history and Fig. 25b shows the updated model after a sequence of updates. The loss of structure in the permeability ﬁeld is quite apparent here and considerable parameter overshooting is observed following data assimilation resulting in localized zones of extreme values of permeability. The log-permeability histogram also tends to become Gaussian and differs signiﬁcantly from the initial multimodal distribution. Repeating the water-cut assimilation with spatial localization using streamline-derived sensitivities, we can signiﬁcantly improve the quality of the solution. This is shown in Fig. 25c. The

3.5. Streamline-based covariance localization: ﬁeld-scale example

Fig. 26. Un-localized cross-covariance estimates between water-cut data (at the indicated well) and grid block permeabilities are plotted at a selected assimilation step in (a) and the corresponding values using streamline-based covariance localization are shown in (b).

144

D. Devegowda et al. / Advances in Water Resources 33 (2010) 129–145

Fig. 27. Matches to water-cut information at three selected wells with the initial predictions on top and the ﬁnal matches in the bottom row.

ﬁgure demonstrates that both the spatial continuity of the predominant geological features and the multimodal nature of the permeability histogram are maintained. This is achieved without creating physically inconsistent regions of extreme permeability values. The adjustments to the initial models are minimal and localized to the water swept zones which are physically intuitive and reasonable. This is the preferred mode of comparison to check the validity of the updates in the absence of a reference or ‘truth’ model. Additionally, the log-permeability histogram of the updated model clearly indicates that the relative proportion of various reservoir rock types is preserved. We can also interpret these results by analyzing the ensemblederived cross-covariance estimates obtained from the un-localized and the localized EnKF variants. In Fig. 26, we plot the cross-covariance at a selected assimilation step between water-cut at one of the wells and the permeabilities. The un-localized estimate shows large values even for distant locations where it is reasonable to expect that there is no correlation between the permeability and the water-cut data. However, with streamline-derived sensitivitybased localization, distant data points are not included in the analysis and the cross-covariance is shown to be more localized in the vicinity of the well. This is to be expected because the water arriving at producers travels mostly across the reservoir volume surrounded by the injectors and the producing wells and the crosscovariance estimates should be restricted to these regions of the reservoir. Consequently, the predicted changes should also be preferentially found in this area. This is clearly not the case with the un-localized cross-covariance estimate used in the standard EnKF formulation. Therefore, although the producers are completed in the center of the ﬁeld, the largely indiscriminate updates contribute to the loss of structure in the permeability ﬁeld during data assimilation when no localization is applied as shown in Fig. 25b. In Fig. 27 we show the matches to production history for three arbitrarily selected wells. Because of new wells and changing well conditions in the ﬁeld, the streamlines as well as the regions of localization at each assimilation step can vary and this is easily accounted in our proposed approach. At each assimilation step, we simply repeat the procedure of tracing the streamlines and then identify the corresponding region of inﬂuence for each well in order to localize the crosscovariance estimates [5].

4. Summary and conclusions The adverse impacts of inaccurate reservoir or groundwater model forecasts can compromise the economic feasibility of a project or signiﬁcantly alter the conclusions of comprehensive engineering studies. Accurate subsurface characterization therefore becomes vital for evaluating various resource management strategies. The EnKF is a viable and promising approach to sequential model updating to reconcile geologic models to dynamic data. The primary advantages of the EnKF are computational ease, portability across diverse numerical ﬂow simulation software and the ability to simultaneously assimilate data and quantify model uncertainties. However, the EnKF is very sensitive to ensemble size and the noise resulting from random statistical variability is shown to lead to degraded model updates and geologically implausible subsurface models thereby severely compromising model forecasts, especially with modest ensemble sizes. In this paper, we propose eliminating or reducing the severity of the statistical noise by conditioning the cross-covariance matrix between the data and the model unknowns and this is achieved through the use of ﬂow relevant streamline-derived information. We compare the proposed approach with the traditional implementation of the EnKF and demonstrate that the use of covariance localization via the use of the streamlines signiﬁcantly improves EnKF performance and leads to signiﬁcant savings in computational effort by reducing the need for large ensemble sizes. The primary conclusions to be drawn from this work are: Spatial localization of the ensemble-derived correlations is a necessary step prior to the process of assimilation. The selected covariance localization methodology should be related to the underlying ﬂow patterns and the geologic structure so that updates are consistent with the prior static information. Streamline trajectory based and streamline-derived sensitivitybased localization schemes are efﬁcient and effective methods to regionalize the cross-covariance to restrict the inﬂuence of the observations on the parameters. Distance-dependent localization may not be appropriate for highly heterogeneous reservoir mapping and the effectiveness is likely to be applicationdependent.

D. Devegowda et al. / Advances in Water Resources 33 (2010) 129–145

Considerable speed-up of the assimilation process is achievable using information derived from streamlines, because the accuracy of the solution is similar to an EnKF assimilation sequence using a much larger ensemble. The extra step of conditioning the cross-covariance far outweighs the need to maintain and update a much larger ensemble thereby making ﬂow relevant conditioning techniques practical enough to use for operational problems. In general, the use of sensitivities allows us to reduce the ensemble sizes further. However, for some data types, for e.g. head measurements, the sensitivity-based approach is difﬁcult to implement because of the formulation of the pressure sensitivities requires signiﬁcantly extra effort.

References [1] Akella S, Efendiev Y, Datta-Gupta A. A coarse-scale constrained ensemble Kalman ﬁlter for subsurface characterization. Water Resour Res; submitted for publication. [2] Anderson JL. An ensemble adjustment ﬁlter for data assimilation. Monthly Weather Review 2001;129(12):2903–94. [3] Anderson JL. An adaptive covariance inﬂation error correction algorithm for ensemble ﬁlters. Tellus 2007;59(2):210–24. [4] Anderson JL, Anderson SL. A Monte-Carlo implementation of the nonlinear ﬁltering problem to produce ensemble assimilations and forecasts. Mon Weather Rev 1999;127(12):2741–58. [5] Arroyo E, Devegowda D, Datta-Gupta A, Choe J. Streamline assisted ensemble Kalman ﬁlter for rapid and continuous reservoir model updating. SPERE 2008;11(6):1046–60. [6] Bishop CH, Etherton BJ, Majumdar SJ. Adaptive sampling with the ensemble transform ﬁlter. Part I: Theoretical aspects. Mon Weather Rev 2001;129(3):420–36. [7] Cheng H, Wen XH, Milliken WJ, Datta-Gupta A. Field experiences with assisted and automatic history matching using streamline models. In: Paper SPE 89857 presented at the SPE annual technical conference and exhibition, Houston, 24– 27 September; 2004. [8] Cheng H, Khargoria A, He Z, Datta-Gupta A. Fast history-matching of ﬁnite difference models using streamline-derived sensitivities. SPERE 2005;8(5):426–36. [9] Cheng H, Oyerinde AS, Datta-Gupta A. Compressible streamlines and threephase history matching. SPEJ 2007;12(4):375–85. [10] Datta-Gupta A, King MJ. 2007. Streamline simulation: theory and practice. Textbook series #11, Society of Petroleum Engineers, Richardson, TX. [11] Deutsch CV, Journel A. GSLIB geostatistical software library and user’s guide. Oxford University Press; 1992. [12] Devegowda D, Arroyo E, Datta-Gupta A, Douma SG. Efﬁcient and robust reservoir model updating using ensemble Kalman ﬁlter with sensitivity based covariance localization. In: Paper SPE 106144 presented at the SPE reservoir simulation symposium, Woodlands, Texas, 24–26 February; 2007. [13] Devegowda D. Streamline assisted ensemble Kalman ﬁlter – formulation and ﬁeld application. PhD dissertation, Texas A&M University; 2008.

145

[14] Emanuel AS, Milliken WJ. History matching ﬁnite-difference models with 3D streamlines. In: Paper SPE 4900 presented at the SPE annual technical conference and exhibition, New Orleans, 27–30 September; 1998. [15] Evensen G. The ensemble Kalman ﬁlter: theoretical formulation and practical implementation. Ocean Dyn 2003;53:353–67. [16] Evensen G. Data assimilation, the ensemble Kalman ﬁlter. Springer Publishers; 2006. [17] Furrer R, Bengtsson T. Estimation of high-dimensional prior and posterior covariance matrices in Kalman ﬁlter variants. J Multivariate Anal 2004;98:227–55. [18] Gaspari G, Cohn SE. Construction of correlation functions in two and three dimensions. Ofﬁce Note Series on Global Modeling and Data Assimilation, DAO Ofﬁce Note 96-03R1; 1996. [19] Gu Y, Oliver DS. The ensemble Kalman ﬁlter for continuous reservoir model updating. J Energy Resour Technol 2006;128:79–91. [20] Hamill TM, Whitaker JS. Distance-dependent ﬁltering of background error covariance estimates in an ensemble Kalman ﬁlter. Mon Weather Rev 2001;129:2776–93. [21] He Z, Yoon S, Datta-Gupta A. Streamline based production data integration under changing ﬁeld conditions. SPEJ 2002;7(4):423–36. [22] Houtekamer PL, Mitchell HL. Data assimilation using ensemble Kalman technique. Mon Weather Rev 1998;126(3):796–811. [23] Houtekamer PL, Mitchell HL. A sequential ensemble Kalman ﬁlter for atmospheric data assimilation. Mon Weather Rev 2001;129(1):123–37. [24] Jimenez E, Sabir K, Datta-Gupta A, King MJ. Spatial error and convergence in streamline simulation. SPERE 2007;10(3):221–32. [25] Johns CJ, Mandel J. A two-stage ensemble Kalman ﬁlter for smooth data assimilation. Environ Ecol Stat; in print. [26] Nævdal G, Johnsen LM, Aanonsen NI, Vefring EH. Reservoir monitoring and continuous model updating using ensemble Kalman ﬁlter. SPEJ 2005;10(1):66–74. [27] Oke PR, Sakov P, Corney SP. Impacts of localization in the EnKF and EnOI: experiments with a small model. Ocean Dyn 2007;57(1):32–45. [28] Oyerinde AS. A composite tracer analysis approach to reservoir characterization. MS thesis, Texas A&M University; 2004. [29] Oyerinde AS. Streamline based three phase history matching. PhD dissertation, Texas A&M University; 2008. [30] Pollock DW. Semi-analytical computation of path lines for ﬁnite difference models. Ground Water 1998;26(6):743–52. [31] Tarantola A. Inverse problem theory. SIAM; 2004. [32] Vasco DW, Yoon S, Datta-Gupta A. Integrating dynamic data into high resolution reservoir models using streamline-based analytic sensitivity coefﬁcients. SPEJ 1999;4(3):389–99. [33] Vega L, Rojas D, Datta-Gupta A. Scalability of the deterministic and bayesian approaches to production-data integration into reservoir models. SPEJ 2004;9(3):330–8. [34] Whitaker JS, Hamill TM. Ensemble data assimilation without perturbed observations. Mon Weather Rev 2002;130(7):1913–24. [35] Yoon S, Barman I, Datta-Gupta A, Pope GA. In-situ characterization of residual NAPL distribution using streamline based inversion of partitioning tracer tests. In: Paper SPE 52729 presented at the SPE/EPA exploration and production environmental conference, Austin, 1–3 April; 1999. [36] Zhang F, Zhang M, Hansen JA. Coupling ensemble Kalman ﬁlter with four dimensional variational data assimilation. Mon Weather Rev; submitted for publication.

Copyright © 2020 KUNDOC.COM. All rights reserved.