Imaging a subsurface geological model at a past intermediate restoration time

Abstract

A system and method is provided for restoring a 3D tomographic model of the Earth's subsurface geology from the present-day to a past restoration time. Whereas at the present time all faults represent active discontinuities, at a past restoration time some faults have not yet formed. Accordingly, the restored model divides the fault network into τ-active faults (discontinuous surfaces for faults that intersect the layer deposited at the past restoration time) and τ-inactive faults (continuous surfaces for faults that do not intersect the layer deposited at the past restoration time). A new 3D restoration transformation is also provided that uses linear geological constraints to process the restoration model in less time and generate more accurate geological images.

Claims

1. A system for decompacting a 3D model of the subsurface geology of the Earth at an intermediate restoration time in the past τ, the system comprising: one or more processors configured to: receive a 3D model of present-day geometry of the subsurface geology imaging geological data produced in a three-dimensional field by tomographic scanning using a transmitter emitting a series of incident waves that reflect at geological discontinuities throughout the subsurface geology of the Earth and a receiver receiving reflections of those waves to image the subsurface geology of the Earth, receive a measure of present-day porosity experimentally measured within the subsurface geology of the Earth, select a value of a restoration time in the past τ before the present-day and after a time an oldest horizon surface in the 3D model of the subsurface was deposited, restore the 3D model from the present: day measured geometry to the predicted past geometry at the restoration time in the past τ using a 3D transformation, and decompact the vertical dimension of the restored 3D model to elongate vertical lengths of geological layers below a horizon layer deposited at the restoration time in the past τ, wherein the vertical lengths are elongated based on a relationship between a depositional porosity of the geological layers at the time sediment in those layers was deposited, restoration porosity of the geological layers at the restoration time in the past τ, and the present-day porosity of the geological layers experimentally measured in the present-day.

2. The system of claim 1, wherein the one or more processors are configured to decompact the vertical dimension of the restored 3D model by a combination of total decompaction corresponding to an increase in porosity from the present day porosity to the depositional porosity and partial recompaction corresponding to a partial decrease in the porosity from the depositional porosity to the restored porosity.

3. The system of claim 1, wherein the restoration porosity is greater than the present-day porosity and less than the depositional porosity.

4. The system of claim 1, wherein the one or more processors are configured to elongate the vertical length as a measure of height in the vertical dimension.

5. The system of claim 1, wherein the one or more processors are configured to elongate the vertical length as a measure of geological-time when the particles of sediment were originally deposited on the Earth's surface.

6. The system of claim 1, wherein the one or more processors are configured to compute the relationship between the depositional porosity, the restoration porosity, and the present-day porosity as: $\frac{1-{\overline{\varphi }}_{\tau }^{\oplus }}{1-{\stackrel{\_}{\varphi }}_{\tau }^{\ominus }},$ where ϕ.sub.τ.sup.⊕=ψ.sub.0−ψ and ϕ.sub.τ.sup.⊖=ψ.sub.0−ψ.sub.p for all points in the restored 3D model.

7. The system of claim 1 comprising a probe configured to extract material from within the subsurface geology of the Earth or one or more wells and experimentally measure the present-day porosity of the extracted material.

8. The system of claim 1, wherein the one or more processors are configured to use a boundary condition that ensures that a top horizon deposited at the restoration time τ is a horizontal plane in the restored 3D model.

9. The system of claim 1, wherein the one or more processors are configured to use a boundary condition that ensures that a direction of change of geological-time when the particles of sediment were originally deposited on the Earth's surface is vertical in the restored 3D model.

10. The system of claim 1, wherein the one or more processors are configured to use a boundary condition that ensures that, for any pair of collocated points on opposite sides of a fault, the two collocated points are decompacted to have the same coordinate.

11. The system of claim 1, wherein the one or more processors are configured to decompact the vertical dimension of the restored 3D model by: computing an elongated geological-time in the restored 3D model, transforming the elongated geological-time from the restored 3D model to generate a 1D geological-time in the present-day 3D model, computing 2D paleo-depositional coordinates based on the transformed geological-time in the present-day 3D model, and performing a 3D transformation comprising the 1D geological-time and 2D paleo-depositional coordinates from the present-day 3D model to the restored 3D model that is decompacted based on the elongated geological-time.

12. The system of claim 1, wherein the one or more processors are configured to iteratively decompact the subsurface layer-by-layer, starting at the top horizon deposited at the restoration time τ and ending at the bottom horizon deposited at the depositional time.

13. The system of claim 12, wherein the depositional porosity and the present-day porosity is independently determined for each geological layer of the subsurface.

14. The system of claim 1, wherein, when the depositional porosity and the present-day porosity are substantially constant throughout the subsurface geology, the one or more processors are configured to decompact in one operation over the entire domain of the restored 3D model.

15. A method for decompacting a 3D model of the subsurface geology of the Earth at an intermediate restoration time in the past τ, the method comprising: receiving a 3D model of present-day geometry of the subsurface geology imaging geological data produced in a three-dimensional field by tomographic scanning using a transmitter emitting a series of incident waves that reflect at geological discontinuities throughout the subsurface geology of the Earth and a receiver receiving reflections of those waves to image the subsurface geology of the Earth; receiving a measure of present-day porosity experimentally measured within the subsurface geology of the Earth; selecting a value of a restoration time in the past τ before the present-day and after a time an oldest horizon surface in the 3D model of the subsurface was deposited; restoring the 3D model from the present day measured geometry to the predicted past geometry at the restoration time in the past τ using a 3D transformation; and decompacting the vertical dimension of the restored 3D model to elongate vertical lengths of geological layers below a horizon layer deposited at the restoration time in the past τ, wherein the vertical lengths are elongated based on a relationship between a depositional porosity of the geological layers at the time sediment in those layers was deposited, restoration porosity of the geological layers at the restoration time in the past τ, and the present-day porosity of the geological layers experimentally measured in the present-day.

16. The method of claim 15 comprising decompacting the vertical dimension of the restored 3D model by a combination of total decompaction corresponding to an increase in porosity from the present-day porosity to the depositional porosity and partial recompaction corresponding to a partial decrease in the porosity from the depositional porosity to the restored porosity.

17. The method of claim 15, wherein the restoration porosity is greater than the present-day porosity and less than the depositional porosity.

18. The method of claim 15 comprising elongating the vertical length as a measure of height in the vertical dimension.

19. The method of claim 15 comprising elongating the vertical length as a measure of geological-time when the particles of sediment were originally deposited on the Earth's surface.

20. The method of claim 15 comprising computing the relationship between the depositional porosity, the restoration porosity, and the present-day porosity as: $\frac{1-{\overline{\varphi }}_{\tau }^{\oplus }}{1-{\stackrel{\_}{\varphi }}_{\tau }^{\ominus }},$ where ϕ.sub.τ.sup.⊕=ψ.sub.0−ψ and ϕ.sub.τ.sup.⊖=ψ.sub.0−ψ.sub.p for all points in the restored 3D model.

21. The method of claim 15 comprising extracting material from within the subsurface geology of the Earth or one or more wells and experimentally measuring the present-day porosity of the extracted material.

22. The method of claim 15 comprising applying a boundary condition that ensures that a top horizon deposited at the restoration time τ is a horizontal plane in the restored 3D model.

23. The method of claim 15 comprising applying a boundary condition that ensures that a direction of change of geological-time when the particles of sediment were originally deposited on the Earth's surface is vertical in the restored 3D model.

24. The method of claim 15 comprising applying a boundary condition that ensures that, for any pair of collocated points on opposite sides of a fault, the two collocated points are decompacted to have the same coordinate.

25. The method of claim 15 comprising decompacting the vertical dimension of the restored 3D model by: computing an elongated geological-time in the restored 3D model; transforming the elongated geological-time from the restored 3D model to generate a 1D geological-time in the present-day 3D model; computing 2D paleo-depositional coordinates based on the transformed geological-time in the present-day 3D model; and performing a 3D transformation comprising the 1D geological-time and 2D paleo-depositional coordinates from the present-day 3D model to the restored 3D model that is decompacted based on the elongated geological-time.

26. The method of claim 15 comprising iteratively decompacting the subsurface layer-by-layer, starting at the top horizon deposited at the restoration time τ and ending at the bottom horizon deposited at the depositional time.

27. The method of claim 26 comprising independently determining the depositional porosity and the present-day porosity for each geological layer of the subsurface.

28. The method of claim 15 comprising, when the depositional porosity and the present-day porosity are substantially constant throughout the subsurface geology, decompacting in one operation over the entire domain of the restored 3D model.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) The principles and operation of the system, apparatus, and method according to embodiments of the present invention may be better understood with reference to the drawings, and the following description, it being understood that these drawings are given for illustrative purposes only and are not meant to be limiting.

(2) FIG. 1 schematically illustrates an exploded view of a 3D present day geological model of the subsurface of the Earth according to an embodiment of the invention. The 3D geological model may comprise a faulted 3D grid Γ 100. Cell edges 106 of the grid are constrained to never cross faults 105. During restoration, twin faces F.sup.− 104 and F.sup.+ 103 on opposite sides of a fault F 105 may slide along one another only in a direction tangential to the surface of the fault F 105. Points (r.sub.F.sup.+, r.sub.f.sup.−) (101,102) are twin-points. “Twin” points or faces may refer to points or faces that were collocated at the time of their deposition, but which may have separated at a later time.

(3) FIG. 2 schematically illustrates a vertical cross section of a 3D subsurface model in a volume deformation: the u.sub.τ v.sub.τ t.sub.τ—transform 201 restores the volume G.sub.τ 202 in a present day space G 220 to a past volume G.sub.τ 203 in a past restored space 219 as the subsurface was configured at restoration time τ according to an embodiment of the invention.

(4) FIG. 3 schematically illustrates a vertical cross section of the present day domain G.sub.τ 202 according to an embodiment of the invention. Horizons 210 are shown in bold lines and level sets 208 of the vertical restoration coordinate t.sub.τ(r) are shown in dashed lines. τ-active faults 105 (depicted as bold black lines) cut the horizon H.sub.τ 210 deposited at restoration time τ whilst τ-inactive faults 300 (depicted as bold gray lines) do not cut horizon H.sub.τ 210. Level sets of the vertical restoration coordinate t.sub.τ(r) (dashed lines 208) are continuous across τ-inactive faults 300 and only separated across τ-active faults 105 because the sedimentary layers of those level sets were not divided when the faults were inactive (at a time before the faults formed).

(5) FIG. 4 schematically illustrates an example 3D geological model used to test the accuracy of the vertical restoration coordinate t.sub.τ(r) according to an embodiment of the invention. This example 3D geological model has a “ramp” geological structure with a restored horizon H.sub.τ 210 as the central, sigmoid surface. This seemingly “simple” test comprises highly complex calculations because there are large variations in layer thickness which make numerical computation of the geological time of deposition t(r) more challenging.

(6) FIG. 5 shows a comparison of multiple histograms 501, 502 and 503 of ∥grad t.sub.τ(r)∥ each using a different method to compute t.sub.τ(r) in the test case domain G.sub.τ of FIG. 4. Depending on the method used to compute t.sub.τ(r), the resulting magnitude of grad t.sub.τ(r) may severely deviate from the ideal value of “1,” as required by equation (10). Histogram 503, which is generated according to an inventive embodiment, is the closest approximation of (10), yielding the most accurate restoration model.

(7) FIG. 6 schematically illustrates a vertical cross section of a 3D subsurface model in which a pair of τ-twin-points (r.sub.F.sup.+, r.sub.F.sup.−).sub.τ (601,602) are deduced from a pair of present day twin-points (r.sub.F.sup.+, r.sub.F.sup.−) (101,102), using a depositional (e.g., GeoChron) model as input, according to an embodiment of the invention.

(8) FIG. 7 schematically illustrates a direct uvt-transform 700 and inverse uvt-transform 701 that transform a 3D model between a present day geological space G 220 and a depositional geological space G 719 according to an embodiment of the invention.

(9) FIG. 8 schematically illustrates a u.sub.τ v.sub.τ t.sub.τ—transform 201 of a vertical cross-section of a 3D subsurface model showing the geological impact between fault 300 being erroneously considered a τ-active fault (top image of FIG. 8) vs. correctly considered a τ-inactive fault with respect to restored horizon H.sub.τ 210 (bottom image of FIG. 8). In the top image of FIG. 8, when fault 300 is erroneously considered a τ-active fault, in order to preserve geological volume, present day fault block 800 is transformed to a restored fault block 801 that intersects τ-active fault 805, which contradicts geological rules. In contrast, in the bottom image of FIG. 8, when fault 300 is correctly considered a τ-inactive fault, a volume-preserving transformation maps the present day fault block 800 correctly to a restored fault block that stays within (and does not cross) τ-active fault 805, according to geological rules.

(10) FIG. 9 shows a comparison of multiple histograms 901, 902 and 903 of ΔV/V each using a different method to compute t.sub.τ(r) in the test case domain G.sub.τ of FIG. 4. Depending on the method used to compute t.sub.τ(r), the resulting magnitude of ΔV/V may severely deviate from the target value of “0”. Histogram 903, which is generated according to an embodiment of the invention, has the least amount of volume variation, yielding the most accurate restoration model.

(11) FIG. 10 schematically illustrates an example of a 3D u.sub.τ v.sub.τ t.sub.τ—transform 201 of a horizon H.sub.τ 210 (depicted as a white layer) from a present day domain G 202 to a restored domain G.sub.τ 203 at restoration time τ according to an embodiment of the invention.

(12) FIG. 11 schematically illustrates an example of a sequence of restorations at geological times {τ.sub.1<τ.sub.2< . . . <τ.sub.n} of vertical cross sections of the 3D model of FIG. 10 according to an embodiment of the invention.

(13) FIG. 12 schematically illustrates an example of fault striae induced on faults by paleo-geographic coordinates u(r) and v(r) of a depositional model provided as input to the restoration according to an embodiment of the invention.

(14) FIG. 13 schematically illustrates an example display of a sequential chronological transformation of a vertical cross section of a 3D subsurface model according to an embodiment of the invention. Chronological transformation may correspond to a plurality of restoration times τi ordered in a sequence either “forward modeling” according to a forward passage of time (e.g., with sequentially ascending or later values of time) or “restoration” according to a reverse passage of time (e.g., with sequentially descending or earlier values of time). In the example shown in FIG. 13, the sequence of chronological transformation progresses according to a forward passage of time, from the start of deposition at earliest geological time τ.sub.1, to one or more subsequent intermediate restorations times τ.sub.2<τ.sub.3, to a latest present day geological time τ.sub.4 (although the sequence may progress according to the reverse passage of time). Further, any other number or orders of restorations or times may be used).

(15) FIG. 14 schematically illustrates a geological tomography technique in which a series of incident and reflected waves are propagated through a subsurface region of the Earth to image the subsurface according to an embodiment of the invention.

(16) FIG. 15 schematically illustrates a system for restoring a present day geological model to an intermediate restoration time τ, according to an embodiment of the invention.

(17) FIG. 16 is a flowchart of a method to restore a geological model with τ-active and τ-inactive faults, according to an embodiment of the invention.

(18) FIG. 17 is a flowchart of a method to restore a geological model with improved accuracy using a new thickness-preserving constraint, according to an embodiment of the invention.

(19) FIG. 18 schematically illustrates an example 3D geological volume of a model that is compacted (right image) and that is decompacted (left image) at an intermediate restoration time in the past τ, according to an embodiment of the invention.

(20) FIG. 19 is a flowchart of a method for decompacting a 3D model of the subsurface geology of the Earth at an intermediate restoration time in the past τ, according to an embodiment of the invention.

(21) For simplicity and clarity of illustration, elements shown in the drawings have not necessarily been drawn to scale. For example, the dimensions of some of the elements may be exaggerated relative to other elements for clarity. Further, where considered appropriate, reference numerals may be repeated among the drawings to indicate corresponding or analogous elements throughout the serial views.

DETAILED DESCRIPTION OF EMBODIMENTS OF THE INVENTION

(22) Embodiments of the invention improve conventional restoration techniques for imaging restored geological models as follows:

(23) “τ-Active” Faults Vs. “τ-Inactive” Faults:

(24) In conventional restoration models, all faults are active (as discontinuous surfaces) at all times. However, in reality, certain faults have not yet formed or activated at various intermediate restoration times T. Accordingly, conventional restoration models generate false or “phantom” faults that erroneously divide geology that has not yet fractured, leading to geological inaccuracies in subsurface images.

(25) Embodiments of the invention solve this problem by selectively activating and deactivating individual fault surfaces to be discontinuous or continuous, respectively, depending on the specific restoration geological-time τ. For each intermediate restoration time in the past τ, embodiments of the invention split faults into two complementary subsets of “τ-active” faults and “τ-inactive” faults. τ-active faults are activated at restoration time τ (e.g., a discontinuous fault surface along which fault blocks slide tangentially), whereas τ-inactive faults are deactivated at restoration time τ (e.g., a continuous surface that does not behave as a fault).

(26) As faults form and evolve over time, they behave differently at different geological times in the past. For example, a fault that forms at an intermediate geological-time τ, where τ.sub.1<τ<τ.sub.2, is τ-active in a restored model at later time τ.sub.2 (after the fault has formed), but τ-inactive in a restored model at earlier time τ.sub.1 (before the fault has formed). This fault classification allows faults to be modelled differently at each restoration time τ in a geologically consistent way, thereby preventing unrealistic deformations from being generated in the neighborhood of these faults.

(27) FIG. 8 shows the problem of a fault 300 being erroneously considered active at a restoration time before it formed (top image of FIG. 8) and the solution of modeling the fault as a τ-inactive fault to correctly deactivate the fault at restoration time τ according to embodiments of the invention (bottom image of FIG. 8). In the top image of FIG. 8, when a horizon H.sub.τ 210 is restored using u.sub.τ v.sub.τ t.sub.τ-transform 201, fault block 800 (shaded region in the top-left image of FIG. 8) is bounded by an active fault 105 and an inactive fault 300. If however fault 300 is erroneously considered as an active fault then, after applying restoration 201: erroneous τ-twin points (803,823) are transformed into a pair of collocated points 813=833, τ-twin points (804,824) are transformed into collocated points 814=834. It is clear that, if d(a, b) denotes the distance between any arbitrary pair of points (a, b), then:
d(803,824)=d(833,835) d(833,834)  (2) This observation shows that erroneously considering fault 300 as a τ-active fault inevitably generates unrealistic deformations.

