INTRODUCTION

The prediction of subseismic scale faults is important for the management and exploration of reservoirs and their overburden. Such faults have significant effects on permeability and on fluid flow (Antonellini and Aydin, 1994). Many authors have shown that the quantification of the extension from seismic interpretation underestimates the true extension because many faults and fractures are likely to be too small to be resolved on a seismic scale (e.g., Marrett and Allmendinger, 1992; Reston and McDermott, 2014). Although seismic-scale faults have displacements between 1 km (3281 ft) and 10 m (33 ft) and borehole core-scale faults have displacements between 1 mm (0.05 in.) and 10 cm (4 in.) (Needham et al., 1996), the scales from 10 cm (4 in.) to 10 m (33 ft) are not resolvable by geophysical exploration methods.

In the framework of geologic storage of CO₂, our project, Prediction of Deformation to Ensure Carbon Traps (PROTECT), aims to predict subseismic scale faults and fractures within the reservoir and the overburden (Krawczyk et al., 2011, 2015). Studies on the geologic storage of CO₂ have either concentrated on the reservoir itself (e.g., Chadwick et al., 2004; Riding, 2006) or aquifers close to surface (Norden, 2011), but evaluation of the complete overburden and reservoir is rare (Keating et al., 2014; Krawczyk et al., 2015). For instance, the In-Salah project in central Algeria integrated numerous methods, such as three-dimensional (3-D) seismic data, satellite imaging, and seismologic studies (e.g., Mathieson et al., 2010; Stork et al., 2014), as well as various modeling techniques, such as static geologic models, multiphase flow models, fracture-flow models, and geomechanical models (e.g., Ringrose et al., 2013), to determine small-scale faulting. These results were used to support operational decisions, but none of these models included subseismic scale fractures. In the Sleipner field (Norway), well log and cuttings, seismic attribute analysis, geomechanical modeling, and time-lapse gravity surveys have been carried out to make sure that the CO₂ plume is not migrating into the overburden (Chadwick et al., 2004; Eiken et al., 2011). In the International Energy Agency Weyburn Project (central Canada), which is one of the largest enhanced oil recovery–carbon capture and storage (CCS) projects, numerous methods, such as time-lapse 3-D seismic, reservoir simulations, geochemical modeling, regional geoscience investigations, and wellbore integrity studies, were applied to characterize the reservoir and overburden (Whittaker et al., 2011). However, as Eiken et al. (2011) noted, dynamic modeling is a challenge and there is still room for further model improvement.

Project PROTECT combined seismic multiattributes (Trappe and Hellmich, 2003; Krawczyk et al., 2015), forward modeling (Aruffo et al., 2014), and 3-D retrodeformation results. The project cooperates with Australia’s first demonstration site of deep geologic storage of CO₂ in the Otway Basin (southwestern Victoria, Australia; Cook, 2014). This paper is primarily concerned with the 3-D retrodeformation of the interpretation of high-quality 3-D seismic data (Ziesch et al., 2017) and the internal strain that occurred as a consequence of restoration. The calculated strain results are used as a proxy to predict the location of subseismic scale fractures and therefore potential CO₂ migration pathways. The insights pertaining to the 3-D retrodeformation method used here can be potentially implemented in a number of applications, especially subsurface storage or exploitation.

GEOLOGIC SETTING

The study area (∼8 km [∼5 mi] ∼7 km [∼4 mi]) is located in southwestern Victoria, 200 km (125 mi) southwest of Melbourne, Australia (Figure 1). It is part of the onshore western Otway Basin, an extensional basin that developed along the southern margin of Australia during the Late Jurassic because of the separation of Australia and Antarctica (Watson et al., 2004). In the mid-Late Cretaceous the tectonic setting changed to compressional deformation with associated uplift and fault inversion (Norvick and Smith, 2001). After a 2.5-m.y. hiatus, the last deformation phase began with fast continental spreading in the Southern Ocean during the Lutetian (41 Ma; Lawver et al., 2011) coupled with thermal subsidence and marine transgression. Previous faults were reactivated under northeast-southwest extension (Norvick and Smith, 2001), in combination with syntectonic deposition. The faults imaged in the seismic data were all formed during the first phase of deformation and were still active during the third phase (Ziesch et al., 2017).

Figure 1. Overview map of the regional tectonics, showing major normal faults in the Otway Basin (compiled from Douglas and Ferguson, 1988; Alley and Lindsay, 1995; Perincek and Cockshell, 1995; Moore et al., 2000). Cretaceous sedimentary rocks marked in gray. Close-up of the study area, showing borehole locations (Cooperative Research Centre [CRC]-1, -2) and the traces of the profiles (0–9) used for the extension versus age diagram (Figure 11). Major fault is marked.

The stratigraphy of the study area can be divided into eight supersequences that are combined here as five lithostratigraphic groups; they are as follows, from oldest to youngest: the Otway, Sherbrook, Wangerrip, Nirranda, and Heytesbury Groups (Figure 2; Tupper et al., 1993). The stratigraphic horizons interpreted in the seismic volume belong to the youngest four of the five groups (Figure 2). The Otway Group consists of fluvial-lacustrine sedimentary rocks that accumulated during the first rifting phase and are overlain by the Sherbrook Group, which contains an alternating sequence of sandstones and mudstones. The Sherbrook Group includes the reservoir (Waarre Formation [Fm.]) and seals (Flaxman Fm., Belfast Mudstone (Mst.), and Skull Creek Mst.) for CO₂ storage. The beginning of the passive margin sedimentation is represented by the Wangerrip Group, which consists of 700 m (2297 ft) of sandstones and mudstones (Norvick and Smith, 2001; Krassay et al., 2004). After a major hiatus of 8.5 m.y., the Nirranda Group was deposited in proximal and distal environments (Norvick and Smith, 2001). Fine-grained facies of the Narrawaturk Marl overlay sandstones of the thin Mepunga Fm. The Miocene Heytesbury Group was deposited under full marine conditions and is characterized by limestones and marls (Krassay et al., 2004).

