The 11th Society of Petroleum Engineers Comparative Solution Project: Problem Definition

This document contains the description of, and call for participation in, the 11th Society of Petroleum Engineers Comparative Solution Project (the 11 th SPE CSP, https://spe.org/csp). It is motivated by the simulation challenges associated with CO 2 storage operations in geological settings of realistic complexity. The 11 th SPE CSP contains three versions: Version 11A is a 2D geometry at the laboratory scale, inspired by a recent CO 2 storage forecasting and validation study. For Version 11B, the 2D geometry and operational conditions from 11A are rescaled to field conditions characteristic of the Norwegian Continental Shelf. Finally, for Version 11C, the geometry of Version 11B is extruded to a full 3D field model. The CSP has a two-year timeline, being launched at the 2023 SPE Reservoir Simulation Conference, and culminating at the 2025 SPE Reservoir Simulation Conference. A community effort is run in parallel to develop utility scripts and input files for common simulators to lower the threshold of participation; see link to supplementary material on the CSP website. At the time of writing, complete input decks for one simulator are already ready for all three versions.


Introduction
Safe and efficient implementation of geological carbon storage (GCS) necessarily relies on reservoir simulators applied to uncertain geological data.While the strengths and limitations of reservoir simulation are well appreciated within petroleum production, GCS raises new challenges both in terms of physical processes and timescales.As an example, the enhancement of dissolution from a CO2-rich supercritical phase to the aqueous phase through convective mixing ensures important long-term storage security, relevant on timescales from decades to centuries.
One consequence of the relative youth of the GCS industry, combined with the long timescales and new physical processes of interest, is that available field data for validation of simulation technology is still rare.This increases the importance of validation against proxy systems, and code verification through comprehensive benchmarking efforts among simulators.

Background and motivation
During 2021-2022, three of the present organizers lead a forecasting and validation study within the academic GCS community (Nordbotten et al., 2022;Flemisch et al., 2023), as illustrated in Figure 1.The primary intent was to validate the long-term performance of numerical simulators for GCS, with particular emphasis on the post-injection period, and to assess the ability to state accurately wellcalibrated forecasting intervals.The study also revealed several numerical challenges, both in terms of numerical accuracy when resolving the reservoir dynamics, and in terms of obtaining good solver performance (see e.g., Flemisch et al (2023), Salo-Salgado et al (2023), Green et al (2023) and Wapperom et al (2023)).
Separately, the development of numerical simulation capabilities for subsurface applications has historically benefited substantially from common reference simulation cases, notably the series of 10 Comparative Solution Projects (CSP) organized within the Society of Petroleum Engineers (SPE) between 1981 and 2001 (see Islam and Sepehrnoori (2013) for a review).These observations provided the initial motivation for developing a set of benchmark cases for CO2 storage within the concept of a new SPE Comparative Solution Project (CSP).
In developing this 11 th SPE CSP (https://spe.org/csp),we hope to provide a common platform and reference case for numerical simulation of GCS.Specifically, we anticipate that the following topics are discussed relative to this baseline.
• Development and verification of accurate and efficient discretization methods for multiphase, multicomponent flow and transport.• Development and verification of space-time adaptive gridding and domain decomposition methods.• Development of upscaling methods for convective mixing and dispersion of in the context of CO2 dissolution into water.• Development and verification of robust and efficient linear and non-linear solvers and solution and time-stepping strategies for 2D and 3D at laboratory and field conditions.• Assessment of the importance of physical processes omitted from this study, including (but not limited to) geochemical reactions, mechanical response, and more realistic boundary conditions.
Furthermore, as of the date of launching this CSP, we do not anticipate that a fully converged solution (in the sense of grid refinement) will be achievable for any of the three versions of the CSP by means of standard numerical methods on desktop hardware.This anticipation is justified in part based on the experiences from the previous academic benchmark cited above, and in part due to the known challenges associated with accurately simulating the dissolution and convective mixing of CO2 into water, in particular outside the context of trivial geometries, see e.g., Pau et al. (2010).As such, providing a reference solution to this CSP is a challenge in itself, that will likely require advanced methods and high-performance computing.

Comparative Solution Project overview
This 11 th SPE CSP (https://spe.org/csp) is based around a synthetic geological cross-section, representative of the structures found in the Norwegian Continental Shelf.The cross-section is defined in terms of seven facies and is used to realize three versions of the CSP: Version 11A is defined as a 2D experiment at laboratory scale and surface conditions.Version 11B is defined as a 2D transect at field scale and conditions.Version 11C is a 3D extrusion of version 11B, thus corresponding to a full synthetic field study.The three versions 11A, 11B and 11C all contain different challenges, and should not be considered to be of increasing difficulty.Indeed, many may find Version 11B to be the easiest, and we recommend starting with this case.
A common 2D transect geometry is utilized for all three versions of the CSP and is designed to strike a balance between a) simplicity of definition, b) geological and operational realism, and c) providing several key computational challenges associated with numerical simulation of CO2 injection, migration, and long-term storage in water-filled porous media.
We recognize that to provide a successful CSP, a balance needs to be attained between including a maximum amount of geological, operational, and physical realism versus a CSP that can be well addressed by a broad spectrum of the community.Such considerations permeate the whole development of this CSP, and some of our main considerations are highlighted as follows.
• Geologic realism: Our goal has been to provide a geometry that is easy to understand and implement, while at the same time conceptually corresponds to a realistic geological formation in the Norwegian Continental Shelf.We have chosen to address this by using a common 2D transect in all three versions of the CSP, that itself was designed to balance complexity with computational addressability.This 2D transect is then stretched, extruded, and deformed between the three variants of the CSP.• Heterogeneity: We have chosen to emphasize structural heterogeneity (facies and faults) over local-scale heterogeneity.This choice is justified by two considerations.Firstly, local-scale heterogeneity is covered by the 10 th SPE CSP.Secondly, we wish to have a CSP that in principle has a computable solution and can be reasonably well approximated by the majority of current simulators, that would be prohibitive with fine-scale (or stochastic) heterogeneities.• Fluid complexity and geochemistry: We have chosen to specify two-phase, two-component flow with thermal effects (Versions 11B and 11C only).This corresponds to pure water in the reservoir, with pure CO2 injection.The omission of salts and other minerals in the system is in part due to the relatively poor characterization of the properties of brine-CO2 mixtures in the literature, and to avoid modeling permeability as a function of salt precipitation and other geochemical processes.These choices substantially simplify the problem definition.The inclusion of thermal effects is an acknowledgment that CO2 may frequently be injected at a temperature significantly below the reservoir temperature.• Geomechanics: Geomechanics is not included in the CSP based on two arguments.Firstly, the fluid-mechanical coupling is primarily a total pressure coupling, and we thus consider the multi-phase flow processes to be a somewhat orthogonal challenge to poromechanics.As such, we believe they can be to a large extent studied separately.Secondly, we acknowledge that a substantial fraction of available reservoir simulators do not include a native coupling to mechanical response, and wish that the CSP should also be open also to these participants.• Petrophysical processes: After careful consideration, we have chosen to include capillary forces and dispersion in the description, but not relative permeability hysteresis.We justify this choice as we consider capillary forces and dispersion to have greater conceptual importance as they represent important terms in the governing equation during injection and post-injection, respectively, and as they are necessary for a grid-converged solution to exist.
On the other hand, hysteresis is of lesser conceptual importance as it is a modification of existing constitutive laws.
The choices just described will invariably exclude active research and development interests of many members of the community, and we therefore emphasize that this CSP is meant to be used as a common baseline, from which the incorporation of additional complexities, not included in the CSP, is more than welcome.
Both the injection protocol, as well as the target quantities for the simulation, are defined with operational realism in mind.In particular, simulation results are to be reported both in terms of spatial maps of the field variables (pressure, saturation, phase composition, and for versions 11B and 11C also temperature).Furthermore, we request target quantities that represent proxies for assessing the longterm security of carbon storage.These are motivated by the following key questions and variables, made precise in Sections 2.5.and 2.6: P1.As a proxy for assessing risk of mechanical disturbance of the overburden: Time evolution of pressure at two observation points.
P2.As a proxy for when leakage risk starts declining and our ability to simulate accurately phase partitioning: Time evolution of phase partitioning of CO2 within a region covering the primary storage dome.
P3.As a proxy for our ability to handle more complex geological features: Time evolution of phase partitioning of CO2 within a region covering a secondary storage dome.
P4.As a proxy for our ability to capture onset of convective mixing: Time evolution of a measure of density-driven fingers in a region below the gas-water interface of the primary storage dome.
P5.As a proxy for our ability to capture migration into low-permeable seals: Time evolution of the total mass of CO2 in the seal facies.

