# Identification of Geological Forms Using Production Data

M. Masmoudi, P. F. Edoa, D. Rahon

<figure><img src="image_1.png" /></figure>

## To cite this version:

M. Masmoudi, P. F. Edoa, D. Rahon. Identification of Geological Forms Using Production Data. Oil & Gas Science and Technology - Revue d'IFP Energies nouvelles, 1999, 54 (1), pp.5-21. 10.2516/ogst:1999001. hal-02075776

HAL Id: hal-02075776

https://ifp.hal.science/hal-02075776

Submitted on 21 Mar 2019

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
# Identification of Geological Forms Using Production Data

P. F. Edoa¹, D. Rahon¹, M. Masmoudi²

1 Institut français du pétrole, Hélioparc Pau-Pyrénées 2, avenue Pierre Angot, 64053 Pau Cedex 9 - France

2 Centre européen de recherche et de formation avancée en calcul scientifique, 42, avenue Gaspard Coriolis, 31057 Toulouse Cedex - France

**Résumé** — *Identification des formes géologiques en utilisant les données de production* — L'identification des formes géologiques de chenaux, de failles ou de limites de réservoirs à partir des données de tests de puits ou d'historiques de production constitue un problème délicat en ingénierie de réservoirs pétroliers. Des modèles analytiques d'interprétation de tests de puits sont utilisés dans des cas simples. Dans les cas multiphasiques ou de géométries complexes, il est nécessaire de recourir à des modèles de simulation d'écoulement de fluides.

Dans cet article, nous nous proposons de résoudre un problème d'inversion de forme associé à un modèle d'écoulement diphasique eau-huile et notre objectif est d'identifier la forme et la position de corps géologiques dans un réservoir d'hydrocarbures par calage des données de production.

En identification de forme, deux problèmes se posent généralement : le choix de la représentation des formes géométriques et le calcul des sensibilités. Pour calculer les gradients, des techniques de perturbation de maillages existent dans des domaines où les éléments finis sont d'utilisation courante. Dans un contexte de maillages strukturés fréquemment utilisés en simulation de réservoirs, ces techniques de perturbation de maillages sont difficilement applicables.

La méthode proposée est basée sur le calcul des sensibilités sur le problème continu par rapport à la géométrie des différents corps géologiques. Cette géométrie est définie par une triangulation. Les paramètres de calage sont les nœuds de la triangulation, et les gradients de la fonction objectif sont calculés par rapport aux déplacements de ces nœuds.

Un algorithme d'optimisation couplé avec un simulateur d'écoulements polyphasiques a été développé. Il permet la prise en compte de contraintes géométriques sur les paramètres de calage et assure la régularité des formes obtenues. Plusieurs applications de type académique ont été effectuées. Ce qui conduit, à l'aide de calages des données de production, à une meilleure caractérisation de la forme, de la taille et de la position des corps sémentaires, plus particulièrement des limites de réservoirs, de la position et de la taille des failles ainsi que de l'épaisseur et de la largeur des chenaux.

Mots-clés : écoulements multiphasiques, historique de production, inversion de forme, problème adjoint, gradient, optimisation.

**Abstract** — *Identification of Geological Forms Using Production Data* — Identifying the geological forms of channels, faults and boundaries of reservoirs on the basis of well-tests or production history is a tricky problem for oil reservoir engineering. Analytical interpretation models of well-tests are used in simple cases. In multi-phase cases or those which present a certain degree of geometrical complexity, it becomes necessary to use fluid flow simulation models.

In this article, we intend to solve a problem of form inversion associated with a two-phase oil-water flow model in which the aim is to identify the form and the position of geological bodies in a hydrocarbon reservoir with production data match.

In identifying the form, two problems generally occur: the choice of representation of geometrical forms and the calculation of sensitivities. To calculate the gradients, mesh perturbation techniques exist in fields where finite elements are commonly used. In the context of structured meshes often used in reservoir simulation, the application of such mesh perturbation techniques is difficult.
The method suggested is based on the calculation of sensitivities on the continuous problem in relation to the geometry of the various geological bodies. Such a geometry is defined through triangulation. The adjustment parameters are the triangulation nodes and the gradients of the objective function are calculated in relation to the displacement of such nodes.

An optimization algorithm coupled with a polyphase flow simulator has been developed. It takes into account the geometrical constraints of the adjustment parameters and guarantees the regularity of the forms obtained. Several laboratory applications have been carried out. With the help of calibrated production data this leads to a better characterization of the form, the size and the position of the sedimentary bodies in particular: reservoir boundaries, position and size of the faults and thickness and width of the channels.

Keywords: multiphase flows, production history, shape inversion, adjoint problem, gradient, optimization.

## NOTATIONS

<table><tr><td>P</td><td>pressure in the reservoir</td></tr><tr><td>P<sub>0</sub></td><td>initial reservoir pressure</td></tr><tr><td>S</td><td>water saturation in the reservoir</td></tr><tr><td>S<sub>0</sub></td><td>initial water saturation</td></tr><tr><td>u,v</td><td>solutions of the adjoint problem</td></tr><tr><td>K</td><td>permeability tensor</td></tr><tr><td>K<sub>h</sub></td><td>horizontal permeability</td></tr><tr><td>K<sub>z</sub></td><td>vertical permeability</td></tr><tr><td>φ</td><td>porosity in the reservoir</td></tr><tr><td>φ<sup>0</sup></td><td>porosity at initial pressure</td></tr><tr><td>c<sub>R</sub></td><td>rock reservoir compressibility</td></tr><tr><td>μ</td><td>viscosity</td></tr><tr><td>kr</td><td>relative permeability</td></tr><tr><td>ρ</td><td>density</td></tr><tr><td>q</td><td>injected mass volume flow or mass volume flow produced at well</td></tr><tr><td>t</td><td>time variable</td></tr><tr><td>(x, y, z)</td><td>point in three dimensional space</td></tr><tr><td>Ω</td><td>spatial field representing the reservoir</td></tr><tr><td>Ω<sub>m</sub></td><td>region of measurements in Ω</td></tr><tr><td>χ<sub>Ωm</sub></td><td>characteristic function on Ω<sub>m</sub></td></tr><tr><td>Γ</td><td>forms of geological bodies</td></tr><tr><td>V</td><td>domain perturbation</td></tr><tr><td><sup>1</sup>DV</td><td>DV transposed</td></tr><tr><td>j</td><td>objective function</td></tr><tr><td>N</td><td>number of geometrical parameters</td></tr><tr><td>NP</td><td>number of production wells in the reservoir</td></tr><tr><td>l</td><td>production well index</td></tr><tr><td>NT</td><td>number of measurements made on wells</td></tr><tr><td>n</td><td>index of an instant of measurements</td></tr><tr><td>d<sub>Γ</sub>j(Γ)·V</td><td>derivative of j in relation to Γ in the direction V</td></tr><tr><td>div</td><td>operator divergence</td></tr><tr><td>grad</td><td>operator gradient</td></tr><tr><td>a·b</td><td>Euclidian scalar product of vectors a and b</td></tr><tr><td>L<sup>2</sup>(Ω<sub>m</sub>)</td><td>space of functions that squares can be integrated on Ω<sub>m</sub></td></tr></table>

## INDICES

w water phase index

o oil phase index

i geometrical node index.

## INTRODUCTION

In order to characterize hydrocarbon reservoirs, calibration techniques of well-tests and production history have been developed over the last few years. Most of these are based on effective gradient methods for solving a problem of inversion. The problem of solving an inversion boils down to the minimization of an objective function and leads to history matching. In such methods, numerical simulations are very often used.

Most of the procedures merely enable identification of petrophysical parameters or characterization of wells. The case of geometrical parameters such as the position of a fault, reservoir boundaries, breadth, thickness or channel forms has been treated in recent publications [1]. Such studies have been concentrated on the definition of geological bodies using mathematical functions describing the permeability and the porosity of the rock reservoir. The sensitivities of the objective function in relation to the petrophysical parameters (permeability, porosity) or geometrical parameters are calculated through derivation of direct state discrete equations [2-6] or through a so-called adjoint state method [7]. The traditional gradient methods can then be used, but generalizing such solutions to complex geological models seems difficult.

The aim of research in the work presented here is to develop a shape inversion method which enables the identification of geometrical parameters which describe geological bodies using numerical simulation flows and match of production data. For this, two separate meshes are used, one supposedly fixed for simulating flows, and the other subject to modifications to represent the geological forms. The approach described in this article has been first applied to single phase flows and then adapted to well-test interpretation [8]. The encouraging results obtained have
contributed to extending its application to polyphase flows to match of production history.

After describing the method, some examples of applications in two-phase water-oil reservoirs are presented. The first two models are synthetic cases. The third is a more realistic, heterogeneous geological model comprising several geological layers and a significant group of faults. The efficiency of the method is clearly illustrated through these examples. Finally, a description of sensitivity calculations is given in the Appendix.

# 1 PROPOSED METHOD

The traditional interpretation of wells production data on the basis of analytical solutions only allows the identification of average values of parameters characterizing simple geometry hydrocarbon reservoirs. For more complex geological models and heterogeneous media, it becomes necessary to use numerical flow simulation models.

In order to achieve efficiency, the optimization or shape inversion process requires that the gradients of an objective function be calculated in relation to the controlled geometrical parameters. There are at least two ways of calculating such gradients.