Figure 2. Stratigraphy of the study area, modified after Ziesch et al. (2017), compiled from Perincek and Cockshell (1995) and Geoscience Australia and Australian Stratigraphy Commission (2017). The key horizons picked in the three-dimensional seismic data set are shown, as are the ranges of thicknesses between the horizons in parentheses. Fm. = Formation; Mst. = Mudstone; Quat. = Quaternary; Sst. = Sandstone.

Ziesch et al. (2017) presented depth and isochore maps of all the horizons in the study area. Figure 3 shows concise summary maps of the reservoir and seal horizons (i.e., Turonion top Waarre Fm. and Santonian top Skull Creek Mst., respectively). Vertical displacements on faults increase from northeast to southwest, where throws up to 800 m (2625 ft) can be observed (Figure 3A, D). Relay ramps act as soft-linking elements between the major faults. Fault Buttress Southeast and fault 6 form a horst structure; between these two faults, a major change in stratigraphic thickness occurs. In general, the thicknesses increase for both stratigraphic horizons toward the southwest and the normal faults were strongly synsedimentary (Figure 3B, E).

Figure 3. Maps of the two stratigraphic horizons: top Waarre Formation (Fm.) and top Skull Creek Mudstone (Mst.) (modified after Ziesch et al., 2017, and used with permission of Wiley): (A) Depth map of top Waarre Fm. (B) Isochore map of top Waarre Fm. to top Skull Creek Mst. It shows strong fault growth and an increase in thickness toward the southwest. (C) Structural map of top Waarre Fm.; fault heave polygons are shown in white. Black arrows represent the average dip direction, and gray arrows indicate the directions of cylindricity vectors of the faults. See text for details. (D) Depth map of Skull Creek Mst. (E) Isochore map of Skull Creek Mst. to top Timboon Sst. Similar to the underlaying formation, the faults show strong growth thicknesses and the sedimentary thickness increases toward the southwest. (F) Structural map of top Skull Creek Mst. Fault heave polygons are shown in white. Grid coordinate system: Map Grid of Australia, zone 54, Transverse Mercator projection.

DATABASE AND METHOD

The database is a seismic volume (32.3 km [20.1 mi] 14.35 km [8.92 mi] 4.1 s two-way time) with a 24-fold bin size of 20 20 m (66 66 ft). The data were acquired and processed in 2000 (Cooperative Research Centre for Greenhouse Gas Technologies, 2013, personal communication). The used signal frequency range was 5–90 Hz, and the recorded offset range was 28 to 2280 m (92 to 7480 ft). The 3-D poststack finite-difference time migration was used for imaging.

Because the target depth of the data acquisition geometry was the Turonian Waarre Fm. reservoir (∼1.6-km [∼1-mi] depth), the uppermost 400 m (1312 ft) of the data set are less well-imaged compared with the rest of the volume. The reflections of the seal and the overburden formations are well imaged (Figure 4). In total, 8 stratigraphic horizons (including the reservoir, seal, and overburden) and 24 major faults were interpreted in the time domain (Ziesch et al., 2017; Figure 4).

Figure 4. Perspective view of the three-dimensional seismic data with two inlines, one crossline, and a variance depth slice. The rear inline shows interpretation of eight stratigraphic horizons; the crossline shows interpretation of the faults; the front inline section is not interpreted. Fm. = Formation.

For time–depth conversion, interval velocities were determined from sonic velocities measured in the injection well Cooperative Research Centre-1. This information was used to generate velocity gradients for each stratigraphic formation. Quality control of the velocity model relied on nine wells in the larger seismic volume (for more information, see Ziesch et al., 2017).

Structural Modeling Workflow

The deformation analysis required several model building and modeling steps (Figure 5). After importing the depth-converted horizons and faults into the structural modeling software, individual fault blocks were defined (Figure 5). This involved separating the stratigraphic horizons along the major faults. This step was important because it was necessary to be able to move each fault block separately.

Figure 5. Workflow used in this study. After seismic interpretation and model building, the uppermost horizon was decompacted and the next top formation was restored. This procedure was repeated until the restoration of the lowest formation. This generates values of incremental elongation, the paleogeomorphology, and distribution of the internal strain within the fault blocks. 3-D = three-dimensional.

Tetrahedral volumes were built for each fault block between stratigraphic horizons (Figure 6A). The magnitude of strain was calculated by the change in position of all vertices in the tetrahedron relative to their initial, undeformed state (La Pointe et al., 2002; Figure 6B, C).

Figure 6. Example of tetrahedral volumes. (A) Excerpt from the volume model with one fault between two fault blocks; (B) unstrained tetrahedron with a sphere of unity radius; and (C) strained tetrahedron with a strain ellipsoid and its principal strain axes , where e is the strain value.

A two-step approach was used for retrodeformation. First, the uppermost horizon was backstripped; after this, the next top formation was restored, from southwest to northeast. In the following, the methods of decompaction and restoration are explained in detail.

Backstripping and Decompaction

Deposition of younger sedimentary rocks on older strata causes compaction. Thus, the first step during retrodeformation is the sequential removal of the overlying stratigraphic formation and decompacting the underlying sedimentary rocks (Figure 5). Backstripping results in each being compensated for their new position by using the porosity of the formations at depth and the constant c (Table 1). This is derived using Athy’s equation