Community resources
The current document constitutes the official CSP description.The development of this description has been paired with an ongoing community effort to develop utilities that lower the threshold for participation.This community effort is hosted at https://github.com/Simulation-Benchmarks/11thSPE-CSP,and at the time of writing contains: • Geometry files and scripts for producing structured meshes for all versions of the CSP.
• Scripts for generating thermodynamic tables based on the NIST database.
Sample input files for common simulators are in development, and corresponding input decks for one simulator are already available for all three versions of the CSP.We emphasize that the input decks are intended to support the participants in setting up their computational workflows.Especially, the contained grids are expected to be optimized by the participants for generating their respective CSP results.

Call for participation and timeline
Participation in this 11 th SPE CSP is open to all interested parties, subject to the condition that the agreement of participation must be completed, signed, and submitted by December 1 st , 2023.
The CSP is planned according to a 2-year cycle, with tentative milestones as follows:

Early Access Team
The authors would like to thank and acknowledge the contributions of the "Early Access Team": The early access team has supported the development of the CSP through giving feedback on the CSP description, quality control of the proposed versions of the CSP, and providing supporting material in the form of example computational grids and run files.This team also formed the initial core of the community developing the resources mentioned in Section 1.3.

Technical description CSP 11A
The CSP 11A is a 2D geometry, set at laboratory scale and conditions.Version 11A is closely inspired by an actual CO2 storage lab experiment, and as such there is an actual physical solution that provides insight into the general system behavior (Fernø et al, 2023).The distribution of gas-phase CO2 and dissolved CO2 at the end of injection can be seen in Figure 2.
Attempting to reproduce actual CO2 storage experiments in silico is important to assess the predictive capabilities of simulators.The experiment underlying CSP 11A was conducted at atmospheric conditions and this introduces challenges that will not necessarily be seen in simulations conducted at reservoir conditions, including large density differences between CO2 and brine and diffusion being a determining factor for the onset and wavelength of the viscous fingers, as discussed in more detail by Wapperom et al. (2023).Version 11A will likely prove to be more difficult to simulate than Version 11B, and even though we describe most of the general setup of the CSP for Version 11A, we recommend Version 11B as a starting point for most participants.

