Estimating Caprock Impact on CO2 Migration in the Gassum Formation Using 2D Seismic Line Data

Realizable CO2 storage potential for saline formations without closed lateral boundaries depends on the combined effects of physical and chemical trapping mechanisms to prevent long-term migration out of the defined storage area. One such mechanism is the topography of the caprock surface, which may retain CO2 in structural pockets along the migration path. Past theoretical and modeling studies suggest that even traps too small to be accurately described by seismic data may play a significant role. In this study, we use real but scarce seismic data from the Gassum Formation of the Norwegian Continental shelf to estimate the impact of topographical features of the top seal in limiting CO2 migration. We seek to estimate the amount of macro- and sub-scale trapping potential of the formation based on a few dozen interpreted 2D seismic lines and identified faults. We generate multiple high-resolution realizations of the top surface, constructed to be faithful to both large-scale topography and small-scale statistical properties. The structural trapping and plume retardation potential of these top surfaces is subsequently estimated using spill-point (static) analysis and dynamical flow simulation. By applying these techniques on a large ensemble of top surface realizations generated using a combination of stochastic realizations and systematic variation of key model parameters, we explore the range of possible impacts on plume advancement, physical trapping and migration direction. The stochastic analysis of trapping capacity and retardation efficiency in statistically generated, sub-seismic resolution features may also be applied for surfaces generated from 3D data.