The first consists of using a finite element formulation and non structured meshes to represent complex geological bodies. In such a case, the response sensitivities to the wells, evaluated through the objective function, in relation to the perturbations of mesh nodes can be calculated by deriving discrete equations. Such a solution requires the use of meshes that can be perturbed, such as the finite elements for simulating flows in the reservoir [9, 10]. However, such meshes are not often used in the petroleum sector.

With the classic fluid flow simulation using regular meshes, the methods of gradient calculation of the objective function on discrete equations are difficult to apply, particularly when the control parameters are not directly involved in the equations of the model. Moreover, the discrete problem cannot be differentiated in relation to perturbations of geological forms when the interface between two geological bodies crosses a mesh.

How does one proceed with structured meshes? The solution proposed is to determine the expression of the sensitivities of the objective function on the continuous problem using mathematical theory and then calculate gradients through the discrete form of the expression obtained.

The main stages of the method may thus be summed up as under:

- determination of sensitivities of the objective function on the continuous model on the basis of flow equations and so called adjoint state equations;

- choice of the representation of geological forms and gradient calculation;

- and lastly, use of gradients in an appropriate optimization algorithm.

# 2 DERIVATIVE IN RELATION TO THE GEOMETRY

All the calculations that follow are formal and the authors will not try to justify them.

## 2.1 Direct State Equations

Let $\Omega$ be a spatial field of two or three dimensions representing a hydrocarbon reservoir of regular boundary $\Sigma$ and $(0,T)$ a time interval covering the production history period. The evolution of the pressure $P$ of the fluids and the saturation $S$ in water corresponding to an oil water two-phase flow without capillary pressure in the reservoir is governed by the partial derivatives equations:

$$ \frac{\partial}{\partial t}(\phi(P)S) - \text{div} \left( K \frac{kr(S)}{\mu_w} \text{grad}(P + \rho_w gz) \right) - \frac{q_w}{\rho_w} = 0 \quad (1) $$
$(x, y, z) \in \Omega, t \in [0, T]$

and:

$$ \frac{\partial \phi(P)}{\partial t} - \text{div} \left( K \frac{kr(S)}{\mu_w} \text{grad}(P + \rho_w gz) \right) - \text{div} \left( K \frac{kr(S)}{\mu_o} \text{grad}(P + \rho_o gz) \right) - \frac{q_w}{\rho_w} - \frac{q_o}{\rho_o} = 0 \quad (x, y, z) \in \Omega, t \in [0, T] $$

with the boundary conditions:

$$ \left( K \frac{kr(S)}{\mu_w} \text{grad}(P + \rho_w gz) + K \frac{kr}{\mu_o} \text{grad}(P + \rho_o gz) \right) \cdot n = 0 \quad (3) $$
$(x, y, z) \in \Sigma, t \in [0, T]$

and the following initial conditions:

$$ P(x, y, z, 0) = P_0(x, y, z) \quad \text{and} \quad S(x, y, z, 0) = S_0(x, y, z) \quad (4) $$
$(x, y, z) \in \Omega$

The boundary conditions (3) express the fact that the external boundary of the reservoir $\Sigma$ is impermeable to the flow; this means that the total fluid flux across this boundary is nil in the course of time.

## 2.2 The Objective Function

The forms of the geological bodies contained in the reservoir $\Omega$ are noted $\Gamma$, and these constitute the unknowns of the
problem to be solved. The aim is to minimize a functional $j$ which depends on $\Gamma$ and whose expression is defined by:

$$j(\Gamma) = \frac{1}{2} \int_{0}^{T} \beta \| P_{\Gamma} - P^{m} \|_{L^{2}(\Omega_{m})}^{2} + \gamma \| q_{w}^{\Gamma}(S_{\Gamma}) - q_{w}^{m} \|_{L^{2}(\Omega_{m})}^{2} \, dt \\ + \frac{1}{2} \int_{0}^{T} \sigma \| q_{o}^{\Gamma}(S_{\Gamma}) - q_{o}^{m} \|_{L^{2}(\Omega_{m})}^{2} \, dt \quad (5)$$

$P^m$, $q_w^m$ and $q_o^m$ represent respectively the pressure measurements as well as those of water and oil rates at the production wells. $P_\Gamma$, $S_\Gamma$, $q_\Gamma^\Gamma$ and $q_o^\Gamma$ represent respectively the pressure, the water saturation and the rates of water and oil produced during the time interval $(0, T)$ corresponding to the position $\Gamma$. $\Omega_m$ is a part of $\Omega$ that is identified with the reunion of the establishment zones of production wells in $\Omega$. It is here that the measurements are carried out. $\beta$, $\gamma$ and $\sigma$ are positive constants which serve as the weighting coefficients. $\beta$, $\gamma$ and $\sigma$ are determined in such a way that the terms evaluated in the cost function are homogeneous.

Optimal shape depends exclusively on the criterion chosen. The objective function chosen enables taking into account certain physical phenomena at the production wells (arrival of water, fall in pressure, closing of a well, etc.).

## 2.3 Derivative of the Objective Function

By using the transport technique and the Lagrangian method introduced in particular by Céa [11], it is possible to show that the derivative of the objective function $j$ in relation to $\Gamma$ is given by:

$$\begin{aligned} d_{\Gamma} j(\Gamma) \cdot V = & \int_{\Omega_{\Gamma}} \left( \left( A \frac{\partial S_{\Gamma}}{\partial t} + B S_{\Gamma} \frac{\partial P_{\Gamma}}{\partial t} + B P_{\Gamma} \frac{\partial S_{\Gamma}}{\partial t} \right) u_{\Gamma} \right) \mathrm{div} V \\ & + \int_{\Omega_{\Gamma}} \left( B \frac{\partial P_{\Gamma}}{\partial t} v_{\Gamma} \right) \mathrm{div} V \\ & + \int_{\Omega_{\Gamma}} \frac{kr}{\mu_{w}} K \mathrm{grad} (P_{\Gamma} + \rho_{w} gz) \cdot \mathrm{grad} u_{\Gamma} \mathrm{div} V \\ & - \int_{\Omega_{\Gamma}} \frac{kr}{\mu_{w}} (K^{t} DV + DV K) \mathrm{grad} (P_{\Gamma} + \rho_{w} gz) \cdot \mathrm{grad} u_{\Gamma} \\ & + \int_{\Omega_{\Gamma}} \frac{kr}{\mu_{w}} K \mathrm{grad} (P_{\Gamma} + \rho_{w} gz) \cdot \mathrm{grad} v_{\Gamma} \mathrm{div} V \\ & + \int_{\Omega_{\Gamma}} \left( \frac{kr}{\mu_{o}} K \mathrm{grad} (P_{\Gamma} + \rho_{o} gz) \right) \cdot \mathrm{grad} v_{\Gamma} \mathrm{div} V \\ & - \int_{\Omega_{\Gamma}} \frac{kr}{\mu_{w}} (K^{t} DV + DV K) \mathrm{grad} (P_{\Gamma} + \rho_{w} gz) \cdot \mathrm{grad} v_{\Gamma} \\ & - \int_{\Omega_{\Gamma}} \frac{kr}{\mu_{o}} (K^{t} DV + DV K) \mathrm{grad} (P_{\Gamma} + \rho_{o} gz) \cdot \mathrm{grad} v_{\Gamma} \end{aligned} \quad (6)$$

where:

$$- \Omega_{T} = \Omega \times (0, T);$$

$$- A = \phi^{0} (1 - c_{R} P_{0}) \text{ et } B = \phi_{0} c_{R}; c_{R} \text{ is the rock compressibility and } \phi_{0} \text{ the porosity of the rock at reference pressure } P_{0};$$

$$- V \text{ is a function called perturbation of the field } \Omega, \text{ in fact of } \Gamma; \text{ this means that } V(\Gamma) \text{ is used to determine a new position of } \Gamma;$$

$$- P_{\Gamma} \text{ and } S_{\Gamma} \text{ are the pressure and the water saturation in the field corresponding to the position of the interface } \Gamma;$$

$$- u_{\Gamma} \text{ and } v_{\Gamma} \text{ are called the adjoint state [7, 12]. These are the solutions to the evolution problem defined by the following partial derivatives equations:}$$

$$\begin{aligned} & - S\phi'(P) \frac{\partial u}{\partial t} - \phi'(P) \frac{\partial v}{\partial t} \\ & - \mathrm{div} \left( \left( \frac{kr(S)}{\mu_w} + \frac{kr(S)}{\mu_o} \right) K \mathrm{grad} v \right) \\ & - \mathrm{div} \left( K \frac{kr(S)}{\mu_w} \mathrm{grad} u \right) + \beta(P_{\Gamma} - P^m) \chi_{\Omega_m} = 0 \end{aligned} \quad (7)$$

$$(x, y, z) \in \Omega, \ t \in [0, T]$$