(28) This problem is solved according to embodiments of the invention, e.g., as shown in the bottom image of FIG. 8. In this image, fault 105 is modeled as a τ-active fault (activating the fault), but fault 300 is modeled as a τ-inactive fault (deactivating the fault). Accordingly, when u.sub.τ v.sub.τ t.sub.τ-transform 201 is applied to fault block 800 (shaded region in the bottom-left image of FIG. 8), restored fault block 801 (shaded region in the bottom-right image of FIG. 8) is no longer bounded by an active fault (τ-inactive fault 300 is inactive at restored time τ). Accordingly, the restored fault block 801 preserves volume and stays within (and does not cross) τ-active fault 805 (because the deactivated boundary transformed from τ-inactive fault 300 may shift to accommodate a shift in the restored τ-active fault 805).

(29) Contrary to conventional methods, the use of τ-active and τ-inactive faults produces more accurate results, e.g., even if there is no continuous path between (no way to continuously connect) a given fault block (e.g., 800) and the horizon Hτ (e.g., 210) deposited at geological time τ, which typically requires additional processing that may induce errors. By selectively activating and inactivating faults at the various restoration times according to when they form, embodiments of the invention eliminate erroneous phantom faults and more accurately represent the faulted geology.

(30) Reference is made to FIG. 16, which is a flowchart of a method to restore a geological model using τ-active and τ-inactive faults, according to an embodiment of the invention.

(31) In operation 1610, a processor may receive a 3D model of the present day measured geometry comprising a network of faults (e.g., present day model 202). The present day model may be measured tomographically by scanning the Earth's subsurface e.g., as described in reference to FIGS. 14 and 15. In the present day geology, all faults in the model have already formed and so, represent active discontinuities that divide fault blocks which slide in opposite directions tangential to the surface of the fault as time approaches a modeled time.

(32) In operation 1620, a processor may select or receive a past restoration time τ that is “intermediate” or prior to the present time and after the start of the subsurface's deposition (the time period when an oldest horizon surface in the 3D model was originally deposited).

(33) In operation 1630, a processor may divide the network of faults into a subset of τ-active faults and a subset of τ-inactive faults. τ-active faults may be faults that are active at the past restoration time τ and τ-inactive faults are faults that are inactive at the past restoration time τ. A fault is determined to be τ-active when the fault intersects a horizon H.sub.τ that was originally deposited at the past restoration time τ (e.g., see τ-active faults 105 of FIG. 3) and a fault is determined to be τ-inactive when the fault does not intersect the horizon H.sub.τ that was originally deposited at the past restoration time τ (e.g., see τ-inactive faults 300 of FIG. 3). Because different faults activate to fracture the subsurface at different geological times, the processor may divide the fault network differently at different geological times. Thus, a fault may be τ-active at a first restoration time τ′ (e.g., a time period during which the fault has formed) and τ-inactive at a second restoration time τ′ (e.g., a time period different than that during which the fault has formed). In one embodiment, iso-value surfaces (e.g., 208 of FIG. 3) of each restoration coordinate (e.g., u.sub.τ, v.sub.τ, and t.sub.τ) are continuous across τ-inactive faults (e.g., 300 of FIG. 3) and discontinuous across τ-active faults (e.g., 105 of FIG. 3).

(34) In operation 1640, a processor may restore the 3D model from the present day measured geometry to the predicted past geometry at the past restoration time τ. During restoration, the processor may flatten a horizon H.sub.τ (e.g., 210 of FIG. 4) that was originally deposited at time τ to a substantially planar surface of approximately constant depth. For horizons older (e.g., deposited deeper in the subsurface) than horizon H.sub.τ, the processor may restore the horizons to non-planar surfaces, e.g., when the thickness of the layers is not constant. Because the region of the subsurface deposited after the restoration time τ (e.g., deposited shallower in the subsurface) did not yet exist at the time of the restored model, restoring the 3D model to a past restoration time τ may eliminate (e.g., removing or not displaying) all relatively shallower horizon surfaces that were originally deposited after the past restoration time τ. During restoration, the processor may treat τ-active and τ-inactive faults differently in operations 1650 and 1660, respectively.

(35) In operation 1650, for each τ-active fault, a processor may model the τ-active fault as an active discontinuous fault surface and restore the horizon surface by removing or omitting the fault surface at the time of restoration. The processor may eliminate the τ-active fault during restoration by sliding its adjacent fault blocks together. This may join end points of a horizon H.sub.τ separated on opposite sides of the fault in the present day model to merge into the same position in the restored model by sliding the end points towards each other in a direction tangential to the surface of the τ-active fault.

(36) In operation 1660, for each τ-inactive fault, a processor may model the τ-inactive fault, not as a discontinuous fault surface, but as a continuous non-fault surface in the restoration transformation. The τ-inactive fault may be modeled as a surface in which the discontinuity induced by the fault has been deactivated to prevent fault blocks from sliding in directions tangential to the surface of the fault as time approaches the restoration time τ. The processor may model the τ-inactive fault during restoration by keeping collocated points on opposite sides of the fault in the present day model together in the restored model.

(37) After the geological model has been restored for a first past restoration time r (operations 1620-1660), the process may repeat to restore the model for a second different past restoration time τ′. In some embodiments, the geological model may be sequentially restored to a sequence of multiple past restoration times τ.sub.1, τ.sub.2, . . . , τ.sub.n. In multiple (all or not all) of the past restoration times τ.sub.1, τ.sub.2, . . . , τ.sub.n, the fault network may be divided into a different subset of τ-active and τ-inactive faults, e.g., because different faults fracture the subsurface at different geological times. In some embodiments, a processor may play a moving image sequence in which the 3D model is iteratively restored in a forward or reverse order of the sequence of past restoration times τ.sub.1, τ.sub.2, . . . , τ.sub.n to visualize changes in the subsurface geology over the passage of time.

(38) In operation 1670, a processor may display a visualization of an image of the subsurface geology of the Earth overlaid with τ-active faults and τ-inactive faults in the restored model at past restoration time τ. The processor may display the τ-active faults and the τ-inactive faults with different visual identifiers, such as, different levels of translucency, different colors, different patterns, etc.

(39) New Restoration Transformation u.sub.τ, v.sub.τ, and t.sub.τ:

(40) A restoration transformation may transform a geological image of the subsurface of the Earth from a present day space (e.g., x,y,z coordinates) to a restoration space (e.g., u.sub.τ, v.sub.τ, and t.sub.τ coordinates) as it was formed at an intermediate restoration time in the past r (before the present-day but after the start of the subsurface deposition). An ideal restoration should transform the vertical coordinate t.sub.τ in a manner that strictly honors the thickness of layers, to preserve areas and volumes of the Earth, so that terrains are not stretched or squeezed over time in the vertical dimension. However, conventional restoration transformations typically deform the vertical coordinates, forcing terrains to stretch and squeeze, resulting in errors in the restoration model.

(41) Embodiments of the invention improve the accuracy of the restoration model by establishing a vertical restoration coordinate t.sub.τ that preserves layer thickness. This may be achieved by implementing a thickness-preserving constraint that sets a variation in the vertical restoration coordinate t.sub.τ between any two points separated by an infinitesimal difference in the direction of maximal variation of the vertical coordinate t.sub.τ to be approximately equal to the distance between the points in the direction of maximal variation. An example of this constraint may be modeled by ∥grad t.sub.τ(x,y,z)∥=1. This constraint, however, is non-linear and highly complex and time-consuming to solve. Due to its complexity, this constraint is rarely used in conventional restoration models, and instead replaced by over-simplifications, such as equations (33) and (34), that result in model errors as shown in histograms 501 and 502 of FIG. 5, and histograms 901 and 902 of FIG. 9, respectively.

(42) Embodiments of the invention improve the accuracy of the restored model by establishing a new thickness-preserving constraint that introduces an error correction term ε.sub.τ. The new thickness-preserving constraint sets the restored vertical coordinate t.sub.τ to be equal to a sum of a first approximation t′.sub.τ of the vertical coordinate and an error correction term ε.sub.τ, wherein the error correction term ε.sub.τ is computed by solving a relationship in which a variation in the sum of the first approximation t′.sub.τ of the vertical coordinate and the error correction term ε.sub.τ between any two points separated by an infinitesimal difference in the direction of maximal variation of the sum is approximately equal to the distance between the points in the direction of maximal variation. An example of this constraint may be modeled by ∥grad (t′.sub.τ+ε.sub.τ)∥=1. The new thickness-preserving constraint preserves layer thickness with greater accuracy as shown in histogram 503 of FIG. 5 as compared to conventional approximations shown in histograms 501 and 502 of FIG. 5 and minimizes volume variation with greater accuracy as shown in histogram 903 of FIG. 9 as compared to conventional approximations shown in histograms 901 and 902 of FIG. 9, respectively.

(43) Embodiments of the invention further improve the performance and computational speed of the computer generating the restored model by linearizing the new thickness-preserving constraint. As an example, the new thickness-preserving constraint may be linearized as follows. ∥grad (t′.sub.τ+ε.sub.τ)∥=1 may be squared to obtain ∥grad t′.sub.τ∥.sup.2+∥grad ε.sub.τ∥.sup.2+∥2.Math.grad t′.sub.τ.Math.grad ε.sub.τ∥=1. The error correction term ε.sub.τ may be generated such that the square of its spatial variation, ∥grad ε.sub.τ∥.sup.2, is negligible. Accordingly, the thickness-preserving constraint simplifies to a new linear thickness-preserving constraint of grad ε.sub.τ.Math.grad t′.sub.τ≅½ {1−∥grad t′.sub.τ∥.sup.2} (eqn. (37)). This thickness-preserving constraint is linear because t′.sub.τ is already known, so the constraint is a relationship between the gradient of the error ε.sub.τ and the gradient of the known first approximation of the vertical coordinate t′.sub.τ. The computer may therefore compute the new thickness-preserving constraint in linear time, which is significantly faster than computing the non-linear constraints ∥grad t.sub.τ∥=1 or ∥grad (t′.sub.τ+ε.sub.τ)∥=1.

(44) Contrary to conventional methods, the computational complexity for performing the restoration transformation according to embodiments of the invention is significantly reduced compared to classical methods that are based on the mechanics of continuous media. As a consequence, the modeling computer uses significantly less computational time and storage space to generate the inventive restoration model.

(45) Contrary to conventional methods that allow variations of geological volumes and deformations, embodiments of the invention implement a new set of geometrical constraints and boundary conditions that preserve geological volumes and deformations while adhering to geological boundaries.

(46) Contrary to conventional methods, embodiments of the invention restore faults along fault striae (e.g., see FIG. 12) induced by the twin points associated with the paleo-geographic coordinates of a depositional (e.g., GeoChron) model, given as input of the restoration method.

(47) An ideal restoration should also transform the horizontal coordinates u.sub.τ and v.sub.τ in a manner that strictly honors lateral spatial distribution, to preserve areas and volumes of the Earth, so that terrains are not stretched or squeezed over time in the horizontal dimensions. However, conventional restoration transformations based on depositional coordinates (e.g., paleo-geographic coordinates u and v) typically deform the horizontal coordinates, forcing terrains to stretch and squeeze, resulting in errors in the restoration model.

(48) Embodiments of the invention improve the accuracy of the restoration model at time τ by establishing horizontal restoration coordinates u.sub.τ and v.sub.τ that restore the horizon surface H.sub.τ deposited at time τ consistently with horizontal depositional coordinates u and v whilst minimizing deformations. In one embodiment, on the horizon surface H.sub.τ only, the horizontal restoration coordinates u.sub.τ and v.sub.τ are equal to the depositional coordinates u and v (see e.g., equation (20)) and the spatial variations of the horizontal restoration coordinates u.sub.τ and y.sub.τ are preserved with respect to the horizontal depositional coordinates u and v (see e.g., equation (21)). Thus, each restoration model at time τ, presents a horizon surface H.sub.τ, as it was configured at that time τ when it was originally deposited. Additionally or alternatively, horizontal restoration coordinates u.sub.τ and v.sub.τ are modeled in a tectonic style (e.g., using constraints (22) or (23)) that is consistent with that of the horizontal coordinates u and v of the depositional model, which makes the restoration more accurate because the geological context is taken into account. Additionally or alternatively, horizontal restoration coordinates u.sub.τ and v.sub.τ are modeled to minimize deformations induced by the restoration of horizon H.sub.τ, rather than minimizing deformations in the whole volume G. This may be achieved by implementing constraints (41) and (42) that only enforce orthogonality of gradients of u.sub.τ and v.sub.τ with local axes b.sub.τ and a.sub.τ, but which do not constrain the norm of grad u.sub.τ and grad v.sub.τ, as is typically constrained for horizontal depositional coordinates u and v consistent with the depositional time model. Horizontal restoration coordinates u.sub.τ and v.sub.τ may also be constrained only in G.sub.τ, thereby only taking into account the part of the subsurface to be restored, not the entire model G. Additionally or alternatively, horizontal restoration coordinates u.sub.τ and v.sub.τ may be constrained to be equal on opposite sides of τ-active faults at twin point locations, where the twin points are computed from fault striae, which also ensures consistency with the depositional model (see e.g., equation (43)). Additionally or alternatively, horizontal restoration coordinates u.sub.τ and v.sub.τ are constrained to be equal on opposite sides of τ-inactive faults at mate point locations to cancel the effect of inactive faults on the restoration model (see e.g., equation (43)).

(49) Reference is made to FIG. 17, which is a flowchart of a method to restore a geological model with improved accuracy using a new thickness-preserving constraint, according to an embodiment of the invention.

(50) In operation 1710, a processor may receive a 3D model of the present day measured geometry (e.g., present day model 202) comprising one or more folded (e.g., curvilinear or non-planar) geological horizon surfaces (e.g., 210). The present day model may be measured tomographically by scanning the Earth's subsurface e.g., as described in reference to FIGS. 14 and 15.

(51) In operation 1720, a processor may select or receive a past restoration time τ that is “intermediate” or prior to the present time and after the start of the subsurface's deposition (the time period when an oldest horizon surface in the 3D model was originally deposited).

(52) In operation 1730, a processor may restore the 3D model from the present day measured geometry (e.g., present day model G.sub.τ 202 in xyz-space G 220) to the predicted past geometry at the restoration time in the past τ (e.g., restored model G.sub.τ 203 in u.sub.τ v.sub.τ t.sub.τ-space 219) using a 3D restoration transformation (e.g., u.sub.τ v.sub.τ t.sub.τ-transform 201). At the restored time in the past τ, the geological layers above H.sub.τ (e.g., H.sub.τ+1 . . . H.sub.n) did not yet exist, so the subregion above H.sub.τ in the present day space G 220 is eliminated or omitted, and only the subregion G.sub.τ 202 below and aligned with H.sub.τ (e.g., H.sub.1 . . . H.sub.τ) in the present day space G 220 is restored. The 3D restoration transformation includes a vertical component that restores the geometry to the vertical coordinate t.sub.τ and two lateral or horizontal components that restore the geometry to the horizontal coordinates u.sub.τ and v.sub.τ. The restored vertical coordinate t.sub.τ and horizontal coordinates u.sub.τ and v.sub.τ represent the predicted vertical and horizontal positions, respectively, where particles in the subsurface were located in the Earth at the restoration time in the past τ. Because the region of the subsurface deposited after the restoration time τ (e.g., deposited shallower in the subsurface than H.sub.τ) did not yet exist at the time of the restored model, the processor may restore and compute coordinates for the part or subregion G.sub.τ of the subsurface G that was deposited at a geological time of deposition t prior to or during the past restoration time τ (e.g., deposited deeper than, or at the same layer in the subsurface as, H.sub.τ). Accordingly, the restored model eliminates or omits all relatively shallower or younger horizon surfaces or layers that were originally deposited after the past restoration time τ.

(53) The processor may restore the vertical coordinate t.sub.τ such that points along a horizon surface H.sub.τ (e.g., 210) modeling sediment that was deposited at the selected restoration time τ have a substantially constant value for the restored vertical coordinate t.sub.τ (see e.g., eqn. (19)). Further, the processor may restore the vertical coordinate t.sub.τ such that at any location in the 3D model, the restored vertical coordinate t.sub.τ is equal to a sum of a first approximation t′.sub.τ of the vertical coordinate and an error correction term ε.sub.τ, wherein the error correction term ε.sub.τ is computed by solving a relationship in which a variation in the sum of the first approximation t′.sub.τ of the vertical coordinate and the error correction term ε.sub.τ between any two points separated by an infinitesimal difference in the direction of maximal variation of the sum is approximately equal to the distance between the points in the direction of maximal variation. The error correction term ε.sub.τ may correct errors in the first approximation t′.sub.τ of the vertical coordinate. This constraint may be represented by a linear second order approximation (see e.g., eqn. (37)).

(54) In some embodiments, the processor computes the first approximation t′.sub.τ of the vertical coordinate by solving a relationship in which the spatial variation of the vertical coordinate t′.sub.1 is locally approximately proportional to the spatial variation of a geological time of deposition t. In some embodiments, the coefficient of proportionality is locally equal to the inverse of the magnitude of the maximal spatial variation of the geological time of deposition (see e.g., eqn. (34)-(1)). This relationship may give the vertical restoration coordinate t.sub.τ the shape of the horizon H.sub.τ because, on the horizon, the gradient of depositional time t is normal to the horizon surface. Thus, the ratio grad t/∥grad t∥ follows the shape of the horizon.

(55) In some embodiments, the processor computes the first approximation t′.sub.τ of the vertical coordinate by solving a relationship in which any infinitesimal displacement in the direction orthogonal to horizon surface H.sub.τ results in a variation of the vertical coordinate t′.sub.τ approximately equal to the length of the infinitesimal displacement for points on the horizon surface H.sub.τ (see e.g., eqn. (33)-1)).

(56) In some embodiments, the processor computes the restored vertical coordinate t.sub.τ in parts of the subsurface which are older than restoration time τ such that iso-value surfaces of the restored vertical coordinate t.sub.τ are parallel to the horizon surface H.sub.τ and the difference in the restored vertical coordinate t.sub.τ between two arbitrary iso-values is equal to the distance between the corresponding iso-surfaces (see e.g., eqn. (31)). Parallel surfaces may be planar parallel in the restored model, and curved parallel (e.g., having parallel tangent surfaces) in present day model, such that the surfaces are non-intersecting at limits.

(57) In some embodiments, the error correction term ε.sub.τ is null at points along the horizon surface H.sub.τ that was deposited at the selected restoration time in the past τ so that the restored horizon surface H.sub.τ is flat (see e.g., eqn. (36)).

