From Legacy Gas Fields to Hydrogen Storage: 3D Seismic-Driven Geological Modelling and Dynamic Simulation in the Northern Upper Rhine Graben

1. Introduction

The transition toward sustainable and low-carbon energy systems is gaining momentum worldwide, driven by the environmental impacts of conventional energy sources and the need to reduce greenhouse gas emissions. In this context, hydrogen (H₂) has emerged as a promising energy carrier; however, large-scale and efficient storage remains a key challenge. Geological storage options such as salt caverns, depleted hydrocarbon fields, saline aquifers, and lined rock caverns are considered potential solutions [1], yet their applicability strongly depends on regional geology and integration within the hydrogen supply chain.

Depleted hydrocarbon reservoirs were the first and most widely used storage sites due to their proven sealing capacity and existing infrastructure, while saline aquifers were later introduced as an alternative, albeit with higher uncertainty and operational complexity [1]. Salt caverns were increasingly developed as a storage option characterised by high deliverability and rapid injection–withdrawal capability, thereby complementing the large working gas (WG) capacities of porous reservoirs such as depleted hydrocarbon fields and aquifers, which are typically used for seasonal storage, whereas salt caverns provide high deliverability and short-cycle flexibility [2,3]. There are four known examples of underground subsurface storage for pure hydrogen, often mentioned in scientific literature: (i) Teesside in Yorkshire, UK, (ii) Clemens, USA, (iii) Moss Bluff, USA and (iv) Spindletop, USA [4]. Though storing hydrogen in salt caverns is already an established technology, porous reservoirs need further assessment. Recently, more options have been subjected to demonstrations or feasibility studies to establish the capabilities of subsurface hydrogen storage. Several pilot projects worldwide, including HyStock (Netherlands), Underground Sun Storage (Austria), HyStorage (Germany), Hychico (Argentina), Hybrit (Sweden) and TH2ECO (Germany), have demonstrated the technical feasibility of underground hydrogen storage (UHS) in salt caverns, depleted gas reservoirs, and lined rock caverns, respectively, highlighting the complementary roles of different storage types in future hydrogen energy systems [1,5,6,7,8,9,10].

In Germany, the northern parts of the country benefit from extensive salt formations and mature underground gas storage (UGS) infrastructure, well connected to offshore wind generation and pipeline networks. The Zechstein salt formations found here are thick and extensive, featuring salt domes and diapirs. As these formations extend southward, they become thinner and consist of stratified evaporite layers that are generally unsuitable for cavern development [11,12,13]. Central and southern Germany lack suitable salt deposits but host numerous legacy hydrocarbon fields and porous reservoirs with strategic storage potential for H₂.

This study focuses on the northern Upper Rhine Graben (URG), where suitable geological salt bodies are absent but porous reservoirs are widespread. In the absence of suitable salt formations, porous reservoirs—particularly depleted hydrocarbon fields and deep saline aquifers—are regarded as the most prospective sites for underground hydrogen storage. Furthermore, depleted gas fields converted into UGS facilities within these formations provide operational experience and field-scale data to evaluate their suitability for hydrogen storage. The Stockstadt UGS represents a depleted hydrocarbon field, whereas the Hähnlein UGS is developed in a saline aquifer. For both sites, a comprehensive set of publicly available geological, petrophysical, and operational data allows the construction and calibration of a high-quality base reservoir model. Part of the data (seismic, wells, and operational data) is available, and the rest is proprietary to the industry partner. In contrast, data coverage for other legacy fields in the region is sparse or less reliable, limiting direct site-specific assessment. The primary objective of this paper is therefore to develop and validate an analogue base model that captures the essential geological and dynamic characteristics of porous reservoirs in this part of the URG and assess the potential of these UGS for future conversion into UHS. The calibrated model is intended to serve as a transferable reference framework for assessing UHS performance in nearby legacy fields with similar geological settings, where detailed data are currently unavailable.

2. Geology and Operational History of the Study Area

The study area is located in southern Hessen, Germany, within the URG. URG has also been a site for developing energy sources and carriers, such as hydrocarbons, geothermal energy, subsurface gas storage, lithium extraction and natural hydrogen. Several former hydrocarbon fields, such as Stockstadt, Frankenthal, and saline aquifer in Hähnlein, are currently utilised for UGS, making them promising candidates for future UHS. Due to the tectonic setting and depositional environment of the URG, most of these fields feature porous reservoirs at different stratigraphic zones that serve as structural traps, making them candidates for repurposing as contingent for UGS/UHS in central Germany to support future sustainability.