where φ is the porosity at depth y (km), c is a constant (km⁻¹), and is the initial porosity (Athy, 1930). In this way, the compaction of sedimentary rocks can be attributed to only two parameters, c and (Table 1).

Porosity was derived from a shale volume log . The average value for each stratigraphic horizon was determined, assuming the percentage of material that is not shale is coarser, that is, silt and sand. Parameters for all the interpreted units were calculated using standard values of c and for shale and sandstone (Sclater and Christie, 1980).

Decompaction is a purely vertical process. Faults were included in this process and retrodeformed to their uncompacted dip; that is, in general, the fault dip becomes steeper after each backstripping and decompaction stage (Durand-Riard et al., 2011).

Restoration

The term “restoration” means that a cross section or map is restored back to its undeformed state. This requires accounting for deformation related to tectonic events (Durand-Riard et al., 2011; Chauvin et al., 2015). Since the last century, structural balancing of cross sections has been used to check the consistency of the geologic interpretation (e.g., Chamberlin, 1910; Dahlstrom, 1969; Gibbs, 1983). However, the transfer of structural balancing into the third dimension has been practiced only for the last two decades (e.g., Rouby et al., 2000; Guzofski et al., 2009; Durand-Riard et al., 2011) and will become increasingly important with the increased availability of 3-D seismic data. Retrodeformation of a geologic 3-D model offers new analytical options: the differing amounts of extension throughout the model, the 3-D paleogeomorphology, and the 3-D internal strain that occurred within the retrodeformed fault block volumes (Figure 5).

The complete 3-D model of the study area, including 24 faults and 8 horizons, was restored back into the undeformed state, from southwest to northeast that is, the northeastern edge of the study area was kept fixed during restoration. The reason for this strategy was that the study area is part of the passive margin between southern Australia and Antarctica (Williamson et al., 1990): major faults dip and become younger toward the southwest. All underlying formations were also passively restored during the structural restoration of the top horizon.

The 3-D retrodeformation raises some new challenges with respect to two-dimensional (2-D) restoration.

• Fault shape: in three dimensions, the fault surface can vary laterally (i.e., fault corrugations) and with respect to depth (e.g., listric faults).

• Displacement: faults show variable displacements along strike until the tip point is reached and the displacement approaches zero.

• Fault kinematics: the 3-D movement vector on the fault surface, can be difficult to extract in three dimensions (Ziesch et al., 2015). We calculated mean 3-D kinematic vectors for the faults by using the cylindricity of corrugations of the fault surface (for more information, see Ziesch et al., 2017), but in reality, the faults can have slightly variable kinematic vectors along strike. Only a few degrees difference can lead to an uneven retrodeformed surface.

Compaction Strain versus Deformation Strain

Strain is caused by compaction and structural deformation. With retrodeformation, the strain caused by decompaction and restoration is the reciprocal of the forward strain tensor and therefore a proxy for subseismic scale fractures. This retrodeformation strain can be differentiated into strain by compaction and strain by deformation. It can be expressed by a finite strain ellipsoid (Figure 6) with three principal axes of finite strain with semiaxis lengths: (Ramsay and Huber, 1983).

An initially circular (spherical) marker is transformed during purely vertical decompaction and compaction to a finite strain ellipse (ellipsoid) (Figure 7). Compaction causes uniaxial compression so the resulting finite strain ellipse (ellipsoid) has a shorter axis but a constant axis. For instance, there is a reduction of 53% magnitude after a 53% bulk volume loss (Figure 7A).

Figure 7. Schematic diagram showing the development of strain during compaction and decompaction. (A) Compaction history over time during which a circle with a unit radius undergoes deformation. The resulting finite strain ellipse is characterized by a shorter axis and a constant axis, where are the principal strain axes and e is the strain value. The graph shows the relationship of magnitude over time. (B) The reverse of (A) to show the effect of decompaction; from the finite compaction stage is the reciprocal of from the finite decompaction stage.

The seismic interpretation represents the compacted, present-day state. The model was backstripped and decompacted to the original predeformation state to analyze the compaction-related strain. As in the finite compaction stage is equal to in the finite decompaction stage, decompaction strain is directly related to compaction strain.

Strain caused by faulting can be differentiated into rotation and shearing. Forward deformation of an initially circular (spherical) marker with a shear angle of 90° in the hanging wall of a listric fault will lead to a sheared and rotated finite strain ellipse (ellipsoid) (Figure 8A). However, the geologic interpretation represents the deformed present-day state. By structural restoring of the hanging wall, it is possible to analyze the strain tensors in all tetrahedra of the hanging wall (Figure 6), thereby investigating the distribution of volumetric strain. With respect to the forward deformation, the ellipsoid’s axes caused by structural restoring are mirrored around the vertical axis, but the -strain magnitude correlates with the equivalent forward deformation -strain magnitude (Figure 8B).

Figure 8. Development of internal strain during structural restoration. (A) Forward deformation with a shear angle of 90° in the hanging wall (HW) of a listric fault using inclined shear. The circle undergoes both rotation and shear. (B) The geologic model of the study area represents the deformed present-day state. By restoring the HW, one is able to analyze the strain magnitude, where is the principal major strain axis and e is the strain value. The strain ellipse caused by structural restoration has the same magnitudes, but as in forward deformation, their axes, which are mirrored around a vertical axis, include rotation. The angle β is the shear angle. FW = footwall.