$$\begin{aligned} & - \phi(P) \frac{\partial u}{\partial t} - \frac{1}{\rho_w} \frac{\partial q_w}{\partial S} u - \left( \frac{1}{\rho_w} \frac{\partial q_w}{\partial S} + \frac{1}{\rho_o} \frac{\partial q_o}{\partial S} \right) v \\ & + \left( \frac{K}{\mu_w} \frac{\partial kr(S)}{\partial S} \mathrm{grad} (P + \rho_w gz) \right) \cdot \mathrm{grad} u \\ & + \left( \frac{K}{\mu_w} \frac{\partial kr(S)}{\partial S} \mathrm{grad} (P + \rho_w gz) \right) \cdot \mathrm{grad} v \\ & + \left( \frac{K}{\mu_o} \frac{\partial kr(S)}{\partial S} \mathrm{grad} (P + \rho_o gz) \right) \cdot \mathrm{grad} v \\ & + \left( \gamma (q_w^\Gamma(S_\Gamma) - q_w^m) \frac{\partial q_w^\Gamma(S_\Gamma)}{\partial S} \right) \chi_{\Omega_m} \\ & + \left( \sigma (q_o^\Gamma(S_\Gamma) - q_o^m) \frac{\partial q_o^\Gamma(S_\Gamma)}{\partial S} \right) \chi_{\Omega_m} = 0 \end{aligned} \quad (8)$$

$$(x, y, z) \in \Omega, \ t \in [0, T]$$

verifying the boundary conditions:

$$u = 0; \quad (x, y, z) \in \Sigma, \ t \in [0, T] \quad (9)$$

and:

$$\left( K \frac{kr(S)}{\mu_w} \mathrm{grad} u \right) \cdot n + \left( \left( \frac{kr(S)}{\mu_w} + \frac{kr(S)}{\mu_o} \right) K \mathrm{grad} v \right) \cdot n = 0 \quad (10)$$

$$(x, y, z) \in \Sigma, \ t \in [0, T]$$
and at the final moment T the condition here below:

$$u(x, y, z, T) = 0 \text{ and } v(x, y, z, T) = 0; (x, y, z) \in \Omega \quad (11)$$

The integration of this adjoint problem is done retrogressively. It may be noted that the direct state equations system describing the two-phase water oil flow comprises two non linear partial derivatives Equations (1) and (2) of the unknowns P and S. The partial derivatives Equations (7) and (8) of the unknowns u and v relative to the adjoint problem are linear. A complete description of the calculation of the derivative of j is given in the Appendix.

# 3 CALCULATION OF GRADIENTS

## 3.1 Representation of Geological Forms

Once the expression of the derivative of the objective function on the continuous problem has been determined (6), the calculation of the gradients can now be carried out.

As specified earlier, the main aim is to develop a shape inversion technique enabling the characterization of oil reservoirs using classic fluid flow simulation. Consequently, the first thing to take into account is to consider a representation of forms of geological bodies which are independent of the fixed Cartesian mesh used for flow calculations.

The form Γ of a geological body is represented by a finite elements mesh with either two or three dimensions. The control parameters are the Cartesian coordinates of the nodes of such triangulation. Each geological body is associated with a set of petro physical parameters: a permeability tensor and a porosity value or a porosity multiplier coefficient. The permeability and porosity maps pertaining to the Cartesian meshes of the reservoir are then generated by determining the meshes located within the body and through homogenization of the meshes cut by the surface of the body. In case of the presence of faults, a change of linkage transmissibility in all the meshes cut by a fault may be carried out.

## 3.2 Choice of Perturbation V

The following step consists of defining a perturbation V of the domain Ω whose role is to deform the geometry. Its support is limited to the neighbourhood of Γ. The control parameters to be identified being the positions on the Cartesian coordinates of the triangulation nodes of Γ, for each index node i, a local perturbation V_i is defined in the following way:

- the support of V_i is determined by the triangulation elements of Γ with a common node i;

- V_i varies from 1 at the index node i to 0 in its support and is nil everywhere else.

The perturbation V is thus equal to:

$$V = \sum_{i=1}^{N} \alpha_{i} V_{i} \quad (12)$$

the α_i are the coefficients to be determined and N the total number of nodes or control parameters. Searching for the values of α_i, is then done through an appropriate optimization algorithm, for example the algorithm of the conjugated gradient of Polak-Ribière [13] with constraints.

The discrete form of the objective function j is given by:

$$j(\Gamma) = \frac{1}{2} \sum_{l=1}^{NP} \sum_{n=1}^{NT} \beta \left( P_{\Gamma}(W_l, t_n) - P^m(W_l, t_n) \right)^2 \\ + \frac{1}{2} \sum_{l=1}^{NP} \sum_{n=1}^{NT} \gamma \left( q_w^{\Gamma}(W_l, t_n) - q_w^m(W_l, t_n) \right)^2 \\ + \frac{1}{2} \sum_{l=1}^{NP} \sum_{n=1}^{NT} \sigma \left( q_o^{\Gamma}(W_l, t_n) - q_o^m(W_l, t_n) \right)^2 \quad (13)$$

where NP is the number of production wells in the reservoir, W_l the position of a production well in Ω, NT the number of measurements carried out at wells and t_n a moment of such measurements.

The expression (6) then allows to calculate for each index i the value of d_Γj(Γ)·V_i, which determine the gradient of j. The functional j not being a priori convex, it would be interesting to look for local minima of j.

# 4 THE OPTIMIZATION PROCESS

After the calculation of the gradients of the objective function in relation to the perturbations of the nodes representing geological forms, such gradients can be introduced in a minimization process.

Most of the time, the identification of a geometrical form brings into play a set of constraints that can be either geometrical or esthetic. Such stresses are often written as inequalities limiting the control parameter values. Thus the problem of real optimization boils down to looking for an optimal form Γ which minimizes the objective function j while respecting the given geometrical constraints. These are defined, for example, in accordance with the dimensions of the reservoir and the data pertaining to well positions.

The process of optimization developed is summarized schematically in the following paragraphs. Each time the optimization process is repeated, the descent of the objective function j is obtained in a direction projected on the constraints space. Such a projection ensures at each step the generation of a new set of parameters respecting the stresses and is very efficient in the case of linear constraints [14, 15].
The optimization process schema used to change the geometry of the model is thus the following:

- Construction of a structured (or non structured) gridded model of the reservoir.

- Representation of the geometry of geological bodies through triangulation defining a 2D or 3D surface.

- Attribution of petro physical data ($\phi$, $K$) to each mesh.

- Attribution of petro physical parameters to volume bodies and modification of transmissivities for faults.

- Beginning of the optimization loop;

• updating of maps of $\phi$ and $K$ by the projection of volume bodies or modification of transmissivities for the faults;

• resolution of the direct problem (1)-(4): flow simulation, calculation of pressure and saturation maps;

• resolution of the adjoint problem (7)-(11);

• calculation of gradient for perturbations of mesh nodes of geological bodies using direct and adjoint state solutions (see Section 3);

• determination of the direction of displacement of nodes with the help of a gradient algorithm: conjugated gradient algorithm of Polak-Ribière with constraints;

• modification of the geometry of the model by displacement of nodes: new positions of geological bodies;

• test of stoppage criteria and verification of conditions of optimality.

- End of the optimization loop.

# 5 NUMERICAL RESULTS

Some results obtained with the method are presented in this section. They are purely academic, but are repre-sentative of certain simple realistic models; real cases still being extremely difficult to study.

The first two examples are synthetic models and the third is the result of a geological model which uses static data (geometrical and petrophysical). The history production was generated using a numerical simulator. In these examples the objective is to identify the form and the position of the geological bodies (faults, channels) with the help of digital simulations and by production history matching. For this, the optimization algorithm used (the projected gradient method) is coupled with a numerical multiphase flow simulator Sarip [16].

## 5.1 Elaborated Schema for the Different Tests

Starting from an initial position of the nodes describing the discrete form of $\Gamma$. The minimization process determines at each iteration a new position of $\Gamma$ reducing the objective function while respecting the constraints. This process stops when one can no longer reduce the objective function or

when the Euclidean norm of the gradient is small, i.e. less than a given positive number (a local minimum is found).

The stages which have been followed in the different applications are:

- Construction of a mesh grid model containing a geological body in reference, to which we give a value of porosity and a map of permeability or multiplicative coefficients of transmissivity for the faults.

- Calculate the pressures and the rates of water and oil in the wells corresponding to this reference that are used as measurements. Only in the third example, there is in addition a random noise on the output data.

- Definition of an initial geometry by the modification of the shape of the geological body in reference; the launching of the optimization process so as to recover measurements by minimization of the cost function.

## 5.2 Example 1: Determination of the Position of a Fault

The first model shown is a single-layer closed 2D reservoir. It is a heterogeneous reservoir with dimensions 1250 m x 750 m, and having a watertight non permeable fault whose reference position is given in the Figure 1. The Cartesian mesh grid of the reservoir consists of 25 x 15 meshes of 50 m side and the roof of the reservoir is situated at a depth of 4700 m. Numerically, the introduction of the fault in the reservoir is done by multiplying with adequate coefficients (equal to 0 here) the transmissivities of the connection in all the meshes cut by the fault. The position of the fault in the reservoir is shown by a 2D line segment (or by a 3D plan) connecting 2 nodes whose Cartesian coordinates represent the parameters of the calibration. For well-tests, the required geometrical parameters of the fault are its distance from the well, its length and its inclination.

<figure><img src="image_3.png" /><figcaption>Figure 1</figcaption></figure>

Example 1: reference position of the fault.
There are two wells drilled into the reservoir, one for production and one for injection. Each well has a radius of 7.85 cm. A production history is simulated over a period of 2000 days. The cumulative rate of flow of oil and water is fixed at 300 m³/d in the production well. The flow of water injected into the reservoir is equal to 150 m³/d in the injection well. The measurements of the pressure and the rates of water and oil relative to the reference fault are taken every 20 days during a period of production spread out over 2000 days.

