Assessment of CO2 storage capacity based on sparse data: Skade Formation

Large North Sea aquifers of high quality are the likely major target for 12 Gt of European CO2 emissions that should be stored in the subsurface by 2050. This involves an upscaling of the present combined injection rate from all European projects, which requires careful examination of the storage feasibility. Many aquifers are closed or semi-closed, with storage capacity mainly constrained by consideration of caprock failure criteria. Because the induced overpressure can propagate to sensitive regions far from the injector, the risk of caprock failure must be examined in terms of large volumes and for long times. This poses challenges with respect to fluid-flow simulation in the presence of highly uncertain aquifer properties. In this work, experts on geology, geophysics, geomechanics and simulation technique collaborate to optimize the use of existing data in an efficient simulation framework. The workflow is applied to the large North Sea Skade Formation, with potential for secondary storage and pressure dissipation in the overlying Utsira Formation. Injection at three locations in Skade gives an overall practical capacity of 1–6 Gt CO2 injected over a 50-year period, depending on the reservoir permeability and compressibility. The capacity is limited by local pressure buildup around wells for the lowest estimated reservoir permeability, and otherwise by regional pressure buildup in shallow zones far from the injection sites. Local deformation of clay due to viscoelastoplastic effects do not have an impact on leakage from the aquifer, but these effects may modify the properties of clay layers within the aquifer, which reduces the risk of lateral compartmentalization. Uplift of the seafloor does not impose constraints on the capacity beyond those set by pressure buildup.


Introduction
Carbon Capture and Storage (CCS) plays a significant role for necessary CO 2 emission reductions in Europe. According to the energy roadmap prepared by the EU Commission (European Commission, 2011), emission reduction by CCS amounts to 12 Gt until 2050, out of a total of 80 Gt that should be kept out of the atmosphere in this time frame compared to "business as usual". The average storage rate from 2025 to 2050 would then be 500 Mt/year. Large offshore aquifers in the North Sea may account for much of these emission reductions. Several approaches of different complexity have been proposed when evaluating the suitability and capacity of target formations. The Norwegian Petroleum Directorate (NPD) has compiled data and performed estimates of the storage capacity of several formations in the North Sea and other parts of the Norwegian Continental Shelf in a "storage atlas" (Halland et al., 2014). The evaluations are based on assumptions on the fraction of the aquifer volume that can be used for storage, so-called storage efficiency, and is part of a world wide effort to map storage potential (e.g. U.S. DoE, 2015;Bradshaw et al., 2011;Bentham et al., 2014). In most situations, dynamic effects such as the temporal pressure buildup need to be evaluated using large-scale fluid-flow simulations. 'These simulations are challenging as they rely on availability of good input data over very large regions, but many aquifers are untested and data is therefore sparse.
T One of the largest aquifer systems in the North Sea is the Skade-Utsira aquifer, spanning over both Norwegian and UK sectors. This aquifer has been considered as one option for large-scale storage (Halland et al., 2014;Pham et al., 2013). Several studies have focused on the Sleipner storage project in the Utsira where CO 2 has been injected since 1996 (e.g. Singh et al., 2010) at a rate of approximately 1 Mt/year. Extensive monitoring and research related to this operation has contributed to improved understanding of geological CO 2 storage. Comparatively limited research has targeted the large injection rates needed to fulfill the roadmap as described above (Wangen et al., 2016;Gasda et al., 2005;Lindeberg et al., 2009), or injection into the Skade Formation below Utsira (Pham et al., 2013;Lie et al., 2016;Allen et al., 2017). In particular, these previous simulations on Skade were performed for a limited part of the formation.
To date, despite the many capacity studies that exist, very little consideration has been given to the interaction between geological characterization and numerical modeling. For instance, geological characterization should take into account the dynamic pressure and plume evolution provided by preliminary computational studies in order to identify and quantify the controlling factors. On the other hand, computational studies that fail to include accurate information about large-scale geology or realistic boundary conditions can produce unreliable predictions of dynamic pressure buildup and plume migration. This disconnect can be solved through close collaboration and interaction between geological and computational studies.
In this study, we use an integrated approach to investigate if large quantities of CO 2 can be safely stored in the Skade Formation, using the Utsira as a potential secondary storage unit. We have performed a complete reevaluation of the Skade/Utsira aquifer system that takes into account the potential for pressure migration over hundreds of kilometers, far from any potential injection well. The new geological characterization improves our understanding of the hydraulic connection between Skade and Utsira by including the proximal sands, time equivalent to Skade and Utsira, in the UK sector as part of the Utsira Formation. (The proximal sand is sometimes referred to as the Hutton Sand, a terminology choice we avoid for clarity.) In this way, we present a more reliable representation of the system boundaries. Given these improved data, we perform a comprehensive capacity evaluation that accounts for CO 2 migration, caprock integrity and injection-induced uplift. Numerical models of varying complexity were used in this evaluation. This includes two-phase flow models to examine the largescale CO 2 and pressure migration, as well as hydromechanical models to study uplift and detailed aspects of caprock integrity. Finally, we make a comparison between analytical and numerical methods for capacity evaluation in order to understand the utility of simple estimates in different geological regimes. The Lower Miocene Skade Formation consists of marine sandstones (mainly turbidites) deposited over a large area of the Viking Graben, see Fig. 1a. It extends over an area approximately between 58°N-61°N and 1°E-3°E, which is considerably larger than the previously mapped portion of the aquifer (Halland et al., 2014). In the Miocene, a deltaic system had developed from the Shetland Platform towards the  Fig. 3 . Also, exploration wells biostratigraphy data used to delineate the Skade Formation (Eidvin and Rundberg, 2007;Eidvin et al., 2013) are shown in black, and the coastlines of Norway and the Shetland Islands are outlined. The blue box is the location of a previous simulation study in the Skade/Utsira system by the Norwegian Petroleum Directorate (Pham et al., 2013). (b) Northwestsoutheast cross section showing the Oligocene to Pleistocene stratigraphy (position given in a). The vertical axis shows the "depth" in time domain. Formation tops are indicated using meters depth (md) below seal level. Source: The Norwegian Petroleum Directorate.
Norwegian sector of the North Sea, and is represented by the Skade and Utsira formations shown in Fig. 1b. These formations constitute the main part of this sandy system and are considered as a single aquifer. There has been no CO 2 storage operation in Skade and very little is known of its properties. As the need for CO 2 storage increases, storage feasibility must be more thoroughly evaluated for all prospective aquifers, including the vast majority that have not experienced injection, and thus Skade is an important and representative formation to study.
An important part of this work is dedicated to the description and reduction of data uncertainty. We focus particularly on the formation extent and depth, the stress regime, permeability and compressibility. The estimated storage capacity resulting from the fluid-flow simulations is presented in terms of a sensitivity matrix that explores the impact of the parameters with the greatest uncertainty. In this case, the most uncertainty is in permeability and compressibility.