(58) In some embodiments, the restored horizontal coordinates u.sub.τ and v.sub.τ are constrained such that for each point along the horizon surface H.sub.τ that was deposited at the selected restoration time in the past τ: the restored horizontal coordinates u.sub.τ and v.sub.τ are equal to depositional horizontal coordinates u and v, respectively, and the spatial variations of the restored horizontal coordinates u.sub.τ and v.sub.τ are equal to the spatial variations of the depositional horizontal coordinates u and v, respectively (see e.g., eqns. (20)(21)). On average, globally over the entire model, the processor may compute ∥grad u∥=1 and ∥grad v∥=1. However, locally, this is not necessarily true e.g., on horizon Hr. So, while the processor sets grad u.sub.τ=grad u and grad v.sub.τ=grad v on Hτ, the processor may not constrain ∥grad u.sub.τ∥=1 and ∥grad v.sub.τ∥=1 on HT. Moreover, the processor may not constrain grad u.sub.τ to be orthogonal to grad t.sub.τ. This results from the boundary condition on Hτ and propagation through its constant gradient.

(59) In some embodiments, the restored horizontal coordinates u.sub.τ and v.sub.τ are constrained in parts of the subsurface which are older than restoration time τ such that directions of maximal change of the restored horizontal coordinates u.sub.τ and v.sub.τ are linearly constrained by a local co-axis vector b.sub.τ and a local axis vector a.sub.τ, respectively (see e.g., eqn. (41)).

(60) In some embodiments, the local axis vector a.sub.τ is oriented approximately in the direction of maximal change of depositional horizontal coordinate u and orthogonal to the direction of maximal change of the vertical restoration coordinate t.sub.τ, and the local co-axis vector b.sub.τ is oriented orthogonal to the direction of the local axis vector a.sub.τ and orthogonal to the direction of maximal change of the vertical restoration coordinate t.sub.τ (see e.g., eqn. (40)).

(61) In some embodiments, if the tectonic style of the 3D model is minimal deformation, the restored horizontal coordinates u.sub.τ and v.sub.τ are computed over the part of the 3D model of the subsurface which is older than restoration time τ such that the directions of maximal change of u.sub.τ and v.sub.τ are approximately orthogonal to the local co-axis vector b.sub.τ and the local axis vector a.sub.τ, respectively. For example, equation (40) constrains the local axis vector a.sub.t-to be parallel to the gradient of u and the local co-axis vector b.sub.τ to be orthogonal to the local axis vector a.sub.τ, which means that the gradient of u is orthogonal to the local co-axis vector b.sub.τ. Equation (41) further constrains the gradient of U.sub.τ to be approximately orthogonal to the local co-axis vector b.sub.τ. Accordingly, the gradient of u.sub.τ is approximately parallel to the gradient of u. The same logic implies the gradient of v.sub.τ is approximately parallel to the gradient of v.

(62) In some embodiments, if the tectonic style of the 3D model is flexural slip, the restored horizontal coordinates u.sub.τ and v.sub.τ are computed over the part of the 3D model of the subsurface which is older than restoration time τ such that projections of their directions of maximal change over the iso-value surfaces of the restored vertical coordinate t.sub.τ are approximately orthogonal to local co-axis vector b.sub.τ and the local axis vector a.sub.τ, respectively (see e.g., eqn. (42)).

(63) In some embodiments, the values of the restored horizontal coordinates u.sub.τ and v.sub.τ are constrained in parts of the subsurface which are older than the restoration time τ to be respectively equal on twin points on τ-active faults, wherein twin points are points on opposite sides of a τ-active fault that were collocated at the restoration time τ and are located on the same fault stria in the present day model, to merge the twin points into the same position in the restored model by sliding the twin points towards each other in a direction tangential to the surface of the τ-active fault (see e.g., eqn. (43)).

(64) In some embodiments, the values of the restored horizontal coordinates u.sub.τ and v.sub.τ are constrained in parts of the subsurface which are older than the restoration time τ to be respectively equal on mate points on τ-inactive faults, wherein mate points are points on opposite sides of a τ-inactive fault that are collocated at present day time, to move mate points together on opposite sides of τ-inactive faults (see e.g., eqn. (43)).

(65) In operation 1740, a processor may display an image of the restored 3D model of the subsurface geology of the Earth such that each point in the 3D model is positioned at the restored coordinates u.sub.τ, v.sub.τ, t.sub.τ defining the location that a piece of sediment represented by the point was located at the restoration time in the past τ.

(66) In some embodiments, the processor may receive an increasing chronological sequence of past restoration times τ.sub.1, τ.sub.2, . . . , τ.sub.n. For each restoration time τ. in sequence τ.sub.1, τ.sub.2, . . . , τ.sub.n) the processor may repeat operations 1720-1730 to compute a corresponding 3D restoration transformation Rτ.sub.i. 3D restoration transformation Rτ.sub.i restores the part of the subsurface older than horizon H.sub.τi to its predicted past geometry at time τ.sub.i, e.g., to 3D restored coordinates u.sub.τi, v.sub.τi, and t.sub.τi.

(67) In operation 1750, in some embodiments, a processor may play a moving image sequence in which the 3D model is iteratively restored in a forward or reverse order of the sequence of past restoration times τ.sub.1, τ.sub.2, . . . , τ.sub.n to visualize changes in the subsurface geology over the passage of time.

(68) In some embodiments, the processor may edit the model in the restoration space and then reverse the restoration transformation to apply those edits in the present day space. For example, the processor may edit the depositional values u, v, and t associated with the restored 3D model, and then reverse transform the restored 3D model forward in time from the predicted past geometry at the restoration time in the past τ to the present day measured geometry using an inverse of the 3D restoration transformation 200 to incorporate the edits from the restored model into the present day model.

(69) Decompaction at Intermediate Restoration Time τ:

(70) Compaction may refer to the pore space reduction in sediment within the Earth's subsurface. Compaction is typically caused by an increase in load weight of overlying geological layers as they are deposited over time. As sediment accumulates, compaction typically increases, as time and depth increase. Conversely, porosity typically decreases, as time and depth increase. For example, at a depositional time to when a layer is deposited with no overlaying geology, the depositional model has minimal or no compaction and maximum depositional porosity ψ.sub.0. At an intermediate restoration time τ, when there is an intermediate load of overlying deposited layers, the restored model has an intermediate level of compaction and an intermediate level of porosity ψ.sub.τ (or simply ψ). At the present-day time t.sub.p, when the present-day model has the most deposited layers, the present-day model typically has a maximal level of compaction and minimum porosity ψ.sub.p. Accordingly, the depositional porosity ψ.sub.0 is greater than the intermediate time porosity ψ.sub.τ, which in turn is greater than the present-day porosity ψ.sub.p, i.e., ψ.sub.0>ψ.sub.τ>ψ.sub.p. Further, because deeper layers are typically deposited at relatively earlier times than are shallower layers, within each model at the same time τ, a relatively deeper geological layer typically experiences a relatively greater load than does a relatively shallower geological layer, resulting in greater compaction and lesser porosity.

(71) Whereas compaction is a result of deposition over the forward passage of time, the process of restoration reverses the passage of time to visualize geology at an intermediate time in the past τ (before the present day and after the start of deposition of the oldest subsurface layer). Accordingly, embodiments of the invention generate a restoration model by reversing the effects of compaction in a process referred to as “decompaction” to more accurately depict how the geometry of geological layers change as their depths increase. Whereas compaction compresses the geological layers, decompaction reverses those effects, decompressing and uplifting terrains, resulting in increased layer thicknesses and increased intermediate time porosity ψ.sub.τ (or simply ψ) in the restored domain as compared with the compacted present-day domain ψ.sub.p. Decompaction decompresses the geology by a greater amount the earlier the intermediate restoration time τ is in the past and the deeper the layer is underneath the Earth's surface.

(72) Conventional decompaction techniques, however, are notoriously unreliable. Laboratory experiments on rock samples show that, during burial when sediments contained in a volume V(r.sub.τ) compact under their own weight, their porosity ψ(r.sub.τ) exponentially decreases according to Athy's law:
Ψ(r.sub.τ)≅Ψ.sub.o(r.sub.τ).Math.exp{−κ(r.sub.τ).Math.δ(r.sub.τ)} ∀r.sub.τ∈G.sub.T  (52)
where V(r.sub.τ) represents an infinitely small volume of sediment centered on a point r.sub.τ∈G.sub.T underneath the sea floor S.sub.τ(0)≡H.sub.τ, δ(r.sub.τ) is the absolute distance, or depth, from point r.sub.τ∈Gτ to sea floor S.sub.τ(0) measured at restoration time τ, and ψ.sub.0(r.sub.τ)<1 and κ(r.sub.τ) are known non-negative coefficients which depend only on rock type at location r.sub.τ. ψ.sub.0(r.sub.τ) is the porosity of the rock type with approximately no (zero) compaction, i.e., the porosity at its depositional time t.sub.0 before any layers were deposited to compress from above. κ is an experimental measurement derived from compression experiments of Athy's law performed in laboratory tests. As an example, assuming that geological depth δ(r.sub.τ) is expressed in meters, the following average coefficients for sedimentary terrains were observed in southern Morocco:

(73) TABLE-US-00001 Rock Type Ψ.sub.o κ Siltstone 0.62 0.57 × 10.sup.−3 Clay 0.71 0.77 × 10.sup.−3 Sandstone 0.35 0.60 × 10.sup.−3 Carbonates 0.46 0.23 × 10.sup.−3 Dolomites 0.21 0.61 × 10.sup.−3
Because, in the restored Gτdr-space, −t.sub.τ(r.sub.τ) measures the vertical distance from point r.sub.τ to the sea floor Sτ(0), the depth δ(r.sup.τ) in equation (52) may be equivalently expressed as:
δ(r.sub.τ)=−t.sub.τ(r.sub.τ) ∀r.sub.τ∈G.sub.T  (53)
Accordingly, in the context of embodiments of the invention, Athy's law may be reformulated as:

(74) $\begin{array}{cc}\stackrel{\_}{\Psi }\left({\stackrel{\_}{r}}_{\tau }\right)\simeq {\stackrel{\_}{\Psi }}_o\left({\stackrel{\_}{r}}_{\tau }\right).Math.\exp \left\{\stackrel{\_}{\kappa }\left({\stackrel{\_}{r}}_{\tau }\right).Math.t_{\tau }\left({\stackrel{\_}{r}}_{\tau }\right)\right\}\forall {\stackrel{\_}{r}}_{\tau }\in {\stackrel{\_}{G}}_{\tau }& \left(54\right)\end{array}$
Athy's law alone, however, incorrectly models porosity ψ and therefore often models decompaction inaccurately. Under Athy's law, restoration porosity ψ depends only on predictions extrapolated based on rock properties (ψ.sub.0 and κ), but does not actually measure real-world porosity. Because Athy's law is not rooted in the real-world geology, it often leads to inaccurate overestimated or underestimated compaction. Further, Athy's law models compaction based on porosity only at the time of earliest deposition, ψ.sub.0, but not porosity that occurs in the present-day, ψ.sub.9. Once the model is transformed from the present-day to restored time z, but prior to decompaction, the restored model still erroneously exhibits present-day compaction ψ.sub.p. Because Athy's law does not eliminate present-day compaction, which erroneously over-compresses terrains compared to restoration porosity, the resulting model is incorrectly decompacted at the restoration time τ.

(75) Embodiments of the invention improve decompaction techniques by modeling decompaction at an intermediate restoration time in the past τ based on real-world measurements of present-day compaction ψ.sub.p experimentally observed within the subsurface geology of the Earth. Modeling decompaction based on present-day compaction measurements accounts for the many real-world geological variables, such as those in the above example scenarios, that Athy's law misses.

(76) Some embodiments accurately decompact the restoration model by simultaneously (1) removing the impact of present-day compaction affecting terrains in (incorrectly) restored version at time τ (e.g., “total” decompaction, such as, defined in equations (58)); and (2) recompacting these terrains according to their depth in the restored model (e.g., “partial” recompaction, such as, defined in equations (59)). Embodiments of the invention solve the difficult problem of performing these two operations (decompaction and recompaction) simultaneously.

(77) Reference is made to FIG. 18, which schematically illustrates an example 3D geological volume of a compacted model 1810 representing the porosity of a subsurface region before decompaction (right image) and a corresponding decompacted model 1800 representing the porosity of the region after decompaction (left image) in the restored Gτ-space at an intermediate restoration time in the past τ, according to an embodiment of the invention. Embodiments of the invention replace original restoration coordinates of the compacted model 1810 with new restoration coordinates {u.sub.τ, v.sub.τ, t.sub.τ}.sub.r.sub.τ of the decompacted model 1800. Decompacted model 1800 may represent a new u.sub.τ, v.sub.τ, t.sub.τ-transform from the present-day model in Gτ-space to the restored Gτ-space that restores the terrains and induces thickness variations as a consequence of decompaction. This decompaction transformation is modeled to be the inverse (“reversing time”) of the compaction that occurred over the forward passage of time between geological-time τ and the present geological-time. Some embodiments may start with a region of the compacted model 1810 under the horizon Hτ (geology deposited before time τ with a present-day level of compaction) and restore the region to the decompacted model 1800 in 6T-space (geology deposited before time τ with a level of compaction at intermediate time τ). Because the compacted model 1810 has not yet been decompacted, its low porosity is similar to the present-day porosity, yielding vertical lengths dh(r.sub.τ) that are too short and compressed for the restored time τ. Accordingly, decompaction vertically stretches the lengths dh.sup.⊕(r.sub.τ) of the decompacted model 1800 to yield a greater porosity predicted at the time in the past τ. This process may repeat iteratively, layer-by-layer, starting at the top horizon Hτ deposited at the restoration time τ and ending at the bottom horizon deposited at the depositional time t.sub.0.

(78) Elasto-plastic mechanical frameworks developed to model compaction rely on a number of input parameters which may be difficult for a geologist or geomodeler to assess and are solved using a complex system of equations. Isostasic approaches are typically simpler to parameterize and still provide useful information on basin evolution. Therefore, compaction may be considered a primarily one-dimensional vertical compression induced by gravity which mainly occurs in the early stages of sediment burial when horizons are still roughly horizontal surfaces close to the sea floor.

(79) At any point r.sub.τ∈Gτ within a geological layer, the decompacted thickness dh.sup.⊕(r.sub.τ) e.g., of a vertical probe of infinitely small volume V(r.sub.τ) comprising an infinitely short column of sediment roughly orthogonal to the restored horizon passing through r.sub.τ is linked to the thickness dh(r.sub.τ) of the shorter, compacted vertical column by, for example, the following relationship:

(80) $\begin{array}{cc}\forall {\stackrel{\_}{r}}_{\tau }\in {\stackrel{\_}{G}}_{\tau }\text{:}.Math.\begin{array}{c}d{\stackrel{\_}{h}}^{\oplus }\left({\stackrel{\_}{r}}_{\tau }\right)=\frac{1}{1-{\stackrel{\_}{\varphi }}_{\tau }\left({\stackrel{\_}{r}}_{\tau }\right)}.Math.d\stackrel{\_}{h}\left({\stackrel{\_}{r}}_{\tau }\right)\\ \mathrm{with}\text{:}{\stackrel{\_}{\varphi }}_{\tau }\left({\stackrel{\_}{r}}_{\tau }\right)={\stackrel{\_}{\Psi }}_o\left({\stackrel{\_}{r}}_{\tau }\right)-{\stackrel{\_}{\Psi }}_o\left({\stackrel{\_}{r}}_{\tau }\right)\in \left[0,1\right]\end{array}& \left(55\right)\end{array}$
In this equation, ϕ(r.sub.τ) denotes the “compaction coefficient” which characterizes the vertical shortening of the probe at restored location r.sub.τ∈Gτ. As an example, FIG. 18 shows the same infinitely short vertical column of sediment where average porosity is equal to (ψ.sub.0=1/3) before compaction and (ψ=1/6) after compaction. The compaction coefficient (ψ.sub.0−ψ) is then equal to (ϕ=1/6) and column shortening (1−ϕ) is (5/6).

(81) Taking Present Day Compaction into Account to Decompact the Restored Model in Gτ-Space:

(82) Compacted model 1810, built assuming there is no compaction, incorrectly ignores the compaction characterized by present-day porosity ψ.sub.p(r.sub.τ). Compacted model 1810 thus results in geology with greater compaction and smaller porosity than occurred at intermediate restoration time τ. Embodiments of the invention correct the restored model by decompacting compacted model 1810. The decompaction process involves decompressing the vertical dimension's compacted height dh(r.sub.τ) or compacted time dt(r.sub.τ) (relatively shorter) to elongate the vertical dimension with a decompacted height dh.sup.⊕(r.sub.τ) or decompacted time dt.sup.⊕(r.sub.τ) (relatively longer) (see e.g., equation (60) and/or (64)). This decompaction of height (e.g., in equation (60)) or time (e.g., in equation (64)) is elongated based on compaction coefficient ϕ.sub.τ.sup.⊕(f.sub.t), which is a function of the present-day porosity ψ.sub.p(r.sub.τ) (see e.g., equation (56)). Because the present-day porosity ψ.sub.p(r.sub.τ) is less than the restoration porosity ψ.sub.τ(r.sub.τ), the ratio term in equations (60) and (64) is >1. Accordingly, the decompacted length dh.sup.⊕(r.sub.τ) and time dt.sup.⊕(r.sub.τ) are greater than the compacted length dh(r.sub.τ) and time dt(r.sub.τ), respectively, resulting in an elongation of the vertical dimension after decompaction. This elongation is thus defined based on real-world measurements of the present-day porosity ψ.sub.p(r.sub.τ), which yields more accurate decompaction than conventional simulations that ignore real-world porosity and compaction, such as Athy's law.

(83) Present-day porosity ψ.sub.p (r.sub.τ) may be measured by direct inspection of the Earth's subsurface material composition. In one example, porosity may be measured by directly analyzing core samples of the Earth's subsurface, for example, using a variety of methods to compare bulk rock volume and total sample volume. In one example, porosity may be derived from well logs, which are measurements performed on rock inside wells. Samples may be collected and porosity measured at regularly or irregularly spaced intervals within the Earth (e.g., bored into the Earth or along well paths). After porosity measurements are taken at those discrete locations, porosity may be extrapolated throughout the entire studied domain. In one example, at least one (and preferably multiple) samples are collected at each distinct depositional layer or depth (e.g., deposited at each distinct period of time).

(84) Example decompaction processes may proceed as follows:

(85) Let ϕ.sub.τ.sup.⊖(r.sub.τ) be a total compaction coefficient (representing a total compaction as a difference between the minimum present-day porosity and maximum depositional porosity) and let ϕ.sub.τ.sup.⊕(r.sub.τ) be an intermediate compaction coefficient (representing a partial compaction as a difference between the intermediate restoration porosity and maximum depositional porosity). The pair of compaction coefficients, ϕ.sub.τ.sup.⊕(r.sub.τ) and ϕ.sub.τ.sup.⊕(r.sub.τ), may be defined, for example, as:

(86) $\begin{array}{cc}\forall {\stackrel{\_}{r}}_{\tau }\in {\stackrel{\_}{G}}_{\tau }\text{:}.Math.\begin{array}{c}{\stackrel{\_}{\varphi }}_{\tau }^{\ominus }\left({\stackrel{\_}{r}}_{\tau }\right)={\stackrel{\_}{\Psi }}_o\left({\stackrel{\_}{r}}_{\tau }\right)-{\stackrel{\_}{\Psi }}_p\left({\stackrel{\_}{r}}_{\tau }\right)\\ {\stackrel{\_}{\varphi }}_{\tau }^{\oplus }\left({\stackrel{\_}{r}}_{\tau }\right)={\stackrel{\_}{\Psi }}_o\left({\stackrel{\_}{r}}_{\tau }\right)-\stackrel{\_}{\Psi }\left({\stackrel{\_}{r}}_{\tau }\right)\end{array}& \left(56\right)\end{array}$
Because compaction typically increases over time, the present day porosity ψ.sub.p(r.sub.τ) may be assumed to be less than the restored time porosity ψ(r.sub.τ):
ψ.sub.p(r.sub.τ)<ψ(r.sub.τ) ∀r.sub.τ∈Gτ  (57)
This inequality implies that intermediate compaction coefficient ϕ.sub.τ.sup.⊕(r.sub.τ)≤total compaction coefficient ϕ.sub.τ.sup.⊖(r.sub.τ), and so, the ratio in equations (60) and (64) is greater than 1, resulting in a vertical elongation in height to dh.sup.⊕(r.sub.τ) and/or time to dt.sup.⊕(f.sub.τ) in the decompacted model 1800 relative to the compacted model 1810.

(87) Considering once again the vertical probe introduced above in restored space Gτ, decompaction may proceed by using equation (55) twice, once in a forward and then in a backward transformation, for example, as follows: 1. First, to completely cancel out the compaction characterized by given, present day porosity ψ.sub.p(r.sub.τ), a “total” vertical decompaction is applied by updating dh(r.sub.τ) as follows:

(88) $\begin{array}{cc}d{\stackrel{\_}{h}}_o\left({\stackrel{\_}{r}}_{\tau }\right)=\frac{1}{1-{\stackrel{\_}{\varphi }}_{\tau }^{\ominus }\left({\stackrel{\_}{r}}_{\tau }\right)}.Math.d\stackrel{\_}{h}\left({\stackrel{\_}{r}}_{\tau }\right)& \left(58\right)\end{array}$ After canceling the total compaction in the first operation, the probe porosity is equal to the depositional porosity ψ.sub.0(r.sub.τ) having no or negligible compaction. 2. Next, a “partial” recompaction is applied as a function of the actual restoration porosity ψ(r.sub.τ) approximated by equation (54) at geological-time τ, for example, as:
dh.sup.⊕(r.sub.τ)={1−ϕ.sub.τ.sup.⊕(r.sub.τ)}.Math.dh.sub.o(r.sub.τ)  (59) After this second operation, the probe porosity is equal to the intermediate time restoration porosity ψ(r.sub.τ).

(89) Therefore, to take present-day compaction into account, equation (55) may be replaced, for example, by:

(90) $\begin{array}{cc}d{\stackrel{\_}{h}}^{\oplus }\left({\stackrel{\_}{r}}_{\tau }\right)=\frac{1-{\stackrel{\_}{\varphi }}_{\tau }^{\oplus }\left({\stackrel{\_}{r}}_{\tau }\right)}{1-{\stackrel{\_}{\varphi }}_{\tau }^{\ominus }\left({\stackrel{\_}{r}}_{\tau }\right)}.Math.d\stackrel{\_}{h}\left({\stackrel{\_}{r}}_{\tau }\right)\forall {\stackrel{\_}{r}}_{\tau }\in {\stackrel{\_}{G}}_{\tau }& \left(60\right)\end{array}$
where compaction coefficient ϕ.sub.τ.sup.⊖(r.sub.4) is based on the measured present-day porosity ψ.sub.p(r.sub.τ), for example, as defined in equation (56). Accordingly, the decompacted vertical thickness dh.sup.⊕(r.sub.τ) at intermediate restoration time τ is elongated based on real-world measurements of the present-day compaction ψ.sub.p (r.sub.t) experimentally observed within the subsurface of the Earth.

(91) Decompaction in GeoChron Based Restoration

(92) In the restored Gτ-space, the geological time of deposition t.sub.1(r.sub.τ) may be interpreted as an arc-length abscissa s(r.sub.τ) along the vertical straight line passing through r.sub.τ oriented in the same direction as the vertical unit frame vector {r.sub.t.sub.τ=r.sub.z}. Therefore, in the Gτ-space,
dt.sub.τ(r.sub.τ)=ds(r.sub.τ)=dh(r.sub.τ)  (61)
may represent the height of an infinitely short vertical column of restored sediment located at point r.sub.τ∈Gτ, subject to present-day compaction. As a consequence, to take compaction into account in the restored Gτ-space, according to equations (60) and (61), geological-time t.sub.τ(r.sub.τ) may be replaced by a “decompacted” geological-time t.sub.τ.sup.⊕(r.sub.τ) such that, for example:

(93) $\begin{array}{cc}\frac{{\mathrm{dt}}_{\tau }^{\ominus }}{{\mathrm{dt}}_{\tau }}{.Math.}_{{\stackrel{\_}{r}}_{\tau }}=\frac{d{\stackrel{\_}{h}}^{\oplus }\left({\stackrel{\_}{r}}_{\tau }\right)}{d\stackrel{\_}{h}\left({\stackrel{\_}{r}}_{\tau }\right)}=\frac{1-{\stackrel{\_}{\varphi }}_{\tau }^{\oplus }\left({\stackrel{\_}{r}}_{\tau }\right)}{1-{\stackrel{\_}{\varphi }}_{\tau }^{\ominus }\left({\stackrel{\_}{r}}_{\tau }\right)}& \left(62\right)\end{array}$

(94) Assuming that {r.sub.t.sub.τ=r.sub.z} is the unit vertical frame vector of the Gτ-space, it follows, for example, that:

(95) $\begin{array}{cc}\mathrm{grad}{t}_{\tau }^{\oplus }\left({\stackrel{\_}{r}}_{\tau }\right).Math.{\stackrel{\_}{r}}_{t_{\tau }}=\frac{{\mathrm{dt}}_{\tau }^{\ominus }\left({\stackrel{\_}{r}}_{\tau }+s.Math.{\stackrel{\_}{r}}_{t_{\tau }}\right)}{\mathrm{ds}}{.Math.}_{s=0}=\frac{{\mathrm{dt}}_{\tau }^{\ominus }}{{\mathrm{dt}}_{\tau }}{.Math.}_{{\stackrel{\_}{r}}_{\tau }}& \left(63\right)\end{array}$
From this, it can be concluded that the compacted geological-time t.sub.τ(r.sub.τ) of point r.sub.τ∈Gτ should be transformed into a decompacted geological-time t.sub.τ.sup.⊕(r.sub.τ), for example, honoring the following differential equation:

(96) $\begin{array}{cc}\mathrm{grad}{t}_{\tau }^{\oplus }\left({\stackrel{\_}{r}}_{\tau }\right).Math.{\stackrel{\_}{r}}_{t_{\tau }}=\frac{1-{\varphi }_{\tau }^{\oplus }\left({\stackrel{\_}{r}}_{\tau }\right)}{1-{\varphi }_{\tau }^{\ominus }\left({\stackrel{\_}{r}}_{\tau }\right)}\forall {\stackrel{\_}{r}}_{\tau }\in {\stackrel{\_}{G}}_{\tau }& \left(64\right)\\ \mathrm{with}\text{:}{\stackrel{\_}{\varphi }}_{\tau }^{\oplus }\left({\stackrel{\_}{r}}_{\tau }\right)={\stackrel{\_}{\Psi }}_o\left({\stackrel{\_}{r}}_{\tau }\right)-\stackrel{\_}{\Psi }\left({\stackrel{\_}{r}}_{\tau }\right)\&{\stackrel{\_}{\varphi }}_{\tau }^{\oplus }\left({\stackrel{\_}{r}}_{\tau }\right)={\stackrel{\_}{\Psi }}_o\left({\stackrel{\_}{r}}_{\tau }\right)-{\stackrel{\_}{\Psi }}_p\left({\stackrel{\_}{r}}_{\tau }\right)& \end{array}$

(97) Due to the vertical nature of compaction, on the top restored horizon {S(0)≡H.sub.τ}, geological-time t.sub.τ.sup.⊕(r.sub.τ) should vanish or reduce to zero and its gradient should be vertical. In other words, in addition to the constraint of equation (64), geological-time t.sub.τ.sup.⊕(r.sub.i) may also honor the following example boundary conditions where f.sub.u.sub.τ and r.sub.v.sub.τ may represent the unit horizontal frame vectors of the Gτ-space:

(98) $\begin{array}{cc}\forall {\stackrel{\_}{r}}_{\tau }^o\in \left\{{\stackrel{\_}{\mathrm{��}}}_{\tau }\left(0\right)\equiv {\stackrel{\_}{H}}_{\tau }\right\}\text{:}.Math.\begin{array}{cc}1)& t_{\tau }^{\oplus }\left({\stackrel{\_}{r}}_{\tau }^o\right)=0\\ 2)& \mathrm{grad}{t}_{\tau }^{\oplus }\left({\stackrel{\_}{r}}_{\tau }^o\right).Math.{\stackrel{\_}{r}}_{{\stackrel{\_}{u}}_{\tau }}=0\\ 3)& \mathrm{grad}{t}_{\tau }^{\oplus }\left({\stackrel{\_}{r}}_{\tau }^o\right).Math.{\stackrel{\_}{r}}_{{\stackrel{\_}{v}}_{\tau }}=0\end{array}& \left(65\right)\end{array}$
Boundary condition (65)(1) may ensure that the top restored horizon H.sub.τ is flat and planar at intermediate restoration time τ when it was deposited. Boundary conditions (65)(2) and (65)(3) may ensure that the direction of change (gradient) of the geological-time t.sub.τ(r.sub.τ) is vertical in the Gτ-space.

(99) As compaction is a continuous process, geological-time t.sub.τ.sup.⊕(r.sub.τ) may be continuous (e.g., C.sup.0-continuous) across all faults affecting Gτ. As a consequence, in addition to the constraints in equations (64) and (65), for any fault F in Gτ, geological-time t.sub.τ.sup.⊕(r.sub.τ) may also honor the following boundary conditions where (r.sub.F.sup.⊕, r.sub.F.sup.⊖) are pairs of “τ-mate-points” defined as collocated points respectively lying on the positive face F.sup.+ and negative face F.sup.− (opposite sides of fault F) at geological time τ:
t.sub.τ(r.sub.F.sup.⊕)=t.sub.τ(r.sub.F.sup.⊖)  (66)
∀F∈G.sub.τ& ∀(r.sub.F.sup.⊕,r.sub.F.sup.⊖).sub.τ∈F
Boundary condition (66) may ensure that, for any pair of collocated points on opposite sides of the fault, the two points have the same decompacted geological-time coordinate t.sub.τ.sup.⊕(r.sub.τ). This ensures there are no (or reduced) gaps or overlaps along the fault in the restored Gτ-space.

(100) Using an appropriate numerical method, t.sub.τ.sup.⊕(r.sup.τ) may be computed in Gτ whilst ensuring that differential equation (64) and boundary conditions (65) and (66) are honored. To ensure smoothness and uniqueness of t.sub.τ.sup.⊕(r.sub.τ), the following constraint may also be added:

(101) 0$\begin{array}{cc}\underset{\left(a,b\right)\in {\left\{u_{\tau },v_{\tau },t_{\tau }\right\}}^2}{.Math.}{\int }_{{\stackrel{\_}{G}}_{\tau }}^{}{\left\{{\partial }_a{\partial }_b{t}_{\tau }^{\oplus }\left({\stackrel{\_}{r}}_{\tau }\right)\right\}}^2.Math.d{\stackrel{\_}{r}}_{\tau }\mathrm{minimum}& \left(67\right)\end{array}$

(102) In summary, the following GeoChron Based Restoration technique may be used to take compaction into account: 1. Compute a numerical approximation of the elongated geological-time t.sub.τ.sup.⊕(r.sub.τ) in GT and use the reverse u.sub.τ, v.sub.τ, t.sub.i-transform to update t.sub.τ (r.sub.τ) in GT:
t.sub.τ(r.sub.τ)←t.sub.τ.sup.⊕(r.sub.τ) ∀r.sub.τ∈G.sub.τ;  (68) 2. Recompute numerical approximations of restoration functions u.sub.τ(r.sub.τ) and v.sub.τ(r.sub.τ) in Gτ to prevent voids and overlaps from being generated in the restored space, as, according to equations (22) and (23), u.sub.τ(r.sub.τ) and v.sub.τ(r.sub.τ) depend on t.sub.τ(r.sub.τ); 3. Build the “decompacted” restored space Gτ as the new, direct u.sub.τ, v.sub.t, t.sub.τ-transform of geological space GT observed today.

(103) This approach to decompaction may be seamlessly integrated into the GeoChron Based Restoration framework according to embodiments of the invention and is wholly dissimilar to the sequential decompaction following Athy's law along IPG-lines. In particular, embodiments of the invention perform decompaction based on real-world present-day porosity, a quantity that is accurately measured and extrapolated for any type of rock without having to make assumptions. Additionally, embodiments of the invention allow decompaction in the restored Gτ space representing the Earth's subsurface at an intermediate restoration time in the past τ, before the present day and after the start of deposition of the oldest subsurface layer being imaged.

(104) An Analytical Solution

(105) In the general case, the system of equations (64), (65) and (66) is typically too complex to be solved analytically and may be approximated using numerical methods. However, in a specific case where κ(r.sub.τ), ψ.sub.0(r.sub.τ) and ψ.sub.p(r.sub.τ) are all constant ∀r.sub.τ∈Gτ, the compaction ratio may be integrated at once over the entire domain, and there is no need to iteratively and independently decompact one layer at a time. This special case allows an analytical solution to the system of equations, for example, as follows.

(106) In this special case, in Gτ, terrain porosity is homogeneous and characterized as for example follows where κ, ψ.sub.0 and ψ.sub.p are known constants:

(107) $\begin{array}{cc}\forall {\stackrel{\_}{r}}_{\tau }\in {\stackrel{\_}{G}}_{\tau }\text{:}.Math.\begin{array}{c}\stackrel{\_}{\kappa }\left({\stackrel{\_}{r}}_{\tau }\right)=\kappa \\ {\stackrel{\_}{\Psi }}_o\left({\stackrel{\_}{r}}_{\tau }\right)={\Psi }_o\\ {\stackrel{\_}{\Psi }}_p\left({\stackrel{\_}{r}}_{\tau }\right)={\Psi }_p\end{array}& \left(69\right)\end{array}$

(108) Due to its homogeneity, Gτ may be considered continuous and the intrinsic, vertical nature of compaction implies that any function φ(r.sub.τ) defined in Gτ associated to compaction may only depend on the vertical component t.sub.τ(r.sub.τ) of r.sub.τ. Therefore, it follows, for example, that:
φ(r.sub.τ)=φ(t.sub.τ(r.sub.τ))=φ(t.sub.τ)  (70)

(109) Let constants A and B be defined, for example, by:

(110) $\begin{array}{cc}A=\frac{1-{\Psi }_o}{1-{\Psi }_o+{\Psi }_p};B=\frac{{\Psi }_o}{1-{\Psi }_o+{\Psi }_p}& \left(71\right)\end{array}$

(111) Let the following example functions be derived from Athy's law in equation (54) and equations (56) and (69):
ϕ.sup.⊖(t.sub.τ)=Ψ.sub.o−Ψ.sub.p; ϕ.sub.τ.sup.⊕(t.sub.τ)=ψ.sub.o−ψ.sub.o.Math.exp(κ.Math.t.sub.τ) ∀t.sub.τ≤0   (72)

(112) On the one hand, the following example indefinite integral holds true:

(113) $\begin{array}{cc}\int \frac{1-{\stackrel{\_}{\varphi }}_{\tau }^{\oplus }\left(s\right)}{1-{\stackrel{\_}{\varphi }}_{\tau }^{\ominus }\left(s\right)}.Math.\mathrm{ds}=\int \left\{A+B.Math.\exp \left(\kappa s\right)\right\}.Math.\mathrm{ds}=A.Math.s+\frac{B}{\kappa }.Math.\exp \left(\kappa s\right)& \left(73\right)\end{array}$

(114) On the other hand, according to equation (62):

(115) $\begin{array}{cc}{t}_{\tau }^{\oplus }\left(t_{\tau }\right)={\int }_0^{t_{\tau }}\frac{1-{\stackrel{\_}{\varphi }}_{\tau }^{\oplus }\left(s\right)}{1-{\stackrel{\_}{\varphi }}_{\tau }^{\ominus }\left(s\right)}.Math.\mathrm{ds}=A.Math.t_{\tau }+\frac{B}{\kappa }.Math.\left\{\exp \left(k.Math.t_{\tau }\right)-1\right\}& \left(74\right)\end{array}$

(116) Therefore, for any {t.sub.τ≤0}, the decompacted restoration function t.sub.τ.sup.⊕(t.sub.τ) may be analytically defined, for example, by:

(117) $\begin{array}{cc}{t}_{\tau }^{\oplus }\left(t_{\tau }\right)=A.Math.t_{\tau }+\frac{B}{\kappa }.Math.\left\{\exp \left(k.Math.t_{\tau }\right)-1\right\}& \left(75\right)\\ \mathrm{with}\text{:}A=\frac{1-{\Psi }_o}{1-{\Psi }_o+{\Psi }_p};B=\frac{{\Psi }_o}{1-{\Psi }_o+{\Psi }_p}& \end{array}$

(118) Other equations or permutations of these equations or terms may also be used.

(119) Reference is made to FIG. 19, which is a flowchart of a method for decompacting a 3D model of the subsurface geology of the Earth at an intermediate restoration time in the past τ, according to an embodiment of the invention.

(120) In operation 1910, a processor may receive a 3D model of present-day geometry of the subsurface geology and a measure of present-day porosity experimentally measured within the subsurface geology of the Earth. The present day model may be measured tomographically by scanning the Earth's subsurface e.g., as described in reference to FIGS. 14 and 15. To obtain the measure of present-day porosity, a probe may burrow into the Earth's subsurface or into one or more wells to collect and/or analyze material from within the subsurface geology of the Earth. Samples of subsurface materials are collected at spaced intervals, from which porosity is extrapolated throughout the studied domain.