The map of absolute permeability of the reservoir is obtained by log-normal simulation. The variance is equal to 1 and the average permeability is 50 mD. The initial porosity in the reservoir is equal to 0.3. The general characteristics of the model are summarized in Table 1.

<caption>TABLE 1</caption>

<table><caption>Example 1: parameters of the fault model</caption><thead><tr><td>Parameters</td><td>Values</td></tr></thead><tbody><tr><td>Number of horizontal meshes</td><td>25 x 15</td></tr><tr><td>Horizontal dimensions</td><td>1250 m x 750 m</td></tr><tr><td>Thickness of the reservoir</td><td>10 m</td></tr><tr><td>Depth of the reservoir</td><td>4700 m</td></tr><tr><td>Reference depth</td><td>4700 m</td></tr><tr><td>Rock compressibility</td><td>5·e-4 bar<sup>-1</sup></td></tr><tr><td>Water compressibility</td><td>4·36e-4 bar<sup>-1</sup></td></tr><tr><td>Oil compressibility</td><td>1·e-4 bar<sup>-1</sup></td></tr><tr><td>Density of water</td><td>1 kg/m<sup>3</sup></td></tr><tr><td>Density of oil</td><td>0.625 kg/m<sup>3</sup></td></tr><tr><td>Viscosity of water</td><td>0.31 Pa·s</td></tr><tr><td>Viscosity of oil</td><td>0.34 Pa·s</td></tr><tr><td>Initial pressure</td><td>300 bar</td></tr><tr><td>Initial saturation of water</td><td>0.3</td></tr><tr><td>Number of injection wells</td><td>1</td></tr><tr><td>Injection rate (water)</td><td>150 m<sup>3</sup>/d</td></tr><tr><td>Number of production wells</td><td>1</td></tr><tr><td>Production rate (water and oil)</td><td>300 m<sup>3</sup>/d</td></tr></tbody></table>

## Results

The reference fault has a length of 360 m and is situated at a distance of 650 m from the production well and at 250 m from the injection well. The initial fault in the beginning of inversion is 350 m long. It is situated at 230 m from the production well and at 640 m from the injection well. The optimal position of the fault was found at 700 m from the production well and at 200 m of the injection well, and it measures 352 m.

Figure 2 represents from right to left, the successive positions of the fault during the process of optimization. Figure 3 shows the graph of water rates of the production well for the initial position of the fault and the setting obtained after 6 iterations.

<figure><img src="image_4.png" /><figcaption>Figure 2</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 3</figcaption></figure>

<caption>Example 1: water rates match at the production well.</caption>

## 5.3 Example 2: Identification of the Limits of a Channel

This model, like the earlier one, is a closed single-layer 2D reservoir with a thickness of 10 m. It is designed in the form of a channel oriented along the *Ox* axis and its reference position is indicated in Figure 4. The dimensions of this reservoir are 2500 m x 1500 m with a thickness of 10 m. The roof of the reservoir is located at a depth of 4700 m. The Cartesian mesh grid of the model consists of 25 x 15 meshes of side 100 m. The channel is described by 13 nodes (6 for one of the limits and 7 for the other) and their Cartesian coordinates represent the calibration parameters.

Eleven wells are drilled in the channel and each one has a radius of 7.85 cm. They are distributed as follows: 5 production wells (Prod. 1 to Prod. 5) and 6 injection wells (Inj. 1 to Inj. 6). Their positions are indicated in the Figure 4. The simulated production history consists in recovering oil and water from the production wells by injecting water
through the injection wells over a period of 1000 days. In each production well the cumulative rate of water and oil produced is fixed at 300 m³/d and the quantity of water produced is limited to 150 m³/d. The injected water rate is equal to 150 m³/d in each injection well. Measurements of pressure and rates of water and oil produced relating to the reference channel were done every 20 days over a simulation period of 1000 days.

<figure><img src="image_6.png" /><figcaption>Figure 4
Example 2: reference channel permeability map.</figcaption></figure>

The absolute permeability is equal to 500 mD along the *Ox* and *Oy* axes inside the channel and 0.0001 mD outside the channel. When the banks of the channel cut a mesh, the permeability in the mesh is calculated using the Cardwell-Parsons homogenization formulas [17]. The initial porosity is uniform and has a value of 0.3 in the reservoir. The general characteristics of the model are summarized in the Table 2.

<table><caption>TABLE 2
Example 2: parameters of the channel model</caption><thead><tr><td>Parameters</td><td>Values</td></tr></thead><tbody><tr><td>Number of horizontal meshes</td><td>25 x 15</td></tr><tr><td>Horizontal dimensions</td><td>2500 m x 1500 m</td></tr><tr><td>Thickness of the reservoir</td><td>10 m</td></tr><tr><td>Depth of the reservoir</td><td>4700 m</td></tr><tr><td>Reference depth</td><td>4700 m</td></tr><tr><td>Rock compressibility</td><td>5·e-4 bar<sup>-1</sup></td></tr><tr><td>Water compressibility</td><td>4·36e-4 bar<sup>-1</sup></td></tr><tr><td>Oil compressibility</td><td>1·e-4 bar<sup>-1</sup></td></tr><tr><td>Density of water</td><td>1 kg/m<sup>3</sup></td></tr><tr><td>Density of oil</td><td>0.625 kg/m<sup>3</sup></td></tr><tr><td>Viscosity of water</td><td>0.31 Pa·s</td></tr><tr><td>Viscosity of oil</td><td>0.34 Pa·s</td></tr><tr><td>Initial pressure</td><td>400 bar</td></tr><tr><td>Initial saturation of water</td><td>0.3</td></tr><tr><td>Number of injection wells</td><td>6</td></tr><tr><td>Injection rate (water)</td><td>300 m³/day</td></tr><tr><td>Number of production wells</td><td>5</td></tr><tr><td>Production rate (water and oil)</td><td>300 m³/day</td></tr></tbody></table>

## Results

Permeability maps corresponding to the initial position of the channel and to the solution obtained after inversion are indicated in Figures 5 and 7. It may be noted that the position obtained is very close to the reference position and enables obtaining a satisfactory match of the measurements generated from the reference over all the production wells. In the aggregate, the volume of the channel is recovered.

<figure><img src="image_7.png" /><figcaption>Figure 5
Example 2: initial guess permeability map.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6
Example 2: water rates match at well Prod. 2.</figcaption></figure>

Figure 8 shows the evolution of the objective function during the 10 iterations of the optimization process. Figures 6 and 9 represent respectively the graphs of water and pressure output at the production well Prod. 2 for the initial position and the setting obtained. Well Prod. 2 is characteristic of the behaviour of production wells for the various positions of the channel. In particular, one can observe that the moments of massive inflow of water to the well are correctly restored.
<figure><img src="image_9.png" /><figcaption>Figure 7

Example 2: inversion permeability map obtained.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 8

Example 2: evolution of the objective function.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 9

Example 2: pressure match at well Prod. 2.</figcaption></figure>

This has a certain importance because these moments enable determining at what moment should be carried out the possible closure of a given well, if its productivity becomes too low.

## 5.4 Example 3: Identification of Several Faults

This third example is based on a more realistic geological model. The geometry and petrophysics of the reservoir are generated with the help of a Heresim 3D, an integrated tool for geo-statistical modelling developed at the Institut français du pétrole (IFP).

### 5.4.1 Construction of the Geological Model

The first stage was to construct the structure of the model. Twenty exploratory wells are used in the model. These wells contain lithofacies information and the reference depth is obtained through a log analysis of the wells. This analysis enabled to define two distinct units in the reservoir. High resolution models of these two units are generated separately with Heresim 3D from the proportion of the graphs representing the vertical distribution of the various lithofacies. The grid used consists of 132 x 76 meshes of 25 m side. The thickness of the fine meshes is 3 m. The average thickness of the highest unit is 40 m and that of the lowest one is 30 m.

Each unit contains three lithofacies with different petrophysical characteristics. The porosity and permeability values in each facies follow respectively gaussian and log normal laws. Figure 10 shows the vertical proportion graphs used to simulate the distributions of the lithofacies in the low and high units of the reservoir. Figure 11 represents a cross-section of the reservoir after combining the two units.

### 5.4.2 Simulation of the Reservoir Model

Only the western half of the geological model was considered to construct a model of reservoir simulation. This part contains sixteen drilling wells. The roof of the reservoir is located at an average depth of 3600 m and its total average thickness is equal to 70 m. The fine meshes are grouped horizontally according to a 4 x 4 size to form a simulation grid of 21 x 19 meshes of side 100 m. Vertically, the geological layers of simulation are 3 in number in the lower part and 2 in the upper part. The permeability values are brought to scale with Heresim 3D by numerical resolution of a steady-state equation in each mesh.

Several large faults have been identified in the model and only two have been taken into consideration for the inversion (Fig. 12). These two faults are supposed to be impervious, vertical and of negligible thickness. From the positions of the two faults proposed after a seismic interpretation and geological studies, the following production diagram was adopted:

- The 10 wells—Prod. 1 to Prod. 10—drilled on the roof of the reservoir are considered as production wells and pass through the reservoir over 2 or 3 layers, from top to bottom.
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>Figure 10

