FAULTS AND SEDIMENTS

FIELD OF THE INVENTION

[0001 ] The present application relates to geological modelling and in particular to methods of determining permeability distributions in physical media using geological models.

BACKGROUND

[0002] The ability to predict fluid flow pathways in sedimentary basins or any other similar geological zones is relevant in many areas of geoscience, including mineral exploration, petroleum exploration and production, assessment of geothermal and groundwater resources, contaminant dispersal, nuclear waste disposal, and CO2 sequestration. In all of the above examples, the distribution of permeability in the basin fill and faults is a key factor controlling fluid flow.

[0003] Permeability can vary by orders of magnitude between different rock/sediment types, and between faults and their host rocks. It can also be highly anisotropic, with the direction of anisotropy typically following stratigraphic and structural features.

[0004] Fault permeability in sedimentary basins is particularly complex, with faults acting as fluid pathways, barriers or combined conduit-barrier systems depending on the nature of the host rocks and their deformation history.

[0005] Consequently, a detailed knowledge of the subsurface geology, along with simulation tools capable of representing the heterogeneous and anisotropic permeability distribution in faults and their host rocks/sediments, are required in order to simulate fluid flow in sedimentary basins.

[0006] Presently, with the targeting of deeper mineral resources, the mineral exploration industry is shifting towards the development of advanced 3D geological models adequately representing both the structural framework (Le., faults and fractures) and the distribution of rock properties. However, one of the drawbacks is that the integration of faults and their damage zones in these models rarely captures the corresponding permeability features, while the distribution of rock properties is commonly estimated using geostatistical approaches which are unreliable in environments where data are sparse and spatially clustered. [0007] The distribution of rock properties in 3D models is commonly estimated using geostatistical methods, however such methods are unreliable where outcrops are lacking and drillhole data are sparse or spatially clustered, as is often the case in mineral exploration. One of the major challenges in fluid flow simulations is handling the complexity of natural fault systems. Examples of such complexities include wavy fault traces, overprinting relationships of different fault generations, and, resulting from this, non-trivial fault intersection geometries.

[0008] Modelling such systems is especially demanding in 3D, not only due to the challenge of building the geometry and associated mesh, but also due to the computational cost of the resulting simulations.

[0009] Volumetric meshing of fault zones, where the small but finite width of the fault is discretised into multiple mesh elements results in numerous very small elements, results in a large computational cost due to the large number of elements and the small timestep that must be used to capture transient fluid flow through such a fine mesh. In addition, it is challenging to generate appropriate meshes for systems that include numerous strongly curved and intersecting fault zones when discretising volumetrically.

[0010] The present applicant has recognized the above-mentioned or related shortcomings and the need to provide an improved method for creating a model with anisotropic permeability in both faults and host rocks.

[001 1 ] Any discussion of the background art or a patent document or other matter throughout the specification should in no way be considered as an admission that such art is widely known or forms part of common general knowledge in the field.

SUMMARY OF THE INVENTION

[0012] The present invention comprises a method for determining permeability distributions in a geological region as defined in the attached claims.

[0013] In one aspect, we describe a computer-implemented method for determining permeability distributions in a geological region, the method comprising: providing a first mesh which comprises at least one cell; wherein each cell of the at least one cell comprises stratigraphic information of the geological region; wherein the stratigraphic information comprises porosity information; building a second mesh; resampling and homogenizing the stratigraphic information from the first mesh to derive anisotropic permeability in the second mesh; building a third mesh, wherein the third mesh comprises at least one geological discontinuity; setting anisotropic permeability in the third mesh from the derived anisotropic permeability in the second mesh by interpolating from the second mesh to the third mesh; determining permeability distributions in the geological region based on the third mesh with set anisotropic permeability.

[0014] In one example, the method further comprises providing the first mesh by stratigraphic forward modelling.

[0015] In one example, the method further comprises using the third mesh with set anisotropic permeability in a fluid flow simulation.

[0016] In one example, the stratigraphic information further comprises information about at least one sediment class.

[0017] In one example, the at least one sediment class is sand and/or shale and/or silt and/or carbonate and/or mud and/or clay and/or any other sediment type.

[0018] In one example, wherein the second mesh has the same horizontal resolution but different vertical resolution than the first mesh.

[0019] In one example, the resampling and homogenizing of the stratigraphic information from the first mesh to derive anisotropic permeability in the second mesh comprises: computing values of longitudinal and transverse permeability from stratigraphic information present in the each cell of the at least one cell in the first mesh; homogenizing and resampling stratigraphic information onto the second mesh by calculating the weighted arithmetic mean of stratigraphic information in each cell of the at least one cell of the first mesh that intersects each cell of the second mesh, wherein weights in the weighted arithmetic mean are proportional to the average height of each cell in the first mesh; homogenizing and resampling longitudinal permeability onto the second mesh by calculating the weighted arithmetic mean of longitudinal permeability in each cell of the at least one cell of the first mesh that intersects each cell of the second mesh, wherein weights in the weighted arithmetic mean are proportional to the average height of each cell in the first mesh; homogenizing and resampling transverse permeability onto the second mesh by calculating the weighted harmonic mean of transverse permeability in each cell of the at least one cell of the first mesh that intersects each cell of the second mesh, wherein weights in the weighted harmonic mean are proportional to the average height of each cell in the first mesh; calculating the average orientation of each cell of the at least one cell of the first mesh that intersects each cell of the second mesh; obtaining a homogenised permeability tensor for each cell in the second mesh by rotation of a diagonal tensor formed by the homogenised and resampled longitudinal permeability and transverse permeability and their average orientation into the global coordinate system.

[0020] In one example, building the third mesh comprises using traces of the at least one geological discontinuity and associated meta data to generate the at least one geological discontinuity in the third mesh.

[0021 ] In one example, the third mesh is a two-dimensional (2D) or a three-dimensional (3D) mesh.

[0022] In one example, the at least one geological discontinuity is at least one geological fault and/or at least one geological fracture.

[0023] In one example, the associated meta data is thickness and/or dip and/or dip direction for each geological discontinuity of the at least one geological discontinuity.

[0024] In one example, the at least one geological discontinuity is meshed as volumetric features or lower-dimensional features.

[0025] In one example, the at least one geological discontinuity and a surrounding volume of the at least one geological discontinuity are discretised into an unstructured mesh.

[0026] In one example, the resampled and homogenised stratigraphic information in the second mesh is interpolated onto the third mesh.

[0027] In one example, the longitudinal permeability, transverse permeability and the homogenised permeability tensor on the second mesh are interpolated onto the third mesh outside of the at least one geological discontinuity.

[0028] In one example, the permeability of the at least one geological discontinuity in the third mesh is derived from the longitudinal permeability and/or the transverse permeability and/or the at least one homogenised permeability tensor and/or at least one sediment class fraction in the adjacent mesh cells of the third mesh using a user-defined relationship.

[0029] In one example, the method further comprises: outputting the determined permeability distributions in the geological region to a processor configured to perform various functions and/or control hardware.

[0030] In one example, the method further comprising: inputting the third mesh with set anisotropic permeability to a fluid flow simulator to determine permeability distributions in the geological region.

[0031 ] In another aspect, we describe a method of determining permeability distributions in a geological region, the method including: receiving stratigraphic data of the geological region including fault locations; receiving sample material properties of the geological region; constraining a stratigraphic forward model of the geological region based on the stratigraphic data and material properties to produce a four-dimensional distribution of porosity and sediment class fractions of the geological region; and estimating transverse and longitudinal permeability from the four-dimensional distribution to generate a permeability distribution of the geological region.

[0032] In one example, the method including the step: inputting the permeability distribution to a fluid flow simulation to estimate fluid flow within the geological region.

[0033] We also describe a computer-readable storage medium comprising instructions which, when executed by a computer, cause the computer to carry out the method of any one of the preceding claims.

[0034] We also describe a computer program comprising instructions which, when the program is executed by a computer, cause the computer to carry out the method of any one of the preceding claims.

BRIEF DESCRIPTION OF THE FIGURES

[0035] Example embodiments of the disclosure will now be described, by way of example only, with reference to the accompanying drawings in which:

[0036] Figure 1 A is a flowchart of the method of the present techniques; [0037] Figure 1 B is an overview diagram of the method of the present techniques;

[0038] Figure 2 is a 2D schematic representation of the first mesh and the second mesh;

[0039] Figure 3 is a 2D example of permeability homogenisation, with the first mesh on the left hand side and the second mesh on the right hand side;

[0040] Figure 4 illustrates discretisation methods for faults. Panel (A) illustrates a volumetric approach where the fault thickness is meshed explicitly with small elements (shaded in grey). Panel (B) illustrates single lower dimensional element (conduit). Panel (C) illustrates an interface created by breaking the mesh along the fault trace and superposing duplicated nodes (barrier). Panel (D) illustrates a combination of B and C where lower dimensional elements are located on either side of the interface (conduit barrier);

[0041 ] Figure 5 illustrates a hypothetical model with two intersecting faults, one with higher permeability and one with lower permeability than the host rock. Fluid flow is driven from left to right by a fluid pressure gradient. Panels (A) and (B) show fluid pressure contours obtained using the volumetric and lower-dimensional meshing approaches, respectively;

[0042] Figure 6 is a 3D example of the third mesh generation from geological discontinuities defined by their traces at the surface and their meta-data;

[0043] Figure 7 is a schematic diagram of a computing device suitable for implementing the method of the present techniques;

[0044] Figure 8A schematically illustrates the location of a case study area in the southern McArthur Basin, showing extent of the stratigraphic forward model (red box/ 802) and fluid flow models (black box/804). Structural framework and sedimentary package distribution are from Blaikie and Kunzmann (2020). “Modelled faults” are those that were included in the SFM and fluid flow model. Stars indicate locations of the Teena (T) and McArthur River (M) Zn-Pb mineral deposits.

[0045] Figure 8B schematically illustrates a lithostratigraphic column for the study area with 3 ^rd order sequence stratigraphic interpretation (Kunzmann et al., 2022). MFS: Maximum Flooding Surface, MRS: Maximum Regressive Surface;

[0046] Figure 9 schematically illustrates a process from well correlation to stratigraphic model. Panel (A) illustrates facies correlation between wells in the Glyde sub-basin. Panel (B) illustrates paleogeographic and thickness maps for the study area. Panel (C) illustrates shale distribution in the Glyde sub-basin area from stratigraphic forward modelling. MFS: Maximum flooding surface, MRS: Maximum regressive surface;

[0047] Figure 10 illustrates a section through the stratigraphic forward model in the Glyde sub-basin area showing: in (A) the original grid (i.e the first mesh) and in (B) the resampled grid (i.e the second mesh). Vertical exaggeration x 5;

