Method and system for evaluating a storage of fluid in an aquifer structure

Abstract

The present invention operates within the technological domain of subsurface fluid storage resource evaluation, particularly focused on aquifers. The (computer-implemented) methods and systems described in the present disclosure can be applied for evaluating the storage of a fluid in an aquifer structure in diverse industries and applications, such as environmental engineering, energy production, natural resource management, geoscience research, and others. An aspect of the present invention relates to a method for evaluating a storage of fluid in an aquifer structure. In broad terms, the method comprises two main components: firstly, acquiring a geological model of the aquifer structure, the geological model defining a plurality of cells, each cell being provided with properties that define the fluid behavior as a function of spatial position within the geological model; and secondly, setting up a dynamic pressure disturbance numerical solver, which evaluates a propagating pressure field based on various simulation parameters representing an injection of fluid into the aquifer structure. Subsequently, the method outputs data that allows for determining the fluid storage of the aquifer structure using analytical techniques.

Description

[0001] METHOD AND SYSTEM FOR EVALUATING A STORAGE OF FLUID IN AN AQUIFER STRUCTURE

[0002] FIELD OF THE INVENTION

[0003] The present invention operates within the technological domain of subsurface fluid storage resource evaluation, particularly focused on aquifers. The methods and systems described in the present disclosure can be applied for evaluating the storage of a fluid in an aquifer structure in diverse industries and applications, such as environmental engineering, energy production, natural resource management, geoscience research, and others.

[0004] BACKGROUND

[0005] Evaluation of the fluid storage of an aquifer is important for effective water resource management, sustainable development, environmental protection, and resilience to future challenges like climate change and increasing water demand. Within the exploration and assessment of fluid storage resources, such as carbon dioxide (CO2) within saline aquifers, precise evaluation methods are crucial for informed decision-making and resource optimization.

[0006] Presently, two primary workflows are employed in this field. One approach relies on simple analytical pressure modelling, offering speed but often necessitating significant reductions in geological heterogeneity, thus limiting its accuracy. The other approach utilizes conventional geomodelling and multi-physics reservoir simulation workflows based on a virtual reservoir. While this method provides precision, its complexity and time-consuming nature make it less suitable for screening or exploratory ranking purposes.

[0007] Alternatively, analytical models, though simpler and faster, struggle to represent critical geological features such as channel systems or irregular boundary shapes. These models often require extensive human input to accommodate reductions in heterogeneity, leading to inaccuracies in fluid storage resource assessments. Conversely, conventional reservoir simulation workflows offer comprehensive assessments but demand specific skills and extensive computational resources. Setting up flow simulations and ensuring inputs validation and convergence are tedious processes requiring expertise to obtain reliable results.

[0008] Hence, there exists a need for a solution that addresses the limitations of the current state of the art, particularly concerning the challenges faced in achieving efficient and reliable results while reducing human input and computational complexity, particularly in preliminary aquifer screening studies conducted under operational time constraints. SUMMARY OF THE INVENTION

[0009] The present invention operates within the technological domain of subsurface fluid storage resource evaluation, with a particular focus on aquifers. The (computer-implemented) methods and systems described in this disclosure are applicable to the evaluation of fluid storage within aquifer structures across diverse industries and applications, including environmental engineering, energy production, natural resource management, geoscience research, and others.

[0010] One objective of the invention is to streamline the process of fluid storage resource evaluation. As disclosed herein, a workflow is proposed that bridges the gap between simplistic analytical models and more complex 3D gridded full physics models. This workflow aims to enhance the accuracy of evaluations by accommodating reservoir heterogeneities, thus improving pressure forecasting capabilities.

[0011] Importantly, this enhancement can be achieved without significantly increasing computational complexity and / or relying excessively on human input, thereby enhancing overall productivity. This workflow is therefore particularly advantageous in the context of preliminary aquifer screening and exploration studies conducted under operational time constraints. However, its versatility extends to various other applications within the field.

[0012] Another objective of the invention is to provide an easy-to-setup, fast, and reliable tool for fluid storage resource evaluation from multiple realizations of geological models. As disclosed herein, the above- mentioned workflow can be integrated into a tool designed for fast and autonomous operation, free from the constraints typically associated with reservoir simulation, such as convergence checks. Unlike traditional tools that rely on strict development concepts, such as the number and positioning of wells, this tool focuses on the concept of pressurizing an area of the aquifer. By extending the analytical view to accommodate complex geometries and heterogeneities within aquifer structures, it facilitates the evaluation of the maximum resource potential.

[0013] The below summary is provided to introduce a selection of key embodiments of the invention in a simplified form. These embodiments are described in further detail in the detailed description of the disclosure below. This summary is not intended to identify key features or essential features of the claimed subject matter, nor is it intended to be used to limit the scope of the claimed subject matter.

[0014] An aspect of the present invention relates to a (computer-implemented) method for evaluating a storage of fluid in an aquifer structure, the method comprising the steps of:

[0015] - acquiring a geological model representing the aquifer structure; the geological model defining a plurality of cells, each cell being provided with properties that define the fluid behavior as a function of spatial position within the geological model;

[0016] - selecting at least one cell forming a source for fluid injection into the geological model; - implementing a numerical simulation describing an injection of predetermined amount of fluid into the geological model; wherein the numerical simulation is configured to initiate an outward propagation of an incident pressure front from the source across the cells of the geological model until a stop condition is met, thereby assigning a set of incident arrival times to each cell as the incident pressure front propagates outward; wherein the numerical simulation is further configured to identify at least one zero-flow boundary where fluid behavior is substantially prevented, prompting the numerical simulation to discontinue the propagation of the incident pressure front, and / or at least one transmissive boundary where fluid behavior is substantially altered based on a predetermined threshold, prompting the numerical simulation to adapt the outward propagation of the incident pressure front to the altered fluid behavior, wherein the numerical simulation is further configured to initiate an inward propagation of a reflected pressure front from at least one reflection source at the zero-flow boundary and / or transmissive boundary across the cells until a stop condition is met, thereby assigning a set of reflected arrival times to each cell as the reflected pressure front propagates inward;

[0017] - determining a set of pressure amplitudes for each group of cells with substantially similar arrival times by relating the amount of injected fluid to the assigned arrival times; and

[0018] - determining a fluid storage of the aquifer from the sets of pressure amplitudes and arrival times.

[0019] In some embodiments, the method comprises the steps of receiving as input a set of geological model parameters indicative of the behavior of fluid in the aquifer structure; and partitioning the geological model into a plurality of discrete cells, and attributing each cell with properties derived from the set of geological model parameters that define the fluid behavior as a function of spatial position within the geological model.

[0020] In some embodiments the attribution of cell properties comprises determining a Mobility M that quantifies the flow of fluid through the cell under pressure, a Storativity S that quantifies the ability of the cell to store fluid under changes to the pressure; and a Diffusivity D that quantifies the rate at which the injected fluid propagates pressure through the cell, whereby the Diffusivity D is determined from a ratio between the Mobility and the Storativity, D=M / S.

[0021] In some embodiments the numerical simulation is configured to identify the transmissive boundary based on a Mobility and / or Diffusivity ratio between an inward cell and an outward cell along the trajectory of the propagating pressure front, and further adapt the outward propagation of the incident pressure front based on said Mobility and / or Diffusivity ratio. In some embodiments the numerical simulation is configured to group cells of a transmissive boundary and / or a zero-flow boundary with similar incident arrival times, identify a reflection source from the group of cells, and initiate a common reflected pressure front from the reflection source.

[0022] In some embodiments the reflected pressure front is initiated with a substantially similar amplitude to the amplitude of the incident pressure front.

[0023] In some embodiments, the method comprises the steps of receiving as input a set of simulation parameters; and implementing a numerical simulation describing an injection of an amount of fluid into the geological model based on the set of simulation parameters.

[0024] In some embodiments the set of simulation parameters comprises at least a fluid injection zone selected from at least one cell or a subset of cells capable of forming the source, a fluid injection time defined as the duration over which fluid is injected, and / or a fluid injection target representing the maximum pressure range at the conclusion of the fluid injection.

[0025] In some embodiments the set of simulation parameters comprises at least a fluid injection zone comprising a plurality of cells, wherein the method further comprises determining the source by applying a numerical simulation is configured to initiate an inward propagation of an incident pressure front from a plurality of positions selected along the perimeter of the fluid injection zone, thereby assigning a set of incident arrival times to each cell as the pressure front propagates inward, and selecting the source from one or more cells with substantially similar arrival times from each pressure front.

[0026] In some embodiments the wherein the amplitude of the propagating pressure front is adjusted based on the ratio of diffusivity between the cells through which the propagating pressure progresses upon encountering the transmissive boundary.

[0027] In some embodiments solving the numerical simulation comprises selecting a numerical solver configured to solve the governing equations for each cell as the pressure front progresses across the cells of the geological model; wherein the selected numerical solver comprises the Fast-Marching Method (FMM) configured to solve the Eikonal equation as the governing equation.

[0028] In some embodiments solving the geological model is segmented into a plurality of zones demarcated by the one or more boundaries; and wherein the numerical model numerically determines a connection between the plurality of zones of the geological model as the numerical simulation progresses through each zone.

[0029] In some embodiments the pressure amplitude determination comprises grouping cells of substantially similar incident and / or reflected arrival times, and relating the pressure amplitude to the group of cells to the volume of injected fluid through the set of geological model parameters including the porosity and compressibility. In some embodiments the propagating pressure front is configured to continue progressing outward from the injection point with a variable amplitude until a stop condition is met, wherein the amplitude of the propagating pressure front is dependent upon the properties of the cells of the geological model.

[0030] In some embodiments the fluid amount determination comprises the step of relating the mass of injected fluid to its volume through the set of geological model parameters including the temperature and density of injected fluid.

[0031] In some embodiments the fluid comprises water and / or air, preferably CO2.

[0032] In some embodiments the method comprises the step of generating of output data, the output data comprising at least one selected from an arrival times map depicting the arrival times of the propagating pressure front at determined locations within the aquifer structure; a pressure map depicting pressure variations at determined locations within the aquifer structure; and / or a histogram depicting the utilization of resources within the geological model during the numerical simulation, said resources preferably comprising at least one of fluid injection rate, fluid injection time, changes in pressure, changes in saturation level, storage capacity, trapping mechanics.