The URG is part of the European Cenozoic Rift System (ECRIS), which stretches for 300 km from Basel to Frankfurt and is linked to the Alpine orogeny (refer to Figure 1a). From Eocene to present, the URG has undergone a complex evolution with multiple tectonic phases [14,15].

The stratigraphic column (modified after [16], shown in Figure 1b) shows the lithostratigraphic units that host hydrocarbon reservoirs and the corresponding source rocks. The lithostratigraphic descriptions of the stratigraphic units in the study area are detailed in previous publications [16,17,18]. The gas accumulations in this study area, mainly found in Hydrobia Beds, Groß-Rohrheim and Weiterstadt/Iffezheim Fm, serve as gas reservoirs with biogenic natural gas origins. However, the Hydrobia Beds are minor as reservoirs, but their bituminous marl clays are crucial as the source of gas. The Miocene sands (Upper Tertiary layers) are also commonly known locally in German as JTI (Jungtertiär I) and JTII (Jungtertiär II). Hydrocarbon fields produced between the 1950s and 1980s are shown in Figure 1a, along with the cumulative production volumes in Table 1.

This paper focuses mainly on the reservoir sands of JTI (Groß-Rohrheim Fm) and JTII (Weiterstadt Fm/Iffezheim Fm) of Stockstadt, and Hähnlein UGS. Figure 2 illustrates a detailed top-depth structure map of the Upper Hydrobia Beds, drawn during the production period (1962) for this area, along with a semi-schematic cross-section (Figure 2a) that reflects the geological understanding at that time across different production sites and a log correlation panel (Figure 2b) from the Stockstadt area. The cross-section AA illustrates the two stratigraphic zones in Stockstadt that produced hydrocarbons. The deeper stratigraphic zone (Pechelbronn Beds) produced oil, and the shallow JTII sands produced biogenic gas. The log panel in Figure 2b, showing the correlation of the various sands, clarifies the thin sand beds and their intercalation with shale layers. According to operational data, well reports, and published literature [17,18,19], the water-bearing Sand 8 of the Hähnlein structure has been used as a storage facility since 1960, with coal gas (town gas) stored. Later in 1973, the reservoir was rebuilt to store natural gas. At approximately the same time, the eastern structure of the Stockstadt gas field was also converted to UGS. In Stockstadt, Sand 7 was converted to storage in 1963, followed by Sand 8 in 1970. Stratigraphically, Sand 7 and Sand 8 are part of JTII.

3. Methodology and Database

3.1. Methodology

The methodology follows a modified reservoir characterisation and simulation workflow adapted from conventional hydrocarbon field development and UGS studies. The workflow was designed to first establish a calibrated base reservoir model using historical production and UGS data and then adapt it for hypothetical UHS feasibility simulations. The overall workflow is illustrated in Figure 3 and comprises three main stages: reservoir model construction, dynamic model calibration through history matching, and UHS scenario simulation.

The first stage involves developing the reservoir model, which provides the geological and structural basis for the subsequent dynamic simulations. Depending on data availability, this model may be generated as either a simplified field-scale tank model or a detailed static geological model based on subsurface data. In this study, a comprehensive dataset was available for the Stockstadt and Hähnlein UGS sites. Well data, including logs and stratigraphic markers, were correlated and calibrated with seismic data to interpret the main faults and horizons. The interpreted structural elements were used to build a regional-scale geological model, from which a field-scale model was extracted. To improve computational efficiency, this field-scale model was converted into a tartan grid, thereby reducing the number of active grid cells while preserving the main structural and reservoir features. While static model building uses the software Petrel™ (SLB) and SKUA-GOCAD™ (AspenTech), subsequent dynamic reservoir modelling has been carried out with the GEM™ (Computer Modelling Group, CMG).

The second stage involves constructing and calibrating the dynamic flow simulation model for the historical production and UGS phases. The dynamic model was initialised by assigning reservoir pressure and temperature, gas–water contact, fluid composition, relative permeability functions, and operational schedules for the active wells. Instead of a black-oil simulator, a compositional model was used to explicitly include the properties of hydrogen and other gas components in the dynamic simulation. The model was first set up to reproduce the historical production phase, including gas production and pressure depletion. Subsequently, the UGS phase was simulated by incorporating injection and withdrawal operations and the associated pressure response. As shown in Figure 3, model calibration was carried out through an iterative history-matching process. The simulated production volumes and pressure response were compared against historical field data. The matching criteria included the actual production and pressure behaviour during the production phase, as well as injection and production volumes and pressure changes during the UGS phase. If the simulated results did not satisfactorily match the historical data, sensitivity runs were performed to identify uncertain reservoir parameters that were controlling the mismatch. These parameters were then adjusted during tuning runs until an acceptable match between simulated and observed rates and pressures was achieved. This iterative loop was essential for constraining uncertain properties, including aquifer extent, reservoir connectivity, permeability distribution, and pressure-support behaviour. Once the model reproduced the historical production and UGS performance with sufficient accuracy, it was considered a calibrated base case model for further UHS assessment.