As in all geologic modeling software programs that use volumetric cells, extremely acute angles within the cells, which can be produced in the process of retrodeformation (e.g., during decompaction, behind faults), can cause unnaturally high strain values. In the 3-D model, the faults are corrugated (Ziesch et al., 2015) and the horizon surfaces are tightly faulted; because interpolation and smoothing were avoided, the surfaces appear in the model with an unsmoothed geometry. Because the areas with abnormally high strains can be clearly seen because of their peculiar shape (e.g., circular shape) and location (e.g., directly behind faults), the strain maps are presented here without any modification of this issue and it is addressed in the text so that the readers can judge for themselves.

Inclined Shear Algorithm

The inclined shear algorithm (also known as simple shear) was chosen for structural restoration (Gibbs, 1983; White et al., 1986) because the faults in the study area are either listric (the dip of the listric faults in the south changes from 60°–70° in the upper parts to 30°–40° in the lower parts, e.g., fault Naylor South or fault Naylor East) or subplanar with a corrugated morphology (e.g., fault Buttress Northwest, Ziesch et al., 2017). The hanging wall of a fault block is moved along a 3-D movement vector (see Table 2). An angle β is chosen, on which shearing will occur in the hanging wall during restoration (Hauge and Gray, 1996). If the dip of the fault in the transport direction changes, simple shear is used to compensate the necessary change in the geometry of the hanging wall. The amount of shearing depends on the change in the dip of the fault in the transport direction. The inclined shear algorithm was first developed for 2-D restoration but is available for 3-D restoration in the structural modeling software (Move 2016.2, Midland Valley Exploration Ltd., 2016). However, the inclined shear algorithm should be considered as a 2.5-dimensional or pseudo 3-D solution, because although a 3-D movement vector was used, the inclined shear algorithm divides the 3-D model into 2-D sections parallel to the 3-D vector (Midland Valley Exploration Ltd., 2017, personal communication). In the transport direction, all particles of the hanging wall are displaced with constant heave, but along the fault strike, that is, in the third dimension, the heave amount varies in proportion to the cutoffs of the faulted strata. The method maintains the volume of the hanging wall in 3-D, but it does not maintain bed thickness if shearing occurs (Williams and Vann, 1987).

Seismic interpretation and variance analysis proved that the faults of the study area possess a distinct morphology and that these subvertical corrugations are real (Ziesch et al., 2017). To quantify the amount and orientations of the corrugations, curvature and cylindricity attributes of the fault surfaces were calculated. The resulting vectors of the cylindricity analysis are shown in Figure 3C. We postulated that the corrugations were generated and maintained during the growth of the synsedimentary faults (Ziesch et al., 2017). Therefore, the corrugation vectors could be used as 3-D movement vectors for each fault in the model (see Table 2; Ziesch et al., 2017, 2015). It was determined that the 3-D movement vectors were either pure dip slip or slightly dextrally oblique slip (gray arrows, Figure 3C). During each restoration step, the same 3-D movement vector was used for each fault.

The optimum β angle is typically steep, antithetic to the main fault (Schultz-Ela, 1992; Withjack and Peterson, 1993; Hauge and Gray, 1996). To demonstrate the influence of the shear angle on the final restoration, a structural restoration with three different shear angles is shown for one fault block (Figure 9). The amount of strain within the hanging wall and the correctness of restoration depend mostly on the β angle (the precondition is of course a correct seismic interpretation). The highest priority is the resulting morphology of the surface, that is, the restored hanging-wall surface should be as smooth as possible and should match the original depositional surface (i.e., the regional). In this example, a shear angle of 55° represents the optimum result with a smooth surface of the hanging wall. This also inherently causes the smallest strain magnitudes within the hanging wall (Figure 9D). Because the dips of the normal faults vary in this model, antithetic shear angles between 45° and 80° were used. The shear angles found to work best for the faults in the study area are listed in Table 2.

Figure 9. Examples of oblique inclined shear restoration. The hanging wall (HW) is deformed by penetrative slip planes, using an inclined shear vector with the shear angle β. (A) Unrestored state, (B) β = 85°, (C) β = 75°, (D) β = 55°. The result obtained in (D) clearly shows a smoother top surface and the smallest strain values, where is the principal major strain axis and e is the strain value. β = 55° is the best option for restoration in this example. FW = footwall.

DEFORMATION QUANTIFICATION

Extension Over Time

The 2-D cross section (extracted from the 3-D model, Figure 10A) shows that the Skull Creek Mst. is highly faulted by southwest-dipping normal faults and northeast-dipping, antithetic faults (Ziesch et al., 2017). The first restoration step begins with the last movement, that is, the youngest fault zone. The study area is part of a passive margin rift basin, that is, the sediment depocenters and fault activity both migrated toward the sea-floor spreading center in the south (see also Finlayson et al., 1998). Therefore, it is highly probable that the order of faulting, even in this small area, was from north to south. The restoration of the fault sequence was carried out in the opposite sense. After complete restoration, the length of the retrodeformed Skull Creek Mst. is reduced and the entire restored formation is tilted northward by approximately 5° (Figure 10B).

Figure 10. The interpreted and retrodeformed Skull Creek Mudstone (Mst.) horizon. (A) The two-dimensional cross section of the deformed Skull Creek Mst.; (B) the same horizon restored and additionally decompacted to the preslip stage. The detailed view shows two restoration steps (vertical exaggeration = 1.5): (C) interpreted horizon with two normal faults, (D) restoration of the hanging wall of the younger, antithetic fault, and (E) restoration of the older synthetic fault. The dashed lines indicate the active objects during movement. Δ_l = increment of elongation; FP = fixed pin; LP = loose pin.

