Joint Inversion of Fracture Model Properties for CO2 Storage Monitoring or Oil Recovery History Matching

Résumé — Inversion conjointe des propriétés d’un modèle de fractures pour le monitoring d’un stockage de CO 2 ou le calage d’un historique de production — Que ce soit pour la production de pétrole ou le stockage de CO 2 , le terme réservoir est utilisé régulièrement. Un réservoir désigne une structure géologique hétérogène en types de roches ou en termes de propriétés de perméabilité et de porosité. Dans le domaine pétrolier, de nombreux réservoirs présentent des fractures ou des failles. Les réservoirs fracturés occupent une part importante de la production de pétrole dans le monde (Moyen Orient, Golfe du Mexique, etc.). Malgré la présence de fractures, ces réservoirs n’en restent pas moins de bons réservoirs. Des études de réservoirs dédiés au stockage de CO 2 ont montré de plus la présence de fractures diffuses ou de failles ainsi qu’un fort impact de celles-ci dans le transfert du CO 2 . Un facteur clé de la connaissance des réservoirs fracturés est la compréhension de la géométrie et de la conductivité hydraulique du réseau formé par les fractures. Cette compréhension nécessite la construction d’un modèle de Abstract — Joint Inversion of Fracture Model Properties for CO 2 Storage Monitoring or Oil Recovery History Matching — For oil recovery or CO 2 storage, “reservoirs” are commonly used to designate geological structures where oil can be found or CO 2 can be stored. All reservoirs present a heterogeneity in terms of rock type and properties (such as porosity and permeability). In addition, some of these reservoirs present fractures and faults. Fractured reservoirs are an important part of the oil reserves in


INTRODUCTION
The concept of reservoir is used in geosciences to represent any porous and permeable structure. This structure is associated with an impermeable layer that stops fluid migration and keeps it inside the reservoir formation. Because of their geological past, some of these reservoirs contain faults and fracture networks. Modeling the complex geometries and fluid flows inside these reservoirs represents a challenging task and has been an active research topic in the field of oil recovery [1]. Another growing application of fracture characterization and modeling deals with the impact of faults and fractures on the safety of underground storage of CO 2 . These works involve natural analogs studies [2] or the characterization of reservoir sites [3,4]. They underline the importance of the characterization and modeling of faults and fractures. The presence of faults can seriously affect the the world (Middle East, Gulf of Mexico, etc.) and some of them are important reservoirs in terms of oil volume and productivity in spite of the fractures. In addition, studies of reservoirs for geologic storage of CO 2