Example 2: the two vertical proportions curves of the model.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 12

Example 3: initial positions of the 2 faults and water saturation map at 3000 days (middle layer).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 11

Example 2: cross section of the reservoir.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 13

Example 3: pressure match at well Prod. 1.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 14

Example 3: oil rates match at production well Prod. 6.</figcaption></figure>
- The wells Inj. 1 to Inj. 6 are considered as injection wells and pass through the reservoir over 3 layers, from bottom to top. The general simulation data and the petrophysical characteristics of the lithofacies are summarized in the Tables 3 and 4.

<table><caption>TABLE 3
Example 3: simulation parameters of the 3D model</caption><thead><tr><th>Parameters</th><th>Values</th></tr></thead><tbody><tr><td>Number of horizontal meshes</td><td>21 x 19 x 5</td></tr><tr><td>Horizontal dimensions</td><td>2100 m x 1900 m</td></tr><tr><td>Thickness of the reservoir</td><td>70 m</td></tr><tr><td>Depth of the reservoir</td><td>3630 m</td></tr><tr><td>Reference depth</td><td>3600 m</td></tr><tr><td>Reference pressure</td><td>300 bar</td></tr><tr><td>Water-oil contact</td><td>3690 m</td></tr><tr><td>Rock compressibility</td><td>5.e-4 bar<sup>-1</sup></td></tr><tr><td>Water compressibility</td><td>4.36e-4 bar<sup>-1</sup></td></tr><tr><td>Oil compressibility</td><td>1.e-4 bar<sup>-1</sup></td></tr><tr><td>Density of water</td><td>1.05 kg/m<sup>3</sup></td></tr><tr><td>Density of oil</td><td>0.625 kg/m<sup>3</sup></td></tr><tr><td>Viscosity of water</td><td>0.381 Pa.s</td></tr><tr><td>Viscosity of oil</td><td>0.715 Pa.s</td></tr><tr><td>Initial pressure</td><td>300 bar</td></tr><tr><td>Irreducible water saturation</td><td>0.15</td></tr><tr><td>Residual oil saturation</td><td>0.20</td></tr><tr><td>k<sub>r</sub> maximum-water phase</td><td>0.65</td></tr><tr><td>k<sub>r</sub> maximum-oil phase</td><td>1</td></tr><tr><td>Number of injection wells</td><td>6</td></tr><tr><td>Injection rate (water)</td><td>300 m<sup>3</sup>/d</td></tr><tr><td>Number of production wells</td><td>10</td></tr><tr><td>Production rate (water and oil)</td><td>250 m<sup>3</sup>/d</td></tr></tbody></table>

### 5.4.3 Inversion Data

Production history was generated with a numerical simulator of two-phase flows for a reference position corresponding to the 2 faults. A random noise was then added on the pressure data as well as on those of oil and water rates on the 10 production wells. Then, the production data of the first

three years were considered as being the measurements. These have been integrated in the objective function.

The calibration parameters are the Cartesian coordinates of the four nodes defining the position and the length of the two faults. These are supposed to be vertical and pass through the entire reservoir from top to bottom.

### 5.4.4 Results

Figure 12 represents the initial positions of the faults and the water saturation maps in the middle layer of the reservoir. Figure 15 shows the fault positions obtained after five iterations of the inversion process. For these positions, the objective function was divided by 100 and a proper match was obtained on the pressure and rates. Figures 13 and 16 correspond to the setting of the pressure data at the most significant wells (Prod. 1 and Prod. 5). Figures 14 and 17 show the match of the oil rate data at the wells Prod. 6 and Prod. 7 on which the impact of the position of the faults was the greatest. It may be noted that these settings were carried out on the first three years of the production history over a total period of a little more than 8 years (3000 days). The structure of the model and the location of the wells are shown on the Figure 18. The results thus obtained through numerical simulation proved to be satisfactory.

## CONCLUSION

The work presented in this article is based on the use and development of several efficient tools adapted to the form inversion. The principal originality of this work lies in the calculation of the gradient. The introduction of an adjoint problem and the use of methods of mathematical analysis enable obtaining the expression of the derivative of the objective function with respect to the geometry on the continuous problem in an analytical manner.

The representation of the geological forms is independent of the Cartesian grid of the reservoir used for the simulation of flows. Which enables to work with structured grids. The calculation of the gradient is afterwards done and effective optimization algorithms can be used.

<table><caption>TABLE 4
Example 3: porosity and permeability values for each lithofacies</caption><thead><tr><th>Units</th><th colspan="3">Lower unit</th><th colspan="3">Upper unit</th></tr><tr><th>Facies</th><th>Facies 1</th><th>Facies 2</th><th>Facies 3</th><th>Facies 4</th><th>Facies 5</th><th>Facies 6</th></tr></thead><tbody><tr><td>Average porosity</td><td>0.25</td><td>0.12</td><td>0.19</td><td>0.2</td><td>0.11</td><td>0.25</td></tr><tr><td>Standard deviation</td><td>0</td><td>3</td><td>3.5</td><td>0</td><td>3</td><td>3.5</td></tr><tr><td>K<sub>h</sub> average (mD)</td><td>100</td><td>9</td><td>350</td><td>100</td><td>7</td><td>235</td></tr><tr><td>Standard deviation</td><td>0</td><td>20</td><td>400</td><td>0</td><td>10</td><td>300</td></tr><tr><td>K<sub>z</sub>/K<sub>h</sub></td><td>1</td><td>0.1</td><td>0.2</td><td>0.2</td><td>0.5</td><td>1</td></tr></tbody></table>
<figure><img src="image_18.png" /><figcaption>Figure 15

Example 3: optimum positions of the 2 faults after 5 iterations.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 18

Example 3: initial water saturation in the 3D simulation model and well positions.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 16

Example 3: pressure match at well Prod. 5.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Figure 19

Example 3: up-scaled $K_s$ permeability in the model - Top layer and two cross sections.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Figure 17

Example 3: oil rates match at production well Prod. 7.</figcaption></figure>

The technique developed was applied in the first instance to mono-phase flows. The encouraging results obtained motivated its generalization to two-phase flows as illustrated in this article. It was coupled with a polyphase flow simulator. The optimization process developed enables to identify 2D and 3D geological forms in heterogeneous media. The geometrical forms obtained are regular.

The academic applications presented concern synthetic cases and a more realistic geological model for history matching of production data. They demonstrate perfectly the efficiency of the method and its validity for several types of geological bodies. Real cases are much more difficult to solve, for example the case of fissured reservoirs. The advantage of the method concerns
particularly the identification of the roof of a reservoir or of
position of large faults which influence the flows. The tools
developed open the way to more and more complex 2D
and 3D applications. Their contribution may be important
in the characterization of reservoirs. They may be
particularly useful to define the location of new drilling
wells and to reduce the uncertainties on the geometry of
the reservoir.

ACKNOWLEDGMENTS

The authors thank the Institut français du pétrole (IFP) and
the Centre européen de recherche et de formation avancée
en calcul scientifique (Cerfacs) which made the publication
of this article possible.

REFERENCES

1 Landa, J.L., Kamal, M.M., Jenkins, C.D. and Horne R. (1996) Reservoir Characterization Constrained to Well Test Data: a Field Example. *SPE 36511*.

2 Anterion, F., Eymard, R. and Karcher, B. (1989) Use of Parameter Gradients for Reservoir History Matching. *SPE 18433*.

3 Bissel, R., Killough, J.E. and Sharma, Y. (1992) Reservoir History Matching Using the Method of Gradients on a Workstation. *SPE 24265*.

4 Rahon, D., Blanc, G. and Guerillot, D. (1996) Gradients Method Constrained by Geological Bodies for History Matching. *SPE 36568*.

5 Roggero, F. and Guerillot, D. (1996) Gradients Method and Bayesian Formalism. Application to Petrophysical Parameter Characterization. *ECMOR V*, Leoben.

6 Chu, L., Reynolds, A.C. and Olover, D.S. (1995) Compu-
tation of Sensitivity Coefficients for Conditioning the
Permeability Field to Well-Test Pressure Data. *In situ*, **19**,
179-223.

7 Chavent, G., Dupuy, M. and Lemonnier, P. (1973) History Matching by Use of Optimal Control Theory. *SPE 4627*, *SPE-AIME*, Las Vegas.

8 Rahon, D., Edoa, P.F. and Masmoudi, M. (1997) Inversion of Geological Shapes in Reservoir Engineering Using Well-Tests and History Matching of Production Data. *SPE 38656*, Texas, 5-8 October 1997.

9 Mohammadi, B. (1996) Optimal Shape Design, Reverse Mode of Automatic Differentiation and Turbulence. *AIAA Paper 97-0099*.

10 Guillaume, P. and Masmoudi, M. (1993) Computation of High Order Derivatives in Optimal Shape Design. *Numerische Mathematik*.

11 Cea, J. (1986) Conception optimale ou identification de formes, calcul rapide de la dérivée directionnelle de la fonction coût, **20**, 3, **MAAN**.

12 Guillaume, P. (1994) Dérivées d'ordre supérieur en conception optimale de forme. *Thesis*, Toulouse.

13 Fletcher, R. and Powell, M.J.D. (1963) A Rapidly Conver-
gent Descent Method for Minimization. *The Computer
Journal*, **6**, 163-168.