(121) In operation 1920, a processor may select or receive a past restoration time τ that is “intermediate” or prior to the present time and after the start of the subsurface's deposition (the time period when an oldest horizon surface in the 3D model was originally deposited).

(122) In operation 1930, a processor may restore the 3D model from the present day measured geometry (e.g., present day model G.sub.τ 202 in xyz-space G 220) to the predicted past geometry at the restoration time in the past τ (e.g., restored model G.sub.γ 203 in u.sub.τ v.sub.τ t.sub.τ-space 219) using a 3D restoration transformation (e.g., u.sub.τ v.sub.τ t.sub.τ-transform 201). The 3D model may be restored, for example, as described in reference to FIG. 17. Prior to decompaction, the restored model may be a compacted model (e.g., 1810 of FIG. 18).

(123) In operation 1940, a processor may decompact the vertical dimension of the restored 3D model. This may expand, stretch or elongate compacted vertical lengths in the compacted model (e.g., 1810 of FIG. 18) to relatively longer vertical lengths in a decompacted model (e.g., 1800 of FIG. 18). In one embodiment, the vertical length may be a measure of height and the vertical dimension may be expanded from relatively shorter heights dh(r.sub.τ) in the compacted model to relatively longer heights dh.sup.⊕(r.sub.τ) in a decompacted model (e.g., as defined in equation (60)). In another embodiment, the vertical length may be a measure of geological-time when the particles of sediment were originally deposited on the Earth's surface and the vertical dimension may be expanded from relatively shorter times dt(r.sub.τ) in the compacted model to relatively longer times dt.sup.⊕(r.sub.t) in a decompacted model (e.g., as defined in equation (64)). The vertical lengths may be elongated based on a relationship between a depositional porosity (e.g., ψ.sub.0(r.sub.τ)) of the geological layers at the time sediment in those layers was deposited, restoration porosity (e.g., ψ(r.sub.τ)) of the geological layers at the restoration time in the past τ, and the present-day porosity (e.g., ψ.sub.p (r.sub.t)) of the geological layers experimentally measured in the present-day. In some embodiments, the relationship between the depositional porosity, the restoration porosity, and the present-day porosity may be, for example:

(124) $\frac{1-{\stackrel{\_}{\varphi }}_{\tau }^{\oplus }}{1-{\stackrel{\_}{\varphi }}_{\tau }^{\ominus }},$
where compaction coefficients ϕ.sub.τ.sup.⊕=ψ.sub.0−ψ and ϕ.sub.τ.sup.⊖=ψ.sub.0−ψ.sub.p for all points in the restored 3D model, e.g., as defined in equations (56), (60) and (64). Since porosity decreases over time, the restoration porosity is typically greater than the present-day porosity (e.g., equation (57)) and typically less than the depositional porosity. Accordingly, the compaction coefficients have a relationship ϕ.sub.τ.sup.⊕≤ϕ.sub.τ.sup.⊖, and the relationship between the depositional, restoration, and present-day porosities, e.g.,

(125) $\frac{1-{\overline{\varphi }}_{\tau }^{\oplus }}{1-{\stackrel{\_}{\varphi }}_{\tau }^{\ominus }},$
is greater than 1, resulting in a stretching or elongating effect to increase the vertical lengths when they are decompacted.

(126) In some embodiments, a processor may decompact the vertical dimension of the restored 3D model by a combination (e.g., equation (60)) of total decompaction corresponding to an increase in porosity from the present day porosity to the depositional porosity (e.g., equation (58)) and partial recompaction corresponding to a partial decrease in the porosity from the depositional porosity to the restored porosity (e.g., equation (59)).

(127) At the restored time in the past τ, the geological layers above H.sub.τ (e.g., H.sub.τ+1 . . . H.sub.n) did not yet exist, so decompaction may elongate lengths of geological layers below the horizon layer H.sub.τ (e.g., H.sub.0 . . . H.sub.τ) deposited at the restoration time in the past τ. In some embodiments, decompaction may be performed by iteratively decompacting the subsurface layer-by-layer, starting at the top horizon H.sub.τ deposited at the restoration time τ and ending at the bottom horizon H.sub.0 deposited at the depositional time. In some embodiments, the depositional porosity and the present-day porosity may be independently determined for each geological layer of the subsurface. In other embodiments, when the depositional porosity and the present-day porosity are substantially constant throughout the subsurface geology, decompaction may occur in one operation over the entire domain of the restored 3D model (e.g., as in equation (76)).

(128) Some embodiments may implement a boundary condition that ensures that a top horizon H.sub.τ deposited at the restoration time τ is a horizontal plane in the restored 3D model (e.g., equation (65)(1)). Additionally or alternatively, some embodiments may implement a boundary condition that ensures that a direction of change of geological-time when the particles of sediment were originally deposited on the Earth's surface is vertical in the restored 3D model (e.g., equation (65)(2) and (65)(3)). Additionally or alternatively, some embodiments may implement a boundary condition that ensures that, for any pair of collocated points on opposite sides of a fault, the two collocated points are decompacted to have the same coordinate (e.g., equation (66)).

(129) In some embodiments, for example, implemented in a past-time model, such as the GeoChron model, a processor may decompact the vertical dimension of the restored 3D model by: computing an elongated geological-time (e.g., dt.sup.⊕(r.sub.τ)) in the restored 3D model (e.g., by solving equation (64)), transforming the elongated geological-time from the restored 3D model to generate a 1D geological-time (e.g., t.sub.τ(r.sub.τ)) in the present-day 3D model (e.g., equation (68)), computing 2D paleo-depositional coordinates (e.g., u.sub.τ(r.sub.τ) and v.sub.τ(r.sub.τ)) based on the transformed geological-time (e.g., t.sub.τ(r.sub.τ)) in the present-day 3D model, and performing a 3D transformation (e.g., a u.sub.τ,v.sub.τ,t.sub.τ-transformation)) comprising the 1D geological-time and 2D paleo-depositional coordinates from the present-day 3D model (e.g., G.sub.τ) to the restored 3D model (e.g., G.sub.τ) that is decompacted based on the elongated geological-time.

(130) Operations of FIGS. 16, 17 and 19 may be performed for example using system 1505 of FIG. 15, e.g., by one or more processor(s) 140 of FIG. 15, or another suitable computing system. The embodiments disclosed in reference to FIGS. 16, 17 and 19 may be performed using other operations or orders of the operations, and the exact set of steps shown in the figures may be varied.

(131) In the past 30 years, many methods have been proposed to build geological models of sedimentary terrains having layers that are both folded and faulted. For any given geological-time τ, checking geological model consistency is considered both simpler and more accurate if terrains have previously been “restored” to their pre-deformational, unfolded and unfaulted state, as they were at geological-time τ.

(132) Embodiments of the invention provide a new, purely geometrical 3D restoration method based on the input of a depositional (e.g., GeoChron model). Embodiments of the invention are able to handle depositional models of any degree of geometrical and topological complexity, with both small and large deformations, do not assume elastic mechanical behavior, and do not require any prior knowledge of geo-mechanical properties. Embodiments of the invention further reduce or eliminate gaps and overlaps along faults as part of the restoration transformation and do not resort to any post-processing to minimize such gaps and overlaps. Compared to other conventional methods, embodiments of the invention minimize deformations and volume variations induced by geological restoration with a higher degree of precision, unequaled so far (see e.g., FIG. 5 and FIG. 9). Embodiments of the invention further ensure that 2D deformations of horizon surfaces induced by the uvt-transform are kept coherent with 3D deformations of volumes induced by the new proposed 3D restoration method.

(133) Referring to FIG. 3, for a given restoration time τ, the set of faults is split into τ-active and τ-inactive subsets. Such a distinction allows: deformations along faults 105 to be minimized, restoration to work correctly even though there are regions of G.sub.τ not continuously connected to H.sub.τ, gaps and overlaps along faults and the geometry of fault striae 600 are minimized by the restoration transformation, so no post-processing is needed to correct gaps or overlaps.

(134) Embodiments of the invention input a 3D model of sedimentary terrains in the subsurface. In one example, the input model may be the GeoChron™ model generated by SKUA® software for use in mining and oil and gas industries. Embodiments of the invention may build a 3D restoration transformation of this model in such a way that, after transformation, the new model represents terrains as they were at a given intermediate restoration-time τ (where τ.sub.1<τ<τ.sub.2, before the present day τ.sub.2 and after the time of the deposition of the oldest layer τ.sub.1).

(135) For example, G may represent the present day 3D geological domain of the region of the subsurface being modeled and G.sub.τ 202 may represent the subset of G containing particles of sediment that were deposited at a time prior to or equal to τ. In some embodiments, for all points r∈G, a geologic restoration transformation may move a particle of sediment observed today at location r to a new restored location f.sub.τ(r), e.g., defined as follows:
r.sub.τ(r)=r+R.sub.τ(r) ∀r∈G.sub.τ  (1)
where R.sub.τ(r) represents a 3D field of restoration vectors, e.g., generated to minimize deformations in G.sub.τ.

Depositional Model

(136) A depositional model may be generated by inputting a tomographic model of the present day subsurface geology of the Earth and transforming that geology to a past depositional time as each particle was configured when originally deposited in the Earth. Sedimentary particles are deposited over time in layers from deepest to shallowest from the earliest to the most recent geological time periods. Since various layers of terrain are deposited at different geological times, a depositional model does not image the geology at any one particular time period, but across many times periods, each layer modeled at the geological time when the layer was deposited. Accordingly, the vertical axis or depth in the depositional model may be a time dimension representing the time period of deposition, progressing from oldest to newest geological time as the model progresses vertically from deepest to shallowest layers.

(137) In one embodiment, the depositional model may be the GeoChron™ model, which is generated by SKUA™ software, that is routinely used by many oil & gas companies to build models of geologic reservoirs which help optimize hydrocarbon production and exploration. An example implementation of the GeoChron model is discussed in U.S. Pat. No. 8,600,708, which is incorporated by reference herein in its entirety. The depositional model is described in reference to the GeoChron model only for example, though any other depositional model may be used.

(138) Reference is made to FIG. 7, which schematically illustrates an example transformation from a present day model (upper-left image) to a depositional GeoChron model (bottom-right image), according to an embodiment of the invention. The transformation may be referred to as a “uvt-transform” 700 that transforms a particle of sediment observed today at location r=r(x,y,z) in the present day geological domain G (also referred to as “G-space”) 220 to be moved to a new depositional location r(r)=r(u,v,t) in the depositional geological domain G (also referred to as “G-space”). The new depositional location f has a vertical coordinate that is the geological time t(r) when the particle at location r was deposited and has horizontal or paleo-geographic coordinates {u(r), v(r)} equal to the lateral spatial location where the particle at r was located at its depositional time t(r). The paleo-geographic coordinates {u(r), v(r)} may be linked to the vertical time coordinate t(r) by different relationships (e.g., constrained by different systems of differential equations) depending on the structural style of their deposition (e.g., minimal deformation or flexural slip).

(139) In the example uvt-transform 700 shown in FIG. 7, when the geological time coordinate t(r) is equal to the curvilinear distance to the top horizon H.sub.τ 210 along curvilinear axis 222, the uvt-transform is a valid technique for imaging the depositional model. In other words, the uvt-transform is a valid depositional rendering technique if the module of its gradient grad t(r) honors the following constraint:
∥grad t(r)∥=1 ∀r∈G  (3)

(140) Embodiments of the invention observe that when ∥grad t(r)∥ differs from “1,” replacing the depositional coordinates {u(r), v(r), t(r)} of the uvt-transform 700 by new restoration coordinates {u.sub.τ(r), v.sub.τ(r), t.sub.τ(r)} where ∥grad t.sub.τ∥=1 allows the uvt-transform to be replaced by a u.sub.τ v.sub.τ t.sub.τ-transform that generates a valid restoration model at restoration time.sub.τ.

(141) In some embodiments, the depositional (e.g., GeoChron) model includes the following data structures stored in a memory (e.g., memory 150 of FIG. 15) (see FIGS. 1, 3, 6, and 7): A network of geological faults 105 within the present day domain G 220. A 3D corner-point grid F 100 that fills the G-space 220 with 3D polyhedral cells 108 (e.g., tetrahedra or hexahedra), without any gaps or overlaps in the studied domain, in such a way that no cell edge crosses any fault. The location of each node α 107 of grid Γ 100 in the G-space is denoted r(α). For each geological fault F 105, two disconnected and independently meshed, collocated surfaces F.sup.+ 103 and F.sup.− 104 on opposite sides of the fault 105. Surfaces F.sup.+ 103 and F.sup.− 104 may be composed of 2D facets from the 3D polyhedral cells of grid Γ 100 bordering F 105 on either side of the fault 105. Fault surfaces F.sup.+ 103 and F.sup.− 104 that are collocated in the present day model may, during the restoration process of transforming the model backwards in time, typically slide against one another, without generating gaps or overlaps between adjacent fault surfaces or fault blocks. Referring to FIGS. 1, 3 and 6, for each fault F 105, a set of pairs of points (r.sub.F.sup.+, r.sub.F.sup.−) (101,102) called “twin-points,” such that: 1. The two twin points in each pair are located on opposite sides of a corresponding pair of twin fault surfaces, r.sub.F.sup.+∈F.sup.+ and r.sub.F.sup.−∈F.sup.−. 2. At geological times before fault F 105 formed in the subsurface, particles of sediment were collocated which are observed today at locations r.sub.F.sup.+ and r.sub.F.sup.−. During the activation of fault F 105, particles of sediment initially located on F are assumed to slide along fault-striae (e.g., see FIG. 12), which are the shortest paths, on F, between pairs of twin points (r.sub.F.sup.+,r.sub.F.sup.−) (101,102). A tectonic style which may be either a “minimal deformation” style or a “flexural slip” style. Choosing this tectonic style is a model decision assumed to have been made by a structural geologist. A triplet (e.g., {u(r), v(r), t(r)}) of discrete coordinates defined on a 3D grid Γ 100 of the depositional G-space, such that, for a particle of sediment observed today at location r, the coordinate values {u(r), v(r)} represent the paleo-geographic coordinates of the particle at geological-time t(r) during the time period when it was deposited. According to the depositional (e.g., GeoChron) model, the paleo-geographic coordinates {u(r), v(r)} may honor different differential equations depending on the tectonic style.

(142) Moreover, referring to FIG. 1 and FIG. 2, the depositional model may have the following properties: Within the present day domain G, each geological horizon H.sub.τ 210 may be defined by a set of particles of sediment which were deposited at geological time r:
r∈H.sub.τ⇐.Math.t(r)=τ  (4) In other words, each horizon H.sub.τ 210 is a level-set (constant value) surface of the geological-time t. Paleo-geographic coordinates {u(r), v(r)} and twin-points (101,102) given as input are linked e.g. by the following equations:

(143) $\hspace{1em}.Math.\begin{array}{c}\{\left(r_F^{+},r_F^{-}\right)\text{is a pair of twin-points}}\\ \iff \{\begin{array}{c}\begin{array}{c}\begin{array}{c}{r}_F^{+}\in F^{+}\&r_F^{-}\in F^{-}\\ u\left(r_F^{-}\right)=u\left(r_F^{+}\right)\end{array}\\ v\left(r_F^{-}\right)=v\left(r_F^{+}\right)\end{array}\\ t\left(r_F^{-}\right)=t\left(r_F^{+}\right)\end{array}\begin{array}{c}\left(5\right)\\ \left(6\right)\\ \left(7\right)\\ \left(8\right)\end{array}\end{array}$ Additionally or alternatively, each pair of twin-points (r.sub.F.sup.+, r.sub.F.sup.−) (101,102) may be the intersection of a level set 210 of vertical depositional coordinate t(r) with a “fault stria” σ(r.sub.F.sup.−) 600 comprising a curved surface passing through point r.sub.F.sup.− 102 whose geometry is defined by geological rules, e.g., defining fault blocks sliding against one another according to tectonic forces and geological context. As a consequence of constraints defined by equations (6), (7), and (8) above, fault-striae (e.g., see FIG. 12) may characterize the paleo-geographic coordinates {u(r), v(r)}, and vice versa. Each point r∈G 214 may be characterized by its present day coordinates (e.g., {x(r), y(r), z(r)}) with respect to a present day coordinate system {r.sub.x, r.sub.y, r.sub.z} 220 comprising three mutually orthogonal unit vectors, e.g., where r.sub.z is assumed to be oriented upward.

(144) It would be appreciated by a person of ordinary skill in the art that the GeoChron model and its features described herein are discussed only as an example of a depositional model and that these elements may differ in other models or implementations without changing the essence of the invention.

u.SUP.τ v.SUB.τ .t.SUB.τ.-Transformation

(145) Referring to the volume deformation of FIG. 2, the restoration timer may be a given geological time in the past and subdomain G.sub.τ 202 may be a part of a 3D present day geological domain G that has terrains older than (deposited at a time prior to) or equal to restoration time τ and defined by a depositional model. Embodiments of the invention provide a new unfolding technique that replaces the uvt-transform (converting the present day model to a depositional model rendering all layers at their many respective times of deposition) by a u.sub.τ v.sub.τ t.sub.τ-transform 201 (converting the present day model to a restored model at a single restoration time τ before the present day but after the earliest times of deposition of the deepest model layer):

(146) $\begin{array}{cc}\left(x,y,z\right)\stackrel{u_{\tau }{v}_{\tau }{t}_{\tau }}{.fwdarw.}\left\{u_{\tau }\left(x,y,z\right),v_{\tau }\left(x,y,z\right),t_{\tau }\left(x,y,z\right)\right\}& \left(9\right)\end{array}$
Accordingly, present day geological space G.sub.τ 202 is transformed into a restored geological space G.sub.τ 203, such that: t.sub.τ(r) is a vertical spatial coordinate of the subsurface at the past restoration time τ, and is derived from, but different than, the geological time of deposition t(r). The vertical restoration coordinate t.sub.τ(r) honors the following constraint:
∥grad t.sub.τ(r)∥=1 ∀r∈G.sub.τ  (10) {u.sub.τ(r), v.sub.τ(r)} are lateral restoration coordinates derived from, but different than, the paleogeographic coordinates {u(r), v(r)} of the depositional model. restoration coordinates {u.sub.τ(r), v.sub.τ(r), t.sub.τ(r)} honor specific inventive constraints described below in such a way that, using the u.sub.τ v.sub.τ t.sub.τ-transform as a restoration operator minimizes deformations in the present day domain G.sub.τ. for each point r∈G.sub.τ 202, the restoration vector field R.sub.τ may be defined e.g. by:
R.sub.τ(r)=u.sub.τ(r).Math.r.sub.x+v.sub.τ(r).Math.r.sub.y+t.sub.τ(r).Math.r.sub.z−r  (11)

Volume Deformation

(147) Compaction may be handled in pre and post-restoration stages, as is known in the art. Thus, the model may be restored without taking compaction into account.

(148) Some embodiments of the invention provide an inventive volume deformation with a new set of inventive geometric constraints on the depositional model to allow geologic layers to be restored at a given geological time τ with a precision that has never before been reached. As shown in FIG. 2, in this volume deformation, paleo-geographic coordinates [u(r) and v(r)] and the geological time coordinate t(r) are replaced by new restoration coordinates respectively denoted {u.sub.τ(r), v.sub.τ(r)} and t.sub.τ(r).

(149) As shown in FIG. 2: a restored volume 203, denoted G.sub.τ, contains a (e.g., direct or “right-handed”) coordinate space 219 having orthogonal coordinate unit vectors {r.sub.u.sub.τ, r.sub.v.sub.τ, r.sub.t.sub.τ} and a family of horizontal planes {S.sub.τ(d):d≥0} 207 parallel to horizontal coordinate vectors {r.sub.u.sub.τ, r.sub.v.sub.τ}; a deformed version G.sub.τ of G.sub.τ contains a (e.g., direct or “right-handed”) coordinate space 220 having orthogonal coordinate unit vectors {r.sub.x, r.sub.y, r.sub.z} and a family of curved surfaces {S.sub.τ(d):d≥0} 208 parallel to horizon {H.sub.τ≡S.sub.τ(0)} (210=208).

(150) For simplicity and without loss of generality, the coordinate frame unit vectors {r.sub.u.sub.τ, r.sub.v.sub.τ, r.sub.t.sub.τ} 219 of the G.sub.τ-space and its origin Ō.sub.u.sub.τ .sub.v.sub.τ .sub.t.sub.τ may be equal to the coordinate frame unit vectors {r.sub.x, r.sub.y, r.sub.z} 220 of the G-space and its origin O.sub.xyz:
r.sub.u.sub.τ≡r.sub.x;r.sub.v.sub.τ≡r.sub.y;r.sub.t.sub.τ≡r.sub.z
Ō.sub.u.sub.τ .sub.v.sub.τ .sub.t.sub.γ≡O.sub.xyz  (12)

(151) Referring to FIG. 2, the following notation is used: each point r∈G.sub.τ 214 is transformed into point r.sub.τ∈G.sub.τ 213 and vice versa:
r.sub.τ∈G.sub.τ←.fwdarw.r∈G.sub.τ  (13) [u.sub.τ(f.sub.τ), ii.sub.τ(f.sub.τ)) represent the horizontal restoration coordinates of f.sub.τ 213 with respect to {r.sub.u.sub.τ, r.sub.v.sub.τ} 219 whilst t.sub.τ(r) represents the vertical restoration coordinate or altitude 204 of r.sub.τ 213 with respect to r.sub.t, oriented upward; for each point r E G.sub.τ 214: {x(r), y(r)} represent the horizontal present day coordinates of r with respect to {r.sub.x, r.sub.y} whilst z(r) represents the vertical present day coordinate or altitude of r with respect to the vertical unit frame vector r.sub.z oriented upward; and {u.sub.τ(r), v.sub.τ(r), t.sub.τ(r)} represent the restoration coordinates of r.sub.τ with respect to the restoration coordinate frame unit vectors {r.sub.u.sub.τ, r.sub.v.sub.τ, r.sub.t.sub.τ} 219 of the restored volume G.sub.τ 203.

(152) Equivalently to equations (12) and in accordance with equation (1), during restoration of G.sub.τ, a particle of sediment observed today at location r 214 is moved to a new location r(r) 213 defined e.g., as follows, where R.sub.τ(r) is a restoration vector field:
r(r)=r+R.sub.τ(r)  (14)
with, in matrix notation:

(153) 0$\begin{array}{cc}{R}_{\tau }\left(r\right)=\left[r_x,r_y,r_z\right].Math.\left[\begin{array}{c}{u}_{\tau }\left(r\right)-x\left(r\right)\\ v_{\tau }\left(r\right)-y\left(r\right)\\ t_{\tau }\left(r\right)-z\left(r\right)\end{array}\right]& \left(15\right)\end{array}$

(154) Referring to FIG. 2, surface S.sub.τ(0) 209 is located at an altitude (t.sub.τ=z.sub.τ.sup.O) with respect to the vertical unit vector r.sub.t.sub.τ oriented upward. In the frame of the presentation of the volume deformation and without loss of generality, z.sub.τ.sup.O may be assumed to be constant, e.g., equal to zero.

(155) Referring to FIG. 2, surface S.sub.τ(d) 207 is located at a distance (d) from S.sub.τ(0) 209, implying that:
t.sub.τ(r.sub.τ)=d−z.sub.τ.sup.O ∀r.sub.τ∈S.sub.τ(d)  (16)
such that:

(156) $.Math.\begin{array}{cc}d<0⟺{\stackrel{\_}{S}}_{\tau }\left(d\right)\mathrm{is}\mathrm{located}\mathrm{above}{\stackrel{\_}{S}}_{\tau }\left(0\right)d>0⟺{\stackrel{\_}{S}}_{\tau }\left(d\right)\mathrm{is}\mathrm{located}\mathrm{below}{\stackrel{\_}{S}}_{\tau }\left(0\right)& \left(17\right)\end{array}$

(157) FIG. 2 shows the folded present-day volume G.sub.τ 202 resulting from the deformation of restored volume G.sub.τ 203 under tectonic forces following either a “minimal deformation” or a “flexural slip” tectonic style:
G.sub.τ.fwdarw.tectonic forces.fwdarw.G.sub.τ.Math.G  (18)

(158) Referring to FIG. 2, the following notation is used: each horizontal plane S.sub.τ(d) 207 is transformed into a curved surface S.sub.τ(d) 208 “parallel” (e.g., this notion of “parallelism” may be characterized by equation (10)) to surface H.sub.τ 210=208 and each surface S.sub.τ(d) 208 is a level set of vertical restoration coordinate t.sub.τ(r); the images in G.sub.τ 202 of the (u.sub.τ), (v.sub.τ) 205, 206 and (t.sub.τ) 204 coordinate axes initially rectilinear and contained in volume G.sub.τ 203 now consist of curvilinear coordinate axes (223, 224) and 222.

(159) As shown in FIG. 2, the part of the subsurface observed today stratigraphically below H.sub.τ 210 may be identified with the deformed volume G.sub.τ 202, e.g., such that: S.sub.τ(0) is assumed to be identical to the horizon H.sub.τ 210 to be restored:
S.sub.τ(0)≡H.sub.τ  (19) which is equivalent to defining that, on horizon H.sub.τ 216, restored vertical coordinate t.sub.τ(r) is equal to z.sub.τ.sup.O; for any t<τ, the actual geologic horizon H.sub.t 216 is included (212) into the deformed volume G.sub.τ 202; note that, contrary to surfaces {S.sub.τ(d):d≥0} 208, horizons {H.sub.t:t <τ} 216 may be non-parallel to {H.sub.τ≡S.sub.τ(0)} 210=208; after restoration of the volume G.sub.τ 202 to its initial, unfolded state G.sub.τ 203: all horizons, faults and geological objects included in G.sub.τ 202 are dragged up by the embedding volume deformation, H.sub.τ≡S.sub.τ(0)} 209=207 may be defined as the restored sea floor as it was at geological time τ.

Minimizing Deformations and Volume Variation

(160) With compaction handled separately in pre and post restoration steps, leaving aside the very particular case of clay and salt layers, tectonic forces generally induce no or negligible variations in volume. Therefore, restoration coordinates {u.sub.τ(r), v.sub.τ(r), t.sub.τ(r)} may be chosen in such a way that the u.sub.τ v.sub.τ t.sub.τ-transform 201 of the present-day volume G.sub.τ 202 into the restored volume G.sub.τ 203 minimizes deformations and volume variations. This is achieved by constraining restoration coordinates {u.sub.τ(r), v.sub.τ(r), t.sub.τ(r)} to honor the two following conditions in the present day G.sub.τ domain: Surfaces {S.sub.τ(d):d≥0} 208 are level sets of the vertical restoration coordinate t.sub.τ (r) and, for any infinitely small increment ε, the thickness of the thin slice of the volume bounded by S.sub.τ(d) and S.sub.τ(d+ε) are, as much as possible, constant and equal to ε. In other words, S.sub.τ(d) and S.sub.τ(d+ε) are as parallel as possible. This is equivalent to honoring equation (10) as precisely as possible. In the frame of this invention, the consistency between the depositional (e.g., GeoChron) model provided as input and its restored version at geological time τ is preserved. Such a consistency is preserved if, and only if, the uvt-transform and the u.sub.τ v.sub.τ t.sub.τ-transform of H.sub.τ are identical. This is achieved by honoring the following inventive boundary conditions, referred to as the (u.sub.τ, v.sub.τ) boundary constraints:

(161) $\begin{array}{cc}\begin{array}{cc}\forall r_{\tau }^o\in H_{\tau }|\begin{array}{cc}1)& u_{\tau }\left(r_{\tau }^o\right)=u\left(r_{\tau }^o\right)\\ 2)& v_{\tau }\left(r_{\tau }^o\right)=v\left(r_{\tau }^o\right)\end{array}& \left(20\right)\\ \forall r_{\tau }^o\in H_{\tau }|\begin{array}{cc}1)& \mathrm{grad}{u}_{\tau }\left(r_{\tau }^o\right)\simeq \mathrm{grad}u\left(r_{\tau }^o\right)\\ 2)& \mathrm{grad}{v}_{\tau }\left(r_{\tau }^o\right)\simeq \mathrm{grad}v\left(r_{\tau }^o\right)\end{array}& \left(21\right)\end{array}& \end{array}$ Whilst taking the same given tectonic style into account (minimal deformation or flexural slip) as the one honored by paleo-geographic coordinates {u(r), v(r)}, lateral restoration coordinates {u.sub.τ(r), v.sub.τ(r)} may be defined so that their associated restoration deformations are minimized. To preserve consistency with boundary conditions (20) and (21), this is achieved by honoring the following inventive constraints: in a minimal deformation tectonic style context, the following “minimal deformation constraints” may be honored by coordinates {u.sub.τ, v.sub.τ}.sub.τ where t.sub.τ(r) is given, e.g., as follows:

(162) $\begin{array}{cc}\forall r\in G_{\tau }|\begin{array}{cc}1)& {\left\{\mathrm{grad}{u}_{\tau }.Math.\mathrm{grad}{v}_{\tau }\right\}}_r\simeq 0\\ 2)& {\left\{\mathrm{grad}{t}_{\tau }.Math.\mathrm{grad}{u}_{\tau }\right\}}_r\simeq 0\\ 3)& {\left\{\mathrm{grad}{t}_{\tau }.Math.\mathrm{grad}{v}_{\tau }\right\}}_r\simeq 0\end{array}& \left(22\right)\end{array}$ in a flexural slip tectonic style context, the following “flexural slip constraint” is coupled (containing both lateral restoration coordinates u.sub.τ and v.sub.τ) and may be honored by coordinates {u.sub.τ, v.sub.τ}.sub.r, e.g., as follows:
∀rϵG.sub.r:{grad.sub.su.sub.τ.Math.grad.sub.sv.sub.τ}.sub.r≅0  (23) where subscript “S” refers to a projection of the directions of maximal change over iso-value surfaces of the restored vertical coordinate tr. So as not to conflict with equations (20) and (21), and contrary to conventional depositional coordinates u and v (e.g., in the GeoChron model), new constraints (22) and (23) do not constrain ∥grad u.sub.τ∥, ∥grad v.sub.τ∥, ∥grad.sub.S u.sub.τ∥, or ∥grad.sub.S v.sub.τ∥ to be equal to “1”.

Restoration

(163) Referring to FIG. 1 and FIG. 2, at geological-time τ, the horizon H.sub.τ 210 to be restored was coincident with a given surface S.sub.τ(0) (208=210) considered as the sea-floor. The task of restoration includes: restoring horizon H.sub.τ 210 to its initial, unfaulted and unfolded state (e.g., mapping horizon H.sub.τ onto the sea floor S.sub.τ(0)) 209 and shifting all sedimentary terrains in such a way that, for each point r∈G: the particle of sediment currently located at point r moves to its former, “restored” location, where the particle was located at geological time τ, no overlaps or voids/gaps are created in the subsurface.

(164) At geological time τ, the sea floor S.sub.τ(0) (209) is assumed to be a contiguous, unfaulted surface whose altitude z.sub.τ.sup.O is a given function z.sub.τ.sup.O(u, v). In practice, S.sub.τ(0) (209=207) is typically a flat, horizontal plane whose altitude z.sub.τ.sup.O(u, v) at geological time τ is constant. Accordingly, for concision, z.sub.τ.sup.O may refer to a given function z.sub.τ.sup.O(u(r), v(r)) which may or may not be constant:
∀r∈H.sub.τ:z.sub.τ.sup.o stands for z.sub.τ.sup.o(u(r), v(r))  (24)

Compaction

(165) Deformation of sedimentary terrains is typically induced both by tectonic forces and terrain compaction. In order to model separately the effects of these phenomena, the restoration process may proceed as follows: First, in a preprocessing phase, a total decompaction may be applied to the terrains to cancel the impact of compaction as it is observed today, at the present day or current geological time; Next, the effects of compaction being canceled, a depositional-based restoration process taking only tectonic deformations into account (and not compaction) is applied to restore the geometry of the subsurface as it would have been observed at geological time τ; Finally, in a post-processing phase, a partial recompaction is applied to the restored terrains in order to take compaction into account, as it could have been observed at geological time τ.

Depositional Based Restoration

(166) As an input to the restoration process, a given depositional (e.g., GeoChron) model may be received from storage in a digital device (e.g., from memory 150 of FIG. 15).

(167) Referring to FIG. 2, a geological time τ may be selected that is associated with the given horizon H.sub.τ 210 to be restored and the given altitude z.sub.τ.sup.O of the surface S.sub.τ(0) 209 onto which the horizon Hτ should be restored.

(168) The region G.sub.τ 202 may be retrieved as the part of the depositional model where geological time of deposition t(r) is less than or equal to τ (subsurface regions deposited in a layer deeper than or equal to the layer deposited at time τ).

(169) The set of faults may be split into a subset of τ-active faults cutting H.sub.τ 210 and a subset of τ-inactive faults which do not cut H.sub.τ.

(170) A geologist or other user may decide to manually transfer some faults from the τ-inactive fault set to the τ-active set, or vice versa, which typically causes greater restoration deformations. For example, manually altering the set of automatically computed τ-active and τ-inactive faults typically makes the restoration less accurate.

(171) Four new 3D piecewise continuous discrete functions {u.sub.τ, v.sub.τ, t.sub.τ, ε.sub.τ}.sub.r may be created that are defined on grid Γ 100 whose discontinuities occur only across τ-active faults;

(172) Referring to FIG. 3, to remove discontinuities of discrete functions {u.sub.τ, v.sub.τ, t.sub.τ, ε.sub.τ}.sub.r across τ-inactive faults, for all τ-inactive faults F 300, one or more of the following inventive (e.g., DSI) constraints may be installed on Γ 100, e.g., as:

(173) $\begin{array}{cc}\begin{array}{cc}1)& u_{\tau }\left(r_F^{\oplus }\right)=u_{\tau }\left(r_F^{\ominus }\right)\\ 2)& v_{\tau }\left(r_F^{\oplus }\right)=v_{\tau }\left(r_F^{\ominus }\right)\\ 3)& t_{\tau }\left(r_F^{\oplus }\right)=t_{\tau }\left(r_F^{\ominus }\right)\end{array}\}\forall {\left(r_F^{\oplus },r_F^{\ominus }\right)}_{\tau }& \left(25\right)\\ \begin{array}{cc}1)& \mathrm{grad}{u}_{\tau }\left(r_F^{\oplus }\right)=\mathrm{grad}{u}_{\tau }\left(r_F^{\ominus }\right)\\ 2)& \mathrm{grad}{v}_{\tau }\left(r_F^{\oplus }\right)=\mathrm{grad}{v}_{\tau }\left(r_F^{\ominus }\right)\\ 3)& \mathrm{grad}{t}_{\tau }\left(r_F^{\oplus }\right)=\mathrm{grad}{t}_{\tau }\left(r_F^{\ominus }\right)\end{array}\}\forall {\left(r_F^{\oplus },r_F^{\ominus }\right)}_{\tau }& \left(26\right)\\ \begin{array}{cc}1)& {\mathrm{.Math.}}_{\tau }\left(r_F^{\oplus }\right)={\mathrm{.Math.}}_{\tau }\left(r_F^{\ominus }\right)\\ 2)& \mathrm{grad}{\mathrm{.Math.}}_{\tau }\left(r_F^{\oplus }\right)=\mathrm{grad}{\mathrm{.Math.}}_{\tau }\left(r_F^{\ominus }\right)\end{array}\}\forall {\left(r_F^{\oplus },r_F^{\ominus }\right)}_{\tau }& \begin{array}{c}\left(27\right)\\ \left(28\right)\end{array}\end{array}$
where (r.sub.F.sup.⊕,r.sub.F.sup.⊖), (304,306) represents a pair of “mate-points” collocated on both sides of F 300 and assigned to F.sup.+ 103 and F.sup.− 104, respectively, and ε.sub.τ(r) represents an error correction constraint. Constraints (25), (26), (27) and (28) may be referred to collectively as “fault transparency constraints.”

(174) Assuming that TH.sub.min>0 is a given threshold chosen by a geologist or other user, fault transparency constraints (25), (26), (27) and (28) may be locally installed at any point r.sub.F on a τ-active fault F where fault throw is lower than TH.sub.min. As a non-limitative example, TH.sub.min may be equal to 1 meter.

(175) Two new discrete vector fields r* and R.sub.τ may be defined on 3D grid Γ 100.

(176) For each node α∈Γ 107: r*(α) may be initialized as the initial location of α:
r*(α)=r(α)  (29) a decompaction transformation C.sup.−(r) known in the art may be used to move a vertically downward from its current compacted altitude z(α) to a new decompacted (e.g., deeper) altitude:
r(α)←C.sup.−(r(α))  (30)

Vertical Restoration Coordinate t.SUB.τ.(r)