[0033] DESCRIPTION OF THE FIGURES

[0034] The following description of the figures relate to specific embodiments of the disclosure which are merely exemplary in nature and not intended to limit the present teachings, their application or uses.

[0035] Throughout the drawings, the corresponding reference numerals indicate the following parts and features: aquifer 10; geological model 20; fluid boundary 25; injection zone 30; source 34; incident pressure front 35; transmissive boundary 40; transmitted pressure front 45; zero-flow boundary 50; reflected pressure front 55; isochrones 60-65 of the pressure front.

[0036] FIG. 1 is a flowchart illustrating a method 100 for evaluating a storage of fluid in an aquifer structure in accordance with certain embodiments of the present disclosure.

[0037] FIG. 2 shows a geological model 20 of an exemplary aquifer with an injection zone 30 in accordance with certain embodiments of the present disclosure.

[0038] FIG. 3 schematically illustrates the diffusivity as a function of spatial position in the aquifer 10 represented by the geological model 20 shown in FIG. 2.

[0039] FIG.4 depicts the geological model 20 previously discussed of FIG. 2 with the presence of a fluid boundary 25. The fluid boundary delineates a zone of the geological model 20 where fluid behavior is altered as determined from the diffusivity as illustrated in FIG. 3, in accordance with certain embodiments of the present disclosure. FIG. 5 illustrates the identification of a plurality of fluid boundaries 25 based on areas with contrasting diffusivity values in accordance with certain embodiments of the present disclosure.

[0040] FIG.6 illustrates the response of a propagation pressure front to a transmissive boundary 40 in accordance with certain embodiments of the present disclosure.

[0041] FIG. 7 illustrates the response of a propagation pressure front to a zero-flow boundary 50 in accordance with certain embodiments of the present disclosure.

[0042] FIG. 8 illustrates the determination of a source 34 within a fluid injection zone 30 that can function as source for the simulation in accordance with certain embodiments of the present disclosure.

[0043] FIG. 9 illustrates the determination of a plurality of potential sources 34 within a fluid injection zone 30 from which at least one can be selected as source for the simulation in accordance with certain embodiments of the present disclosure.

[0044] FIG. 10 is a schematic representation of a system 500 configured for assessing fluid storage in an aquifer structure in accordance with certain embodiments of the present disclosure.

[0045] FIG. 11 illustrates the propagation of a pressure front within a geological model 20 without the presence of a zero-flow boundary or a transmissive boundary as described herein.

[0046] FIG. 12 presents a visualization of Figure 11 as an equivalent radial system.

[0047] FIG. 13 is a pressure distribution graph derived from Figure 11 showcasing the pressure amplitudes AP (bars) relative to the distance travelled (in meters) along the normal to isochrones (years). Further, the correspondence between the investigated Radius (R) and the isochrones (T) is presented as a dotted line (R-T).

[0048] FIG. 14 illustrates a propagating pressure front in a geological model 20 with a zero-flow boundary 50 as described herein.

[0049] FIG. 15 is representation of Figure 14 as an equivalent radial system illustrating the presence of a reflected pressure front 55 in accordance with certain embodiments of the present disclosure.

[0050] FIG. 16 is a pressure distribution graph of Figure 14 with the pressure amplitudes AP (bars) relative to the distance travelled (in meters) along the normal to isochrones (years). Further, the correspondence between the investigated Radius (R) and the isochrones (T) is presented as a dotted line (R-T).

[0051] FIG. 17 illustrates the propagation of a pressure front within a geological model 20 with the presence of a transmissive boundary 40 as described herein.

[0052] FIG. 18 presents a visualization of Figure 17 as an equivalent radial system illustrating the presence of a transmitted pressure front 45 and a reflected pressure front 55.

[0053] FIG. 19 is a pressure distribution graph derived from Figure 17 showcasing the pressure amplitudes AP (bars) relative to the distance travelled (in meters) along the normal to isochrones (years). Further, the correspondence between the investigated Radius (R) and the isochrones (T) is presented as a dotted line (R-T).

[0054] FIG. 20 is an example of a contrasting mobility M ratio between an inward cell (Min) and an outward cell (Mout) along a transmissive boundary in accordance with certain embodiments of the present disclosure.

[0055] FIG. 21 is an example of a determination of a Volume V corresponding to an isochrone T of a transmitted segment (VpT) and a reflected segment (VpR) of the pressure front in accordance with certain embodiments of the present disclosure.

[0056] FIG. 22 is an illustration of elementary pressure amplitude profile AP (bars) relative to isochrones (years) along the radius of investigation (meters) from the source, to be used in the calculation of a pressure profile in accordance with certain embodiments of the present disclosure.

[0057] FIG. 23 presents a visualization of the arrival times for an incident propagating pressure front (IPF) within a geological model in accordance with certain embodiments of the present disclosure.

[0058] FIG. 24 presents a visualization of the arrival times for a first reflected pressure front (RPF) originating from the zero-flow boundary labelled as (A) in Figure 23 in accordance with certain embodiments of the present disclosure.

[0059] FIG. 25 presents a visualization of the arrival times for a second reflected pressure front (RPF) originating from the zero-flow boundary labelled as (B) in Figure 23 in accordance with certain embodiments of the present disclosure.

[0060] FIG. 26 presents a visualization of the arrival times for a third reflected pressure front (RPF) originating from the zero-flow boundary labelled as (C) in Figure 23 in accordance with certain embodiments of the present disclosure.

[0061] Figure 27 is a pressure distribution map exhibiting the spatial distribution of pressure amplitudes AP (bars) within a geological model in accordance with certain embodiments of the present disclosure.

[0062] Figure 28 is a histogram representing fluid storage resources (metric tons) within an aquifer as measured in metric tons (MT) in accordance with certain embodiments of the present disclosure.

[0063] DETAILED DESCRIPTION

[0064] Although certain embodiments and examples are disclosed below, it will be understood by those in the art that the present disclosure extends beyond the specifically disclosed embodiments and / or uses of the present disclosure and obvious modifications and equivalents thereof. Thus, it is intended that the scope of the present disclosure should not be limited by the particular disclosed embodiments described below. In the present disclosure, technology is described that relates to the evaluation of fluid storage in a subterranean geological reservoir, specifically, an aquifer. The technology described herein aims to address the limitations of the current state of the art in efficiently evaluating the total volume or capacity of fluid that the aquifer can hold within its geological formation.

[0065] Evaluating the aquifer storage, in accordance with aspects of the present disclosure, involves analysing various geological parameters and features to obtain a comprehensive geological model representing the aquifer. These geological parameters may provide information about fluid flow behaviour, aquifer characteristics, and hydrogeological processes.

[0066] As used herein, "geological model” refers to a computer-generated, virtual representation of the aquifer and optionally any number of substructures forming part of the aquifer structure. Advantageously, this virtual representation is primarily three-dimensional (3D) in nature, although a two-dimensional (2D) model is also contemplated, with the option to project a 2D model from a 3D base for computational efficiency. The numerical model may take various forms, including a volume mesh, a polygon mesh with implicit volume definition, or a collection of surfaces and curves outlining the boundaries of the aquifer structure.

[0067] In specific embodiments, the numerical model may adopt a parametric approach, utilizing parametric curves and surfaces or parametrized mesh deformations to define the shape boundary of the aquifer. Importantly, the model integrates geological model data associated with the volume of the virtually represented aquifer, incorporating information about materials or substructures within the aquifer materials. This data may be embedded within the geological model, appended as metadata, or organized as a separate structure indexing the modelled components.

[0068] Consequently, "geological model parameters" refers to characteristics or properties used to define the geological model of the aquifer structure. These parameters allow for accurately representing the behaviour of fluids within the aquifer structure. Some non-exhaustive exemplary parameters may include permeability (hydraulic conductivity), porosity, lithology, fluid, and / or various structure-fluid parameters, such as pressure, temperature, water density, and pore compressibility. Adjusting the values of the geological model parameters appropriately can improve the accuracy of the geological model in representing the aquifer structure.

[0069] As used herein, "numerical simulation" refers to a computational technique used to model the injection of fluid into the geological model as described herein in order to predict the resulting pressure changes over time. The numerical simulation generates a dynamic pressure disturbance profile, offering valuable information on the spatiotemporal evolution of pressure changes within the geological model due to fluid injection. The simulation involves implementing one or more mathematical algorithms to solve the governing equations that describe the pressure propagation within the geological model. Consequently, "numerical simulation parameters" refers to the variables and settings used to configure and control the behaviour of the numerical simulation as described herein. These parameters define the conditions under which the simulation operates and influence its accuracy, efficiency, and reliability. Some non-exhaustive exemplary parameters may include various initial conditions defining the starting state, various solver settings used to solve the governing equations such as convergence criteria, solution tolerance, and iterative methods, and / or various fluid properties, such as fluid density, viscosity, and compressibility. Adjusting the values of the numerical simulation parameters appropriately can improve the reliability and accuracy of output data from the simulation.

[0070] In the following section, the technology underlying the present disclosure will be described by means of different aspects thereof. It will be readily understood that the aspects of the present disclosure, as generally described herein, and illustrated in the figures, can be arranged, substituted, combined, and designed in a wide variety of different configurations, all of which are explicitly contemplated and make part of this disclosure. This description is meant to aid the reader in understanding the technological concepts more easily, but it is not meant to limit the scope of the present disclosure, which is limited only by the claims. Hence, the description below is to be regarded as illustrative in nature, and not as restrictive.

[0071] An aspect of the present invention relates to a method for evaluating a storage of fluid in an aquifer structure. In broad terms, the method comprises two main components: firstly, setting up a geological model of the aquifer structure as described herein, including a gridded model structure incorporating various geological model parameters that define the behaviour of fluid in the aquifer structure; and secondly, setting up a dynamic pressure disturbance numerical solver as described herein, which evaluates a propagating pressure field based on various simulation parameters representing an injection of fluid into the aquifer structure. Subsequently, the method outputs data that allows for determining the fluid storage of the aquifer structure using known techniques.

[0072] The present method and embodiment thereof will be discussed further with reference to FIG 1, which shows a flow diagram of the fluid storage evaluation method 100 as described herein. It is understood that additional steps can be provided before, during, and after the steps of the method, and that some of the steps described can be replaced or eliminated for other implementations of the method. For example, although not shown, the method may be preceded by one or more steps for determining the quantitative values of one or more input parameters to the improve modelling accuracy.