Methods for effective storage capacity estimates
Effective storage capacity estimates have been previously evaluated for open (U.S. DoE, 2008) and closed (Zhou et al., 2008) systems. These methods do not consider practical limitations such as dynamic pressure buildup, and their use is therefore limited to screening purposes (Bachu, 2015). The effective storage capacity [kg CO 2 ] is Here, E is an efficiency factor, the ratio of potentially stored CO 2 volume to total aquifer pore volume, ρ n is CO 2 density at reservoir conditions and V f is the aquifer pore volume. Estimates of the effective storage capacity in the Skade/Utsira aquifer and many other formations in the North Sea have been reported in the North Sea storage atlas (Halland et al., 2014) and in an investigation by the Geological Survey of Norway (NGU) (Bøe et al., 2002). The NGU survey considered the part of Skade that is more than 800 m below sea level and treated it as an open system due to its connection to Utsira. A storage efficiency coefficient of 6% was applied and they arrived at M CO2 = 15 Gt in Skade, in addition to 42 Gt in Utsira. In the storage atlas, a smaller portion of the Skade-Utsira aquifer was considered suitable. Using the pore volume of this region with 4% storage efficiency, they arrived at an estimate of 16 Gt. In comparison, the range for E in closed aquifers as suggested by Bachu (2015) and Goodman et al. (2013) is only 0.3-1.2%. Here, we improve on these effective storage capacity estimates by extended aquifer characterization and examination of E. Since previous interest was focused on storage capability in the Norwegian sector of the North Sea (Halland et al., 2014), the UK (western) proximal equivalent to the Utsira was poorly constrained. It has therefore been part of our objective to assess whether the western subcrop aquifer boundary is open or closed. The investigations performed in our work assume that the system is closed or semi-closed. For closed systems, Zhou et al. (2008) provide where c T is the total compressibility (see definition in Section 3.2) and Δp is the pressure increase that the caprock can tolerate. The shallowest, most sensitive, region is used when evaluating Δp. Dissolution of CO 2 into brine can also have a positive impact on geological CO 2 storage capacity. First, dissolution reduces the leakage risk by removing CO 2 from the buoyant CO 2 phase into the brine, where it is no longer buoyant. Second, because of the overall density increase, the pressure buildup is reduced, which reduces the risk of caprock failure. Dissolution is facilitated when convective mixing occurs below the CO 2 plume (Riaz et al., 2006;). In the system of interest here, convection is enhanced due to high permeability and low salinity. The rate of dissolution F [kg/m 2 s] and the mass dissolved M diss [kg] can be estimated as, Here, k z is the vertical permeability, c CO2 is the mass concentration of CO 2 in brine at the solubility limit, w CO 2 is the density increase of brine with dissolved CO 2 , µ w the brine viscosity, and g the gravitational constant. The parameter λ accounts for interaction with the so-called capillary transition zone in the lower parts of the plume (Elenius et al., 2014). The value has previously been estimated at λ = 0.04 in Utsira (Elenius et al., 2014), and the same value is used here for Skade. In addition, estimates are needed for the density increase of brine with CO 2 relative to brine without CO 2 , w CO 2 , the brine viscosity, µ w , and the average footprint area of the CO 2 plume, A CO2 , up to the time t of interest. This estimate of dissolution does not account for the limited volume of brine accessible for CO 2 dissolution and should therefore be seen as an upper limit.

Methods for practical storage capacity estimates
As discussed by Bachu (2015), the practical storage capacity is further constrained from the effective capacity by taking into account technical, economic and regulatory considerations. Here, we examine the impact of restrictions on pressure buildup, CO 2 migration, uplift of the seafloor and small-scale deformations. The overpressure is examined over the entire aquifer system during the injection and postinjection periods and realistic bounds are determined for the pressure tolerance of the caprock. In addition, we impose an injection scenario with three well sites in Skade and injection over 50 years, which we considered to be technically and economically realistic.
Recently, Allen et al. (2017) performed simulations to study the practical storage capacity of several formations in the Norwegian Continental Shelf. Injection scenarios were chosen to maximize injected amounts subject to constraints on CO 2 migration and overpressure. The latter was constrained to 90% of the lithostatic pressure, which we consider to be too generous (Section 5.5), and restrictions on the pressure applied only in the presence of CO 2 . The storage efficiency E of Skade and Utsira was less than 1%. Earlier, Pham et al. (2013) estimated a practical capacity of 0.17 Gt in a southern segment of Utsira-Skade using 4 wells and a constant overpressure restriction of 10 bars.
We improve on previous estimates by improved mapping of the reservoir, by performing sensitivity studies related to permeability and compressibility and by using more realistic pressure tolerance for the caprock. In the following, we describe four models that were used in our work to evaluate the practical storage capacity. The models are presented in the order of increasing complexity. First, an analytical Theis model is presented that is used to produce a rough estimate of pressure buildup. Second, a two-phase flow model is used to simulate CO 2 migration and pressure buildup during the injection and post-injection periods. Third, a regional poroelastic model is used to investigate the amount of uplift that can occur during the injection period. Fourth, a viscoelastoplastic model is employed to investigate dynamic processes in clay due to its interaction with CO 2 that can increase its permeability. Modified clay properties may have an impact on flow compartmentalization and connectivity between Skade and Utsira.

Theis solution
In some cases, dynamic pressure buildup close to wells may limit the capacity. Analytical solutions to the pressure change resulting from well operations were first obtained by Theis (1935) for single-phase flow from/to one well of negligible radius (no wellbore storage) that fully penetrates the formation. It applies to laterally unbounded aquifers of uniform properties and thickness, with no seepage through the seals above and below the formation. The change in hydraulic head (Δh) caused by injection/production of rate Q is: Δh = s(Q/4πT)W(r 2 S/4Tt), in an aquifer of transmissivity T and storativity S, at a distance r from the well and time t after commencement of well operations. The exponential integral, W, can be evaluated directly or further simplified. The method is used here in addition to performing simulations. Although several extensions have been developed, such as consideration of two-phase flow (Mathias et al., 2009), we use the simplest version.