have shown the existence of diffuse fractures and faults and their strong impacts on flow.
A key point in fractured reservoirs is to understand the geometry and hydraulic conductivity of the network formed by the fractures. This requires the construction of a reservoir model that integrates all available conceptual knowledge and quantitative data. The topic of the present paper deals with a new methodology able to perform the history matching of a fractured reservoir model by adapting the sub-seismic fault properties and positions. The main difficulty of this work is to generate a sub-seismic fault network whose fault positions can be easily modified while respecting the statistical fault model. The sub-seismic fault model we have chosen allows us to obtain a sub-seismic fault network that is consistent with the seismic fault network and that succeeds in capturing the specific spatial organization of the faults. In a first step, the geometry of the seismic fault network is characterized using fractal methods. Sub-seismic faults are then generated according to a stochastic algorithm. Finally, the geometry of this discrete fracture network is optimized in order to match the hydrodynamic data about the reservoir. The optimization algorithm modifies the sub-seismic fault positions, leading to the history matching of the reservoir model. Fractal properties are preserved during the deformation process. These different steps are demonstrated on a realistic synthetic case.
hydrodynamic behavior of the reservoir [5,3]. Conductive faults will tend to create preferential flow paths to the wells. These can lead to early water breakthrough. At the opposite, sealing faults will cause compartmentalization problems [6]. An accurate description of the fractured network is hence necessary to understand and predict the hydrodynamic behavior of the reservoir.
The first question in fractured medium characterization is to figure out the working. A fractured network can be divided into different sets having their specific characteristics. Usually, fractures are divided in two main types: the background fractures, mainly composed of joints (tens of meters), and the faults (hundreds of meters to kilometers).
Often, the density of background fractures is large, which facilitates the identification of fracture sets and the characterization of their geometrical parameters from observations on borehole images or on analogs [7]. Due to the large number of observed fractures, background fracture sets are well constrained. This is not the case with seismic faults (several hundred meters to kilometers).
Fault network geometry is determined from seismic images [5]. The critical point is that all faults cannot be detected because of the resolution limit of seismic methods. For example, a vertical fault throw may be too small to be visible on seismic images. These undetectable faults are called sub-seismic faults (hundreds of meters to kilometers). In addition, the seismic fault density is not high enough to build a statistical fault model in the same way as diffuse fracture models. Nevertheless, it is possible to capture the fault organization using fractal theory and to generate sub-seismic faults constrained by a fractal or correlation dimension.
Fractal theory is used to describe objects -here, fracture networks -that display some similarity whatever the observation scale. Many natural processes exhibit this fractal behavior [8,9], and are characterized using the non-integer fractal dimension, D f . For a 2D network in which fractures are represented using lineaments, D f varies between 1 and 2. When D f is close to 1, the fracture network is typically made of distinct linear objects. When it is close to 2, the fractures tend to uniformly fill the 2D space.
An extension of fractal theory, the multifractal theory, is classically used to describe quantities such as fault density. Multifractals are characterized by a set of generalized dimensions D q [10]. These dimensions can be easily measured using the multifractal box-counting algorithm.
The study domain is subdivided in n(r) square cells of size r. A quantity F i is measured in each box (cumulated fault length or number of fracture barycenters) and normalized weights p i are computed: The q-order moments are computed: If the property is multifractal, these moments scale according to: for -∞< q < ∞. The dimensions D q form the multifractal spectrum of generalized dimensions. D 0 is equivalent to the fractal dimension D f . D 2 is the correlation dimension, D c , which will be the main focus in this paper.
Another method for inferring D c is to consider the faults midpoints and to compute the function representing the number of points whose mutual distance is less than or equal to r. If the population of points is a scaling structure with a fixed correlation dimension, this function has a power law behavior with exponent D c .
Stochastic methods based on fractal geometry have been extensively used to generate sub-seismic fault models [11][12][13][14]. These methods assume that a fracture network is similar at different observation scales. They are convenient because they propose to analyze the available seismic faults and extend their properties to smaller observation scales.
Another question deals with the choice of the fluid flow modeling approach. The flow simulation cannot be carried out on the geological grid because of too large computational requirements. Computations of equivalent flow parameters and upscaling of the properties are therefore used to convert object-based models into reservoir simulation models. Because of the uncertainties involved in the fracture modeling process, these models are not able to correctly reproduce the past hydrodynamic behavior of the reservoir. As a consequence, these models have poor production predictability and need to be adjusted to match the historical data.
Traditional history matching methods, in which petrophysical parameters (permeabilities and porosities) are adjusted, are not efficient with sub-seismic faults. The reason is that the objects are large and scarcely distributed. The hydrodynamic behavior is therefore predominantly influenced by fault positions rather than by the petrophysical parameters. When the geometrical parameters of the faults are also adjusted, traditional assisted history matching methods directly optimize the uncertain parameters, which leads to two major problems. First, the modification of individual fault parameters (position, orientation or length) can lead to incoherent geological models; in such a case, a model matching the historical data can be found but the predictability of forecast simulations can be poor. Secondly, optimizing individual fault parameters is not practical in terms of computation cost. Recent work presented an attractive methodology to perform the history matching through an optimization of fault positions [15]. Nevertheless, the sub-seismic fault model used in this method was difficult to constrain to observed data. Obtained realizations were not realistic and this method is therefore difficult to apply to real field cases. The concepts of this paper were then revisited [16,17], leading to a new model that is able to capture the spatial organization of the seismic faults while constraining the sub-seismic faults positions to dynamic data. This paper presents the latest developments of this methodology based on a new sub-seismic fault generation and on an optimized approach of history matching. Special emphasis is put on the parameterization of fault positions and on the history matching algorithms.
In the first section, the new stochastic sub-seismic fault generator is presented. The algorithms are based on the preservation of a spatially correlated structure and designed to be compatible with our history matching methodology. The second part deals with a parameterization of the displacement of the faults. This method is able to overcome the limitations of traditional history matching methods. The fractal properties of the fault model are preserved when the geometry of the fault network are modified, thus ensuring a geologically coherent model. Considering the global fault organization and not each sub-seismic fault independently, the number of parameters is greatly reduced in order to bring the optimization down to reasonable times. Lastly, the validity of the methodology is tested on a semi-synthetic case study.