Introduction
Geological sequestration of CO 2 is essential in achieving rapid and efficient reductions in anthropogenic greenhouse gas emissions, as recommended by the IPCC. Deep saline aquifers are considered to hold the largest potential storage volumes by far globally, while depleted oil and gas fields may add additional geological storage potential (Benson et al. 2005). Saline aquifers are not otherwise considered resources (i.e., uneconomical), and therefore there are generally much less data available compared to subsurface areas with prospective or depleted hydrocarbon resources. Thus, storage quality and security assessments in CO 2 storage reservoirs are often based on sparse data (i.e., seismic surveys, geophysical well logs, rock samples).
Nevertheless, practical storage of CO 2 in saline aquifers has already been demonstrated at industrial scale, both onshore (Duong et al. 2019) and offshore (Ringrose 2018). A large share of brine-filled reservoirs relevant for CO 2 storage are sloping, (semi-)open aquifers. These are inclined, porous and permeable reservoir layers (e.g., sandstone), overlain by impermeable, confining units (e.g., shale), but with no lateral confinement, meaning that injected CO 2 can migrate gradually up dip within the reservoir along the caprock. Several of the saline aquifers off the Norwegian coast are sloping and with open or semi-open boundaries (Riis and Halland 2014). In some settings, there could be a risk for buoyant CO 2 to eventually reach the surface or seabed, at least in principle. However, the long-term risk of leakage is likely low in practical terms, both because CO 2 movement is slow relative to the long time scales involved, and because physical and chemical trapping mechanisms gradually immobilize CO 2 and retard the plume migration. In microscale (mm-m), droplets of CO 2 will be capillary trapped at narrow pore throats, some will dissolve and dissociate in water, and may eventually even mineralize, ensuring absolute immobilization of CO 2 in-situ (Bachu et al. 1994). Further, significant volumes of free-phase, buoyant CO 2 may be structurally trapped beneath the caprock surface due to meso-scale (ms to 10s of ms) irregularities (i.e., faults, undulations or folds). Predicting how this trapping mechanism may prevent CO 2 leakage or migration outside of a defined storage license is essential in risk assessment for sloping aquifer storage prospects. Open or semi-open boundary conditions may reduce the risk of pressure buildup during CO 2 injection, which is advantageous both in terms of injectivity, storage capacity and reservoir integrity.
The goal of the CO 2 Upslope project has been to investigate methods for risk assessment of sloping, (semi-)open saline aquifers where available data are scarce for practical storage of CO 2 , using the Gassum Formation as a case study. The Gassum Formation is a sandstone of Late Triassic age and vast aerial extent offshore Norway and Denmark (Skagerrak) and onshore Denmark (Nielsen 2003), see Fig. 1. Reservoir grade, brine-filled layers in the Norwegian sector are located at burial depths down to ca. 3500 m below the seafloor, with an upwards sloping trend toward the north (Baig et al. 2013;Olivarius et al. 2019). The overlying Fjerritslev Formation is sealing, comprising a thick succession of mudstones and shales (Nielsen 2003). The Gassum saline aquifer is considered to be a promising candidate for CO 2 storage (Bergmo et al. 2017), although available data remain limited. While the CO 2 Upslope project has considered how several trapping mechanisms are limiting long-time migration, as well as investigated mechanisms behind near-well salt precipitation (Parvin et al. 2020), the study presented in this article focuses on the potential for faulted and undulating caprock to slow down the advancement of a CO 2 plume and trap CO 2 . This is relevant for proving low risk of leakage to the seabed. The main challenge is data scarcity, with interpretations of top reservoir and cap rock topography based solely on spatially spread-out 2D seismic line data.
The effect of caprock topography on CO 2 migration and trapping has been extensively studied in the past, both in terms of large-scale structural/stratigraphic trapping (Riis and Halland 2014;Nilsen et al. 2015a, b;Goater et al. 2013) and in terms of impact of smallscale undulations ("rugosity") on plume advancement. Both conceptual (Gasda et al. 2012(Gasda et al. , 2013Nilsen et al. 2016) and synthetic Shariatipour et al. 2016) studies have shown how caprock rugosity has the practical effect of slowing down the advancing plume front. In the conceptual framework of vertical equilibrium modeling (cf. Sect. 3.2), the upscaled effect is shown to be analogous to that of a modified CO 2 relative permeability curve, reducing propagation speed for thin plume tips (Gasda et al. 2013;Nilsen et al. 2016).
Studies of the impact of caprock topography on plume migration for a hypothetical storage scenario in the Gassum Formation is complicated by the lack of 3D data. Contrary to other well-studied reservoirs such as Sleipner (Singh et al. 2010;Ringrose 2018), for which there exist a large amount of 3D seismic and other monitoring data (Equinor: Sleipner 4d seismic dataset 2019; Equinor: Sleipner 2019, benchmark model 2019), the data coverage concerning the Gassum Formation in the Norwegian sector of Skagerrak is relatively scarce, including some few 2D seismic surveys and one research well confirming the top formation. There is no hydrocarbon production in the study area, and as saline aquifers lack economic potential, data scarcity is a common challenge in reservoir characterization for CO 2 storage. Whereas it is possible to interpolate a surface between the interpreted 2D lines, such a surface will generally be smoothed and lack the fine-scale structure that is nevertheless important to evaluate plume advancement speed and depletion from smallscale structural trapping (Fig. 1, middle plot).
In addition, a large number of faults with varying displacements can be inferred from the seismic lines. Some faults are sufficiently large to be correlated between several seismic lines (Fig. 1, middle plot). Most faults are smaller scale, however, and die out between far-spaced seismic lines. Little can be known regarding the lateral extent and orientation of these faults. For clarity of discussion, we will in this article make an explicit distinction between (I) large-scale faults, being the regional faults that can be correlated between multiple seismic lines and (II) intermediate scale faults, being the faults that can be identified on a single seismic line only, but still have large enough throws (displacements) to be recognized within seismic resolution. Similarly, we use the term (III) small-scale faults to refer to faults too small to be recognized as such on the seismic data. Alongside with other small-scale undulations in top surface (which may be hard to distinguish from seismic noise), the latter is considered as a component of caprock rugosity, and will be modeled separately. For the purpose of this paper, we seek to infer as much information about the Gassum top surface as we can from the seismic lines. This includes statistical analysis to recreate small-scale topographical features through stochastic simulation, combined with a statistical model of fault sizes, distributions and shapes. By varying several unknown parameters, we are able to generate a large set of detailed surfaces which remain overall consistent with the observed data. In a second phase, we numerically simulate buoyant fluid flow up along and below these surfaces, assessing the impact of small-scale undulations and faults on long-term CO 2 migration in comparison with migration along smooth surfaces. We then analyze the result by looking at the spread of several aggregate quantities, including total structural trapping potential, plume tip advancement and plume center of mass, as well as how different choice of parameters impact these aggregates. For In the flow simulation phase, several simplifying assumptions are made. In particular, only the top surface topography is considered in terms of plume retention; other trapping mechanisms are ignored. Moreover, faults are considered only in their geometric effect on the caprock topography, and we do not consider the reduction in permeability across fault planes, nor the risk of upward leakage through faults. Reservoir properties and cap rock integrity are major controls on long term CO 2 storage, and have been addressed in other studies (Olivarius et al. 2019;Nielsen 2003;Bruno et al. 2014;Springer et al. 2020). In this study, faults are considered transmissive across sand-sand contacts (Fossen and Bale 2007), and sealing along mudstone intervals (Bruno et al. 2014). While the work presented here focuses on caprock topography alone, we emphasize the importance of these other effects. Although our present study suggests that the presence of small-scale detail and a large number of small faults has the potential to slow down plume advancement by some 10-30 percent, the cumulative effect of all known trapping mechanisms and reduced fault permeability may plausibly reduce plume migration significantly more. An important riskreducing factor at Gassum is that the sloping aquifer is open, and thus hazardous pressure build-up and fracturing is unlikely (Chadwick et al. 2008).