The elongation of the faulted horizons was determined by analyzing 10 2-D cross sections through the 3-D model (for location, see Figure 1), before and after restoration. The extension versus age diagram (Figure 11) illustrates the elongation of the stratigraphic horizons over time. The rate of extension decreased abruptly from top Waarre Fm. to top Narrawaturk Fm. Although the median extension of the Waarre Fm. is 9.2% and that of Skull Creek Mst. is 7.1%, the variation in the extension for the 10 cross sections is up to 12%, i.e., almost twice the median value. After Timboon Sandstone (Sst.) deposition (66 Ma), the total elongation after this point was only approximately 1%. Most extension occurred during the deposition of Skull Creek Mst. and the underlying Waarre Fm., that is, most faults terminated at approximately 66 m.y.a.

Figure 11. Extension versus age diagram showing the incremental elongation of the 10 cross sections (see locations in Figure 1). The box represents the 25%–75% percentile, and the whiskers represent the minimum and maximum extent. During Late Cretaceous (top Waarre Formation [Fm.] and top Skull Creek Mudstone [Mst.]), the cross sections show median values of 7%–10% extension, but there is a very wide range (see text for details). The deformation activity dies out rapidly after Timboon Sandstone (Sst.) deposition. In all 10 cross sections, 99% of extension occurred before deposition of top Timboon Sst. (66 Ma). After 66 Ma, the elongation was almost only 1%, that is, most faults had died out by this time. Gr. = Group; Mst. = Mudstone.

Strain Distribution

To better understand the finite retrodeformation strain maps, the different components of compactional strain and deformation strain for the Skull Creek Mst. were separately analyzed. The development of the compaction strain within the Skull Creek Mst. over time is presented in Figure 12. The backstripping of the overburden strata between the topography and the top Mepunga Fm. cause compaction strain magnitudes up to 10% within the Skull Creek Mst. (Figure 12A). The northwestern part of the study area displays slightly higher deformation than the remainder of the study area (Figure 12A). After the backstripping of the Mepunga Fm., the compaction strain increases up to 25% across the area (Figure 12B). Note that because of software limitations, extreme strain values that occur near the faults should not be considered meaningful (the software uses vertical decompaction and calculates therefore unnatural strain in the shadow of the faults; Midland Valley Exploration Ltd., 2016, personal communication). The next backstripping step caused an increase of the strain values by approximately 5% in large areas (Figure 12C). Following the backstripping of the Pebble Point Fm. and Timboon Sst., the compaction strain increases significantly. The southern area now has the highest compaction strain values (Figure 12D).

Figure 12. Cumulative strain caused by compaction within the Skull Creek Mudstone (Mst.) calculated by stepwise backstripping of the following overburden strata: (A) between topographical surface and top Mepunga Fm., (B) between top Mepunga and top Dilwyn Fm., (C) between top Dilwyn and top Pebble Point Fm., (D) between top Pebble Point Fm. and top Skull Creek Mst. (see Figure 2). See text for more information about strain artifacts. Grid coordinate system: Map Grid of Australia, zone 54, Transverse Mercator projection.

The deformation strain within the Skull Creek Mst. is shown at different restoration steps (Figure 13). The color-coded maximum strain magnitude increases southward and ranges from 0% to approximately 75% (Figure 13). A total of 16 restoration steps were performed (Figure 13D), which corresponds to the number of faults that were restored. Four faults are a major factor; they are as follows: fault 16, fault Naylor South, fault Buttress Northwest, and fault 2, because the strain is significantly unequally distributed between these structures. In general, the Skull Creek Mst. can be divided into two parts: the southern part underwent much more deformation than the northern area (Figure 13D). In Figure 13C and D, stripe features can be observed. These features are related to the fault corrugations (see Discussion section). They do not reflect the strain orientation; rather, they only reflect zones of high strain magnitudes.

Figure 13. Maps of the deformation strain ( strain magnitudes, where e is the strain value and is the principal major strain axis) after three-dimensional (3-D) structural restoration of the Skull Creek Mudstone (Mst.). The faults causing the tripartition are labeled. (A) Unrestored Skull Creek Mst.; fault heave polygons are shown in white. (B) First intermediate step up to fault 16 and fault Naylor South. (C) second intermediate step after restoration up to fault Buttress Northwest (NW) and fault 2. (D) complete 3-D structural restoration of the Skull Creek Mst. See text for more information about stripe artifacts. Dashed line shows profile location of section in Figure 16. Grid coordinate system: Map Grid of Australia, zone 54, Transverse Mercator projection.

Decompaction and restoration together determine the finite strain. The results after backstripping of the overburden and subsequent 3-D structural restoration of the Skull Creek Mst. are shown in Figure 14. The background strain is the compaction strain that occurred within the Skull Creek Mst. over the complete approximately 80 m.y. The northern area has strain values between 12% and 25%, whereas the central and southern part has slightly higher pre-Skull Creek Mst. deformation strain. After the restoring of all fault blocks, the strain is significantly inhomogeneously distributed (Figure 14, similar to Figure 13).

Figure 14. Maps of the strain magnitudes, where e is the strain value and is the principal major strain axis, after decompaction and subsequent three-dimensional (3-D) structural restoration of the Skull Creek Mudstone (Mst.). The strain magnitude is color-coded and represents the cumulative strain from both decompaction and structural restoration. The faults causing the tripartition are labeled. (A) Unrestored Skull Creek Mst. after decompaction and structural restoration of the overburden. (B) First intermediate step up to fault 16 and fault Naylor South. (C) Second intermediate step after retrodeformation up to fault Buttress Northwest (NW) and fault 2. (D) Complete 3-D retrodeformation of the Skull Creek Mst. grid coordinate system: Map Grid of Australia, zone 54, Transverse Mercator projection.