WORKFLOW TO CREATE A CONSISTENT FAULT NETWORK
Because the sub-seismic faults are by definition not observable, we must introduce some assumption on their organization.
Here, they are assumed to have the same organization as the seismic faults. The correlation dimension of the seismic faults has to be established using the box counting algorithm, for example. In this part, a method is proposed to distribute a given number of sub-seismic faults under the constraint of a given correlation dimension. The spatial distribution will be controlled by a density map. This map is computed in two steps. First, a fault center density map covering the whole fractured domain is computed. The cell size of this map is set equal to r 0 . In a second step, the multifractal behavior is extended below r 0 using the multiplicative cascade algorithm [18].

Seismic Fault Density Map
The multiplicative cascade algorithm generates multifractal images [18][19][20]. Its principle is to generate a scaling structure by recursively replicating a given pattern at different scales.
In the particular case of a 2D image, the pattern can be constructed using four weights. These are constrained by the system: where D 2 is the correlation dimension at which the scaling structure should be constrained. In this application, D 2 should be equal to the correlation dimension measured on the seismic faults. The weights are assigned to the cells of a 2 × 2 grid (Fig. 1a). Each cell is then subdivided in 4 parts and new weights are assigned to the subdivided cells. Either new weights are recomputed using the previous equation or the same values are used at all iterations (as is the case in the figure). At each stage, weights are then multiplied by the value of the parent cell (Fig. 1b). The process of subdividing cells and multiplying them by the previous iterations is repeated until the targeted resolution is obtained (Fig. 1c). The multiplicative cascade process yields the targeted correlation dimension D 2 because each subdivision honors the 2nd order box counting condition defined by Equations (1) and (2).

Fault Geometrical Properties
The density map computed with the cascade algorithm is used as a discrete probability distribution function to draw a population of sub-seismic fault centers. The two-point correlation dimension D c of these points is equal to the correlation dimension used to constrain the cascade algorithm.
Fault lineaments are generated according to a power-law length distribution and randomly distributed on the simulated centers. The orientation of each lineament is given by the average of the surrounding seismic faults. To pass from a 2D Multiplicative cascade algorithm.
model to a 3D model, lineaments are then extended along the dipping angle and clipped to the geological units of the model.

FULL FIELD UPSCALING
The next step is to model the multiphase flow in the reservoir. At this step, an obvious choice would be to take the faults explicitly into account. This approach is not applied here because it requires the construction of a new reservoir grid every time a new fracture model is built. Furthermore, gridding individual faults tends to generate very large grids requiring long simulation times. Thus, the object-based seismic and sub-seismic fault network cannot be used to directly simulate the flow and must therefore be converted into a simulation model. The chosen workflow is therefore based on an upscaling approach [21]. First, conductivity values are associated to each fault. These can be deterministic values for all the faults or can be drawn in a given distribution. An equivalent permeability tensor is then computed for each cell of the reservoir grid.
For fractured reservoir simulation, it is important to choose the adequate flow model. The reservoir engineer may choose between three models: -a single-medium porosity model, where the fractured reservoir is represented by an equivalent porous medium. There are one porosity value and one permeability value per grid cell [22]; -a dual-porosity single-permeability model where two porous media are defined. An exchange flux is introduced to model the communication between the fractures and the matrix. There are two porosity values, two permeability values associated with the matrix and the fractured medium [23,24]. Matrix block size values are used to evaluate the exchange flux between matrix and fracture cells. The use of "single-permeability" is to indicate that there is no transfer between the matrix blocks; -a dual-porosity dual-permeability model which is similar to the dual-porosity single-permeability model. The use of "dual-permeability" means that there are transfers between the matrix blocks. The model choice depends on the geometrical and physical properties of the fractured media [21]. Equivalent permeability and porosity values have to be determined in each cell of the reservoir grid. This can be achieved using one of the following two upscaling methods:

Numerical Upscaling
Numerical uspcaling is based on a numerical model. A discrete fracture model (DFN) is loaded or generated thanks to a statistical fracture model. A mesh is built based on the DFN geometry and used for flow simulations. Stationary single phase flow is computed using varying boundary conditions. The results of these simulations are used to derive the equivalent permeability tensor [21,25].