Top Surface Modeling
The first part of our study focuses on generating top surface realizations that are consistent with the input data. The data available for the study consist of 13 seismic 2D lines from 3 different seismic surveys (IKU88, FSB88 and SKAGRE96), as well as 41 fault sticks corresponding to interpreted regional faults (intersecting more than one seismic line). Regional faults have throws of more than 50 m, and are generally striking WSW-ENE. The domain of study, a four-sided polygon with an extent of 45.7 km east-west and 50.3 km north-south, was chosen to tightly fit the region of a prospective CO 2 storage site (Baig et al. 2013;Bergmo et al. 2017) covered by the seismic lines.
Regional faults (I) were imported as planes from Gregersen et al. (2018); Olivarius et al. (2019), as shown in the middle plot in Fig. 1. Intermediate scale faults (II) were interpreted as paired (displaced) points along the seismic reflector representing an acoustic impedance contrast between the sandy reservoir and shaly seal. The reflector (line data) and faults (point data) were subsequently converted from the seismic time domain (i.e., two-way travel time in m/s) to depth, using a model for the seismic velocity distribution in the overburden by Baig et al. (2013).  Baig et al. (2013). Model area with the available 2D lines (yellow) shown in polygon (pink). Middle: Top surface of the Gassum Formation inside the model area, shown in TWT and interpolated (smoothed) between 2D seismic lines. Regional faults, as interpreted by Gregersen et al. (2018) and displacing the top Gassum reservoir are displayed as intersecting planes. Note that small scale undulations/faults are only represented along line data traces. Bottom: Example of reservoir (Gassum, yellow) and seal (Fjerritslev, gray) formation geometries shown in North-South-oriented 2D seismic line through the study area (from the SKAGRE96 survey, location marked in the regional map top left). Suggested location of CO 2 injection well shown. Buoyant fluid will migrate upwards within the reservoir, and upslope along the top reservoir/base seal, gradually becoming trapped in structural traps (faults), as well as capillary and chemically immobilized along the migration path ▸ As can be seen in Fig. 2, the seismic lines criss-cross the domain, but there are large areas for which no data exists. The reflector interpreted to correspond to the top Gassum Formation is irregular in all seismic lines, and cross-cut by numerous faults, i.e., fracture zones with downward displacement, i.e., normal faults, (mostly) due to large-scale, structural extension. The seismic lines can be seen on the right plot of Fig. 2, as well as in Fig. 3, where a single seismic line from the FSB88 survey is plotted in profile.
The software algorithms used for top surface modeling in this section (spline surface approximation, local correction, stochastic simulation and fault generation) have all been implemented in MRST, and are scheduled to be included as a separate module in a future release, although this has not yet happened at the time of writing of this article.

Constructing the Trend Surface
The first step is to generate a trend surface, which estimates the overall shape of the top seal topography. To this end, we use bi-cubic B-spline tensor product penalized least square surface approximation (Floater 1998;von Golitschek and Schumaker 2002). To specify the spline surface, we use a grid of 40 × 40 control points (degrees of freedom) and a thin-plate smoothing term weight of 10 −8 . This choice of parameters allows us to define a surface that provides a good approximation of the input data without introducing the smallscale detail locally around the seismic lines that is seen in Fig. 1. As input data, we used all points on the seismic lines (a total of 25,703 points) after depth conversion, as well as a set of densely sampled points along the interpolated regional faults (18600 points). The resulting surface is shown on the top plot of Fig. 4. While the surface shown in the upper plot of Fig. 4 generally provides a good approximation to the input data while avoiding local artifacts around the seismic lines, it does not resolve the sharp edge and discontinuities around regional faults. In order to include this, we generate a correction spline surface where only the sampled points from the fault polygons are used as input data. Instead of least square interpolation, we here use a multilevel B-spline algorithm better suited to local corrections (Lee et al. 1997), where the approximating surface remains constant (zero) away from input data. After adding the correction surface to the trend surface, we obtain the surface shown on the bottom plot of Fig. 4, where regional fault discontinuities are much better resolved.