The combination of compaction and deformation strain offers a more detailed subdivision: the strain north of fault Buttress Northwest and fault 2 is, in comparison to the other areas, mostly from deposition of the overburden and not from tectonics. In contrast, the southern part is strongly deformed, mostly because of tectonic movement. The central area between fault Naylor South and fault 16 in the south and fault Buttress Northwest and fault 2 in the north underwent intermediate amounts of deformation.

DISCUSSION

From 3-D retrodeformation analysis, it is possible to infer paleogeomorphology, elongation rate, and subseismic scale strain that occurred within the faulted blocks of geologic formations. Comparative studies are rare (e.g., Mohr et al., 2005; Lohr et al., 2008; Freeman et al., 2015). Most case studies analyzed the structural style, fault growth, and fault population of seismic-scale faults (e.g., Needham et al., 1996; Mansfield and Cartwright, 2001; Nicol et al., 2005; Tvedt et al., 2013), but subseismic scale faults and fractures are commonly not considered. For the first time, both the seismically visible faults (Ziesch et al., 2017) and, in this paper, the prediction of the location of subseismic scale structures are evaluated in the context of a CCS project.

Paleogeomorphology

After decompaction and structural restoration of the 3-D geologic model, the complete Skull Creek Mst. is tilted northward (Figure 10B). Typically, a stratigraphic horizon is tilted by faulting, for instance, if the fault is nonplanar. Extensional movement on a listric fault can generate a roll-over anticline (Williams and Vann, 1987) or block rotation can be caused by domino-style faulting on a planar fault. Therefore, it is proposed that the study area was probably part of the hanging wall of a much larger normal fault (Figure 15). This major fault is located in the north of the study area and counts as one of the biggest faults of the onshore Otway Basin (for location, see Figure 1). However, thickening of any sedimentary units toward the northeast was not observed, in fact, the opposite is the case. The study area is too far from the major fault (∼20 km [∼12 mi]) to observe this feature. Activity on this major fault probably led to tilting of the Skull Creek Mst. The fault was only active prior to the end of Skull Creek Mst. deposition because, after this time, the retrodeformed beds are no longer tilted.

Figure 15. Schematic diagram showing an interpretation to explain the northward-tilt of the restored Skull Creek Mudstone. The area was probably affected by a major fault (see Figure 1 for location) outside the study area, to the north.

Figure 16. Two-dimensional cross sections of the retrodeformed Skull Creek Mudstone (Mst.) and Timboon Sandstone (Sst.) with their respective strain distributions and interpretations (for location of cross section, see Figure 13D). The strain magnitude, where e is the strain value and is the principal major strain axis, is color-coded. (A) Timboon Sst.: the strain magnitude is close to zero, except for the southernmost tip. (B) Sketch showing the amount of subseismic scale fracturing is minimal. (C) Skull Creek Mst.: the strain magnitudes are higher in the southwestern and central parts. (D) Sketch showing the relative distribution of the subseismic scale faults in the Skull Creek Mst. Areas with patches of higher strain magnitudes (violet). These areas are potential pathways for fluid or gas migration, but subseismic scale fractures of the underlain Skull Creek Mst. do not continue into the Timboon Sst. above.

Elongation Rate

The analysis of extension over time shows that the largest variation in the amount of total extension is between top Waarre and top Timboon (Figure 11). The average horizontal extension rate within the study area during the Late Cretaceous was 0.36% per m.y., which corresponds to approximately 0.028 mm/yr (∼0.001 in./yr). This is very small, compared with the extension rates of other passive margins (>3.5 mm/yr [0.1 in./yr]; Brune et al., 2014). The extension rate of the study area does not, however, reflect the extension rate of the whole Otway Basin, which would require the analysis of the complete passive margin.

Most deformation occurred during the Cretaceous and stopped at the beginning of the Paleocene. The abrupt stop of extension at 66 Ma, which is associated with the absence of faults, represents the transition from synrift to postrift in this area and is probably caused by the start of the production of true oceanic crust in the Southern Ocean (Blevin and Cathro, 2008). This marks the relaxation of the extensional stress field in the Otway Basin and the northern passive margin as a whole (see discussion in White et al. (2013)).

Subseismic Scale Strain

The distribution of strain within a 3-D geologic model is of special interest because this is important for the prediction of the location of potential leakage pathways (e.g., Lohr et al., 2008; Krawczyk and Tanner, 2010; Watkins et al., 2015, 2018). Walsh and Watterson (1992) estimated the number of subseismic scale faults using a fractal method. However, this method assumes that subseismic scale fractures exist homogeneously throughout a volume. Other authors have attempted to predict fracture density from outcrop or seismic surfaces by using curvature or the azimuthal amplitude versus offset method, albeit sometimes with conflicting results (Lisle, 1994; Hunt et al., 2010; Suo et al., 2012). The approach taken in this paper is to kinematically restore each fault and thereby assess the strain that occurs within each hanging-wall volume. We postulate that the magnitude and distribution of the strain in the hanging wall after restoration is a proxy for the amount and distribution of subseismic scale deformation (Plesch et al., 2007; Lohr et al., 2008). This deformation, above a certain magnitude, is probably manifested as subseismic scale fracturing. According to van der Pluijm and Marshak (2004), brittle fracturing occurs at finite strains of approximately 5%. However, this model could not be calibrated with, for example, formation microimager (FMI) data from boreholes (as carried out by, for instance, Lohr et al., 2008), and therefore the amount of strain calculated in the 3-D retrodeformation is treated here in a relative manner, that is, we postulate that high amounts of modeled deformation correlate with a large density of subseismic scale fracturing. The amount of subseismic scale fracturing can be evaluated by interrogating the strain magnitude of the tetrahedra that make up the hanging wall (e.g., Lohr et al., 2008; Krawczyk et al., 2015). The finite strain tensor within a stratigraphic volume is composed of strain caused by differential compaction and strain caused by tectonic deformation.