The third stage consists of adapting the calibrated dynamic model for hypothetical UHS feasibility simulations. H₂ is added to the gas composition to include the hydrodynamic, thermodynamic, and compositional properties of hydrogen along with other gas components present in the system. In the UHS simulations, the remaining natural gas in the reservoir was treated as cushion gas (CG), or the choice of gas that can be injected as CG, while different WG compositions were tested to evaluate their influence on storage performance. The scenarios were designed to assess operational conditions, recovery efficiency, pressure response, water production, hydrogen flow behaviour, and the interaction between hydrogen and other reservoir components. The scenario simulations provide insight into the feasibility of using the studied reservoir system for hydrogen storage under shallow reservoir conditions. In addition to flow behaviour and production performance, the workflow allows further extensions to include geochemical reactions, microbial activity, and gas–brine–rock interactions where sufficient data are available. Based on the outcomes of the UHS scenario simulations, the model can also be used to support future assessments of infrastructure requirements, surface facilities, and economic feasibility. Thus, the workflow provides a structured approach for moving from a calibrated hydrocarbon/UGS reservoir model toward a field-scale evaluation of UHS potential.

3.2. Database

A comprehensive subsurface database was compiled for the selected study area by integrating well, seismic, production, and archival data from institutional and industrial sources. The selected study area is located in the northern URG in southern Hessen and includes several legacy hydrocarbon fields, such as Wolfskehlen, Stockstadt, Pfungstadt, and Darmstadt-SW. The location was chosen for this feasibility study based on comprehensive data, infrastructure, and UGS experience. Over 150 wells have been selected, with drilling reports from Hessian Agency for Nature Conservation, Environment and Geology (Hessisches Landesamt für Naturschutz, Umwelt und Geologie, short form as HLNUG) providing stratigraphic and lithological data. Additionally, 34 2D seismic lines were provided by MND Energy Storage Germany GmbH. 3D seismic data covering approximately 10 km x 30 km from Wolfskehlen to Worms, along with a velocity model, was provided by HLNUG. Production histories from the State Authority for Mining, Energy and Geology (Landesamt für Bergbau, Energie und Geologie, short form as LBEG) were added to the database, which provide insights into monthly production but lack specific details on reservoir pressure, production rates, or water cut during the production phase (1950s-1970s). Some archived reports about the area have been used to characterise the reservoir and reconstruct the production phase [17,18,19].

Lithological descriptions with well markers for the Stockstadt area have been loaded from the drilling reports and correlated with the well markers from Wolfskehlen and Hähnlein. The industry partner provided the top markers of the reservoir sands in the Stockstadt and Hähnlein areas, which were used to delineate the sands in the UGS areas. The lithological description of the reservoir section was transformed into a facies log for sand correlation and facies modelling. Basic log sets of gamma ray, resistivity, and spontaneous potential were used to identify the lithological character of the reservoir. Most of the wells were drilled during the 1950s-1980s period and lack a modern log suite. Table 2 summarises the stratigraphic names used in the HLNUG drilling reports, which provide detailed lithologic descriptions and the corresponding seismic horizons interpreted for the study area.

Markers were correlated with seismic reflectors to identify changes in lithology as indicated by the seismic data. The stratigraphic zones in the study area are primarily comprised of clastic sediments, which can be confirmed from well reports and geological descriptions from published literature. Nine seismic horizons were identified; refer to Table 2: JTII Top, JTI Top, OHY Top, UHY Top, CBS Top, BNS Top, RT Top, Tertiary Base, and Basement Top. The German nomenclature for the horizons was used for simplification because these names are widely used in the scientific community. JTII Top, JTI Top, and OHY Top were locally refined in Stockstadt and Hähnlein based on MND Energy Storage Germany GmbH’s 2D data interpretation. The fault interpretation in the 3D survey was conducted and compared with the 2D interpretation. Time horizons were converted to depth for static modelling, ensuring reliable depth accuracy for shallow horizons, which were tied to the well markers. In contrast, deeper horizons have reduced accuracy due to lateral velocity variations and limited calibration points, though the error remains acceptable.

3.3. Structural Model