Analytical Upscaling
Equivalent permeability tensors are computed using analytical expressions based on the Oda's method [26]. If the statistical properties of the DFN are known, these properties can be directly integrated to obtain the equivalent permeability tensor [27]. In the contrary case, each fracture is studied separately. The equivalent permeability tensor is derived from the dips, strikes and volumes of the fractures [26]. Basic assumptions for this model are: -a linear pressure gradient in the fracture network; -a complete connectivity of the fracture network; -independent flow in each fracture. This approach can be used when there are no preferential flow paths (it is typically the case of diffuse fracture networks with a large density) or when there are few faults whose position has to be explicitly taken into account (it is the case of seismic faults).
Matrix and fracture porosity are given by the reservoir model while the matrix block size (a dual-porosity model parameter) is derived from the fracture set orientation and spacing.
All these options are available in our tools and have to be chosen depending on fractured reservoir properties. In this paper, the synthetic case is studied using a dual-porosity single-permeability model and an analytical upscaling. This choice is made, considering that there is a strong density of diffuse fractures. The matrix blocks are isolated and transfers occur only through the fractures. In addition, analytical upscaling allows to reduce computation time. This is very useful because a history matching loop will be run and the model will have to be recomputed for each new DFN.

INVERSE MODELING
The construction of the correlated sub-seismic fault model followed by the upscaling and the fluid flow simulation can be considered as the resolution of the direct problem. Because of the high amount of uncertainty involved in the sub-seismic fault modeling this model does not reproduce the past hydrodynamic behavior measured on the field. The objective of this paper is to present a method that is able to integrate the dynamic data into the fault model while preserving the scaling properties of the sub-seismic model.
The history matching of the sub-seismic fault realizations is formulated as an optimization problem. An Objective Function (OF), is defined to measure the difference between field (d obs ) and simulated (d sim ) production data: Minimizing OF is a difficult task because the number of parameters (size and location of the faults) is large. Furthermore, the optimized realization cannot be arbitrary and must honor the structure of the stochastic model. These difficulties can be overcome using a statistical parameterization method like the object-based gradual deformation method [15]. We first recall how random numbers can be continuously deformed without changing their statistical properties. Then, we explain how the gradual deformation method can be used to change the fault positions (without changing their spatial distribution).

The Gradual Deformation Method
The gradual deformation method proposes to deform an initial realization of a Gaussian random function in a continuous and coherent way in order to improve its match to dynamic data [28]. It is based on the fact that linear combinations of Gaussian random functions are Gaussian random functions. Let Z(x) be a standard normal random function. A new realization z of Z is constructed as a linear combination of N + 1 realizations z i of Z: (6) with weights ρ i satisfying: The constraint can be satisfied using spherical coordinates. In the case of two realizations this leads to: z t (x) = z 0 (x) cos(πt) + z 1 (x) sin(πt) The variation of t in the interval [-1, 1] yields a progressive variation of the function z t (x). The correlation between z 0 (t) and z(t) is strong if t is close to -1, 0, or 1. With t = ± 0.5, the new realization is independent of z 0 (t).
This process can also be used to deform a random variable. The gradual deformation of non-Gaussian random variable is achieved through a Gaussian transformation. For example, to deform a uniform random variable U with value in [0, 1], we consider the Gaussian random variable Z = G -1 (U), where G is the standard normal c.d.f. So, starting from an initial value u 0 , we derive z 0 = G -1 (u 0 ) deform it into z(t), and take:

Gradual Deformation of a Poisson Point Process
The presented methodology is to modify the fault positions by deforming the random numbers used to compute these fault positions. Therefore, simple non stationary point process sampling algorithms based on acceptation and reject methods cannot be used. The sequential method is used instead.