[0048] Figure 11 illustrates porosity and sediment class fractions on cross-sections through the stratigraphic forward model in the Glyde sub-basin area. Results are shown at the end of the simulated period (1624 Ma). Vertical exaggeration x 5;

[0049] Figure 12 illustrates porosity on an E-W section through the West side of the Glyde sub-basin in (A) the original stratigraphic forward model (i.e. the first mesh) and in (B) the resampled stratigraphic forward model (i.e. the second mesh). Vertical exaggeration x 5;

[0050] Figure 13 illustrates longitudinal (left) and transverse (right) permeability (m ² ) on cross-sections through the resampled stratigraphic forward model (i.e. the second mesh) in the Glyde sub-basin area. Vertical exaggeration x 5;

[0051 ] Figure 14 illustrates an exemplary third mesh with set anisotropic permeability used in fluid flow simulations. Panel (A) illustrates an entire mesh. Panel (B) illustrates a detail of the mesh around a fault; and

[0052] Figure 15 illustrates fluid flow simulation results for 4 permeability scenarios. Left: Permeability and fluid flow vectors in faults. Middle: Fluid flow vectors and shale content on cross-sections and shale content in faults. Right: Tracer distribution on cross-sections at end of simulation. Vertical exaggeration x 2 in all plots.

DESCRIPTION OF THE INVENTION

System overview

[0053] Described herein is a method that allows building a model which captures the anisotropic permeability distribution from a stratigraphic forward model (SFM) and structural framework, in a mesh (i.e. third mesh with anisotropic permeability) which is suitable to run fluid flow simulations. Thus, beneficially, the method can predict fluid pathways through a geographical/geological zone under various fault permeability scenarios. Another advantage of the present method is that the subtleties of the fluid pathways predicted by the method of the present techniques cannot be easily predicted by models that assign uniform properties to faults and their host rocks.

[0054] The features which are the same or substantially similar retain the same reference numbers in the Figures and throughout the description. Therefore, a repeated description of such features is omitted.

[0055] All examples of the present techniques may be combined together, unless it is explicitly stated otherwise.

[0056] Described herein is a method to generate 3D models of faults/fractures and their host rocks or sediments, using information from process-based models of sediment deposition as opposed to traditional geostatistical methods to populate the rock property values in the model. By using stratigraphic information from process-based models, such as the stratigraphic architecture and rock property distribution, a homogenised anisotropic permeability field is derived based on user-provided porosity-permeability relationships and the orientation of stratigraphic layers. An unstructured 3D grid is then built which accounts both for the structural framework and the intervening host rocks. The porosity and anisotropic permeability fields generated from the process-based model are used to populate the porosity and permeability values in the unstructured 3D grid, both within and outside the faults. Outside the faults, the anisotropic permeability is set directly from the process-based model. Permeability within the faults is set as a function of the corresponding host rock permeability and/or sediment type, depending on the nature of the discontinuity. The method of the present techniques may use a mixed dimensional approach for the discretization of a fault region, whereby faults are represented by elements of a lower dimension: 1 D elements (lines) in 2D models and 2D elements (surfaces) in 3D models. In other words, the method of the present techniques may discretise faults using volumetric elements of the same dimension as the overall mesh (e.g. 3D elements in a 3D mesh) and/or lower-dimensional elements (e.g. 2D elements in a 3D mesh). The method of the present techniques can be applied in simulations of flow in discrete fracture networks and their host rocks, where the fractures are represented as zero-thickness features in conforming meshes. However, the method of the present techniques can also be utilised to represent fault zones if the scale of the model is large enough that the fault zone thickness is negligible compared to the overall size of the domain (e.g., entire sedimentary basins or crustal- scale models). [0057] Figures 1 A and 1 B show an overview of the method according to present techniques. Figure 1 B shows an exemplary diagram overview of the method of the present techniques shown in Figure 1 A. This example shows how the combined elements generate a final 3D model with anisotropic permeability suitable for use in fluid flow simulations, starting from the structural framework (fault traces, thickness, dip, dip direction), geological information of the zone and user-selected porosity-permeability relationships derived from sedimentary rock properties.

[0058] Returning to Figure 1 A, in S100, the method comprises providing a first mesh. Generally speaking, a mesh may be defined as a geometrical component which represents a geometry of a target domain by smaller discrete cells, where the corner points of the discrete cells are called nodes. Accordingly, a mesh partitions space into cells, the corners of which are nodes. Therefore, a mesh is made of at least one cell. A cell may be defined as a control volume into which the target domain is broken up. Equations can be solved per each cell or node, which then approximates the solution to the equations over the target domain. The term “mesh” and “grid” may also be used interchangeably. The term “mesh cell” and “mesh element” can be used interchangeably. The terms: initial mesh, outputted mesh, SFM output and the first mesh may be used interchangeably. The initial mesh may be produced by Stratigraphic Forward Modelling. Thus, the first mesh may be a Stratigraphic Forward Model (SFM) grid.

[0059] Stratigraphic Forward Modelling is the use of mathematical formulae and algorithms to simulate sedimentary processes in order to create synthetic stratigraphy, with the aim of understanding and predicting dynamic sedimentary systems. The input of the stratigraphic forward modelling may be the geological information of the target zone. Geological information of the zone may comprise at least one of the following: age constraints on the studied stratigraphic interval, lithological descriptions from boreholes and/or outcrops, paleoenvironmental interpretations, stratigraphic interpretation and associated thickness and distribution of the sedimentary units in the study zone, eustatic variations, expert knowledge of climatic evolution during the simulated geological time affecting rainfall, evaporation and erosion over the study area, expert knowledge of the main structural and tectonic events and their impact on uplift and subsidence over the study area, and any other similar information.

[0060] The initial mesh comprises a sequence of vertically stacked layers with uniform grid spacing in the horizontal directions (X and Y axes), but varying spacing in the vertical direction (Z axis), with some layers potentially having zero thickness in parts of the domain. The first mesh contains at each cell stratigraphic information including the porosity and, optionally, proportions of one or more sediment classes.

[0061 ] In the absence of any information about sediment classes in the first mesh, it can be assumed that the entire first mesh represents a single sediment class; alternatively, sediment classes may be defined based on other information contained in each cell, such as porosity and/or grain size. The sediment classes may be for instance: sand and/or shale and/or silt and/or carbonate and/or mud and/or clay. However, the skilled person will appreciate that other sediment classes, depending on the geographical zone of interest may be used. Information regarding the sediment classes may also comprise information regarding each sediment class fraction. Sediment class fractions represent the volumetric proportions of each sediment class (numbers adding to 100%) in a cell.

[0062] Stratigraphic information comprises outputs from stratigraphic forward modelling, which include the porosity and optionally one or more sediment class fractions together with the thickness and orientation of all cells in the first mesh. The outputs from stratigraphic forward modelling may also include environmental and hydrodynamic properties (e.g. water discharge, wave energy). The skilled person will appreciate that any or all of these properties could be used to constrain permeability in the first mesh. In other words, stratigraphic information comprises porosity information. Stratigraphic information may further comprise information about at least one sediment class. The information about at least one sediment class may comprise information about sediment class fractions. The information about at least one sediment class may further comprise one or more sediment class fractions together with the thickness and orientation of all cells in the SFM output. In other words, stratigraphic information may further comprise one or more sediment class fractions together with the thickness and orientation of all cells in the SFM output. The stratigraphic information may further comprise environmental and/or hydrodynamic properties.

[0063] In an example, the stratigraphic forward modelling output may be in the form of a structured mesh. The outputted mesh (i.e. the first mesh/the initial mesh) may comprise stacked layers of 1 x 1 km cells with varying thicknesses. The skilled person will appreciate that the stacked layer cells may have different dimensions. These dimensions may be based on user input and/or the goals for a given task. In other words, the horizontal dimensions of the cells in each layer of the initial mesh may be chosen based on user input and/or the goals for a given task. Each layer in the outputted mesh may represent the net result of sediment deposition, erosion and compaction over a given time interval. The thickness of the layers varies over the target domain, and can be zero if no sediment was deposited or all sediment was eroded in a given cell. Each cell of the first mesh contains stratigraphic information, which includes the porosity and optionally also the proportion of each sediment class.

[0064] Therefore, Stratigraphic Forward Modelling may reconstruct the stratigraphic architecture of a geographical or a geological zone by simulating erosion and/or transportation and/or deposition of sediments in response to basement subsidence/uplift and/or sea level fluctuations and variations in sediment supply. A geographical/geological zone may be a basin and/or a rock reservoir. However, the skilled person will appreciate that a geographical zone may be any geographical zone in which is it beneficial to explore fluid flow. The terms geographical zone, geological zone, target/ modelled domain and geological region may be used interchangeably.

[0065] In one example, which is further described in the “Exemplary case study” section below, a structural framework comprising the location and orientation of major structural elements (e.g. faults) in a geological zone may be used to constrain a Stratigraphic Forward Model of that geological zone. The skilled person will appreciate that information about the location and orientation of major structural elements (e.g. faults) could be obtained in various forms, such as geological maps and 3D models. By way of example, the DionisosFlow software, developed by Beicip-Franlab, may be used to perform stratigraphic forward modelling simulation and thus provide the first mesh. The skilled person will appreciate that any other software package that is similar to DionisosFlow may be used to provide the first mesh.

[0066] Information about faults may be used to constrain subsidence at the base of the stratigraphic forward model, however, the faults are not represented explicitly within the stratigraphic forward model. Furthermore, the relatively coarse grid size of the stratigraphic forward model, which may span at least tens of meters, means that displacement across faults is not well-resolved.

[0067] The skilled person will appreciate that the stratigraphic forward model (SFM) may be calibrated against data from the area represented by the model, such as observed thicknesses and interpreted paleobathymetries of the main stratigraphic units. [0068] Therefore, the SFM may produce a 4D distribution of porosity and sediment class fractions over the modelled domain, comprising a series of meshes and corresponding porosity and sediment class distributions representing user-specified output times in the SFM simulation period. This provides enough detail to estimate the permeability distribution in the modelled area at any of the user-specified output times. Using the SFM results corresponding to the desired output time, longitudinal permeability (parallel to the depositional layers within the SFM) may be estimated from porosity using a Kozeny-Carman type relationship parameterised for each sediment class. Transverse permeability (perpendicular to the depositional layers within the SFM) may be assumed to be a fraction of longitudinal permeability, with the multiplier varying between sediment classes. The skilled person will appreciate that setting the multiplier equal to one will result in isotropic permeability.

[0069] Figure 2 is a 2D schematic representation of the first mesh 200 and the second mesh 210. In practice, the first mesh is usually 3D or 4D, however, for ease of explanation, the first mesh is shown as a 2D mesh in Figure 2. The skilled person will appreciate that 4D representation of the first mesh may comprise a series of 3D meshes representing points in time. Figure 2 illustrates the first mesh 200 with regular horizontal meshing and irregular vertical meshing. Successive layers 202, 204,206 (represented in blue 202, orange 204 and green 206) have varying thicknesses, zero in some places (represented by dashed lines).

