Improved CRM Model for Inter-Well Connectivity Estimation and Production Optimization : Case Study for Karst Reservoirs

Due to the coexistence of multiple types of reservoir bodies and widely distributed aquifer support in karst carbonate reservoirs, it remains a great challenge to understand the reservoir flow dynamics based on traditional capacitance–resistance (CRM) models and Darcy’s percolation theory. To solve this issue, an improved injector–producer-pair-based CRM model coupling the effect of active aquifer support was first developed and combined with the newly-developed Stochastic Simplex Approximate Gradient (StoSAG) optimization algorithm for accurate inter-well connectivity estimation in a waterflood operation. The improved CRM–StoSAG workflow was further applied for real-time production optimization to find the optimal water injection rate at each control step by maximizing the net present value of production. The case study conducted for a typical karst reservoir indicated that the proposed workflow can provide good insight into complex multi-phase flow behaviors in karst carbonate reservoirs. Low connectivity coefficient and time delay constant most likely refer to active aquifer support through a high-permeable flow channel. Moreover, the injector–producer pair may be interconnected by complex fissure zones when both the connectivity coefficient and time delay constant are relatively large.


Introduction
The hydrocarbon resources stored in carbonate reservoirs play an important role in new proven reserves worldwide over the last decade.Based on the type of reservoir bodies and geological origin and production response, carbonate reservoirs can be divided into three major categories, i.e., porous type, fractured-porous type, and karst type.So far, many karst carbonate reservoirs have been found in the northern Tarim basin, China.Compared with porous or fractured-porous carbonate reservoirs which are mainly distributed in the Middle East, Central Asia, and North America, domestic karst carbonate reservoirs are usually subjected to multiscale geological structures and strong heterogeneity.The coexistence of matrix, fracture, and karst caves usually leads to extremely complicated multiphase fluid dynamics in this type of reservoir.On the other hand, the lack of powerful techniques to accurately evaluate reservoir dynamic connectivity while considering the impact of aquifer support, further increases the uncertainty for finding an optimal waterflood development scheme.Therefore, it is urgent to establish a practical methodology for inter-well connectivity estimation while coupling the effect of aquifer support and subsequent real-time injection production optimization in order to improve waterflooding recovery as much as possible.
A typical workflow for analyzing waterflood performance involves a combination of reservoir characterization, geological models, and numerical simulation engines for modeling multiphase flow in porous media.This ab initio integrated type of model is useful for evaluating strategies at the larger scale, but is not suited to discrete, localized waterflood operations in karst carbonate reservoirs.Over the last decades, data-driven physics-based models have emerged as an attempt to provide a simple picture of reservoir fluid dynamics, providing some vital information, e.g., inter-well connectivity and advanced decision-making processes by scarifying part of the exact description (e.g., evolution of the saturation field over time).Production history, well logs, and interventions are common examples of hard data that are used to train a data-driven reservoir model to answer specific questions.The data-driven models for inter-well connectivity estimation mainly include capacitance-resistance (CRM) models [1][2][3][4], flow-network models [5], and inter-well numerical simulation models [6][7][8].
As a powerful approach for effective analogy between source/sink terms in hydrocarbon reservoirs and electrical conductors, the CRM only requires the most readily available production history and producers' bottom hole pressure (BHP) when available.In addition, the CRM can easily handle dynamic boundary conditions such as new wells or shut-ins.Full-field history matching and channeling detection can be carried out in minutes making it a suitable tool for real-time monitoring.The CRM concept was initiated by the work of Albertoni and Lake [9].Yousef et al. [1] advanced the model by incorporating the tank material balance concept and derived the model to multiple producers and injectors by applying superposition in space.Sayarpour et al. [2] further developed the CRM equations to several types of reservoir control volumes: tank model (CRMT), producer-based model (CRMP), and injector-producer pair based model (CRMIP).The CRM has also been integrated with other analytical tools (e.g., rate-transient analysis, (RTA)) for screening or monitoring the performance of various enhanced oil recovery (EOR) processes [10][11][12][13][14][15].Recently, Mamghaderi and Pourafshary [16] established an improved CRM to investigate the effect of layers using data from production logging tools (PLTs).Zhang et al. [17,18] presented a system of multi-layer CRM models for layered reservoirs by considering both BHP and crossflow.Holanda et al. [19] further introduced the state-space (SS) theory to describe the dynamic behavior of CRMs as a multi-input/multi-output system.Capacitance-resistance models were also incorporated with fractional flow models [20,21] to allow the prediction of oil rates.It indicated that the Koval model proposed by Cao et al. [21] may not be a good choice for mature waterfloods, but it is effective enough to revisit abrupt breakthroughs caused by active aquifer support, which widely exist in karst carbonate reservoirs.During production of karst reservoirs, large-scale fractured-vuggy units are usually treated as isolated targets for continuous waterflood, thereby the previous CRM models should be modified to supply the aquifer influx rate as needed with all contributing injection rates.
Owing to the fact that understanding reservoir fluid dynamics to achieve optimal decision-making by grid-based reservoir models is computationally expensive and frequently require large volumes of uncertain data to estimate petrophysical properties, full field-scale models are not easy for rapid reservoir analysis to make reliable decisions.As a newly developed technique for optimal decision, production optimization where a grid-based model is substituted by a useful yet tractable data-driven model for waterflood systems was developed to find the optimal well controls that maximize cumulative oil production or net present value (NPV).The constrained optimization we consider here is essentially a non-linear, least-squares problem.The efficient optimization mostly depends on the computations of gradients of objective function to control variables.Theoretically, gradient of objective function can be accurately obtained using the adjoint method that requires only two simulations for this calculation, no matter how many variables [22,23].However, calculation of the adjoint gradient by traditional reservoir simulation is inconvenient and code intrusive.Besides, it is not feasible to compute finite difference gradient when the number of control variables is large.
A viable alternative is approximation of the gradient by an ensemble optimization technique known as EnOpt, which has received appealing attention over the past years by reservoir engineers after the pioneering work of Chen et al. [24,25].Using the deterministic or standard EnOpt, hybrid constrained optimization problems were resolved by generalized non-linear programming algorithms [7,[26][27][28].Do and Reynolds [29] analyzed the deterministic EnOpt and demonstrated its close connection with other stochastic gradient algorithms.Stordal et al. [30] confirmed that deterministic EnOpt can be treated as a special example of Gaussian mutation.Jafroodi and Zhang [31] utilized the CRM model as an underlying reservoir proxy for ensemble-based production optimization and used a power law relationship to predict oil productivity.In their work, The CRM equation was only used to detect faults or low-permeability areas between wells, rather than acting as a proxy for the simulator, with the advantage of using actual production and injection data.Hong et al. [32] presented a proxy-model workflow where a grid-based model was supplemented by the useful yet tractable CRM model as a proxy for waterflood operations.The results indicated that CRM models have high potential to serve as a cogent proxy model for waterflooding related decision-making and obtain robust results that result in a near-optimal solution.Recently, a novel ensemble-based technique, stochastic-simplex-approximate-gradient (StoSAG), was developed by Fonseca et al. [33,34].The StoSAG deals with reservoir simulator as a black box and approximates gradient through the inputs and outputs of all the ensemble runs.Various theoretical analyses [35][36][37][38][39] showed that StoSAG can yield a significantly higher NPV than that obtained with the standard EnOpt.However, there are few proposals using the newly developed StoSAG for estimation of inter-well connectivity.We will explore the possibility of using StoSAG for the case where two objectives are to assimilate the production history data to infer inter-well connectivity when active aquifer support exists and to perform real-time production optimization by maximizing the NPV or cumulative oil production under hybrid non-linear constraints, respectively.
In this paper, an improved CRMIP model by coupling the effect of active aquifer support is first proposed and integrated with the newly developed StoSAG optimization algorithm for better understanding of inter-well connectivity in a waterflood operation.The improved CRM-StoSAG workflow is further employed for real-time production optimization to find the optimal injection rate at each control step by maximizing the objective function, i.e., net present value (NPV) of waterfloods with regard to typical karst carbonate reservoirs.Finally, the conclusions of this work will be provided.

Methodology
This section discusses in detail the integrated workflow for inter-well connectivity estimation and production optimization in a waterflood reservoir.

The Improved CRM-Koval Model
With regard to different CRM representations (e.g., CRMT, CRMP, and CRMIP), the two main parameter of a CRM are connectivity coefficient and time delay constant if BHP data of producers are not available.For an oil-water system, the time delay constant denotes how long a pressure wave from injectors takes to reach a producer.Due to reservoir heterogeneity, time delays of displacing fluid from adjacent injectors to a certain producer may differ from each other significantly.So, assuming only a time delay constant for each producer, as in CRMP, may be unreliable.In such a case, it is necessary to develop one continuity equation for every injector-producer pair, i.e., the CRMIP model.However, the traditional CRMIP model is not suited to investigate the impact of aquifer support on inter-well connectivity estimates, which are widely distributed in karst carbonate reservoirs.As shown in Figure 1, when aquifer support is inevitable, the governing differential equation should be modified, as follows where q ij (t) is the liquid production rate of an injector-producer pair at time t, m 3 /d; τ ij is the time delay constant, d; e wij is the water influx rate, m 3 /d; I i (t) is the injection rate of injector i, m 3 /d; J ij is the liquid production index of an injector-producer pair, m 3 /(MPa•d); P w f ,j is the bottom hole pressure of producer j at time t, MPa; f ij is the inter-well connectivity coefficient between injector i and producer j, ∈ [0,1].The semi-analytical solution of Equation (1) will be further derived using superposition in space, which is stated as According to superposition in time, Equation (2) has the following form For injector-producer pairs, the liquid production rate  ( ) of producer j at time  can be given by When aquifer support is active and BHP change of producer with time is negligible, there are four unknown parameters for each injector-producer-pair based control volume, i.e.,  ,  ( ),  and  .The total number of unknown parameters for waterflood reservoir is equal to 4 ×  ×  .Obviously, when the waterflood reservoir is not affected by aquifer support, the semi-analytical solution Equation ( 3) can be simplified as the traditional CRMIP model.In order to guarantee a reasonable balance between injection and production, the inter-well connectivity coefficient  should satisfy inequality constraints.In many cases, the liquid production rate  ( ) must be subjected to some equality constraints as follows The semi-analytical solution of Equation (1) will be further derived using superposition in space, which is stated as According to superposition in time, Equation (2) has the following form For injector-producer pairs, the liquid production rate q j (t k ) of producer j at time t k can be given by When aquifer support is active and BHP change of producer with time is negligible, there are four unknown parameters for each injector-producer-pair based control volume, i.e., f ij , q ij (t 0 ), τ ij and e wij .The total number of unknown parameters for waterflood reservoir is equal to 4 × N pro × N inj .Obviously, when the waterflood reservoir is not affected by aquifer support, the semi-analytical solution Equation (3) can be simplified as the traditional CRMIP model.In order to guarantee a reasonable balance between injection and production, the inter-well connectivity coefficient f ij should satisfy inequality constraints.In many cases, the liquid production rate q ij (t 0 ) must be subjected to some equality constraints as follows If the two-phase flow in each control volume can be regarded as steady-state flow, the Koval model proposed by Cao et al. [21] will be used to compute the fractional flow of water in porous media by considering the effect of local heterogeneity and viscosity ratio, which is given by where K val is the Koval factor, reflecting both reservoir heterogeneity and fluid-viscosity contrast.
A large Koval factor (greater than unity) usually implies either a high degree of reservoir heterogeneity or a large oil-water viscosity ratio.t D is a dimensionless time denoting the cumulative water injection in control volume.
where f ij is the interwell connectivity coefficient, which can be estimated using the improved CRMIP model to assimilate the production history; V pj is the drainage volume of a injector-producer pair, m 3 ; I i is the injection contribution to the producer at timestep t k , m 3 /d.For oil-water two-phase system, the oil or water production rate of jth producer at timestep t k can be easily obtained based on the physical meaning of fraction-flow equation, which can be written as [21] q wj (t k ) = q j (t k ) f wj (t k ) Using the improved CRM-Koval model for inter-well connectivity estimation, the six unknown parameters for each injector-producer pair, i.e., connectivity coefficient, time delay constant, water influx rate, liquid production rate at timestep t 0 , Koval factor, and drainage volume were estimated with non-linear multivariate regression, the required least-squares objective function were finally described as minmize Except for Equations ( 5) and ( 6), the objective function was also constrained by where the superscript obs and cal denote the observed and predicted production data, respectively; V pField denotes the total pore (drainage) volume of a reservoir or block, m 3 .

Waterflood Production Optimization
When the improved CRM-Koval model was used for a better understanding of reservoir fluid dynamics including inter-well connectivity, aquifer influx rate, etc., the NPV of production was defined as the objective function by having the proxy model serve as a precursor of a grid-based reservoir model and finally maximized in order to find the optimal well controls (rate or pressure) under hybrid non-linear constraints, which can be given as follows [22] where u is a N u -dimensional column vector containing all the well controls over the production lifetime.
For the problem of interest here, we treated the injection rate of each injector as the control variable; r o is the oil revenue, USD/STB; c w and c wi , respectively, are the disposal cost of produced water and the cost of water injection, USD/STB; b is annual discount rate; t n denotes the time at end of the nth time step of the proxy model; ∆t n is the size of the nth time step; N t is the total number of time steps; P and I denote the number of producers and injectors, respectively; q n o,j and q n w,j respectively, denote the average oil and water production rate at the jth producer, STB/day; q n wi,k is the average water injection rate at the kth injector, STB/day.
The generalized waterflood production optimization problem can be written as Satisfying the following constraints where N(u) is the objective function for waterflood optimization; Equations ( 16)-( 18) denote the boundary constraint, inequality, and equality constraints, respectively; u low i and u up i are the lower and upper limits of ith control variable u i , respectively; n i and n e denote the number of inequality and equality constraints, respectively.

Ensemble-Based Optimization Method
The key to sequential data assimilation or robust optimization in large-scale non-linear dynamics is to obtain the gradients of objective function and handle the hybrid inequality and equality constraints.Due to the complexity of fluid flow in heterogeneous porous media, accurate computation of gradient with the adjoint or finite difference method in a traditional reservoir simulation code is usually time consuming and code intrusive.To overcome the drawbacks, an ensemble-based method, StoSAG algorithm, was used to obtain the approximate gradient of least-squares objective function.The StoSAG algorithm can treat the reservoir simulator or proxy model as a black box and approximates the gradient of objective function through the inputs and outputs of all the ensemble runs efficiently.Moreover, the augmented Lagrange method and log-transformation method will be integrated to handle the hybrid nonlinear constraints.

Augmented Lagrange Objective Function
For the kth control variable u k , an unbounded optimization problem can be obtained by applying a log transformation method [23] to transform the original bound-constrained optimization problem, which takes the form of where u low k and u up k are the lower and upper bounds of the kth control variable u k , respectively; v k is the log-domain kth control variable varying from −∞ to ∞.
The robust optimization is performed in log-domain, but the control variable u k can be obtained with the inverse log-transformation after each iteration, which can be illustrated as Using the log-transformation to eliminate the bound constraints, the following augmented Lagrange objective function [15] will be ultimately established where λ e,j and λ c,i , respectively, denote the Lagrangian multiplier of equality and inequality constraints; µ denotes the penalty parameter; s e,j and s c,i , respectively, denote the scaling factors for equality and inequality constraints.Note that the term N(u) in Equation ( 14) is substituted with L a (u, λ, µ) for a hybrid constrained optimization problem.

StoSAG Gradient Computation
To address the data assimilation problem of Equation (11) and waterflood optimization problem of Equation ( 15), we used the steepest ascent algorithm [24] to compute the estimate of the optimal control vector at the (k + 1)th iteration, which is given by where u 0 is the initial guess; u k is the optimal control vector at the kth iteration.a k is the step size; d k denotes the search direction vector.For our study, it was a stochastic search direction.The StoSAG algorithm proposed by Fonseca et al. [33] with the stochastic search direction, described as following, was used to maximize the augmented Lagrange objective function L a (u, λ, µ): where δ U k j = u k j − u k , and δL a | k j = L a u k j − L a u k .The superscript "+" on a matrix denotes the Moore-Penrose pseudo-inverse.the superscript "T" denotes the transpose process for a vector or matrix.Note that, u k j denotes N e ensembles of Gaussian random vector, where X~N u k , I (i.e., the mean of u k j equals to u k and its covariance matrix I), and can be written as where k is the iteration number of inner loop; u k is the optimal control vector at the kth iteration; L is a lower triangular matrix obtained by Cholesky decomposition of the covariance matrix C U ; Z j is the vector satisfying a Gaussian distribution N(0, I), and I is unit matrix with dimension of N u × N u .Therefore, L•Z j denotes a (N u × 1) Gaussian random vector, namely, L•Z j ∼ N(0, C U ).The following spherical covariance function is further adopted to compute the (i, j) entry of covariance matrix C U by C i,j , Figure 2 displays the integrated CRMIP-StoSAG workflow for inter-well connectivity estimation and waterflood optimization.This workflow was divided into two parts: data assimilation by the CRMIP-Koval model (shaded in orange) and real-time waterflood optimization based on good understandings of inter-well connectivity relationships (shaded in blue).
Energies 2019, 12, x 8 of 15 is the vector satisfying a Gaussian distribution (0, ), and  is unit matrix with dimension of  ×  .Therefore,  ⋅  denotes a ( × 1) Gaussian random vector, namely,  ⋅  ∼ (0,  ).The following spherical covariance function is further adopted to compute the (i, j) entry of covariance matrix  by  , , Figure 2 displays the integrated CRMIP-StoSAG workflow for inter-well connectivity estimation and waterflood optimization.This workflow was divided into two parts: data assimilation by the CRMIP-Koval model (shaded in orange) and real-time waterflood optimization based on good understandings of inter-well connectivity relationships (shaded in blue).

Case Study for Karst Carbonate Reservoir
As previously mentioned, when the impact of aquifer support is negligible, the improved CRMIP model will be simplified to the traditional CRM representation as described in Sayarpour et al. [2].Accuracy of the improved CRM-StoSAG workflow and its feasibility for strongly-heterogeneous waterflood reservoir was validated in our previous work [15], which will not be revisited here.The following focuses on using the improved CRM-StoSAG workflow to infer the inter-well connectivity relationship of karst carbonate reservoirs and subsequent waterflood optimization.
The typical karst carbonate reservoir we consider here is located in northern Tarim Basin, which is subject to strong heterogeneity and contains multiple types of reservoir bodies in many cases, mainly including karst caves, fractured-vuggy bodies, etc. [40].As a matter of fact, karst caves are usually regarded as the major reservoir spaces.Rapid drilling time drops, drilling rigs ventilation, overflow and large mud loss often occurs.In addition, strong bead reflection, as shown in Figure 3a, can be utilized to determine the karst cave systems.As for the fractured-vuggy reservoir bodies, the main components are composed of dissolved pores and small caves.In addition to communicating the dissolved pores and caves, a relatively high hydrocarbon storage is founded in high-angle tectonic fractures with varying conductivities.The main seismic amplitude for this type of reservoir body is strong flake and weak messy reflection, see Figure 3b.As previously mentioned, when the impact of aquifer support is negligible, the improved CRMIP model will be simplified to the traditional CRM representation as described in Sayarpour et al. [2].Accuracy of the improved CRM-StoSAG workflow and its feasibility for stronglyheterogeneous waterflood reservoir was validated in our previous work [15], which will not be revisited here.The following focuses on using the improved CRM-StoSAG workflow to infer the inter-well connectivity relationship of karst carbonate reservoirs and subsequent waterflood optimization.
The typical karst carbonate reservoir we consider here is located in northern Tarim Basin, which is subject to strong heterogeneity and contains multiple types of reservoir bodies in many cases, mainly including karst caves, fractured-vuggy bodies, etc. [40].As a matter of fact, karst caves are usually regarded as the major reservoir spaces.Rapid drilling time drops, drilling rigs ventilation, overflow and large mud loss often occurs.In addition, strong bead reflection, as shown in Figure 3a, can be utilized to determine the karst cave systems.As for the fractured-vuggy reservoir bodies, the main components are composed of dissolved pores and small caves.In addition to communicating the dissolved pores and caves, a relatively high hydrocarbon storage is founded in high-angle tectonic fractures with varying conductivities.The main seismic amplitude for this type of reservoir body is strong flake and weak messy reflection, see Figure 3b.In addition, due to the complex distribution patterns of active aquifer support in karst carbonate reservoir, water influx behaviors observed mainly conform to abrupt watered-out and fluctuation or staircase rise characteristics, as shown in Figure 4.Moreover, the lack of efficient techniques to replenish formation void age may further result in a relatively low recovery efficiency and high natural decline rate of oil production.On the other hand, the majority of previously-used mathematical models are based on Darcy's flow, and the coexistence of porous and free-flow domains over a wide range of scales usually lead to sophisticated fluid flow dynamics in karst carbonate reservoirs.The traditional percolation theory and reservoir simulators are not adapted to represent the coupled flow dynamics in this type of waterflood reservoir.In addition, due to the complex distribution patterns of active aquifer support in karst carbonate reservoir, water influx behaviors observed mainly conform to abrupt watered-out and fluctuation or staircase rise characteristics, as shown in Figure 4.Moreover, the lack of efficient techniques to replenish formation void age may further result in a relatively low recovery efficiency and high natural decline rate of oil production.On the other hand, the majority of previously-used mathematical models are based on Darcy's flow, and the coexistence of porous and free-flow domains over a wide range of scales usually lead to sophisticated fluid flow dynamics in karst carbonate reservoirs.The traditional percolation theory and reservoir simulators are not adapted to represent the coupled flow dynamics in this type of waterflood reservoir.According to the problems of rapid production decline and low recovery rate caused by weak natural energy and rapid rise of water cut, pilot tests of sequential waterflooding have been carried out for some potentially interconnected injector-producer pairs, which are sorted out according to tracer surveillance and interference well testing.But with the amount of water injection increasing, oil yield effect gradually becomes worse, and much of the remaining oil is still unexploited around producers.For this study, typical injector-producer-pairs, Pro_1, Inj_1, and Pro_2, corresponding to the same large-scale fractured-vuggy unit, was screened for application for the improved CRM-StoSAG workflow.Using the liquid and oil production data of Pro_1, Inj_1, and Pro_2 with a relatively stable production scheme for history matching, the control variables for inter-well connectivity including connectivity coefficient, time delay constant, water influx rate, Koval factor, and drainage volume of each injector-producer pair were ultimately estimated and summarized in Table 1.The CRMIP oil production match results of Pro_1 and Pro_2 are shown in Figures 5-6, respectively.The iteration procedure of the least-squares objective function using the improved CRM-StoSAG workflow for inter-well connectivity estimation is illustrated in Figure 7.According to the problems of rapid production decline and low recovery rate caused by weak natural energy and rapid rise of water cut, pilot tests of sequential waterflooding have been carried out for some potentially interconnected injector-producer pairs, which are sorted out according to tracer surveillance and interference well testing.But with the amount of water injection increasing, oil yield effect gradually becomes worse, and much of the remaining oil is still unexploited around producers.For this study, typical injector-producer-pairs, Pro_1, Inj_1, and Pro_2, corresponding to the same large-scale fractured-vuggy unit, was screened for application for the improved CRM-StoSAG workflow.Using the liquid and oil production data of Pro_1, Inj_1, and Pro_2 with a relatively stable production scheme for history matching, the control variables for inter-well connectivity including connectivity coefficient, time delay constant, water influx rate, Koval factor, and drainage volume of each injector-producer pair were ultimately estimated and summarized in Table 1.The CRMIP oil production match results of Pro_1 and Pro_2 are shown in Figures 5 and 6, respectively.The iteration procedure of the least-squares objective function using the improved CRM-StoSAG workflow for inter-well connectivity estimation is illustrated in Figure 7.It demonstrates that, the estimated results of connectivity coefficient and time delay constant between Pro_2 and Inj_1 was far less than those between Pro_1 and Inj_1, suggesting that the injected water makes little contribution to oil production, which is mainly influenced by active aquifer support, and the most likely geological structure between Pro_2 and Inj_1 was high-permeable flow channel.Moreover, the estimated water influx rate of Pro_1 was lower than 1.0 m 3 /d, indicating no active aquifer support from Inj_1.The remarkable oil yielding effect was mainly dependent on water injection.It is inferred that Pro_1 and Inj_1 are interconnected with a complex fissure zone, which is in good agreement with the understandings obtained from interference well testing.When performing production optimization for Pro_1, Inj_1, and Pro_2 pairs, water injection for Pro_1 needs to be properly increased along with real-time adjustment of operation parameters, e.g., nozzle size, shut-off, etc., in order to prohibit too-early water breakthroughs from nearby injectors.
Based on the accurate estimation of inter-well connectivity relationships, water injection rates of Inj_1 at each control step were regarded as unknown control variables, and the NPV of oil production with regard to this large-scale fractured-vuggy unit was maximized using the StoSAG optimization method in order to find the optimal well control variables.Figure 8 displays the optimized distribution of water injection rate at each timestep.Comparison of well group cumulative oil production prior to and after optimization is shown in Figure 9.The result indicates that, when cumulative volume of water injection is assumed constant, the cumulative oil increase of this well group after real-time injection-production optimization in karst carbonate reservoir was equal to 1290.2 m 3 , indicating a remarkable oil yielding effect.Obviously, the proposed workflow can provide good insights into accurate estimation of dynamic connectivity and subsequent waterflood optimization in karst carbonate reservoirs with aquifer support widely existing.It demonstrates that, the estimated results of connectivity coefficient and time delay constant between Pro_2 and Inj_1 was far less than those between Pro_1 and Inj_1, suggesting that the injected water makes little contribution to oil production, which is mainly influenced by active aquifer support, and the most likely geological structure between Pro_2 and Inj_1 was high-permeable flow channel.Moreover, the estimated water influx rate of Pro_1 was lower than 1.0 m 3 /d, indicating no active aquifer support from Inj_1.The remarkable oil yielding effect was mainly dependent on water injection.It is inferred that Pro_1 and Inj_1 are interconnected with a complex fissure zone, which is in good agreement with the understandings obtained from interference well testing.When performing production optimization for Pro_1, Inj_1, and Pro_2 pairs, water injection for Pro_1 needs to be properly increased along with real-time adjustment of operation parameters, e.g., nozzle size, shut-off, etc., in order to prohibit too-early water breakthroughs from nearby injectors.
Based on the accurate estimation of inter-well connectivity relationships, water injection rates of Inj_1 at each control step were regarded as unknown control variables, and the NPV of oil production with regard to this large-scale fractured-vuggy unit was maximized using the StoSAG optimization method in order to find the optimal well control variables.Figure 8 displays the optimized distribution of water injection rate at each timestep.Comparison of well group cumulative oil production prior to and after optimization is shown in Figure 9.The result indicates that, when cumulative volume of water injection is assumed constant, the cumulative oil increase of this well group after real-time injection-production optimization in karst carbonate reservoir was equal to 1290.2 m 3 , indicating a remarkable oil yielding effect.Obviously, the proposed workflow can provide good insights into accurate estimation of dynamic connectivity and subsequent waterflood optimization in karst carbonate reservoirs with aquifer support widely existing.

Conclusions
We develop an improved CRMIP model by coupling the effect of active aquifer support and integrated it with the newly developed StoSAG optimization algorithm for accurate evaluation of reservoir dynamic connectivity in a waterflood operation.Then, the improved CRM-StoSAG workflow was employed for real-time waterflood production optimization in order to find the optimal water injection rate at each control step by maximizing the objective function, i.e., net present value (NPV) of production.
Case studies showed that the proposed workflow can provide good insights into accurate estimation of inter-well connectivity and subsequent waterflood optimization.With regard to typical karst carbonate reservoirs, low connectivity coefficient and time delay constant most likely refer to active aquifer support through a high-permeable flow channel.The injector-producer pair may also be interconnected by complex fissure zones when both the estimated connectivity coefficient and time delay constant are relatively large.

Conclusions
We develop an improved CRMIP model by coupling the effect of active aquifer support and integrated it with the newly developed StoSAG optimization algorithm for accurate evaluation of reservoir dynamic connectivity in a waterflood operation.Then, the improved CRM-StoSAG workflow was employed for real-time waterflood production optimization in order to find the optimal water injection rate at each control step by maximizing the objective function, i.e., net present value (NPV) of production.
Case studies showed that the proposed workflow can provide good insights into accurate estimation of inter-well connectivity and subsequent waterflood optimization.With regard to typical karst carbonate reservoirs, low connectivity coefficient and time delay constant most likely refer to active aquifer support through a high-permeable flow channel.The injector-producer pair may also be interconnected by complex fissure zones when both the estimated connectivity coefficient and time delay constant are relatively large.

Figure 1 .
Figure 1.Schematic diagram of an injector-producer-pair-based control volume that is widely distributed in karst waterflood carbonate reservoirs.

Figure 1 .
Figure 1.Schematic diagram of an injector-producer-pair-based control volume that is widely distributed in karst waterflood carbonate reservoirs.

Figure 2 .
Figure 2. Workflow for inter-well connectivity estimation and waterflood optimization.

Figure 2 .
Figure 2. Workflow for inter-well connectivity estimation and waterflood optimization.

Figure 3 .
Figure 3. Seismic amplitude of different reservoir bodies in the Tarim Basin: (a) strong bead-like reflection for karst caves; (b) strong flake-like reflection for fractured-vuggy bodies.

Figure 3 .
Figure 3. Seismic amplitude of different reservoir bodies in the Tarim Basin: (a) strong bead-like reflection for karst caves; (b) strong flake-like reflection for fractured-vuggy bodies.

Figure 4 .
Figure 4. Schematic diagram of typical water influx behavior patterns: (a) abrupt watered-out; (b) fluctuation or staircase rise in water cut.

Figure 4 .
Figure 4. Schematic diagram of typical water influx behavior patterns: (a) abrupt watered-out; (b) fluctuation or staircase rise in water cut.

Figure 7 .
Figure 7. Iteration procedure of objective function based on the CRM-Koval model.

Figure 7 .
Figure 7. Iteration procedure of objective function based on the CRM-Koval model.

Figure 7 .
Figure 7. Iteration procedure of objective function based on the CRM-Koval model.

Figure 7 .
Figure 7. Iteration procedure of objective function based on the CRM-Koval model.

Figure 8 .
Figure 8. Distribution of water injection rate obtained by waterflood optimization.

Figure 8 .
Figure 8. Distribution of water injection rate obtained by waterflood optimization.

Figure 9 .
Figure 9.Comparison of cumulative oil production before and after optimization.

Figure 9 .
Figure 9.Comparison of cumulative oil production before and after optimization.

Table 1 .
The estimated control variables of CRMIP-Koval model.

Table 1 .
The estimated control variables of CRMIP-Koval model.

Table 1 .
The estimated control variables of CRMIP-Koval model.

Table 1 .
The estimated control variables of CRMIP-Koval model.