A regional model was created using the nine interpreted seismic horizons in the depth domain, covering the areas of Wolfskehlen, Stockstadt, and Hähnlein. The structural framework considers faults that originate from the basement and affect the Quaternary sediments, including those that impact the reservoir zone. This regional model allows for the extraction of field-scale models for dynamic simulation.

A field-scale simulation grid which includes Stockstadt and Hähnlein was extracted from the regional model for history matching of the production and, UGS phases as well as, for the feasibility study of potential future UHS conversion scenarios. The cell size in the reservoir zone was 100 m × 100 m × 2 m, totalling 328,860 cells. Local refinement was applied to the UGS areas to capture changes in the dynamic results in fine resolution. In X and Y directions the cells were divided by 2, and in Z direction by 4, resulting in the resolution of 50 m × 50 m × 0.5 m with ~ 1.5 × 10⁶ cells in total.

A simplified approach for the facies model was considered, and two facies types were defined: reservoir facies (Sand) and non-reservoir facies (Shale). Though Sand 7 and Sand 8 were deposited in channel systems, hence the reservoir characteristics vary laterally and vertically, also the connectivity between the two sands is difficult to quantify. Lithological descriptions from wells in the JTII reservoir were converted into facies logs and upscaled to the 3D grid using the ‘Largest proportion’ approach. Facies continuity was assessed using vertical proportion curves (VPCs). The VPCs calculated layer-wise facies proportions across the grid using upscaled logs at well locations. Vertical variograms for sand were generated from well data; however, for the lateral variogram, assumptions were made considering the continuity and heterogeneity of the sand bodies based on the depositional environment. The facies proportion from the upscaled logs was used as hard data and the VPCs were used as secondary data to distribute the facies in 3D space using the Sequential Indicator Simulation (SIS) algorithm. The reservoir facies was assigned a constant porosity of 25-30% and a permeability of 1000 mD.

3.4. Dynamic Model for History Match of Production and Ugs Phase

The dynamic model includes the Stockstadt gas field and UGS Hähnlein, initialised on January 1, 1955, utilising LBEG annual reports for the production phase. Sand 7.2 was identified as gas-filled in the UGS Stockstadt, while sands 8.1 and 8.2 were water-filled. However, all the respective sands in the other fault blocks in the Stockstadt area were gas-filled. Sand 8 in the Hähnlein area was water-bearing on the initialisation date. To effectively manage the configuration of sands, the reservoir grid is divided into three regions to assign the initial conditions for reservoir pressure and gas-water contact (Figure 4). The initial reservoir temperature was 29 °C, with an initial pressure of 5100 kPa (51 bar) at a depth of 430 meters below MSL. The gas accumulation in the Stockstadt area was of biogenic origin, and it is believed to be a localised and small gas accumulation. Hence, the assumption was made that the gas saturation was quite low. Default relative permeability curves for sand from Petrel were used.

A five-component model was used, reflecting gas compositions (94.6% CH₄, 3.0% C₂₊, 2.2% N₂, and 0.2% CO₂) reported in [21] during production, with H₂ included for hypothetical UHS scenarios. Volume multiplier (VM) was applied to the cells below the gas-water contact to define the aquifer dimension. Data from LBEG spanning 1955 to 1978 indicated 26 producing wells in the Stockstadt area, with gas production rates estimated from their yearly outputs. The individual rates for the wells are not available due to data archival limitations. The LBEG report does not specify the water production rate; however, some older literature cites an average water production rate of approximately 100 m³/day at the end of the production phase [18]. The maximum gas production rate per well was assumed to be 6,000 m³/day, with a minimum bottomhole pressure of 500 kPa. Based on the assumed water production rate, the VM was adjusted to determine the aquifer dimensions, and three sizes were simulated for sensitivity analysis.

The UGS conversion of the eastern block of Stockstadt and Hähnlein started in 1960s. Sand 7.2 was transformed into a gas storage in 1963. Gas injection into sand 8.1 on the Stockstadt structure began in 1970, shortly followed by the start of injection into sand 8.2. The Hähnlein facility was first utilised for storage in 1960 when coal gas was injected into the structure. In 1973, the reservoir was modified to store natural gas [19].The model simulates UGS operations from March 2016 to April 2025. At the Stockstadt UGS, 12 storage wells are operational, while in the Hähnlein UGS, 15 are operational. The injection and withdrawal rates for this time duration were downloaded from the company’s official website. These rates were assumed to simulate the approximate conditions of the UGS to the present day. Due to the scope of the work and data confidentiality, actual pressure data, injection, and withdrawal rates for the UGS phase were not simulated; however, an attempt was made to replicate conditions as similar as possible to understand the lifetime of a depleted gas reservoir or saline aquifer.