[0073] In the shown embodiment, the method 100 begins at step 101 by receiving as input a set of geological model parameters indicative of the behaviour of fluid in the aquifer structure. As previously described, these geological model parameters can serve several key functions relevant to characterization of aquifer properties and spatial distribution of the aquifer structure for subsequent modelling and analysis. These geological model parameters may be acquired through various methods known in the art that determine or estimate the quantitative values of the relevant geological model parameters. For example, the input can include seismic or geologic models containing a description of the aquifer's subsurface geology.

[0074] In certain embodiments the geological model parameters may comprise at least one of a spatial map (such as a top-base map) depicting the spatial extent and physical boundaries of the aquifer, a thickness map depicting the thickness variations of the aquifer across an area, a permeability map indicating the spatial distribution of permeability values within the aquifer, a porosity map showing the spatial distribution of porosity values within the aquifer, a materials or substructures map (such as a net to gross map) depicting the ratio of net reservoir rock thickness to the total thickness of the aquifer, a thermal / temperature map indicating the spatial variation of temperature within the aquifer, initial pressure map illustrating the initial pressure conditions within the aquifer (before any injection activity), pore compressibility map displaying the spatial distribution of compressibility values within the aquifer, pressure distribution map (such as a cap rock frac pressure map) maps indicating the fracture pressure distribution within one or more sealing layers (cap rock) of the aquifer. It may be appreciated that advantageously the input can comprise combinations of various maps and data types to improve the modelling accuracy of the aquifer. In certain embodiments the geological model parameters may comprise one or more parameters describing the viscosity and compressibility of the fluid as functions of pressure and temperature variations. These properties of the fluid enable modelling the fluid behaviour within the aquifer. Further parameters may include information about the fluid, such as salinity, chemical composition, presence of additives or contaminants capable of altering the viscosity and compressibility, and the like. For example, the input may comprise tables, graphs, or mathematical equations that correlate water viscosity and compressibility with changes in pressure and temperature. These parameters allow for more accurately estimating the fluid properties at specific pressure and temperature conditions encountered within the aquifer.

[0075] In a specific embodiment, the geological model parameters may comprise at least a permeability, a porosity, a viscosity, and a compressibility of the aquifer structure and / or a substructure of the aquifer. As used herein, permeability determines how easily fluids can flow through the aquifer, while porosity quantifies the volume of pore space available for fluid storage. Compressibility is a measure of how much the aquifer materials can be compressed or compacted under pressure to undergo volume changes in response to changes in fluid pressure.

[0076] Returning to Figure 1, at step 102, the method 100 may include a step of acquiring at least one geological model, as a virtual representation of the aquifer structure, to initiate the fluid storage evaluation process. As previously described, a geological model may be generated based on the set of geological model parameters provided at step 101.

[0077] Figure 2 depicts an example of a geological model 20 representing an exemplary aquifer. This geological model has a surface area defined by the perimeter of the aquifer. Additionally, the presence of a fluid injection zone 30 is indicated, from which the injection of a fluid into the geological model 20 will be simulated as described in later steps of the method 100.

[0078] Figure 3 depicts a visualization of a geological model parameter as a function of spatial position in an aquifer 10 represented by the geological model 20 shown in figure 2 Specifically, it shows a labelled map containing information about the diffusivity (m2 / s) distribution within the aquifer 10 as defined herein. The color coding (in grayscale) corresponds to specific values of the diffusivity. Alternatively or in combination, other geological model parameters may be included in the same or complementary maps, such as the permeability, porosity, water viscosity, and pore compressibility of the aquifer.

[0079] In certain embodiments, the numerical model is generated at least partially from 'a priori' knowledge of the aquifer structure. The generation of the geological model from a priori knowledge implies that this model is generated based on predefined specifications or known characteristics of the aquifer in question. This a-priori knowledge can originate from an estimation of similar or standardised aquifer structures, including a number of anticipated features and structural details of the aquifer without relying on any actual input data. For instance, certain parameters relevant to porosity may be derived from information regarding the expected subsurface materials contained within the aquifer.

[0080] In some embodiment the method may include a step of partitioning the geological model into a plurality of grid cells, thereby obtaining a gridded geological model. The grid generation can be implemented using modelling software or computational tools configured to divide the geological model into smaller, advantageously discrete, components of a selected geometric shape, each representing a portion of the aquifer. In exemplary embodiments the selected shape may include cubes, rectangular prisms, hexahedral prisms, tetrahedrons pyramids, prisms, wedges, and other shapes. The choice of grid cell shape depends on factors such as the complexity of the geological formation, computational efficiency requirements, and the accuracy needed for the simulation.

[0081] In a specific embodiment, the gridded geological model can be generated by delineating the upper and lower surface boundaries of the geological model, which preferably correspond to the upper and lower surfaces of the aquifer, and dividing the volume between these boundaries into grid cells. The surfaces representing the top and base boundaries of the aquifer can be delineated using interpolation techniques or by digitizing geological interpretations from available data sources. Optionally, further partitioning can be performed along the height of the geological model, corresponding to the distance between the upper and lower surface boundaries, if relevant for the simulation. Advantageously, this gridded geological model spans the entirety of the aquifer's extent, ensuring adequate resolution to capture geological features accurately.

[0082] In a specific embodiment, the method may further include one or more refinement steps aimed at optimising the grid resolution. These refinement can be applied for instance in areas of interest such as near faults or regions exhibiting significant geological heterogeneity. Additionally, mesh optimization methods can be employed to ensure that the grid conforms well to the geological features and boundary conditions. For instance, these techniques may include smoothing or alignment algorithms. Other refinement criteria may be included to guide the refinement process to areas where it will have the most significant impact on simulation accuracy. This refinement process can enhance the representation of fine-scale features, thereby augmenting the model's accuracy.

[0083] In some embodiments the method may further include a step of attributing each cell with properties derived from the set of geological model parameters to define the fluid behaviour as a function of spatial position within the geological model. Each grid cell may contain information about the corresponding portion's physical or chemical properties, such as density, porosity, lithologic facies, and so on, as previously described. The combination of information assigned to each will determine the amplitude of propagating pressure front in the simulation described below. The amplitude is defined to be a proxy to represent how the pressure front propagates through the model. In certain embodiments multiple functions and modifiers can be defined, enabling more complex geology and geometries to be incorporated into this technique of quantifying the connective quality.

[0084] In an embodiment, the attribution of cell properties comprises determining a Storativity S based on S = <pct, wherein <p is the porosity of the aquifer, which represents the volume of void space in the rock or sediment and indicates how much water it can hold, and ctis the total compressibility of the system corresponding to a summation of the water compressibility in the aquifer and pore pressure compressibility, which measures how much the volume of water stored in the aquifer changes with changes in hydraulic head. Consequently, a higher porosity or compressibility will result in a higher storativity, indicating that the aquifer can store more water relative to changes in hydraulic head. Conversely, lower values of porosity or compressibility will result in lower storativity.

[0085] In an embodiment, the attribution of cell properties comprises determining a Mobility M based on M = K / / i, wherein K is the hydraulic conductivity of the porous medium, which represents its ability to transmit fluids, measured in units like meters per second (m / s) or centimeters per second (cm / s), and / i is the dynamic viscosity of the fluid, which measures its resistance to flow, measured in units like pascal- seconds (Pa-s) or centipoise (cP). Consequently, a higher hydraulic conductivity or lower dynamic viscosity will result in a higher mobility, indicating that the fluid can move more easily through the porous medium. Conversely, lower hydraulic conductivity or higher dynamic viscosity will result in lower mobility, indicating more resistance to fluid flow through the medium.

[0086] In an embodiment, the attribution of cell properties comprises determining a Diffusivity D based on D = M / S, which is the ratio of the above defined Mobility and Storativity. Consequently, a higher mobility or lower storativity will result in a higher diffusivity, indicating that the pressure front can spread more rapidly through the porous medium. Conversely, lower mobility or higher storativity will result in lower diffusivity, indicating slower spreading or diffusion of the pressure.

[0087] Figure 4 depicts the geological model 20 previously discussed in Figure 2, with the presence of a fluid boundary 25 indicated by the dashed line. This fluid boundary represents an interface with a significant contrast in the behaviour of the fluid, specifically in terms of diffusivity (m2 / s) as illustrated in Figure 3. In this example, the fluid boundary delineates two fluid zones where the cells of the geological model 20 have substantially similar properties. Consequently, the behaviour of fluid passing through the fluid boundary 25, as indicated by the direction of the arrow, will be significantly altered.

[0088] In an alternative to any of the above embodiments, the method may include utilizing of a pre-generated geological model comprising a plurality of cells provided with properties that define the fluid behavior as a function of spatial position within the geological model. This approach allows for the subsequent execution of the method steps on a pre-existing geological model. The preceding steps primarily serve to provide guidance to the skilled person in preparing the requisite groundwork for performing the simulation as described below.

[0089] The cell properties as defined herein can be used to detect fluid boundaries in the geological model representing regions or sections of the aquifer where the behavior of the fluid is substantially altered. The boundaries can be identified by analyzing the previously discussed parameters, such as permeability, porosity, lithology, and fluid properties, which determine the behavior of fluid within the aquifer structure.

[0090] In an embodiment, the identification of fluid boundaries comprises comparing the cell properties of adjacent, preferably neighboring, cells. Alternatively, a plurality of cells having substantially similar cell properties can be grouped together into fluid zone, and the identification of fluid boundaries comparing the cell properties of neighboring or adjacent fluid zones. By analyzing the attributes across the grid cells, patterns and contrasts in properties can be identified, indicating potential zones where fluid behavior may be altered. For example, abrupt changes in permeability values or lithological properties may signify the presence of faults or lithological boundaries that can influence fluid flow patterns. In an embodiment, identification of fluid boundaries comprises comparing the above defined Storativity S, Mobility M and Diffusivity D determined for each cell and / or fluid zones. These parameters can improve the identification of regions or clusters where fluid behavior deviates from surrounding areas. These regions may correspond to areas of hydraulic contrasts or preferential flow paths within the geological model.

[0091] In the above embodiments, identification of fluid boundaries may further comprise identifying a plurality of adjacent cells with substantially different values of diffusivity D as described herein; and designating the boundary as a transmissive boundary and / or zero-flow boundary based on one or more predetermined thresholds. For example, the identified cell can be designated as a transmissive boundary if the contrasting diffusivity D is above a predetermined threshold, and as a zero-flow boundary if the contrasting diffusivity D is below the predetermined threshold.