[0070] Figure 2 also illustrates the second mesh 210, represented by black lines is created with the same horizontal discretisation, but different vertical spacing so that all six cells of the second mesh 210, in this non-limiting example, have a non-zero thickness. This may be done by selecting a series of geological units to maintain (in this example the whole model), with their top and bottom horizons being preserved and new mesh cells being created in between, regularly distributed vertically between the top and bottom surfaces from a user-defined number of layers and/or cells. The skilled user will appreciate that the cells of the second mesh could have non-uniform vertical spacing. In this example, two horizontal layers are shown. However, the skilled person will appreciate that the user may define more than two layers.

[0071 ] It is not possible to import information from the SFM output (i.e. the first mesh) directly into some fluid flow simulators (including finite element codes) due to the existence of zerothickness cells. Moreover, the fine vertical resolution of the SFM is likely to create problems with numerical stability and convergence in the fluid flow simulations. Beneficially, the method of the present techniques overcomes these problems by resampling and homogenizing the SFM output (i.e. the first mesh) onto a new structured mesh (i.e. a second mesh) with the same horizontal resolution but coarser vertical resolution as described in S102 of the method of the present techniques. Fluid flow simulations may use either the second mesh with resampled stratigraphic information and permeability S102 if geological discontinuities are not required, or the third mesh with set anisotropic permeability S106.

[0072] Returning to Figure 1 A, the method of the present techniques comprises building a second mesh S101.

[0073] The step S101 building a second mesh comprises creating a second mesh by specifying a vertical distribution of non-intersecting laterally continuous surfaces that match the horizontal extent of the SFM output (i.e. the first mesh) and span all or part of the vertical extent of the SFM output (i.e. the first mesh) and, optionally, the regions above and/or below it, and using these surfaces to define grid cells with the same horizontal discretization as the SFM output.

[0074] In other words, the second mesh has the same horizontal resolution but different vertical resolution than the first mesh. The second mesh has the same horizontal resolution as the first mesh, but the second mesh has a vertical resolution that is suitable for fluid flow simulations. The vertical resolution of the second mesh is usually within an order of magnitude of the horizontal resolution, which oftentimes results in coarser vertical resolution of the second mesh when compared with the first mesh.

[0075] The second mesh may keep the same resolution as the first mesh along X and Y axes, but the thickness of cells along the Z axis in the second mesh may be larger than that of the first mesh. Thus, each cell of the second mesh may cover multiple cells of the first mesh and result in smaller horizontal resolution: thickness ratio (aspect ratio) of the cells, for instance the aspect ratio of cells in the second mesh may be 5:1 instead of original ratio in the first mesh 100000:1 . Thus, the new thickness of cells along the Z axis is coarser than the original thickness of cells along the Z axis. Accordingly, the second mesh may have the same horizontal resolution but coarser vertical resolution than the first mesh. The coarser vertical resolution of the second mesh means that the second mesh has smaller horizontal resolution: thickness ratio (aspect ratio) than the first mesh. For some types of fluid flow simulations (e.g. simulations performed using a finite element code), it is desirable to maintain the aspect ratio of mesh cells as close as possible to 1 :1 . For example, the aspect ratio around 5 is suitable. However, the cells of the second mesh may have a larger aspect ratio if the second mesh exists only to provide input to the third mesh; that is, if the second mesh will not be used in fluid flow simulations.

[0076] In S102 of the method of the present techniques, the method comprises resampling and homogenizing the first mesh to derive anisotropic permeability in the second mesh. Deriving anisotropic permeability tensors in all elements of a second mesh is achieved by resampling and homogenizing the SFM output (i.e. the first mesh) onto a second structured mesh (i.e. a second mesh) with the same horizontal resolution but with a vertical resolution adapted for fluid flow simulations, usually within an order of magnitude of the horizontal resolution, which oftentimes results in coarser vertical resolution of the grid. The skilled person will appreciate that the second mesh may not have a coarser vertical resolution grid than the first mesh. A permeability tensor is defined for each cell of the second mesh by identifying all layers of the first mesh present in each cell of the second mesh and applying a homogenisation procedure.

[0077] For fluid flow simulations, it is desirable to maintain the aspect ratio of mesh cells as close as possible to 1 :1 . Therefore, as an example, it is beneficial to avoid the ratio of horizontal resolution: thickness exceeding 5:1. In a non-limiting example, the SFM output (i.e. the first mesh) may have cells of 1 km along X and Y axes but with only approximately 1 cm thickness (Z axis), resulting in horizontal resolution: thickness ratio of 100000:1. After resampling and homogenizing the SFM output (i.e. first mesh) to derive anisotropic permeability in a second mesh, the second mesh may keep the same resolution along X and Y axes (i.e. 1 km in this example), but cell thicknesses along the Z axis in the second mesh may increase to 200m which may cover multiple cells of the first mesh and result in horizontal resolution: thickness ratio of 5:1. Thus, in this example, the new thickness of cells along the Z axis (i.e. 200m) is coarser than the original thickness along the Z axis (i.e. 1 cm).

[0078] Homogenizing and resampling are common terms used in statistical analysis. In the examples of the present techniques, homogenizing means bringing and averaging all data from at least one cell of the first mesh into at least one cell of the second mesh. Resampling refers to interpolation from at least one node or cell of the first mesh to at least one modified node or cell of the second mesh. The step S102 of the method of the present techniques: resampling and homogenizing the first mesh to derive anisotropic permeability in a second mesh also comprises converting the SFM output (i.e. the first mesh) into a permeability distribution to estimate the longitudinal permeability (parallel to the layers) and transverse permeability (perpendicular to the layers) of each sediment class in each cell of the first mesh. Longitudinal permeability is assumed to follow a porosity-permeability relationship, e.g. Kozeny-Carman equation, with transverse permeability being calculated from longitudinal permeability by using a user-provided longitudinal/transverse permeability ratio. In 3D, the longitudinal permeability is assumed to be the same in all directions within the plane aligned with the corresponding layers. Therefore, resampling and homogenizing the first mesh to derive anisotropic permeability in a second mesh S102 also comprises computing values of longitudinal and transverse permeability for each sediment class present in each cell of the first mesh, from user-provided porosity-permeability relationships (e.g. Kozeny-Carman equation) and ratios of longitudinal to transverse permeability specific to each sediment class. In one example, parameter values derived from published porosity-permeability data for the rock types corresponding to the sediment classes in the SFM (e.g. Ehrenberg and Nadeau, 2005; Yang and Aplin, 2010) may be used to calculate permeability.

[0079] Porosity and sediment class fractions are homogenized and resampled onto the new mesh (i.e. the second mesh) by calculating the weighted arithmetic mean of their values in each of the first mesh cells that intersect each cell of the second mesh, where the weights are proportional to the height (thickness) of each first mesh cell.

[0080] Homogenizing and resampling permeability is complex because it is anisotropic. Anisotropy within each cell in the second mesh arises from: The anisotropy of individual sediment classes; and/or the layering of sediment classes within each new cell of the second mesh.

[0081 ] Layering can result in very strong anisotropy, with the permeability potentially being orders of magnitude higher parallel to layering than perpendicular to it. For infinite parallel layers of uniform thickness and permeability, the overall layer-parallel (longitudinal) permeability equals the weighted arithmetic mean of the layer permeabilities, while the overall transverse permeability is the weighted harmonic mean, where the weights are proportional to the thickness of each layer (e.g. Renard and de Marsily, 1997; Matheron 1967). In reality, the longitudinal and transverse permeabilities lie between these bounds due to the heterogeneous and discontinuous nature of the layers. These principles are applied to derive homogenised longitudinal and transverse permeabilities on the second mesh by calculating appropriate weighted means of the longitudinal and transverse permeabilities of the first mesh cells that lie within each cell of the second mesh, where the weights are proportional to the average thickness of each cell of the first mesh. In other words, the cells of the first mesh that lie within each cell of the second mesh are approximated as a series of parallel layers, such that the method described in Matheron, G. (1967). Elements pour une theorie des milieux poreux. France: Eds. Masson et Cie, which is known in the art, can be applied to calculate homogenised longitudinal and transverse permeability. Accordingly, Matheron (1967) is hereby incorporated by reference.

[0082] Accordingly, the equations for calculating homogenised longitudinal and transverse permeabilities of a cell in the second mesh that spans n cells of the first mesh are: where k _long and k _trans are the longitudinal and transverse permeabilities, respectively, of the cell in the second mesh, k _tong is the weighted arithmetic mean of the longitudinal permeabilities of the cells of the first mesh within the cell of the second mesh, k _trans is the weighted harmonic mean of the transverse permeabilities of the cells of the first mesh within the cell of the second mesh, hi is the average thickness of the ith cell of the first mesh in the cell of the second mesh, iong, i and k _{trans i} are the longitudinal and transverse permeabilities, respectively, of the ith cell of the first mesh in the cell of the second mesh, a < 1 is a weighting factor for the homogenised longitudinal permeability and p > 1 is a weighting factor for the homogenised transverse permeability. The values of the weighting factors a and p are chosen by the user to reflect the heterogeneous and discontinuous nature of the layers, such that k _iong is less than or equal to kiong ’ and k _trans is greater than or equal to ktrans- The values of a and p should be chosen such that k _iong > k _rans .

[0083] Subsequently, resampling and homogenizing the first mesh to derive anisotropic permeability in a second mesh S102 further comprises calculating the average orientation of the cells of the first mesh that fall within each cell of the second mesh. In a 3D example the average orientation of the cells of the first mesh that fall within a cell of the second mesh is calculated by using three nodes from the top surface of each cell of the first mesh to define a plane, then using this plane to define a unit normal vector for each cell of the first mesh and calculating the average of these unit normal vectors for all cells of the first mesh in a cell of the second mesh. The skilled person will appreciate that the average orientation calculated in this way will differ slightly depending on which three nodes are chosen to define the plane. Moreover, a slightly different average orientation would be obtained if the three nodes were chosen from the bottom surface of each cell rather than the top surface. However, the skilled person will also appreciate that these differences would have negligible impact on the average orientation given the typically flat shape of the cells of the first mesh. The skilled person will also appreciate that alternative methods could be used to derive the average orientation.