Governing equations and constitutive laws
As governing equations, we state the standard isothermal two-phase, two-component extension of Darcy's law.For a detailed description, refer to Nordbotten & Celia (2011) or Lake et al. (2014); the main equations are summarized below.We emphasize that all parameters and constitutive functions are defined as constant within each facies (see the geometry description in the next section and parameters in Section 2.3).
Multi-phase Darcy's law for phases  =  (CO2-rich non-wetting gas phase), and  =  (H2O-rich wetting liquid phase): Here,   is the volumetric flux of phase ,   is the phase pressure,  , and  are the relative and intrinsic permeabilities, respectively,   is the phase viscosity,   is the phase density, and  is the gravitational force, defined with three significant digits, || = 9.81 ms −2 , pointing "down".(2.5) In addition to the governing equations provided above, the following constitutive laws are considered.
Brooks-Corey type relative permeability and capillary pressure: For primary drainage (initial period of injected gas displacing water) we consider normalized saturations: , 0) and  , = max ( where  ,imm is the saturation below which the phase is immobile.The relative permeability is given as a unique function of saturation by:  , (  ) = ( , ) ,1 (2.7) Here, the exponents  ,1 determine the non-linearity of the relative permeability.Similarly, the basic Brooks-Corey capillary pressure is given by: Here,   is the entry pressure for the rock.We remark that capillary pressure is only physically meaningful when the phases are connected.For non-wetting saturations less than  , , the phase will exist as disconnected bubbles of various radii, and thus a unique non-wetting pressure may not exist.This has no impact on the flow calculations (because the non-wetting phase relative permeability is zero), but this will impact thermodynamical calculations.For this CSP, we resolve this issue by defining an extended capillary pressure, valid for all saturations, as: Here,  cap,max defines a maximum capillary pressure, while the error function ensures a smooth transition of the capillary pressure function to the maximum value.For completeness, we emphasize that as the error function is defined as the integral of the normal distribution, it is well-defined at infinity: erf(∞) = 1.Capillary pressure-saturation relationships are illustrated in Section 2.3 and discussed in more detail in Section 5.1 and 5.4.
Thermodynamics: The description of the thermodynamics is split into three parts: Phase partitioning, pure-phase properties, and mixture properties.
1) We define the phase partitioning of the pure CO2-H2O system according to Spycher, Pruess and Ennis-King (2003); see in particular section 4.3 and Tables 1 and 2. To be precise, we emphasize that for the purpose of this CSP, the solubility limits of each phase (i.e.CO2 solubility limit in the water-rich wetting phase and H2O solubility limit in the CO2-rich non-wetting phase) should be calculated based on the pressure of the same phase.Mass transfer between phases is assumed to be instantaneous at each point, immediately dissolving/vaporizing available CO2 and H2O until the solubility limit is reached.
2) The pure-phase CO2 and H2O properties are defined according to the NIST database, see https://webbook.nist.gov/chemistry/fluid/,(Lemmon et al., 2023).The respective phase pressure is used when evaluating the pure phase properties.3) For the range of conditions considered herein, the mutual solubilities are quite small.
Rock compressibility: The rock is considered incompressible, thus  does not vary over time.

Geometry, boundary, and initial conditions
The geometric description is motivated by laboratory experiments relevant for North Sea storage formations and has been developed in consultation with faculty and researchers at the Department of Earth Science, University of Bergen 6 The geometry presented here is a simplification of the geometry defined in Nordbotten et al. (2022).
The length of the porous medium is 2.8 m and height is 1.2 m.Version 11A of the CSP is 2D, but in order to present quantities in common units of "mass per volume", we consider a uniform depth in the third dimension of 0.01 m.The porous domain consists of the full porous medium, except for the two injection wells, as specified below.Figure 3 provides a sketch of the geometry; a precise definition is given in the gmesh-compatible file spe11a.geo,available as part of the benchmark description7 .As evident from the figure, the geometry contains seven facies, all of which are considered to be internally homogeneous.The facies properties are given in the next section.
The figure contains three boxes that are used for reporting.Their (, ) coordinates are measured relative to the lower left corner of the domain and are specified in terms of their bottom left and top right corners: Boundary conditions: The left, right, and bottom boundaries are impermeable, i.e.: ⋅  = 0 and    ⋅  = 0 (2.12) where  is the (outward) normal vector to the boundary.
The top boundary is considered a constant pressure boundary in contact with pure water, i.e.: Here,  ,0 is the boundary pressure, defined as  ,0 = 1.1 ⋅ 10 5 Pa.
The two injection wells for pure CO2 injection,  = 1,2, are also defined as (internal) boundary conditions.Both wells are defined to have a radius of  = 9 ⋅ 10 −4 m, centered at their respective coordinates.We refer to the respective well-to-reservoir circular boundaries as Γ  , where we impose the boundary conditions: Here   () are the injection rates per unit depth (injection rate per well length for 11C), specified in the section "Operational conditions".The injection is pure CO2, thus Eq. (2.14) is complemented by: Initial conditions: The initial conditions represent a water-filled medium at rest, compatible with the boundary conditions, i.e.: = 1,   H 2 O = 1 and   = 0 (2.16)

Facies properties
The geometry contains seven facies: One seal (facies 1), five permeable reservoir sands (facies 2-6) and one impermeable (facies 7).These provide the definition of the material properties, as given in Tables 2.1 and 2.2.These properties are consistent with laboratory measured parameters for unconsolidated sands (Nordbotten et al., 2022;Fernø et al, 2023).We remark that the maximum capillary pressure is chosen to avoid any potential for phase transition from water to vapor.
As examples, capillary pressure and relative permeability functions are given in Figure 4.
We emphasize that to obtain 3D rates, one must multiply by the nominal thickness of the geometry, thus the injection rates correspond to 1.7 ⋅ 10 −7 kg s .Exemplifying these injection conditions, the specified mass rate corresponds to a volumetric rate of about 5 cubic centimeters per minute at 20 °C and 1.1 ⋅ 10 5 Pa.

Measurables
In addition to spatial maps of the field variables, the following measurables shall be reported, details of which, including time resolution, are described in Section 2.6.These measurables correspond to proxies for the motivating questions in Section 1.2.

Pressure (Proxy P1)
Pressure shall be reported at each of the two pressure observation points, POP 1 and POP 2 in [Pa].

Phase composition (Proxies P2 and P3)
The distribution of the CO2-phase shall be reported within boxes labeled A and B in Figure 2. The phase distribution (in kg) shall be reported in the following categories: 1. Mobile free phase (CO2 at saturations for which the non-wetting relative permeability exceeds 0); 2. immobile free phase (CO2 at saturations for which the non-wetting relative permeability equals 0); 3. dissolved (CO2 in water phase) and 4. seal (CO2 in any form in facies 1).The sum of the three first categories (mobile, immobile, and dissolved) shall equal the total mass of CO2 in the respective box.

Convection (Proxy P4)
For the box labeled Box C in Figure 2, the extent of convective mixing shall be reported as the integral of the magnitude of the gradient in relative concentration of dissolved CO2.In other words, given the mass-fraction of CO2 in water    , and the solubility limit is denoted  ,max  , then the following quantity  shall be reported: (2.17) This quantity corresponds to the normalized total variation of the concentration field within Box C, and (taking into account the nominal depth of the geometry) has units of meters squared.

CO2 in sealing units (Proxy P5)
Similar as in Section 2.5.2, the total mass of CO2 in all sealing units (CO2 in any form in facies 1) shall be reported in kg.