Along with the remaining gas, 50 MMm³ was injected into Sand 7 as CG for UGS operations. Approximately 38 MMm³ and 80 MMm³ of natural gas, with the same composition [21], are injected as CG during the initial setup of Sand 8 at the Stockstadt and Hähnlein UGS sites over 1 year, from March 2014 to February 2015. Following this, a similar volume of WG (Figure 5) was considered in the schedule, along with similar operational conditions (refer to Table 3), which included the number of wells, bottom-hole pressures, withdrawal and injection rates, to simulate the interval from March 2016 to March 2025. Based on the net rock volume (NRV) of Stockstadt UGS, it is assumed that approximately 60% of the WG is from Stockstadt. This volume assumption was considered in the schedule file. The average injection rates for Sand 7 and Sand 8 in Stockstadt UGS are 183,000 m³/day and 61,200 m³/day per well, respectively. The average withdrawal rates for the respective sands are 193,000 m³/day and 65,000 m³/day. For Hähnlein UGS, the average injection/withdrawal rates are 59,000/89,500 m³/day, respectively.

3.5. Dynamic Model for Hypothetical Uhs Phase

After this UGS simulation stage, the model was considered as the base case model to test different UHS scenarios. UGS Stockstadt was simulated for 10 years for two cases: 1) 5% H₂ + 95% CH₄ (Case X) and 2) 100% H₂ (Case Y), where the remaining natural gas is considered as CG. The two scenarios were conducted over a 10-year period, from April 2025 to March 2034. The injection duration was derived from typical UGS use, with a main period from around April to late October or early November. Production occurs from November through the end of March. A total of 12 storage wells were used for operations (four for Sand 7 and eight for Sand 8). A group target for injection volume of 650,000 m³/day was set for seven months, resulting in a total injection of 139 million cubic meters (MMm³). During the annual cycle, a withdrawal volume of 925,000 m³/day was planned for five months, also totalling 139 MMm³. For Sand 7, the average injection rate is 90,000 m³/day per well, while the average withdrawal rate is 125,000 m³/day per well. In contrast, Sand 8 has an average injection rate of 37,500 m³/day per well and an average withdrawal rate of 73,000 m³/day per well.

4. Results

4.1. Seismic Interpretation and Structural Model

The detailed reinterpretation of the 3D seismic survey allowed for the generation of higher-resolution depth structure maps for the area of interest. This process also enhanced the understanding of the structural configuration of faults within the reservoir zone and in the surrounding regional setup. The shallow horizons (JTII Top, JTI Top and OHY Top) have minimum time-depth conversion error. The interval velocity, calculated from well markers and time picks at well locations, ranges from 2,000 to 2,200 m/s, which is applicable for shallow depth zones with clastic sediments. Data from checkshot analysis validate the observation regarding Two-Way Traveltime (TWT). Consequently, the marker calibration of shallow horizons in the velocity model resulted in an improved definition of reservoir structure and a more accurate depth prediction of the sand tops.

Figure 6a depicts the reinterpretation of the data, enabling more accurate identification of the 3D configuration of fault surfaces. Figure 6b displays a seismic section across Wolfskehlen-Stockstadt-Hähnlein, highlighting the interpreted horizons and faults. The Basement Top picking is incomplete due to a poorly defined seismic reflector laterally. Analysis of sediment thickness maps shows a general increase from NW to SE, with local variations driven by fault activity. The interval between OHY Top (18 Ma) and UHY Top (21 Ma) has been analysed and shows significant thickening in the eastern fault block, aligned with periods of increased fault activity. Fault-related depocentres confirm that differential fault movement has caused heterogeneity in sediment thickness and a complex relationship between tectonics and sedimentation (explained in detail [22]).

The nine depth-converted seismic surfaces were used to create a general structural model with the most important faults. Furthermore, faults were added to the reservoir grid to improve the compartmentalisation of the zone of interest. The faults are represented as stair-stepped faults to capture more complex structural configurations between two faults and facilitate faster simulation. Figure 7 shows lithological cross-sections through the Stockstadt and Hähnlein reservoirs extracted from the reservoir grid, illustrating the stratigraphic arrangement of the main reservoir units (Sand 7 and Sand 8) and the positions of the storage wells. The facies distribution represents a geologically reasonable realisation of a channelised sandstone system constructed during static modelling and constrained by well data and regional depositional understanding. However, because the channel sands are sub-seismic, their internal geometry and lateral continuity cannot be uniquely resolved from the available seismic data; therefore, the facies architecture shown should not be interpreted as a deterministic description of the subsurface. Although multiple stochastic realisations could in principle be generated, this was not pursued, as dynamic model calibration shows that the upscaled facies-dependent properties, together with the selected aquifer representation, reproduce the observed pressure evolution and injection–withdrawal volumes during both the production and UGS operation periods. This indicates that the model validation and predictive behaviours are not controlled by a specific facies realisation at the scale relevant for storage performance.