Two-phase flow model
The flow simulations were performed using a vertical equilibrium model (VESA) (Gasda et al., , 2013 with two phases, CO 2 and water. This model assumes an instantaneous gravity segregation of CO 2 and brine at the time scale of the simulation and compares well with 3D simulations in our previous work on Skade (Elenius et al., 2016) and in other comparison studies, e.g. Nordbotten et al. (2012). The simulations performed here include residual and structural trapping, while solubility and mineral trapping are neglected. The VESA model simulates the pressure and saturation change in the reservoir in time and space. Saturation is a depth-averaged quantity, which is used to reconstruct the height of the CO 2 plume. Also the porosity and permeability are depthaveraged properties. The relative permeability for each phase is a saturation and permeability weighted average. It is especially important to properly account for permeability differences from the average in regions where CO 2 is flowing. Several fluid properties used in the VESA simulations are chosen according to the Sleipner benchmark dataset (Singh et al., 2010), see Table 1.
The small capillary pressure in the sand is neglected, i.e. the simulations are performed in sharp-interface mode, but residual trapping is included. In the VESA model, compressibility is defined by the storativity, S, which is related to fluid compressibility, c f , and compressibility of the rock, c s . The latter assumes negligible grain compressibility and uniaxial deformation, thus Here, ρ is the density of the fluid, p is pressure, ϕ is porosity and H is the so-called uniaxial, or P-wave, modulus. The total compressibility is defined as c T = S/ϕ = c f + c p .

Coupled geomechanical-flow models
Poroelastic regional model. A regional Biot poroelasticity model coupled to single-phase fluid flow (Wangen et al., 2016) was applied to the storage unit and its over-and underburden to investigate uplift caused by the injection. The simulator has been benchmarked against the analytical Theis solution for transient fluid pressure and 1D and 2D analytical solutions for deformations caused by pore fluid pressure variations. It solves for deformation and stress from a near-well scale to the basin scale. The stress state is divided into two parts: (1) the initial stress state and (2) a poroelastic stress change that is due to a change in the fluid pressure. The equilibrium equations with respect to poroelastic stress changes are coupled to an equation for pressure buildup and fluid flow. The poroelastic fluid pressure equation has the source term expressed with stress, instead of the volumetric strain (fixed stress split).
Poroviscoelastoplastic local model. Local deformations in the reservoir and caprock may be significant also where averaged large-scale regional deformations are small. Laboratory experiments on sandy rocks and shales show that the geomechanical response of such rocks is nonlinear at small pressures and displays time-and stress-dependent behavior (Chang and Zoback, 2009;Tsai et al., 2008;van Noort and Yarushina, 2016). The failure limit of these rocks obtained in laboratory tests is fairly low (Mohamadi and Wan, 2016;Skurtveit et al., 2013). Therefore, the rheological behavior can be described as viscoelastoplastic with porosity and pressure sensitive elastic and viscous moduli. CO 2 injection may be affected by local creep deformation allowing vertical fluid flow through an a priori impermeable formation (Räss et al., 2014;, see Fig. 2.
Such flow is affected by viscous rock deformation, which introduces characteristic length and time scales to the fluid flow via the so-called compaction length δ, and viscous compaction time, Here, k is the permeability, η is the bulk viscosity, g is the constant of gravity and v is the difference between the rock and CO 2 densities, i.e. the fluid is assumed to be CO 2 . The compaction length and compaction time were used to render the problem non-dimensional. The model for viscoelastoplastic deformations of porous rocks extends Biot poroelasticity to large irreversible deformations in such a way that the model still remains relatively simple and yet consistent with experimental observations on rock mechanical behavior . The governing equations, as given in the  (Zweigel et al., 2004) with c f = 0.44 GPa −1 (Freeze and Cherry, 1979) CO 2 viscosity μ n μ n (T) (Vesovic et al., 1990) Flow localization. A buoyant pocket of higher porosity filled with CO 2 is overlain by clay, which represents the caprock or a layer within the aquifer. Buoyancy forces acting on the porous pocket lead to its upward propagation.
references above, describe how porosity, and therefore permeability, change as a result of differences in total and fluid pressures.

Reservoir and caprock
Solutions to the model problems defined above are pending on the geological setting, injection rates, initial conditions and boundary conditions. These choices are motivated below, including a description of uncertainties.
Because our previous simulations (Elenius et al., 2016) based on the storage atlas dataset (Halland et al., 2014) revealed a very large range in CO 2 storage capacity depending on the uncertain aquifer extent, considerable efforts were made in the present study to reduce the uncertainty in the horizontal and vertical extent of the Skade Formation. These results are presented below and in more detail in Baig (2016). In addition, because the seismostratigraphic mapping performed here strongly indicates that Skade is connected to Utsira via the proximal sands, these are included in the mapping and are grouped into a single unit referred to as Utsira. However, the lateral extent of the mapping was limited to the extent of Skade. (Later, we merged with a dataset (Kirby et al., 2001) from the British Geological Survey (BGS) as described below.) To complement the seismostratigraphic mapping, we also discuss and make conclusions based on the structural setting.

Seismostratigraphic mapping
Seismic mapping of the Skade Formation and its overburden succession is developed based on the interpretation of 2D seismic reflection and well data, including well logs and well-completion information, see Fig. 1.
Seismic interpretation was carried out in the IHS TM Kingdom software (www.kingdom.ihs.com) and depth conversion together with geological reservoir modeling were completed in the Petrel Software (www.software.slb.com). More than 250 wells have penetrated Skade in the Viking Graben area. Data from these wells is available in the NPD database. However, the database for Skade is inconsistent. In some wells, the top of Skade is a lithological boundary, while in others it marks the sequence stratigraphic boundary. This gives rise to a correlation conflict between the Utsira and Skade formations. As one example, the lower part of Utsira Formation in well 16/1-1 correlates with the Skade Formation in well 24/12-1, cf. Fig. 6 in Rundberg and Eidvin (2005). In order to avoid this correlation conflict, we used data from wells that have a better age constraint from biostratigraphy (Eidvin and Rundberg, 2007;Eidvin et al., 2013). First, the top and base of Skade were taken from these wells, to start seismic interpretation in blocks 16/1, 16/2, 24/12 and 25/10 to the south, and in blocks 25/1 and 25/2 to the north. Then, the top and base were traced outside of this region. The formation tops from NPD factpages were also cross checked with the wire line logs using gamma-ray transit time and gamma-ray deepresistivity log templates. The formation tops were tied to seismic data and updated when necessary according to seismic and wireline log correlation. The Skade Formation is here divided into two zones in the vertical direction, based on log-to-seismic correlation, see Fig. 3. Top Skade 1 is recognized as a continuous, high amplitude peak reflection. Wireline logs show a clean sand with minor shales. The maximum thickness exceeds 300 m at the depocentre between 59°N and 60°N. The  Fig. 1a shows the location of these sections. lower zone, Skade 2, has a more chaotic seismic reflection. Wireline logs show mainly silty shales with some minor sand units. It is more than 300 m thick in block 25/10. The overburden to Skade includes a mid-Miocene mud unit, the Utsira sands and Quaternary prograding units (Fig. 1b). The mid-Miocene mud unit separates the Skade and Utsira formations within part of the Norwegian sector, as shown in Figs. 1b and 3 a. The thickness of the mud unit exceeds 300 m towards the northeast and southwest, whereas in the center, it is mostly less than 100 m thick. The Utsira Formation is overlaid by thick Quaternary prograding units (> 500 m) within the Norwegian sector but the Quaternary succession is less than 200 m thick towards the west in the UK sector.

Depth conversion and uncertainty
A velocity model was built in the Petrel software from interval velocities taken from checkshot surveys in approximately 80 wells. We used averaged interval velocity from the seabed to Base Skade 2 in our velocity model. The calculation of interval velocities was based on checkshot data. Interval velocities at shallow depth between the seabed and top Utsira are somewhat uncertain due to lack of data but for the lower sections, the interval velocities are better constrained. These interval velocities were then used to build the velocity model. Interval velocities were estimated as shown in Fig. 4. The interval velocity for water was set to 1480 m/s. The resulting top surfaces of Utsira and  (Kirby et al., 2001). In the BGS data, there was no mapping of proximal sands in the UK sector.
After converting the seismic profiles to depth, we compared the top of Skade 1 with the detailed mapping presented in Eidvin and Rundberg (2007) and Eidvin et al. (2013). The difference in depth was only ± 20 m, which is considerably less than in the storage atlas, where deviations between surface horizons and well data were an order of magnitude larger (Halland et al., 2014). Variability in the overburden thickness and relative composition (mud/sand) may affect accuracy in the west and east. In the seismic time domain, the mid Miocene mud unit is observed to onlap top Skade 1. However, there is a small gap between the corresponding depth-converted surfaces, which is probably caused by varying velocity. To address this inconsistency in the depth conversion, we constructed a depth-converted model by defining the mudstone wedge as the difference between the bottom Utsira and top Skade 1 horizons, which is depicted in gray in Fig. 6.
The thickness and extent of the mudstone wedge may influence the extent to which Utsira is used as a secondary storage unit. There are, however, some uncertainties to consider. First, uncertainty arises due to the interpolation of interpreted 2D seismic lines. Second, the exact pinch out of the mud unit may be further west as the mud layer might be thinner and below seismic resolution. Third, there are few wells for confirmation of pinch out, especially to the north. Fourth, the varying accuracy in depth conversion gives rise to a resulting uncertainty in mudstone estimated dip angle and western extent. The largest difference in western extent between the seismic and depth-converted maps is approximately 25 km in the north. The minimum western extent of the mud unit is estimated to be approximately the extent as interpreted in the seismic time domain, see the right sub-figure of Fig. 6.

Structural setting
The Lower Miocene Skade Formation has been interpreted as sandrich gravity deposits, sourced from the Shetland Platform in the West (Eidvin et al., 2013). Further basinwards, in the east, fine-grained deposition from suspension fall-out dominated. The Skade sand is most likely unconsolidated given the shallow burial depths in the order of 1 km, and paleo-temperature less than 70°C (Bjørlykke and Egeberg, 1993). Although the sands in Skade 1 mostly appear thickly bedded and amalgamated there may be mud layers present that are below log-scale resolution (< 0.5 m). Turbidites deposited in a deep shelf setting, as interpreted in a depositional environment for Skade by Eidvin et al. (2013), would expectedly comprise aerially extensive sand sheets with draping layers of mud. Compartmentalization in turbidite reservoirs The purple stars mark the location of the injection wells used in simulations and the black circle is the location of well 16/4-1 discussed below. The mid-Miocene mud unit is limited to the eastern part of Skade as shown by a black contour. The outline of Skade is also shown in black.

Fig. 6.
The three figures to the left display the thickness in meters of the mid-Miocene mud unit with different scales. The thin purple contours show the eastern limit of conditions where CO 2 exists as a gas in Skade 1 (left contour) and Utsira (right contour). The right figure shows the extent of the mid-Miocene mud unit derived in two ways. First, the region in gray depicts the depth-converted boundary derived from the distance between the bottom Utsira and top Skade horizon surfaces. Secondly, the red contour depicts the time-domain boundary interpreted from seismic data, which is considered to be the minimum extent.
can cause difficulties in hydrocarbon production and is a recognized issue (e.g. Fugitt et al., 2000). Lateral compartmentalization is dependent on mud layer geometries, their thickness and frequency, and the degree of consolidation. If significant compartmentalization exists, it can cause local increases in pressure that limit the ability to practically achieve regional CO 2 storage capacity.
The geometry of mud layers depends on the nature of the depositional system. With the depo-center located at approximately 60°N (Halland et al., 2014;Eidvin and Rundberg, 2007) and given the semicircular shape of the sandy deposits, the frequency and individual thickness of flow baffles would be expected to increase with distance from the delta-front, towards the east and south. In turbidite complexes, the relative bed thickness, the degree of amalgamation and the channelization relate to the location within the depositional environment and the nature of the transport processes (Bouma, 2000;Carlson and Grotzinger, 2001). With this in mind, erosive turbidite channels would effectively reduce the likelihood of compartmentalization. In more distal settings, however, mud-layers draping sandy lobes may be extensive (10s of kilometers). Skade 2, representing an early phase of deposition, displays thinner-bedded sands and appears to represent more distal deposits compared with Skade 1. Additionally, experience from Sleipner shows that mudstones are acting as flow baffles but do not noticeably block pressure migration. Therefore, the risk for reservoir compartmentalization would therefore be greater in Skade 2.
The potential for intra-formational mud layers to act as efficient fluid and pressure seals also depends on the degree of consolidation of mud to shale. In the Skade Formation, muds are likely still highly porous and unconsolidated, which would allow for pressure dissipation and potentially also sediment-fluidization with formation of chimneys (Räss et al., 2014;. In conclusion, the risk of severe pressure buildup due to intra-formational mud layers in at least Skade 1 is considered to be low. The mid-Miocene mud unit (basal part of the Nordland Group) forms a seal above the eastern parts of Skade. Haugvaldstad (2014) showed that sands in the lower part of the Utsira Formation were likely extruded from the underlying Skade Formation. This is documented only in the far south-west of the Skade deposits, and sand extrusions are not proven to penetrate the mid-Miocene mud unit in the east. If the mid-Miocene mud unit is not intact, or if CO 2 migrates to the western extent of this unit, some of the injected CO 2 may escape to Utsira, which can be considered as a secondary storage unit. This may or may not lead to CO 2 forming a gas, depending on the pressure and temperature at the location where a CO 2 plume forms, cf. Figs. 6 and 7 . CO 2 as a gas occupies a larger volume and this might lead to higher overpressures. On the other hand, there are large volumes for this pressure to dissipate, and if mainly brine escapes to Utsira, this would be beneficial as it would reduce pressure buildup in Skade around the injection wells.
In addition to evaluating the structural setting of the storage units and the mid-Miocene mud unit, the overburden is of interest due to the risk of leakage. Gas pockets (unrelated to CO 2 storage) have been observed at approximately 500 m depth in the overburden of Utsira, but there are no gas pockets present below a specific fluid-expelling seafloor fracture called the Hugin fracture (Pedersen et al., 2013). This could indicate leakage of gas through a complicated network of channels and fractures. Fractures may have formed during glacial loading and unloading with pathways starting at the interface between clay and sand-dominated fill in shallow buried channels (Landschulze et al., 2014). The offset between the so-called Hugin fracture observed at the seafloor and a buried glacifluvial channel at 50 m depth below the sediment-water interface is consistent with results of geomechanical simulations employing glacial loading and unloading. Other fractures similar to the Hugin fracture may have formed preferentially in vertical Quaternary deposits as a result of glacial loading and unloading in the presence of sand channels. It is less likely that these structures would occur in deeper formations such as the mid-Miocene mud unit between Skade and Utsira, or in the deeper parts of the overburden of Utsira. However, there are uncertainties with respect to the western extent of the system considered here, where the caprock is more shallow and the western boundary is not mapped.

Reservoir model
In this section, we first describe the geological facies model, which consists of one geostatistical realization tied to wells. Then, rock properties applied in simulations are given, including an estimate of parameter uncertainty. The uncertainty in the most important parameters, permeability and compressibility, is explicitly accounted for in a sensitivity matrix. Then, the averaging used in regional models is explained, followed by an evaluation of the pressure that can be tolerated by the caprock without inducing mechanical failure, the boundary conditions and the injection strategy.

Geological model
The new stratigraphic mapping performed here was merged with the existing BGS Utsira dataset (Kirby et al., 2001). A geological model was developed for the area mapped here, assuming uniform conditions in the remainder of the system. The model also includes the mid-Miocene mud unit that separates Utsira from Skade in part of the domain. The mid-Miocene sand unit in the west (Fig. 1b) was merged with Utsira, because their properties are similar and they cannot be separated in the seismic data. In addition, the large-scale geomechanical simulations include the mud/shale units above and below this model domain. Cross-sections of Utsira and Skade are illustrated in Fig. 7 and the top surfaces and thicknesses are shown in Fig. 8.
The facies model was created with the aim of achieving overall netto-gross ratios of sand distributed in the latitudinal and longitudinal directions (as opposed to being used directly as a 3D model), see Fig. 8. Discrete lithology logs in wells were generated based on volume of clay (calculated from porosity and gamma ray logs) and divided in sand, fine sand and clay. The overall fraction of each litho-facies is presented in Table 2.
The lithology logs were upscaled to model resolution, which is 1000 × 1000 m in northing/easting and variable in the vertical direction with a constant number of grid cell layers (10 in Utsira, 20 in Skade 1 and 10 in Skade 2). The distribution of lithologies was generated using sequential indicator simulation in Petrel, using the overall lithology fractions (%) from well data ( Table 2) as input and conditioned to the upscaled logs. A spherical semivariogram model was used with ranges 57,500 m and 30 m in the lateral and vertical directions respectively, as obtained from the well data. The indicator simulation honors the local upscaled lithologies in wells, in as much as the grid resolution allows. Uncertainty in the distribution between sand and clay is larger in the UK sector due to lack of well data. The mudstone wedge was assigned a uniform mudstone lithology.
Next, we review the modeling choices for porosity, permeability, and rock-mechanical properties. A summary is given in Table 3. Very little is known about Skade rock-mechanical properties. However, because Skade is unconsolidated similar to the nearby Utsira, several measured Utsira parameters were taken to representative of Skade.

Porosity and permeability
Although the porosity may vary with pressure, the variability measured in the Utsira sand by Springer et al. (2002) was negligible in the context of other uncertainties. The estimated range was 0.36-0.37 at 100 bars and 0.38-0.39 at 70 bar. We therefore apply the value at 100 bars to the entire aquifer system. We also apply the sand porosity for clay within Utsira and Skade. The clay above Utsira has porosity 0.35 as given by Pillitteri et al. (2003). We assigned the same porosity for the mid-Miocene mud unit and the shale below Skade 2. Springer et al. (2002) also measured the permeability of Utsira sand in several core sections and conclude that the gas and Klinkenberg corrected gas permeability were within 1.5-2.5 D. The liquid permeability was lower at 1-1.5 D. The reduction was mainly due to fines migration but it was not clear if that originated from the sand itself or from drilling mud. We use a uniform distribution between 0.5 and 3 D to account for the measured range with an additional uncertainty of 0.5 D in either direction. Clay permeability was set to 1 × 10 −9 D (Springer and Lindgren, 2006).
When considering local effects of flow through the caprock and clay lenses in our viscoelastoplastic model, we account for a strong coupling between geomechanical deformation and fluid flow. This results in a porosity-dependent permeability (Dong et al., 2010;van Noort and Yarushina, 2016) and we use k = k 0 (ϕ/ϕ 0 ) n with initial values represented by subscript 0, as given in Table 3. The exponent n can vary between 3 and 20 in shales and clays (Dong et al., 2010). The exact values for Mid-Miocene mud and Utsira caprock are unknown, so we use n = 3 as it provides the lower bound. Heterogeneity of the rock is represented by vertical variation in rock bulk viscosity in the range η = 10 17 − 10 19 cP.

Rock-mechanical properties
The geomechanical response to large-scale CO 2 injection has a significant impact on the amount of CO 2 that can be stored safely, and the related parameters are therefore described here in detail. First, note that poroelasticity allows for large number of different moduli that are dependent on each other. For example, the P-wave modulus H, also known as the longitudinal or uniaxial modulus, is related to the Poisson ratio ν, Young's modulus E and the drained bulk modulus K via, Pillitteri et al. (2003) measured E and ν in two samples from the Nordland shale above the Utsira. The measured values for each sample were E = 0.19 and 0.29 GPa, and ν = 0.39 and 0.18. The average values are used here; E = 0.24 GPa, ν = 0.29. We then compute the clay compressibility 1/H = 3.2 GPa −1 based on Eq. (8). Because pressure buildup is less sensitive to the overburden and underburden properties compared to the aquifer properties, we do not account for uncertainty in the former. On the other hand, as the underburden will be likely more consolidated than the overburden, we represent it as a stiffer clay (shale) by increasing Young's modulus by a factor of ten.
In the sand, there are measurements of compressibility, but not of Young's modulus or Poisson's ratio. Springer et al. (2002) measured 1/ (Kϕ) = 6.0 × 10 −5 psi −1 = 8.7 GPa −1 in Utsira sand, at hydrostatic (i.e. isotropic) confining pressure s = 70 bars and porosity ϕ = 0.39. This gives a bulk compressibility 1/K = 3.4 GPa −1 . Measurements of compressibility were also performed at other confining pressures as shown in Fig. 9a (red circles), and are seen to plateau at large confining , net-to-gross sand ratio (N/G) and limiting pressure (plim) [bar] for the four horizons, here limited to the region currently mapped. Eir represents the mid-Miocene mud unit here. The lower sub-figures show maximum sustainable pressure increase, plim, which is defined as the difference between the initial hydrostatic pressure and the pressure that the caprock can tolerate [bar], cf. Section 5.5. The shallowest part of Utsira is only 238 m below sea level, corresponding to 8 bar allowable overpressure as the most critical value. pressures. The last measurement is therefore a reasonable approximation to our conditions of 150-200 bar confining stress at injection. Mavko et al. (2009) give a range in compressibility of poorly consolidated sandstone that varies by a factor of two. Because their variability was for one given effective stress, we account for larger uncertainty and vary 1/K a factor of 3 below and above the fitted model as shown in Fig. 9. The Poisson ratio ν is needed to convert between different moduli K, H and E, cf. Section 3.3. There is no measurement of ν in the Utsira or Skade sands. Mavko et al. (2009) give the range 0.30-0.34 in poorly consolidated sand at 30 MPa effective pressure. Accounting for additional uncertainty with respect to differences in pressure regimes, we apply a range between 0.25 and 0.39. The resulting variability in 1/H and E is illustrated in Fig. 9. The viscoelastoplastic simulations additionally employ the rock facies viscosities. As we do not have data on viscosity of clay, shale or sandstone from Utsira and Skade, we use typical parameters as reported in Räss et al. (2014), Räss et al. (2017).

Parameter upscaling for regional models
The 3D geological model was developed to provide the basis for upscaled parameters used in the regional simulations (Sections 3.2 and 3.3). For example, the modeling in 3D gives important insight into the proportion of clay present in the different horizons, however, we presume that the lateral extent of individual clay layers in the 3D model is overestimated since erosional contacts (i.e. sand-sand) are not modeled. As such, the 3D model overpredicts compartmentalization at any given point. Instead, we have used the 3D model to produce depth-averaged values of net-to-gross ratio (per area of aquifer), accounting for different levels of compartmentalization in Utsira, Skade 1 and Skade 2 according to the estimated extent of erosional contact in each horizon, described below.
The upscaled values in Utsira Skade 1 are based on the assumption of no compartmentalization, but very low permeability of clay. The upscaled porosity is therefore taken as the intrinsic porosity multiplied with the net-to-gross ratio, f s (x, y). The same upscaling is used for the horizontal permeability, where the sand permeability k s is multiplied with the same factor f s . The vertical permeability is assumed equal to the horizontal permeability. We note that this value is not used in the VESA simulations due to the assumption of vertical equilibrium, and it has only a small impact in the regional poroelastic simulations because the flow is dominated by the lateral permeability for pressure plumes that spread out with a lateral extent that is much larger than the thickness of the reservoir. Due to the large fraction of clay and risk of compartmentalization in Skade 2, the flow field is likely more complex, and the upscaled permeability is further reduced. Using subscripts s for sand and up for upscaled quantities, The VESA model takes the regionally varying f s as the local input to a heterogeneous model, whereas the poroelastic model applies homogeneous properties within each zone, i.e. f s is calculated from the first two columns of Table 2.

In-situ stress field and failure criteria
Pressure buildup in the storage formation must be controlled to avoid failure of the caprock and leakage of CO 2 . In fact, our previous simulations showed that pressure buildup may limit the storage capacity during large-scale injection (Gasda et al., 2005). In the flow simulations (Section 3.2), a regional pressure criteria is therefore used to  Values below 0.24 GPa at high pressure are not realistic, because it is not likely that the sand is softer than the unconsolidated clay above.
limit the injection rate [Mt/well region/year] and avoid failure. The poroelastic regional model (Section 3.3) applies the same volumetric injection rates, where a constant density 700 kg/m 3 was used in the conversion from mass to volume. As a first step towards establishing the tolerable pressure increase, the initial stress state of the rock, including the faulting regime, is evaluated. The latter defines the relation between the vertical stress s v , maximum horizontal stress s H and minimum horizontal stress s h . Regional databases on faulting regimes based on analysis of earthquake focal mechanisms often focus on the brittle crust (Lindholm et al., 1995(Lindholm et al., , 2000, which below Skade has strike-slip component (Zoback, 1992 Here, z is depth below sea level, h = 117 m (NPD, www.factpages.npd. no) is mean seawater depth in the Skade region, 1022 w = kg/m 3 (Zweigel et al., 2004) is brine density (Utsira pore water, similar composition as sea water), and ρ b = 2100 kg/m 3 is wet bulk density as measured for the Utsira caprock (Pillitteri et al., 2003). This value is a constant value assumed for the entire overburden, which was chosen for convenience in lieu of integrating the density logs.
We then approximate s h based on leakoff tests. There are 387 leakoff pressures (LOP) measured from 207 wells penetrating Skade (www.factpages.npd.no). Among these, we restrict the dataset to the 202 tests that are performed at representative depths shallower than 1500 m, see Fig. 10a. When evaluating s h from leakoff tests, LOP is usually estimated to vary between s h and the formation breakdown pressure, FBP. On the Norwegian Continental Shelf, an upper bound for FBP, and therefore also for LOP, is generally set by Kirsch's equation, which gives LOP ≤2s h − p h (Andrews et al., 2016). This corresponds to a lower limit for s h , for which the regression based on LOP gives with λ = 0.134 bar/m (13.4 MPa/km), see Fig. 10a. (We also tried adding a constant in the fitting procedure to account for sea depth but this gave negligible improvement due to shallow sea water depth.) As seen in the figure, deeper parts of the formation can withstand higher overpressures. Our estimate is slightly lower than a typical range of 15-20 MPa/km discussed by Thibeau and Mucha (Thibeau and Mucha, 2011), probably reflecting the conservative choices we made, in addition to addressing different regions. For well 16/-1, the pressure limit has been estimated to be 13.5-19 MPa based on well log measurements and compaction of reconstituted drill cuttings sample (Nooraiepour et al., 2017), giving a good agreement with the estimates used in our model.
Our regression can also be compared with that of Breckels and van Eekelen (Breckels and van Eekelen, 1982), which was based on a dataset collected from a larger extent of the North Sea (Fig. 10b), here presented in units of bars and meters At greater depths, this gives a higher estimate of s h . These authors also recognized the uncertainty in the interpretation of leakoff tests and chose to represent s h by a fit to the lowest leakoff pressures of their examined dataset, at any depth. Recent evaluation of different approaches for fracture pressure determination (Bohloli et al., 2017) suggest that fracture pressure might deviate from predictions and should be determined with greater accuracy during operation. There are two common failure scenarios related to increasing fluid pressure. The most critical shear failure is for cohesionless material, e.g. faults, whereas tensile failure will be most critical in intact material. Tensile failure occurs when the pore pressure exceeds the sum of the least horizontal stress and the tensile strength. The tensile strength can be quite low. Therefore as a conservative choice, we neglect the tensile strength and arrive at the following tensile limit for the pore pressure, The areal distribution over the Skade and Utsira formations is shown in the lower sub figures of Fig. 8. Occasionally, estimates of the tensile pressure limit have instead been performed using a fraction of the lithostatic stress, e.g. Wangen, 2002) (Fig. 10b). Andrews et al. (2016) had access to a large database of XLOT measurements on the Norwegian Continental Shelf. XLOT data measure the fracture closure pressure compared with a LOT which measures the fracture opening pressure. Fracture opening requires overcoming both s h and tensile strength, while fracture closure requires only s h , and therefore XLOT tests are preferred over LOT when available. Analysis of XLOT data in Andrews et al. (2016) found that the minimum stress was greater than 84% of the maximum stress in 50% of cases, which would allow for higher overpressures than we have used. Since this was not limited to our region, we have not used it. In light of this, and because we neglected tensile strength and used the lower limit when interpreting s h from LOP, our estimate for the tensile pressure limit can be viewed as conservative. Detailed analysis of faults penetrating Skade and into the overlaying caprock has not been performed. The clay-rich caprock is expected to be uncemented with a low degree of consolidation and its behavior is expected to be weak and ductile (Pillitteri et al., 2003). In addition, no faults are observed in the seismic data. Therefore, cohesionless faults are not considered a high risk for this system.
For the relatively small extent of the caprock that comes in contact with CO 2 , other potential pathways are seepage of CO 2 into the caprock if the pressure exceeds the entry pressure, and formation of localized flow pathways through creep effects. Given the low permeability of the caprock and significant vertical distance to the seafloor, seepage into the caprock would occur at exceedingly slow rates, taking hundreds or thousands of years to reach the seafloor. Therefore, we conclude that seepage would pose no additional risk for CO 2 leakage to the seafloor. We study creep effects in Section 7.2.

Boundary conditions and injection strategy
Because Skade is connected to Utsira, the dataset developed here was merged with the BGS dataset for Utsira as described above. This accounts for the entire aquifer system, which is therefore considered to be closed in the flow simulations. An exception is the western extent of the Utsira sand, which is here treated as closed although it may outcrop to the seafloor. As shown by the simulation results, this does not significantly impact storage capacity in the majority of cases considered. In coupled geomechanical-flow simulations, the system up to the watersediment surface is included, and the surface is free and open, while other boundaries are closed.
Potential leakage through abandoned wells or other high-permeability pathways such as existing sand injectites are not considered in large-scale simulations. Nonetheless, leakage due to adverse pressure buildup is considered when setting constraints on the injection rate.
The targeted simulations of viscoelastoplastic effects are performed on small domains extending 30-90 compaction lengths, δ (Eq. (6)), which is on the order of meters to tens of meters depending on formation and fluid properties as explained below. Along all outer boundaries of the domain, the rock matrix is stagnant and the normal component of the fluid flux is set to zero. CO 2 injection is simulated for 50 years from three well regions, see Fig. 5 and Table 4 . In each well region, the injection interval is from the top of Skade 1 to the bottom of Skade 2 in the injector cell. Given the large cell sizes, this can be considered a region with several wells. The locations were chosen based on a few trial simulations to reduce pressure buildup while increasing the fraction of CO 2 remaining in the primary storage target, Skade. Favorable factors for reduced pressure are large aquifer thickness and distance from closed boundaries. As explained before, the constant flow rate of each case was optimized to maximize storage capacity while keeping the caprock intact. The optimized rates produce pressures between 96 and 100% of the limiting pressure. Simulations with VESA continue for 450 years after injection stops and the pressure constraint applies at all times and all points along the top surface of Skade and Utsira. Table 5 shows the estimated rock and pore volumes. We used net-togross ratio and porosity from Tables 2 and 3 respectively. Despite the smaller areal extent compared to Utsira, the Skade is thicker and thus has a larger rock volume. However, the total pore volume of Skade is smaller due to the low net-to-gross ratio of Skade 2 Based on the sensitivity matrix, the total compressibility of cases considered is c T = 2.1, 12 and 21 GPa −1 , cf. Tables 1 and 3 . The least tolerable pressure increase in our domain is 8.1 bar. This gives (volumetric) efficiency factors of 0.17, 0.94 and 1.7%, which is a slightly larger range than 0.3-1.2% suggested as a "rule-of-thumb" by Goodman et al. (2013) and Bachu (2015) and smaller than 4% used in the NPD storage atlas (Halland et al., 2014). The effective storage capacity with our efficiency factors is approximately 1, 5 respectively 10 Gt, see Table 5.

Effective storage capacity
Storage capacity was also evaluated with the Theis solution subject to the tolerable pressure increase 500 m from the wells at 50 years. We used the thickness of the well regions (well cells), the upscaled permeability and the total bulk compressibility S = c T ϕ. For our 9 combinations of low, medium and high permeability versus compressibility, the Theis solution gives capacity (1.2, 1.5, 1.7) Gt for k = 0.5 D, (3.7, 4.5, 4.9) Gt for k = 1.75 D, and (6.0, 7.2, 7.8) Gt for k = 3.0 D.
As discussed before, some CO 2 would dissolve into the formation brine. An upper estimate for the amount of dissolved CO 2 after 500 years by Eq. (4) is between 0.2 Gt (lowest permeability) and 2.5 Gt (highest permeability). Here, we used approximate CO 2 plume areas 1500 km 2 respectively 3000 km 2 , based on simulation results. Parameters not given above were estimated from Elenius and Johannsen (2012). For large permeability the upper estimate of dissolution is therefore large compared to the storage capacity. For small permeability, dissolution is not important during the first 500 years but would continue to act over longer times. Dissolution was not included in the simulations below.

Large-scale behavior
Here, we first describe VESA results of the system response to the optimized injection rate of 17 Mt/well region/year (2.6 Gt in total) using the medium permeability k = 1.75 D and the lowest compressibility 1/H = 0.6 GPa −1 . Although the entire domain as shown in Fig. 5 was simulated, a zoomed region close to injection is shown here and in subsequent figures. Figs. 11 and 12 show the CO 2 migration and pressure buildup. The CO 2 plume reaches Utsira. The injection rate was limited by the pressure buildup in Utsira at the end of simulation, i.e. 500 years. At this time, all of Utsira is exposed to the pressure plume, which shows the importance of including the larger aquifer system in simulations.

Viscoelastoplastic effects
We have evaluated the effect of rock viscosity on CO 2 seepage through the caprock and through clay lenses present in the sandstone formation. The results are first presented in non-dimensional units and then interpreted in terms of the propagation in caprock and lenses. The non-dimensional time and length shown here can easily be converted to dimensional units for rocks of different properties by multiplication with the scales given by Eqs. (6) and (7).
Due to strong coupling between the flow and viscous deformation, fluid overpressures develop locally in a clay barrier. This results in an effective weakening of the solid matrix and generation of shear stresses locally. This further leads to the opening of the pore space in regions Table 4 CO 2 injection wells. The depth (z 1 ) is at the top of Skade 1 at the well and the local thickness is given for Skade 1 (h 1 ) and Skade 2 (h 2 ). We used geothermal gradient 0.0373°C/m with reference temperature 12.5°C at 400 m depth as reported for Utsira (Lindeberg et al., 2009 Elenius et al. International Journal of Greenhouse Gas Control 79 (2018) 252-271 where the fluid pressure is the highest and facilitates formation of localized vertical channels. Increased fluid pressure at the top of the channel leads to further opening of a pore space ahead of the channel, see Fig. 13. This is accompanied with the compacting behavior at the bottom of the channel. Thus, fluid-filled channels can propagate upward in a self-sustainable solitary-wave-like manner. Presence of heterogeneities in the rock disturbs channel formation slightly causing their diffusion and nonhomogeneous distribution of porosity, fluid pressure and stresses within and around the channel, see Fig. 14.
However, the momentum that these channels acquired during propagation is enough for their upward propagation as long as the permeability does not vary significantly. If such channels would meet another layer of yet more impermeable but still viscous rock, they would be halted temporarily but then resume their propagation on another time scale.
Channel formation exhibits self-organizing features in that the number of channels and their spatial distribution are controlled by the compaction length, δ, which depends on the permeability, fluid viscosity and compaction viscosity of the rock. Assuming typical values Table 5 Volumes and capacity. For Skade 1 and Skade 2, values are given as part underlying the mid-Miocene mud unit + part outside of this unit. Only the part below the mud unit can be considered as the primary storage unit of Skade, since CO 2 in the remaining area will quickly reach the secondary storage unit, Utsira. We used ρ n = 700 kg/m 3 for Skade below the mid-Miocene mud unit and otherwise ρ n = 100 kg/m 3 , which is a representative CO 2 density (gas) where CO 2 reaches west of the mud unit (the lowest density is 55 kg/m 3 ).    Fig. 11. μ = 0.06 cP for CO 2 at near-well conditions, η = 10 18 cP and clay permeability in the range 10 −9 -10 −7 D (Springer and Lindgren, 2006;van Noort and Yarushina, 2016), the compaction length will vary between 0.1 and 1 m and the velocity of channel propagation varies between 0.01 and 1 m/year. The width of one channel is up to a few compaction lengths. The time for propagation through clay lenses of  Evolving permeability increase (k/k 0 ) and associated porosity in a domain with initial heterogeneity consisting of layers and a square inclusion with higher initial viscosity. The domain size is here 30 × 80 compaction lengths. Comparison with Fig. 13 shows that heterogeneities only causes small disturbances in rock properties in the immediate vicinity of these heterogeneities. thickness 1 m is therefore 1-100 years which in most cases is quick enough to reduce the risks related to compartmentalization. The time for propagation through 100 m caprock is 100-10,000 years, and this propagation can in most cases be neglected compared to propagation around the western extent of the mid-Miocene mud unit.

Practical storage capacity
Based on VESA flow simulations (Section 3.2) and regional poroelastic simulations (Section 3.3), the sensitivity of the storage capacity is investigated with respect to the permeability and the uniaxial bulk compressibility of the sand, in total 9 cases. Upscaled values were used in simulations, according to Section 5. The constant injection rate of flow simulations was optimized in all cases until the storage capacity in terms of pressure buildup was reached. Due to the nonlinear nature of the problem including impacts of CO 2 migration on pressure buildup especially at late times, this required up to four simulations per case, and these had to be run for the entire period of investigation, which is 500 years. This was possible to achieve because we used the VESA model, requiring only between 20 min and 2 h execution time for a single simulation on one processor. The regional poroelastic simulations then applied the same rate to simulate the injection period of 50 years. Fig. 15a and b shows the CO 2 plume in Skade at the end of injection and 450 years later. We show the results of the 9 simulations with optimized injection rates, in which the end-point and the average parameter values of compressibility and permeability were applied. For these and following figures, the permeability and compressibility are indicated for each row and column, respectively. In addition, the injection rate Q [Mt/well region/year] is printed out in each sub-figure.
In all but one case, the advancing tip of the CO 2 plume reaches west of the mud unit and migrates to the top of Utsira during the 500 year simulation period, see Fig. 15b. When the CO 2 reaches Utsira, it spreads towards the western extent of the model domain (Fig. 15c). Fig. 16a depicts the overpressure in Skade at the end of injection.
The pressure is concentrated around the injection wells mainly in simulations with low permeability and high compressibility. Fig. 16b and c shows the overpressure in Skade at end of injection and in Utsira at end of simulation. The pressure is here shown as a percentage of the allowable overpressure. Note that the shallow western extent of Utsira has a very low tolerance to overpressure, down to 8 bars. The distribution of CO 2 mass is shown in red in Fig. 17. The cases with permeability of 500 mD have negligible CO 2 mass in Utsira at the end of simulation. The largest fraction of CO 2 that migrates to Utsira is 18%, but migration has not ceased after 500 years. Approximately 20% of the CO 2 is residually trapped in Skade. If dissolution or sub-seismic caprock roughness had been included, the amount migrating west of the mud unit would be lower. Fig. 17 further displays the most critical overpressure in Skade and Utsira during the 500 year simulation period. With the lowest permeability, 0.5 D, the local pressure around injectors limits the practical storage capacity. When also the compressibility is at the lowest value, the capacity is 9 Mt/well region/year, i.e. a total of 1.4 Gt. With higher permeability, the storage capacity is instead limited by pressure buildup in Utsira. The greater permeability gives rise to faster CO 2 migration to Utsira where it becomes a gas. The expansion leads to higher pressures. In addition, the buoyancy pressure of the thick CO 2 gas column in Utsira causes additional buoyancy pressure on the caprock. The maximum injection rate, 37 Mt/well region/year (5.6 Gt) is therefore not at the highest permeability, but at the intermediate permeability. Here, limited amounts of CO 2 has reached Utsira, and in addition, the local injectivity around the wells is good. If the simulation period was extended, more CO 2 would reach Utsira also for the intermediate permeability. However, at late times, also other processes not included here would become important and reduce the overpressure. In particular, a portion of the CO 2 would dissolve, pressure would dissipate through the over-and underburden, and pressure would dissipate to sands that are not mapped to the west of the simulation domain. In summary, the total storage capacity over the 50 year injection period varies between 1 and 6 Gt with our assumptions (3 well regions, etc.). for the nine simulation cases (VESA model). In red, we show the distribution of CO 2 mass as percent of total injected CO 2 : solid line for mobile CO 2 in Skade, dash-dotted line for residually trapped CO 2 in Skade and dashed line for CO 2 in Utsira (all in Utsira is mobile). The total capacity is labeled in each sub-figure. In blue, we show the most critical overpressure as percentage of the allowable overpressure: solid line for Skade and dashed line for Utsira.
The regional coupled flow-poroelastic model was used to examine if large-scale uplift of the seafloor would likely pose an additional constraint on storage capacity. First, the overpressures produced with this model (Fig. 18a) are compared to those of the VESA flow simulations (Fig. 16a). Overall, the overpressures show the same trends, but with differences due to neglecting CO 2 viscosity and compressibility in the poroelastic model (Section 3.3), both of which increase pressure in this single-phase model. The seafloor uplift obtained with the coupled model is shown in Fig. 18b. It is at most 10 m, which is significantly greater than uplift observed in other CCS projects. However, those injection rates are significantly smaller than those used in this study, so larger deformation is to be expected. We may also compare a 10-m uplift with deformation associated with offshore petroleum production, such as the 9-m subsidence observed at the Ekofisk field (Jones et al., 2014), although the latter is for carbonate rocks that do not compact poroelastically. The sea depth is over 100 m in the Skade region and we therefore do not anticipate that uplift would cause additional constraints on storage capacity.

Discussion and conclusions
Uncertainty in input data is often very large for untested aquifers. Here, uncertainties related to the large North Sea Skade/Utsira storage complex were reduced to a large degree based on interpretations of existing data by experts in various disciplines. The effective storage capacity was then estimated at 1, 5 or 10 Gt, depending on the compressibility, using the updated pore volumes and pressure limits. After this, we evaluated the practical storage capacity of a scenario with 50 years of injection in three injection zones in Skade. Regional uplift of the seafloor is limited to 10 m and is not expected to pose constraints on the practical storage capacity. The capacity is instead limited by pressure buildup and varies between 1 and 6 Gt depending on the compressibility and permeability of the aquifer system. This shows that CO 2 storage in the Skade Formation can contribute significantly to the EU target of reduced emissions (12 Gt until 2050). If XLOT measurements of stresses in the Utsira/Skade region were made available, arguments could probably be made to allow larger overpressures, which would enhance the effective as well as the practical storage capacity. Fig. 19 shows the practical (simulated) storage capacity together with the effective capacity. In addition, we compare with the capacity obtained by the Theis solution. Only a few cases have good "prediction" of the storage capacity from the effective capacity or the Theis solution.

Comparison of capacity estimation methods
The practical storage capacity with different compressibility and permeability can be classified as below:

Lowest permeability (500 mD):
The practical capacity is limited by local pressure buildup around wells and increases with compressibility. The wells are not interfering and the practical capacity could be increased with further wells in Skade (and Utsira). The Theis solution underestimates the capacity, probably because it assumes a constant aquifer thickness, while the actual thickness increases to the west of the injectors.

Medium and highest permeability (1750-3000 mD):
The practical capacity is limited by pressure buildup in Utsira, and it may be preferable from an economic standpoint to reduce the number of wells in Skade. The capacity increases with enhanced compressibility. It is larger with the intermediate permeability than with the highest permeability, because of reduced CO 2 migration to Utsira. Despite the fact that far-field pressure buildup limits the capacity, the effective capacity does not provide good estimates. For low compressibility, the effective capacity is lower than the simulated capacity. This is probably due to greater expansion of the aquifer near the wells where overpressures are large. This expansion leads to reduced pressures. With higher compressibility, the effective capacity overestimates the capacity. Further investigations using these permeabilities would require mapping of the western extent of Utsira and analysis of the impacts of dissolution and pressure dissipation through the over-and underburden. In addition, the quality of the shallow overburden to the west should be evaluated, including mapping possible sand channels and connectivity up to the Quaternary sediments.
From this we conclude that analytical estimates give the correct order of magnitude of the storage capacity simply because compressibility is the dominating factor for closed systems (even very large basin such as the Skade/Utsira). More detailed estimates of the true capacity needed for extended development of large storage resources (through multi-site injection) must be made via a combined approach of smart characterization and stochastic-driven dynamic simulation.
Brine production was not included here. By comparison with oil and gas terminology, we can think of our approach in terms of "primary recovery" (Thibeau et al., 2014;Bachu, 2015). Already this primary recovery is very large in some cases, but would be further enhanced with brine production. Other considerations such as local viscoelastoplastic effects may increase vertical connectivity within Skade but are not expected to cause leakage risk through the caprock.

Implications for sparse data in untested aquifers
The main goal of this study was to improve capacity estimation of untested aquifers where data are sparse. We recommend a workflow starting from optimal use of existing data through to uncertainty analysis by stochastic simulations. Below we list the key guidelines to aid future studies of this type: • Cross-disciplinary collaborations to guide characterization and focus resources on the most critical elements of the basin, particularly with regard to large-scale pressure communication internal and external to the storage aquifer.
• Geological mapping of features relevant to large-scale CO 2 migration and pressure dissipation. These include determination of aquifer boundaries, evaluation of caprock quality and evaluation of risk for compartmentalization (and loss of capacity) due to lowpermeability zones.
• Use of existing data, preferably including XLOT measurements, core logs and seismic surveys, to constrain parameters relevant for pressure and plume propagation.
• Simulations performed using reduced modeling tools, such as vertical equilibrium simulators, that resolve relevant processes but are still sufficiently fast to allow for stochastic analysis of the uncertainty.