(177) Referring to FIG. 2 and FIG. 7, the depositional uvt-transform 700 of G.sub.τ 202 is typically correct when equation (3) is honored. Based on this observation, embodiments of the present invention adapt equation (3) for the inventive restoration technique, replacing the vertical depositional coordinate t(r) by a vertical restoration coordinate t.sub.τ(r) and replacing equation (3) by the following inventive thickness-preserving constraint to ensure layer thickness is preserved and surfaces {S.sub.τ(d):d≥0} are parallel:
∥grad t.sub.τ(r)∥≅1 ∀r∈G.sub.τ  (31)
In addition, to allow H.sub.τ 210 to be restored on surface S.sub.τ(0) 209, the vertical restoration coordinate t.sub.τ(r) may honor the following boundary condition, e.g., as a DSI constraint on grid Γ 100, referred to as the “paleo-sea-floor constraint”:
t.sub.τ(r.sub.H)=z.sub.τ.sup.O ∀r.sub.H∈H.sub.τ  (32)

(178) Due to its non-linearity, thickness-preserving equation (31) cannot be implemented as a DSI constraint, which must be linear. In order to incorporate the thickness-preserving equation into the restoration model using the DSI method, various linear surrogates of equation (31) may be used to approximate t.sub.1(r) as follows: Referring to FIG. 1, to approximate thickness-preserving equation (31), as a non-limitative example, install the following DSI constraints on the grid Γ 100 where r.sub.T⋄ and r.sub.T* are arbitrary points belonging to any pair (T.sup.⋄, T*) of adjacent cells 108 of grid Γ 100 and where N(r.sub.h) is the field of unit vectors defined on H.sub.τ, orthogonal to H.sub.τ and oriented in the direction of younger terrains:
1) grad t.sub.τ(r.sub.H)=N(r.sub.H) ∀r.sub.H∈H.sub.τ  (33)
2) grad t.sub.τ(r.sub.T⋄)≅grad t.sub.τ(r.sub.T*) ∀(T.sup.⋄,T*) Referring to FIG. 1, to approximate thickness-preserving equation (31), alternatively as a non-limitative example, install the inventive DSI constraints on grid Γ 100 as follows:

(179) $\begin{array}{cc}\begin{array}{ccc}1)& \mathrm{grad}{t}_{\tau }\left(r\right)=\frac{\mathrm{grad}t\left(r\right)}{.Math.\mathrm{grad}t\left(r\right).Math.}& \forall r\in G_{\tau }\\ 2)& \mathrm{grad}{t}_{\tau }\left(r_{T\mathrm{.diamond-solid.}}\right)\simeq \mathrm{grad}{t}_{\tau }\left(r_{T^{*}}\right)& \forall \left(T^{\mathrm{.diamond-solid.}},T^{*}\right)\end{array}& \left(34\right)\end{array}$
where r.sub.T⋄ and r.sub.T* are arbitrary points belonging to any pair (T.sup.⋄, T*) of adjacent cells of grid Γ 100 (e.g., the centers of T.sup.⋄ and T*, respectively).

(180) Constraints (33) and (34) are only examples of possible surrogate-thickness-preserving constraints. Other examples of such surrogate thickness-preserving constraints may be used.

(181) Referring to FIG. 3, contrary to constraint (33)(1), new inventive constraint (34)-(1) benefits from the geologic observation that, throughout the entire domain G.sub.τ 202, surfaces {S.sub.τ(d):d≥0} 208 generally have a shape roughly similar to the level sets of the geologic time of deposition t(r);

(182) Assuming that constraints (32) and (33) or (34) are installed on grid Γ 100, a first approximation of vertical restoration coordinate t′.sub.τ(r) may be computed by running the DSI method on grid Γ 100.

(183) Honoring constraint (31) significantly increases the accuracy of the restoration model and a violation of this constraint not only degrades the accuracy of the vertical restoration coordinate t.sub.τ(r) but also the horizontal restoration coordinates {u.sub.τ(r), v.sub.τ(r)} as they are linked to t.sub.τ(r) (e.g., by equations (22) and (23)). Accordingly, there is a great need for validating any approximation technique used to compute t.sub.τ(r).

(184) To test the accuracy of the various approximations of t.sub.τ(r), an example geological terrain is provided in FIG. 4. Despite the apparent simplicity of this terrain, because the thicknesses of the layers vary, this test example is challenging and useful in comparing the accuracy of inventive embodiments with other conventional techniques.

(185) FIG. 5 shows histograms 501 and 502 of the distributions of ∥grad t.sub.τ∥, where t.sub.τ is approximated using constraints (33) or (34), respectively, in the example geological terrain G.sub.τ 202 shown in FIG. 4. FIG. 5 shows that when t.sub.τ is approximated by constraints (33) or (34), ∥grad t.sub.τ∥ significantly differs from “1”. Therefore, while constraints (34) provide a better approximation of the thickness-preserving equation (31) than constraints (33), both of these approximations are inaccurate.

(186) Similarly, FIG. 9 shows histograms 901 and 902 of relative variations of volume ΔV/V induced by the restoration of H.sub.τ 210 over G.sub.τ 202 shown in FIG. 4, where t.sub.τ is approximated using constraints (33) or (34), respectively. Ideally, a restoration transformation should minimize variations in volume ΔV/V from the present day to the restored model. FIG. 9 however shows that a restoration based on constraints (33) or (34) results in a volume variation ΔV/V that significantly differs from the ideal value of “0”. While constraints (34) result in a smaller volume variation ΔV/V than constraints (33), both of these approximations induce a significant volume variation ΔV/V and are inaccurate.

Improving Vertical Restoration Coordinate t.SUB.τ.(r)

(187) An approximation of the vertical restoration coordinate t′.sub.τ(r) may be improved by a “t.sub.τ-incremental improvement” constraint, e.g., as follows:
t.sub.τ(r)=t′.sub.τ(r)+ε.sub.τ(r) ∀r∈G.sub.τ  (35)
where ε.sub.τ(r) is an error correction term, e.g., as characterized below.

(188) Accordingly, assuming that an initial approximation t′.sub.τ(r) has already been obtained, to compute an improved version of t.sub.τ(r), the following inventive incremental procedure may be executed: For each point r.sub.H∈H.sub.τ, set the following equation as an inventive sea-floor-error constraint e.g., using the DSI method:
ε.sub.τ(r.sub.H)=0 ∀r.sub.H∈H.sub.τ  (36) this constraint assumes that constraint (32) remains honored. For each cell T∈Γ∩G.sub.τ 108, choose a point r.sub.τ in the cell T (e.g., its center) and install the new linear thickness-preserving constraint, e.g., using the DSI method as follows:
grad ε.sub.τ(r.sub.τ).Math.grad t′.sub.τ(r.sub.τ)≈½ {1−∥grad t′.sub.τ(r.sub.τ)∥.sup.2}  (37) This constraint is linear, deduced from a linear second order approximation of equation (31). Further, this constraint ensures that, after applying the t.sub.τ-incremental improvement correction constraint (35), the local value of ∥grad t.sub.τ(r)∥ at any point r∈G.sub.τ is as close as possible to “1.” For each sampling point r located on a τ-active fault, install for ε.sub.τ(r) the following inventive DSI constraint referred to as the “t.sub.τ-incremental boundary” constraint:
grad ε.sub.τ(r)×grad t.sub.τ(r)≅0  (38) This constraint specifies that, after applying correction constraint (35), in the close neighborhood of τ-active faults, the shape of level sets of t.sub.τ(r) remains roughly unchanged. To ensure piecewise continuity of the error correction ε.sub.τ(r), install DSI gradient smoothness constraints, known in the art, for the error correction ε.sub.τ(r). Assuming that constraints (36), (37) and (38) are installed on grid Γ 100, to interpolate the error correction ε.sub.τ(r), run DSI on grid Γ 100. For each node α∈Γ 107, update the vertical restoration coordinate t.sub.τ(α) as follows:
t.sub.τ(α)=t′.sub.τ(α)+ε.sub.τ(α)  (39) In the test case represented by FIG. 4 and FIG. 5, the histogram 503 of the distribution of ∥grad t.sub.τ∥, where t.sub.τ is approximated by constraints (37) over G.sub.τ 202 is now considerably better than histograms 501 and 502 obtained with constraints (33) or (34), respectively. In particular: As specified by equation (10), distribution 503 is now centered on value “1”. This condition is of paramount importance to minimize deformations during the restoration process generated by a u.sub.τ v.sub.τ t.sub.τ-transform. The standard deviation of distribution 503 is considerably reduced as compared to the relatively wider standard deviation of distributions 501 and 502. Moreover, in the test case represented in FIG. 4 and FIG. 9, the histogram 903 of the distribution of relative variations of volume ΔV/V induced by a restoration of H.sub.τ 210 over G.sub.τ 202, where t.sub.τ is obtained using constraints (37) is significantly more accurate than histograms 901 and 902, where t.sub.τ is obtained using constraints (33) or (34), respectively. The center of histogram 903 of volume variation ΔV/V is closer to the ideal value “0” than histograms 901 and 902, which indicates that variations of volume are better minimized after applying second order constraints (37) than constraints (33) or (34).

Horizontal Restoration Coordinates {u.SUB.τ.(r), v.SUB.τ.(r)}

(189) Referring to FIG. 2, with respect to surfaces {S.sub.τ(d):d≥0} 208, horizontal restoration coordinates {u.sub.τ(r), v.sub.τ(r)} play a role similar to the one played by paleo-geographic coordinates {u(r), v(r)} with respect to horizons {H.sub.t:t≥0} 216 of the depositional model provided as input. Based on this similarity, horizontal restoration coordinates {u.sub.τ(r), v.sub.τ(r)} may be generated as follows: install equations (20) and (21) as inventive boundary constraints. for all points r∈G.sub.τ 214, define as follows inventive vectors fields a.sub.τ(r) and b.sub.τ(r) respectively, referred to as the “τ-axe” and “τ-coaxe” vector fields:
a.sub.τ(r)=grad t.sub.τ(r)×grad u(r)×grad t.sub.τ(r)
b.sub.τ(r)=grad t.sub.τ(r)×a.sub.τ(r)  (40)
The τ-axe and τ-coaxe vector fields a.sub.τ(r) and b.sub.τ(r) differ considerably from the local axe and co-axe vectors fields a(r) and b(r), e.g., as discussed in U.S. Pat. No. 8,711,140, which is incorporated by reference herein in its entirety. These new τ-axe and τ-coaxe vectors a.sub.τ(r) and b.sub.τ(r) strongly depend on the new vertical restoration coordinate t.sub.τ(r) (e.g., already computed as above) and also take into account the gradient of the paleo-geographic coordinate u(r) (e.g., associated to the depositional model provided as input). if the tectonic style is minimal deformation then, to approximate equations (22), install the following inventive “surrogate minimal-deformation” constraints e.g., using the DSI method:

(190) $\begin{array}{cc}\forall r\in G_{\tau }:|\begin{array}{c}\mathrm{grad}{u}_{\tau }\left(r\right).Math.b_{\tau }\left(r\right)\simeq 0\\ \mathrm{grad}{v}_{\tau }\left(r\right).Math.{\alpha }_{\tau }\left(r\right)\simeq 0\end{array}& \left(41\right)\end{array}$ if the tectonic style is flexural slip then, to approximate equations (23), install the following inventive “surrogate flexural-slip” constraints e.g., using the DSI method:

(191) $\begin{array}{cc}\forall r\in G_{\tau }:|\begin{array}{c}{\mathrm{grad}}_s{u}_{\tau }\left(r\right).Math.b_{\tau }\left(r\right)\simeq 0\\ {\mathrm{grad}}_s{v}_{\tau }\left(r\right).Math.{\alpha }_{\tau }\left(r\right)\simeq 0\end{array}& \left(42\right)\end{array}$ where subscript “S” refers to a projection of the directions of maximal change over iso-value surfaces of the restored vertical coordinate tr. Referring to FIG. 1 and FIG. 6, to prevent the u.sub.t v.sub.τ t.sub.τ-transform used as a restoration operator from generating gaps and overlaps along τ-active faults, specific constraints may be added along fault striae induced by twin-points of the depositional model provided as input. For that purpose, for each pair of twin-points (r.sub.F.sup.+,r.sub.F.sup.−) (101,102) located on faces F.sup.+ 103 and F.sup.− 104 of a τ-active fault F 105, respectively, a process may proceed according to as follows: retrieve the fault stria σ(r.sub.F.sup.−) 600 passing through twin points (r.sub.F.sup.+, r.sub.F.sup.−) (101,102), and on curve σ(r.sub.F.sup.−) 600, retrieve a point {tilde over (r)}.sub.F.sup.+ 601 located on F.sup.+ 103 and such that t.sub.τ({tilde over (r)}.sub.F.sup.+) is approximately equal to t.sub.τ(r.sub.F.sup.−). install the following inventive “τ-fault-striae” constraints e.g., using the DSI method:

(192) $\begin{array}{cc}|\begin{array}{cc}1)& u_{\tau }\left({\tilde{r}}_F^{+}\right)\simeq u_{\tau }\left(r_F^{-}\right)\\ 2)& v_{\tau }\left({\tilde{r}}_F^{+}\right)\simeq v_{\tau }\left(r_F^{-}\right)\end{array}& \left(43\right)\end{array}$ To ensure piecewise continuity of horizontal restoration coordinates u.sub.τ(r) and v.sub.τ(r), install gradient smoothness constraints e.g., using the DSI method. To compute the pair of horizontal restoration coordinates u.sub.τ(r), v.sub.τ(r)) honoring constraints (20), (41 or 42) and (43), run DSI on grid Γ 100.

Computing the Restoration R.SUB.τ.(r)

(193) The restoration vector field R.sub.τ(r) represents the field of deformation vectors from the present day (e.g., xyz) space to the restoration (e.g., u.sub.t v.sub.τ t.sub.τ) space, e.g., computed from the u.sub.τ v.sub.τ t.sub.τ-transform.

(194) Referring to FIG. 1, for each node α 107 of 3D grid Γ 100, move a to restored location r(α), e.g., defined as follows:
r(α)=u.sub.τ(α).Math.r.sub.x+v.sub.τ(α).Math.r.sub.y+t.sub.τ(α).Math.r.sub.z  (44)

(195) For each node α 107 of 3D grid Γ 100: if, to compute vertical restoration coordinate t.sub.τ(r), compaction was taken into account, then, using a recompaction operator C.sup.+(r) known in the art, move α vertically upward from its current decompacted altitude z(α) to a new recompacted (shallower) altitude:
r(α)←C.sup.+(r(α))  (45) save the restoration vector Rτ(α) on a digital device:
Rτ(α)=r(α)−r*(α)  (46) where r*(α) is defined e.g., in equation (29). reset location r(α) of a to its initial location before restoration:
r(α)←r*(α)  (47) stop.

Scanning the Subsurface Through Time

(196) Consider a series of geological restoration times {τ.sub.1<r.sub.2< . . . <r.sub.n} associated with reference horizons H.sub.τ.sub.1, H.sub.τ.sub.2, . . . , H.sub.τ.sub.n, respectively. Using the restoration method described herein, for each (τ.sub.i=τ), build and store on a digital device a restoration vector field R.sub.τ.sub.i(r)=R.sub.τ(r), e.g., as:

(197) $\begin{array}{cc}\begin{array}{ccccccc}{\tau }_1& <& {\tau }_2& <& \mathrm{.Math.}& <& {\tau }_n\\ \updownarrow & & \updownarrow & & & & \updownarrow \\ R_{{\tau }_1}& & R_{{\tau }_2}& & & & R_{{\tau }_n}\end{array}& \left(48\right)\end{array}$

(198) In addition to these reference restoration times, an additional restoration time τ.sub.n+1 may be added to be associated with the horizontal plane H.sub.t.sub.n+1 located at a constant altitude z.sub.τ.sub.n+1.sup.0 of the sea level. Time τ.sub.n+1 may be the present day geological time and, provided that τ.sub.n+1 is greater than τ.sub.n, any arbitrary value may be chosen for τ.sub.n+1. As a non-limitative example, τ.sub.n+1 may be defined as:
τ.sub.n+1=τ.sub.n+1  (49)

(199) Because τ.sub.n+1 is the present day, applying the restoration vector field Rτ.sub.n+1(r) to the present day horizon H.sub.t.sub.n+1 should leave H.sub.t.sub.n+1 unchanged e.g., as follows:
Rτ.sub.n+1(r)=0 ∀r∈G  (50)

(200) To explore subsurface evolution throughout geological times, a process may proceed as follows: as input, read a depositional (e.g., GeoChron) model and associated series of restoration vector fields {R.sub.τ.sub.1, R.sub.τ.sub.2, . . . , R.sub.τ.sub.n+1} stored on a digital device; using an input device such as, in a non-limitative example, the keyboard of a computer or the wheel of a computer mouse, select a geological time τ.sub.i in the given list of geological times {τ.sub.1<τ.sub.2< . . . <τ.sub.n+1}; for each vertex α∈Γ 107, save a copy r*(α) of the location of this node in the depositional model given as input; apply the restoration vector field R.sub.τ.sub.i (r) to the depositional model given as input; display the restored model on a device such as, in a non-limitative example, a display (e.g., display 180 of FIG. 15), such as, a screen, a hologram or a 3D printer; repeat the previous operations as long as desired. optionally, to modify the geometry of the horizons at geological time τ.sub.i, use a computerized tool known in the art to edit the geological time of deposition t(r); for each vertex a 107 of 3D grid Γ 100, use copy r*(α) to restore r(α) to its present day location:
r(α)←r*(α) ∀αℑΓ  (51) such an operation implicitly and automatically propagates the modifications of the geometry of horizons optionally performed above; return to the first step above.

Geological Tomography

(201) Geological models are generated using geological or seismic tomography technology. Geological tomography generates an image of the interior subsurface of the Earth based on geological data collected by transmitting a series of incident waves and receiving reflections of those waves across discontinuities in the subsurface. A transmitter may transmit signals, for example, acoustic waves, compression waves or other energy rays or waves, that may travel through subsurface structures. The transmitted signals may become incident signals that are incident to subsurface structures. The incident signals may reflect at various transition zones or geological discontinuities throughout the subsurface structures, such as, faults or horizons. The reflected signals may include seismic events. A receiver may collect data, for example, reflected seismic events. The data may be sent to a modeling mechanism that may include, for example, a data processing mechanism and an imaging mechanism.

(202) Reference is made to FIG. 14, which is a schematic illustration of a geological tomography technique in which a series of incident rays 111 and reflected rays 121 are propagated through a subsurface region of the Earth 30 to image the subsurface, according to an embodiment of the invention.