[0092] Figure 5 illustrates an example of an aquifer comprising a plurality of zones with significantly altered fluid behaviour. In present example, diffusivity is calculated from the provided values of mobility M and storativity S. When reviewing the diffusivity along direction of the arrow differences can be identified based on the delineated zones. Specifically, it is shown that the diffusivity starts with a value of 3.0 m2 / s, but as the first fluid boundary 25 is crossed, the diffusivity increases from 3.0 m2 / s to 5.0 m2 / s, corresponding to a zone with reduced fluid flow resistance indicative of a lower amplitude pressure increase. Subsequently, when crossing through a second fluid boundary 25', the diffusivity decreases from 5.0 m2 / s to 1.0 m2 / s, indicative of a zone with increased fluid flow resistance indicative of a higher amplitude pressure increase. Finally, a third fluid boundary 25" is reached where passage of fluid is prevented resulting in zero fluid flow. It should be noted that the provided diffusivity values are exemplary for illustrating the classification of fluid boundaries as described herein. In real-world applications, the contrast difference may be dynamically adjusted based on the desired accuracy and computational complexity of the model.

[0093] As will be later described, at least one identified fluid boundary can be designated as a transmissive boundary characterized in that the fluid movement is permitted to a significant extent across the fluid boundary. This type of boundary condition typically represents areas where fluid flow is not restricted by the boundary of the aquifer but represent a change in one or more geological model parameters determining the fluid behaviour. For example, the transmissive boundary may correspond to a region of the aquifer that is characterised by a different properties of the subsurface materials, for example, passing from porous material having a relatively high porosity to rocks having a relatively low porosity, resulting into a higher pressure increase for the same amount of fluid passing through, or vice versa. Figure 6 illustrates an example of such a transmissive boundaries 40, 40' that are identified within the fluid zone depicted in the previously discussed Figure 5. Additionally, the presence of a fluid injection zone 30 is indicated, from which the injection of fluid into the geological model 20 can be simulated, as described in subsequent steps of the method 100. As depicted, the fluid injection initiates a propagating incident pressure front that is capable of passing through each one of the transmissive boundaries 40, 40', thereby adjusting the associated parameters with the pressure front, as it continues its propagation in the geological model.

[0094] As will be later described, at least one identified fluid boundary can be designated as a zero-flow boundary characterized by fluid movement being prevented to a significant extent. This type of boundary condition typically represents areas where fluid flow is constrained or blocked, such as the boundary of an aquifer or a geological barrier. For instance, if the aquifer contains impermeable layers on its boundaries such as rocks or other materials with very low to no porosity, these boundaries may be defined as zero-flow boundaries to simulate the lack of flow across them. Consequently, the fluid flow is restricted, and the fluid is expected not to pass across the boundary, or at least not in an amount that would permit a significant flow of fluid.

[0095] Figure 7 illustrates an example of such a zero-flow boundary 50 that is situated along the perimeter of the fluid zone depicted in the previously discussed Figure 5. Additionally, the presence of a fluid injection zone 30 is indicated, from which the injection of fluid into the geological model 20 will be simulated, as described in subsequent steps of the method 100. As depicted, the fluid injection initiates a propagating pressure front 60 that cannot pass through the zero-flow boundary 50, thereby preventing the pressure front 60 from advancing further and halting its propagation in the geological model.

[0096] Returning to Figure 1, at step 103, the method 100 may include a step of receiving as input a set of simulation parameters relevant to setting up a dynamic pressure disturbance simulation in the geological model. As previously described, these geological model parameters serve several key functions relevant to the properties and spatial positioning of the propagating pressure front during subsequent modelling. In an embodiment, the set of simulation parameters comprises at least a fluid injection point or a fluid injection zone consisting of group of spatial points, a fluid injection time, and / or a fluid injection target. As used herein, the fluid injection point or zone determines the location within the geological model where fluid is introduced into the system and can be specified in terms of its coordinates within the model’s coordinate system. The fluid injection time determines the duration or length of time over which fluid will be injected into the geological model. It specifies the period during which the injection process occurs and influences the temporal aspect of the simulation. The fluid injection target is defined as the maximum pressure range that is desired or acceptable at the end of the fluid injection process within the geological model. It serves as a target or threshold for pressure buildup resulting from the fluid injection and helps to guide the simulation towards achieving the desired pressure conditions within the model.

[0097] In certain embodiments, the method may comprise a step of selecting one or more fluid injection points within a designated fluid injection zone received as a simulation parameter. Specifically, as will be later described, the simulation optimally initiates the propagating pressure front at a single point. Consequently, as the propagating pressure front will continue propagating for a predetermined amount of time, the fluid injection point can be denotated as the starting time 0, referred to hereafter To. However, choosing an optimal starting point based solely on available aquifer information can pose challenges due to its complexity.

[0098] Therefore, in a further embodiment, the method may further comprise selecting a fluid zone consisting of a plurality of spatial points from which one or more fluid injection points can be selected. The selection may comprise an inward simulation, whereby a pressure front is initiated at various locations along the perimeter of the selected fluid injection zone. In this simulation the pressure front is configured to continue progressing inward towards a central region of the fluid injection zone. Consequently, the simulation may produce one or more sets of arrival times as a function of spatial position within the fluid injection zone. From these sets, at least one point representing substantially similar arrival times within the fluid injection zone can be designated as the fluid injection point. Optionally, multiple injection points may be selected, which can then be ranked based on predetermined criteria to determine the preferred fluid injection point.

[0099] Figure 8 illustrates an exemplary selection of a single point forming a source 34 within a fluid injection zone 30, which could correspond to the zone 30 highlighted in the previously discussed geological map 20 of Figure 2. In this example, the fluid injection point is determined by running a numerical solver from the edges of the fluid injection zone 30 towards its centre. The numerical solver may for example include FMM as described herein. The point with the highest arrival time, representing the longest time needed to reach the central point, is designated as the source 34 from which the below discussed simulation can be initiated.

[0100] Figure 9 illustrates another example wherein a plurality of points are selected forming a potential source 34 within a fluid injection zone 30. Similar to the previous example, the fluid injection point is determined by running a numerical solver from the edges of the fluid injection zone 30 towards its centre. However, depending on the shape or characteristics of the map, multiple sources 34, 34 can be identified, either corresponding to substantially similar arrival times or by establishing sets of local maximums.4 Advantageously, the plurality of identified injection points can be sorted based on the corresponding arrival times. For instance, in the present example, the first point 34 (labelled as 1) may have an arrival time of 3 years, while the second point 34 (labelled as 2) may have an arrival time of 2.5 years. Therefore, depending on the simulation preferences, a selection can be made either automatically based on a predetermined criteria (such as selecting the point corresponding to the longest arrival time) or based on further user input (such as selecting a more accessible position). Optionally, the below described simulation may be iteratively repeated for each one of the selected fluid injection points to facilitate comparative analysis.

[0101] In an alternative to any of the above embodiments, the method may include utilizing of a set of predetermined simulation parameters. For example, the predetermined simulation parameters may include a list of commonly implemented values based on previous or commonly implemented simulations. The preceding steps primarily serve to provide guidance to the skilled person in implementing the simulation described below.

[0102] Returning to Figure 1, at step 104, the method 100 may further include a step of simulating a dynamic pressure disturbance in the geological model as a function of time. Specifically, the pressure disturbance can be simulated by implementing numerical techniques that describe how the pressure dynamically evolves in response to an injection of a predetermined amount of fluid over a predetermined amount of time, defined by the provided simulation parameters as described herein. In essence, this numerical simulation aims to predict the temporal evolution of pressure changes within the geological model by tracking the changes in pressure magnitude and distribution at different spatial locations within the geological model.

[0103] The propagation of the pressure front can be started from a single cell or a group of cells, referred to herein as the "source", its location depending on the simulation parameters described herein. The amplitude of the propagating front varies across the gridded geological model is determined by the arrival times of the pressure front. Accordingly, as the front propagates the travel times from the source to other cells in the reservoir can be calculated, with each cell possessing its unique travel time. These times, referred to as arrival times, signify the duration taken for the propagating front to reach a specific cell from the source, indicating the connectivity quality of cells in the reservoir. Longer travel times imply lower connectivity quality, while shorter times indicate higher connectivity quality. The arrival times will determine the amplitude of the pressure at connected cell within the geological model.

[0104] The calculation of travel times can be performed through a numerical solver configured to computationally solve the governing equations of each cell of the geological map as the propagating pressure front progresses through the geological map. As time progresses, this numerical solver utilizes algorithms and numerical techniques to solve these equations and obtain numerical solutions approximating the behavior of the injected fluid over time and space. The numerical solver continues computing the arrival times at the cells of the propagating front until meeting termination criteria, such as reaching a maximum arrival time, maximum distance, or the front encountering a target object.

[0105] The governing equations of the cells can be often complex, being nonlinear, of mixed character, and coupled partial differential equations with spatially varying coefficients. Consequently, the numerical solver can be configured to provide an asymptotic approximation of the pressure field to simplify the complex governing equations, as will be described below. This approach is particularly useful for dealing with such complicated equations, allowing for the derivation of simplified solutions or expressions that still capture essential features of the original problem. Nevertheless, the techniques described here are general and may be applied to any system of nonlinear partial differential equations. Thus, one may consider coupled processes involving thermal effects. However, as additional variables and equations are added, the resulting expressions for the relevant parameters become increasingly complicated and may require increased computational efforts.

[0106] In an embodiment, the configuration of the numerical solver involves setting parameters such as the discretization method (e.g., finite difference, finite element, finite volume), spatial and temporal discretization schemes (e.g., grid size, time step), boundary conditions, and convergence criteria. Selection of these parameters determines the accuracy, stability, and efficiency in solving the governing equations numerically based on the available computational resources.

[0107] In an embodiment, the parameters primarily influencing the propagation of the pressure front include the Mobility M, Storativity S and the Diffusivity D as described herein. It is notably understood that parameters influencing the respective values of the Mobility M, Storativity S and the Diffusivity D inherently affect the propagation of the pressure front.