Statistical Description of Small-Scale Detail
Since we in this study are interested in the impact of detailed caprock topography on CO 2 migration, we need to add smaller-scale detail and a representation of the intermediate scale faults (II) which cannot be correlated between seismic lines. A smooth top reservoir/ seal surface is not representative of reality. The frequency, size and orientation of faults, folds or facies features (e.g., channels) should form an integral part of risk assessment, also in areas with sparse data, as the caprock topography can accelerate or decelerate plume propagation, and direct preferential migration pathways. We seek to generate plausible representations based on statistical analysis of the input data.
In Fig. 5, we compare the seismic lines to the corresponding lines drawn from the trend surface. As can be seen from the upper plot, the trend surface follows the Gassum seismic reflectors closely, but does not capture small-scale detail. The residuals are shown on the lower plot, where we note a certain spatial coherence. We consider these small-scale variations to be a combination of small-scale faults (III), other small-scale variations and seismic noise, without making any attempt to separately identify these components. We seek to extend these variations ("rugosity") to the whole surface, so that any arbitrarily sampled seismic line would exhibit similar small-scale variations as those seen on the figure. On the other hand, we are not particularly concerned about exact interpolation of the input data, since it is the statistical nature of the small-scale topography that impacts CO 2 migration. We use the residuals seen in Fig. 5 to compute a semivariogram (Fig. 6, left), which we can subsequently use to generate stochastic realizations of a 2D Gaussian random field with the same spatial correlations and mean deviations as those observed in the input. One such realization is shown on the right plot of Fig. 6.
By adding a realization of this Gaussian random field to the trend surface, we obtain a surface like the one seen in Fig. 7. The small-scale detail in the seismic lines thus now extends to the whole surface.

Physical Fault Model
In order to account for the roughness induced by the large number of intermediate scale faults (II) identified in the top reservoir reflector (cf. Fig. 2), we generate a number of synthetic faults and compute their impact on top surface shape. Since little can be Trend surface generated by least square bi-cubic spline approximation of the seismic lines and densely sampled fault polygons. The seismic lines are also shown in the plots. Top: basic trend surface with no local correction for regional fault edges; Bottom: trend surface with local correction for regional fault edges. The regional fault polygons (black) are superposed known about each individual fault, we seek to estimate the total number and size distribution of faults using statistical analysis.
To generate synthetic faults, we use a conceptual model proposed in Kim and Sanderson (2005), illustrated in Figs. 8 and 9. In this model, the fault slip surface is described as an ellipse with major axis length L (fault length, i.e., strike slip component) and  where r is relative distance from fault center (Fig. 8c, d).
In general, the Gassum top surface will not intersect the fault directly at the fault center, but at some arbitrary vertical offset h (Fig. 8a). We consider that the fault will impact the geometry of the top surface up to a lateral distance D away from the fault plane intersection (Figs. 8b,9). Before using this model to stochastically generate a plausible set of faults intersecting the Gassum top surface, several parameters must be chosen: where we use n = 1 (self-similarity) and c = 0.016 , roughly in line with the statistical data presented in Kim and Sanderson (2005) for normal faults. As this is a highly uncertain parameter (Kim and Sanderson 2005;Schultz et al. 2013), we also perform a simple sensitivity analysis when running flow simulations in the second part of the paper. Sensitivity analyses are also used for item 4 and 5. For item 5 (fault strike), we use the general direction of the regional faults, which we estimate from the data to be N62 • E as the base case. For item 4, we have limited information, but based on the known vertical and plane parallel displacements, we test a range of ratios D/L from 1 to 4. For fault dip, we use an average dip angle of 75 • . Since the overall slope of the trend surface is about 14 • , this corresponds to an average angle of 61 • between fault and trend surfaces.