Data reporting
All result data will be uploaded by the participants to the dedicated website provided by SPE.An account at spe.org is required for the upload process which is described in detail on the website.At the end of the CSP, corresponding to the release of the final report, all submitted data will be turned public.
The reported data will be analyzed in two respects: Both in terms of an intercomparison of general numerical simulation capability (global spatial maps), but also in terms of our ability to correctly assess key properties of the system (the measurables outlined in Section 2.5).Consequently, we establish both a "dense" and a "sparse" reporting protocol.
Groups are encouraged to submit up to four results, if they deem it interesting, of which at least one should be representative of a computation that is reasonable within common reservoir engineering practice.

Sparse data
All measurables identified in Section 2.5 shall be reported at 600s (10-minute) intervals starting at the initial injection and lasting 432000s (5 days).The data is expected in csv format in a file spe11a_time_series.csv of the form according to the measurables defined in Section 2.5.

Dense data
A spatial map of all field variables (pressure, saturation, phase composition) shall be reported for each hour from injection start.While the computational grids are generated by the participants individually and should be chosen by each group as they find most appropriate, for cross-group comparison the spatial maps shall be reported relative to a uniform Cartesian grid of 280 by 120 cells (.01 m by .01m grid cells from the bottom of the domain).For each temporal snapshot indicated by H hours, H = 0, 1, The origin of the coordinate system should be located in the lower-left corner with the x-axis positively oriented towards the right and the z-axis positively oriented towards the top.(The reported x and y values refer to the lower-left corners of each cell in the uniform report grid.)Moreover, note that intensive variables (pressure, saturation, and mass fractions) should be reported as cell-center values, while extensive variables (total mass) should be reported as integral/average values for the cell.

Performance data
Reporting of performance data is strongly encouraged to the extent possible, but not mandatory.
Note that several reporting quantities may not be relevant for certain participating groups, depending on their choice of numerical method and solution strategy.Performance data should be reported in three categories: As time-series (similar to sparse data in Section 2.6.1), as spatial maps (similar to dense data in Section 2.6.2) and as a questionnaire.These are detailed below.The participants should stick to the prescribed format, particularly the number of columns in the csv files, even if not all requested values can be reported.A missing value should be indicated by n/a.
The participants are encouraged to provide time-series at 600s (10-minute) intervals starting at the initial injection and lasting 432000s (5 days).The data is expected in csv format in the repository in a file spe11a_performance_time_series.csv of the form The reporting quantities provided as time-series are defined as: • tstep: Average time-step size over the last 600s.
• fsteps: Number of failed time steps over the last 600s.
• mass: Mass balance, i.e., total mass of CO2 in the domain plus any mass that has crossed the boundaries.• dof: Average number of degrees of freedom per time-step over the last 600s.
• nliter: Total number of non-linear iterations for all time-steps (both failed and converged) used to advance the solution the last 600s.• nres: Number of local residual evaluations (e.g., per cell) for all time-steps the last 600s.
• liniter: Total number of linear iterations spent to advance the solution over the last 600s.
• runtime: Total runtime for advancing the solution the last 600s.Details of the computational platform should be reported in the questionnaire, i.e., type of hardware (CPU/GPU model and make, amount and type of memory, memory bandwidth), how this was used (number of threads/cores utilized, GPU offloading, etc.) • tlinsol: Total time spent in the linear solver if applicable.
It is encouraged to include an extra file spe11a_performance_time_series_detailed.csv that reports the same quantities for all the time steps and not just for the 600s report steps.
The participants are encouraged to provide spatial maps for each hour from injection starts.For crossgroup comparison the spatial maps shall be reported on a uniform Cartesian grid of 280 by 120 cells (.01 m by .01m grid cells from the bottom of the domain).For each temporal snapshot indicated by H hours, ( = 0, 1, 2, … ), cell values should be provided in csv format in a file spe11a_performance_spatial_map_<H>h.csv of the form: The reporting quantities provided as spatial maps are defined as: • cvol: Average geometric cell volume of the computational grid inside the given 1cm x 1cm cell of the reporting grid at latest time-step.• arat: Average aspect ratio (vertical/lateral) of the computational grid used inside the 1cm x 1cm reporting cell.and the surface density   of the respective component.If this quantity cannot easily be computed for your numerical method, please specify how you define the mass-balance error in the questionnaire.
• post_est: A posteriori estimate of the error between the computed and the true solution within the 1cm x 1cm reporting box at convergence of latest time-step.Please specify how you define the error estimate in the questionnaire.
For both sparse and dense data, the participants are encouraged to add columns with extra performance metrics that are relevant for their solution approach.As an example, if adaptive implicit methods (AIM) are used, extra performance data may be spatial maps of implicit/explicit variables.Moreover, we strongly encourage all participants to provide as much additional data as possible that is helpful for the reproduction of the submitted results.This can be input decks for the performed simulations, source code, container images etc.The information should be provided either in form of persistent links to data (and software) published elsewhere or by uploading the data (and software) as supplementary material (maximum 100MB) following instructions given on the CSP website (https://spe.org/csp).

Technical description CSP 11B
The CSP 11B is a 2D geometry, set at field scale and conditions.No real analogue exists, but the setup is typical of the Norwegian Continental Shelf and the scenario that the downscaled experiment underlying Version 11A sought to mimic.Many of the challenges from Version 11A carry over to Version 11B, although the balance of timescales between different processes are altered as one goes from lab conditions to field conditions, see discussion in Kovscek et al. (2023).A significant simplification of 11B relative to 11A lies in the higher density and lower compressibility of the CO2rich phase that make the system easier to solve numerically.We recommend that participants start with CSP 11B.Version 11B nonetheless contains separate challenges.Firstly, thermal effects are included.Furthermore, we emphasize again the challenge of capturing the onset and development of convective mixing accurately.Moreover, the injection of a relatively cool CO2 super-critical phase will lead to potentially significant thermal effects in the near-well regions.Another challenge arises due to the boundary conditions, as Version 11B is essentially a closed system.The 2D nature of the domain, seen together with the lack of geomechanical response, therefore implies that some pressure buildup is to be expected.

Relationship to CSP 11A
The CSP 11B is a geometric scaling of CSP 11A, as illustrated in Figure 5.We summarize the commonalities and differences between CSP 11A and 11B as follows: 1) The governing equations for CSP 11A apply to 11B, with the extension to thermal effects.Thus Section 2.1 applies to CSP 11B, including the consideration of pure water (no salts).2) A thermal equation is introduced for CSP 11B, as detailed below.
3) The geometry of CSP 11A is reused in CSP 11B, with the following scaling: a) The horizontal scale (x-axis) is scaled 1:3000.b) The vertical scale (z-axis) is scaled 1:1000.Thus, the overall dimensions of the CSP 11B is a vertical cross-section measuring 8.4km horizontally and 1.2km vertically.We emphasize that we keep the coordinate system oriented with the vertical direction pointing "up" and the origin in the "lower-left" corner of the domain.4) The geometric scaling applies also to the definition of Boxes A, B and C, and the well placements.5) The only exception to the pure scaling of the geometry are the injection wells, which are kept circular, with a radius of 0.15m.6) As with CSP 11A, we assign a nominal depth to Version 11B of 1m to allow us to work with volumetric quantities.7) The initial and boundary conditions, together with the injection schedule, are updated to be consistent with field conditions, see specification below.
8) Facies properties are updated to be representative of field conditions, see specification below.
We define the three boxes in terms of their bottom left and top right corners, stated as (, ) (measured relative to the lower left corner of the domain): Box A: Bottom left (3300, 0), top right (8300, 600) Box B: Bottom left (100, 600), top right (3300, 1200) Box C: Bottom left (3300, 100), top right (7800, 400) Note that in contrast to CSP 11A, Box A and Box B do not extend to the boundary.The reason for this will be clear in Section 3.3.
The two injection wells have (, ) coordinates: Well 1: (2700, 300) Well 2: (5100, 700) The two pressure observation points have (, ) coordinates: POP 1: (4500, 500) POP 2: (5100, 1100) Furthermore, to avoid ambiguity, we define the following quantity: One year is defined as exactly 365 days, containing 31,536,000 seconds.In this equation, we denote the solid phase by  = .As convention set the "saturation" of the solid phase such that   = 1 − , and its "Darcy flux" as zero  α = 0. Furthermore, in Eq. (3.1)   is the internal energy per mass, while   is the thermal conductivity.As in Eq. (2.2),  denotes reservoir volume per domain volume, and will be detailed in Section 3.3.