4.2. Dynamic Simulation

Table 4 summarises the NRV and corresponding gas in place (GIIP) for the Stockstadt gas field, the Stockstadt UGS, and the Hähnlein UGS as derived from the structural model and compared with the dynamic model initialisation. The Stockstadt gas field volume includes the gas-filled Sand 7 interval within the UGS area, while GIIP for Sand 8 in both UGS sites was not considered during initialisation because these intervals were water-bearing. The comparison between static and dynamic volumes shows a deviation of approximately 10–12%, which lies well within industry-typical volumetric uncertainty ranges of ±10–20% for GIIP estimates based on limited well control and sub-seismic heterogeneity.

This consistency is further supported by the production match: the dynamic simulation yields 63 MMm³ of cumulative production from the Stockstadt UGS compared with 57 MMm³ of actual production (+10%), and 500 MMm³ simulated versus 525 MMm³ actual cumulative production for the Stockstadt gas field. These results demonstrate that, despite the non-uniqueness of the facies realisation, the static–dynamic model linkage is reliable and that the model reproduces volumetrics and production behaviour within an acceptable uncertainty range.

During the production phase, field operation was controlled by constraints on maximum daily gas production, minimum bottom-hole pressure, and maximum allowable water production rate. Under these operational controls, the dynamic model closely reproduces the reported number of producing wells and the overall production duration of the field, consistent with the published LBEG report. This further supports the reliability of the model in representing both volumetric and operational aspects of the historical field development.

The aquifer strength played an important role in the calibration of volumes, water production and pressure development. Figure 8a presents the results of the aquifer sensitivity analysis conducted for the production phase from 1955 to 1978. Simulated gas and water production rates for three aquifer strength scenarios—weak (1,000 VM), medium (10,000 VM), and strong (100,000 VM)—are compared with observed gas production data. Aquifer strength was varied by applying a VM to the aquifer layers below the gas–water contact (GWC), representing different effective aquifer sizes and degrees of hydraulic support. Among the tested cases, the strong aquifer scenario provides the best agreement with the observed production behaviour and was therefore selected for subsequent modelling steps.

Figure 8b shows the simulated reservoir pressure evolution for the three aquifer strength scenarios during the production phase. While all cases exhibit a pressure decline associated with gas production, the magnitude and timing of the decline depend strongly on aquifer strength. In the weak aquifer case (1,000 VM), pressure decreases rapidly and reaches limiting levels at an earlier stage, leading to premature termination of production. The medium aquifer scenario (10,000 VM) delays this pressure depletion but still results in an earlier production stop compared with the strong aquifer case. In contrast, the strong aquifer scenario (100,000 VM) maintains higher pressure levels over a longer period, allowing production to continue for an extended duration. This behaviour is consistent with Figure 8a, where production rates are broadly similar across all aquifer scenarios during active production, but the weaker aquifer cases cease production earlier (red and orange curves). Although direct pressure measurements are not available for validation, the smoother pressure evolution associated with the strong aquifer case is consistent with its improved ability to sustain production and storage volumes in the dynamic simulations, supporting its selection as the base case for subsequent UGS and UHS modelling.

The calibrated model was subsequently applied to the UGS operation period from 2016 to 2025. Figure 8c compares simulated and observed storage capacity, injected volumes, and withdrawn volumes for the Stockstadt and Hähnlein UGS sites. The dynamic simulation was calibrated against injected and produced gas volumes and demonstrates a close match to the observed operational data published by MND Energy Storage Germany GmbH. The improved agreement obtained with the strong aquifer scenario confirms that enhanced hydraulic support from the aquifer is essential for reliably reproducing UGS performance under cyclic injection and withdrawal conditions.

4.3. Hypothetical Uhs Simulation

After calibration of the production and UGS phases, the model was initialised for UHS scenarios. Figure 9 compares the simulation results for the two UHS cases at the Stockstadt UGS over a ten-year cyclic operation. In both cases, injection and production volumes are of comparable magnitude, with injection volumes of approximately 100 MMm³ and production volumes of approximately 80 MMm³ per cycle. However, Case Y (100% H₂) exhibits slightly lower gas production volumes with sharp peaks (Figure 9a) and higher water production rates (Figure 9b) compared with Case X (5% H₂ + 95% CH₄). This behaviour can be attributed to differences in fluid properties of WG and resulting reservoir pressure support. In Case Y, the higher hydrogen content results in lower gas density, high mobility, less stable gas-water displacement behaviour, and reduced volumetric storage efficiency under reservoir conditions. Consequently, the injected gas provides less effective pressure support to counteract water encroachment, allowing increased water production and limiting the recoverable gas volume during withdrawal. In contrast, Case X benefits from stronger effective pressure support, resulting in lower water production and slightly higher gas recovery per cycle as the WG composition is similar to the CG (natural gas).