14 Bonnans, J.F., Gilbert, J.C., Lemarechal, C. and Sagastizabal, C. (1995) Methodes numeriques d'optimisation. *INRIA, Projet Promath*, Rocquencourt, Sept. 1995.

15 Minoux, M. (1980) *Programmation mathématique : Théorie et algorithmes*, tomes 1 and 2, Dunod.

16 Le Thiez, P., Lemonnier, P. and Le Gallo, Y. (1993) Notice d'utilisation du modèle SARIP<sup>CH</sup>. *IFP Report*, ref. 40927, November 1993.

17 Cardwell, W.T. and Parsons, R.L. (1945) Average Permeability of Heterogeneous Oil Sands. Trans. Am. Inst. Min. Metall. Pet. Eng., **160**, 34-42.

18 Simon, J. (1980) Differentiation with Respect to the Domain, Numerical Functional Analysis and Optimization, **2**, 7-8.

19 Raviart, P.A. and Thomas, J.M. (1989) *Introduction à l'analyse numérique des équations aux dérivées partielles*, Masson.

20 Murat, F. and Simon, J. (1976) Sur le contrôle par un domaine géométrique. *Prépublication du laboratoire d'analyse numérique n° 76015*, Paris 6.

21 Rahon, D., Edoa, P.F. and Masmoudi, M. (1998) Identification of Geological Shapes in Reservoir Engineering by History Matching Production Data. *SPE 48969*, New Orleans, Louisiana, 27-30 September 1998.
# APPENDIX

## Calculation of the Derivative Relative to the Geometry

Here, we show the calculation of the derivative of the objective function $j$ relative to $\Gamma$. We note that actual unknown in this problem of optimization is $\Gamma$. To simplify and to use the formalism introduced by Céa in [11], where one can find the justification for the calculation which will follow, we will actually determine the derivative of the objective function relative to the domain $\Omega$. The techniques used will be the transport method [18] called the parametrization of the domain and the Lagrangian method.

The idea involves replacing the search for the optimal form $\Gamma$ contained in a reference domain $\Omega$, by the search for a transformation $F$ such that $F(\Omega)$, implying $F(\Gamma)$, is optimal.