Thermal equations
Thermodynamics: The general thermodynamics defined in Section 2.1 also apply to CSP 11B, with the addendum that the internal energy and thermal conductivity of the wetting and non-wetting phase is considered independent of composition, i.e.,   =   (, ), where the dependencies follow the same reference (NIST) as given in Section 2.1.For the solid phase, the internal energy is given in terms of temperature only, and is assumed to correspond to a constant specific heat capacity, i.e.: The specific heat capacity and the thermal conductivity are given in section 3.4.
We remark that since we consider the internal energy of each phase as a function of pressure and temperature, Eq. (3.1) can be linearized to provide a temperature equation.

Initial and boundary conditions
We define initial and boundary conditions consistent with a depth of about 2000 m to the top of the defined geometry.Based on a presumed geothermal gradient of 25 °C per km, we define the function   () for  = (, ) as: Boundary conditions: All boundaries are given no-flow boundary conditions for the fluid as stated in Eq. (2.12).For the energy equation we consider an insulating boundary condition for the left and right boundary: and constant temperature for the top and bottom boundaries; thus, for any point any point   on the top or bottom boundary: To avoid an unphysical increase in reservoir pressure, we introduce additional volume at the horizontal boundaries.This is a variant of what is commonly known as pore volume multipliers, that are typically implemented by giving cells at the boundary of the domain an elevated volume content.Herein, we specify these in rigorous mathematical terms as follows.
We define by ℓ  (  ) (units of meters) the volume per area of the boundary.For any  > 0, we can extend this function to the whole domain, denoted    , by using the notation   () to identify the closest point on the boundary to : In particular, the limit   0 () = lim ϵ→0    () implies that   0 is a Dirac-type distribution on the boundary with weight ℓ  , where the extra volume is interpreted to lie "just inside" the boundary.We now formally augment the volume of the domain to include this extra volume, i.e.: This has the effect of introducing an extra volume per area of ℓ  immediately inside the boundary.The effect is that the boundary has the capacity to buffer fluids and energy, while all other aspects of the governing equations remain essentially unchanged.Importantly, the definition of the effective reservoir volume given in Eq. (3.6) can be approximated by using any finite value of , so that () ≈ 1 +    (), and then implemented using a single layer of grid cells with this width.
The boundary volumes per area are specified as • ℓ  (  ) = 5 ⋅ 10 4 m for   on the left and right boundaries, within facies 2 to 5.
• ℓ  (  ) = 0m for   on the left and right boundaries, within facies 1 and 7.
• ℓ  (  ) = 0m for   on the top and bottom boundaries.
Additionally, the two injection wells ( = 1,2) are equipped with fixed temperature boundary conditions during injection, thereafter zero heat transfer, i.e., for   ∈ Γ  : Initial conditions: The CSP 11B is initialized at  = −3.1536⋅ 10 10 s (1000 years before injection begins).The initial condition is given by an initially stagnant water-filled reservoir following the geothermal gradient: To make the initial condition well-posed, we specify a pressure at the center of Well 1 of 3.0 ⋅ 10 7 Pa.
Due to the possibility of natural thermal convection, this initial condition is discussed in more detail in Section 5.5.

Facies properties
The geometry contains seven facies, six permeable and one impermeable.These provide the definition of the material properties, as given in Tables 3.1 and 3.2.The material properties are synthetic, and are chosen to be characteristic of storage reservoirs at the Norwegian Continental Shelf.

Measurables
The same measurables as for CSP 11A (defined in Section 2.5) also apply to CSP 11B.In addition, we request CO2 accumulated in boundary volumes:

CO2 in boundary volumes (SPE 11B and 11C only)
The total mass of CO2 in all boundary volumes (CO2 in any form within the region where   0 ≠ 0 in the sense defined in Section 3.3) shall be reported in kg.