The compositional response of the produced gas is shown in Figure 9c and Figure 9d. In both cases, the hydrogen mole fraction in the production stream increases rapidly and remains high from the third cycle onward (Figure 9c and Figure 9d), indicating efficient displacement of the initial gas inventory, cyclic behaviour and achieving quasi-steady state. Despite the observed differences in volumetric performance and water production, both scenarios demonstrate the technical feasibility of sustained hydrogen withdrawal with high hydrogen purity after the initial conditioning cycles. The pressure plot (Figure 9e) shows that both systems remain within a similar cyclic pressure envelope because both still rely on the remaining natural gas as CG. However, the pure-H₂ case (Case Y) exhibits more mobile, less damped dynamic behaviour.

4.4. Map-Based Estimation of Storage Potential

Based on hypothetical UHS simulation results for the Stockstadt UGS, a probabilistic, map-based volume estimation was conducted for the other known former gas fields in the area, i.e., the Wolfskehlen, Pfungstadt, Eich, and Darmstadt fields. Field areas were derived from published structural maps (Figure 2), while thickness, porosity, and gas saturation were assigned representative mean values from the Stockstadt–Hähnlein analogue. Uncertainty was propagated via Monte Carlo simulation using a uniform distribution for area (±20%), a triangular distribution for thickness (5–15 m, mode 10 m), and normal distributions for porosity (μ = 0.25, σ = 0.05) and gas saturation (μ = 0.70, σ = 0.0875).

Figure 10 combines input-property variability (box-and-whisker plots) with resulting cumulative distribution functions (CDFs) of WG volume. The boxplots demonstrate that area (Figure 10a) exhibits the largest inter-field variability, whereas thickness (Figure 10b), porosity (Figure 10c), and gas saturation (Figure 10d) show broadly comparable distributions across all fields. Consistently, the CDFs (Figure 10e) indicate that differences in GIIP are primarily controlled by areal extent, with Pfungstadt displaying the highest volumes and widest uncertainty range, and Darmstadt the lowest and most constrained.

The aggregate probabilistic volumes are P90 = 575 × 10⁶ m³, P50 = 910 × 10⁶ m³, and P10 = 1.39 × 10⁹ m³, corresponding to approximately 5.61, 8.89, and 13.62 TWh of CH₄ (1TWh = 102,359,965.88 m³), respectively (equivalent to 1.87, 2.96, and 4.54 TWh of H₂) [23]. The relatively steep CDF slopes for Eich and Wolfskehlen reflect lower uncertainty, whereas the broader spread for Pfungstadt highlights higher structural uncertainty associated with area. However, the CDF slope for Darmstadt is the steepest, possibly because the area distribution has a very narrow spread, which contributes to the steep CDF curve. The Wolfskehlen gas field site can be considered a potential site for UHS conversion or a pilot project, as the data is more complete than at other sites. Wolfskehlen lies within the 3D seismic survey area and has a few wells with wireline logs, and the sand correlates with the Stockstadt and Hähnlein areas. The volume estimates (216 × 10⁶ m³ ≈ 2.11 TWh of CH₄ ≈ 0.70 TWh of H₂) fall within the range of the analogue sites and have comparatively less uncertainty than the other sites.

5. Discussion of Model Uncertainties

Geological uncertainty in the model is primarily related to subsurface interpretation and property representation. Earlier 2D seismic interpretations were refined using the available 3D seismic data, resulting in an improved structural framework and reduced uncertainty in horizon and fault interpretation. However, sub-seismic features such as the internal geometry and lateral continuity of channelised sand bodies cannot be uniquely resolved, and the facies distribution therefore represents a single geologically plausible realisation rather than a deterministic description. Petrophysical property heterogeneity is further simplified by assigning constant porosity and permeability values to the reservoir units, which does not capture small-scale spatial variability but is considered appropriate for the analogue-model objective of this study. Although hydraulic communication between the Stockstadt and Hähnlein reservoirs was reported in the 1980s, this connectivity could not be reproduced in the current model, likely due to limitations in resolving subtle structural or stratigraphic pathways at the available data resolution. These uncertainties are characteristic of legacy reservoirs and are explicitly acknowledged in the static model construction.