For this, we consider the internal transformations of the domain $\Omega$ of form $F = I + V$. These transformations belong to an appropriate functional space which we formally denote as $U$ [12]. $V$ denotes a perturbation of the domain $\Omega$. The support of $V$ is defined only in the neighbourhood of $\Gamma$. We denote $\Lambda = F(\Omega) = \Omega + V(\Omega)$ the perturbed domain, $\partial\Lambda$ its border, $(x', y', z') = F(x, y, z)$, $(x, y, z) \in \Omega$, $\Omega_T = \Omega \times (0, T)$ and $\Sigma_T = \Sigma \times (0, T)$.

Let $P_\Lambda$ and $S_\Lambda$ be the solutions of the partial derivatives equations:

$$ \frac{\partial}{\partial t}(\phi(P_\Lambda)S_\Lambda) - \operatorname{div} \left( K \frac{kr}{\mu_w}(S_\Lambda) \operatorname{grad}(P_\Lambda + \rho_w gz') \right) - \frac{q_w}{\rho_w} = 0 \quad (14) \\ (x', y', z') \in \Lambda, t \in [0, T] $$

and:

$$ \frac{\partial \phi(P_\Lambda)}{\partial t} - \operatorname{div} \left( K \frac{kr}{\mu_w}(S_\Lambda) \operatorname{grad}(P_\Lambda + \rho_w gz') \right) - \operatorname{div} \left( K \frac{kr}{\mu_o}(S_\Lambda) \operatorname{grad}(P_\Lambda + \rho_o gz') \right) - \frac{q_w}{\rho_w} - \frac{q_o}{\rho_o} = 0 \quad (15) \\ (x', y', z') \in \Lambda, t \in [0, T] $$

verifying the boundary conditions:

$$ \left( K \frac{kr}{\mu_w} \operatorname{grad}(P_\Lambda + \rho_w gz') + K \frac{kr}{\mu_o} \operatorname{grad}(P_\Lambda + \rho_o gz') \right) \cdot n = 0 \quad (16) \\ (x', y', z') \in \partial\Lambda, t \in [0, T] $$

and the initial conditions:

$$ P_{\Lambda}(x', y', z', 0) = P_0(x', y', z'); (x', y', z') \in \Lambda \quad (17) $$

$$ S_{\Lambda}(x', y', z', 0) = S_0(x', y', z'); (x', y', z') \in \Lambda \quad (18) $$

$P_\Lambda$ and $S_\Lambda$ are the pressure and the water saturation in the “perturbed” reservoir in the presence of a two-phase oil-water flow without any capillary pressure. $P_\Lambda$ and $S_\Lambda$ belong to a functional space which we formally denote as $H(\Lambda)$.

The variational formulation of the problem (14)-(18) is obtained by assuming that $P_\Lambda$ and $S_\Lambda$ are sufficiently regular, next by multiplying the Equation (14) by a test function $u_\Lambda \in H(\Lambda)$ and the Equation (15) by a test function $v_\Lambda \in H(\Lambda)$, and by integrating over $\Lambda \times (0, T)$. Hence, we find that by taking the boundary conditions (16) and using the classic Green formula [19]:

$$ E_1(\Lambda, P_\Lambda, S_\Lambda, u_\Lambda) = 0, \quad \forall u_\Lambda \in H(\Lambda) \quad (19) $$

$$ E_2(\Lambda, P_\Lambda, S_\Lambda, v_\Lambda) = 0, \quad \forall v_\Lambda \in H(\Lambda) \quad (20) $$

with:

$$ E_1(\Lambda, P_\Lambda, S_\Lambda, u_\Lambda) = \int_{\Lambda \times (0, T)} \frac{\partial}{\partial t} (\phi(P_\Lambda)S_\Lambda) u_\Lambda \\ + \int_{\Lambda \times (0, T)} K \frac{kr}{\mu_w} \operatorname{grad}(P_\Lambda + \rho_w gz') \cdot \operatorname{grad} u_\Lambda - \frac{q_w}{\rho_w} u_\Lambda \\ + \int_{\partial\Lambda \times (0, T)} \left( K \frac{kr}{\mu_w} \operatorname{grad}(P_\Lambda + \rho_w gz') \right) \cdot n u_\Lambda \quad (21) $$

and:

$$ E_2(\Lambda, P_\Lambda, S_\Lambda, v_\Lambda) = \int_{\Lambda \times (0, T)} \frac{\partial \phi(P_\Lambda)}{\partial t} v_\Lambda \\ + \int_{\Lambda \times (0, T)} \left( K \frac{kr}{\rho_w} \operatorname{grad}(P_\Lambda + \rho_w gz') \right) \cdot \operatorname{grad} v_\Lambda \\ + \int_{\Lambda \times (0, T)} \left( K \frac{kr}{\mu_o} \operatorname{grad}(P_\Lambda + \rho_o gz') \right) \cdot \operatorname{grad} v_\Lambda \quad (22) $$

Under conditions of sufficient regularity, we assume that $P_F$, $S_F$, $u_F$ and $v_F \in H(\Omega) = \{z \in H(\Lambda), F \in U\}$ are such that:

$$ P_{\Lambda} = P_{F} o F^{-1}, \quad S_{\Lambda} = S_{F} o F^{-1} $$

$$ u_{\Lambda} = u_{F} o F^{-1}, \quad v_{\Lambda} = v_{F} o F^{-1} $$

For this, it is sufficient to show that the application:

$$ f \in H(\Omega) \rightarrow f \circ o F^{-1} \in H(\Lambda), \Lambda = F(\Omega) $$

is an isomorphism as soon as $F$ is in an appropriate space $U$.

The transport method suggests that the transported solutions $P_F = P_\Lambda o F$ and $S_F = S_\Lambda o F$ verify the variational equations after changing of the variables:

$$ \tilde{E}_1(F, P_F, S_F, u) = 0, \quad \forall u \in H(\Omega) \quad (23) $$

$$ \tilde{E}_2(F, P_F, S_F, v) = 0, \quad \forall v \in H(\Omega) \quad (24) $$
with:

$$
\begin{align*}
\tilde{E}_1(F, P_F, S_F, u) &= E_1(\Lambda, P_\Lambda, S_\Lambda, u_\Lambda) \\
&= \int_{\Omega_T} \left( A \frac{\partial S_F}{\partial t} + B S_F \frac{\partial P_F}{\partial t} + B P_F \frac{\partial S_F}{\partial t} - \frac{q_w}{\rho_w} \right) u J F \\
&\quad + \int_{\Omega_T} \left( K \frac{kr_w}{\mu_w} \right) \left( \begin{aligned}[t] &'DF^{-1}\operatorname{grad}(P_F + \rho_w g z) \\ &\cdot DF^{-1}\operatorname{grad} u J F \end{aligned} \right) \\
&\quad + \int_{\partial\Omega_T} \left( K \frac{kr_w}{\mu_w} \right) \left( \begin{aligned}[t] &'DF^{-1}\operatorname{grad}(P_F + \rho_w g z) \\ &\cdot DF^{-1}\operatorname{grad} n J F \end{aligned} \right)
\end{align*}
$$

and:

$$
\begin{align*}
\tilde{E}_2(F, P_F, S_F, v) &= E_2(\Lambda, P_\Lambda, S_\Lambda, v_\Lambda) = \int_{\Omega_T} B \frac{\partial P_F}{\partial t} v J F \\
&\quad + \int_{\Omega_T} \left( K \frac{kr_w}{\mu_w} \right) \left( \begin{aligned}[t] &'DF^{-1}\operatorname{grad}(P_F + \rho_w g z) \\ &\cdot DF^{-1}\operatorname{grad} v J F \end{aligned} \right) \\
&\quad + \int_{\Omega_T} \left( K \frac{kr_o}{\mu_o} \right) \left( \begin{aligned}[t] &'DF^{-1}\operatorname{grad}(P_F + \rho_o g z) \\ &\cdot DF^{-1}\operatorname{grad} v J F \end{aligned} \right)
\end{align*}
$$

where DF is the jacobian of F, DF⁻¹ is the inverse of DF and
JF = |det DF| the determinant of DF; A = φ⁰ (1 - cₛRₚ₀) and
B = φ⁰ cₛR.

We are thus taken into the domain of reference Ω on
which the objective function can be defined on U by:

$$
\begin{align}
\tilde{j}(F) = j(F(\Omega)) &= \frac{1}{2} \int_0^T \beta \left\| P_{F(\Omega)} - P^m \right\|_{L^2(\Omega_m)}^2 \mathrm{d}t \nonumber \\
&\quad + \frac{1}{2} \int_0^T \gamma \left\| q_w^{F(\Omega)} (S_{F(\Omega)}) - q_w^m \right\|_{L^2(\Omega_m)}^2 \mathrm{d}t \tag{25} \\
&\quad + \frac{1}{2} \int_0^T \sigma \left\| q_o^{F(\Omega)} (S_{F(\Omega)}) - q_o^m \right\|_{L^2(\Omega_m)}^2 \mathrm{d}t \nonumber
\end{align}
$$

Here onwards we have to find the transformation F which
minimize the “new definition” j of the objective function j.
For this, one must determine the directional derivative

d<sub>Fj</sub>(F).V of j̃ relative to F in the direction V. To calculate the
first derivative of j̃, we use the Lagrangian method [11].

The Lagrangian method involves defining the following
function on U × H(Ω)⁴:

$$
\begin{equation}
\begin{split}
L(F, P, S, u, v) = \tilde{j}(F) + \tilde{E}_1(F, P, S, u) + \tilde{E}_2(F, P, S, v) \quad &(26) \\
\qquad \forall u \text{ and } v \in H(\Omega)
\end{split}
\end{equation}
$$

in which the direct state equations in the variational form
(23) and (24) appear as constraints.

Since, P<sub>F</sub> and S<sub>F</sub> are the solutions of the Equations (23)
and (24), we precisely have:

$$
L(F, P_F, S_F, u, v) = \tilde{j}(F), \quad \forall u \text{ and } v \in H(\Omega) \quad (27)
$$

It is possible to obtain an expression for the total
derivative of the function j̃ which is nothing but a partial
derivative of the Lagrangian. Thus, having u and v as
constants and taking into consideration (27), we have:

$$
\begin{align*}
\frac{d\tilde{j}}{dF}(F) \cdot V &= \lim_{\epsilon \to 0} \frac{\tilde{j}(F + \epsilon V) - \tilde{j}(F)}{\epsilon} = \frac{dL}{dF}(F, P_F, S_F, u, v) \cdot V \\
&= \frac{\partial L}{\partial F} \cdot V + \frac{\partial L}{\partial P} \left( \frac{\partial P}{\partial F} \cdot V \right) + \frac{\partial L}{\partial S} \left( \frac{\partial S}{\partial F} \cdot V \right)
\end{align*}
$$

The functions u and v can be chosen such that the
derivatives $\frac{\partial L}{\partial P} \left( \frac{\partial P}{\partial F} \cdot V \right)$ and $\frac{\partial L}{\partial S} \left( \frac{\partial S}{\partial F} \cdot V \right)$ can be elimi-
nated. These introduce us to the following equations called
the associated adjoint state equations.

$$
\frac{\partial L}{\partial P} (F, P_F, S_F, u, v) \cdot \psi = 0, \quad \forall \psi \in H(\Omega) \quad (28)
$$

$$
\frac{\partial L}{\partial S} (F, P_F, S_F, u, v) \cdot \varphi = 0, \quad \forall \varphi \in H(\Omega) \quad (29)
$$

The functions u and v are dependent on the parametrization
of the domain Ω. Proceeding from the Equations (27) and
(28) for F = I—these are the ones relevant to us—and by
observing that P = P_I = P_Γ and S = S_I = S_Γ, we get:

$$
\begin{align}
& \int_{\Omega_T} \left( -S\phi'(P) \frac{\partial u}{\partial t} - \phi'(P) \frac{\partial v}{\partial t} - \operatorname{div} \left( K \frac{kr_w(S)}{\mu_w} \operatorname{grad} u \right) \right) \psi + \int_{\Omega_T} \left( -\operatorname{div} \left( \left( K \frac{kr_w(S)}{\mu_w} + K \frac{kr_o(S)}{\mu_o} \right) K \operatorname{grad} v \right) + \beta (P_\Gamma - P^m) \chi_{\Omega_m} \right) \psi \nonumber \\
& \qquad + \int_{\Omega} \left( (S\phi'(P)u\psi)(x, y, z, T) - (S\phi'(P)u\psi)(x, y, z, 0) \right) + \int_{\Omega} \left( (\phi'(P)v\psi)(x, y, z, T) - (\phi'(P)v\psi)(x, y, z, 0) \right) \tag{30} \\
& \qquad + \int_{\Sigma_T} \left( \left( K \frac{kr_w(S)}{\mu_w} \operatorname{grad} u \right) \cdot n \right) \psi + \int_{\Sigma_T} \left( \left( \left( \left( \frac{kr_w(S)}{\mu_w} + \frac{kr_o(S)}{\mu_o} \right) K \operatorname{grad} v \right) \cdot n \right) \psi + \int_{\Sigma_T} \left( K \frac{kr_w(S)}{\mu_w} \operatorname{grad} \psi \right) \cdot nu = 0, \quad \forall \psi \in H(\Omega) \nonumber
\end{align}
$$
and:

$$
\begin{align}
& \int_{\Omega_T} \left( -\phi(P) \frac{\partial u}{\partial t} - \frac{1}{\rho_w} \frac{\partial q_w}{\partial S} u - \left( \frac{1}{\rho_w} \frac{\partial q_w}{\partial S} + \frac{1}{\rho_o} \frac{\partial q_o}{\partial S} \right) v \right) \varphi + \int_{\Omega_T} \left( \frac{K}{\mu_w} \frac{\partial kr_w(S)}{\partial S} \operatorname{grad}(P + \rho_w gz) \cdot \operatorname{grad} u \right) \varphi \nonumber \\
& \qquad + \int_{\Omega_T} \left( \gamma(q_w^\Gamma(S_\Gamma) - q_w^m) \frac{\partial q_w^\Gamma(S_\Gamma)}{\partial S} \right) \chi_{\Omega_m} \varphi + \int_{\Omega_T} \left( \sigma(q_o^\Gamma(S_\Gamma) - q_o^m) \frac{\partial q_o^\Gamma(S_\Gamma)}{\partial S} \right) \chi_{\Omega_m} \varphi + \int_{\Omega_T} \left( \frac{K}{\mu_w} \frac{\partial kr_w(S)}{\partial S} \operatorname{grad}(P + \rho_w gz) \right) \cdot \operatorname{grad} v \varphi \nonumber \\
& \qquad + \int_{\Omega_T} \left( \frac{K}{\mu_o} \frac{\partial kr_o(S)}{\partial S} \operatorname{grad}(P + \rho_o gz) \right) \cdot \operatorname{grad} v \varphi + \int_{\Omega} \left( (\phi(P)u\varphi)(x, y, z, T) - (\phi(P)u\varphi)(x, y, z, 0) \right) \nonumber \\
& \qquad + \int_{\Sigma_T} \left( \frac{K}{\mu_w} \frac{\partial kr_w(S)}{\partial \mu_w} \varphi \operatorname{grad}(P + \rho_w gz) \right) \cdot nu = 0, \quad \forall \varphi \in H(\Omega) \tag{31}
\end{align}
$$

where $\chi_m$ is the characteristic function on $\Omega_m$, the region in the domain where the measurements were carried out, in other words in the production wells.

Next, by making a right choice of functions $\psi$ and $\varphi$, we formally deduce that $u = u_\Gamma$ and $v = v_\Gamma$ verify the problem of evolution defined by the following equations with partial derivatives:

$$
\begin{align}
-S\phi'(P)\frac{\partial u}{\partial t} - \phi'(P)\frac{\partial v}{\partial t} - \operatorname{div}\left(\left(\frac{kr(S)}{\mu_w} + \frac{kr(S)}{\mu_o}\right)K\operatorname{grad}v\right) & \nonumber \\
-\operatorname{div}\left(K\frac{kr(S)}{\mu_w}\operatorname{grad}u\right) + \beta(P_\Gamma - P^m)\chi_{\Omega_m} & = 0 \tag{32} \\
(x, y, z) \in \Omega, \quad t \in [0, T] & \nonumber
\end{align}
$$

$$
\begin{align}
& -\phi(P) \frac{\partial u}{\partial t} - \frac{1}{\rho_w} \frac{\partial q_w}{\partial S} u - \left( \frac{1}{\rho_w} \frac{\partial q_w}{\partial S} + \frac{1}{\rho_o} \frac{\partial q_o}{\partial S} \right) v \nonumber \\
& \quad + \left( \frac{K}{\mu_w} \frac{\partial kr_w(S)}{\partial S} \operatorname{grad}(P + \rho_w gz) \right) \cdot \operatorname{grad} u \nonumber \\
& \quad + \left( \frac{K}{\mu_w} \frac{\partial kr_w(S)}{\partial S} \operatorname{grad}(P + \rho_w gz) \right) \cdot \operatorname{grad} v \nonumber \\
& \quad + \left( \frac{K}{\mu_o} \frac{\partial kr_o(S)}{\partial S} \operatorname{grad}(P + \rho_o gz) \right) \cdot \operatorname{grad} v \tag{33} \\
& \quad + \left( \gamma(q_w^\Gamma(S_\Gamma) - q_w^m) \frac{\partial q_w^\Gamma(S_\Gamma)}{\partial S} \right) \chi_{\Omega_m} \nonumber \\
& \quad + \sigma \left( q_o^\Gamma(S_\Gamma) - q_o^m \right) \frac{\partial q_o^\Gamma(S_\Gamma)}{\partial S} \right) \chi_{\Omega_m} = 0 \nonumber \\
& \quad (x, y, z) \in \Omega, \quad t \in [0, T] \nonumber
\end{align}
$$

the boundary conditions:

$$ u = 0; \ (x, y, z) \in \Sigma, \ t \in [0, T] \tag{34} $$

and:

$$
\begin{align}
& \left( K \frac{kr(S)}{\mu_w} \operatorname{grad} u \right) \cdot n + \left( \left( \frac{kr(S)}{\mu_w} + \frac{kr(S)}{\mu_o} \right) K \operatorname{grad} v \right) \cdot n = 0; \tag{35} \\
& \qquad (x, y, z) \in \Sigma, \ t \in [0, T] \nonumber
\end{align}
$$

and at the instant $T$ with the conditions as follows:

$$ u(x, y, z, T) = 0 \quad \text{and} \quad v(x, y, z, T) = 0; \ (x, y, z) \in \Omega \tag{36} $$

We thus have the result which follows: for $P = P_F$ and $S = S_F$ solutions of the variational direct state Equations (23) and (24) and for $u = u_F$ and $v = v_F$ solutions of the adjoint state Equations (28) and (29), the derivative of $\tilde{j}$ relative to $F$ in the direction $V$ is given by:

$$ \frac{d\tilde{j}}{dF}(F) \cdot V = \frac{\partial L}{\partial F}(F, P_F, S_F, u_F, v_F) \cdot V \tag{37} $$

If we consider the applications $\delta$: $F \in U \to JF$ and $\theta$: $F \in U \to ^tDF^{-1}$, we can show that they are of the class $C^\infty$ for $F$ in an adequate space $U$ and that their derivatives are:

$$ \frac{\partial \theta}{\partial F}(F) \cdot V = {}^tDF^{-1} {}^tDV {}^tDF^{-1} \tag{38} $$

and:

$$ \frac{\partial \delta}{\partial F}(F) \cdot V = DF_c : DV = \operatorname{tr} [ {}^tDF_c DV ] \tag{39} $$

${}^tDF_c$ is the co-matrix of $DF$ and is defined by the relation $DF {}^tDF_c = (\det DF)I$; tr is the trace operator.

Taking into consideration the formulas (38) and (39), and by making the hypothesis that the perturbation $V$ is null on $\Omega_m$, the region where the source terms relative to the wells are situated, we have:
$$
\begin{aligned}
\frac{d\tilde{j}}{dF}(F) \cdot V &= \frac{\partial L}{\partial F}(F, P_F, S_F, u_F, v_F) \cdot V = \frac{d\tilde{j}}{dF}(F) \cdot V = \frac{\partial \tilde{E}}{\partial F}(F, P_F, S_F, u_F) \cdot V + \frac{\partial \tilde{E}_2}{\partial F}(F, P_F, S_F, v_F) \cdot V \\
&= \int_{\Omega_T} \left( \left( A \frac{\partial S_F}{\partial t} + B S_F \frac{\partial P_F}{\partial t} + B P_F \frac{\partial S_F}{\partial t} \right) u_F \right) DF_c : DV + \int_{\Omega_T} \left( B \frac{\partial P_F}{\partial t} v_F \right) DF_c : DV \\
&\quad - \int_{\Omega_T} K \frac{kr}{\mu_w} {}^tDF^{-1} {}^tDV {}^tDF^{-1} \text{grad}(P_F + \rho_w gz) {}^tDF^{-1} \text{grad} u_F JF - \int_{\Omega_T} K \frac{kr}{\mu_w} {}^tDF^{-1} \text{grad}(P_F + \rho_w gz) {}^tDF^{-1} {}^tDV {}^tDF^{-1} \text{grad} u_F JF \\
&\quad + \int_{\Omega_T} K \frac{kr}{\mu_w} {}^tDF^{-1} \text{grad}(P_F + \rho_w gz) {}^tDF^{-1} \text{grad} u_F DF_c : DV - \int_{\Omega_T} K \frac{kr}{\mu_w} {}^tDF^{-1} {}^tDV {}^tDF^{-1} \text{grad}(P_F + \rho_w gz) {}^tDF^{-1} \text{grad} v_F JF \\
&\quad - \int_{\Omega_T} K \frac{kr}{\mu_o} {}^tDF^{-1} {}^tDV {}^tDF^{-1} \text{grad}(P_F + \rho_o gz) {}^tDF^{-1} \text{grad} v_F JF + \int_{\Omega_T} K \frac{kr}{\mu_w} {}^tDF^{-1} \text{grad}(P_F + \rho_w gz) {}^tDF^{-1} \text{grad} v_F DF_c : DV \\
&\quad + \int_{\Omega_T} K \frac{kr}{\mu_o} {}^tDF^{-1} \text{grad}(P_F + \rho_o gz) {}^tDF^{-1} \text{grad} v_F DF_c : DV - \int_{\Omega_T} K \frac{kr}{\mu_w} {}^tDF^{-1} \text{grad}(P_F + \rho_w gz) {}^tDF^{-1} {}^tDV {}^tDF^{-1} \text{grad} v_F JF \\
&\quad - \int_{\Omega_T} K \frac{kr}{\mu_o} {}^tDF^{-1} \text{grad}(P_F + \rho_o gz) {}^tDF^{-1} {}^tDV {}^tDF^{-1} \text{grad} v_F JF
\end{aligned}
$$

with: $A = \phi^0(1 - c_R P_0)$ and $B = \phi^0 c_R$.

We note the following fundamental result [11]: the partial derivative of j relative to F at the point F = I coincides with

the total derivative of j relative to Ω (or Γ) for every perturbation V of the domain. Finally, given that $P_Γ = P_p$ $S_Γ = S_p$, $u_Γ = u_p$, and $v_Γ = v_p$, we get:

$$
\begin{aligned}
d_{\Gamma} j(\Gamma) \cdot V = d_{\Omega} j(\Omega) \cdot V = d_{\tilde{F}} \tilde{j}(I) \cdot V = & \int_{\Omega_{\Gamma}} \left( A \frac{\partial S_{\Gamma}}{\partial t} + B S_{\Gamma} \frac{\partial P_{\Gamma}}{\partial t} + B P_{\Gamma} \frac{\partial S_{\Gamma}}{\partial t} \right) u_{\Gamma} \mathrm{div} V + \int_{\Omega_{\Gamma}} \left( B \frac{\partial P_{\Gamma}}{\partial t} v_{\Gamma} \right) \mathrm{div} V \\
& + \int_{\Omega_{\Gamma}} \frac{kr}{\mu_w} K \mathrm{grad}(P_{\Gamma} + \rho_w gz) \cdot \mathrm{grad} u_{\Gamma} \mathrm{div} V - \int_{\Omega_{\Gamma}} \frac{kr}{\mu_w} (K {}^t DV + DV K) \mathrm{grad}(P_{\Gamma} + \rho_w gz) \cdot \mathrm{grad} u_{\Gamma} \\
& + \int_{\Omega_{\Gamma}} \left( \frac{kr}{\mu_w} K \mathrm{grad}(P_{\Gamma} + \rho_w gz) \right) \cdot \mathrm{grad} v_{\Gamma} \mathrm{div} V + \int_{\Omega_{\Gamma}} \left( \frac{kr}{\mu_o} K \mathrm{grad}(P_{\Gamma} + \rho_o gz) \right) \cdot \mathrm{grad} v_{\Gamma} \mathrm{div} V \\
& - \int_{\Omega_{\Gamma}} \frac{kr}{\mu_w} (K {}^t DV + DV K) \mathrm{grad}(P_{\Gamma} + \rho_w gz) \cdot \mathrm{grad} v_{\Gamma} - \int_{\Omega_{\Gamma}} \frac{kr}{\mu_o} (K {}^t DV + DV K) \mathrm{grad}(P_{\Gamma} + \rho_o gz) \cdot \mathrm{grad} v_{\Gamma}
\end{aligned}
\qquad (40)
$$