Data reporting
Data reporting shall follow the same structure for CSP 11A (detailed in Section 2.6), with the following changes: 1) Sparse data: Shall be reported at 3.1536 ⋅ 10 6 s intervals (10 data points per year) starting at the injection in a file spe11b_time_series.csv, on the same form as in Section 2.6.1, with the addition of a final column with values of boundary CO2 (defined in Section 3.6.1).2) Dense data: A spatial map of all field variables (pressure, saturation, phase composition and temperature) shall be reported for every five years from injection start.The spatial maps shall be reported on a uniform Cartesian grid of 840 by 120 cells (10 m by 10 m grid cells from the bottom of the domain).For each temporal snapshot indicated by Y years, ( = 0, 5, 10, …), cell values should be provided in csv format in a file spe11b_spatial_map_<Y>y.csv on the same form as in Section 2.6.2, with the addition of a final column with temperature values in Celsius.
3) The performance data shall be reported with the same change in temporal and spatial resolution as indicated in points 1) and 2) above.

Technical description CSP 11C
The CSP 11C is a 3D geometry, set at field scale and conditions, typical of the Norwegian Continental Shelf.Version 11C, due to its 3D structure, induces significant computational overhead relative to 11B.While for 11A and 11B it is feasible to consider grids that fully resolve the development of convective mixing, this may be impossible for 11C.Furthermore, the existing challenges associated with correctly capturing the dynamics near the gas-water contact will be exacerbated.On the other hand, the 3D nature of the domain accommodates the injection volume better, and the pressure buildup seen in Version 11B should be lessened.

Relationship to CSP 11B
The CSP 11C is an arched 3D extrusion of CSP 11B, as illustrated in Figure 6.We summarize the commonalities and differences between CSP 11B and 11C as follows: 1) The governing equations and facies properties for CSP 11B apply to 11C.
2) The geometry of CSP 11B is reused in CSP 11C, with the following modification: The depth (along the y-axis) is extended to 5000 meters and deformed such that the top surface forms a parabola elevating the central part of the domain by 150m.See section 4.2 for precise definition.
We emphasize that we keep the coordinate system oriented with the vertical direction pointing "up" and the origin in the "lower-left" corner of the domain, using the convention that the two horizontal axes are enumerated first, i.e., with reference to Figure 5 coordinate triplets are given as (, , ).
3) The geometric scaling affects the definition of Boxes A, B and C, and the well placements, as detailed in Section 4.2.4) The initial and boundary conditions, together with the injection schedule, are updated to be consistent with the 3D extension; see specification below.

Definition of geometry, well placement and reporting boxes
For a precise definition of the geometry, we consider a reference and a physical configuration, with a mapping between them.Thus, the domain has a baseline gradient of 2 meters per km in the -direction, so that the back boundary (at  =  = 5000) is 10 meters higher than the front boundary (at  =  = 0).Furthermore, the domain is arched in the -direction following a parabolic shape with maximum elevation difference (relative to the baseline gradient) of 150 meters at the central ridge of the domain ( =  = 2500).
The absolute maximum elevation of the ridge is at  =  = 2541 + 2/3.
We note that Φ(, , ) is an invertible function for the domain of interest  = Φ(), and we denote the inverse as Ψ = Φ −1 .For any physical coordinate triplet  = (, , ) ∈  ⊂ ℝ 3 , we can thus recover the position in reference configuration as: The two injection wells are given as follows.Injection well 1 is considered as straight (horizontal) in physical space, and is open/perforated between the near and far points with (, , ) coordinates: Well 1: Near end (2700, 1000, 300), far end (2700, 4000, 300) Injection well 2 is considered as curved following the geology in physical space, and is thus a straight well in reference space, which is open/perforated between the near and far points with (, , ) coordinates: Well 2: Near end (5100, 1000, 700), far end (5100, 4000, 700) Note in particular that this implies that Well 2 is slightly longer than Well 1, which must be taken into account as the injection rates are given per well-length.

Initial and boundary conditions
We define initial and boundary conditions based on a depth of about 2000 m to the top of the defined geometry.
Boundary conditions: The left, right, front, back, boundaries are defined to have boundary conditions as defined in Section 3.3 for the left and right boundaries.The top and bottom boundaries are defined to have the same boundary conditions as defined in Section 3.3.
We point out that as the top and bottom boundaries are non-planar, the temperature boundary condition given in Eq. (3.3) implies that the temperature on the top boundary varies between a minimum temperature of about 36.12 °C at the top of the arch and a maximum of 40 °C where the top boundary meets the front.Similar variation applies to the bottom boundary.
The two injection wells ( = 1,2) are equipped with fixed temperature boundary conditions during injection, thereafter zero heat transfer, as stated in Eq. (3.7).
Initial conditions: The CSP 11C is initialized at  = −3.1536⋅ 10 10 s (1000 years before injection).The initial condition is given by an initially stagnant water-filled reservoir following the geothermal gradient, as stated in Eq. (3.8).
To make the initial condition well-posed, we specify a pressure at the centerline of Well 1 of 3.0 ⋅ 10 7 Pa.

Facies properties
The geological model contains seven facies; six permeable and one impermeable.These have properties as given in Tables 3 and 4, and the entry pressure follows the scaling given in Eq. (3.10).
From the horizontal permeability, the full facies permeability is defined based on a 10:1 horizontal to vertical anisotropy ratio in the reference configuration as: The permeability at a point  = Φ() in the physical configuration is then given by the standard transformation rules as: This definition of the permeability ensures that the anisotropy follows the layering in the y-direction but intersects the layering in the -direction.This captures, albeit only conceptually, the computational difficulties associated with both the geological situations of large-scale deformation and erosion surfaces.

Measurables
The same measurables as for CSP 11B (defined in Section 3.6) also apply to CSP 11C.

Data reporting
Data reporting shall follow the same structure for CSP 11B (detailed in Section 3.7), with the following changes: 4) Sparse data: Shall be reported at 3.1536 ⋅ 10 6 s intervals (10 data points per year) starting at the injection in a file spe11c_time_series.csv.5) Dense data: A spatial map of all field variables (pressure, saturation, phase composition and temperature) shall be reported for every five years from injection start.The spatial maps shall be reported on a uniform Cartesian grid relative to the reference configuration, of 168 by 100 by 120 cells (50 m by 50 m by 10 m grid cells from the bottom of the domain).For each temporal snapshot indicated by Y years,  = { 0, 5,10,15,20,25,30,35,40,45,50,75,100,150,200,250,300,350,400,450,500,600,700,800,900, 1000 } (26 reporting times) cell values should be provided in csv format in a file spe11c_spatial_map_<Y>y.csv on the same form as in Section 2.6, with the addition of a final column with temperature values in Celsius.6) The performance data shall be reported with the same change in temporal and spatial resolution as indicated in points 1) and 2) above.