[0084] Resampling and homogenizing the first mesh to derive anisotropic permeability in a second mesh S102 further comprises obtaining a homogenised permeability tensor of each cell in the second mesh by rotation of the diagonal tensor (formed by the homogenised longitudinal and transverse permeabilities in their corresponding average directions) into the global coordinate system. In other words, the average orientation of the cells in the first mesh that fall within each cell of the second mesh is used to rotate the diagonal tensor formed by the longitudinal and transverse permeabilities into the global coordinate system of the second mesh. This process is repeated for each cell in the second mesh. This permeability tensor reflects the local rock composition (i.e. the proportions of each sediment class and their corresponding properties) and orientation from the SFM output (i.e. from the first mesh) which is beneficial for fluid flow simulations. An example of deriving anisotropic permeability tensors in all elements of the second mesh S102 is shown in Figure 3.

[0085] Figure 3 shows a 2D example of permeability homogenisation and its effect on fluid flow. Figure 3 illustrates that fluid flow is guided by the homogenised permeability anisotropy even on a very coarse mesh (the second mesh 210), capturing the anisotropy and geometry implied by the initial/first mesh 200. More specifically, the geometry of the initial mesh 200 is captured by the permeability anisotropy of the second mesh 210 where the permeability anisotropy allows the resulting computed flow 230, which is shown with arrows in Figure 3, to follow the stratigraphy although the second mesh 210 itself does not follow the stratigraphy. This directional fluid flow field represents the solution of the partial differential equations that describe fluid flow in a porous medium with anisotropic permeability, where the permeability tensor is defined in the global (model) coordinate system and is independent of the simulation mesh. [0086] Returning to Figure 1 A, in S104, the method further comprises building a third mesh. The second mesh does not include faults. Moreover, horizontal resolution of the second mesh is insufficient to capture the detail of fluid movement between the faults and sediments. Faults may be used to constrain subsidence at the base of the SFM model and thus the first mesh. However the faults are not represented explicitly within the first mesh, therefore the location of the faults and displacement across faults are not well-resolved in either one of the first mesh or the second mesh. Therefore, a third mesh which represents the faults and other geological discontinuities appropriately is required for the fluid flow simulations.

[0087] In the method of the present techniques, fluid flow simulations are carried out using a fluid flow simulator in which the permeability tensor is defined in the global (model) coordinate system, such that the calculated fluid flow directions are independent of the shape and orientation of the mesh cells. Therefore, it is not necessary for the third mesh to follow the stratigraphic horizons of the SFM output or the layers of the second mesh. Consequently, an unstructured tetrahedral mesh may be used to represent the geological zone of interest. This greatly simplifies the process of generating the third mesh, especially in areas where the intersections of faults and stratigraphic horizons create irregular 3D shapes, and where the stratigraphic layers thin on the margins of the geological zone. For example, MOOSE is an open-source finite element simulator which can be used to import the permeability distribution from the second mesh onto the unstructured third mesh and subsequently used to run fluid flow simulations. The ability to use a mesh that is independent of stratigraphic layering represents a significant advantage over other finite difference codes in which permeability is commonly defined in the local coordinate system of the mesh cells, meaning that the flow direction is meshdependent and the mesh must therefore follow the stratigraphic horizons. The method of the present techniques avoids this limitation by using a fluid flow simulator in which the permeability tensor is defined in the global (model) coordinate system, such that the calculated fluid flow directions are independent of the shape and orientation of the mesh cells.

[0088] In order to provide more context, depending on host rock lithology, the character of a single fault (i.e. a geological discontinuity) can vary with depth and along strike, ranging from conduit to barrier and any combination in between. Representing open flow conduits as lower dimensional elements (see Fig. 4B) can be an efficient method for the discretisation of large numbers of intersecting structural features. Representing low permeability barriers requires a slightly different approach (i.e discontinuous Galerkin approach) as described by Poulet et al. (2021 ). In this discontinuous Galerkin approach, the conforming mesh is split along the fault plane where nodes are duplicated and superposed, with two nodes having different fluid pressure values on either side of the fault (see Fig. 4C). The thickness of the fault is not meshed explicitly but is accounted for when solving for fluid flow and transport. The strategy is modular: all parts can be used separately or simultaneously. For example, a mixed conduit-barrier system can be represented with two lower dimensional surfaces representing conduits (damage zones) surrounding a central interface (fault core) that acts as a barrier to flow (see Fig. 4D). Fig. 5A and Fig. 5B illustrate the equivalence of the volumetric and lower dimensional meshing approaches with a 2D example representing a body of rock containing two intersecting faults 501 , 502, with the more permeable fault 501 overprinting the less permeable fault 502. Figure 5B shows that the resulting fluid pressure distribution is essentially identical, confirming the validity of the lower-dimensional meshing approach.

[0089] Examples of the present invention aim to address at least the following three issues with respect to geological discontinuities modelling: (i) realistic discretisation of complex fault geometries; (ii) easily varying fault behaviour within faults; and (iii) removing the need to use very small elements to discretise fault zones, thus reducing the computational cost of the simulations.

[0090] Returning back to Figure 1 A, S104, building the third mesh further comprises using traces of geological discontinuities and associated meta data (such as dip, dip direction and thickness) to generate corresponding geological discontinuities in the mesh. In certain examples, the third mesh represents at least one geological discontinuity. However, the skilled person will appreciate that the third mesh may still be built without any geological discontinuities if there are no geological discontinuities present. The third mesh may be a 2D or a 3D mesh. The geological discontinuities may be, for example, faults and/or fractures. The associated meta data may be, for example, thickness and/or dip and/or dip direction for each fault or fracture. The geological discontinuities may be meshed as volumetric features of the same dimension as the second mesh. Alternatively, the geological discontinuities may be meshed as lowerdimensional features, such as 2D surfaces within a 3D mesh. The treatment of faults as lowerdimensional features (e.g. 2D surfaces within a 3D model) is advantageous because it reduces significantly the number of elements that are required to represent the faults, thus reducing simulation time and memory requirements. This improves the functioning of the processor of a computer which implements the method of the present techniques. The skilled person will appreciate that the geological discontinuities and their associated meta data may be meshed as lower-dimensional features or as volumetric features of the same dimension as the mesh, based on the type of simulations that will be run on the model created and/or required application. The skilled person will appreciate that the choice of mesh dimension for the geological discontinuities will depend on the type of simulations that will be run on the model created and/or the required application.

[0091 ] In one example, geological discontinuity surfaces are generated from the surface traces of faults within the model domain, with dips and dip directions informed by the user, e.g. interpreted by a skilled person from geophysical data (e.g. gravity, magnetics) and geological data (e.g. surface geology maps). An example may be found following the work of Blaikie and Kunzmann (2020). Therefore, in one example, the third mesh is generated in which the geological discontinuities are represented as 2D surfaces meshed with triangular elements, and the intervening volumes are meshed with tetrahedral elements that are refined in the areas of interest. In one example, a computer script includes parameters to refine the mesh as a function of the distance to faults. Mesh elements closer to the fault will be smaller and gradually enlarge away from the fault, up to a set threshold distance. The threshold distance may be set by a user. All parameters involved can be varied and set by a user.

[0092] The mesh resolution of the third mesh is sufficient to capture spatial variations in permeability, porosity and sediment class fractions within a geological zone even though the third mesh does not follow the layering of the second mesh.

[0093] For instance, the third mesh may be built using a meshing package such as Gmsh techniques described in Geuzaine, C. & Remade, J. -F (2009). Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities. International Journal for Numerical Methods in Engineering, Wiley, 79, 1309-1331 to programmatically generate the third mesh from the user input. In an example, the third mesh is generated by a Python script that uses the programmatic interface to Gmsh. The script performs the following steps: (1 ) use fault traces from a structural analysis to define splines on the top surface of the modelled domain; (2) extrude the splines in the appropriate direction (given by dip and dip direction provided by the user or in the metadata for the faults) to create 2D surfaces extending from the top to the bottom of the modelled domain; (3) optionally, if the user wishes to mesh the faults volumetrically, the 2D surface defining each fault is translated by half the fault thickness (which is specified by the user or in the metadata for the fault) in an average direction normal to the fault, then the translated surface is extruded in the opposite direction over the full thickness of the fault to create a corresponding volume centred on the original 2D surface; (4) perform topological operations to compute all intersections between all surfaces and volumes; (5) tag all surfaces and volumes such that each region of the mesh can be identified for the purpose of assigning properties (e.g. permeability); (6) generate an unstructured mesh comprising triangular cells on surfaces and tetrahedral cells in volumes; (7) save the mesh as a file suitable for reading into a fluid flow simulation code.

[0094] Optionally, the method of the present techniques may further comprise testing the third mesh to ensure that errors due to mesh discretization did not affect the results. The testing may be done by generating a series of meshes of varying resolutions and comparing the results of fluid flow simulations run on each of these meshes.

[0095] In one example, in fluid flow simulations, given the uncertainty regarding fault permeability in unconsolidated sediments, the method of the present techniques may explore the possibility of the faults being either more or less permeable than the surrounding sediments.

[0096] Figure 6 shows a 3D example of building the third mesh 406 from geological discontinuities 400 (here faults) defined by their traces at the surface and their meta-data (thickness, dip, dip direction). In this example, traces of geological discontinuities and associated meta data are provided in 1 D in block 400. In this example, the geological discontinuities are meshed as volumetric features of the same dimension as the third mesh in block 402. Alternatively or additionally, the geological discontinuities may be meshed as lowerdimensional features, such as 2D surfaces within a 3D mesh as shown in a block 404. Subsequently, the third mesh is obtained as shown in blocks 404 and 406. Block 406 illustrates the third mesh for volumetric faults and block 404 illustrates the third mesh for the case where faults are modelled as lower dimensional features.

[0097] Returning back to Figure 1 , in S106, the method further comprises setting anisotropic permeability in the third mesh. Therefore, the method comprises setting anisotropic permeability in the third mesh from the derived anisotropic permeability in the second mesh by interpolating from the second mesh to the third mesh and thus creating the third mesh with set anisotropic permeability. In an example, the method comprises setting anisotropic permeability in the third mesh from the permeability field computed in the second mesh. In this example, this may be achieved by interpolation from the second mesh to the third mesh.

[0098] In the geological discontinuities in the third mesh, the permeability is a user-specified function of the properties of the adjacent geological zone (e.g. host rock or host sediment) which are obtained from the second mesh. In an example, a multiple of the longitudinal permeability of the host rock (or host sediment) is used to provide the permeability of the geological discontinuity, and the anisotropy may be either retained from the host rock (or host sediment) or rotated to match the orientation of the geological discontinuity. The permeability of the geological discontinuity may also depend on the corresponding sediment class fractions at that location. Geological discontinuity is at least one geological fault and/or at least one geological fracture.

[0099] In a case where the anisotropic permeability tensor is rotated to match the orientation of the geological discontinuity, the longitudinal and transverse directions used to rotate the permeability tensor are obtained from the gradient of a field that can be obtained (a) either by solving numerically a diffusion problem for a fake property that is set to a constant value on a surface at the centre of the geological discontinuity, or (b) by using third-party implicit modelling software producing such potential field. The gradient of the resulting field is used to determine the transverse direction of the geological discontinuity, with the longitudinal directions set in two perpendicular directions.