[0108] In some preferred embodiments, the governing equations include a continuity equation that is transformed into an Eikonal equation, which describes the propagation of fronts in a medium with a known velocity distribution, especially distribution of a diffusivity square root. Specifically, the Eikonal equation is a known partial differential equation used to model the propagation of fronts in various fields such as fluid dynamics. It describes the evolution of the front by specifying the speed at which information propagates through a medium.

[0109] In some preferred embodiments, the numerical solver may include the Fast Marching Method (FMM) as a numerical algorithm used for solving the governing equations, advantageously, the Eikonal equation as described herein. The FMM is particularly efficient in computing the arrival times of pressure fronts at discrete grid points in a domain, given the speed of propagation at each point. It does so by iteratively updating the arrival times based on the known arrival times at neighboring points, prioritizing the points with the smallest arrival times to update next. This iterative process continues until all points in the domain have been visited and their arrival times computed.

[0110] Further at step 106 of Figure 1, after completion, the numerical solver may provide as solution a set of arrival times of the propagating pressure front at determined locations in the geological model. These arrival times can be represented as an "isochrone", which herein refers to a line or contour representing the arrival of the pressure front at different points within the model at a specific time. In essence, the isochrones are a graphical representation showing where the pressure front has reached in the geological model at a given moment in time. Each isochrone typically corresponds to a specific time interval, indicating the spatial extent of the pressure disturbance at that particular time. These isochrones may help visualize the propagation of pressure over time and provide insights into how the pressure front evolves within the geological model during the simulation.

[0111] In an embodiment, the isochrones can be primarily treated as a moving no-flow boundary in a steady state. In another embodiment, the isochrones can also be interpreted as a moving radius or distance of investigation of the transient. It is assumed that the variation in diffusivity is smooth enough to directly map the drained volume isochrones onto an equivalent radial system.

[0112] Further still at step 107 of Figure 1, the relevant pressure profiles can be derived from the solution comprising the arrival times or sets of arrival times. Specifically, these pressure profiles can be derived by interpolating or extrapolating the arrival times to estimate the pressure distribution throughout the geological model at different points in time. This process may involve applying mathematical or numerical techniques to relate the arrival times to pressure changes, taking into account factors such as fluid properties, the below discussed boundaries, and the geometry of the geological model.

[0113] In an embodiment, the boundary condition at the edge of the injection is defined as the maximum increase in pressure at the boundary. The injection duration can be defined as the duration of the pressure increase at the boundary. By representing the model as a radial analytic, analytical solutions used in hydrogeology can be implemented to describe the behaviour of fluid flow in porous media. In some embodiments the Theis solution and / or the Thiem solution as known in the art can be considered, or any combination thereof, with the Thiem solution being of particular interest due to reflecting a zero increase of pressure condition at the edge.

[0114] In an embodiment, the numerical solver results in solutions that are calculated along outward trajectories from the source, which can include the injection point of zone, towards a point of interest, which may correspond to an end point or zone, such as the perimeter of the geological model. When using an asymptotic approximation of the pressure field, trajectories refer to the paths followed by pressure changes within the geological model over time. These trajectories represent the evolution of pressure distributions as fluid is injected into the geological model.

[0115] In a further embodiment, pressure mapping the 'radial' pressure field onto the real 'isochrones' enables the determination of the local increase in pressure and thus the volume of water expelled from the injection zone. Inside the injection zone, and to be conservative, the compressible volume is determined by applying the maximum increase in pressure at its boundary to the entire injection zone. The summation of these procedures may provide the total injected volume at downhole conditions, which can later be converted to fluid mass, facilitated by the end of injection condition within the injected area.

[0116] In an embodiment, the evaluation process can initially employ water as reference fluid, which is later converted into mass stored of a fluid of choice. For example, the selection may include CO2 mass stored to allow determination of the CO2 storage resource. This conversion can be justified by the fact that the compressed volume in the injection area is generally lower than the flow outside of the injection area. When the injection area covers the overall aquifer, the fluid storage preferably being CO2 storage resource is directly determined as the allowable compressible volume of water.

[0117] Returning to Figure 1, at step 105, the method may include a step of identifying one or more fluid boundaries characterized by substantial changes to fluid behavior as the pressure front propagates outwards across the cells of the geological model. These fluid boundaries may include a transmissive fluid boundary, where fluid behavior is substantially altered based on predetermined criteria, and / or a reflective fluid boundary, where fluid behavior is substantially prevented. It should be appreciated that the method described herein can be implemented based on the identification of only a transmissive fluid boundary, only a reflective fluid boundary, or a combination of both, depending on the aquifer properties. Examples of these boundaries have been discussed previously. Additionally, the identification of these boundaries may occur during the simulation, or alternatively, may happen in a preceding processing step before the simulation initiates to optimize computational resources.

[0118] In an embodiment, the numerical simulation's propagating pressure front may be configured to discontinue outward progression upon encountering the zero-flow boundary, thereby triggering propagation of a reflected pressure front in the opposite direction, preferably inward back towards the source. The initiation of this reflected pressure front can be adapted based on the desired injection period, determining the 'primary' reflection point from which the reflected pressure front is simulated. Specifically, at each instance when an isochrone intersects a boundary, the cells associated with the boundary from this juncture to the last isochrone can be set to generate a "reflected segment", i.e., a reflected pressure front initiated from the zero-flow boundary that propagates in a direction substantially opposite to the incident pressure front. In further embodiments, a reflected propagating pressure front can be reinitiated from this reflected segment inward within the geological model, ensuring continuity. The time of first arrival from the source is enforced as an initial condition for this new front propagation. Subsequently, the pressure field can be computed similarly to the aforementioned embodiments from the primary injection point, integrating reflections from the boundary while maintaining consistency in mapping.

[0119] This latter embodiment may mitigate the generation of numerous reflected pressure fronts, particularly when the maximal arrival time falls below a predefined threshold. This enhancement may enable the method to better address numerical artifacts associated with grid cells, such as intersecting polygons or other relevant shapes. Additionally, if a reflected segment generates secondary isochrones intersecting another boundary, a similar recursive process as described earlier can be applied, continuing until the last isochrone observed by numerical solver from the primary injection point.

[0120] The concept of the reflected segment addresses a specific scenario of a pressure front in complete reflection (no-flow), necessitating consideration of significant changes in diffusivity to refine pressure estimation. Further, a transmission coefficient can be employed to calculate a reflected pressure front and a transmitted pressure front. These distinct processes generate different isochrones, which can be remapped in a similar manner as described above.

[0121] In an embodiment, the numerical simulation's propagating pressure front may be configured to continue progressing outward upon encountering the transmissive boundary. The continuation of this pressure front can be adapted based on the desired injection period, determining the main reflection point from which the pressure front continues. Specifically, at each instance when an isochrone intersects a boundary, the cells associated with the boundary from this juncture to the last isochrone can be set to generate a "transmitted segment", i.e., a continuation of the incident pressure front that is adapted to the altered fluid behavior following the transmissive boundary.

[0122] In further embodiment, a reflected propagating pressure front can be reinitiated from this "transmitted segment" outward within the geological model, ensuring continuity from the primary propagating pressure front. The time of first arrival from the source is enforced as an initial condition for this new front propagation.

[0123] In a complementary embodiment, the numerical simulation's propagating pressure front may be configured to trigger backward reflection, thereby initiating a reflected pressure front in the opposite direction, upon encountering the transmissive boundary. The initiation of this reflected pressure front can be implemented similarly to the above embodiment of the backward reflection initiated upon encountering the zero-flow boundary. Consequently, the method provides the possibility for a transmissive boundary to include both a transmitted pressure front and a reflected pressure front as described herein. Such an implementation can improve the modelling accuracy at the cost of increased complexity.

[0124] In some embodiments the dynamic pressure disturbance simulation may provide as a solution the arrival times of the propagating pressure front at determined locations in the geological model. In line with the above embodiments, the amplitude of the propagating pressure front depends upon the arrival times and can be determined by relating the amount of injected fluid to the assigned incident and reflected arrival times.

[0125] Once the travel times are determined, they can be used to estimate the pressure amplitudes experienced at different locations within the aquifer at specific points in time. This estimation involves analyzing how the pressure changes over time as the pressure front propagates through the aquifer, influencing factors such as fluid storage, hydraulic conductivity, and porosity.

[0126] By mapping the pressure solution to travel times, it becomes possible to quantify the pressure amplitudes at various locations within the aquifer structure, providing valuable insights into fluid behavior, aquifer dynamics, and the potential for fluid storage and transport within the subsurface environment. The mapping can be implemented by applying a mathematical model that analyses the flow of fluid within aquifer structures.

[0127] In some embodiments the mathematical model may include the Thiem model and / or the Theis model as known in the art to describes the hydraulic conductivity of the aquifer and various transient conditions. Both the Thiem and Theis models are known tools in groundwater hydrology that provide insights into the behavior of aquifer structures. It should be appreciated that the disclosed method can be implemented through a different mathematical model if relevant for the accuracy of the fluid storage evaluation.

[0128] In some embodiments the isochrones of are mapped out onto an equivalent radial system, applying a simple transformation which ties the radial isochrones to the real isochrones. In such embodiments the variation in diffusivity D as described can be regarded as smooth enough to consider the isochrones of the drained volume to be mapped out directly, but the use of correction may be considered. The implementation of a radial system allows the solutions of a constant rate drawdown to be applied in transient infinite (Theis model), although it may result in a more complex implementation. Alternatively, the solutions of a constant rate drawdown to be applied as a system with imposes pressure condition at its edge (Thiem model). This is of particular interest because the drained volume obtained by the Eikonal equation reflects a zero increase of pressure condition at the edge. Consequently, variable rate can be used by convolution (superposition theorem) and is directly a summation of the next steps with a shift in time. Further at step 111, the method can include a step of generating output data based on a user selection. The output data may be used for the evaluating a storage of fluid in an aquifer structure using known analytical techniques and methods.

[0129] In some embodiments, the output data may comprise a map of isochrones generated by providing a map of arrival times by associating selected spatial points in the geological model with a corresponding arrival time of the propagating pressure front; assigning contour lines representing a plurality of isochrones based on a predetermined interval of the arrival times; plotting the isochrones on the arrival times to visually depict regions with substantially similar travel times from the predetermined injection point; and, providing the plotted isochrones as output data.