Remarks
We do not expect that it is possible at the time of issuing this CSP to solve any of the three versions of this CSP in the classical mathematical sense (grid-converged numerical approximation within computable error bounds).Moreover, based on our experience, we expect that even getting the various versions of the CSP to run will require careful choices in terms of how to approximate the spatial variable (i.e., grids), the constitutive functions, and parameters in the computational algorithms.In this section, we highlight some experiences that may be useful to those attempting to simulate the CSP for the first time.

Considerations for version 11A
Several quantitative observations can be made that guide the numerical approximation choices, as illustrated by the results and discussion presented in the numerical studies of the physical experiment described in the introduction to Section 2 (see e.g., Flemisch et al (2023), Salo-Salgado et al (2023), Green et al (2023) and Wapperom et al (2023)).Some of these are: • Grid resolution: A large part of the domain will not see substantial gradients in the variables and can be represented by relatively coarse grids.On the other hand, the gas-water interface (Figure 7, sub-figure B) is quite sharp, and will require a very fine grid to resolve accurately.
Similarly, the onset of convective mixing (Figure 7, sub-figure D) happens at a cm-scale, suggesting a grid resolution on the order of millimeters to resolve accurately.• Constitutive functions: The capillary pressure in the coarser sands (Facies 2-5) is relatively small, and gravity segregation dominates even during the injection phase.Thus, the saturation is almost binary (discontinuous).If these discontinuities are not resolved by the grid, it may be advantageous to upscale the constitutive functions to the scale of the grid (likely leading to less non-linear relative permeabilities) • Maximum capillary pressure: Gas-phase CO2 does not enter any of the seal facies (Facies 1).
Considering the elevation of the spill-point, the gas column thus does not exceed 10-12 cm height.Correspondingly, it is not expected that capillary pressures exceeding 10-12 cm water column (1000-1200 Pascal) will appear in the domain.As such, modifying the capillary pressure for values above this value should not change substantially the solution (e.g., we expect that using a value of  cap,max = 2500 Pascal will not substantially alter the solution).
• Pressure boundary: Because the CO2 is in gas phase the compressibility is high, and the volume is therefore quite sensitive to pressure variations.As such, care should be taken if approximating the pressure boundary conditions.

Considerations for Versions 11B and 11C
Additional to the considerations presented in Sections 3 and 4, some more detailed comments can be made based on theoretical considerations.
• Grid resolution: Gravity-driven convective mixing is also expected to be important for 11B and 11C.Using theoretical estimates (see e.g., Riaz et al., 2006 andElenius et al. 2012) and parameters from Facies 5, we expect an onset time of a few years, with critical wavelength of less than five meters.This indicates that a grid-converged solution may require grid cells in the order of a meter or less.On the other hand, suitable scale separation may exist for upscaled models of convection to be applicable for coarser grids.• Maximum capillary pressure: As for 11A, there is a finite limit on the possible height of a gas column, and therefore it is likely that a significantly lower value of cut-off value for the capillary pressure can be used without altering the solution.
• Bottom boundary: The lowest facies is impermeable; however, it does still conduct heat.This leads to two alternatives: Either including it as part of the domain or asserting the boundary conditions above the facies.In this latter case, care must be taken to implement thermal boundary conditions that are consistent.

Dispersion
If coarse grids are applied, it is commonly considered that the numerical dispersion will dominate the physical diffusion and dispersion, and that the terms related to physical diffusion and dispersion are not necessary to include explicitly.
On the other hand, physical diffusion is needed when considering grid convergence studies, otherwise there is no lowest wavelength for convective mixing, and a grid-converged solution cannot be computed.
The threshold between what is considered "coarse" and "fine" depends on the particular numerical method, but as a rule of thumb this is proportional to the dispersivity (about 1cm for 11A, and about 10m for 11B and 11C).

Vertical equilibrium pressure and saturation distribution
Capillary barriers are important in excluding gas from the seal facies and also for the migration patterns in the upper reservoir where the facies are fining upwards.To give some intuition to this for those who are unfamiliar with capillary barriers in these systems, we exemplify the typical pressure and saturation under a capillary barrier for a 1D system in equilibrium in Figure 8.This example is based on the exact solution to the hydrostatic system when density variations are negligible.As can be seen, the saturation is nearly independent of vertical elevation near the discontinuity when the lower material properties correspond to Facies 4 and 5.
We note that the sharpness in transition of the saturation depends completely on the shape of the capillary pressure function and must be seen in relation to the vertical extent.Indeed, for finer sands, such as Facies 3, we see more gradual saturation changes.
A detail plotted in Fig. 8 (left) is the "virtual" gas pressure in regions with no gas saturation.If an additional amount of gas was introduced in the system, this would shift the gas-water contact downwards, until the point where the real gas pressure immediately below the discontinuity (solid red line at negative heights) would equal the "virtual" gas pressure immediately above the discontinuity (dashed red line at positive heights).When pressure continuity is thus achieved, the capillary barrier will no longer be able to sustain the gas column, and gas will begin to migrate into the facies above the discontinuity.
While the above discussion is in the context of 11A, entry pressure effects are also present in 11B and 11C.