[00100] Therefore, the third mesh may be populated with porosity and permeability information from the second mesh, using an interpolation algorithm to map values from the second mesh onto the third mesh. The third mesh with set anisotropic permeability may be a tetrahedral mesh.

[00101 ] Returning back to Figure 1 , the method of the present techniques may further comprise step S108 of determining permeability distributions in a geological region based on the third mesh with set anisotropic permeability. Therefore, the skilled person will appreciate that the third mesh with set anisotropic permeability is suitable for use in fluid flow simulations.

[00102] For instance, fluid flow simulations can be performed using the PorousFlow module of the open-source finite element code MOOSE (Multiphysics Object-Oriented Simulation Environment). The simulated domain may be assumed to be a continuous porous medium that is fully saturated with water, with fluid flow obeying Darcy’s law with anisotropic permeability. The skilled person will appreciate that any other similar package may be used.

[00103] The skilled person will appreciate that each one of the first mesh, second mesh, third mesh and the third mesh with set anisotropic permeability have at least one cell.

[00104] Beneficially, the method of the present techniques allows building a 3D model (i.e. the third mesh with set anisotropic permeability) which captures the permeability anisotropy distribution from a stratigraphic forward model (first mesh) and structural framework, in a third mesh with set anisotropic permeability which is suitable to run fluid flow simulations. One of the advantages of the method of the present techniques is that by using the SFM to constrain permeability of both faults and sediments, the method of the present techniques predicts fluid pathways through a geological zone under various fault permeability scenarios. The subtleties of the fluid pathways predicted by the method of the present techniques cannot be predicted by models that assign uniform properties, including permeability, to faults and their host rocks. Advantageously, the method of the present techniques may be used in the early stages of mineral and hydrocarbon exploration, where borehole data tend to be quite sparse.

[00105] Further exemplary benefits of the method of the present techniques over existing solutions, including over the separate components used here in a combined manner are as follows: Stratigraphic Forward Modelling, a type of process-based modelling used in the method of the present techniques usually only includes major elements of the structural framework with corresponding limitations about the structural framework itself (e.g. vertical faults only, smaller faults ignored) and without the ability to set material properties within the generated geological discontinuities, even when those features are thick enough (e.g., fault zone) to be represented in the corresponding mesh. Moreover, 3D geological models known in the art generally use geostatistical distributions of rock properties. In other words, the 3D models known in the art do not account for the physics of the processes responsible for heterogeneous property distributions. This can result in geologically unrealistic property distributions. Furthermore, permeability information within geological faults known in the art is usually set per fault as a constant, sometimes as a function of depth, but rarely as a function of the material properties of the corresponding host rock and rarely with anisotropy following the orientation of the fault. Moreover, in many geological fluid flow simulators known in the art, the longitudinal and transverse permeability directions are assumed to be parallel and perpendicular to the local coordinate axes of the mesh; that is, anisotropy must follow the layers in the mesh, and the mesh layers must therefore follow stratigraphic layering in order to generate realistic fluid flow fields. Accordingly, it is difficult to construct such meshes for complex 3D architectures, especially those involving folded layers and faults. In some other geological fluid flow simulators known in the art, the longitudinal and transverse permeability directions are assumed to be horizontal and vertical, resulting in fluid flow fields that do not reflect the orientation of stratigraphic layers and their anisotropy. These problems may be resolved by simulators that provide full tensor representation of permeability (as in the method of the present techniques), as the direction of fluid flow is then independent of the mesh, however this functionality is not commonly used because it is difficult to determine the orientation of the permeability tensor throughout the modelled domain. This difficulty is resolved by the method of the present techniques.

[00106] The skilled person will appreciate that the method of the present techniques may be applied by geological modellers trying to understand the impact of the structural framework and physics-based stratigraphic information on the movements of fluids underground, which has applications in a variety of fields, including, the following, non-limiting examples: (a) mineral exploration, where fluid flow simulations on 3D geological models generated by the present techniques could be used to understand mineralisation patterns; (b) hydrocarbon reservoir engineering (conventional or unconventional), where the method of the present techniques could be used to predict more accurately the formation of hydrocarbon reservoirs and/or their behaviour during production; (c) geothermal energy, where the method of the present techniques could be used to simulate underground circulation of fluids advecting heat; (d) safe waste disposal, where the method of the present techniques could be used to predict the movement of potentially contaminated fluids around waste buried in fractured or faulted geological environments; (e) groundwater contamination, where the method of the present techniques could be used to capture more accurately the plume propagation of contaminants underground.

[00107] Therefore, the computer-implemented method of the present techniques comprises features representing an interaction with an external physical reality by building a third mesh with set anisotropic permeability which can be used for simulating fluid movement in a real geological zone. Therefore, the computer-implemented method of the present techniques predicts the physical state of existing real objects (e.g. the state of fluid flow accounting for the geological discontinuities) in a geographical zone and thus makes a technical contribution regardless of what use is made of the results.

[00108] Returning back to Figure 1A, optionally, the method of the present techniques may further comprise step S1 10 of applying the output of the determined fluid flow to a processor configured to perform various functions and/or control hardware. For instance, the output of the present method may be used in a hardware which can determine a drilling target zone. For instance, the output of the present method may be used in a hardware which can determine a drilling target zone for mineral and/or hydrocarbon exploration. In another example, the output of the present method may be used in a hardware which can determine a drilling target zone for operation of a drilling machine for mineral and hydrocarbon exploration. In an example, the drilling machine may be configured to drill a geographical zone according to the output of the method of the present techniques.

[00109] In one example, the method of the present techniques may comprise step S1 10 of applying the output of the determined fluid flow to a processor configured to perform various functions and/or control hardware. For instance, the output of the present method may determine a mineralisation distribution (e.g. iron ore distribution) in a target geological zone. Beneficially, this result may be input into a hardware which can coordinate drilling and mineral exploration based on the result of the method of the present techniques.

[001 10] The skilled person will appreciate that alternatively, the method of the present techniques may comprise a step S1 10 of applying the output of the determined fluid flow to a processor configured to map mineralisation patterns in a geological zone. In this example, machinery may be operated based on the results of the mapped mineralisation patterns.

[001 1 1 ] Alternatively, the method of the present techniques may comprise a step S1 10 of applying the output of the determined fluid flow to a processor configured to predict the formation of hydrocarbon reservoirs and/or their behaviour during production. In this example, a piece of machinery may be operated based on the result of the predicted formation of hydrocarbon reservoirs and/or their behaviour during production.

[001 12] Alternatively, the method of the present techniques may comprise step S1 10 of applying the output of the determined fluid flow to a processor configured to simulate the underground circulation of fluids advecting heat in a geological zone of interest. In this example, a drilling machine may be operated based on the simulated result of the underground circulation of fluids advecting in a geological zone of interest.

[001 13] Alternatively, the method of the present techniques may further comprise step S1 10 of applying the output of the determined fluid flow to a processor configured to predict the movement of potentially contaminated fluids around waste buried in fractured or faulted geological environments. In this example, a piece of machinery may be operated based on the result of the predicted movement of potentially contaminated fluids.

[001 14] Alternatively, the method of the present techniques may further comprise step S1 10 of applying the output of the determined fluid flow to a processor configured to capture the plume propagation of contamination underground. In this example, a machinery may be operated based on the result of the predicted plume propagation of contamination underground.

[001 15] Figure 7 is a schematic illustration of the system which can be used to implement the method of the present techniques. The method of the present techniques as described in Figure 1 A is a computer-implemented method. This method may be implemented on a computing device 500.

[001 16] Figure 7 comprises a computing device 700, which may be any suitable electronic device, e.g. a personal computer or computing device, a laptop, a tablet, a smart phone etc. It will be understood that this is a non-exhaustive and non-limiting list of example devices. Figure 7 shows some of the components of the computing device 700, it will be appreciated that there may additionally be other standard components which are not shown. The apparatus comprises at least one processor 702 coupled to memory 704. The at least one processor 702 may comprise one or more of: a microprocessor, a microcontroller, and an integrated circuit. The memory 704 may comprise volatile memory, such as random access memory (RAM), for use as temporary memory, and/or non-volatile memory such as Flash, read only memory (ROM), or electrically erasable programmable ROM (EEPROM), for storing data, programs, or instructions, for example. The computing device 700 typically comprises at least one user interface 706 for a user to input parameters of the geological zone of interest. The at least one user interface 706 may take any appropriate form, e.g. a keyboard, a mouse, a touchpad or other input device or a computer screen. Similarly, the computing device 700 typically comprises at least one display 708 for providing the results and/or data generated during the method described above. Exemplary case study

[001 17] The method of the present techniques, as described in Figure 1A, will now be illustrated on an exemplary case study which focuses on the Glyde sub-basin of the McArthur Basin in northern Australia.

[001 18] Determining the permeability distribution in the basin fill and faults is particularly challenging where outcrops are lacking and drillhole data are sparse. This approach is illustrated using a case study in the McArthur Basin (northern Australia) where drillholes provide sufficient data to constrain a stratigraphic forward model (SFM) of the Palaeo-proterozoic Barney Creek Formation (BCF), host to several economically important Zn-Pb-Ag deposits. The SFM is used to derive a permeability distribution in the BCF and related faults, which is then used in fluid flow simulations to demonstrate the impact of the sedimentary facies distribution on fluid flow.

[001 19] The Paleo- to early Meso-proterozoic McArthur Basin formed in an intra-cratonic setting on the North Australian Craton. Known sediment-hosted mineralisation in the basin includes the world-class McArthur River Zn-Pb-Ag deposit and the smaller Teena deposit, which occur in the BCF of the McArthur Group. The BCF comprises dolomitic, organic-rich and pyritic shales and siltstones with minor mass-flow breccias and other carbonate rocks. Deposition of the BCF occurred primarily in sub-basins formed by sinistral strike-slip movement on the Emu Fault.

[00120] The timing of mineralisation relative to deposition of the BCF has long been debated, with some authors arguing for a syngenetic, SEDEX-style origin on the seafloor, and others arguing for diagenetic or epigenetic mineralisation within the sediments. These contrasting genetic models have different implications for fault permeability: Syngenetic mineralisation implies that the faults that carried the mineralizing fluid were more permeable than the surrounding rocks/sediments all the way to the seafloor, whereas diagenetic/epigenetic mineralisation implies that fluid deviated from the faults into the adjacent sediments, which in turn implies that the faults changed from fluid conduits to barriers at some depth. A detailed knowledge of the permeability distribution within the faults and sediments is necessary in order to test these conceptual models, or to show where different styles of mineralisation may have occurred. [00121 ] Faults that cut through sedimentary sequences typically display large variations in permeability as they intersect layers with different mechanical properties, grain size and mineralogy. To predict these variations requires knowledge of the distribution of sediment/rock types and properties along the fault. Here, it is demonstrated how this may be achieved in the absence of outcrops or dense drillhole data, using the Glyde sub-basin of the McArthur Basin, northern Australia to illustrate the method of the present techniques.