The calculated, uniaxial, compaction-related deformation is a proxy for the amount of volume loss (Isaacs, 1981). The results show that the presently 2.2-km (7217-ft)-deep Skull Creek Mst. has minimum compaction strain values of 15%. The southern part is slightly more compacted (20%–22%) than the northeastern area (Figure 12). This could have resulted in higher compaction-related flattening of the fault surfaces, which are more listric in the south. This is a similar observation to that of the theoretical discussion in Jones and Addis (1984).

The strain results of the structural restoration within the Skull Creek Mst. show, at 30%–50%, higher values in the southwest than in the northeast (Figure 13D). The differences in strain along fault strike are probably caused by fault morphology. In addition, the fault surfaces are irregularly corrugated. The corrugations are subparallel to the dip direction of the faults, and they have widths of 100–500 m (328–1640 ft) and lengths of 500–2750 m (1640–9022 ft), that is, above seismic resolution (Ziesch et al., 2017). They are mostly not continuous along the length of the faults, although the breaks are not related to a particular strata. Therefore, when restoring, even parallel to the corrugation direction of the faults, the irregularity of the fault surface reveals internal strain in the hanging wall. The amount of subseismic scale deformation is postulated to be dependent on the amplitude of the corrugations because regions with higher strain correlate with more corrugations of the fault surfaces (Figure 13B–D). All geologic formations underwent more deformation in areas where the fault corrugations are more prominent. Similar to these results, Lohr et al. (2008) observed in the Northwest German Basin that the deformation of the hanging wall is strongly controlled by the fault morphology. The highest strain areas were connected to strong changes in fault strike or fault dip.

The Appendix presents forward modeling of a hanging-wall volume on fault Buttress Southeast (Figures 17, 18). The forward and retrodeformation results are comparable because the strain tensor is mirrored (Figure 8B). Even by reducing the resolution of the fault grid (Figure 18A), the strain distribution is still comparable. The strain maximums and the strain minimums are located in the same areas. Therefore, we interpret that the fault possesses a basic morphology, with topographic variation above 300 m (984 ft), that causes a basic strain pattern. The corrugations below 300 m (984 ft) produce additional strain that is superimposed on this basic pattern. The distribution of strain does not give information about the orientation of subseismic scale fractures. Figure 18B shows the calculated strike directions of subseismic scale fractures based on the orientation of the strain tensor plane). Note that lineations of high strain magnitude are not parallel to the fracture orientations.

Figure 17. (A) Map showing the location of fault Buttress Southeast (SE). Black box shows the area modeled in Figure 18. (B) Surface of fault Buttress SE; color-coded by curvature (above) and cylindricity (below); insert shows stereographic projection of poles to the surface triangles. (C) Two triangulated fault surfaces of Buttress SE. The lighter blue surface represents the 75-m (246-ft) gridded surface used in this work. The darker green surface represents the 150-m (492-ft) gridded surface. Grid coordinate system: Map Grid of Australia, zone 54, Transverse Mercator projection.

Figure 18. (A) Forward modeling of a planar volume (75 m [246 ft] mesh) along the fault Buttress Southeast at different amounts of movement (100 m [328 ft], 300 m [984 ft], 500 m [1640 ft], and 700 m [2297 ft]). The strain magnitude, where e is the strain value and is the principal major strain axis, is color-coded and represents the cumulative strain. The left column represents the 150-m (492-ft) fault grid and the right column represents the 300-m (984-ft) fault grid. (B) Example for the strain map and the resultant strike direction (black lines) of subseismic scale fractures.

Because of the complex fault surface morphologies, the best shear angle value for the retrodeformation of each fault was determined separately (Figure 9; Appendix, Table 2). The appropriate shear angle for inclined shear has been the subject of discussion since the early 1980s. Although Groshong and Richard (1989) recommend shear angles that are exactly antithetic to a major fault, Xiao and Suppe (1992) suggested shear orientations for normal faults of 67°, which approximate the Coulomb collapse of deforming materials. The shear angles used in this work fit very well with the suggestion of Xiao and Suppe (1992), especially considering that lateral variation of the fault shape can change the shear angle by up to ±15° (Yamada and McClay, 2003).

The strain occurs as stripe features during restoration. It can be excluded that these are caused by a numerical error because of the mesh density. All surfaces and volumes were manually constructed and checked for integrity (e.g., bad triangles in surfaces and acute angles in tetrahedra; see the quality of the triangulation in the 3-D supplementary PDF in Ziesch et al., 2017). The cause of the strain magnitude stripes are the fault corrugations and their impact on the inclined shear method. Other authors, for instance Lohr et al. (2008) and Watkins et al. (2015, 2018), have shown similar features occur during 3-D restoration using the same software (Move, Midland Valley Exploration Ltd., 2016). However, as Brandes and Tanner (2014) remark, all kinematic models of fault restoration in the present-day are limited because of the use of the same kinematic behavior for every part of the fault (instead of a combination of different algorithms) and the still rudimentary application of originally 2-D algorithms in the third dimension.