[0130] In some embodiments, the output data may comprise a map of isobars generated by providing a map of arrival times by associating selected spatial points in the geological model with a corresponding amplitude of the propagating pressure front; assigning contour lines representing a plurality of isobars based on a predetermined interval of the amplitudes; plotting the isobars on the amplitude map to visually depict regions with substantially similar amplitude from the predetermined injection point; and, providing the plotted isobars as output data.

[0131] In some embodiments, the output data may comprise at least one selected from an arrival times map depicting the arrival times of the propagating pressure front at determined locations within the aquifer structure; a pressure map depicting pressure variations at determined locations within the aquifer structure; and / or a histogram depicting the utilization of resources within the geological model during the numerical simulation, said resources preferably comprising at least one of fluid injection rate, fluid injection time, changes in pressure, changes in saturation level, storage capacity, trapping mechanics Optionally, at step 110, the preceding steps of the method 100 may be reiterated for any number of times. As previously described, during the implementation of the method 100, the simulation parameters may be determined based on a balance between computational complexity and accuracy of results. Consequently, this implies that the output data may comprise inaccuracies or estimations.

[0132] As such, the output data can be analysed to adapt one or more simulation parameters to obtain improved results, depending on the desired accuracy of the method. Subsequently, a comparison of the further output data(s) to the earlier output data(s) can be performed to evaluate the results.

[0133] At step 111, the method 100 may include the step of analysing the output data to determine the fluid storage of the aquifer structure. In specific embodiments, the fluid storage can be determined from the arrival times to determine an amount of fluid transported through the geological model that is correlated to the fluid storage. In some embodiment the evaluation can be performed with a fluid, which can be converted at a subsequent step to another fluid storage amount or mass. For example, water can be selected as the fluid that can converted into CO2 mass stored. The latter can be justified by the observation that compressed volume in the injection area is generally lower than the flow outside of the injection area. Accordingly, when the injection area covers the overall aquifer, the CO2 storage resource is directly known as the allowable compressible volume of water.

[0134] It is important to note that beyond step 111, the determined fluid storage can be utilized in subsequent stages to determine an amount of fluid that can be injected into the evaluated aquifer. In an exemplary application, the present technique can be used to inject an amount of fluid into the into aquifer structure based on the determined fluid storage by performing the method as described herein.

[0135] In certain embodiment, the system is configured for performing the method as described herein in accordance with any embodiments thereof. In a particular embodiment, the system may comprise a processing unit, which when receiving the computer program product as described herein, performs the method as described herein. It is understood that the system as described herein also includes a use of the system as described herein for inspecting an item with a feature by performing the method as described herein in accordance with any embodiments thereof.

[0136] Another aspect of the present invention relates to a processing system for assessing fluid storage in an aquifer structure, comprising at least one memory and at least one processor in communication with the memory, wherein the memory is configured for receiving sets of data as described herein as input, said data included a geological model of the aquifer; and, wherein the processor is configured for performing the (computer-implemented) method as described herein. In particular, the memory can be configured to store computer-executable instructions, and the at least one processor can be configured to access the memory and execute the computer-executable instructions to perform the (computer-implemented) method as described herein. Consequently, it should be appreciated that any embodiments of the method as described herein from embodiments of the system.

[0137] The present system will be discussed with reference to Figure 10, which is a block diagram illustrating a processing system for assessing fluid storage in an aquifer structure 500 as described herein. It is understood that additional components can be included. For example, although not shown, the method 500 may comprise a housing, various electronic, mechanical, safety, and other components that ensure proper and safe functioning of the system.

[0138] The processing system 500 may comprise at least one input means designed to receive input from a user and an output means 520 configured to provide output to the user. The input means may include various input devices such as keyboards, mice, touchscreens, voice recognition systems, or any other suitable input mechanism capable of receiving user input. Similarly, the output means may encompass display screens, speakers, printers, or any other appropriate output devices capable of conveying information to the user. These input and output means facilitate interaction between the user and the processing system, allowing for the input of commands, data, or instructions and the output of results, feedback, or visualizations.

[0139] The processing system 500 may further comprise at least one processing unit 550 configured to receive input from the at least one input means 520 and provide output to the output means 520. The processing unit 550 may further comprise a processor 570 in communication with the memory 560, configured to receive and store data on the memory 560, in order to perform the steps of the (computer-implemented) method as described herein. These processing steps may include generating or obtaining the models and simulations as described herein.

[0140] Another aspect of the present invention relates to computer program comprising instructions which, when the program is executed by a computer, performs the (computer-implemented) method as described herein. The computer program may be stored on a non-transitory computer-readable memory medium, which when loaded by at least one processor cause the at least one processor to perform the (computer-implemented) method as described herein. Consequently, it should be appreciated that any embodiments of the method as described herein from embodiments of the computer program and / or the non-transitory computer-readable memory medium.

[0141] Reference throughout this specification to "one embodiment" 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" or "in an embodiment" in various places throughout this specification are not necessarily all referring to the same embodiment.

[0142] As used herein, words of inclusion, such as "comprising", "comprises" and "comprised of" are synonymous with "including", "includes" or "containing", "contains". Unless specified otherwise, these words are open-ended without limitation and do not exclude additional, non-recited members, elements or method steps. The terms "comprising", "comprises" and "comprised of" when referring to recited members, elements or method steps also include embodiments which "consist of" said recited members, elements or method steps. Similarly, the singular forms "a", "an", and "the" include both singular and plural referents unless the context clearly dictates otherwise.

[0143] As used herein, words of approximation, such as "about," "almost," "substantially," "approximately," and the like, can be used herein to mean "at," "near," "nearly at," "within acceptable manufacturing tolerances of," or any logical combination thereof. They are intended to provide flexibility to a numerical value or range endpoint by providing that a given value may be "a little above" or "a little below" said value or endpoint, depending on the specific context. For example, the recitation of "about 30" should be construed as not only providing support for values a little above and a little below 30, but also for the actual numerical value of 30 as well. Similarly, terms "vertical" or "horizontal" are intended to additionally include "near," "nearly at," "within acceptable manufacturing tolerances of," a vertical or horizontal orientation, respectively.

[0144] As used herein, words of direction, such as "top," "bottom," "left," "right," "above," "below," "over," "under," and so on, are used for descriptive purposes and not necessarily for describing permanent relative positions. It is to be understood that such terms are interchangeable under appropriate circumstances and are intended to relate to the equivalent direction as depicted in a reference illustration; as understood contextually from the object(s) or element(s) being referenced, such as from a commonly used position for the object(s) or element(s); or as otherwise described herein.

[0145] Objects described herein as being "adjacent" to each other reflect a functional relationship between the described objects, that is, the term indicates the described objects must be adjacent in a way to perform a designated function which may be a direct ( / .e. physical) or indirect ( / .e. close to or near) contact, as appropriate for the context in which the phrase is used.

[0146] Objects described herein as being "connected" or "coupled" reflect a functional relationship between the described objects, that is, the terms indicate the described objects must be connected in a way to perform a designated function which may be a direct or indirect connection in an electrical or nonelectrical ( / .e. physical) manner, as appropriate for the context in which the term is used.

[0147] The recitation of numerical ranges by endpoints includes all numbers and fractions subsumed within the respective ranges, as well as the recited endpoints. Furthermore, the terms first, second, third and the like in the description and in the claims, are used for distinguishing between similar elements and not necessarily for describing a sequential or chronological order, unless specified. It is to be understood that the terms so used are interchangeable under appropriate circumstances and that the embodiments of the disclosure described herein are capable of operation in other sequences than described or illustrated herein.

[0148] Reference in this specification may be made to devices, structures, systems, or methods that provide "improved" performance (e.g. increased or decreased results, depending on the context). It is to be understood that unless otherwise stated, such "improvement" is a measure of a benefit obtained based on a comparison to devices, structures, systems or methods in the prior art. Furthermore, it is to be understood that the degree of improved performance may vary between disclosed embodiments and that no equality or consistency in the amount, degree, or realization of improved performance is to be assumed as universally applicable.

[0149] EXAMPLE

[0150] To better illustrate the implementation of the simulation of a dynamic pressure disturbance as described herein, an exemplary workflow will be outlined below with reference to a number of drawings. The below example is meant to aid the reader in understanding the technological concepts more easily, but it is not meant to limit the scope of the present disclosure. It is understood that additional steps can be provided before, during, and after the steps of the exemplary workflow, and that some of the steps described can be replaced or eliminated for other implementations of the dynamic pressure disturbance simulation.

[0151] The exemplary method commences with the assumption that the geological model representing the aquifer has already been provided. An example of such a model is shown in Figure 3. As a user selects an injection zone to be pressurized by DPmax, the method may initiate by finding the source using a selected computational solver from the edge of the injection zone toward its centre. In the present example the Fast-Marching Method (FMM) is selected to solve the Eikonal equation defining the cell properties.

[0152] The cell with the highest arrival time computed by the FMM is then designated as the source point and the time needed to reach the source defines the source To. While multiple suitable sources can exist depending on the shape of the injection zone, for example, by defining the same To or alternatively defining a set of local maximums {TO,}, the present example is simplified by selecting only one, although extension by superposition is feasible. An example of source selection is shown in Figure 8.

[0153] Once the source is selected, the FMM propagates across the entire grid of the geological model from the source, assigning an arrival time to each cell until the end of the injection is reached (T=Tend). The injection zone is then updated to be the isochrone To from the source point.

[0154] Based on the selected equations, such as the Eikonal solution, the radial isochrone T is mapped onto the real model T, and the pressure increase t P(t To,Tend) for To < t < Tend is calculated using constraints at isochrones To and Tend. This calculation involves a logarithmic function and a maximum pressure increase

[0155] (DPmax) defined

[0156] Alternatively, to implement the Theis transient model, the natural logarithm in the calculations can be z x z x °o e~y dy replaced by the Theis W(u) function defined as W(u) = Ju— - — . Additionally, the isochrones T can be extended beyond Tend to capture the compressed volume more accurately.

[0157] As the FMM propagates across the entire grid of the geological model, it may encounter a cell designated as a zero-flow boundary and / or a transmissive boundary as defined herein. Depending on the characteristics of the cell, three scenarios may possibly occur within the present example. These scenarios will be described separately below. However, it should be noted that in an implementation of the method, all three outcomes may occur, sequentially or simultaneously, within a single simulation as different cells with different characteristics are encountered while the FMM propagates through the entire geological model.

[0158] No boundary

