A method for building a depositional space corresponding to a geological domain includes the steps of (i) partitioning the present day geological domain with at least one conformal mesh sensibly matching the boundaries of said geological domain, (ii) calculating depositional coordinates defining a depositional space, wherein the depositional coordinates calculations includes calculations of fields of displacement including simulations of mechanical deformations in the geological domain, using a solid material deformation model. A computer program implementing the method is also described.
|
1. A method for building a depositional space corresponding to a geological domain of a stratified terrain that comprises horizons, the method comprising:
partitioning said geological domain with at least one conformal mesh comprising element coordinates, said at least one conformal mesh matching boundaries of said geological domain, wherein the boundaries comprise the horizons,
calculating depositional coordinates defining a depositional space to represent the horizons as substantially planar and parallel,
wherein said depositional coordinates calculation comprises computer calculations of fields of displacement along the element coordinates with specified boundary conditions for the horizons using a solid material deformation model for simulating mechanical deformations in the geological domain,
wherein said depositional coordinates calculation provides for attribution of point coordinates in the geological domain to point coordinates in the depositional space, allowing points to be mapped therebetween, and
wherein said depositional space provides for computing physical properties of the stratified terrain based on measured data for the stratified terrain.
25. A computer-readable non-transitory storage medium comprising computer-executable instructions to instruct a computer to:
build a depositional space corresponding to a geological domain of a stratified terrain that comprises horizons by
partitioning said geological domain with at least one conformal mesh comprising element coordinates, said at least one conformal mesh matching boundaries of said geological domain, wherein the boundaries comprise the horizons
calculating depositional coordinates defining a depositional space to represent the horizons as substantially planar and parallel,
wherein said depositional coordinates calculation comprises calculations of fields of displacement along the element coordinates with specified boundary conditions for the horizons using a solid material deformation model for simulating mechanical deformations in the geological domain,
wherein said depositional coordinates calculation provides for attribution of point coordinates in the geological domain to point coordinates in the depositional space, allowing points to be mapped therebetween, and
wherein said depositional space provides for computing physical properties of the stratified terrain based on measured data for the stratified terrain.
2. A method according to
u(x,y,z)=x+dx(x,y,z), v(x,y,z)=y+dy(x,y,z), w(x,y,z)=z+dz(x,y,z), wherein the calculation of the depositional coordinates (u, v, w) comprises:
boundary matching, comprising an operation of associating the boundaries of the geological domain each with a matching set of mesh facets or elements, said set including at least one facet of a 3D conformal mesh or element of a 2D conformal mesh,
mechanical simulation setting up, comprising a definition of the specified boundary conditions aiming at constraining the calculated geometry of the geological domain in the depositional space,
mechanical simulation, comprising a calculation in the geological domain (x, y, z) of a field of displacements (dx, dy, dz) matching said specified boundary conditions.
3. A method according to
4. A method according to
an output of the mechanical simulation step, and
a value which is explicitly specified.
5. A method according to
6. A method according to
7. A method according to
8. A method according to
9. A method according to
10. A method according to
11. A method according to
12. A method according to
13. A method according to
14. A method according to
15. A FEM method according to
16. A FEM method according to
17. A FEM method according to
18. A FEM method according to
building a global stiffness matrix, and
iteratively building and updating local stiffness matrices defined at the nodes of the 3D mesh representing the geometry of the geological domain.
19. A FEM method according to
retrieving the mesh element containing the location (x, y, z),
computing the depositional coordinates (u, v, w) of said location (x, y, z) using an analytical interpolation scheme based on the values of the depositional coordinates (u, v, w) corresponding to the nodes coordinates (x, y, z) of said container mesh element,
retrieving the value of said terrain property corresponding to coordinates (u, v, w) in the depositional space.
20. A method according to
21. A BEM method according to
22. A BEM method according to
23. A BEM method according to
24. A BEM method according to
computing the depositional coordinates (u, v, w) of said location (x, y, z) used as an observation point, and
retrieving the value of said terrain property corresponding to coordinates (u, v, w) in the depositional space.
|
The general field of the invention is, not limited to, the modeling of stratified terrains in the subsurface, in oil reservoir geosciences notably. The invention pertains to a method for building a virtual and physically-reliable depositional space matching the chronostratigraphic environment at the time of deposition of geological terrains, said depositional space being then used as a computational space where many applications can be run advantageously.
Such applications include accurate modeling and interpolation of geological terrains properties such as porosity or permeability, controlling the quality of seismic reflection data interpretation, or building grids fitting the complex geometry of geological structures.
There are potentially many applications that benefit from the use of a computational space that models the chronostratigraphic environment at the time of deposition of the terrains (which is referred to as the “depositional space”). In the following, as a review of the prior art, the specific example of the modeling of the physical properties of geological terrains is put forward in order to provide a clear understanding of the background of the invention.
The physical properties of a 3D geological domain such as porosity or permeability are usually modeled using geostatistical methods. These methods interpolate the ground properties inside a 3D high-resolution mesh on the basis of a usually sparse data, such as measurements obtained from wells. This interpolation process makes an intensive use of distances, such as Euclidean distances between the centers of mesh elements, or of correlation distances given by variograms.
The computed distributions of physical properties must reflect the paleo-environment at the time of deposition of the terrains, so the interpolation between the data obtained from the wells is accurate only if the computed Euclidean distances are close to the equivalent distances at the time of deposition (which are referred to as the “geodesic distances”).
However, rocks have usually been altered since the time of deposition by erosion or faulting and folding caused by applied tectonic stresses. As a consequence, the present-day geometry of a geological domain is usually significantly different from its geometry at the time of deposition.
Thus, geostatistical methods produce inaccurate results if applied directly to the present-day space, described by the Cartesian coordinate system (x, y, z). This problem can be overcome by rather applying geostatistical methods in a “computational space” or “depositional space”, which aims at matching the chronostratigraphic environment at the time of deposition. This computational or depositional space is usually fitted with a curvilinear coordinate system (u, v, w) also called “computational coordinate system” or “depositional coordinate system”.
In all the following, an “isochron” refers to a surface joining points of the present-day space where sediments deposited at the same time. Seismic horizons or the top of geological layers that do not represent a gap in the geologic record are typically isochrons. Boundaries of sequences corresponding to stratigraphic unconformities are not isochrons. “Gaps” refer to situations where points of the depositional space have no corresponding location in the present-day space, and “overlaps” refer to situations where points of the depositional space have more than one corresponding location in the present-day space.
Depositional spaces, as defined in prior art documents, are expected to honor the two following properties:
So far, two main types of solutions can be found in the prior art for defining such depositional space.
A first type of solutions can be referred to as “(i, j, k) indexation” solutions, as they consist in building a 3D structured conformal mesh in the present-day space with a Cartesian coordinate system (x, y, z), which represents the present-day geometry of the geological domain. This mesh is structured, so an (i, j, k) index can be assigned to every node and element of the mesh. This indexation is such that the neighbors of a node or an element can be found by simple transformations of its index, these transformations being the same for all the nodes and elements of the mesh. Examples of such transformations are (i−1, j, k), (i, j−1, k), (i, j, k−1), (i+1, j, k), (i, j+1, k) and (i, j, k+1).
Once the (i, j, k) indexation is completed, the depositional space is very simply defined by a curvilinear coordinate system (u, v, w) such as u(x, y, z)=i(x, y, z), v(x, y, z)=j(x, y, z) and w(x, y, z)=k(x, y, z). As the mesh is conformal, there exists a set of facets in the mesh corresponding to every horizon. As a consequence, if all the nodes of these facets have the same k index, the property “P1” is honored. In other words, the mesh must be what is often called a “stratigraphic grid”.
This solution suffers from several drawbacks and limitations:
A second type of solutions can be referred to as “parametric” solutions, as they consist in computing depositional coordinates with an interpolation of some quantities following a set of geometric rules.
An example of such a geometric solution is given in documents WO 03/050766 by Deny et al. and WO 2005/119304 by Dulac et al., both incorporated herein by reference. Three transfer functions u(x, y, z), v(x, y, z) and w(x, y, z) are interpolated at the nodes of a 3D conformal mesh representing the present-day geometry of the geological domain. The transfer function w(x, y, z) is computed first. In order to meet the requirements described above, the interpolation of the w(x, y, z) transfer function is done so that specific iso-values (iso-surfaces) of the function approximate the geometry of the horizons, which guarantees the property “P1” to be honored. Then, the two other transfer functions u(x, y, z) and v(x, y, z) are computed. The interpolation of these transfer functions is done so that their gradients are orthogonal both to each other and to the gradient of the w(x, y, z) transfer function previously computed, and so that the length of their gradients are equal. Moreover, as an additional boundary condition, the values of these two transfer functions are computed on a reference horizon, specifying that their gradients are orthogonal to each other and that the lengths of their gradients are equal. This reference horizon is chosen so that it intersects a maximum number of faults. WO 2005/119304 discloses running the interpolation on a tetrahedral mesh (a mesh whose elements are all tetrahedra).
This solution suffers from several drawbacks and limitations:
The purpose of the invention is to propose a method for physically computing a curvilinear coordinate system (u, v, w) defining a depositional space.
In order to reach that objective, the invention proposes a method for building a depositional space corresponding to a geological domain, comprising the steps of:
wherein said depositional coordinates calculation comprises calculations of fields of displacement including simulations of mechanical deformations in the geological domain, using a solid material deformation model.
The depositional space as computed by the method according to the invention may advantageously be defined (without any intention to be limitative) so that it honors the two following properties (referred to as the “P1” and “P2” properties):
Advantageously,
The solid material mechanical behavior can be any behavior, such as an elastic behavior, a plastic behavior, a behavior involving the viscosity for instance, or any combination of behaviors.
The calculation of the depositional coordinates (u, v, w) can include:
So the invention makes possible the computation of a depositional coordinate system (u, v, w), advantageously honoring the “P1” and “P2” properties, where the geometry of the 3D geological domain is physically reliable, thereby allowing the definition of a low-distorted computational space where measured Euclidean distances are close to geodesic distances. When applied in that depositional space, geostatistical interpolations or simulation algorithms can advantageously model the physical properties of the 3D geological domain with an improved accuracy. The generation of meshes honoring the complex geometries of geological structures can also be made much easier when performed in that depositional space. As another example, the quality of a 3D seismic reflection interpretation can be more precisely checked if run in that depositional space.
In all the following, the word “horizon” is used in a very general sense, referring to a specific isochron identified within a sequence of the geological domain, without representing a gap in the geologic record. In the depositional coordinate system (u, v, w) computed by the method according to the invention, horizons may be planar and parallel (“P1” property). As a consequence, layers may not be faulted in the depositional space.
Compared to the one produced by a parametric method from prior art such as “Space-time mathematical framework for sedimentary geology” by Mallet, a fundamental difference is that in a depositional space computed with the present invention:
With the present invention, faults (or any other boundaries) are free to slip and deform, provided that mechanical contact is honored (“P2” property). This means that equilibrium of forces has been reached, without overlapping or penetration across the faults. The corresponding slip directions and magnitudes are thus an output of the computation instead of an input. The produced slip vectors result from a simulation that honors the fundamental principles of both fracture mechanics and physics that govern rock deformation. As a consequence, they provide a reliable physically-based geometry of the geological domain in the depositional space. The same remark applies to horizons: with a parametric method such as the one disclosed in document WO 2005/119304, it is not possible to model the bedding slip/detachment (or any other kind of sliding along a horizon) that might have affected some layers of the geological domain, whereas it is very easy to do so with the method according to the invention, just by setting a boundary condition that allows mechanical contact with shearing at layers interfaces.
Moreover, the present invention does not simplify the processed geological model. All the faults are taken into account, whatever the connections between them. Faults can have free borders and do not need to be artificially extended so that to be entirely connected with other boundaries of the model. The present invention is also not limited to a particular tectonic context. For example, extensional structures or compressive structures are handled exactly the same way. So the present invention can advantageously handle the full complexity of geological models.
Advantageously, in order to honor the “P1” property, the mechanical simulation setup step can further comprise the definition of a boundary condition prescribing that the nodes matching a given horizon of the geological domain have substantially the same final w(x, y, z) value,
In order to honor the “P2” property, the mechanical simulation setup step can further comprise the definition of:
Advantageously,
In a first embodiment, referred to as the “FEM” embodiment, the invention further provides a method wherein the field of displacements (dx, dy, dz) is calculated using a Finite Element Method (FEM), the (present-day) geometry of the geological domain being represented by at least one 3D conformal mesh, the elements of said 3D conformal mesh being polyhedra.
Specific aspects of that embodiment are as follows:
The method according to the “FEM” embodiment can further comprise the computation at location (x, y, z) in the geological domain of the value of a terrain property calculated (that is for instance simulated by a geostatistical method performed in the depositional space) at location (u, v, w) in the depositional space, said terrain property computation comprising the following steps:
According to some specific aspects:
In a second embodiment, referred to as the “BEM” embodiment, the invention provides a method wherein the field of displacements (dx, dy, dz) is calculated using a Boundary Element Method (BEM), the (present-day) geometry of the geological domain being represented by a set of 2D conformal meshes embedded in 3D, each of these meshes representing a boundary of said geological domain, the elements of said 2D meshes being polygons.
Specific aspects of that embodiment are as follows:
The method according to the “BEM” embodiment can further comprise the computation of the value at location (x, y, z) of a terrain property calculated (that is for instance simulated by a geostatistical method performed in the depositional space) at location (u, v, w) in the depositional space, said terrain property computation comprising the following steps:
It is a further object of the invention to provide a computer program implementing the method.
The methods according to embodiments of the present invention may be better understood with reference to the drawings, being understood that these drawings are given for illustrative purposes only and are not meant to be limiting. Other aspects, goals and advantages of the invention shall be apparent from the descriptions given hereunder.
With reference to
In the following description, several embodiments of the invention are described and specific examples are put forward in order to provide a thorough understanding of the invention. However, it will also be apparent to someone skilled in the art that other embodiments are possible and the invention is not limited to the examples described herein. Furthermore, features usually well known by someone skilled in the art may be omitted or simplified in order to avoid obscuring the description of the embodiments of the invention.
With reference to
The present time geological domain 1 of
With reference to
There are potentially many possible implementation methods for the boundary matching step, depending on the algorithms and the data structures it uses. Usually, the more complex a data structure, the simpler the algorithm will be, although the simplicity of the algorithm is altered by the data structure update management. On the other hand, a simple data structure is efficient to update but forces the algorithm to perform explicitly a set of operations to recreate some data needed each time. As a consequence, the final choice is often a compromise.
The geological IDs associated with the boundaries of the geological domain 1, B={B1, B2, B3, . . . , Bn}, where B is the restricted set of boundaries 4 involved in the further steps of the method, can be represented by a series of n arbitrary integer numbers {b1, b2, b3, . . . , bn}. Although not mandatory and for the sake of clarity, B may be represented by a linear list of n records, that is a set of n records whose structural properties essentially involve only the relative position of the records. Record positions can be put into a one-to-one correspondence with the set of the first n integers, so that the geological ID bi, associated with the ith record of the list B, may be equal to i.
With reference to
The method assumes that the tetrahedral mesh is conformal, so there is at least one facet that matches each of the boundaries 4 contained in the set B. The values of aj may be set during an initialization step taking place at the beginning of the boundary matching step. Then, retrieving all the facets of the tetrahedral mesh that correspond to a given boundary Bi of the geological domain may involve the following steps:
The list of facets Fi may be further split into two lists Fi+ and Fi−. The list Fi+ contains all the facets of Fi having a normal vector that locally matches (in the most general sense) the normal vector of the boundary Bi, the normal vector of a facet (T, I) being for example defined as the normal vector to the plane defined by the three nodes of the facet pointing towards the exterior of T. The list Fi− contains all the other facets of Fi. These two lists may be used to retrieve all the facets that correspond to a given side (“+” or “−”) of a given boundary of the geological domain.
In that first embodiment or as soon as the present-day geometry of the geological domain 1 is represented by a 3D conformal mesh 9 whose elements 11 are general polyhedra, the unknown displacements (dx, dy, dz) can be advantageously determined using a Finite Element Method (FEM). The FEM method is a classical numerical method, well known from the one skilled in the art, so it will be only partially described here in order not to obscure embodiments of the invention. A detailed description can be for instance found in the book by Zienkiewicz and Taylor, “The Finite Element Method—Volume 1: The Basis”, Butterworth-Heinemann editions (2000).
In that embodiment, referred to as the FEM method, every element 11 of the 3D conformal mesh 9 is associated with at least two parameters representing the mechanical properties of the material contained in the portion of the geological domain corresponding to the inside of that element. The number of parameters depends on the solid material deformation model implemented. An elastic deformation model requires two parameters, while plastic or more advanced deformation models require more parameters. The values of the parameters can differ from element 11 to element 11 through the mesh 9 and also change with time during the mechanical simulation step.
If the material is assumed to deform elastically during the mechanical simulation step, these parameters are usually the Poisson's ratio and the Young's modulus (or the corresponding Lame coefficients) relative to the material contained in the mesh element.
In order to honor the “P1” property, a boundary condition is set to prescribe that all the nodes 10 of the facets 12 corresponding to a given horizon 5 of the geological domain 1 have the same final w(x, y, z) value, so as to appear as a plane in the depositional space as shown in
In order to honor the “P2” property, another boundary condition is set to prescribe that all the nodes 10 of the facets 12 corresponding to a given side of a given fault 7 or horizon 5 of the geological domain 1 are free to move provided that they stay in contact with the facets corresponding to the other side of the fault 7 or horizon 5. Optionally, another boundary condition may be set to prescribe that all the nodes 10 located at the intersection between a given horizon 5 of the geological domain 1 and a given side of a given fault 7 of the geological domain 1 are free to move provided that they stay in contact with the edges corresponding to the intersection between the horizon 5 and the other side of the fault 7.
The computation of the unknown displacements (dx, dy, dz) at the location of each node 10 of the 3D mesh 9 representing the present-day geometry of the geological domain 1 is done with a Finite Element Method based on the boundary conditions defined in the previous mechanical simulation setup step.
Advantageously, each sequence 8 identified in the geological domain 1 can be processed separately, leading to a representation of its own depositional space 2, as shown on
The FEM computation can be done by building a global stiffness matrix as done in any classical implicit Finite Element Method known in the art.
Alternatively, The FEM computation can be done iteratively by building and updating local stiffness matrices defined at the nodes 10 of the 3D mesh 9 representing the present-day geometry of the geological domain 1, allowing forces to be transmitted from node to node (or from a set of nodes to other sets of nodes) until equilibrium is obtained as done in any “dynamic relaxation”, “multi-agent”, Block Gauss-Seidel or Block Jacobi method known in the art. Such methods are described for instance in the book by Saad, “Iterative Methods for Sparse Linear Systems”, Siam editions (2003).
The FEM computation is not limited to the small displacement hypothesis and can handle large deformations, avoiding in that case the undesirable side effects of rotations.
In order to constrain the nodes 10 to stay mechanically in contact with sets of facets 12 or edges, as it might be prescribed by some boundary conditions defined in the mechanical simulation setup step, a “master/slave” formulation is used. That formulation can be for instance based on a classical penalty method or on the use of Lagrange multipliers.
In the following, an example of application of the Finite Element Method to the computation of mechanical deformation of solid materials is given, based on an elastic deformation model. Only a brief overview is given, as the FEM method itself is common knowledge.
Under the small displacement hypothesis, the Green-Lagrange deformation tensor ε is defined by:
where dI,J denotes the Jth partial derivative of the Ith component of the unknown displacement d=(dx, dy, dz). For isotropic elastic solids, the stress tensor σ is given by the Hooke's law:
where the Lamé coefficients λ and μ are defined from the Young's modulus E and Poisson's ratio v:
The momentum conservation equation in the quasi-static case gives:
div σ(d)=0
with usual displacement and boundary conditions. In particular, when nodes 10 are required to stay mechanically in contact with sets of facets 12 or edges, a non-penetration constraint is set, prescribing that the distance of the projection of such nodes on such sets of facets or edges to the local normal plane to the nodes is null (“master/slave” formulation).
According to the Finite Element Method, the unknown displacement field d is searched in the kinematic admissible displacement field D. We also define the space V of “test” admissible virtual displacement field v. The virtual work principle gives the momentum conservation law as an integral form:
where Ωs is the considered solid defined by its spatial coordinates. The latter equation can be considered as a minimization problem of a functional. Then, in order to take into account the non-penetration constraints, the minimization problem is rewritten, replacing the functional by a Lagrangian function L(d, λ), where λ is an unknown Lagrange multiplier. Let h be the mesh used for the discretization of the considered solid. The Finite Element Method approximates the unknowns dh and λh on finite dimensional spaces, and a non-linear vectorial system R(dh, λh) may be obtained by enforcing:
where {g1, . . . , gn} and {Λ1, . . . , Λn} are the basis of given shape functions approximating displacements and Lagrange multipliers respectively. The resulting non-linear system can be solved iteratively in order to nullify the first order residual at step (n+1) according that R(dh, λh) is not null at step n. Defining:
it follows that:
is the tangent matrix of the residual R. This matrix system can be directly inversed in two dimensions with a LU factorization, or solved iteratively in three dimensions with an Uzawa or Minimal Residual Fashion algorithm (see for example Saad, “Iterative Methods for Sparse Linear Systems”, Siam editions (2003)).
Finally, the solved nodal displacements (dx, dy, dz) are used to define the depositional coordinate system (u, v, w) at the location of each node 10 used in the 3D conformal mesh 9. Referring to
The computation of the value at any present-day location (x, y, z) of a property 30 simulated by any geostatistical method performed in the computed depositional coordinate system (u, v, w) can be then done with the following operations:
In a second embodiment, the geological domain 1 can be represented by a set M of 2D independent meshes (for example made of triangles 12), each of these meshes representing a boundary 4 of the set B. The generation of the 2D meshes is not within the scope of the invention and can be done following any classical method known in the art. M can be represented for instance by a linear list M={M1, M2, M3, . . . , Mn}. The data structure representing such a triangular mesh Mj is encoded so that memory slots are reserved to store and retrieve an integer value A(Mj). The value of A(Mj) is set so that if the mesh Mj represents the boundary Bi, then A(Mj) is equal to the corresponding geological ID bi (equal to i for example). The value of the A(Mj) may be set during an initialization step taking place at the beginning of the boundary matching step. Then, retrieving all the triangles that correspond to a given boundary Bi of the geological domain may involve the following steps:
In that second embodiment, the unknown displacements (dx, dy, dz) are advantageously determined using a Boundary Element Method (BEM). The BEM method is a classical numerical method, well known from the one skilled in the art, so it will be only partially described here in order not to obscure embodiments of the invention. A detailed description can be for instance found in the book by Aliabadi, “The Boundary Element Method: Applications in Solids and Structures—Vol. 2”, Wiley editions (2002).
In the context of geomechanics, a Boundary Element representation applies when the partial differential equations (PDEs) governing the deformation of rocks are piecewise homogeneous, and for which a Green's function can be found or approximated. In such cases, the unknown displacement of the system can be reduced to the self-contained displacement of the interfaces separating layers 3 with homogeneous mechanical properties. The Boundary Element Method uses the given boundary conditions to fit boundary values into the integral equation.
The BEM method assumes that the present-day geometry of the geological domain is represented by a set of 2D meshes, each of these meshes representing a boundary 4 of said geological domain (see
Each individual polygon 12 forming the 2D meshes can be associated with mixed traction-displacement boundary conditions, and when such boundary conditions are specified, the method solves the corresponding unknown displacements. The polygons can advantageously be triangles, but common Boundary Element Methods naturally handle other types of polygons.
In that embodiment, referred to as the BEM embodiment, every layer of the geological domain is associated with at least two parameters representing the mechanical properties of the material contained in the portion of the geological domain corresponding to the inside of that layer, thus defining a homogeneous material for that layer. The number of parameters depends on the solid material deformation model implemented. An elastic deformation model requires two parameters, while plastic or more advanced deformation models require more parameters. The values of the parameters can change with time during the mechanical simulation step.
If the material is assumed to deform elastically during the mechanical simulation step, these parameters are usually the Poisson's ratio and the Young's modulus (or the corresponding Lame coefficients) relative to the material contained in the layer.
In order to honor the “P1” property, a boundary condition is set to prescribe that all the nodes 10 of the elements 12 corresponding to a given horizon 5 of the geological domain 1 have the same final w(x, y, z) value, so as to appear as a plane in the depositional space as shown in
In order to honor the “P2” property, another boundary condition is set to prescribe that all the nodes 10 of the elements 12 corresponding to a given side of a given fault 7 or horizon 5 of the geological domain 1 are free to move provided that they stay in contact with the elements corresponding to the other side of that fault 7 or horizon 5. Optionally, another boundary condition can be set to prescribe that all the nodes 10 located at the intersection between a given horizon 5 of the geological domain 1 and a given side of a given fault 7 of the geological domain 1 are free to move provided that they stay in contact with the edges corresponding to the intersection between said horizon 5 and the other side of said fault 7.
The computation of the unknown displacements (dx, dy, dz) at the location of each node 10 of the 2D meshes representing the present-day geometry of the geological domain 1 is done with a Boundary Element Method based on the boundary conditions defined in the previous mechanical simulation setup step. Each sequence 8 identified in the geological domain 1 can be processed separately. The BEM computation is not limited to the small displacement hypothesis and can handle large deformations, avoiding in that case the undesirable side effects of rotations.
Optionally, other unknown displacements (dx, dy, dz) can be computed at arbitrary present-day locations (x, y, z) referred to as “observation points”, using the computed displacements at the location of each node 10 of the 2D meshes representing the present-day geometry of the geological domain 1, as described in any Boundary Element Method known in the art.
In the following, an example of application of the Boundary Element Method to the computation of mechanical deformation of solid materials is given, based on an elastic deformation model. Only a brief overview is given, as the BEM method itself is of common knowledge.
This method assumes that the constitutive and equilibrium equations defining the mechanical problem involve two operators or matrices T(p, q) and U(p, q), where T(p, q) is the traction operator giving the tractions t at “field” point q due to a unit load at “source” point p, and U(p, q) is a Green's function giving the displacement d at “field” point q due to a unit load at “source” point p. For example, these two operators can be given by the fundamental solutions of linear elasticity. The Betti theorem gives the unknown displacement d at point p:
where S is the boundary of the considered solid. Using a discretization of S into polygonal elements, we obtain the system of equations:
MT·d=MU·t
with t the set of prescribed tractions, and MT and MU two dense and non-symmetric matrices. Therefore, the amount of computer memory needed to store the system of equations is in the order of n2 and the complexity of an iterative solver is in the order of k·n2, with n being the number of unknowns and k the number of iterations needed to solve the system. However, specific techniques such as multipole expansion (as described for instance in Geengard and Rokhlin, “A Fast Algorithm for Particle Simulations”, J. Comput. Phys., Vol. 73, 1987) or hierarchical matrices (as described for instance in Hackbusch, “A Sparse Matrix Arithmetic Based on H-Matrices—Part I—Introduction to H-Matrices”, Computing, Vol. 86, No 4, 1999) can be used to reduce the complexity in the order of n·log(n).
Finally, the solved nodal displacements (dx, dy, dz) are used to define the depositional coordinate system (u, v, w) at the location of each node 10 of the 2D meshes representing the present-day geometry of the geological domain 1 and/or at any observation point.
With reference to
More generally, the invention includes also a post-processing method for “painting” the elements of a 3D mesh 9 representing the present-day geometry of a geological domain 1 with properties 30 simulated in the depositional space which has been produced by the method according to the invention. It comprises the following steps:
Illustrations of this process can be found on
Any method embodiment (or portion thereof) described herein may be implemented in terms of program instructions. The program instructions may be stored on any of various kinds of readable computer memory media. The program instructions may by readable and executable by any computer or set of computers to implement the method embodiment (or portion thereof).
While this invention has been described in conjunction with a number of embodiments, it is evident that many alternatives, modifications and variations would be or are apparent to those of ordinary skill in the applicable arts. Accordingly, it is intended to embrace all such alternatives, modifications, equivalents and variations that are within the spirit and scope of this invention, the essential point being to use the theory of mechanics and rock deformation to build a depositional space, advantageously, (i) where isochrons (horizons) identified within a 3D geological domain are substantially planar and parallel, (ii) such that every point of the depositional space located inside a stratigraphic sequence has one and exactly one corresponding location in the present-day space, and (iii) where the geometry of the 3D geological domain is physically reliable, said depositional space being thus used as a computational space for applications like accurate modeling and interpolation of the physical properties of the subsurface, mesh generation for complex geometries of geological structures, or quality checking of a 3D seismic reflection interpretation.
Patent | Priority | Assignee | Title |
10036829, | Sep 28 2012 | ExxonMobil Upstream Research Company | Fault removal in geological models |
10087721, | Jul 29 2010 | Young Living Essential Oils, LC | Methods and systems for machine—learning based simulation of flow |
10088596, | Mar 15 2013 | Schlumberger Technology Corporation | Meshless representation of a geologic environment |
10107938, | Oct 31 2014 | ExxonMobil Upstream Research Company | Managing discontinuities in geologic models |
10319143, | Jul 30 2014 | ExxonMobil Upstream Research Company | Volumetric grid generation in a domain with heterogeneous material properties |
10359523, | Aug 05 2014 | ExxonMobil Upstream Research Company | Exploration and extraction method and system for hydrocarbons |
10803534, | Oct 31 2014 | ExxonMobil Upstream Research Company | Handling domain discontinuity with the help of grid optimization techniques |
10839114, | Dec 23 2016 | ExxonMobil Upstream Research Company | Method and system for stable and efficient reservoir simulation using stability proxies |
11163079, | Oct 07 2016 | ELIIS | Method for producing a geological vector model |
11409023, | Oct 31 2014 | ExxonMobil Upstream Research Company | Methods to handle discontinuity in constructing design space using moving least squares |
11662501, | Apr 03 2018 | EXXONMOBIL TECHNOLOGY AND ENGINEERING COMPANY | Geologic modeling methods and systems having constrained restoration of depositional space |
12099159, | Apr 15 2019 | Schlumberger Technology Corporation | Modeling and simulating faults in subterranean formations |
8552755, | May 12 2010 | Intermolecular, Inc. | High throughput current-voltage combinatorial characterization tool and method for combinatorial solar test substrates |
8674984, | Jul 16 2010 | IFP Energies Nouvelles | Method for generating a hex-dominant mesh of a geometrically complex basin |
9058445, | Jul 29 2010 | ExxonMobil Upstream Research Company | Method and system for reservoir modeling |
9058446, | Sep 20 2010 | ExxonMobil Upstream Research Company | Flexible and adaptive formulations for complex reservoir simulations |
9134454, | Apr 30 2010 | ExxonMobil Upstream Research Company | Method and system for finite volume simulation of flow |
9187984, | Jul 29 2010 | ExxonMobil Upstream Research Company | Methods and systems for machine-learning based simulation of flow |
9361414, | Aug 10 2012 | Schlumberger Technology Corporation | Hybrid local nonmatching method for multiphase flow simulations in heterogeneous fractured media |
9489176, | Sep 15 2011 | ExxonMobil Upstream Research Company | Optimized matrix and vector operations in instruction limited algorithms that perform EOS calculations |
9600608, | Aug 17 2011 | IFP Energies Nouvelles | Method of constructing a geological model comprising setting a depositional position of stratigraphic units |
Patent | Priority | Assignee | Title |
7043410, | Apr 04 2000 | ConocoPhillips Company; NORSKE CONOCO, A NORWEGIAN CORP | Method of modeling of faulting and fracturing in the earth |
7123258, | Dec 10 2001 | Paradigm France | Method, device and program for three-dimensional modeling of a geological volume by 3D parametering of the geological domain |
7424415, | Jun 06 2001 | Schlumberger Technology Corporation | Automated system for modeling faulted multi-valued horizons |
7480205, | Apr 20 2005 | Landmark Graphics Corporation | 3D fast fault restoration |
7603265, | Apr 14 2004 | Institut Francais du Petrole | Method of constructing a geomechanical model of an underground zone intended to be coupled with a reservoir model |
20020091502, | |||
20030216897, | |||
20090265152, | |||
EP2317348, | |||
WO3050766, | |||
WO2005119304, | |||
WO2006113939, |
Executed on | Assignor | Assignee | Conveyance | Frame | Reel | Doc |
Nov 20 2009 | Schlumberger Technology Corporation | (assignment on the face of the patent) | / | |||
Dec 15 2009 | LEPAGE, FRANCOIS | IGEOSS | ASSIGNMENT OF ASSIGNORS INTEREST SEE DOCUMENT FOR DETAILS | 023864 | /0428 | |
Nov 03 2010 | IGEOSS | Schlumberger Technology Corporation | ASSIGNMENT OF ASSIGNORS INTEREST SEE DOCUMENT FOR DETAILS | 025239 | /0014 |
Date | Maintenance Fee Events |
Apr 22 2016 | STOL: Pat Hldr no Longer Claims Small Ent Stat |
May 05 2016 | M1551: Payment of Maintenance Fee, 4th Year, Large Entity. |
May 07 2020 | M1552: Payment of Maintenance Fee, 8th Year, Large Entity. |
May 08 2024 | M1553: Payment of Maintenance Fee, 12th Year, Large Entity. |
Date | Maintenance Schedule |
Nov 20 2015 | 4 years fee payment window open |
May 20 2016 | 6 months grace period start (w surcharge) |
Nov 20 2016 | patent expiry (for year 4) |
Nov 20 2018 | 2 years to revive unintentionally abandoned end. (for year 4) |
Nov 20 2019 | 8 years fee payment window open |
May 20 2020 | 6 months grace period start (w surcharge) |
Nov 20 2020 | patent expiry (for year 8) |
Nov 20 2022 | 2 years to revive unintentionally abandoned end. (for year 8) |
Nov 20 2023 | 12 years fee payment window open |
May 20 2024 | 6 months grace period start (w surcharge) |
Nov 20 2024 | patent expiry (for year 12) |
Nov 20 2026 | 2 years to revive unintentionally abandoned end. (for year 12) |