Dynamic uncertainty is mainly associated with data availability and scope-related limitations. The dynamic model was validated against publicly available operational data from the UGS period. Within the scope of this study, validation focused on reproducing cumulative production and injection volumes together with the overall field production history. Nevertheless, the available data provide a suitable basis for assessing the model’s field-scale performance. As a result, the model was not calibrated to reproduce short-term operational transients or pressure responses. Despite these limitations, the dynamic simulations provide internally consistent behaviour and reproduce volumetric trends within industry-typical uncertainty ranges. The model is therefore considered fit for purpose as an analogue framework to assess general storage behaviour in porous reservoirs of the northern URG, while recognising that site-specific operational studies would require additional data and further calibration.

Further, the map-based volume estimation approach for the nearby legacy fields implicitly assumes spatial uniformity within each field and does not fully capture intra-reservoir heterogeneity, such as facies variations, compartmentalisation, or localised petrophysical anomalies. The use of independent input distributions does not account for possible correlations between parameters (e.g., porosity–permeability relationships or thickness–facies trends), which may bias the resulting volume distributions. Such heterogeneities may lead to deviations from the modelled distributions, particularly affecting effective porosity, net thickness, and fluid distribution. Consequently, while the probabilistic framework captures first-order uncertainty, the results should be interpreted as screening-level estimates, with the understanding that additional subsurface data (e.g., seismic refinement, well control, and facies modelling) could significantly reduce uncertainty, particularly in defining structural closure and reservoir continuity.

6. Conclusions

This study developed and validated an integrated geological and dynamic model workflow to evaluate underground hydrogen storage in shallow, porous reservoirs in the northern URG. To evaluate this potential, an analogue model representing the Miocene sands of the northern URG was developed, incorporating both Stockstadt and Hähnlein UGS. The integrated geological model successfully reconstructed the structural framework and provided a reliable basis for history matching and UHS prediction.

The validated field-scale model provides a robust analogue for evaluating hydrogen storage behaviour in shallow porous reservoirs while accounting for structural complexity, facies heterogeneity, and aquifer support. The dynamic simulation results of production and UGS phases emphasise that the dimension of the aquifer plays a significant role for these shallow reservoirs, in terms of pressure support and also maintaining storage capacity. The hypothetical UHS scenarios of Stockstadt UGS compare two WG options: i) low-H₂ content (5%, natural gas mixture) and ii) pure-H₂ (100%). The prediction for the case with low-H₂ (5%) content shows stable flow with less water production, similar to natural gas; however, the pure-H₂ scenario exhibited greater gas mobility and higher water production than the low-H₂ mixture. This behaviour is primarily attributed to the lower compressibility of hydrogen under the investigated reservoir conditions, resulting in reduced water displacement and weaker pressure support. Despite these differences, both scenarios achieved WG capacities comparable to those of existing natural gas storage operations. These findings demonstrate that depleted hydrocarbon fields and saline aquifers in similar geological settings can be used for seasonal hydrogen storage with only moderate operational modifications, provided that reservoir-specific hydrodynamic behaviour is properly considered.

This work demonstrates the value of integrating legacy hydrocarbon production data, operational UGS data, modern seismic interpretation, and dynamic compositional simulation within a single workflow. Such an integrated approach reduces uncertainty in UHS feasibility assessments and provides a transferable methodology for evaluating similar porous reservoirs elsewhere. As some of the depleted gas fields, such as Wolfskehlen, have similar geological settings, map-based estimation yields a volume capacity comparable to that of the other two UGS sites. Similar field-scale models for other nearby gas fields or potential prospects can be extracted for hypothetical UHS simulations to test various CGs, geochemical reactions, and methanation. Future work should extend the integrated workflow by coupling the dynamic reservoir model with geomechanical simulations to evaluate well integrity, caprock stability, fault reactivation, surface deformation, and operational pressure limits during hydrogen storage.

The integrated workflow developed in this study provides a transferable framework for assessing underground hydrogen storage in shallow porous reservoirs. Operational experience from the Stockstadt and Hähnlein UGS facilities, combined with validated geological and dynamic models, can be applied to evaluate nearby depleted gas fields and saline aquifers with similar geological characteristics. Future research should expand this workflow by exploring CG optimisation, geochemical and microbial processes, methanation, and integrated geomechanical modelling to more effectively minimise uncertainties in UHS deployment. By enabling the identification of suitable storage sites in the northern URG, this approach supports the development of regional hydrogen infrastructure and contributes to the energy transition and long-term sustainability.