[00122] The summary information regarding the McArthur Basin can be seen in Figure 8. The Proterozoic McArthur Basin is part of one of the largest zinc-lead provinces in the world, with deposits such as McArthur River and Teena (Fig. 8A) occurring in the Barney Creek Formation (BCF) of the McArthur Group. The BCF is dominated by dolomitic, carbonaceous, and pyritic shales and siltstones, deposited primarily in subbasins along the Emu Fault, overlying dolomitic carbonates and evaporites of the lower McArthur Group. Faults are believed to have played an important role in transporting metal-bearing fluids from deeply buried metal source rocks to shallow depths, where mineralisation occurred either within the sediment pile or on the seafloor. Of particular interest is the path taken by ascending fluids where the faults intersect the BCF. A detailed knowledge of the permeability distribution within the faults and sediments is necessary in order to test hypotheses about the timing of mineralisation relative to sediment deposition.

[00123] There is no known economic mineralisation in the Glyde sub-basin but there are sufficient boreholes with publicly accessible data (see Fig. 8A) to constrain the sedimentary sequence within the BCF, which in turn can be used to estimate fault permeability. Kunzmann et al. (2019, 2022) used drill core descriptions (Fig. 9A) to conduct a facies analysis of the BCF, from which they developed a sequence stratigraphic framework for the depositional system (Fig. 8B). The case study of the present invention utilised their analysis along with the structural framework proposed by Blaikie and Kunzmann (2020) to constrain a Stratigraphic Forward Model of the BCF in the southern McArthur Basin. Stratigraphic forward modelling reconstructs the stratigraphic architecture of basins by simulating erosion, transportation, and deposition of sediments in response to basement subsidence/uplift, sea level fluctuations and variations in sediment supply and thus creates the first mesh.

[00124] Stratigraphic forward modelling simulations were performed using DionisosFlow. The stratigraphic forward modelling includes four sediment classes: sand, silt, clay and carbonate. Regional scale faults (See Fig. 8) were used to constrain subsidence at the base of the model, however note that faults are not represented explicitly within the SFM/first mesh. [00125] The stratigraphic forward modelling produces a 4D distribution of porosity and sediment class fractions over the modelled domain (See Fig. 9C and Fig. 11 ). Using the final state of the SFM we estimated longitudinal permeability (parallel to the depositional layers within the SFM) from porosity using a Kozeny-Carman type relationship parameterised for each sediment class. Transverse permeability (perpendicular to the layers) was assumed to be a fraction of longitudinal permeability, with the multiplier varying between sediment classes. The resulting porosity-permeability distribution was then resampled onto a mesh (i.e. the third mesh with set anisotropic permeability) suitable for fluid flow simulation in the finite element code MOOSE, which included faults represented as 2D surfaces. The model domain extends down to a horizontal surface at 2 km depth (below the deepest part of the BCF), with the area below the BCF being assigned relatively low permeability representing consolidated rocks of the lower McArthur Group (see Fig. 14).

[00126] In the area below the BCF fault permeability was assumed to be higher than that of the host rock, consistent with the behaviour of faults in consolidated rocks. The faults were assigned a thickness of 100 m.

[00127] In petroleum geology, much attention is focused on predicting the sealing capacity of faults in reservoir rocks. Conventional methods for estimating sealing capacity, such as the Shale Gouge Ratio (SGR), focus on the effect of clay being entrained (smeared) along the fault, which reduces cross-fault permeability. However, such methods do not account for entrainment/injection of sand into the fault plane, which may enhance fault-parallel permeability in faults that intersect alternating layers of unconsolidated sandy and clay-rich sediments. Such faults act as combined conduit-barrier systems due to the combined effects of sand entrainment/injection and clay smearing.

[00128] Given the uncertainty regarding fault permeability in unconsolidated sediments, we explore the possibility of the faults being either more or less permeable than the surrounding sediments where the faults intersect the BCF in our fluid flow simulations.

[00129] To illustrate the present invention, fluid flow through the domain was simulated by imposing a constant fluid flux of 10' ⁸ m ³ /m ² /s at the base of the faults (representing some larger- scale process driving fluid flow through the system) and a fixed hydrostatic fluid pressure on the top boundary, then solving for the steady-state fluid pressure distribution. The corresponding fluid flux was calculated using Darcy’s law. Workflow

Sedimentary and stratigraphic analysis

[00130] The case study to illustrate the present invention is based on the facies analyses of the middle McArthur Group presented in Kunzmann et al. (2019). Based on drill core descriptions as shown in Figure 9A, the inventors identified four main facies associations in the interval from the Tooganinie Formation to the Hot Spring Member of the Lynott Formation (Fig. 8B). The depositional environments range from continental - supratidal to deep subtidal - slope, including intertidal and subtidal facies. Studying facies variations together with their geometric relationships enabled Kunzmann et al. (2022) to propose a detailed sequence stratigraphic framework for this interval, in which five 3 ^rd order depositional sequences were identified (Fig. 8B). In the Glyde sub-basin, the BCF is dominated by sediments deposited in offshore transition to offshore environments, with a small thickness of shallow water sediments (the W-Fold Shale) at the base. Shale content varies within the BCF; one shale-rich interval occurring at the maximum flooding surface of the B1 sequence (Fig. 8B) has been interpreted to play a major role in the formation of Zn-Pb-Ag deposits due to its role as either a barrier to upward fluid flow in a diagenetic/epigenetic mineralisation scenario, or a chemical trap on the seafloor in a syngenetic mineralisation scenario (Kunzmann et al., 2019).

Stratigraphic forward modelling

[00131 ] The facies analysis described above was used along with the structural framework of the southern McArthur Basin proposed by Blaikie and Kunzmann (2020) to constrain a SFM corresponding to depositional sequences B1 and B2 (i.e. the HYC Pyritic Shale and Undifferentiated unit of the BCF;). Stratigraphic forward modelling aims at reconstructing the stratigraphic architecture of basins by simulating the erosion, transportation, and deposition of sediments. The present example uses DionisosFlow, a stratigraphic forward modelling tool developed by IFPen, to simulate deposition of sediments in the southern McArthur Basin. Amongst other processes, DionisosFlow simulates the diffusion of multiple sediment classes with different transport properties. In DionisosFlow, sedimentary material displacement is calculated using a water- and gravity-driven diffusion equation:

Q _s . = -K _s . * h + K _w .Q™S ⁿ (1 ) [00132] where Q _Si is the flux of sediment class i; K _s ., K _w .are the hill-slope creeping transport and water-driven diffusion coefficients for sediment class i (km ² /a); h is elevation (m); Q _w is water flux; S is the slope; and m and n are empirical constants. Equation 1 accounts for linear slope- driven diffusion and non-linear slope- and water-driven diffusion. To accurately represent the stratigraphic evolution of basins, DionisosFlow accounts for the main geological controls on accommodation creation and evolution, namely sea level variations, changes in subsidence/uplift rates, variations in sediment source locations and supply, carbonate precipitation, and compaction. DionisosFlow’s output for this study is a 4D distribution of porosity and sediment class fractions, which provides the necessary detail to estimate the permeability distribution in the modelled area. As such, stratigraphic forward modelling is a useful method for estimating a geologically realistic permeability distribution in the absence of dense drillhole data, as is often the case in mineral exploration.

SFM parameterization and calibration.

[00133] For the present case study, the inventors built a 100 by 200 km SFM (red box in Fig. 8A) with a horizontal cell size of 1 km. The SFM area is large enough to minimize boundary effects in the areas of interest (e.g. the Glyde sub-basin) and spans from proximal to distal parts of the sedimentary system. The SFM has four sediment classes: three detrital sediments (sand, silt, shale), and one carbonate. In this model, carbonate production peaks in environments with 70 m water depth and is limited to areas with limited detrital sediment delivery, to represent the chemically precipitated nature of these carbonates. Detrital sediment supply enters the model by the western and southern edges, where previous works identified the continent to be located. High frequency eustatic variations are introduced by a sinusoidal curve with an amplitude of 10 m and a period of 2 Ma (million years).

[00134] Simulation of the basin evolution is carried out from 1639 Ma to 1624 Ma, which are interpreted to respectively represent the base of the BCF pyritic shale member and top of the undifferentiated BCF (Kunzmann et al., 2022). Sediment transport is parameterized to represent displacement by water, gravity, and waves. The main structural features, as shown in Figure 8A, are used to create subsidence maps, which in turn control the rate of subsidence in the subbasins. Note that faults are not explicitly represented within the SFM, and the relatively coarse horizontal grid size (1 km) of the SFM means that displacement across faults is not well- resolved. [00135] Four intervals, whose thicknesses and final bathymetries are derived from drill core data (Fig. 8B) in the modelled area (including the Glyde sub-basin), are used to calibrate the model. Calibration is carried out iteratively by adjusting the sediment supply, carbonate production rate and subsidence/uplift rates during the four intervals, aiming for an overall misfit less than 5% between measured and modelled thicknesses and bathymetries. It should be noted that, while the SFM extent covers the locations of known mineralisation at McArthur River and Teena, the lack of publicly available drill hole data in those areas means that the model is not well-calibrated there. Hence focus is made on the Glyde sub-basin for the present study. The temporal evolution of the stratigraphic model in this area is illustrated in Fig. 9C, which shows the distribution of the shale sediment class at four output times.

Stratigraphic forward modelling output

[00136] The stratigraphic forward modelling output is in the form of a structured mesh (i.e. the first mesh), comprising stacked layers of 1 x 1 km cells with varying thickness as shown in Fig. 10A. Each layer represents the net result of sediment deposition, erosion and compaction over a given time interval. The thickness of the layers varies over the model domain, and can be zero if no sediment was deposited or all sediment was eroded in a given cell. Each cell contains the porosity and the proportion of each sediment class as shown in Figure 1 1 . There are high porosities (> 30%) throughout the sediment pile, reflecting the unconsolidated nature of the sediments, and an area of high shale content on the east side of the sub-basin; this shale plays a key role in controlling fluid flow through the sediments.

Permeability estimation

[00137] The first step in converting the SFM output (i.e. the first mesh) into a permeability distribution is to estimate the longitudinal permeability (parallel to the layers) and transverse permeability (perpendicular to the layers) of each sediment class in each cell of the first mesh. Longitudinal permeability was assumed to follow a Kozeny-Carman relationship:

[00138] where is the longitudinal permeability of sediment class j in cell i; k _o is the permeability corresponding to porosity < > _Oj for sediment class 7; and 3 < n, < 7 is an empirical constant. Transverse permeability was assumed to be related to longitudinal permeability: [00139] where k^ ^ans is the transverse permeability of sediment class j in cell i and cr _; > 1 is the ratio of longitudinal to transverse permeability for sediment class j (i.e. the anisotropy). Table 1 lists the parameter values that were used to calculate permeability. These parameter values yield permeabilities consistent with published datasets for the relevant sediment types (e.g. Ehrenberg and Nadeau, 2005; Yang and Aplin, 2010).

Table 1. Sediment class parameters used for permeability estimation

Resampling

[00140] It is not possible to import information from the SFM output directly into a fluid flow simulation due to its fine vertical resolution and the existence of zero-thickness cells. This problem is addressed by resampling the SFM output onto a new structured mesh (i.e. the second mesh) with the same horizontal resolution but coarser vertical resolution as shown in Figure 1 1 B. The new mesh (i.e. the second mesh) retains some key horizons from the original SFM output mesh and is extended down to a horizontal surface at 2 km depth which corresponds to the base of the fluid flow simulations. The area below the SFM is assigned a relatively low porosity and permeability, representing the more consolidated sediments underlying the BCF.

[00141 ] Porosity and sediment class fractions are homogenized and resampled onto the new mesh (i.e. the second mesh) by calculating the weighted arithmetic mean of their values in each old mesh cell (i.e. the first mesh) that intersects each new mesh cell (i.e. the second mesh), where the weights are proportional to the height (thickness) of each old cell (i.e. the first mesh). The original and resampled porosity distributions in part of the mesh are illustrated in Fig. 12.

[00142] Resampling permeability is more complex because it is anisotropic. Anisotropy within each cell in the new mesh arises from:

1 . The anisotropy of individual sediment classes;

2. The layering of sediment classes within each new cell. [00143] Layering can result in very strong anisotropy, with the permeability potentially being orders of magnitude higher parallel to layering than perpendicular to it. For infinite layers of uniform thickness and permeability the overall layer-parallel (longitudinal) permeability equals the weighted arithmetic mean of the layer permeabilities, while the overall transverse permeability is the weighted harmonic mean. In reality the longitudinal and transverse permeabilities lie between these bounds due to the heterogeneous and discontinuous nature of the layers. These principles are applied to derive longitudinal and transverse permeabilities on the new mesh (i.e. the second mesh) by calculating appropriate means of the longitudinal and transverse permeabilities of the old grid cells (i.e. cells of the first mesh) that lie within each new grid cell (i.e. cells of the second mesh). In addition the average direction of the cells of the first mesh that lie within each cell of the second mesh is used to convert the average longitudinal and transverse permeabilities to components of the permeability tensor in the global coordinate system (required for the fluid flow simulations). The resampled permeability distribution is shown in Fig. 13.

Fluid flow simulations

[00144] Fluid flow simulations are performed using the PorousFlow module of the open- source finite element code MOOSE (Multiphysics Object-Oriented Simulation Environment). The simulated domain is assumed to be a continuous porous medium that is fully saturated with water, with fluid flow obeying Darcy’s law with anisotropic permeability.

[00145] For the case of fully-saturated flow of a single-phase fluid, MOOSE solves the continuity equation:

[00146] where < > is porosity, p _f the fluid density (kg m ^-3 ), and the Darcy fluid flux (a volumetric flux; m s ^-1 ):

[00147] where k is the permeability tensor (m ² ), p the fluid viscosity (Pa s), P the fluid pressure (Pa) and g (m s ^-2 ) the acceleration due to gravity. [00148] To track movement of fluid through the model we introduce an inert tracer in the pore fluid with mass fraction X. The tracer is transported by advection with the moving pore fluid, represented by the following mass balance equation:

[00149] Further information on the governing equations, validation tests and numerical solution techniques can be found at

Simulation mesh

[00150] The resampled SFM mesh (i.e. the second mesh) does not include the faults, and its horizontal resolution is insufficient to capture the detail of fluid movement between the faults and sediments. An alternative mesh is therefore required for the fluid flow simulations. It is not necessary for the simulation mesh to follow the stratigraphic horizons of the SFM output because the permeability tensor in MOOSE is defined in the global (model) coordinate system, such that the calculated fluid flow directions are independent of the shape and orientation of the mesh cells. Consequently, an unstructured tetrahedral mesh is used to represent the basin fill. This greatly simplifies the process of mesh generation, especially in areas where the intersections of faults and stratigraphic horizons create irregular 3D shapes, and where the stratigraphic layers thin on the margins of the sub-basin. The ability to use a mesh that is independent of stratigraphic layering represents a significant advantage over finite difference codes in which permeability is commonly defined in the local coordinate system of the mesh cells, meaning that the flow direction is mesh-dependent and the mesh must therefore follow the stratigraphic horizons.

[00151 ] Fault surfaces are generated from the surface traces of faults within the model domain (Fig. 8A), with dips and dip directions informed by the work of Blaikie and Kunzmann (2020). A mesh (i.e. the third mesh) is generated in which the faults are represented as 2D surfaces meshed with triangular elements, and the intervening volumes are meshed with tetrahedral elements that are refined in the areas of interest, i.e. close to the faults and within the BCF (Fig. 14).

[00152] The resolution of the third mesh is sufficient to capture spatial variations in permeability, porosity and sediment class fractions within the BCF even though the mesh does not follow the layering of the resampled SFM. The treatment of faults as lower-dimensional features (i.e. 2D surfaces within a 3D model) in the third mesh is advantageous because it reduces significantly the number of elements that are required to represent the faults, thus reducing simulation time and memory requirements. The mesh is generated using Gmsh, and a series of meshes of varying resolution were tested to ensure that errors due to mesh discretization did not affect the results. The final mesh comprises 38.7 million elements and 6.3 million nodes.

Porosity and permeability