[0159] In a first scenario, if none of the isochrones T reach a cell designated as a zero-flow boundary and / or a transmissive boundary as defined herein, the fluid storage resources can be calculated based on cells with arrival times higher than To and those lower than or equal to To. Specifically, if none of the cells has a AP(t; T0, Tend) higher than the fracture pressure, alternatively, if the fracture pressure has not been set outside the injection zone, the calculation can be performed as follows:

[0160] Firstly, a summation over cells with arrival time higher than To their compressed volume Vcomp=Vceu * Ct* AP(t; To, Tend. Next, a summation over cells with arrival time lower than or equal to To to get the Vcomp=Vceu * Ct* DPmax. Based on a combination of these summations, the fluid storage resources are converted in mass with fluid density within the injection area using the final pressure and temperature conditions. In a following step, the cells AP(t; T0,Tend) higher than the fracture pressure are identified, and based on this subset, the critical cell with the small arrival time Tcrit: all AP(t; T0,Tend) are rescaled by a factor equals to ^Pfrac Pini^crit-rhef |ujd resources can then be calculated by repeating the preceding steps.

[0161] Figure 11 illustrates an example of a propagating pressure front effectuated by a simulated injection of fluid from a fluid injection zone 30 into a simplified geological model without a zero-flow boundary or a transmissive boundary as defined herein. The propagating pressure front is visualized through the use of isochrones that indicate the pressure front at a corresponding time interval. In the present example, a total of four isochrones is exemplified. Specifically, the first isochrone 61 at corresponds to the pressure front at a first arrival time 1, the second isochrone 62 corresponds to a second arrival time 2, the third isochrone 63 corresponds to a third arrival time 3, and so on, until the most outward isochrone of interest 60. As can be seen, the propagating pressure front 60 propagates outward across every cell of the grid, forming an ever expanding area ranging along the trajectory arrow from the fluid injection zone 30 (To) until the final isochrone 60 (Tend).

[0162] Figure 12 shows a representation of Figure 11 as an equivalent radial system, whereby the solutions of a constant rate drawdown are known. By representing the model as a radial analytic, the derivation can be simplified as the pressure profile is entirely known, and scaling between the edge of the injection zone 30 (To) and the isochrone of interest (Tend) is sufficient to obtain the pressure field outside the injection zone for any time t > 0.

[0163] Figure 13 shows of the calculation of pressure amplitudes AP (bars) in from Figure 12 relative to the distance travelled (in meters) along the normal to isochrones (years). The results are presented in solid lines for three different injection duration times: 8 years (8y), 16 years (16), and 24 years (24y) from the start of fluid injection. Further, the correspondence between the investigated Radius (R) and the isochrones (T) is presented as a dotted line (R-T). The graph demonstrates that as the duration time increases, there is a corresponding increase in both the pressure amplitudes and the distance travelled. This trend demonstrates that longer injection periods lead to higher pressure amplitudes, resulting in fluid propagation over greater distances from the injection point. It should be noted that the curves have been cut-off to show an 'average' pressure within the injection zone equal to the pressure at its edge.

[0164] Zero-flow boundary

[0165] In a second scenario, if a subset of isochrones T < Tend reach a cell designated as a zero-flow boundary as defined herein, the computation of fluid storage resources further entails considering cells containing reflected pressure fronts generated in response to the FMM encountering the zero-flow boundary. Specifically, these reflected pressure fronts, representing cells or a group of cells with an isochron reaching a zero-flow boundary, are discerned within this subset based on consecutive decreasing arrival times.

[0166] If two reflected segments RSI and RS2 share one or more cells ceRSlnRS2, they are consolidated as a single segment only if the difference in arrival time of the cells falls below a predetermined threshold, considering that shared cells belong to the isochrone

[0167] T1i-ntersect -

[0168] For each reflected segment RSita FMM computation is conducted from it by initializing the cells of the RSLsegment with the arrival time from previous propagation, yielding a set of isochrones T^iective. Employing the image theorem and the Thiem analytical model, the isochrones T^iectivecan be projected onto the radial isochrone T, resulting in a corrective pressure increase being calculated in a similar manner.

[0169] For cells common to the initial front propagation (i.e., belonging to isochrone T) and a set of reflective ones (i.e., belonging to isochrone Tre^lective) emitted from one or several reflected segments RSL, the sum of can be employed in computing the compressed volume. If a subset of isochrones T^iectiveextends to cells belonging to a zero-flow boundary, the process described above is recursively applied until the last isochrone corresponding to Tend does not intersect any zero-flow boundary. Figure 14 illustrates an example of backward reflection initiated upon encountering a zero-flow boundary 50. In a setup similar to that of the previously discussed Figure 11, the propagation of the pressure front is initiated from a selected fluid injection zone 30 forming a source in a simplified geological model 20. As in the previous example, the pressure front continues propagating outward across every cell of the grid along the arrow trajectory resulting in a number of isochrones 60. However, in the present example, the incident pressure front encounters a zero-flow boundary 50 - indicated by a full line. This encounter results in two responses. Firstly, the initial pressure front discontinues propagating outward - indicated by the absence of the dashed line, similar to the other outward regions of the geological model 20. Secondly, a reflected pressure front is initiated from a reflection source at the zero-flow boundary 55 in the opposite direction - indicated by the dashed line directed inwards, back towards the source located at the center of the geological model 20.

[0170] Figure 15 further represents Figure 14 as an equivalent radial system. By implementing the selected analytical model, isochrones T corresponding to the reflected pressure front 55 initiated from the reflective source of the zero-flow boundary 50 can be projected onto the radial isochrone 60 of the incident pressure front from the source 30, thereby resulting in a correction of the computation of the pressure amplitude.

[0171] Figure 16 shows of the calculation of pressure amplitudes AP (bars) in from Figure 14 relative to the distance travelled (in meters) along the normal to isochrones (years). The results are presented as solid lines for three different injection duration times: 8 years (8y), 16 years (16), and 24 years (24y) from the start of fluid injection. Further, the correspondence between the investigated Radius (R) and the isochrones (T) is presented as a dotted line (R-T). The graph depicts that the 16-year duration allows reaching the previously mentioned zero-flow boundary 50; however, the associated pressure is insufficient to initiate a reflected pressure front. Conversely, the graph demonstrates that the longest duration of 24y leads to reaching the zero-flow boundary 50 and, upon contact, initiates a reflected pressure front labeled as '55'. Notably, this reflected pressure front 55 exhibits a decrease in amplitude in a manner akin to the theoretical continuation of the incident pressure front, denoted by the dashed line labeled '56', had the zero-flow boundary 50 not stopped its progress. Consequently, the pressure profile for 24y will result in an increased pressure value corresponding to the summation of the pressure values from the pressure front 55 and the discontinued incident pressure front 56 as indicated by the double-sided arrow. This observation underscores that longer injection periods result in higher reflected pressure fronts.

[0172] Transmissive boundary In a third scenario, if a subset of isochrones T < Tend reach a cell designated as a transmissive boundary as defined herein, the computation of fluid storage resources further entails considering cells containing transmissive generated in response to the FMM encountering the zero-flow boundary. Specifically, these transmitted pressure fronts, representing cells or a group of cells with an isochron reaching a transmissive boundary, are discerned within this subset based on consecutive decreasing arrival times.

[0173] The interface HD corresponding to the transmissive boundary is characterized by cells with a significant contrast in diffusivity Dout / Dinand / or a significant contrast of mobility Mout / Min. Here, Doutdenotes the diffusivity of the downstream isochrone cell (i.e., further away from the source), while Dinsignifies the diffusivity of the upstream isochrone cell (i.e., closer to the source).

[0174] FIG. 20 is an example of a contrasting mobility M ratio between an inward cell (Min) and an outward cell (Mout) along a transmissive boundary.

[0175] To accommodate this condition, the procedure previously outlined in for a zero-flow boundary is adapted to employ a reflection coefficient RD .n Dout(T and a transmission coefficient TD .n Dout(T = 1 + RDinDout(T) along the interface. Consequently, the contribution of each coefficient towards the fluid resources computation is described below.

[0176] For a normal incident pressure front encountering an infinite plane separating porous media in and out, with D and M representing their respective diffusivity and mobility, the reflection coefficient is calculated

[0177] Under the mapping out process to a radial model, invoking the relationship between drained volume and

[0178] K,(t) diffusivity D = — — the reflection coefficient is «

[0179] Two FMM processes, denoted as FMMRand FMMT, are conducted along cells belonging to the interface HD and adjacent cells, corresponding to isochrones Thit<T < Tend. These processes enable the definition of isochrones from the HD interface. At any time from Thit, especially at Tend, isochrones TrHeDflectlveand TtHrDansmitive, respectively, define •

[0180] On each side of interface HD, an average mobility M can be computed as the geometrical mean using cells belonging to the interface and isochrones Thit<T < Tend.

[0181] Similar to the calculations of the reflected segment for the second scenario, the compressed volume per cell comprises P(t; To, Tend) along with sums of reflection and transmission coefficients resulting in

[0182] FIG. 21 is an example of a determination of a Volume V corresponding to an isochrone T of a transmitted segment (VpT) and a reflected segment (VpR) of the pressure front.

[0183] Figure 17 illustrates an example of transmission initiated upon encountering a transmissive boundary 40. In a setup like that of the previously discussed Figure 11, the propagation of the pressure front is initiated from a selected fluid injection zone 30 forming a source in a simplified geological model 20. As in the previous example, the pressure front continues propagating outward across every cell of the grid along the trajectory resulting in a number of isochrones 60. However, in the present example, the incident pressure front at a specific time encounters a transmissive boundary 40 - indicated by a full line. This encounter results in two responses. Firstly, the incident pressure front continues propagating outward - indicated by the absence of the dotted line, similar to the other outward regions of the geological model 20, the pressure front is altered based on the cell properties after the transmissive boundary 40. Secondly, a reflected pressure front is initiated from a reflection source at the transmissive boundary - indicated by the dashed line directed inwards, i.e., in the opposite direction towards the source located at the center of the geological model 20.

[0184] Figure 18 further represents Figure 17 as an equivalent radial system. By implementing the selected analytical model, isochrones corresponding to the reflected pressure front 55 initiated from the reflective source of the transmissive boundary 40 can be projected onto the radial isochrone 60 of the incident pressure front from the source 30, thereby resulting in a correction of the computation of the pressure amplitude. Additionally, the isochrones 45 of the transmitted pressure front of the incident pressure front 60 can be altered based on the cell properties after the transmissive boundary 40.

[0185] Figure 19 shows of the calculation of pressure amplitudes AP (bars) in from Figure 17 relative to the distance travelled (in meters) along the normal to isochrones (years). The results are presented as solid lines for three different injection duration times: 8 years (8y), 16 years (16), and 30 years (30y) from the start of fluid injection. Further, the correspondence between the investigated Radius (R) and the isochrones (T) is presented as a dotted line (R-T). The graph demonstrates that the 16-year duration enables reaching the previously mentioned transmissive boundary 40; however, the associated pressure is insufficient to trigger a transmitted pressure front. In contrast, the graph shows that the longest duration of 30 years allows reaching the transmissive boundary 40 and, upon contact, initiates a transmitted pressure front labeled as '45'. Consequently, the pressure profile for 24y will result in a pressure value that can be determined as a summation of these three pressure profiles as further described in Figure 22 below. Notably, the contrast in diffusivity / mobility across the transmissive boundary 40 results in a lesser decrease in amplitude compared to what would be expected for a hypothetical continuation of the pressure front in the absence of the transmissive boundary 40.

[0186] Figure 22 illustrates the calculation of pressure amplitudes AP (bars) relative to isochrones (years) along the radius of investigation (meters) from the source, represented by the injection point in the geological model. The results are presented for the scenario depicted in Figure 19, where, at 16 years (16y), the pressure front encounters a high diffusivity contrast, transitioning from a value of 0.63 m2 / s to a value of 2.80 m2 / s. The graph demonstrates that the incident propagating pressure front 35 generates a certain pressure amplitude. At the equivalent radius for the isochrone 16y, a reflective pressure front 55 is generated, akin to the scenario with a no-flow boundary, but its amplitude is influenced by a reflective coefficient. Meanwhile, the transmissive pressure front 45 continues from isochrones 16y to 30y using the 2.8 m2 / s diffusivity, with its amplitude affected by a transmissive coefficient T. These three pressure profiles are summed-up over an equivalent investigation radius (R) to generate the 30y pressure profile depicted in Figure 19.

[0187] Figures 23-26 further present an implementation of the above discussed exemplary workflow on an exemplary geological model 20.

[0188] Figure 23 displays the arrival times for an incident pressure front (IPF) initiated from the the selected source 30 across the cells of the geological model, following the arrow trajectories. As the incident pressure front propagates, arrival times are assigned to each cell until the outer perimeter of the geological model is reached. The expanding lines represent isochrones with similar arrival times. Zeroflow boundaries are identified along the outer perimeter of the geological model, labelled as (A), (B), and (C). Backward reflection will be initiated at each of these zero-flow boundaries as described below.

[0189] Figure 24 displays the arrival times for a first reflected pressure front consisting of a reflected pressure front (RPF) initiated from a first selected reflection source at zero-flow boundary labelled as (A) in Figure 22. Notably, as the (A) flow boundary is closest to the source (30), the reflected pressure front will travel the furthest distance with the highest amplitude relative to the reflected pressure fronts of (B) and (C).

[0190] Figure 25 displays the arrival times for a second reflected pressure front consisting of a reflected pressure front (RPF) initiated from a second selected reflection source at the zero-flow boundary labelled as (B) in Figure 22. Notably, as the (B) flow boundary is furthest from the source (30), the reflected reflected pressure front will travel the shortest distance with the lowest amplitude relative to the reflected pressure fronts of (A) and (C).

[0191] Figure 26 displays the arrival times for a third reflected pressure front consisting of a reflected pressure front (RPF) initiated from a third selected reflection source at the zero-flow boundary labelled as (C) in Figure 22. Figure 27 presents a pressure distribution map based on the simulation for Figure 23. This map is generated by mapping the pressure amplitudes from the Theis Model onto the Arrival time FMM Isochrones and then aggregating the respective values according to the computational procedures outlined in the preceding sections. The resulting pressure map offers a visual representation of the spatial distribution of pressure (bars) delineated through isobars within the geological model. Notably, the map illustrates that pressure amplitude values decrease as distance from the source increases.

[0192] Figure 28 shows a histogram of fluid storage resources in metric tons (MT) determined from the simulation for Figure 23 with CO2 as the selected fluid. The values are generated through the computation described in the embodiments above. Each bar in the histogram represents a range of fluid storage capacities, with the height of the bar indicating the frequency or count of occurrences within that range.

[0193] By combining these values, the cumulative volume of CO2 that can be stored within the aquifer structure can be determined.

Claims

CLAIMS1. A computer-implemented method (100) for evaluating a storage of fluid in an aquifer structure, the method comprising the steps of: acquiring a geological model representing the aquifer structure; the geological model defining a plurality of cells, each cell being provided with properties that define the fluid behavior as a function of spatial position within the geological model; selecting at least one cell forming a source for fluid injection into the geological model; implementing a numerical simulation describing an injection of predetermined amount of fluid into the geological model; wherein the numerical simulation is configured to initiate an outward propagation of an incident pressure front from the source across the cells of the geological model until a stop condition is met, thereby assigning a set of incident arrival times to each cell as the incident pressure front propagates outward; wherein the numerical simulation is further configured to identify at least one cell forming a zeroflow boundary where fluid behavior is substantially prevented, prompting the numerical simulation to discontinue the propagation of the incident pressure front, and / or at least one cell forming a transmissive boundary where fluid behavior is substantially altered based on a predetermined threshold, prompting the numerical simulation to adapt the outward propagation of the incident pressure front to the altered fluid behavior, wherein the numerical simulation is further configured to select at least one cell forming a reflection source at the at the zero-flow boundary and / or the transmissive boundary, and initiate an inward propagation of a reflected pressure front from the reflection source across the cells of the geological model until a stop condition is met, thereby assigning a set of reflected arrival times to each cell as the reflected pressure front propagates inward; determining a set of pressure amplitudes for each group of cells with substantially similar arrival times by relating the amount of injected fluid to the assigned arrival times; and, determining a fluid storage of the aquifer structure from the sets of pressure amplitudes and arrival times.

2. The method according to any of the preceding claims, wherein the cell properties include at least a Mobility M quantifying the flow of fluid through the cell under pressure, a Storativity S quantifying the ability of the cell to store fluid under changes to the pressure; and a Diffusivity D quantifying the rate at which the injected fluid propagates the pressure front through the cell, whereby the Diffusivity D is defined as a ratio between the Mobility and the Storativity, D=M / S.

3. The method according to claim 2, wherein the numerical simulation is configured to identify the transmissive boundary based on a ratio of the Mobility and / or Diffusivity between an inward cell and an outward cell along a trajectory of the incident pressure front, and further adapt the outward propagation of the incident pressure front based on said ratio.

4. The method according to any of the preceding claims, wherein the numerical simulation is configured to group a plurality of cells with substantially similar incident arrival times to form the transmissive boundary and / or the zero-flow boundary, and select at least one cell from this group of cells to form the reflection source.

5. The method according to any one of the preceding claims, wherein the pressure amplitude determination comprises grouping cells of substantially similar arrival times, and relating the pressure amplitude to the volume of injected fluid through the porosity and compressibility, and relating the mass of injected fluid to the volume through the temperature and density of injected fluid.

6. The method according to any of the preceding claims, wherein the method comprises receiving a set of simulation parameters as input, the simulation parameters including at least a fluid injection zone selected from at least one cell or a group of cells capable of forming the source for fluid injection into the geological model, a fluid injection time defined as the duration over which fluid is injected into the geological model, and a fluid injection target representing a pressure range to be achieved at the conclusion of the fluid injection into the geological model.

7. The method according to claim 6, wherein the fluid injection zone comprises a plurality of cells; wherein the method further comprises determining the source by applying a numerical simulation configured to initiate an inward propagation of an pressure front from a plurality of positions selected along the perimeter of the fluid injection zone, thereby assigning a set of incident arrival times to each cell as the pressure front propagates inward, and selecting the source from one or more cells with substantially similar arrival times from each pressure front.

8. The method according to any of the preceding claims, wherein solving the numerical simulation comprises selecting a numerical solver configured to solve the governing equations for each cell as the incident and / or reflected pressure fronts progress across the cells of the geological model; wherein the selected numerical solver includes the Fast-Marching Method (FMM) configured to solve the Eikonal equation included as governing equation.

9. The method according to any one of the preceding claims, wherein determining the set of pressure amplitudes comprises mapping the pressure amplitudes to the arrival times through ananalytical model describing the behavior of fluid within an aquifer; wherein the selected analytical model includes the Thiem and / or the Theis model.

10. The method according to any one of the preceding claims, wherein the fluid comprises water and / or air, preferably CO2.

11. The method according to any one of the preceding claims, further comprising the steps of: generating a map of arrival times by associating selected spatial points in the geological model with a corresponding arrival time and / or amplitude of the propagating pressure fronts; assigning contour lines representing a plurality of isochrones and / or isobars based on a predetermined interval of the arrival times and / or amplitudes; plotting the isochrones and / or isobars on the arrival times and / or amplitude map to visually depict regions with substantially similar travel times and / or amplitude from the predetermined injection point; and, providing the plotted isochrones and / or isobars as output data.

12. The method according to any of the preceding claims, further comprising the step of generating output data, the output data comprising at least one selected from an arrival times map depicting the arrival times of the propagating pressure fronts at determined locations within the aquifer structure; a pressure map depicting pressure variations at determined locations within the aquifer structure; and / or a histogram depicting the utilization of resources within the geological model during the numerical simulation, said resources preferably comprising at least one of fluid injection rate, fluid injection time, changes in pressure, changes in saturation level, storage capacity, trapping mechanics.

13. A method of injecting fluid into an aquifer structure, comprising the steps of determining the fluid storage of the aquifer structure by performing the method according to any one of claims 1 to 12; injecting an amount of fluid into the into aquifer structure based on the determined fluid storage.

14. A system (500) for assessing fluid storage in an aquifer structure, comprising a memory (560) and a processor (570) in communication with the memory, wherein the memory is configured for acquiring a geological model representing the aquifer structure; and, wherein the processor is configured for performing the method according to any one of claims 1 to 12.

15. A computer program comprising instructions which, when the program is executed by a computer, performs the method according to any one of claims 1-12.

Images