The Sequential Method to Draw Poisson Points
We will center the sub-seismic faults on points of a Poisson point process with a highly variable intensity. This is usually done by coupling the simulation of a Poisson point process with constant intensity with an acceptation-rejection technique. Such a method, however, requires a variable number of random numbers to produce a single Poisson point and is therefore not suitable to the gradual deformation method. Therefore, we will simulate the Poisson points with a sequential simulation method [29]. Let f(x i ,y j ) be the bivariate density of the nonstationary Poisson point process, discretized on an M × N grid. The marginal distribution of x i is given by: A coordinate x i is drawn by inverting the cumulative distribution function for a number u selected from a uniform distribution in [0, 1]: (11) Once x i is known, the conditional distribution f x i (y) is computed: (12) The y i coordinate is computed in a similar manner from a uniform number v: (13) Instead of modifying a single point, the same parameter t can be used to modify the coordinates of all Poisson points. As a consequence of the gradual deformation method, the spatial variations of the Poisson density are preserved, and therefore, the imposed scaling density pattern and the twopoint correlation dimension are also preserved. Another important aspect is that the number of parameters needed to deform the point process does not depend on the number of fractures but rather on the number of random realizations that are combined using Equation (8).
When the density map is homogeneous, the gradual deformation applied to Poisson points simulated with the sequential method yields smooth and continuous trajectories. Multifractal density maps generated with the cascade algorithm have a very heterogeneous nature and this approach is not efficient in that case: deformation trajectories are noisy and points tend to jump from one position to another in the y direction (Fig. 2). Indeed, when the gradual deformation method is applied, the abscissa changes from x 1 → x 2 . The cumulated distribution function in the y direction changes The difference between these functions in the heterogeneous case can be quite important. As a consequence, even small variations of u lead to large variations of y (Fig. 3).
The continuity of the gradual deformation process is a crucial point in the success of the optimization process. In order to reduce the number of evaluations of the Objective Function the Gauss-Newton algorithm is used to perform the optimization. In this algorithm, search directions for optimal parameters are deduced from the gradient and the hessian of the Objective Function. These are computed using small perturbations of the input parameters. Therefore, the Objective Function must be continuous and 2 times derivable. This condition cannot be satisfied if fault trajectories are discontinuous because they would lead to sudden modifications of the hydrodynamic response.
To avoid these discontinuities, a smoother density map may be used. At the last step of the multiplicative cascade algorithm and before the gradual deformation, interpolations may be used throughout the iterations to avoid squared features (Fig. 1c). Instead of multiplying new weights with the parent cells, they are multiplied by a bicubic interpolation of the parent map. Figure 4 shows the difference between interpolated and non-interpolated cascades.
This interpolation does not grant the exact conservation of the correlation dimension. The interpolations have a smoothing effect on the final map. Therefore the real correlation dimension tends to be higher than the targeted dimension D 2 used to build the multiplicative cascade. The exact correlation between the interpolation and the final dimension is difficult to establish and has to be verified.  Differences between F x1 -1 (u) and F x2 -1 (u) when x changes from x 1 to x 2 yielding a jump in the vertical direction. Non interpolated a) and interpolated b) multiplicative cascades.
hence be deformed by inversion of the cumulated distribution function: where l i (t) is the length of fault i, P the power law c.d.f. of fault length, and u i (t) a uniform number.

Geological Context
The Bloemendaal reservoir has no reality. Nevertheless, its geological data are based on studied carbonate reservoirs and the associated geological model was built using geological analysis, statistical property generations as described below.
The Bloemendaal field is a 12 × 15 km oil saturated carbonate reservoir. The reservoir has the structure of an anticline with North/South orientation. Mean thickness is equal to 10 m.
The geology of the study area can be described with a single reservoir unit and three facies. Carbonate lithotypes with good reservoir qualities (in terms of porosity and permeability) are merged into the LT1 facies. LT2 represents dolomitized lithotypes with secondary porosity attributed to dissolution of macrofossils in carbonate rocks. Third facies LT3 represents shales and other poor lithotypes.
Past tectonic events have created faults with a N110E mean orientation. 106 seismic sub-vertical faults are identified. Average reservoir pressure is equal to 150 bars.
Statistic distribution of porosity was established from well logs and a correlation between porosity and permeability has been established as log(K X ) = A.ω + B.
A facies realization is generated using a sequential indicator function algorithm, [32], with an exponential variogram. Porosity realizations are generated for each facies with a Gaussian variogram. Mean porosities in all three lithotypes range from 15 to 25% and permeabilities from 5 to 10 mD. The matrix is supposed to be isotropic (K X = K Y ) and vertical permeability is given as a function of the horizontal one K Z = 0.2K X . The matrix is water wet with capillary pressures ranging from 0.35 to 1.38. Relative permeabilities range from 0 to 1 (Fig. 6). Simple crossed permeability curves are assigned to the fractured medium and capillary pressures are neglected.
The production is ensured by a water injection strategy. Four production wells are located on top of the anticline. Eleven injectors build up pressure around the trap (Fig. 7). The injection rate of each well is set to 2 500 m 3 /day. Producers are set to a constant bottom hole pressure of 60 bar. Deformation trajectories with the multiscale sequential method.
Nevertheless, the seismic fault correlation is estimated on field data and a range of correlation dimension values may be established. A variation of the generated correlation dimension may be assumed. In addition, using interpolations the final density map seems to be more natural from a geological point of view. Future work would require to take the smoothing factor into account while modeling the faults or to use a method that directly builds continuous correlated density maps [31]. In this method, the continuity of fault displacements is enhanced using smoother density maps.