After the combination of decompaction and subsequent 3-D restoration of the Skull Creek Mst., the study area has a strongly deformed (55%–80%) southern part, a medium deformed (25%–55%) central area between the four major faults, and the northern area is less deformed (10%–25%; Figure 14D). The tripartition of the study area supports the hypothesis of a hierarchical system of normal faults. Although 24 normal and antithetic faults were defined as major faults (based on fault length and throw), these 4 faults acted as higher-order master faults (see also Needham et al., 1996; Watterson et al., 1996). They are the biggest faults in the study area; they have the largest lengths (in the strike direction) and offsets up to 800 m (2625 ft).

It is noticeable that, particularly in the areas of relays (Figure 3C, F), more subseismic scale deformation occurred (Figure 13A; e.g., fault 10, fault Naylor South). Ziesch et al. (2017) observed torsion of the relays at the depth of the Skull Creek Mst. By retrodeforming the complete model, including the relays, the amount of additional deformation that the relays underwent is made apparent (see Figure 14D).

The detailed analysis of subseismic scale deformation within the complete overburden provides insight into potential pathways. These could significantly influence the flow of CO₂, hydrocarbons, or hydrothermal fluids. To better illustrate this, Figure 16 shows a 2-D section (for location, see Figure 13D) through the decompacted and restored Skull Creek Mst. and Timboon Fm. The strain magnitude within the stratigraphic formations is used as a proxy to evaluate the location of potential migration pathways (Lohr et al., 2008). Strain by compaction is uniaxial and therefore not taken into account. We propose that parts of the study area with strain magnitudes over 2%–3% are potential areas of subseismic scale faulting, and therefore point to pathways of potential fluid or gas migration.

The strain magnitude within the overlying Timboon Fm. is close to zero, which means that within this and also all other overburden formations, the amount of subseismic scale fracturing is minimal. Even if the model was increased to the north, using available data, which is not presented here, this effect would still be observed. Therefore, subseismic scale fracturing does not continue beyond the seal formation into higher stratigraphic formations. Within the higher stratigraphy, only seismically visible faults could be pathways. Thus, the retrodeformation method is able to locate possible critical areas of subseismic scale fracturing (Figure 16). Validation as to whether these fractures and therefore pathways are sealing or not could not be accomplished. However, in a case study of an inverted fault in the Northwest Germany Basin, Lohr et al. (2008) showed that strain shown by 3-D retrodeformation correlates well with the results of FMI data from boreholes within the hanging wall, that is, there was a positive correlation between the amount of strain and the amount of fractures determined by FMI.

The subseismic scale deformation in the Otway Basin occurred coevally with the Cretaceous sedimentation. Isochore maps showed that nearly all faults are growth faults at different stratigraphic levels, that is, strong fault growth occurred during movement on the faults. Therefore, shortly after deposition, the partially consolidated Skull Creek sedimentary rocks were deformed. Because the rocks were still close to the surface, the sediment φ was approximately 40%–45%. It is likely that, because of increased fluid flow at this time, mineral precipitation occurred within the fractures.

SUMMARY AND CONCLUSIONS

For the first time in the Otway Basin, a 3-D structural restoration and strain analysis of a geologic model of a storage reservoir was carried out. This paper shows the following.

• The compaction-related deformation increases with depth up to approximately 30% volume loss.

• The study area is controlled by four major faults, where the highest strain values occur in the southwestern part of the study area. These faults are of a higher hierarchy, that is, they acted as higher-order master faults, with respect to the other 20 major faults.

• Subseismic scale fracturing of the seal horizon (Skull Creek Mst.) occurs and is inhomogeneously distributed.

• Subseismic fracturing does not extend into the overlaying formations.

• Restoration yields a unique understanding of the geologic evolution of the area. The Skull Creek Mst. was northwardly inclined during its deposition because of movement on a major fault outside the study area.

It is concluded that it is possible to retrodeform a complete 3-D model using the workflow presented in this paper to identify the distribution and magnitude of subseismic scale fractures. Strain maps help to predict the relative distribution of fracture density in areas without well data, and thus to predict potential migration pathways. The strain results remain relative until calibration is possible using FMI data or core data.

We recommend that the 3-D retrodeformation method should always be used on baseline data before any subsurface use to determine whether subseismic pathways exist, which are prime structures along which leakage out of a reservoir could occur.

This appendix deals with the theoretical consideration as to what degree fault corrugations influence strain in the hanging wall. Forward modeling is presented for fault Buttress Southeast, which is located in the eastern part of the study area (Figure 17A). In the seismic interpretation, the fault surface is characterized by corrugations, which are strong, parallel linear zones of high curvature, see Figure 17B. Ziesch et al. (2015) demonstrated that these corrugations are above seismic resolution. To determine the influence of the corrugations on the strain distribution within the hanging wall, the original triangulated fault surface was resampled (direct triangulation, mean approximately 50 m [∼164 ft]) at two coarser regular grid sizes (150 m [492 ft] and 300 m [984 ft]), and forward modeling was carried out on the regular-gridded hanging-wall volume on these fault surfaces. During the forward modeling, the hanging-wall volume was moved stepwise downward every 100 m (328 ft), for a total of 700 m (2297 ft) along the fault surface. The strain distribution in the hanging wall is shown after 100 m (328 ft), 300 m (984 ft), 500 m (1640 ft), and 700 m (2297 ft) (Figure 18A). Figure 18B shows the orientation of the strain tensor after 700-m (2297-ft) movement of the hanging wall on the 150-m (492-ft) gridded fault surface.