Hydrogen storage in a heterogeneous sandstone formation: dimensioning and induced hydraulic effects

Large-scale energy storage in the geological subsurface (e.g. by storing hydrogen gas) may help to mitigate effects of a fluctuating energy production arising from the extensive use of renewable energy sources. The applicability of hydrogen (H2) storage in a porous sandstone formation is investigated by defining a usage scenario and a subsequent numerical simulation of a storage operation at an existing anticlinal structure in the North German Basin. A facies modelling approach is used to obtain 25 heterogeneous and realistic parameter sets. The storage operation consists of the initial filling with nitrogen used as cushion gas, the initial filling with H2, and six withdrawal periods with successive refilling and shut-in periods. It is found that, on average, the storage can sustain a continuous power output of 245 MW for 1 week when using five storage wells, while peak performance can be as high as 363 MW, indicating that the storage is mainly limited by the achievable extraction rates. The median of the maximum pressure perturbation caused by this storage is around 3 bars and can be observed at a distance of 5 km from the wells.

One option to reduce carbon dioxide emissions is to increase the share of electricity produced from renewable sources. Currently, 27.4% of the electricity consumed in Germany in 2014 is produced from such sources, with the overall target of the 'Energiewende' being 80% by the year 2050 (BMWi 2015). Locally, in areas of strong wind-energy production, about 100% of the electricity is provided by renewables. Electricity produced from renewable sources, however, is subject to strong fluctuations due to the local weather conditions. Shortages in power production will thus occur which could last up to several days (Klaus et al. 2010). The idea of an economy based on renewable power production and using hydrogen (H 2 ) as a storage medium had already been introduced in the mid-1970s (Sørensen 1975). In recent years, more research has been conducted in this field: for example, Sørensen et al. (2004) investigated several usage scenarios for H 2 as an energy carrier in the context of Denmark's plan at that time to cover 50% of its total energy supply by renewable sources by 2030. Klaus et al. (2010) discussed scenarios for a 100% renewable power supply by 2050 for Germany that included energy storage utilizing H 2 . In fact, powerto-gas concepts using H 2 are being investigated worldwide with 48 projects in 15 countries being either planned or already realized (Gahleitner 2013). Depending on the characteristics of the shortage periods, large storage capacities and high power supply rates are necessary.
The geological subsurface is already being used for bulk energy storage (e.g. natural gas storage) as it can offer very large storage capacities and long withdrawal periods (e.g. Gregory & Pangborn 1976;Carden & Paterson 1979;Ogden 1999;Evans & West 2008). Experience in the field of subsurface H 2 storage, however, is scarce. Up to now, pure H 2 gas has been stored in salt caverns in Texas, USA, and in Tesside, UK (Crotogino et al. 2010). Even less experience has been gained regarding porous-media storage of H 2 . The storage of H 2 -rich town gas is the only application in which H 2 has been stored in porous geological formations. Albeit much more widespread, the aquifer storage of town gas near Beynes, France is the only field application documented in the literature (Carden & Paterson 1979;Foh et al. 1979), and has been referenced in subsequent publications investigating geochemical effects (e.g. Panfilov 2010). Even though technical issues arising from the use of H 2 instead of natural gas are diversesuch as the corrosion of installations, as well as potential impacts on well and formation integrity (Reitenbach et al. 2015)natural gas storage can be used as an analogue for porous-media H 2 storage as the overall principles are comparable (Carden & Paterson 1979;Foh et al. 1979). As more and more types of use in the subsurfaceranging from groundwater abstraction over geothermal applications to large-scale energy storagecompete for the limited suitable space, analyses of induced effects and impacts should be conducted in addition to an assessment of storage capacity and retrieval rates . For this, an integrated concept for site investigation and monitoring can be used, as explained, for example, in Bauer et al. (2012) for the application of CO 2 storage. In this approach, synthetic but realistically parameterized storage sites are used to simulate the impacts and induced effects, such as, for example, the pressure increase in the subsurface. These are then used to test and improve investigation and monitoring methods. This approach has been applied by Benisch & Bauer (2013) for pressure monitoring, and by Benisch et al. (2014) for assessing geophysical monitoring methods for gas-phase detection in the context of CO 2 geological storage. Pfeiffer et al. (2016) could successfully transfer this concept to geophysical monitoring of hydrogen gas storage.
This work therefore aims at assessing the hypothetical application of porous-media H 2 gas storage, utilizing an existing anticline structure in the North German Basin that provides a realistic geological setting and parameters. A likely usage scenario is derived using a projection of typical fluctuations of renewable energy production and the energy demand to define the storage demand, and then simulating the hypothetical storage operation. At the selected storage site, the Rhaetian sandstones are chosen as the designated storage formation, which are overlain by Middle Jurassic and the Lower Cretaceous deposits that represent possible cap rocks (Hese 2011(Hese , 2012. At the storage site, formation thickness is up to 20 m, with the formation depth being about 450 -500 m. A heterogeneous and realistic parameter set consisting of 25 equally likely realizations is used to investigate the effect of the reservoir heterogeneity on the storage operation. Required storage dimensions and retrieval rates are determined for the ensemble. In addition, the hydraulic impact of the storage operation due to induced pressure perturbations is investigated.

Definition of the storage demand
For a scenario analysis of subsurface energy storage, an estimation of the storage demand is necessary. This can be broken down into the required power output that has to be sustained for a given period of time and the frequency with which such periods of demand occur due to the fluctuating nature of the renewable energy production. Storage demand is thus mainly controlled by the fluctuating renewable power production, as well as the deployment of loadbalancing schemes, which reduce the storage requirements by shifting energy demands on smaller timescales of up to hours. Numerical simulations of renewable power production in Germany, based on actual meteorological data in a scenario in which 100% of the electricity is produced by renewable sources, indicate that power shortage periods of up to 14 days in a year may occur frequently (Klaus et al. 2010). However, it can be expected that several storage options or back-up systems, such as conventional power plants, will be employed during such prolonged shortage periods, which essentially reduces the time span over which one storage has to operate at its maximum capacity. Thus, a typical shortage period of 1 week is assumed in this analysis, but the analysis can easily be extended to other shortage scenarios.
For determining the required delivery rate from porous-media H 2 storage, the power delivered is compared to the average weekly demand of electricity for the state of Schleswig-Holstein, northern Germany, with a population of about 2.8 million people. In 2011, a total of 42.82 × 10 6 GJ of electricity was consumed in Schleswig-Holstein (MELUR 2013). Thus, a week-long shortage period would result in a deficit of 0.82 × 10 6 GJ (228 × 10 6 kWh) of electricity, corresponding to an average load of 1356 MW. Klaus et al. (2010) estimated the round-trip efficiency of a H 2 storage system to be around 42%, incorporating conversion of electric power to gaseous H 2 and from H 2 back to electric power. The analysis in this paper is focused on the deliverability of power from the porous-medium hydrogen storage and thus only the re-electrification process is considered, which can be technically accomplished either through fuel cells (Büchi et al. 2014) or gas turbines (Forsberg 2009). The efficiency of the re-electrification process is assumed to be 60%, as this efficiency can be achieved with either method. Assuming an energy density for H 2 of about 124 MJ kg −1 (Carden & Paterson 1979) and the corresponding H 2 density of 0.085207 kg m −3 at surface conditions defined as 1 bar and 15°C, the volume of H 2 that has to be stored to cover for the complete deficit of 0.82 × 10 6 GJ in a week-long shortage period equates to 129 × 10 6 m 3 of H 2 gas at surface conditions (Sm 3 ).

Structural and geological model
Several criteria have to be fulfilled at a given site in order to qualify a geological structure as a potential gas storage formation. Besides providing sufficient reservoir volume, a potential storage site should provide competent sealing formations above and below the storage formation in order to prevent gas migration into adjacent formations. Also, a high intrinsic permeability is required to ensure well deliverability (Bennion et al. 2000). In addition, the depth of the storage formation is important as it affects the operational pressure range of the storage and the well deliverability. While deeper formations potentially allow a broader operating pressure range, reducing the well length by choosing a shallow storage formation reduces the pressure drop occurring within the well during gas flow and thus increases the overall well deliverability (Carden & Paterson 1979;Wang & Economides 2009). H 2 gas is prone to fingering in the porous formation as a result of its specific properties, thus a steeply dipping structure is favourable (Paterson 1983).
The anticline used in this analysis ( Fig. 1) is based on an actual structure found in Schleswig-Holstein, northern Germany. The storage operation, however, is purely hypothetical. The subsurface of northern Germany is strongly affected by movements of the Zechstein salt deposits (Doornenbal & Stevenson 2010). Halokinesis of these salt deposits, which led to the formation of the structure, started in the Triassic period (Baldschuhn et al. 2001). The Rhaetian sandstones of the Upper Triassic (Exter Formation) have been investigated as a potential host formation for CO 2 sequestration (Hese 2011(Hese , 2012. Furthermore, the overlying deposits of the Middle Jurassic and the Lower Cretaceous are regarded as possible cap rocks (Hese 2011(Hese , 2012. Apart from being investigated as potential CO 2 storage formations, the Rhaetian deposits are proven reservoir formations for natural gas exploration in NW Germany (Fahrion & Betz 1991). Several NW -SEorientated faults cut through the centre of the anticline as a result of changes in the local stress state. In this paper, the faults are assumed tight and represented as a hydraulically closed model boundary.
During the Rhaetian, the depositional system of the North German Basin changed from the non-marine environment of the Late Triassic to the marine environment of the Early Jurassic (Doornenbal & Stevenson 2010). As a result, the depositional system changed spatially from a non-marine system in the east through a paralic system to a marine setting in the west of the North German Basin (Doornenbal & Stevenson 2010). The study site is assumed to lie in the paralic facies belt, consisting of shallowmarine, deltaic and coastal sandstones, and separating the non-marine facies in the east from the marine facies to the west (Doornenbal & Stevenson 2010).
In general, the Rhaetian deposits are made up of several distinct sandstone layers of varying thickness with intermediate shale layers (Gaupp 1991;Doornenbal & Stevenson 2010;Hese 2011Hese , 2012. Multiple classifications of the Rhaetian appear in the literature, of which most include a separation into the Lower, Middle and Upper Rhaetian (Battermann 1989;Gaupp 1991;DSK 2005). The Lower Rhaetian consists of cyclic coarsening-upwards successions of mudstones to sandstones (Gaupp 1991). The Middle Rhaetian, also known as the Conterta succession, comprises five units: the Lower Shale, the Main Sandstone, the Middle Shale, the Flaser Sandstone and the Upper Shale. Similar to the Lower Rhaetian, the Upper Rhaetian shows a coarsening-upwards succession from mudstones at its base to sandstones in the upper parts of the formation (Gaupp 1991;DSK 2005;Doornenbal & Stevenson 2010). The Main Sandstone of the Middle Rhaetian has the highest potential for use as an underground gas storage site, as it shows, in general, a thickness of several metres and consists of medium-to coarsegrained relatively pure quartzite (Fahrion & Betz 1991;Gaupp 1991;DSK 2005;Hese 2011Hese , 2012. At the given site, the Main Sandstone exhibits thicknesses of around 15 -20 m near the top of the structure and thicknesses greater than 30 m at the flanks (Fig. 1).
A regional structural model developed by Hese (2012) is used as a basis for the geological site model created for this study. In order to get a more realistic presentation of the storage formation, the internal succession of the Rhaetian was included in the model. For this, the Lower Shale of the Middle Rhaetian, the Main Sandstone of the Middle Rhaetian, and the Upper Shale of the Middle and Upper Rhaetian were integrated into the structural model, based on facies descriptions given by Gaupp (1991) and Hese (2011Hese ( , 2012. At the storage site used in this analysis, the depth of the Rhaetian deposits varies from around 400 m near the top of the structure to more than 3000 m at the flanks of the anticline. Halokinesis, as well Owing to its shallower depth, the eastern flank of the structure is selected to accommodate the hydrogen gas storage site. Assuming that the faults in the centre of the anticline are tight, they represent a closed boundary and the simulation model was hence reduced to the eastern part of the structure. The heterogeneous geological model of the storage site is obtained by facies modelling based on the described structural model, and accounting for lateral and vertical trends. Doornenbal & Stevenson (2010) determined the major direction of sediment influx during the Rhaetian for the study region to be around 70°. While deltaic facies such as those found at the study site can extent several tens of kilometres (Morse 1994), the lateral correlation lengths were set to 2500 m in the major direction and to 500 m perpendicular in the minor direction in the absence of direct site-specific data. The facies descriptions given by Gaupp (1991) and Hese (2011Hese ( , 2012 were used to obtain correlation lengths corresponding to the coarsening-upwards of the individual sub-formations. The data were used to create 25 equally probable facies distributions. The final heterogeneous parameter distributions were obtained in a second step by assigning the parameters based on the facies distributions in combination with the parameter ranges assigned to each individual facies component (Table 1). One distribution of permeability for the geological model is depicted in Figure 2, which shows the vertically evenly distributed values in the storage formation, the Main Sandstone, compared with the stronger coarsening-upwards trends in the remaining sub-formations.
Site-specific data on hydraulic properties of the Rhaetian deposits are scarce. Parameter ranges for porosity and intrinsic permeability are available as on-site (Hese 2011(Hese , 2012 and off-site data from the North German Basin (Gaupp 1991). Capillary pressure and relative phase permeability data for the Raethian, however, are not openly accessible. Hence, the necessary parameters, such as residual saturations, displacement pressures, maximum gas-phase permeabilities and pore distribution indices, were derived from values typically found in the literature for rocks of similar characteristics (e.g. Hildenbrand et al. 2002Hildenbrand et al. , 2004Bachu & Bennion 2008;Wollenweber et al. 2010). Based on this the relative phase permeability data were calculated using Corey-type equations by Burton et al. (2009), while capillary pressure was calculated using the standard Brooks & Corey formulation (Brooks & Corey 1964). To account for difference in phase properties, the capillary entry pressure data taken from the literature were scaled to the phases used in the simulation using the Laplace equation (e.g. Helmig 1997). Assuming temperature and pressure conditions that can be expected during storage (25°C; 50 bars), the interfacial tensions calculated for H 2 and N 2 using the empirical formulation by Massoudi & King (1974) are 0.068 and 0.071 N m −1 , respectively. Consequently, the entry pressures for both gas components are comparable and thus phase composition was not taken into account for calculating the capillary pressure data.
The realizations generated here are used to obtain a first estimate of the effects of geological formation sub-structure and heterogeneity on the storage operation, especially the pressure evolution and the gas-phase distribution, as well as the well deliverability rates. They are not intended as a full-scale uncertainty assessment, as this would require more knowledge on site-specific data, as well as more realizations.

Numerical simulation model
The geological model of the storage site was transferred to a simulation model using a discretization of 50 × 50 m laterally, 0.2 -5 m vertically and a constant number of layers for each subformation. Initially, a hydrostatic pressure distribution and a fully water-saturated pore space is assumed throughout the model domain. The initial pressure distribution in the central plane of the storage formation is shown in Figure 3, where the steeply dipping anticline suitable for H 2 storage can be discerned from the hydrostatic pressure distribution. In order to get a representative pressure response in the far field of the storage, a pore volume multiplier and a permeability factor equal to the reciprocal of the pore volume multiplier were applied to the boundary elements of the simulation grid to represent the lateral extent of the storage formation. Based on the distribution of the Rhaetian deposits given in Baldschuhn et al. (2001), the extent towards the north, east and south were set to 48, 20 and 20 km, respectively. The boundary to the west is assumed to be closed to the fault system apart from the most southern part of the western boundary in which the storage formation continues beyond the model boundary for 14 km. Furthermore, the cap rocks above and below the Rhaetian deposits are assumed impermeable, and are thus omitted from the simulations. The storage operation is simulated using five wells, labelled 1 -5 from north to south (see Fig. 1). Preliminary simulations using homogeneous test cases showed that the individual screen length should be at least 12 m in order to sustain the storage rates required in this scenario. Placement of the storage wells can furthermore affect the theoretical maximum storage capacity of a structure, as shown for CO 2 storage by Mitiku & Bauer (2013), with the highest usable potential being available when the wells are placed at the top of a structure. Thus, the storage wells are placed at the shallowest depth possible near the top of the structure while at the same time fulfilling the minimum reservoir thickness criterion of 12 m. The individual reference depths of the storage wells are 493, 481, 479, 463 and 474 m for wells 1 -5, with Placing the wells further down the anticline would allow for higher injection and extraction rates, but also require a larger gas storage volume.
The necessity of cushion gas injection prior to the storage operation represents a costly one-time loss for porous-medium gas storage (Carden & Paterson 1979;Evans & West 2008). Using a cheaper alternative, such as inert N 2 or CO 2 , instead of the working gas can thus reduce the cost of the storage, and has been discussed and already successfully conducted for natural gas storage (e.g. Laille et al. 1988;Dussaud 1989;Oldenburg 2003), as well as for compressed air energy storage (Oldenburg & Pan 2013). The use of N 2 as a cushion gas in porous-medium H 2 storage also helps to reduce the extremely sharp density contrast between the H 2 and the formation water. On the other side, using a separate cushion gas does introduce the risk of impurities of the produced gas phase. Experiments on the effect of added N 2 on the combustion characteristics of H 2 -rich synthetic gas in gas turbines, however, indicated a relatively low overall impact (Lee et al. 2012).
The model simulations are carried out using the multiphasemulticomponent reservoir simulator ECLIPSE E300 (Schlumberger). The simulation takes into account two-phase immiscible flow, as well as pressure-dependent dissolution of both gas components in the water phase at the reservoir temperature. Evaporation of water into the gas phase is not included in the simulations. The solubility data of H 2 and N 2 was determined prior to the simulation using solubility data from Young (1981) and Battino (1982). Phase densities and viscosities are calculated using a generalized formulation of the Peng-Robinson equation of state (Schlumberger 2014). For a temperature of 25°C and pressures from 30 to 70 bars, which are representative conditions at the simulated storage site, the equation of state yields densities to within 1.9 and 1.3% of the values given in Lemmon et al. (2016) for H 2 and N 2 , respectively. The parameters used are given in Table 2. The calculated viscosities deviate by less than 13.9 and 1.9% from the reference values for H 2 and N 2 in Lemmon et al. (2016), respectively. Nevertheless, the obtained values are within the uncertainty interval specified in Lemmon et al. (2016). Diffusion and dispersion processes between the two immiscible fluids or the gas components are neglected in this study. In a porous geological formation, diffusion is reduced compared to free atmosphere diffusion as a result of high pressure, as well as porosity and tortuosity effects (e.g. see Oldenburg et al. 2004). With the given parameters of the storage site, the diffusion coefficient of H 2 would be about 3.6 × 10 −7 m 2 s −1 , corresponding to a diffusion length of less than 2 m during one storage cycle. As hydrogen gas moves by about 100 m during one cycle due to injection and extraction, we consider this effect small and neglect it. Mechanical dispersion, which represents mixing due to subscale hydraulic heterogeneity not explicitly represented in the model, is a macroscopic, scale-dependent effect. This process is also not included in the simulation, as the major contribution to mixing, being the heterogeneity of the hydraulic permeability of the formation, is explicitly accounted for in the model approach. In the case of a homogeneous reservoir representation, Feldmann et al. (2016) showed that effects of dispersion should be considered.  Critical temperatures (T crit ) and pressures (P crit ) are taken from Lemmon et al. (2016), critical volumes (V crit ) are from Kaye & Laby (2016), and acentric factors (Ω) are from Gasem et al. (2001). The values for the binary interaction coefficients are defaulted to 0.

H 2 storage in a heterogeneous sandstone formation
For each of the 25 realizations, the storage operation is simulated in three stages. In the first stage, the gas storage site is initialized by injecting the N 2 serving as the cushion gas for 710 days with a target injection rate of 55 625 Sm 3 /day/well. The time span of 710 days required for the initial filling of the reservoir compares well with reported times by, for example, Dussaud (1989), and shows that this initial period requires considerable effort. This is followed by the initial filling of the reservoir with H 2 for another 210 days with a target rate of 155 000 Sm 3 /day/well. Subsequent to the first two stages, the actual storage operation consisting of six storage cycles is simulated. Each cycle consists of a production phase of 1 week at a target rate of 1 000 000 Sm 3 /day/well followed by a short 1 day shut-in period. Consequently, the simulated storage should suffice to provide about 370 MW under ideal conditions, which is around 27% of the total demand defined in the scenario. The storage then is replenished with H 2 gas for 50 days with the target rate set to 150 000 Sm 3 /day/well, followed by a 30 day shut-in period. Thus, the frequency of the extraction periods is roughly 3 months, representing a strong withdrawal period every 3 months, with a longer time required to replenish the extracted gas by H 2 produced from surplus power. The BHP of the wells is limited to the range of 30 -65 bars, equating to ±50% of the initial pressure value. If the storage operation results in a violation of these limits in one of the wells, the injection/extraction rate applied in this well is automatically adjusted so that the BHP limit is respected.

Storage initialization
The consecutive injection of N 2 and H 2 gas in the first two stages of the simulation results in a quick pressure increase within the storage formation (Figs 3 and 4), with the maximum overpressure being 20 bars, as governed by the BHP limits of the storage wells. The pressure responses of the individual wells strongly depend on the distribution of the high-and low-permeability zones at and close to the injection wells in the individual heterogeneous realizations (not shown). Ultimately, the upper BHP limit of 65 bars is reached in all wells in every realization, resulting in an automatic reduction of the well flow rates. The pressure continues to increase over the course of the initial filling stage of the simulation at the observation points (see Fig. 1 for observation point positions), with the signal at 500 m reaching a maximum of slightly less than 12.5 bars. With increasing distance to the injection points, the overpressure signal is dampened considerably (Fig. 4). Furthermore, the arrival time of the pressure signal at the observation wells at 2500 and 5000 m is delayed, resulting in the maximum pressure being observed after the initial filling is completed and the first storage cycle already being well under way. The variability in the pressure signal from the 25 realizations is more or less constant throughout the first stage of the simulation, with a range between the 5th and 95th percentiles of about ±1.5 bars around the median at the observation points at a distance of 500 and 1000 m from the well gallery, and of about ±1 bar around the median at the distant observation points.
At a distance of 5000 m from the well gallery, the median of the maximum pressure perturbation caused by the storage reaches 3.1 bars, with the 5th and 95th percentiles being 2.6 and 3.8 bars, respectively.
The injected gas accumulates at the top of the structure due to its lower density compared with the formation water (Figs 5 and 6). The gas-phase composition exhibits a clear zonation as the N 2 is injected prior to the H 2 . While the N 2 is spread throughout the stored gas phase, reaching up to 2 km from the injection wells in a northern direction, the H 2 gas is concentrated around the injection wells (Fig. 5). The distribution of the gas-phase saturations in lateral and vertical directions indicates the strong influence of the formation  W. T. Pfeiffer et al. heterogeneity (Figs 5,6 and 7). As result of the upper BHP limit being reached in all wells during the storage initialization, the median of the injected N 2 and H 2 volume is just 82.75 × 10 6 and 21.75 × 10 6 Sm 3 , respectively. This represents just about 42 and 13% of the intended injection volume of N 2 and H 2 . Consequently, less H 2 gas than intended is available for extraction in the first storage cycle (Table 3). This clearly shows in the composition of the gas phase in the storage prior to the first extraction period, with the only significant H 2 fractions in gas being found in close proximity to the wells (Fig. 5).

Storage operation
The third stage of the simulation consists of the actual storage operation represented by six storage cycles. The production of gas from the storage formation results in sharp and pronounced pressure drops in and around the storage wells by more than 20 bars. This results in pressure levels of more than 10 bars below the initial hydrostatic levels, while local pockets of overpressurized fluid are persistent in some realizations over the course of the extraction periods (Fig. 3). The minimum theoretical pressure that can occur is again governed by the BHP limits of the storage wells and thus equates to 30 bars. The pressure change observed at the observation points 500 and 1000 m distant from the well gallery during the first storage cycle from day 923 to day 930 exhibits a similar characteristic compared to the storage wells with the median pressure dropping by about 14 and 8 bars, respectively (Fig. 4). Consequently, the pressure at 500 m is already below the initial hydrostatic pressure at the end of the extraction period, while the observation well at 1000 m still shows a slight median overpressure of about 1.8 bars. The variability of the pressure change as determined by the spread between the 5th and 95th percentiles around the median varies but stays within ±2.7 and ±1.0 bars at the observation points at 500 and 1000 m distance, respectively. The more distant observation points again show a delayed pressure response to the gas extraction and smaller amplitudes (Fig. 4). The refill of the storage with H 2 from day 932 to day 982, subsequent to the gas extraction, again results in a pressure increase in the vicinity of the storage wells (Fig. 3). Even though the amplitude of the pressure increase is greater than during the initial fill, the maximum value of the median pressure observed is lower than the value reached during the storage initialization, as the pressure level prior to the refill is lower than the initial pressure (Fig. 4). The rate of the pressure change during the storage refill is lower than during the production phase because the gas flow rates applied at the wells, and thus the gradient imposed on the system, is considerably smaller. The variability in the pressure signal during the storage refill periods is nearly identical with the variability observed during the production phases. During the subsequent 30 day shut-in phase, the pressure slowly declines towards the initial

321
H 2 storage in a heterogeneous sandstone formation value in the injection wells (not shown). This behaviour is also visible at the observation points at 500 and 1000 m, with a pressure decrease of about 1 and 0.5 bars, respectively. Contrary to this, a slight pressure increase resulting from the H 2 refill can be observed at 2500 and 5000 m during this period (Fig. 4). The subsequent pressure drop at 1020 days is already a result of the second production period. Thus, the pressure signal of the 30 day shut-in period separating the storage cycles is not visible at the distant observation points as the strong pressure gradient imposed on the system during production outpaces it. The variability of the pressure signal during the shut-in period is similar to the refilling period.
With an increasing number of storage cycles, the pressure signal from the initial filling of the storage dissipates slowly as a result of the large pore volume boundary conditions (Fig. 4). The maximum and minimum pressures observed during the storage cycles also decrease, while the variability of the pressure signal is more or less constant throughout the simulation.
While both gas injection and extraction are clearly visible in the pressure distribution, few or no visible differences are observable in the gas saturation (Figs 5 and 6), indicating a relatively stable gas-water contact during the storage cycles. This shows that the storage operation is mainly supported through gas expansion and compression and is therefore controlled by the gas compressibility. This conclusion is further supported by the observed changes in gas density of around 20% during each storage cycle. The total gas in place increases with the number of storage cycles as slightly more gas is injected than extracted during each cycle (Table 3; Fig. 5). Nevertheless, the strong effect of the formation heterogeneity is persistent throughout the simulation, with the distribution of the gas phase being highly variable between the individual realizations until the last storage cycle (Fig. 7). While few realizations show a relatively homogeneous gas-phase distribution (e.g. run #15 in Fig. 7), most realizations display a more variable pattern. Well placement in a heterogeneous reservoir can thus result in individual wells being more or less separated from the main storage (e.g. centre well 3 in run #3 and #19 in Fig. 7). This well-known effect has to be checked for in the field during well testing.
The composition of the produced gas phase changes considerably during each extraction period (Figs 5 and 7). As more H 2 is injected than extracted in each storage cycle, the H 2 in place increases with Fig. 6. (a) Absolute gas saturations before the sixth production period and (b) the change in gas saturation after the sixth storage cycle; (c) molar fractions of H 2 in gas before the sixth production period and (d) the change in the molar fractions of H 2 in gas after the sixth production period (magnitudes) in the storage formation in run #14. The displayed formation is vertically exaggerated by a factor of 5. Each grid block is 50 × 50 m in lateral direction. the number of cycles (Fig. 5). While the absolute molar fractions of H 2 in gas show high values at the top of the storage formation (Fig. 7), the vertically averaged values are considerably lower in most of the storage site (Fig. 5). Thus, a density-contrast-driven vertical separation of the two gas components is occurring, with H 2 occupying the upper parts of the reservoir and N 2 making up the gas phase in the lower parts of the storage formation.
During the first extraction period, the lower BHP limit is reached in at least one well in almost every realization, resulting in reduced well flow rates and thus reduced gas extraction rates for the whole storage (Table 3; Fig. 8). The median of the achieved gas flow rate is above 4 × 10 6 Sm 3 /day for most of the first extraction period, ultimately dropping down to about 3.5 × 10 6 Sm 3 /day towards the end. The variability in the storage gas flow rate is large, with the range between the 5th and 95th percentiles being about ±1.4 × 10 6 Sm 3 /day around the median value. The composition of the gas phase in storage and also that of the produced gas changes considerably during each production period (Figs 5,7 and 8). While the molar fraction of H 2 in the produced gas phase initially is 0.8, fractions of about only 0.3 are reached at the end of the extraction period (Fig. 8). The variability in the data stays more or less constant and within ±0.05 of the median. Consequently, the median of the achieved H 2 flow rate drops considerably from initially 3.5 × 10 6 Sm 3 /day to slightly less than 1.1 × 10 6 Sm 3 /day at the end of the production phase, equating to about 260 and 78 MW, respectively (Fig. 8). The median of the total power output over the whole production period equates to 13.64 × 10 6 Sm 3 or 24 104 MWh (Table 3), which is about 11% of the required power specified in the storage demand analysis. However, taking the median power output at the end of the production phase of 78 MW as a guideline for the potential continuous power output of the storage, the storage can only supply 5.7% of the continuous power demand of 1356 MW specified in the scenario.
In the subsequent storage cycles, the median of the achieved storage flow rates increases slightly, with the lowest value obtained at the end of the last storage cycle being 4.0 × 10 6 Sm 3 /day (Fig. 8).
With the exception of the last day of the third storage cycle, the 95th percentile is constantly at a maximum value of 5.0 × 10 6 Sm 3 /day throughout the third and sixth storage cycle. Thus, the best 5% of all simulation runs can sustain the envisaged extraction rate of 1 × 10 6 Sm 3 /day/well from the third storage cycle onwards. The variability in the storage flow data decreases slightly with the number of storage cycles. Compared with the first storage cycle, the median value of the molar H 2 fraction in the produced gas, and thus the H 2 flow rate of the storage, increases considerably in the later cycles, with the median being 3.32 × 10 6 Sm 3 /day at the end of the sixth cycle, equating to about 245 MW (Fig. 8). This means that the storage is capable of continuously delivering 18% of the required power in this scenario in the last storage cycles. The corresponding cumulative power output of the storage equates to 52 022 MWh, which is nearly 23% of the total demand specified in the scenario. Obviously, the storage performance is greatly increased compared to the first few storage cycles (Table 3).
The variability of the H 2 flow data is increased compared to the first cycle, with the 5th and 95th percentiles ranging from 2.19 × 10 6 to 4.14 × 10 6 Sm 3 /day in the sixth production period, providing about 162 and 305 MW, respectively. Thus, the minimum and maximum power output that can be sustained over the whole sixth production period equates to 11.9 and 22.5%, respectively, of the required delivery rate of power. The variability in the data of the H 2  323 H 2 storage in a heterogeneous sandstone formation flow rate is mainly due to the variability in the total gas flow rate as the data of the molar H 2 fractions in the produced gas show only a small spread (Fig. 8). Thus, the storage performance strongly depends on the distribution of the gas saturation in the storage formation. An approximately homogeneous gas-phase distribution as observed, for example, in run #15 (Fig. 7) allows for high H 2 gas flow rates and thus a power output twice as large as for a strongly varying phase, such as that in run #3 (Fig. 7). Thus, the performance of the storage site is mainly governed by the distribution of the H 2 gas phase within the storage formation, which, in turn, depends strongly on the formation heterogeneity.

Conclusions
The heterogeneous parameter distribution of the storage system affects the behaviour of the storage in the short and long term. In the first two stages of the simulated storage operation (initial fill with the cushion and working gases), the heterogeneity of the hydraulic parameters results in reduced injection rates as a consequence of the pressures in the injection wells reaching the specified upper BHP limit. Even if good injection properties of the reservoir are given, the volume of H 2 that can be provided for the first storage cycle is lower than intended. This, in part, can be attributed to the injection of inert N 2 in the first stage prior to the first H 2 injection, resulting in already elevated pressures in the storage formation and thus reducing the amount of H 2 that can be injected in the second stage. Changing the injection pattern, such that the N 2 is injected at the side of the structure purely for pressure mitigation after the initial filling with H 2 , could potentially improve the storage performance of the first cycles. The N 2 , however, is also acting as an inert barrier between the formation water and the H 2 , which could help in minimizing chemical and microbial reactions. Other techniques, such as additional shut-in periods during the storage initialization to reduce the formation pressure, could also help to increase the H 2 volume available prior to the first storage cycle. The relatively slow dissipation of the overpressure signal, however, requires long shut-in times during which the storage cannot be used. Starting the storage operation, Fig. 8. (a) Total storage flow rate, (b) volume averaged molar H 2 fraction in the produced gas phase per storage and (c) resulting H 2 flow rate, as well as equivalent power output, for storage cycles 1, 3 and 6. All observation points are located at the central layer of the storage formation. The solid black line depicts the median of all realizations, the dark grey shaded area is the interval spanning between the 25th and 75th percentiles, and the light grey shaded area is the interval spanning between the 5th and 95th percentiles. The dashed lines indicate the absolute minimum and maximum values.

324
W. T. Pfeiffer et al. even if the performance is sub-par in the first cycles, does result in comparable or even better storage performances after the same time span with the added benefit that the storage is, albeit on a smaller scale than intended, already in use. The heterogeneous distribution of the gas phase in combination with the varying composition results in reduced storage flow rates, as well as N 2 gas as an impurity in the produced H 2 gas. Furthermore, the simulation of the storage operation clearly shows that the storage is not operating at its full potential in the first storage cycles, as the extracted H 2 volume and thus the power output increases consistently with the number of completed storage cycles. The increase in storage performance is due to the fact that more H 2 is injected than extracted in each cycle.
Even though the achieved H 2 extraction rates increase with the number of completed cycles, the variability among the heterogeneous realizations stays more or less constant. This is due to the gas-phase distribution being strongly spatially heterogeneous throughout the whole simulation. Consequently, the explicit simulation of the storage initialization as opposed to defining an initial homogeneous gas-phase distribution is very important in order to obtain valid storage flow rates later in the simulation. Monitoring of the gas-phase distribution by, for example, a combination of seismic and hydraulic methods, as reported in Pfeiffer et al. (2016) or Benisch et al. (2014), could thus enhance the understanding of the storage formation and help to optimize the storage operations.
As every simulation run exhibits a decrease in the H 2 fraction in the produced gas in the week-long production period and most simulation runs also show a reduction in storage flow rate, the H 2 flow rate and thus the power output from the storage is not constant. Taking the lowest value as a reference, it can be concluded that the simulated storage is capable of supplying a continuous power output of around 245 MW for 1 week when using five storage wells, while peak performance can be as high as 363 MW. Thus, the storage is sufficient to provide 18% of the total delivery rate required in the defined scenario, which assumes no power production from renewable sources at all during the 1 week period. If favourable parameter distributions and, thus, homogeneous, high saturation gas-phase distributions close to the wells prevail, the storage formation can supply up to 300 MW continuously, equating to about 22.5% of the required rate. Increasing the power output requires more wells to be used, as the storage performance is clearly rate limited. Storage capacity is not limiting, as the spill points of the structure are not yet reached by the injected gas, so that the stored gas volume could be increased easily. Dispersion effects, which are not included in this simulation via dispersivities but rather by explicitly representing formation heterogeneity, can be used to represent mixing processes during gas storage and may have considerable impact on the component distribution in certain settings, as found, for example, by Feldmann et al. (2016). A site-and model-approachspecific evaluation of these mixing effects should thus be conducted, which requires, however, site-specific determination of dispersivities. In cases of higher mixing than represented in this work, fractions of hydrogen in the produced gas may slightly decrease.
The pressure response of the reservoir to the storage operation is far reaching, with a pressure increase of about 3 bars being observed at distances of 5 km and the storage operation being visible in the recorded pressure signals. These signals could thus be used to monitor the storage site, and the injection and extraction periods, as, for example, Benisch & Bauer (2013) did for CO 2 storage. Pressure changes at such distances from the storage wells caused by other usages, however, will most probably not affect the storage operation drastically as the observed storage-induced pressure change is already relatively low.