The Multiscale Sequential Method
The multiscale method proposes to modify the sequential method in order to reduce discontinuities during gradual deformations. First, a low resolution map is computed by averaging densities inside larger cells. This map is generally equal to the initial density map D 0 . Then for each cell the number of points N loc to draw is computed as a function of the total number of points N tot : (14) Then for each cell of D 0 the N loc points are drawn. With this method deformation trajectories are restricted to the cells of D 0 . Vertical displacement is hence limited and trajectories are more continuous (Fig. 5).

Gradual Deformation of Fault Lengths
Fault lengths are sampled from a power law distribution using the inverse transform sampling method. They can

Fault Characterization
Because of the limited resolution range of seismic methods, all Bloemendaal faults are not identified. The seismic fault network is the only available data. A geostatistical model is proposed to generate the missing objects. In this part, a two-step method is proposed to analyze the properties of the seismic faults. First, length distribution is studied. Then, a multifractal method is proposed to analyze and to constrain the sub-seismic fault model to the visible seismic ones. The fault length distribution, [14], is modeled using a power law distributions; n(l)∝l -a where n(l) is the number of faults shorter than l and a is the power law exponent. This choice is justified from outcrop observations and geomechanical models. Power laws allow one to extend the seismic fault distribution to smaller sub-seismic faults. The power law exponent is determined by plotting the cumulated length distribution on a log-log scale. The plot is linear in a given length range [l 0 , l max ]. The slope yields the exponent "-a".
The classical censoring and truncation effect, [33], can be observed on the measured fault length distribution. The deviation at larger fault lengths is due to edge effects: larger faults are incompletely sampled and appear smaller than reality. Deviation for smaller fault lengths is due to the resolution limit of seismic methods. The number of sub-seismic faults can be derived using an extension of the power law and is a function of the minimum simulated fault length l min as shown in Figure 8. This length was estimated at 200 m on a natural analogs.
To constrain the sub-seismic fault center distribution to the seismic faults, the fault centers are supposed to be a scaled structure with a constant correlation dimension. It is the major assumption of the proposed methodology. The analysis of the Capilary pressures and relative permeabilities for the initial model.

Figure 7
Bloemendaal reservoir structure: 11 injection wells and 4 production wells. Seismic faults are represented in black. seismic fault network starts with the computation of these centers. The multifractal box counting algorithm is applied to these centers to compute the 2nd-order moment M 2 . The measured quantity F i is equal to the number of centers per grid block. The 2nd-order moment M 2 is then plotted on a log-log scale. The plot is linear within a limited observation range r min < r < r max . D 2 is the slope of the fitted line. Deviation at box sizes larger than r max is attributed to the finite size of the system. Deviation for box sizes smaller than r min is due to the underestimation of smaller faults. The fault center are hence fractal within the observation scales [r min , r max ].
The seismic fault length distribution and multifractal parameters are analyzed (Fig. 8, 9). About 1 060 faults are below seismic resolution. Minimum subseismic fault length is 200 m and maximum length 1 000 m with a power-law exponent of 1.63. D 2 is set to 1.7 and the initial resolution r 0 of the multiplicative cascade is 2 000 m.

Subseismic Fault Model
The following step describes the method for generating 1 060 subseismic faults while preserving the seismic fault structure. The algorithm is initiated with the computation of the initial density map (Fig. 10a). The multiplicative cascade algorithm is applied to refine the density map (Fig. 10b). Fault centers are drawn using the multiscale sequential algorithm (Fig. 10c).
Fault lengths are drawn and the initial geological model is built (Fig. 11). Fault conductivities and apertures are respectively set to 3 000 mD.m and 0.01 m.