[00153] The tetrahedral elements of the third mesh are populated with porosity and permeability from the resampled SFM (i.e. the second mesh), using an interpolation algorithm within MOOSE to map values from the resampled SFM mesh onto the tetrahedral mesh. Fault permeability within the BCF is assumed to be related to the longitudinal permeability and shale content ( _sha(e ) of the adjacent sediments, such that k _fault = ak _long , where a is a multiplying factor. 4 permeability scenarios were investigated:

1 . Fault permeability > sediment permeability: a = 10

2. Fault permeability < sediment permeability: a = 0.01

3. Fault permeability decreasing with increasing shale content (X _shaie y. toga = 2 — 4 • min(X _sfta(e , 0.6)/0.6.

This relationship results in a = 100 when X _sfiaie = 0, and a = 0.01 when X _sfiaie > 0.6.

4. Fault permeability as in scenario 2, but with principal permeability directions of the host rocks/sediments assumed to be horizontal and vertical:

[00154] In all scenarios we assume a = 100 for fault permeability for the parts of the faults below the BCF, consistent with the relatively high permeability of faults that cut through consolidated sedimentary rocks. The resulting fault permeability distributions are illustrated in Fig. 15 (left column; grey-red colour scale in faults). The skilled person will appreciate that the fault permeability relationships used in the present example are just some of many possible relationships that could be implemented.

[00155] For the present investigation fault permeability is treated as isotropic; while this is an over-simplification, it is unlikely to be significant in the simulated scenarios where the dominant flow pathway is upwards along the faults and out into the host sediments. Initial and boundary conditions

[00156] Fluid pressure is initialized to hydrostatic and fixed at 10 ⁵ Pa on the top boundary. A fluid flux of 10 ^-8 m ³ m ^-2 s ^-1 is applied at the base of the faults, representing fluid being driven upwards through the faults by larger-scale processes, such as thermal convection. The bottom and sides of the model are treated as impermeable boundaries, except where the fluid flux is applied at the base of the faults.

[00157] Equation (4) is solved for steady-state fluid pressure subject to these boundary conditions, and the resulting Darcy fluid flux is calculated using Equation (5). Equation (6) is then used to simulate tracer transport in the steady-state fluid flow field for a duration of 2 x 10 ¹¹ seconds, with a boundary condition of X = 0.1 at the base of the faults and initial concentration X = 0 throughout the model. The simulated duration of tracer transport was chosen to be sufficient to demonstrate movement of fluid through the model; it does not have any geological significance.

Results

[00158] The first two rows in Fig. 15 illustrate the effect of fault permeability being either higher (Scenario 1 ) or lower (Scenario 2) than that of the host sediments. When fault permeability is high, fluid flows uninterrupted through the faults to the top of the model and has only minimal interaction with the sediments of the BCF. The tracer remains largely confined to the faults in this scenario. Conversely, when fault permeability is lower than that of the BCF sediments, fluid deviates from the faults once it reaches the BCF. This is reflected in the tracer distribution, which extends into the sediments following the orientation of the stratigraphic layers. This happens on both sides of the sub-basin. In Scenario 3 (third row in Fig. 15) the fault permeability is a function of shale content, being lowest where shale content is high on the East side of the sub-basin. Consequently, in this scenario fluid flowing up the West fault continues largely uninterrupted to the top of the model, whereas fluid flowing up the East fault deviates into the sediments.

[00159] Scenario 4 (fourth row in Fig. 15) illustrates the importance of considering the direction of permeability anisotropy in the sediments as well as its magnitude. Fault permeability is the same as in Scenario 2, but the principal permeability directions in the BCF are assumed to be horizontal and vertical rather than following the stratigraphic layering. The effect of this simplified treatment of anisotropy is most apparent in the tracer distribution, which shows that the tracer has moved horizontally away from the faults rather than following the dip of the layers as in Scenarios 2 and 3.

[00160] The subtleties of the fluid pathways predicted by the present method cannot be predicted by models that assign uniform properties to faults and their host rocks. The predicted fluid pathways in Scenario 3 of Figure 15 are consistent with the fluid flow regime associated with diagenetic mineralisation elsewhere in the basin, where black shales associated with a maximum flooding surface near the base of the BCF are believed to have trapped ascending fluid and forced it to travel through the sediments, depositing metals there. Conversely, scenarios with high fault permeability (i.e., Scenario 1 of Figure 15 and the West fault in Scenario 3 of Figure 15) predict fluid pathways that are more consistent with syngenetic mineralisation on the seafloor. Exploring the sensitivity of the model to varying fault permeability scenarios enables testing of hypotheses concerning mineralisation or other fluid-related phenomena in the basin.

[00161 ] The importance of having permeability anisotropy following stratigraphic layering (illustrated by comparison of scenarios 2 of Figure 15 and 4 of Figure 15) was noted previously in literature by Yager et al. (2009), who showed that groundwater ages and flow pathways in a folded aquifer system differed significantly between models that used different representations of permeability anisotropy. The ability to predict both the magnitude and orientation of the permeability tensor is therefore an important feature of the SFM.

[00162] Stratigraphic forward modelling is particularly useful in scenarios where there are sufficient boreholes and/or outcrop data to constrain a SFM, but not sufficient to construct a detailed facies model simply by interpolating between data points. As such, it is anticipated that opportunities to apply this technique will arise in the early stages of mineral and hydrocarbon exploration, where borehole data tend to be quite sparse. Stratigraphic forward modelling also provides an alternative to various statistical approaches that are used to provide permeability distributions within aquifers in groundwater models. Other areas where the present approach could be applied include modelling of contaminant dispersal and CO2 sequestration; a realistic permeability distribution is important in both cases.

[00163] The fluid flow simulations presented in this application were conducted using the finite element code MOOSE. However the workflow is not confined to a specific code; it could be implemented using other simulation codes that permit unstructured grids, such as SUTRA. Unstructured grids with tensorial representation of permeability enable modelling of arbitrary 3D geometries with permeability anisotropy independent of the mesh, which is critical for representing the small-scale permeability variations within sediments and the permeability contrasts between faults and their host rocks/sediments.

[00164] Transitioning to codes that use unstructured grids may be particularly beneficial in the field of groundwater modelling, in which the representation of faults has historically been hindered by the use of finite difference codes such as MODFLOW that make it difficult to represent steeply dipping structures intersecting stratigraphic layers. The present approach opens the door to a more realistic representation of permeability in both faults and host rocks in groundwater flow models.

[00165] The permeability of faults and their host rocks/sediments is a critical input for models of fluid flow in sedimentary basins. This study demonstrates the use of stratigraphic forward modelling to estimate permeability of the basin fill, including both its magnitude and the direction of anisotropy. This, in turn, is used to provide the permeability distribution in fluid flow simulations, including both the sediments and associated faults. The simulations reveal subtle aspects of the fluid flow pathways through a basin that could not be predicted by models using uniform permeabilities. The proposed workflow has a range of potential applications, including mineral/hydrocarbon exploration, groundwater studies, contaminant dispersal modelling and CO2 sequestration.

[00166] Therefore, the skilled person will appreciate that embodiments of the present invention can be applied by geological modellers trying to understand the impact of the structural framework and physics-based stratigraphic information on the movements of fluids underground. This has applications in a variety of fields, including: Mineral exploration, where fluid flow simulations on 3D geological models generated by this technology could be used to understand mineralisation patterns; Hydrocarbon reservoir engineering (conventional or unconventional), where the technology could be used to try and predict more accurately the formation of reservoirs and/or their behaviour during production; Geothermal energy, to simulate underground circulation of fluids advecting heat; Safe waste disposal, to model the movement of potentially contaminated fluids around waste buried in fractured or faulted geological environments; and Groundwater contamination, to capture more accurately the plume propagation of contaminants underground.

[00167] The skilled person will also appreciate that embodiments of the present invention described above provide advantages over existing methods including:

• Stratigraphic Forward Modelling, a type of process-based modelling used in our workflow/method, usually only includes major elements of the structural framework that can be accounted for in the simulation of basin generation, with corresponding limitations about the structural framework itself (e.g. vertical faults only, smaller faults ignored) and without the ability to set material properties within the generated geological discontinuities, even when those features are thick enough (e.g., fault zone) to be represented in the corresponding mesh.

• 3D geological models generally use geostatistical distributions of rock properties, i.e. not accounting for the physics of the processes responsible for heterogeneous property distributions. This can result in geologically unrealistic property distributions.

• Permeability information within geological faults in fluid flow models is usually set per fault as a constant, sometimes as a function of depth, but rarely as a function of the material properties of the corresponding host rock and rarely with anisotropy following the orientation of the fault.

• In many geological fluid flow simulators, the longitudinal and transverse permeability directions are assumed to be parallel and perpendicular to the local coordinate axes of the mesh; that is, anisotropy must follow layers in the mesh. It is difficult to construct such meshes for complex 3D architectures, especially those involving faults. This problem is resolved by simulators that provide full tensor representation of permeability, as the direction of fluid flow is then independent of the mesh, however this functionality is not commonly used because it is difficult to determine the orientation of the permeability tensor throughout the modelled domain.

[00168] The tools presented here can be applied in any field in which fluid flow in faults and their sedimentary host rocks is of interest. Stratigraphic forward modelling is useful for predicting fault properties where outcrop data are lacking and drillholes are sparse, as is often the case in mineral exploration. In hydrocarbon exploration there may be sufficient drillhole data to constrain sedimentary sequences without an SFM, however the lower dimensional element approach is useful for exploring the impact of faults that have varying characteristics along their length due to varying sediment properties. This approach may also lead to improvements in the handling of faults in groundwater studies and models of contaminant transport or CO2 sequestration.

References

[00169] The following is a list of references, the contents of which are incorporated herein by way of cross-reference.

[00170] Blaikie, T.N., Kunzmann, M., 2020. Geophysical interpretation and tectonic synthesis of the Proterozoic southern McArthur Basin, northern Australia. Precambrian Res. 343.

[00171 ] Kunzmann, M., Crombez, V., Blaikie, T.N., Catuneanu, O., King, R., Halverson, G.P., Schmid, S., Spinks, S.C., 2022. Sequence stratigraphy of the ca. 1640 Ma Barney Creek Formation, McArthur Basin, Australia. Aust. J. Earth Sci.

[00172] Kunzmann, M., Schmid, S., Blaikie, T.N., Halverson, G.P., 2019. Facies analysis, sequence stratigraphy, and carbon isotope chemostratigraphy of a classic Zn-Pb host succession: the Proterozoic middle McArthur Group, McArthur Basin, Australia. Ore Geol. Rev.

[00173] Poulet, T., Lesueur, M., Kelka, II., 2021. Dynamic modelling of overprinted low- permeability fault cores and surrounding damage zones as lower dimensional interfaces for multiphysics simulations. Comput. Geosci. 150, 104719.

[00174] Ehrenberg, S.N., Nadeau, P.H., 2005. Sandstone vs. carbonate petroleum reservoirs: A global perspective on porosity-depth and porosity-permeability relationships. Am. Assoc. Pet. Geol. Bull. 89, 435-445. https://doi.Org/10.1306/11230404071 .

[00175] de Marsily, G., Delay, F., Gongalves, J., Renard, P., Teles, V., Violette, S., 2005.

Dealing with spatial heterogeneity. Hydrogeol. J. 13, 161-183. hit s; /dgi.prg/10. 004-0432-3. [00176] Renard, P., de Marsily, G., 1997. Calculating equivalent permeability: A review. Adv. Water Resour. 20, 253-278.

[00177] Yang, Y., Aplin, A.C., 2010. A permeability-porosity relationship for mudstones. Mar. Pet. Geol. 27, 1692-1697. https://doi.Org/10.1016/J.MARPETGE0.2009.07.001 .

[00178] Matheron, G. (1967). Elements pour une theorie des milieux poreux. France: Eds. Masson et Cie.,

Interpretation

[00179] 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”, “analyzing” or the like, refer to the action and/or processes of a computer or computing system, or similar electronic computing device, that manipulate and/or transform data represented as physical, such as electronic, quantities into other data similarly represented as physical quantities.

[00180] In a similar manner, the term “controller” or "processor" may refer to any device or portion of a device that processes electronic data, e.g., from registers and/or memory to transform that electronic data into other electronic data that, e.g., may be stored in registers and/or memory. A “computer” or a “computing machine” or a "computing platform" may include one or more processors.

[00181 ] Reference throughout this specification to “one embodiment”, “some embodiments” or “an embodiment” means that a particular feature, structure or characteristic described in connection with the embodiment is included in at least one embodiment of the present disclosure. Thus, appearances of the phrases “in one embodiment”, “in some embodiments” or “in an embodiment” in various places throughout this specification are not necessarily all referring to the same embodiment. Furthermore, the particular features, structures or characteristics may be combined in any suitable manner, as would be apparent to one of ordinary skill in the art from this disclosure, in one or more embodiments.

[00182] As used herein, unless otherwise specified the use of the ordinal adjectives "first", "second", "third", etc., to describe a common object, merely indicate that different instances of like objects are being referred to, and are not intended to imply that the objects so described must be in a given sequence, either temporally, spatially, in ranking, or in any other manner. [00183] In the claims below and the description herein, any one of the terms comprising, comprised of or which comprises is an open term that means including at least the elements/features that follow, but not excluding others. Thus, the term comprising, when used in the claims, should not be interpreted as being limitative to the means or elements or steps listed thereafter. For example, the scope of the expression a device comprising A and B should not be limited to devices consisting only of elements A and B. Any one of the terms including or which includes or that includes as used herein is also an open term that also means including at least the elements/features that follow the term, but not excluding others. Thus, including is synonymous with and means comprising.

[00184] It should be appreciated that in the above description of exemplary embodiments of the disclosure, various features of the disclosure are sometimes grouped together in a single embodiment, Fig., or description thereof for the purpose of streamlining the disclosure and aiding in the understanding of one or more of the various inventive aspects. This method of disclosure, however, is not to be interpreted as reflecting an intention that the claims require more features than are expressly recited in each claim. Rather, as the following claims reflect, inventive aspects lie in less than all features of a single foregoing disclosed embodiment. Thus, the claims following the Detailed Description are hereby expressly incorporated into this Detailed Description, with each claim standing on its own as a separate embodiment of this disclosure.

[00185] Furthermore, while some embodiments described herein include some but not other features included in other embodiments, combinations of features of different embodiments are meant to be within the scope of the disclosure, and form different embodiments, as would be understood by those skilled in the art. For example, in the following claims, any of the claimed embodiments can be used in any combination.

[00186] In the description provided herein, numerous specific details are set forth. However, it is understood that embodiments of the disclosure may be practiced without these specific details. In other instances, well-known methods, structures and techniques have not been shown in detail in order not to obscure an understanding of this description.

[00187] Similarly, it is to be noticed that the term coupled, when used in the claims, should not be interpreted as being limited to direct connections only. The terms "coupled" and "connected", along with their derivatives, may be used. It should be understood that these terms are not intended as synonyms for each other. Thus, the scope of the expression a device A coupled to a device B should not be limited to devices or systems wherein an output of device A is directly connected to an input of device B. It means that there exists a path between an output of A and an input of B which may be a path including other devices or means. "Coupled" may mean that two or more elements are either in direct physical, electrical or optical contact, or that two or more elements are not in direct contact with each other but yet still co-operate or interact with each other.

[00188] Embodiments described herein are intended to cover any adaptations or variations of the present invention. Although the present invention has been described and explained in terms of particular exemplary embodiments, one skilled in the art will realize that additional embodiments can be readily envisioned that are within the scope of the present invention.