(203) One or more transmitter(s) (e.g., 190 of FIG. 15) located at incident location(s) 60 may emit a series of incident rays 111. Incident rays 111 may include for example a plurality of energy rays related to signal waves, e.g., sonic waves, seismic waves, compression waves, etc. Incident rays 111 may be incident on, and reflect off of, a subsurface structure or surface 90 at a reflection point 50. Multiple reflection points 50 may be identified or imaged or displayed in conjunction to display, for example, a horizon.

(204) One or more receiver(s) (e.g., 140 of FIG. 15) located at reflected location(s) 65 may receive the reflection rays 121. Reflection rays 121 may be the reflected images of incident rays 111, for example, after reflecting off of image surface 90 at target point 50. The angle of reflection 55 may be the angle between corresponding incident rays 111 and reflected rays 121 at reflection point 50. An incident rays 111 and a corresponding reflected rays 121 may propagate through a cross-section of a subsurface structure 30. Incident rays 111 may reflect off of a subsurface feature 90 at a reflection point 50, for example, a point on an underground horizon, the seafloor, an underground aquifer, etc.

(205) One or more processor(s) (e.g., 140 of FIG. 15) may reconstitute incident and reflected rays 111 and 121 to generate an image the subsurface 30 using an imaging mechanism. For example, a common reflection angle migration (CRAM) imaging mechanism may image reflection points 50 by aggregating all reflected signals that may correspond to a reflection point, for example, reflected signals that may have the same reflection angle. In other examples, imaging mechanisms may aggregate reflected signals that may have the same reflection offset (distance between transmitter and receiver), travel time, or other suitable conditions.

(206) The processor(s) may compose all of the reflection points 50 to generate an image or model of the present day underground subsurface of the Earth 30. The processor(s) may execute a restoration transformation (e.g., u.sub.τ v.sub.τ t.sub.τ-transform 201) to transform the present day model of subsurface 30 to a restored subsurface image 203 at a restoration time τ. One or more display(s) (e.g., 180 of FIG. 15) may visualize the present day subsurface image 30 and/or the restored subsurface image 203.

System Overview

(207) Reference is made to FIG. 15, which schematically illustrates a system including one or more transmitter(s), one or more receiver(s) and a computing system in accordance with an embodiment of the present invention. Methods disclosed herein may be performed using a system 1505 of FIG. 15.

(208) System 1505 may include one or more transmitter(s) 190, one or more receiver(s) 120, a computing system 130, and a display 180. The aforementioned data, e.g., seismic data used to form intermediate data and finally to model subsurface regions, may be ascertained by processing data generated by transmitter 190 and received by receiver 120. Intermediate data may be stored in memory 150 or other storage units. The aforementioned processes described herein may be performed by software 160 being executed by processor 140 manipulating the data.

(209) Transmitter 190 may transmit signals, for example, acoustic waves, compression waves or other energy rays or waves, that may travel through subsurface (e.g., below land or sea level) structures. The transmitted signals may become incident signals that are incident to subsurface structures. The incident signals may reflect at various transition zones or geological discontinuities throughout the subsurface structures. The reflected signals may include seismic data.

(210) Receiver 120 may accept reflected signal(s) that correspond or relate to incident signals, sent by transmitter 190. Transmitter 190 may transmit output signals. The output of the seismic signals by transmitter 190 may be controlled by a computing system, e.g., computing system 130 or another computing system separate from or internal to transmitter 190. An instruction or command in a computing system may cause transmitter 190 to transmit output signals. The instruction may include directions for signal properties of the transmitted output signals (e.g., such as wavelength and intensity). The instruction to control the output of the seismic signals may be programmed in an external device or program, for example, a computing system, or into transmitter 190 itself.

(211) Computing system 130 may include, for example, any suitable processing system, computing system, computing device, processing device, computer, processor, or the like, and may be implemented using any suitable combination of hardware and/or software. Computing system 130 may include for example one or more processor(s) 140, memory 150 and software 160. Data 155 generated by reflected signals, received by receiver 120, may be transferred, for example, to computing system 130. The data may be stored in the receiver 120 as for example digital information and transferred to computing system 130 by uploading, copying or transmitting the digital information. Processor 140 may communicate with computing system 130 via wired or wireless command and execution signals.

(212) Memory 150 may include cache memory, long term memory such as a hard drive, and/or external memory, for example, including random access memory (RAM), read only memory (ROM), dynamic RAM (DRAM), synchronous DRAM (SD-RAM), flash memory, volatile memory, non-volatile memory, cache memory, buffer, short term memory unit, long term memory unit, or other suitable memory units or storage units. Memory 150 may store instructions (e.g., software 160) and data 155 to execute embodiments of the aforementioned methods, steps and functionality (e.g., in long term memory, such as a hard drive). Data 155 may include, for example, raw seismic data collected by receiver 120, instructions for building a mesh (e.g., 100), instructions for partitioning a mesh, and instructions for processing the collected data to generate a model, or other instructions or data. Memory 150 may also store instructions to divide and model τ-active faults and τ-inactive faults. Memory 150 may generate and store the aforementioned constraints, restoration transformation (e.g., u.sub.τ v.sub.τ t.sub.τ-transform 201), restoration coordinates (e.g., u.sub.τ, v.sub.τ, t.sub.τ), a geological-time and paleo-geographic coordinates (e.g., u, v, t), a model representing a structure when it was originally deposited (e.g., in uvt-space), a model representing a structure at an intermediate restoration time (e.g., in u.sub.τ, v.sub.τ, t.sub.τ-space), and/or a model representing the corresponding present day structure in a current time period (e.g., in xyz-space). Memory 150 may store cells, nodes, voxels, etc., associated with the model and the model mesh. Memory 150 may also store forward and/or reverse u.sub.τ, v.sub.τ, t.sub.τ-transformations to restore present day models (e.g., in xyz-space) to restored models (e.g., in u.sub.t, v.sub.τ, t.sub.τ-space), and vice versa. Memory 150 may also store the three-dimensional restoration vector fields, which when applied to the nodes of the initial present day model, move the nodes of the initial model to generate one of the plurality of restored models. Applying a restoration vector field to corresponding nodes of the present day model may cause the nodes to “move”, “slide”, or “rotate”, thereby transforming modeled geological features represented by nodes and cells of the initial model. Data 155 may also include intermediate data generated by these processes and data to be visualized, such as data representing graphical models to be displayed to a user. Memory 150 may store intermediate data. System 130 may include cache memory which may include data duplicating original values stored elsewhere or computed earlier, where the original data may be relatively more expensive to fetch (e.g., due to longer access time) or to compute, compared to the cost of reading the cache memory. Cache memory may include pages, memory lines, or other suitable structures. Additional or other suitable memory may be used.

(213) Computing system 130 may include a computing module having machine-executable instructions. The instructions may include, for example, a data processing mechanism (including, for example, embodiments of methods described herein) and a modeling mechanism. These instructions may be used to cause processor 140 using associated software 160 modules programmed with the instructions to perform the operations described. Alternatively, the operations may be performed by specific hardware that may contain hardwired logic for performing the operations, or by any combination of programmed computer components and custom hardware components.

(214) Embodiments of the invention may include an article such as a non-transitory computer or processor readable medium, or a computer or processor storage medium, such as for example a memory, a disk drive, or a USB flash memory, encoding, including or storing instructions, e.g., computer-executable instructions, which when executed by a processor or controller, carry out methods disclosed herein.

(215) Display 180 may display data from transmitter 190, receiver 120, or computing system 130 or any other suitable systems, devices, or programs, for example, an imaging program or a transmitter or receiver tracking device. Display 180 may include one or more inputs or outputs for displaying data from multiple data sources or to multiple displays. For example, display 180 may display visualizations of subsurface models including subsurface features, such as faults, horizons and unconformities, as a present day subsurface image (e.g., 202), a restored subsurface image (e.g., 203) and/or a depositional model (e.g., 703). Display 180 may display one or more present day model(s), depositional model(s), restoration model(s), as well as a series of chronologically sequential restoration models associated with a sequence of respective restoration times (e.g., τ.sub.1<τ.sub.2<τ.sub.3<τ.sub.4, as shown in FIG. 13). The models may be displayed one at a time, two at a time, or many at a time (e.g., the number selected by a user or automatically based on the difference between models or the total number of models). Display 180 may display the models in a sequence of adjacent models, through which a user may scan (e.g., by clicking a ‘next’ or ‘previous’ button with a pointing device such as a mouse or by scrolling through the models).

(216) Input device(s) 165 may include a keyboard, pointing device (e.g., mouse, trackball, pen, touch screen), or cursor direction keys, for communicating information and command selections to processor 140. Input device 165 may communicate user direction information and command selections to the processor 140. For example, a user may use input device 165 to select one or more preferred models from among the plurality of perturbed models, recategorize faults as τ-active faults and τ-inactive, or edit, add or delete subsurface structures.

(217) Processor 140 may include, for example, one or more processors, controllers or central processing units (“CPUs”). Software 160 may be stored, for example, in memory 150. Software 160 may include any suitable software, for example, DSI software.

(218) Processor 140 may generate a present day subsurface image (e.g., 202), a restored subsurface image (e.g., 203) and/or a depositional model (e.g., 703), for example, using data 155 from memory 150. In one embodiment, a model may simulate structural, spatial or geological properties of a subsurface region, such as, porosity or permeability through geological terrains.

(219) Processor 140 may initially generate a three dimensional mesh, lattice, grid or collection of nodes (e.g., 100) that spans or covers a domain of interest. The domain may cover a portion or entirety of the three-dimensional subsurface region being modeled. Processor 140 may automatically compute the domain to be modeled and the corresponding mesh based on the collected seismic data so that the mesh covers a portion or the entirety of the three-dimensional subsurface region from which geological data is collected (e.g., the studied subsurface region). Alternatively or additionally, the domain or mesh may be selected or modified by a user, for example, entering coordinates or highlighting regions of a simulated optional domain or mesh. For example, the user may select a domain or mesh to model a region of the Earth that is greater than a user-selected subsurface distance (e.g., 100 meters) below the Earth's surface, a domain that occurs relative to geological features (e.g., to one side of a known fault or riverbed), or a domain that occurs relative to modeled structures (e.g., between modeled horizons H(t.sub.1) and H(t.sub.100)). Processor 140 may execute software 160 to partition the mesh or domain into a plurality of three-dimensional (3D) cells, columns, or other modeled data (e.g., represented by voxels, pixels, data points, bits and bytes, computer code or functions stored in memory 150). The cells or voxels may have hexahedral, tetrahedral, or any other polygonal shapes, and preferably three-dimensional shapes. Alternatively, data may include zero-dimensional nodes, one-dimensional segments, two-dimensional facet and three-dimensional elements of volume, staggered in a three-dimensional space to form three-dimensional data structures, such as cells, columns or voxels. The cells preferably conform to and approximate the orientation of faults and unconformities. Each cell may include faces, edges and/or vertices. Each cell or node may correspond to one or more particles of sediment in the Earth (e.g., a node may include many cubic meters of earth, and thus many particles).

(220) Data collected by receiver 120 after the time of deposition in a current or present time period, include faults and unconformities that have developed since the original time of deposition, e.g., based on tectonic motion, erosion, or other environmental factors, may disrupt the regular structure of the geological domain. Accordingly, an irregular mesh may be used to model current geological structures, for example, so that at least some faces, edges, or surfaces of cells are oriented parallel to faults and unconformities, and are not intersected thereby. In one embodiment, a mesh may be generated based on data collected by receiver 120, alternatively, a generic mesh may be generated to span the domain and the data collected by receiver 120 may be used to modify the structure thereof. For example, the data collected may be used to generate a set of point values at “sampling point”. The values at these points may reorient the nodes or cells of the mesh to generate a model that spatially or otherwise represents the geological data collected from the Earth. Other or different structures, data points, or sequences of steps may be used to process collected geological data to generate a model. The various processes described herein (e.g., restoring a geological model using τ-active and τ-inactive faults, or restoring a geological model using a new thickness-preserving constraint) may be performed by manipulating such modeling data.

(221) Restoration coordinates may be defined at a finite number of nodes or sampling points based on real data corresponding to a subsurface structure, e.g., one or more particles or a volume of particles of Earth. Restoration coordinates may be approximated between nodes to continuously represent the subsurface structure, or alternatively, depending on the resolution in which the data is modeled may represent discrete or periodic subsurface structures, e.g., particles or volumes of Earth that are spaced from each other.

(222) The computing system of FIG. 15 may accept the data used in the operations of FIGS. 16, 17 and 19 as for example a set of data generated by tomographic scanning of a subsurface geological region of the Earth as disclosed in reference to FIG. 14, or such data augmented by another process. The computing system may accept one or more of seismic and well data. The computing device may generate one or more of seismic and well data.

(223) “Restoration” or “intermediate” time τ may refer to a time in the past before the present day and after a time when an oldest or deepest horizon surface in the 3D model was deposited. “Restoration” or “intermediate” transformation or model may refer to a model or image of the surface as it was configured at the “intermediate” time in the past τ. An intermediate horizon may refer to a horizon that was deposited at the “intermediate” time τ, which is located above the deepest horizon and below the shallowest horizon.

(224) “Time” including the present-day, current or present time, the past restoration time τ, and/or the depositional time t, may refer to geological time periods that span a duration of time, such as, periods of thousands or millions of years.

(225) “Geological-time” t(r) may refer to the time of deposition when a particle of sediment represented by point r was originally deposited in the Earth. In some embodiments, the geological-time of the deposition may be replaced, e.g., by any arbitrary monotonic increasing function of the actual geological-time. It is a convention to use an monotonically increasing function, but similarly an arbitrary monotonic decreasing function may be used. The monotonic function may be referred to as the “pseudo-geological-time”.

(226) The geological-time of the deposition and restoration time of particles are predicted approximate positions since past configurations can not typically be verified.

(227) “Current” or “present day” location for a particle (or data structure representing one or more particles) or subsurface feature may refer to the location of the item in the present time, as it is measured.

(228) In stratified terrain, layers, horizons, faults and unconformities may be curvilinear surfaces which may be for example characterized as follows. A horizon, Hτ, may be a surface corresponding to a plurality of particles of sediment which were deposited approximately at substantially the same geological-time, τ. A fault may be a surface of discontinuity of the horizons that may have been induced by a relative displacement of terrains on both sides of such surfaces. In other words, the geological-time of deposition of the sediments is discontinuous across each fault. Faults may cut horizons and may also cut other faults. An unconformity may be a surface of discontinuity of the horizons that may have been induced by for example an erosion of old terrains replaced by new ones. In other words, similarly to faults, the geological-time of deposition of the sediments is discontinuous across each unconformity.

(229) Terrain deformed in the neighborhood of a point r in the G-space may occur according to a “minimal deformation” tectonic style when, in this neighborhood: the strain tensor is approximately equal to the null tensor.

(230) Terrain deformed in the neighborhood of a point r in the G-space may occur according to a “flexural slip” tectonic style when, in this neighborhood: the length of any small increment of distance d(r) on the horizon passing through point r is preserved, e.g., in any direction, and, the volume of the terrains in the neighborhood of point r does not vary.

(231) Discrete-Smooth-Interpolation (DSI) is a method for interpolating or approximating values of a function ƒ(x,y,z) at nodes of a 3D grid or mesh Γ (e.g., 100), while honoring a given set of constraints. The DSI method allows properties of structures to be modeled by embedding data associated therewith in a (e.g., 3D Euclidean) modeled space. The function ƒ(x,y,z) may be defined by values at the nodes of the mesh, F. The DSI method allows the values of ƒ(x,y,z) to be computed at the nodes of the mesh, F, so that a set of one or more (e.g., linear) constraints are satisfied. DSI generally only applies linear constraints on the model.

(232) In some embodiments, bold symbols represent vectors or multi-dimensional (e.g., 3D) functions or data structures.

(233) In some embodiments, the “simeq” symbol “≅” or “≅” may mean approximately equal to, e.g., within 1-10% of, or in a least squares sense. In some embodiments, the symbol “≡” may mean identical to, or defined to be equal to.

(234) While embodiments of the invention describe the input depositional model as the GeoChron model, any other depositional model visualizing the predicted configuration of each particle, region or layer at its respective time of depositional may be used.

(235) While embodiments of the invention describe the present day coordinates as xyz, the restoration coordinates as u.sub.τv.sub.τt.sub.τ, the depositional coordinates as uvt, the restoration transformation as a u.sub.τv.sub.τt.sub.τ-transform, and the depositional transformation as a uvt-transform, any other coordinates or transformations may be used.

(236) In the foregoing description, various aspects of the present invention have been described. For purposes of explanation, specific configurations and details have been set forth in order to provide a thorough understanding of the present invention. However, it will also be apparent to one skilled in the art that the present invention may be practiced without the specific details presented herein. Furthermore, well known features may have been omitted or simplified in order not to obscure the present invention. Unless specifically stated otherwise, as apparent from the following discussions, it is appreciated that throughout the specification discussions utilizing terms such as “processing,” “computing,” “calculating,” “determining,” or the like, refer to the action and/or processes of a computer or computing system, or similar electronic computing device, that manipulates and/or transforms data represented as physical, such as electronic, quantities within the computing system's registers and/or memories into other data similarly represented as physical quantities within the computing system's memories, registers or other such information storage, transmission or display devices. In addition, the term “plurality” may be used throughout the specification to describe two or more components, devices, elements, parameters and the like.

(237) Embodiments of the invention may manipulate data representations of real-world objects and entities such as underground geological features, including faults and other features. The data may be generated by tomographic scanning, as discussed in reference to FIG. 14, e.g., received by for example a receiver receiving waves generated e.g., by an air gun or explosives, that may be manipulated and stored, e.g., in memory 150 of FIG. 15, and data such as images representing underground features may be presented to a user, e.g., as a visualization on display 180 of FIG. 15.

(238) When used herein, a subsurface image or model may refer to a computer-representation or visualization of actual geological features such as horizons and faults that exist in the real world. Some features when represented in a computing device may be approximations or estimates of a real world feature, or a virtual or idealized feature, such as an idealized horizon as produced in a u.sub.τ v.sub.τ t.sub.τ-transform. A model, or a model representing subsurface features or the location of those features, is typically an estimate or a “model”, which may approximate or estimate the physical subsurface structure being modeled with more or less accuracy.

(239) It will thus be seen that the objects set forth above, among those made apparent from the preceding description, are efficiently attained and, because certain changes may be made in carrying out the above method and in the construction(s) set forth without departing from the spirit and scope of the invention, it is intended that all matter contained in the above description and shown in the accompanying drawings shall be interpreted as illustrative and not in a limiting sense.

(240) It is also to be understood that the following claims are intended to cover all of the generic and specific features of the invention herein described and all statements of the scope of the invention which, as a matter of language, might be said to fall therebetween.