History Matching
The equivalent fault permeability and porosity are computed on a 121 × 150 × 1 dual-medium grid. Matrix permeability is set to 50 mD in the X and Y directions. Porosity is set to 0.3. Matrix-fracture exchange terms are defined through block sizes computed according to the fault geometries. A twophases (water-oil) flow simulation is run at the field scale using matrix and fracture relative permeabilities and capillary pressures (previously described). Water-cut levels are considered because they are representative of the connectivity of the network. As expected, simulated production data differs from measured data. The irregular shapes of the curve are the consequence of the presence of different injectors. Indeed, each injector contributes to the water breakthrough of the producers Multifractal box-counting algorithm applied to the fault center distribution.

Figure 10
The seismic fault network is used to calculate an initial low resolution density map a) which is used to initiate the multiplicative cascade algorithm b), the final density map is used to draw a population of fault centers c). but the arrival times are different. Early water breakthroughs can be observed on wells 1, 2 and 3, (Fig. 12). The position of the sub-seismic faults increases the water flow between the wells. Some sub-seismic faults have to be moved away while respecting the statistical fault model. History matching is hence necessary. All the wells cannot be matched simultaneously when the history matching method is applied to the whole fractured domain. If a single parameter is used to control the positions of all faults at the same time, it is highly probable that the match to dynamic data will be improved for some wells and will worsen for some others. It is better to have independent control on the match to each well. The reservoir model is therefore divided in four zones defined according to groups of injectors and producers (Fig. 13). One gradual deformation parameter is assigned to modify the fault positions in each zone. Fault lengths remain constant. Centers are displaced independently in each zone. Interactions can occur between regions because the lineament of a given fault can cross a lineament of another zone. The initial geological model can then be deformed independently in each zone to increase the Objective Function convergence rate. Fault positions are optimized using the Sequential Quadratic Programming and Augmented Lagrangian (SQPAL) gradient-based algorithm Seismic fault (black lines), initial sub-seismic fault realization (gray lines) and elliptic study zones (gray).

Figure 12
Initial water cut levels for P1, P2, P3 and P4. Field production data is in red, and simulated values in yellow. [34,35]. When an optimal realization has been found, the match is further increased by combining it with a new realization. Figure 14 shows the evolution of the Objective Function during the gradual optimization process. The match to production data is satisfying as shown in Figure 15. Figure 16 shows the optimal fault realization. The water breakthrough time are correctly predicted. Comparison of the gray elliptic zones of Figures 11 and 16 shows that the number of faults has been locally lowered. In comparison with the initial model the water flow takes more time to reach P1 and P3 in the optimized model. Nevertheless, the last periods of the watercut curves of wells P2 and P4 are not correctly predicted. As the water breakthrough time is correctly estimated, we can assume we correctly model the flow between one injector and the associated producer but not the influence of the other injectors. More work is needed to perfectly match these two wells.

CONCLUSION
In this paper, we presented a methodology to perform the history matching of sub-seismic fault models. We first presented a stochastic fault generator that was specifically developed for history matching purposes. We then proposed Division of the reservoir into 4 independent optimization zones. The different colors are used to represent the different deformation zones.

Figure 14
Evolution of the Objective Function during the history matching process.
an original history matching method in which the faults are displaced to match the production data. The main conclusions of this work are that: -the proposed stochastic model uses a reduced number of parameters to describe the geometry of the realizations (fractal dimension, length distribution). These parameters are easily inferred from an analysis of the seismic fault network; -the fault generator is based on fractal geometry, which has been widely used in petroleum applications to describe and model fractured networks. Realizations are realistic and constrained to the available seismic faults; -the fault generator is compatible with the gradual deformation method: all fault positions depend on a reduced number of gradual deformation parameters; -the gradual deformation preserves the statistical coherency of the model: length distribution and scaling properties are left unchanged during the optimization; -the use of the multiscale sequential algorithm provides smooth and continuous deformation trajectories. As a consequence, hydrodynamic behavior is continuous with respect to gradual deformation parameters; -the reduced number of deformation parameters and the continuity of the OF enables the use of gradient-based Optimal water cut levels for P1, P2, P3 and P4. Field production data is in red, and simulated values in yellow.
algorithms. This reduces the number of fluid flow simulations used to evaluate the Objective Functions; -the proposed methodology proves to be successful on a semi-synthetic case and should be tested on real case studies. Future work will focus on: -testing global optimization methods: indeed, some fault configurations are difficult to optimize because the Objective Function contains local minima; -take the presence of sealing faults into account. It is a challenging task because compartmentalization effects yield non linear hydrodynamic behavior when the faults are displaced.