Statistical Model for Fault Size Distribution
As explained in the previous section, we consider fault displacement max and L to be proportional. We here assume that maximum displacement length max can be modeled as a random variable that follows an exponential probability distribution with unknown parameter . Its probability density function is written: Although there is no apriori physical argument for using the exponential distribution here, we will note that after estimating that this choice leads to a curve that fits the observed data quite well (c.f. Fig. 10).
In "A", we provide detail on the derivation of formulas and how we use maximum likelihood estimation (MLE) to estimate the unknown parameter from the observed data. By assuming that a seismic line may intersect a fault surface at any arbitrary location, and that a line is more likely to intersect with a larger fault than a smaller one, we end up with the following formula for the probability distribution of measured fault displacement length interpreted from seismic reflection line data: with: However, as noted in Kim and Sanderson (2005), limitations in the resolution of seismic data can cause observed displacements or seismic horizons below a threshold of 10 m or more to remain undetected. This is consistent with the data in our study, where the frequency of observed displacements suddenly drops below 9m, as seen in Fig. 10. For the purpose of estimating , we therefore work with a truncated probability distribution and exclude samples below the threshold min = 9 m. Using * measured to represent "observable" measured displacements larger than min , we obtain the probability distribution plotted in red in Fig. 10, and defined as: where we use F measured to denote the cumulative distribution function (CDF) of measured . Using MLE, we finally estimate to be ̂= 0.0584m −1 . (

Estimating Total Number of Intersecting Faults
In "B" we show that the number of intermediate (II) faults per unit area, N s , can be estimated as: where n l is the number of seismic lines, B l is the number of fault intersections per unit length for seismic line l, ̂ as estimated above, and l the angle between seismic line l and the general fault strike direction. A minor fault is counted as within a unit area if its center is located inside that area. Since we cannot detect all fault intersections for displacements smaller than min , we only count observed fault displacements larger than this threshold when computing B l , call the resulting estimate N measured , and divide its value by 1 − F measured ( min ) (cf. end of Sect. 2.4) to get an estimate of the total number of minor faults per unit area (whether or not within the seismic resolution): This procedure, when applied on our input data, gives us the estimate N s = 1.945 × 10 −6 minor faults/m 2 .

Adding Synthetic Fault Realizations to the Base Surface
Once values have been chosen for the variables in the list 1-6 in Sect. 2.4, stochastic minor fault fields can be generated. One such realization is shown in Fig. 11, generated with parameter values summarized in Table 1. The total number of fault/surface intersections within the rectangular bounding box of the model surface is obtained by multiplying N s by total area ( 3.4 × 10 3 km 2 ), giving a value of 6650 intersections (many of which resulting in displacements below the seismic resolution).
It should be noted that due to lack of more specific information, the spatial distribution of faults was assumed to be uniform. However, one might expect that in reality there could be a clustering of minor faults (i.e., damage zone) around the intermediate (II) and large regional (I) faults.

Impact on CO 2 Trapping and Plume Migration
Having established a procedure to create top surface realizations with small-scale irregularities and synthetic faults, we now seek to study the impact of topography on total structural trapping capacity and long-term retardation during upslope migration of CO 2 injected into the aquifer. Although key parameters (semivariogram, N s , ) have been chosen to be statistically consistent with the input data, uncertainties are still considerable. Our goal is thus limited to get some rough idea of what the impact could be. To this end, we vary the values of c, D/L ratio and fault strike and create a large number of realizations subsequently used for analysis. We define some quantitative measures, and look at the impact of different parameters on their means and standard deviations. Our reference set of parameters is that of Table 1, which summarizes the estimates and choices made in Sect. 2. We then introduce alternative values for c, fault strike and D/L-ratio as listed in Table 2. These are used to construct and analyze alternative sets of top surface realizations.

Impact on Structural Trapping Capacity
By "structural trapping capacity," we refer to the total volume of topographical traps for a given top surface, regardless of trap size, and regardless of whether these traps would actually be reached by the CO 2 plume for a particular migration scenario. In order to identify the traps and compute the total trap volumes for a given reservoir top surface realization, we use the spill-point tool in MRST-co2lab . The result of this purely topographical analysis is shown in Fig. 12, where we have applied it three times: on the base surface, on the base surface with added small-scale rugosity but no intermediate (II) faults, and on the base surface with intermediate (II) faults added. (The regional faults are added beforehand as Fig. 11 Left: stochastically generated fault field using parameters of Table 1, shaded to enhance visibility of detail. Right: Generated fault field superimposed on the Gassum trend surface, also including the regional faults and small-scale detail discussed in Sect. 2.2. C.f. Fig. 7 Table 2 Table 1 was used to generate this surface realization. We see that the base surface with no added detail has only five identified traps, whereas the inclusion of small-scale variation and intermediate (II) faults adds a large number of traps of different sizes. For this particular realization, the total trap volume of the base surface is 0.40 km 3 , whereas that of the middle surface is 0.61 km 3 , and 2.1 km 3 for the right surface. As such, the small-scale variations add a modest but still significant volume of structural trapping, whereas the inclusion of intermediate (II) faults increases the total trapping volume more than five-fold. It should be noted that these figures represent the total volume of the traps; in order to get the equivalent pore volumes, these values should be multiplied by the corresponding average rock porosity (in the order of 0.22 to 0.25 for the top Gassum reservoir, Ref. Olivarius et al. (2019)). These values also depend on the particular realization, although the contribution of small-scale detail remains relatively constant. The leftmost bar in Fig. 13 indicates the variations in trap volume obtained when running 30 realizations of the top surface generated using the reference parameters. We now assess the impact of varying the parameters in Table 2. We generate another 30 realizations for each of the six alternative parameter values in the table, while holding the other parameters at their reference values (Table 1). In the second column of Fig. 12 and the left column of Fig. 14, we see the result of varying c, the ratio between maximum fault displacement max and fault length L. Since the distribution of max remains fixed, smaller values of c correspond to fewer but longer faults, c.f. Eqs. (2) and (7). We can observe this on the left column of Fig. 14, where the upper plot represents a "large" value  For the "small" value of c, the few but large traps yields a high standard deviation, whereas the "large" value of c yields a very small standard deviation-the effect of a large number of small traps means that variations tend to cancel out. As such, the contribution of the fault field becomes more like that of the small-scale variations.
In the reference case, we kept the D/L-ratio to 1, which means that a fault perturbs the top surface up to a distance from the fault equal to the fault length L. This is a somewhat artificial ratio that follows from the definition of the synthetic fault geometry. As such, one might hope that this ratio would not impact the result too much, but from the third column of Fig. 12 and the second column of Fig. 14 we see that increasing this value also increases the estimated total trap volume. This highlights that results from this type of study must be interpreted with proper caution, as particularities of the synthetic fault model can (and do) influence the result.
Another uncertain parameter is the general strike direction of intermediate faults (II). When varying this value ±45 • from its reference value, we get the results seen in the right columns of Figs. 12 and 14. It is clear that this angle plays a significant role with regards to structural trapping capacity, as is to be expected. The more this angle is aligned with the general upslope direction, the less the impact of the intermediate (II) faults on structural trapping becomes.

Impact on Migration
To study the impact of the top surface on long-term plume migration, we run numerical simulations using the vertical equilibrium (VE) flow simulation functionality in MRST-co2lab . The assumption of VE entails that the two phases (brine and CO 2 ) are considered to be in hydrostatic equilibrium in the vertical direction at all times. The known vertical distribution of fluids allows upscaled flow simulations to be run on a 2D domain, with 3D reconstruction as a post-processing step. The elimination of the vertical dimension from the upscaled flow equations drastically reduces computational requirements, allowing for practical simulations on lateral high resolution grids within reasonable computational time even on standard laptops. VE can be considered a reasonable modeling assumption when vertical fluid segregation occurs on a timescale much shorter than the  Table 2 timescale of interest, which is frequently the case when studying CO 2 migration in a saline aquifer. Several past studies have shown the suitability of VE-based simulation for such cases (Bandilla et al. 2019;Nilsen et al. 2011).
We first study how top reservoir topography impacts CO 2 migration injected at three different locations (southwest, south, and southeast within the domain), and for different selections of top surface realizations. We first run the simulations on grids where the top surface is set to the trend surface with no added detail (cf. Fig. 4). We then run the simulations again three times on grids where small-scale detail has been added, based on three different realizations of Gaussian random fields based on the semivariogram estimated in Sect. 2.2. Finally, we run the simulations again another three times on grids where three different realizations of intermediate (II) fault fields have been added, using the parameters of Table 1. We then compare the results across injection points and top surface type.
As we focus on the impact of the top surface itself, we employ simple, homogeneous rock parameters (permeability and porosity) and uniform vertical thickness throughout the simulated domain, chosen to be within the range estimated in Gregersen et al. (2018); Olivarius et al. (2019). It should be noted that in reality, the presence of a large number of small (III) and intermediate (II) faults would probably introduce significant local variations in effective parameters. Dissolution would likely also play an important role at the timescale studied, but is not included in the analysis below. Table 3 lists key simulation parameters. Brine density and viscosity are kept constant, whereas CO 2 properties are functions of local temperature and pressure, following (Span and Wagner 1996;Fenghour et al. 1998).
The shape and position of the free CO 2 plume at year 500 for the different simulations is shown in Fig. 15. From these plots, we see that adding small-scale detail (row 2-4) appears to have some limited impact on the advancement of the plume compared to the smooth base case (row 1), and the plume is visibly more irregular. The actual realization does not seem to matter much; the plumes in row 2-4 are all fairly similar. When minor faults are added in (row 5-7), the plumes become more fragmented, with quite significant variation between realizations.
In order to assess the impact quantitatively, we define the following three indicators: (1) maximum distance from injection point to plume tip; (2) distance from injection point to plume center of mass; and (3) direction of migration (direction of center of mass from injection point, with 90 • being north). We plot the value of these indicators as a function of time for all simulated scenarios. The graphs are presented in Fig. 16. From these plots, a few observations can be made.  • Small-scale detail (red curves) had a very small impact on plume advancement for the south and southeast injection points. For the southwest injection point, the scenarios with small-scale detail saw a reduction in plume advancement (whether measured at tip or at center of mass) of about 15 percent. It is worth noting that migration from the southwest migration point occurs along a narrow sloping ridge (cf. Fig. 15). With the exception of one particular realization, small-scale detail had very minor impact on migration directions for all injection points. • With the exception of the above-mentioned outlier, there was very little difference in outcome between realizations in terms of small-scale detail-the red curves match each other closely. It is not clear why one realization caused a deviation of about 4 • in migration direction. • For the top surfaces with added minor faults (blue curves), plume advancement was reduced by 0 to 30 percent depending on injection location and whether the tip or the center of mass was used as a measure. Deviations in migration direction of up to 18 • While the number of samples is low (three realizations for small-scale detail and another three for minor faults), the results suggest the impact of these features on plume migration distance and direction depends significantly on the injection location. It does not seem possible here to define any globally valid top surface retardation factor. To examine the sensitivity of the outcome to variations in the key parameters listed in Table 2 we run another set of simulations focusing on the middle well. We run a batch with 30 simulations of the reference case, for different top surface realizations. We then run six batches with 30 simulations each, where we vary one parameter at a time, for the alternative values listed in Table 2. In total, we run 210 simulations. We then plot the graphs of the indicators defined above (tip migration distance, center of mass migration distance and migration direction) in terms of batch averages and batch mean deviations, as shown in Fig. 17.
The impact of variations in c is shown in the left column on the figure. We see that the average impact on plume migration distances is not much affected, although the "small" value of c = 5.0 ⋅ 10 −3 (representing average fault lengths 3.2 times longer than and fault orientation (right column). The batch parameter value is indicated by the color, as indicated by the legends. The color black represents the reference case the reference case) show significant variation between outcomes. The migration direction, however, seems to be significantly impacted by variations in c, where longer faults leads to stronger deviations (and more spread in outcomes). Note that the orientation of the faults, N62 • E, corresponds to an angle of 28 • on these plots.
In the middle column of Fig. 17, we see that increasing the value of D/L significantly decreases migration distances and reduces the eastwards deviation of plume migration. As this is an "artificial" parameter used with our synthetic fault model, these results serve as a reminder that our study is sensitive to the specifics of our fault model and should therefore be interpreted with due caution.
The right column of the figure shows the impact of varying the average orientation of minor faults. The plots confirm that minor fault orientation may have significant influence on the plume migration direction, as well as (in the S73 • E case) migration distance.

Conclusion
In this study, we set out to assess how smaller-scale features in top surface topography may impact long-term migration of CO 2 in the Gassum aquifer, based on a set of interpreted 2D seismic lines with small-scale variations and a large number of identified minor faults. When generating top surface realizations, rather than trying to interpolate the lines exactly, we sought to create surfaces whose small-scale detail were statistically consistent with the behavior observed in the lines. With the limited amount of input data, and a number of model assumptions and highly uncertain estimates, it is not possible to say anything definite in terms of potential trapping capacity or quantified impacts on plume migration distances and orientations. Nevertheless, we believe that our results may provide some general understanding in terms of possible magnitude and behavior of these effects. In particular, our study suggests that the topographical effects of smaller-scale features and faults on plume migration depends significantly on the injection point, i.e., its location with respect to larger features in the trend surface. Consequently, it might not be possible to define universal, global "retardation factors" that predicts plume retardation from the statistical nature of smaller-scale features alone. For instance, the impact on migration from the southwest injection point (where CO 2 is funneled northeastward up a narrow ridge) is notably different from the other injection points. Our sensitivity studies also suggest that average length and orientation of minor faults can play a significant role in directing and/or slowing down CO 2 migration. It would therefore be particularly useful to be able to deduce these characteristics from monitoring data. In any case, testing the effect of intermediate scale geological features on structural trapping capacity and potential for plume retardation (e.g., fault strike perpendicular to updip, buoyant migration) or the formation of preferential flow paths (e.g., fault strike parallel with updip, buoyant migration) should be an integral part of risk assessment.
We note that our synthetic fault model introduces effects on its own. This comes from the way faults are inserted in the trend surface as a post hoc modification, which requires specification of a specific "distance of influence," D/L. An alternative would be to include minor faults directly in the same way as regional faults are defined, using local corrections based on the multilevel B-spline algorithm (Lee et al. 1997) or something similar. Computationally, this would be a costlier option and was not attempted in this study, although it may be a worthwhile idea to pursue in the future.
Finally, we would like to reemphasize that our study focuses on the impact of small-scale caprock features, and many other important factors impacting plume development have been left out. Aquifer heterogeneities, dissolution, residual trapping, loss of permeability from salt precipitation and upward leakage through faults are all factors that could likely have first-order effects on plume development. We believe, however, that by keeping the modeling simple, the impact of caprock shape can be more easily discerned and understood. function f � ( max ; ) whose exact form we will proceed to determine. The measured displacement at the intersection is thus the random variable measured = � × Y , whose probability distribution becomes: Although we assume that fault displacements and sizes follow an exponential distribution, the same cannot be said of faults intersecting with a given seismic line l. It is natural to assume that the probability of l intersecting a given fault surface will be directly proportional to the fault surface area (issue 2 in the list above). What the random variable ′ represents is therefore the maximum displacement of a fault, conditioned on the fact that it was intersected by the seismic line. To derive the corresponding probability distribution f � ( max ; ) , we apply Bayes' theorem, and write I(l) to represent the situation where the fault is intersected by l: Since our assumption is that the probability of a fault intersecting l is proportional to its area (and thus proportional to its squared length), the probability that a fault with max displacement max will intersect l can be written on the form c 2 max , where c is some constant. For an arbitrary fault with unknown max therefore have: Likewise, for the conditional probability P(I(l)| < max ) , we have: where F ( max ; ) = ∫ max 0 f ( ; )d is the cumulative distribution function (CDF) of the random variable . Considering that by definition we also have P( < max ) = F ( max ; ) , we insert (14), (15) and (9) into (13), and obtain: Since (16) is the CDF of ′ , we obtain f � ( max ; ) by differentiation: Combining (17) with (12) and (11), we obtain the complete expression for f measured ( ; ). The last adjustment we need to make before estimating is to account for the fact that fault displacements below a certain threshold might not be identified in the seismic data due to limited resolution (item 3 on the list above). This is discussed in Sect. 2.4, where the truncated probability distribution f * measured ( ; ) is introduced in Eq. 6. It is this truncated distribution we use for the MLE estimation, along with all observed fault displacements larger than the lower threshold of min = 9m . The MLE algorithm, which searches for the value of that optimizes the likelihood of the observed data, gives the estimate ̂= 0.0584m −1 .

B Estimating Total Number of Minor Faults
We here aim to show how Eq. (7) was derived. We start by assuming a seismic line L projected on a geological horizon of infinite extent. The horizon is intersected with faults of varying sizes, all having the same strike direction at angle with the seismic line (Fig. 18). We further assume that the size of the faults follow an exponential distribution (9) with smaller faults being far more frequent (they also occur associated with larger faults, in the damage zone), and we have aligned the coordinate system so that L passes through the origin and is aligned with the y axis. We aim to establish a relation between the average number N s of fault-horizon intersections per unit area and the average number N l of faults intersected by L per unit length. The first step is to establish the distribution of fault intersection lengths with the geological horizon. We use the fault model described in Sect. 2.3 and illustrated in Fig. 8. The relationship between fault length and displacement is described by Eq. (2), where we use n = 1 . Although faults sizes in the 3D domain are assumed to follow an exponential distribution, the same cannot be said for the faults that actually intersect with the geological horizon. The probability that a fault intersects with this surface is directly proportional to fault height, and we follow an argument similar to the one used in "A" and developed in Eqs. (13)-(17) to establish the following probability density function for maximum fault displacements I , conditioned on the fault intersecting the horizon: x y L θ Fig. 18 Zenith view of a seismic line l (in red) traced on a geological horizon intersected by a number of minor faults (black diagonal lines), at angle with l. We have here chosen our coordinate system to make l coincide with the y-axis As illustrated in Fig. 8a, we consider that a given fault will intersect with the horizon at some arbitrary vertical distance h from the fault center. If we use h r to represent the relative vertical distance: h r = h H∕2 , then h r ∈ [0, 1] and for an unknown fault this value can be considered a random variable that follows uniform distribution h r ∼ U[0, 1] . If L I (h r ) denotes the length of the intersection between the (elliptical) fault plane and the horizontal plane as a function of h r , we can define: and since h r ∼ U[0, 1] , we can consider I(h r ) a random variable whose CDF becomes: Since we have that L I = Ic I , with c from (2), and with probability distributions for I and I obtained from (18) and (20) respectively, we obtain the CDF of L I by integration: To derive this expression, we keep in mind that by definition P(I < cx max ) = 1 whenever cx max ≥ 1.
The second step is to estimate the average number N l of fault intersections with the seismic line L along the horizon, per unit length. We let N s denote the average number of fault intersections with the top reservoir horizon per unit area, where a fault intersection line is considered to belong to a unit surface square if its center point is contained within the square. From Fig. 18, we see that a fault intersection line with center located at (x, y) will intersect L if its length l I > 2| x sin | . This means that the expected (mean) number of faults intersecting L per unit length is expressed: On line three of the development above, the interval of integration (−∞, ∞) was changed to [− , ] where = max sin 2c , since by definition P(I > x) = 0 for x > 1. Given an estimate of N l , we can estimate N s as: Since N l can be estimated by counting observed fault intersections in each available seismic line, dividing by the line length and averaging over the number of seismic lines, we get the following estimator for N s : where n l is the number of seismic lines, B l the number of fault intersections per unit length for seismic line l, ̂ as estimated in "A" and l the angle between seismic line l and the general fault strike direction.
Funding Open access funding provided by SINTEF AS. The research presented in this article was funded under Grant 268512 in the CLIMIT-programme of the Norwegian Research Council.
Availability of data and materials Generated surface realizations from the study can be made available by the corresponding author upon reasonable request.
Code availability The software algorithms used for top surface modeling in this section (spline surface approximation, local correction, stochastic simulation and fault generation) have all been implemented in the open-source software MRST, and are scheduled to be included as a separate module in a future release. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.