Initial state for Versions 11B and 11C
Both Versions 11B and 11C are initialized at 1000 years before injection, based on a static water column and constant geothermal gradient, as defined in Section 3.3.In particular, this implies that the vertical pressure at the initial state,  ,0 (), is given by a hydrostatic column satisfying: ,0 () = 3 ⋅ 10 7 Pa −  ∫   ( ,0 ( ′ ),  geo ( ′ ))  300m  ′ (5.1) Here 300 appears as the vertical elevation of Well 1, and 3 ⋅ 10 7 Pa is the specified pressure at this well.Figure 9 shows water and CO2 densities based on the initial pressure distribution and geothermal gradient.As seen in the figure, the water density is decreasing with depth due to the thermal expansion being of greater importance than the mechanical compression.The profile of water density with depth indicates that water convection cells may appear in the domain8 .By initializing the CSP 1000 years before injection, we therefore provide some time both for the thermal profile to equilibrate (in response to the facies-dependent thermal conductivity) as well as potential water convection to initiate.On the other hand, we do not expect that 1000 years is sufficient to calculate a true equilibrium state (if such exists), and the pre-injection simulation period should therefore be considered a compromise between computational cost and establishing a reasonable initial state.
Also evident from Fig. 9 is that the injected CO2 is denser than the ambient water but is significantly lighter than water at reservoir temperatures.Therefore, some injected CO2 may sink in as the nearwell region is cooled by the injection, however in general CO2 will rise and accumulate under the lowpermeable facies.
We remark that formation water in real reservoirs will in general be saline brines, leading to higher densities than those used in this CSP.The effects described in this section are therefore not necessarily expected to be common.However, prospective storage sites do exist where injection conditions for CO2 lead to densities close to or above the formation water, and the dynamics are therefore a realistic challenge for reservoir simulation.

Figure 2 :
Figure 2: Image of physical lab-scale CO2 storage experiment after about five hours, with contours of regions with mobile gas (yellow) and dissolved CO2 (green) overlain, from Fernø et al (2023).

Figure 3 :
Figure 3: A sketch of the benchmark geometry.The geometry includes an anticline (right side) where CO2 accumulation is anticipated.There are three fault-like structures in the geometry, with different permeability (two high and one low).The lower fault (left side) is heterogeneous and consists of several facies.The upper-left fault is a homogeneous sealing fault.The upper-right fault is homogeneous consisting of a single facies.There are two CO2 injection wells (blue circles), and pressure observation points (POPs, red circles) with locations given in the text.Box A, B and C indicated in the figure correspond to regions of interest motivated in Section 1.2, and further detailed in Section 2.5, with positions defined in the text.The facies are identified by colors, and further detailed in Section 2.3.

Figure 4 :
Figure 4: Left: Capillary pressure for Facies 1, 3, 4 and 5, as function of wetting saturation.Right: Relative perm and fractional flow for Facies 5, also as function of wetting saturation.

Figure 5 :
Figure 5: A sketch of the benchmark geometry for CSP 11B, which is a scaling of the geometry in CSP 11A.For a detailed description, please refer to the caption of Figure 2.

Figure 6 :
Figure 6: A sketch of the benchmark geometry for CSP 11C.For a detailed description, please refer to the text below and to caption of Figure 2.

Figure 7 :
Figure 7: A sketch of the benchmark geometry for CSP 11C in the reference configuration.
The facies distribution, as well as the location of boundaries, for positions  ∈  in the physical configuration can now be obtained based on the facies distribution, and boundaries, of the reference configuration and the mapping Ψ.The boxes are adjusted to account for the 3D configuration.We state the (, , ) coordinates of the three boxes (measured relative to the lower left corner of the domain) in terms near bottom left and far top right corners: Box A: Near bottom left (3300, 0, 0), far top right (8300, 5000, 750) Box B: Near bottom left (100, 0, 750), far top right (3300, 5000, 1350) Box C: Near bottom left (3300, 0, 250), far top right(7800, 5000, 550)    Note that the definitions of Boxes A and B extend outside the bottom (resp.top) boundary, and the actual boxes can be truncated to the part within the physical domain .

Figure 8 :
Figure 8: Pressure and saturation distribution near hypothetical discontinuities between Facies 1 (above discontinuity) and Facies 3-5 (below the discontinuity) for Version 11A, assuming equilibrium pressure and saturation.The left figure shows the phase pressures corresponding to a gas column is shown with a 10cm height, based on Facies 3 being the lower facies.Dashed lines indicate the "virtual" gas pressures obtained by adding the entry pressure to the water pressure for the two facies.The right figure shows the saturation profile if the lower facies is Facies 3, 4 or 5, respectively.The shape of the curves is exactly the same as if capillary pressure be plotted instead of height, as can be seen by comparison to Figure 4.

Figure 9 :
Figure 9: Density of water and CO2 for the initial state of CSP 11B.The CO2 density is plotted both at reservoir temperature (based on the geothermal gradient) and the injection temperature (10 °C).

•
March 29, 2023: Official announcement of the 11 th SPE CSP at the 2023 SPE Reservoir Simulation Conference, Galveston, Texas.• October 1, 2023: Final date for publication of corrections or amendments to the CSP description (this document).• October 16-18, 2023: Special session at SPE Annual Technical Conference and Exhibition (ATCE).• December 1, 2023: Open call for participation period ends.• March 1, 2024: Deadline for submission of early CSP simulation results.• March 2024: First intercomparison workshop for all CSP participants (virtual).• September 1, 2024: Deadline for submission of final CSP simulation results.• September 2024: Final intercomparison workshop for all CSP participants (hybrid).• December 2024: Completion of draft report on the results of the CSP.• February 2025: Report on the results of the CSP finalized and submitted.• March 2025: Special session at the 2025 SPE Reservoir Simulation Conference.

Table 2
.1: Properties for 11A that vary between facies.

Table 2 .2: Properties
for 11A that are equal in all facies.
2, …, cell values should be provided in csv format in a file spe11a_spatial_map_<H>h.csv of the form # x [m], z [m], pressure [Pa], gas saturation [-], mass fraction of CO2 in liquid [-], mass fraction of H20 in vapor [-], phase mass density gas [kg/m3], phase mass • max_norm_res: Maximum normalized residual within the 1cm x 1cm reporting box at convergence of latest time-step, defined for each component  as namely,the maximum (over all computational cells  inside the reporting box) absolute value of the residual  , (the discrete version of the summand in Eq. (2.2)) divided by the local pore volume   , multiplied by the time step , and divided by the surface density   of the respective component (CO2/H2O).If this quantity cannot easily be computed for your numerical method, please specify how you define a similar measure in the questionnaire.

Table 3 . 2 :
Facies properties for 11B and 11C that are equal in all facies.