Methane occurrence and origin in Dutch groundwater: from shallow aquifers to deep reservoirs

Abstract

Methane is a common constituent of groundwater with multiple possible origins. Elevated methane concentrations may also result from anthropogenically induced pathways between the deep and shallow subsurface caused by oil and gas production. A baseline characterisation of methane occurrence and origin in the subsurface of the Netherlands was made using a large set of methane concentrations in shallow groundwater (n = 12,219, up to 500 mbgs). Additionally, targeted sampling (n = 40) was carried out in (1) the shallow aquifers at locations where the presence of thermogenic methane was deemed most probable, such as above faults and known gas reservoirs, (2) deep groundwater aquifers below the depth of Neogene and Paleogene marine clays that form the hydrogeological base in the country and (3) geothermal formation waters at 1640–2625 mbgs. Median methane concentrations in shallow aquifers are relatively high from an international perspective (0.2 mg L⁻¹). The highest methane concentrations (up to 120 mg L⁻¹) are attributed to reactive organic matter in Holocene deposits and Pleistocene marine and glacial formations. However, elevated concentrations are also found at greater depth (100–160 m bgs) in Pleistocene aquifers in the eastern and southern inland areas of the Netherlands. Isotopic evidence and gas composition of naturally occurring methane indicate that methane in the targeted samples from shallow aquifers was of biogenic origin, and that methanogenesis predominantly occurs via CO₂ reduction. Only trace amounts of methane (<0.2 mg L⁻¹) were observed in the deep groundwater aquifers. A combination of methane and ethane isotopic composition showed that this methane consists of varying fractions of both biogenic and thermogenic origin. Methane in the geothermal reservoirs has an oil associated thermogenic origin. Overall, these findings highlight that future observations of thermogenic methane in Dutch shallow groundwater (post-Paleogene) are most probably linked to anthropogenically induced connections with the deep subsurface.

Introduction

The subsurface occurrence and chemistry of methane has gained interest in the last decade, since it is a strong greenhouse gas with rapidly growing atmospheric concentrations (Nisbet et al., 2019), its introduction in fresh groundwater can adversely affect water quality (Rice et al., 2018), and it can constitute an explosion and asphyxiation hazard when accumulating in confined spaces (Williams and Aitkenhead, 1991). Notably, the massive scale production of oil and gas from shale formations in the United States, enabled by the application of horizontal drilling and hydraulic fracturing techniques, has raised concerns about the contribution of anthropogenic methane leakage from the deep subsurface to these hazards (Jackson et al., 2014). Leakage of methane through impaired wellbore systems into shallow aquifers and the atmosphere has indeed been shown to locally contribute significantly to greenhouse gas emissions (Kang et al., 2014), lead to undesirable changes in groundwater quality (Forde et al., 2019), and in rare cases cause explosions (Miyazaki, 2009). Furthermore, evidence of stray gas migration originating in the deep subsurface could be a precursor for leakage of more hazardous fluids, such as brines or drilling and production fluids (Vengosh et al., 2014) and its presence signals a conduit for transport that can also hamper the implementation of other technologies relying on the deep subsurface, such as CO₂ storage (Lackey et al., 2019).

As a consequence, groundwater methane concentrations and isotopic composition have increasingly been used to assess the presence of anthropogenically induced connections between shallow aquifers and deeper geological formations. However, methane in groundwater may also occur naturally as the result of either local microbial methanogenesis, leakage of gaseous methane through natural methane seeps and mud volcanoes (Etiope et al., 2009), or mixing with deeper, methane containing, groundwater (Warner et al., 2012). Distinguishing between these multiple possible origins of methane and attributing either a natural or anthropogenic source can be challenging. Therefore, determining the pristine, pre-drilling, methane baselines in groundwater can be vital to avoid later conflicts. Indeed, a number of studies carried out in US states overlying the Marcellus shale play came to conflicting conclusions: while some studies contributed elevated methane concentrations to their proximity to shale gas wells (Jackson et al., 2013; Osborn et al., 2011), others concluded that topographic location played a more important role in determining their distribution (Molofsky et al., 2013, 2011; Siegel et al., 2015).

In the shallow groundwater system of the Netherlands, which largely consists of Quaternary sandy aquifers down to depths of around 400 m, the presence of dissolved methane is ubiquitous and is typically assumed to have a local microbial (i.e. ‘biogenic’) origin (Stuyfzand et al., 1994). A dataset analysed by Fortuin and Willemsen (2005) showed that methane concentrations close to solubility are found near the surface (<20 mbgs) in the part of the country where reactive peat layers occur at shallow depth. The fact that large amounts of methane can be found in the shallow subsurface has been known for centuries; however, as reports of gaseous, flammable, blowouts encountered during drilling of shallow water wells go back as far as the 18^th century. From the end of the 19^th century midway through the 20^th century, gas accumulations that formed below shallow clay deposits were even exploited for private use at farms in the coastal provinces, by extracting groundwater and allowing the degassing methane to accumulate below a suspended kettle balanced by counterweights (Bol, 1991).

Methane isotopes of gas collected from two such systems indeed revealed a biogenic (microbial) origin (Buijs and Stuurman, 2003). More recently, methane in nine groundwater wells with known high methane concentrations was found to be of biogenic origin as well, with a sample collected near the site of a catastrophic gas well blowout that occurred in 1965 in the eastern Netherlands being the exception (Cirkel et al., 2015). These findings confirm that methane in shallow groundwater of the Netherlands is predominantly of biogenic origin. However, thermogenic gas locally migrating from depth through either natural or anthropogenic conduits cannot be ruled out, given the presence of numerous oil and gas reservoirs in the deep subsurface (Fig. 1). Since production of these resources started in the 1940s, around 2500 oil and gas wells have been drilled, of which more than 1700 onshore (NLOG, 2020). The majority of produced gas comes from the Groningen gas field, which is by far the country’s largest natural gas accumulation and one of the biggest in the world. However, reservoir depletion here has led to induced seismicity with significant societal impact (van Thienen-Visser and Breunese, 2015). As a result, production from this reservoir has recently been entirely halted and drilling of new onshore wells from other reservoirs has also decreased dramatically. Besides these conventional reservoirs, two extensive shale formations are present (Fig. 1) that are potentially economically viable (Bergen et al., 2013). However, a decree by the Dutch government has banned both shale gas exploration and production due to environmental concerns (MIE and MEAC, 2018).

Figure 1. Map of the Netherlands showing the onshore oil and gas fields, major shale formations, maximum extent of Holocene Marine deposits (after Griffioen et al., 2013), averaged methane concentrations per well in the datamined dataset (A), locations of the samples collected for isotopic analysis (A + B), depth to base of the Neogene Breda formation and the main faults penetrating through to the surficial geology (B). Samples labelled SGW and DGW are from shallow and deep groundwater wells, respectively, indicating whether the sampled well screens were above or below the base of the Breda Formation. Samples labelled GTW were collected from geothermal wells. Geological data retrieved from the Netherlands Oil and Gas Portal (NLOG.nl) and the data and information site of the Geological Survey of the Netherlands (www.DINOloket.nl). Area west and north of the line showing the maximum extent of Holocene marine deposits is referred to as the coastal lowlands.

In spite of the lack of production from shale plays, hydraulic fracturing has been applied to improve production in a number of ‘conventional’ fields. A survey commissioned by the State Supervision of Mines (SodM) tallied a total of 245 frack jobs dating back to 1954, 63 of which did not involve high pressure fluid injection but only an acid treatment (TNO, 2018). Given the depth (900–4000 m below mean sea level) of these frack jobs, the thickness of overlying cap rocks and the stress regime, the risk of leakage to shallow groundwater occurring through the induced fractures is low (TNO, 2018). However, migration of methane to shallow groundwater through impaired wellbores has been shown to occur. A recent survey by the SodM found some form of well barrier failure at 227 (23%) out of 986 gas wells (SodM, 2019). Failure of a single barrier does not necessarily imply that leakage to the environment occurs, since there are typically multiple barriers in place. Nonetheless, observations of thermogenic gas bubbles in flooded well cellars showed that gas migration originating in the deep subsurface occurred at at least 13 of the surveyed gas wells (1.3%). Well integrity failure has also been confirmed by means of gas flux measurements in the vadose zone above a cut and buried, fully decommissioned gas well (Schout et al., 2019). Lastly, migration of thermogenic gas through the fractures opened by the aforementioned blowout site is believed to be ongoing (Schout et al., 2017).

Contrary to these anthropogenic leakages, and in spite of the presence of a number of significant faults and fractures in the country (Fig. 1), no clear evidence exists of gas migration from the deep subsurface to the shallow groundwater environment through natural conduits. Fluid flow towards the shallow aquifers is assumed to be limited by a thick sequence of impermeable Neogene and Paleogene marine clays that directly underlie the aquifer system in nearly the entire country (Griffioen et al., 2016b). There is some evidence that local conduits for fluid flow through these barriers exist on the Dutch continental shelf in the North Sea, as methane with an apparent mixed biogenic and thermogenic origin was observed near an offshore pockmark (Schroot et al., 2005). However, the mixed isotopic signature could also be explained by oxidative fractionation of biogenic methane, and hence the evidence was inconclusive. Other observations of shallow gas encountered in the Dutch part of the North Sea were determined to be of microbial origin (Verweij et al., 2018), with the exception of one well in the northern offshore where thermogenic methane has apparently migrated in to the Quaternary deposits (Ten Veen et al., 2013). Onshore, migration of thermogenic methane through a fault to the atmosphere was observed in a former coal mining area in Germany, around 75 km east of the border with the Netherlands. However, neither the Paleogene marine clays nor other Cenozoic unconsolidated sediments that characterise the Dutch subsurface are present at this location. Hence, gas migration was concentrated in faults penetrating the outcropping Cretaceous consolidated rocks (Thielemann et al., 2000).

Given the abundant presence of methane in shallow groundwater, the long history of oil and gas extraction, and potential for future use of the deep subsurface, there is a clear need to further develop the groundwater methane baseline of the Netherlands. The limited number of observations suggest that shallow methane is purely of microbiological origin, and hence that findings of thermogenic methane in this groundwater compartment would indicate an anthropogenic disturbance. A thorough analysis on the controls on methane concentrations in the country is lacking. Importantly, the occurrence and origin of methane between the shallow aquifer system and the depth of the oil and gas reservoirs is unknown. In this paper, the distribution of methane concentrations in shallow groundwater was assessed using a methane concentration dataset compiled through data-mining and consisting of more than 12,000 individual samples, partly consisting of repeat samples over time in the same well screens. Furthermore, a targeted sampling campaign was conducted, analysing both groundwater below the marine clay barriers underlying the shallow aquifers and above these barriers at locations where thermogenic gas migration was deemed most likely due to presence of faults or fractures. The dissolved gas compositions were then compared to those from oil and gas wells in the country, such that the molecular and isotopic composition of (dissolved) gasses throughout the entire depth interval is analysed.

Hydrogeology of the Netherlands

Quaternary and Neogene shallow groundwater system

The Netherlands is a largely deltaic country where several major European rivers coalesce (Rhine, Meuse and Scheldt). It has a total land surface area (including inland waters) of 37,363 km². It is located at the south-east of the subsiding North Sea sedimentary basin and neighbours Belgium and Germany on its southern and eastern border, respectively (Fig. 1). The climate is temperate with an annual precipitation of around 850 mm. The hydrogeology is determined by the presence of thick successions of Cenozoic unconsolidated sediments of dominantly marine and fluvial origin that were deposited during various stages of sea level highs and lows. At the maximum extent of the Saalian Glaciation, the northern half of the country became covered by ice. Glacial till formations were left behind in the northern half of the country and ice pushed ridges are present in the centre and east of the country that have a maximum elevation of 110 m above sea level (masl). Sand dunes are situated along the coastline that play an important role in flood protection and reach a maximum of around 50 masl.

During the Holocene, encroaching sea levels caused peat to form at a regional scale behind the coastal barriers. The peat area then became inundated, and clayey marine sediments were deposited in the resulting tidal lagoons. This part of the country is referred to as the coastal lowlands (Fig. 1). A significant part of the Holocene coastal area is now formed by polders that are situated below sea level and are drained by a system of ditches and canals. This causes groundwater exfiltration of sometimes saline groundwater in these polders. Another part of the Holocene landscape comprises of riverine deposits. This Holocene riverine area has not been inundated by seawater. The Holocene peat and clay act as a confining layer to the underlying Pleistocene aquifers in both Holocene areas, while regional recharge occurs from the coastal dunes and ice-pushed ridges. In most of the other parts of the country, the surficial geology is defined by Pleistocene fluvial sediments that form phreatic aquifers. The surficial Pleistocene sands in the Centre, East and South constitute the major groundwater recharge areas of the country (Meinardi, 1994). The most south-eastern part of the country has an entirely different geology with Cretaceous formations at shallow depth. The groundwater table is generally very shallow, with the exception of elevated areas in the south-eastern parts of the country and the ice pushed ridges in the centre and East. In these locations, the depth to groundwater can be several 10s of meters.

The units that host active circulation of fresh groundwater are the more coarse-grained Quaternary and Neogene deposits that increase in thickness from roughly 50 m in the South-East to around 400–500 m in the North-West (Dufour, 1998) and are part of the Upper North Sea Group (Table 1). Here, turnover time of groundwater is in the order of days to locally more than 10,000 years (de Vries, 2007). Full stratigraphic details of recognised Cenozoic formations in the Netherlands are given in Fig. S1. Notable aquifers are found in the Pleistocene fluvial deposits of the Kreftenheye, Urk, Sterksel, Waalre, Beegden, Appelscha, Peize and Kieseloolite Formations as well as the sandy units of the marine Early Pleistocene Maassluis Formation and the Neogene Oosterhout and Breda Formations. However, clay layers of regional significance may be present in these units. Low permeability glacial clays of the Peelo (Elsterian) and Drente (Saalian) Formations form barriers to flow in the North-East of the country. Marine clays were also deposited in glacial basins formed during the Saalian ice-age that are assigned to the Eem Formation (de Gans et al., 2000).

Table 1. Simplified stratigraphic overview of the Netherlands showing the geological units relevant to this paper. The dots in columns R, C and S show whether the unit typically serves as a reservoir/aquifer (R), caprock/confining layer (C) or source rock (S). Based on Jager and Geluk, 2007. Detailed lithostratigraphy of Cenozoic formations is shown in Fig. S1

The Early Pleistocene and Neogene marine formations that form the base of the Upper North Sea Group, increasingly consist of thick clay deposits with depth. In the coastal lowlands, they contain saline groundwater that has also intruded parts of overlying fluvial aquifers (de Vries, 2007), and the fresh/brackish interface varies in depth between 50 and 130 mbgs (Mendizabal et al., 2011). In the Pleistocene part of the country, the fresh–salt interface is generally deeper (between 100 and 500 mbgs). However, towards the eastern border it follows the depth of the underlying Paleogene marine formations of the Middle and Lower North Sea Groups (Table 1) which can be just 10 m below the surface. The depth and thickness of Neogene formations, as well as the depth of the fresh–salt interface, are much greater in the Roer Valley Graben in the South of the Netherlands, which is bounded by South-East North-West running faults on either side (Fig. 1). Here, rapid deposition of these formations could occur in the Cenozoic due to strong subsidence associated with rifting. Steps in hydraulic head over these faults show that they form strong barriers to horizontal groundwater flow in the shallow groundwater systems (Bense et al., 2003). Rifting during earlier geological eras associated with the same rift system accounts for the extension of these faults into the west of the Netherlands, where they may also penetrate though to the base of the Neogene Formations (Verweij et al., 2012).

Paleogene and older groundwater and hydrocarbon systems

The Paleogene formations underlying the shallow groundwater system are found in virtually the entire country and dominantly consist of low permeable clays that can be up to several 100 m thick (Verweij et al., 2012). Notable deposits include the Rupel Clay and Ieper Clay Members. The former is known as the Boom Clay in Belgium where it is considered as a host rock for the permanent underground disposal of nuclear waste (Verhoef et al., 2014). However, more permeable units also exist in the Middle and Lower North Sea Groups, such as the Brussel Sand and the Basal Dongen Sand Members. Groundwater in these formations is generally saline or hypersaline, except at the borders where outcrops of these formations may lie relatively nearby (Griffioen et al., 2016b). Locally, gas accumulations have been encountered in the deepest permeable Paleogene units. These have remained mostly undeveloped due to anticipated subsidence problems (Jager and Geluk, 2007). The origin of these shallow gas accumulations is subject to debate, and molecular and isotopic evidence is limited. The three available gas samples from such wells from onshore Netherlands (data from NLOG, 2020) are also inconclusive as isotopic analyses were not carried out. The C₁/(C₂ + C₃) molar ratio measured for these samples range from 422 to 1930 (Table S1), which appears to indicate a mixed biogenic – thermogenic origin.

An analysis of fluid pressures in the Dutch subsurface has shown that overpressurised conditions, resulting from the balance between sediment loading and pressure dissipation by fluid flow, are primarily observed in the Northern Netherlands and northern offshore. Overpressurised conditions are not encountered in the Southern Netherlands, which is most likely linked to aforementioned faulting in this area. Together with a gradual increase in the deposition of sediments with larger grain sizes towards the margins of the North Sea Basin (located near the Southern border of the Netherlands), these faults have likely contributed to the dissipation of overpressures by vertical fluid flow on geological time scales (Verweij et al., 2012). Prominent geological units below which overpressures have been observed are the claystones of the Lower Cretaceous Rijnland Group and the salt deposits of the Triassic Röt Formation and the Permian Zechstein Group (Table 1).

These two salt deposits also form the most important caprocks for the Paleozoic gas system that is by far the largest in the Netherlands by volume. The major source rock for this system comes from Carboniferous coals. In the south-eastern tip of the country, these coals are located much closer to the surface, where they were mined commercially up to the 1970s. Noteworthy reservoir rocks for this system are the Permian Slochteren Formation (which hosts the Groningen gas field) and the Triassic Main Buntsandstein Subgroup (Table 1). The lower Cretaceous claystones form the caprocks for a Mesozoic mixed oil and gas system, where the Jurassic Posidonia Shale Formation serves as the major source rock. Reservoirs are typically found in sandstones belonging to the Schieland and Rijnland Groups (Table 1). Salt diapirs formed by movement of the underlying Zechstein salt deposits have led to faulting (ten Veen et al., 2012) resulting in potential pathways for upward fluid movement. However, these diapirs as well as the associated faults have also been shown to contribute to trapping of hydrocarbons in these formations (Jager and Geluk, 2007). Lastly, the two shale formations in the Netherlands that have been identified to hold potentially economically viable gas contents are the aforementioned Posidonia Shale Formation and the Carboniferous Geverik Member (Bergen et al., 2013).

Methods

Dataset 1: Nationwide methane concentrations

Data sources and data processing

A dataset of methane concentrations in Dutch groundwater was compiled by datamining. The majority of data comes from two main sources: (1) monitoring data from drinking water production companies, (2) monitoring data from thermal energy storage system operators (Kappelhof et al., 2006) and (3) all methane concentration data available in the online database of TNO Geological Survey of the Netherlands (www.dinoloket.nl). Samples in the dataset were collected in the period 1981 to 2013, with the exception of three samples that were collected decades earlier. Sampling methodologies used for individual samples are not traceable; however, the most commonly applied methodology is likely the inverted bucket method, given the time period. More advanced sampling methods were also used; however, including sampling under high pressure in stainless steel vessels using submersible pumps. The inverted bucket methods have been shown to be reliable, at least up to concentrations where effervescence would occur at atmospheric conditions (Molofsky et al., 2016). Given the variety of sources and methods, the minimum reporting limit (MRL) also varied but an MRL of 0.01 mg L⁻¹ was most frequently encountered. To homogenise the dataset concentrations reported as zero or <0.01 mg L⁻¹ were assigned values of 0.01 mg L⁻¹.

This dataset consists of 12,219 datapoints. For quality control considerations, datapoints lacking information regarding either the sampled depth interval or the coordinates were excluded from the further analysis at this point. The samples were collected from 2,063 wells, including 426 multi-level monitoring wells that were sampled at at least 2 depths. The number of unique well screens in the dataset is 3,250. The median, average and maximum number of samples collected from each well screen are 1, 3.6 and 31, respectively. Calculated coefficients of variation (CV) for wells screens with at least 2 samples are less than 0.2 in 52% of measurements. Only 3% of well screens have a CV larger than 1 (Fig. 2). Out of only 4 well screens with a CV greater than 1 and an average methane concentration above 1 mg L⁻¹, the maximum methane concentration was just 1.47 mg L⁻¹. This shows that methane concentrations generally have low levels of natural temporal variability (or limitations in sampling) and hence averaged methane concentration values were used throughout the paper, unless otherwise noted. For each of these well screens the host formation, depositional environment and geological age were determined, based on the depth of the middle of the well screen and extracted from the Digital Geological Model of the Netherlands (DGM; www.dinoloket.nl).

Figure 2. Relative frequency histogram of the coefficient of variation (CV) of methane concentrations for well screens that were sampled more than once. Only well screens with average methane concentrations >0.01 mg L⁻¹ were considered.

Statistical analysis

Methane concentration data were grouped to test whether there were significant differences between subgroups of data based on: the stratigraphic unit at the well screen depth, the sedimentary depositional environment of the formation at the well screen depth, proximity to gas wells, proximity to faults penetrating the surficial Upper North Sea Subgroup, the presence of oil and gas reservoirs at the well location, the presence of the Posidonia Shale, the presence of the Geverik Shale, the presence of the Nieuwkoop Formation and the presence of the Eem Formation. The presence of Holocene marine sediments was also used, and the two categories formed by this subdivision will be referred to as the coastal lowlands versus the inland area. Threshold distances of 1 and 2.5 km were used to test whether proximity to oil and gas wells or faults influenced methane concentrations. The geospatial data used to make these groupings is shown in Fig. S2. Given that distributions of methane concentrations tend to be positively skewed, the non-parametric Mann–Whitney U test was used to test for differences in median methane concentrations between two subgroups and the Kruskal–Wallis test for categories with more than two independent samples, as is typically done in statistical assessments of dissolved methane in groundwater (e.g. Atkins et al., 2015; Christian et al., 2016; Humez et al., 2016; Jackson et al., 2013). Kruskal–Wallis tests were extended by a stepwise stepdown multiple comparison post hoc procedure to determine statistical differences between the subsets. Statistical analyses of the methane concentration dataset were carried out using SPSS v25 and evaluated at p-values of 0.01.

Dataset 2: Targeted groundwater samples

Sampling locations

An additional 40 water samples were collected for more detailed isotopic analysis of dissolved gasses. The aim was to cover the entire depth interval, from shallow groundwater aquifers down to formation waters at reservoir depth. In the following, we consider shallow groundwater to be groundwater sourced from Neogene and younger formations – those where active groundwater circulation takes place in most of the country. Shallow groundwater sampling locations were selected based on likelihood of encountering thermogenic methane or because anomalously high methane concentrations had been observed before. Likelihood of encountering thermogenic methane was deemed high when water wells were located in close proximity to a fault penetrating through to the surficial Upper North Sea Subgroup, a salt diapir and/or were situated above a known gas reservoir. In total, water samples were collected from 17 such groundwater wells (Fig. 1). Two of these 17 were multi-level monitoring wells that were sampled at four different depths each (SGW-02 and SGW-16). The 15 other wells were only sampled at one depth. The depth of the middle of the well screen for wells in this category ranged from 5 to 280 mbgs.

A further 11 samples were collected from deep groundwater wells (Fig. 1), which in this context we consider to be below the Paleogene Rupel Clay Member. The relatively small areas of the country where this formation is not present where therefore excluded. Local geology was determined using DGM. Given the large depth of the Rupel Clay Member in most of the country, the salinity of the water at these depths, and the typically poor productivity of Paleogene formations, no more wells could be sampled. Also, the 11 sampled wells are all situated near the southern and eastern border of the country, given that the depth of the Rupel Formation increases rapidly towards the north-west of the country (following the general pattern of the Breda Formation, Fig. 1B). The depth of the middle of the well screens for the sampled wells in this category ranges from 130 to 871 mbgs (Table 3). The majority of sampled deep and shallow groundwater wells were collected from dedicated groundwater monitoring wells with short well screen lengths (<3 m). A subset of samples was collected from groundwater production wells with well screens ranging from 10 to 30 m.

Besides groundwater wells, molecular and isotopic compositions of dissolved gasses were also analysed for 6 samples retrieved from geothermal wells (Fig. 1), that range in depth from 1640 to 2625 mbgs. Results of these analyses, as well as that of 9 shallow groundwater wells, have been reported before in Dutch literature (Cirkel et al., 2015; Table 3). For comparison, isotopic compositions of natural gas from Dutch oil and gas wells were retrieved from the governmental ‘Netherlands Oil and Gas Portal’ (NLOG, 2020). Lastly, the water composition of 10 samples from oil field formation waters (Eggenkamp, 1994) was re-analysed. These samples were retrieved from oil wells in the west of the Netherlands in 1983 and preserved in glass bottles with a thick oil layer on top of the water. These wells range in depth from 950 to 1766 mbgs. However, due to the age of these samples and lack of preservation measures, they were not re-analyzed for dissolved gas isotopic composition.

Sampling and analysis

Wells were purged prior to sampling until around three times the total volume of water in the well tube had been flushed or after the temperature, pH and electrical conductivity of the water had stabilised and dissolved oxygen indicated stable anoxic conditions (field measurements carried out using WTW probes). Different pumps were used depending on the casing diameter and groundwater level. Submersible pumps were used for wells with casing diameters >2” and peristaltic pumps for diameters <2”. Inertial pumps had to be used in two cases (DGW-08 and DGW-09), as the casing diameters were <2” and groundwater heads were too deep to use a peristaltic pump. Samples were filtered through 0.45 μm pore filters and stored at 5°C until analysis. Major anions and cations were determined using ion chromatography (IC) and inductively coupled plasma optical emission spectrometry (ICP-OES), respectively. IC measurements were carried out on a Dionex ICS-1500 equipped with an IONPAC AS14 anion exchange column and an A SRSs-Ultra 14 mm suppressor (Dionex Corporation, Sunnyvale CA, USA). ICP-OES measurements were carried out on a a Perkin-Elmer Avio 500. Samples for ICP-OES analysis were acidified using HNO₃. Alkalinity was determined using titration with sulphuric acid and ammonium concentrations using a photo spectrometer.

Samples for the analysis of molecular and isotopic composition of dissolved gasses were collected in Isoflasks® that were connected directly to the pump tubing to avoid loss of gas to the atmosphere. However, this was not possible for the two samples collected using an inertial pump; here, water was first pumped into a bucket and subsequently pumped into the Isoflasks® using a peristaltic pump. Analysis of dissolved gasses was carried out by ISOLAB, the Netherlands. Samples were allowed to equilibrate to the room temperature and air pressure of the laboratory. Subsequently, the total headspace gas was extracted using a syringe. The total water volume was determined by weighing the samples. Meanwhile, the pressure and temperature in the laboratory were recorded with which the mass concentration of gas components in the aqueous phase of the original sample were back calculated using Henry’s law. Up to three Isoflasks were used per well to ensure sufficient gas was available to analyse at least the gas composition and δ¹³C-CH₄ ratio reliably. The gas composition was analysed on Agilent 6890N/7890A/7890B Gas Chromatographs (GC). Carbon isotopes of methane were analysed with an Agilent 6890N GC interfaced to a Finigan Delta S isotope-ratio mass spectrometer (IRMS). Carbon isotopes of C_2, C₃ and CO₂, as well as hydrogen isotopes of C₁ were analysed on an Agilent 7890A GC interfaced to a MAT 253 IRMS. Isotopic compositions of carbon and hydrogen are reported relative to PDB and SMOW standards, respectively, using the standard delta notations (e.g. Whiticar, 1999).

Chloride stable isotope ratio (³⁷Cl/³⁵Cl) was determined by mass spectrometry on chloromethane (CH₃Cl) following the procedures in Eggenkamp (1994) and Godon et al. (2004) at the Institut de Physique du Globe de Paris (IPGP). δ³⁷Cl is reported relative to the Standard Mean Ocean Chloride standard (SMOC; Kaufmann, 1984). δ³⁷Cl was only determined for the deep groundwater samples, a subset of the shallow groundwater samples, and the formation waters. For the formation waters, obtained values were in good agreement with reported values in Eggenkamp (1994), who also analysed δ³⁷Cl values for 8 out of the 10 available samples. Only 1 of the 10 samples showed a large increase in chloride concentration over this time due to evaporation caused by improper sealing of the bottle. Nevertheless, the difference in δ³⁷Cl was only 0.2‰. The average difference between reported δ³⁷Cl in 1994 and now was only 0.08‰, which shows that the consistency of the measurement technique is high. δ¹³C of dissolved inorganic carbon (DIC) and subsequently ¹⁴C activity were determined by the Centre for Isotope Research (CIO), University of Groningen, on an Accelerator Mass Spectrometer (AMS). Samples for ¹⁴C determination were collected in dark glass bottles and sterilised using an I₂-KI solution. ¹⁴C activity is reported as percentage modern carbon (pMC). Conversion of ¹⁴C activity to pMC was done according to definitions in Mook and Plicht (1999). ¹⁴C was only determined for the deep groundwater samples and a subset of the shallow groundwater samples.

Results and discussion

Nationwide groundwater methane concentration dataset

Median and average and maximum methane concentrations in the 3,250 well screens in the dataset were 0.20 mg L⁻¹, 3.00 mg L⁻¹ and 96.50 mg L⁻¹, respectively. This shows that methane concentrations are positively skewed, as is typically observed in dissolved methane baseline studies. Concentrations are greater than 0.01 mg L⁻¹ in 79% of the well screens. The 10 mg L⁻¹ and 28 mg L⁻¹ hazard mitigation threshold values as proposed by the US Department of the Interior (Eltschlager et al., 2001) are exceeded in 7.5% and 2.8% of well screens. Median methane concentrations in the Netherlands are considerably higher than those observed in the neighbouring German state of Lower Saxony (∼0.002 mg L⁻¹, Schloemer et al., 2016), which borders the Netherlands in the north-east, and in Great Britain, where the largest median value was observed for South Wales (∼0.032 mg L⁻¹, Bell et al., 2017). While dissolved methane concentrations in the Netherlands are thus relatively high, much higher median methane concentrations have been observed at the regional scale, such as in the Los Angeles Basin (39 mg L⁻¹, Kulongoski et al., 2018).

Methane concentrations and depth

Methane concentrations varied significantly with depth within sampled multilevel wells (Fig. S3). Out of 368 groundwater wells with at least two sampled well screens with methane concentrations >0.01 mg L⁻¹, the coefficient of variation, calculated as the standard deviation of the methane concentrations in the well screens present in the well divided by their total mean concentration, was greater than 1.0 in 20% of wells. For this group, 33 of the 76 wells had a mean methane concentration above 1 mg L⁻¹, showing that such variability does not only appear at the low end of the methane concentration range. This vertical spatial variability is much greater than the variability of methane concentrations within single well screens over time (Fig. 2).

Overall, median methane concentrations per depth interval increase with depth down to 100 - 140 mbgs, where they exceed 1 mg L⁻¹. Below 140 mbgs median concentrations gradually drop off again to 0.02 mg L⁻¹ for depths greater than 240 mbgs. Maxima vary between 40 and 100 mg L^-1 down to 140 mbgs, after which they drop off to around 10 and 20 mg L⁻¹. Only the maximum reported single measurement (before averaging concentrations over time in each well screen) of 120 mg L⁻¹ exceeded the theoretical solubility corresponding to the depth of the middle of the well screen and a 100% methane mole fraction (∼96 mg L⁻¹ at 22.3 mbgs). The pattern of decreasing methane concentrations with depth suggests the high concentrations result from local biogenic methanogenesis rather than diffusion from underlying sources of methane, which would yield an overall pattern of increasing methane concentrations with depth.

When isolating the wells that are in the coastal lowlands (n = 369, Fig. 1A), methane concentrations show a distinctly different pattern. Median concentrations are greatest at shallow depths and peak between 40 and 60 mbgs at 4.0 mg L⁻¹, compared to just 0.2 mg L⁻¹ in the same depth interval in the inland areas (Fig. 3). Furthermore, the variability in methane concentrations down to around 100 mbgs is more significant in this area, as shown by the greater interquartile range. However, from 100 to 160 mbgs, both median and 3^rd quartile concentrations are actually smaller in the Holocene part of the country. At this depth interval, third quartile concentrations in the inland areas are around 4 mg L⁻¹. The lower methane concentrations at these depths in the coastal lowlands versus the inland areas can possibly be explained by the generally shallower depth of the fresh-salt interface in this region (Mendizabal et al., 2011) and hence the prevalence of sulphate, which prevents methanogenesis (Molofsky et al., 2018). These results give a more nuanced picture of methane concentrations in the country than earlier work by Fortuin and Willemsen (2005), who concluded that the highest concentrations occur at shallow depth in the presence of reactive organic matter such as peat. While this is indeed true for the coastal areas, there are other sources of methane in slightly deeper groundwater aquifers in the inland areas, resulting in overall higher methane concentrations there.

Figure 3. Medians and interquartile range of methane concentrations per depth interval of twenty meters. Results are shown for wells in the part of the country with Holocene marine sediments (blue) and other areas (green), as per the red line in Fig. 1A.

Depositional environment and stratigraphic unit

Methane concentrations in water samples from glacial deposits (1.7 mg L⁻¹) are significantly higher than in sediments with other depositional origins (Fig. 4A). This finding appears consistent with observations in US groundwater, where methane concentrations are also higher in unconsolidated glacial sediments than in unconsolidated non-glacial sediments (McMahon et al., 2017). However, median concentrations in glacial sediments in the United States (∼0.001 mg L⁻¹) are much lower overall than in the Netherlands. Concentrations in formations of marine and fluvial origin are lower and not significantly different from one another, with median concentrations of 0.21 and 0.31 mg L⁻¹, respectively. Median concentrations in formations with other depositional origins were again statistically significantly lower with a median of just 0.02 mg L⁻¹. Most of the samples in this category are from formations with an aeolian origin, or from ice pushed ridges. These generally very permeable sandy formations serve as important infiltration areas for groundwater in the Netherlands. Hence, the observed low methane concentrations are not surprising, given the predominant Fe-anoxic redox state of the groundwater (Griffioen et al., 2013).

Figure 4. Box and whisker plots of methane concentrations versus depositional environment (A) and versus Neogene and Quaternary geological units (B), shown in order of increasing age (as per Fig. S1). Concentrations in Holocene (not including the Boxtel Formation) and Paleogene formations are aggregated in plot B and excluded from A. Colours in B indicate the depositional environment as per plot A (Fluvial = Green, Marine = Orange, Glacial = Purple, White = Other). Boxes represent values within the 25^th and 75^th percentile, whiskers the minimum and maximum concentrations. Number labels show the total number of samples per category, capital letters indicate homogeneous subsets following the multiple comparison test procedure.

Multiple comparison tests on the stratigraphic units in which wells were screened revealed a more complex image, with significant differences between 9 homogeneous subsets that are much larger than the differences between the depositional environments (Fig. 4B). The highest median concentrations were observed in the glacial Peelo Formation (4.85 mg L⁻¹), much larger than those in the Drente Formation (0.14 mg L⁻¹), which is the only other glacial formation in the Netherlands. Organic matter content in the Peelo Formation is indeed known to be higher than in the Drente Formation, particularly in the clay deposits that can be more than one hundred meter thick where they have partially infilled deeply incised subglacial tunnels that are not present in the Drente Formation (Griffioen et al., 2016a). The data show that local methanogenesis in these clays is likely a source of dissolved methane in the Netherlands.

Median methane concentrations in marine deposits are highest in the Late Pleistocene Eem Formation (1.50 mg L⁻¹). Third quartile concentrations in the Eem Formation are particularly high at 49.6 mg L⁻¹, indicating that significant methanogenesis likely occurs in this formation. This corresponds to findings of shallow gas in the Danish subsurface, which were contributed to significant methanogenesis in a comparable Late Pleistocene marine formation formed in glacial basins (Laier et al., 1992). Median concentrations in Neogene and Early Pleistocene marine deposits were considerably lower (<0.22 mg L⁻¹). Only 5 of 38 total measurements in Paleogene marine formations were above the MRL, and only 1 is larger than 1 mg L⁻¹. This finding also follows from the low observed methane concentrations at depths greater than roughly 200 mbgs (Fig. 3), where the Paleogene formations typically occur in the Netherlands. Similarly, methanogenesis in these formations does not have to be carbon limited but is more likely the result of unfavourable redox conditions due to the presence of SO₄ rich saline groundwater.

The occurrence of methane in fluvial deposits also varies greatly, with median concentrations in the Appelscha Formation (4.65 mg L⁻¹) and to a lesser extent the combined Peize and Waalre Formations (1.10 mg L⁻¹) being significantly larger than in other fluvial deposits. A local methanogenic source is less likely in these formations, given that they are characterised by medium to coarse grain sands and gravel (Griffioen et al., 2016a) that are not typically associated with high organic matter contents. Since the Appelscha Formation overlies the Peize-Waalre Formation, the occurrence of high concentrations of dissolved methane here is more likely attributable to varying overlying sources, such as the Peelo clays, the Marine Eem Formation or Holocene deposits. One may note that the ice-pushed ridges, which are frequently composed of fluvial sediments from the Urk Formation, have very low methane concentrations, with third quartile concentrations just slightly higher than 0.01 mg L⁻¹.

Surprisingly, the median methane concentration in Holocene deposits is only 0.08 mg L⁻¹. However, the interquartile range is relatively large. The high methane concentration observed at shallow depth in the Holocene marine area (Fig. 3) are therefore likely very local phenomena, linked to reactive carbon in peat of the Nieuwkoop Formation and/or the marine Holocene Naaldwijk Formation.

Spatial position

The spatial distribution of methane concentrations was tested by comparing the average and median concentration in a subset of wells that overlap with certain spatial criteria (Fig. S2) to those that don’t. The biggest difference in average methane concentrations was observed for the subset of wells overlapping with the Holocene Nieuwkoop Formation (10.1 mg L⁻¹), which is the major peat formation in the Netherlands. The difference in median and average methane concentrations for wells in the area with Holocene marine deposits (Fig. 1A) was much smaller (0.5 and 4.7 mg L⁻¹, respectively). This suggests that the Holocene marine deposits might play a smaller role in determining methane concentrations than the presence of peat does. Surprisingly, methane concentrations in wells that intersect the Eem Formation were not significantly elevated (Table 2), in spite of the high observed concentrations in samples directly from the Eem Formation. The discrepancy results from the fact that while 283 samples are from wells that intersect the Eem Formation, only 30 water samples were directly obtained from it (Fig. 4). This highlights a key limitation of using statistical analysis of horizontal spatial relations for determining sources of groundwater methane, which in reality are determined by depth and time dependent processes.

Table 2. Difference in mean and median values for subgroups of well screens using a criterion with respect to a given threshold distance of or intersecting a certain spatial feature (shown in Fig. S2) that is a potential source of methane in shallow groundwater. P-values are for associated Mann–Whitney U tests

Differences in median concentrations were also high for wells that overlap with known oil and gas reservoirs (7.3 mg L⁻¹). To a lesser extent, median concentrations in groundwater wells in close proximity to gas wells were also substantially higher than those further away, for both cut off distances of 1 (2.9 mg L⁻¹) and 2.5 km (2.7 mg L⁻¹), respectively. Whether this is evidence of a deep, thermogenic source at play in controlling the methane concentrations is questionable however. There is strong overlap between the occurrence of peat and that of oil and gas accumulations in the Netherlands (Fig. S2), which cannot be distinguished between here. Although a dominantly biogenic origin of methane remains the most likely factor controlling methane concentrations in the Netherlands, a local thermogenic source near gas wells or above oil and gas fields in general cannot be ruled out on the basis of methane concentrations alone. Hence, isotopic evidence is required to determine whether a substantial thermogenic component contributes to the elevated concentrations above oil and gas fields and near gas wells. The spatial analysis did reveal that other potential deep sources can be ruled out as major sources of methane in shallow groundwater, as concentrations in groundwater wells near faults were not significantly different, and methane concentrations above the Posidonia and Geverik shales differed less than 1 mg L⁻¹ with those that do not overly these formations.

Methane origin in targeted groundwater samples

Methane concentrations in the targeted samples from shallow groundwater ranged up to 79 mg L⁻¹, compared to a maximum methane concentration of just 0.18 mg L⁻¹ in the deep groundwater samples (Fig. 5). Even the smallest concentrations were still considerably higher than methane concentrations of groundwater in contact with air at 10°C (∼0.00006 mg L⁻¹). Methane concentrations above 3 mg.L⁻¹ only occurred in samples with sulphate concentrations below 3 mg L⁻¹ and vice versa (Table 3). Two samples had higher co-occurring concentrations, but these were taken from wells with long well screen lengths, which can allow for mixing of different water types. The low methane concentrations in the deeper groundwater samples are therefore likely attributable to redox conditions that are unfavourable for microbial methanogenesis, as sulphate concentrations ranged from 4 to 1591 mg L⁻¹ in these waters.

Figure 5. Methane concentrations versus carbon isotopic composition (δ¹³C-CH₄) in targeted samples from deep and shallow groundwater_..

Table 3. Molecular and isotopic composition of dissolved gases, ¹⁴C-DIC, chloride stable isotope ratio and chloride and sulphate concentrations for samples that were(1) collected and newly reported for this study, (2) previously collected and published in Dutch literature in Cirkel et al (2015) and (3) partly reported previously in Eggenkamp (1994) and now re-analyzed. bd = below detection limit

Shallow aquifers

All but one of the shallow groundwater samples with methane concentrations greater than 1 mg L⁻¹ had a δ¹³C-CH₄ <−60‰, indicative of a biogenic origin, with the exception of the sample collected near the blowout in Sleen (SGW-09). A follow-up study carried out at this location showed that leakage of natural gas has likely been ongoing for several decades. As such, the sample proved to be clear evidence of an anthropogenically opened connection with the deep subsurface (Schout et al., 2017). The δD-CH₄ values for shallow groundwater samples with biogenic methane in concentrations >1 mg L⁻¹ were around -250‰ or higher, which suggests that CO₂ reduction is the dominant methanogenic pathway (Whiticar, 1999). This was also found to be the dominant methanogenic process in numerous other groundwater methane baseline studies worldwide (e.g. Humez et al., 2016; Kulongoski et al., 2018; McIntosh et al., 2014; Nicot et al., 2017).

Five samples from the shallow aquifers had δ¹³C-CH₄ values above −50‰ (Table 3), typically associated with a thermogenic origin. For two of the samples, the combination of δ¹³C-CH₄ and the light alkane ratio (C₁/(C₂ + C₃) even closely resembles that of SGW-09 from the blowout site, which in turn is similar to that of natural gas reservoirs in the Netherlands (Fig. 6A). However, these fives samples all show extreme deuterium enrichment with δD-CH₄ values between 100 and 300‰ (Fig. 6B). Such highly enriched values have previously been associated with fractionation caused by methane oxidation (Etiope et al., 2011; Schloemer et al., 2016). Rayleigh fractionation modelling with carbon isotope fractionation factors ε_C of 7 and 12 (definitions in Whiticar, 1999) confirms that this shift in isotopic and molecular composition of the dissolved gasses in these samples can indeed be attributed to methane oxidation (Fig. 6A).

Figure 6. Methane carbon isotope ratio versus light alkane ratio (A) and hydrogen isotope ratio (B), for gasses found in Dutch shallow groundwater, deep groundwater, geothermal wells and natural gas accumulations. Dashed black lines show hypothetical mixing lines for biogenic and thermogenic gas, with bold labels in A showing the fraction of thermogenic gas in the mixture. Green and orange lines in A depict the calculated effect of anaerobic methane oxidation (AOM) for two carbon fractionation factors (ε_C). Orange line in B depicts isotopic fractionation up to a residual methane fraction of 1% with fractionation factors ε_C = 12 and ε_D = 96. Calculations according to formulations in Whiticar (1999).

The process by which these anomalous isotopic values came to be is best illustrated by the analyses of the four samples collected in the multi-level monitoring wells SGW-02 and SGW-16. For example, only the sample collected from the greatest depth in SGW-02 has high methane concentrations (SGW-02d, 36 mg L⁻¹), negligible sulphate concentrations and a clear biogenic isotopic composition. The three other samples collected from the overlying well screens all have trace amounts of methane and high sulphate concentrations (Table 3). The extreme deuterium enrichment in these samples is well explained by isotopic fractionation factors ε_C of 12 and ε_D of 96 (Fig. 6B), that fall within the observed range of fractionation factors for microbial methane consumption (Whiticar, 1999), leading to an enrichment ratio (ΔD/Δ¹³C) of 8. Using this method, an estimated 99% of the initial methane would have had to be oxidised to reach such extremely enriched values. Given the low CH₄ concentrations (<0.1 mg L⁻¹) this still corresponds to only a minor amount of oxidation in absolute terms. Hence, rather than a thermogenic origin and connection with the deep subsurface, methane in these wells originates from biogenic sources and has subsequently been transformed isotopically as a result of anaerobic oxidation of methane (AOM).

Paleogene and older deep aquifers

In spite of the low methane concentrations and the presence of sulphate in all deep groundwater samples, only sample DGW-09 displays a similar shift in both isotopic and molecular composition that clearly shows the impact of AOM (Fig. 6B). The other samples (except two; see discussion below) lie somewhere on a mixing line, and hence the small methane mass concentrations appear to consist of varying fractions of biogenic and thermogenic methane (Fig. 6A). Measurements of the δ¹³C of ethane and propane corroborate the conclusion of a mixed origin for these samples, as their δ¹³C-C₂H₆ composition resembles that of Dutch natural gas samples closely, and is notably different from the ethane in the shallow groundwater samples with biogenic methane (Fig. 7). This indicates that some upward seepage of thermogenic gas does occur. The ratio between the δ¹³C of methane and ethane, as incorporated in the isotope factor O_(C1/C2) proposed by Cesar et al. (2021), further suggests that the thermogenic ethane in most of these samples is very early mature. This shows that the source is likely more local than the deeper natural gas reservoirs, which contain more mature thermogenic gas (Fig. S4).

Figure 7. ‘Gas plot’ of carbon isotopic composition versus carbon number for methane (C₁), ethane (C₂) and propane (C₃) in gasses found in Dutch shallow groundwater, deep groundwater and natural gas accumulations. Also shown are potential effects of oxidation and mixing on methane isotopes and ranges of δ¹³C₂H₆ values interpreted to represent either biogenic or thermogenic ethane.

Since the fraction of ethane in thermogenic gas is much larger than in biogenic gas, even mixtures with a small thermogenic component should show a clear thermogenic ethane signature, in conjunction with a mixed isotopic signature of methane which shifts more or less linearly with the fraction of thermogenic methane. Therefore, the lighter ethane signature in the shallow groundwater samples also provides evidence for a lack of even small amounts of thermogenic methane at these depths. This highlights the added value of δ¹³C-C₂H₆ as a tracer for gas migration. However, the amount of sampled ethane is often too small to reliable measure its isotopic composition, which also explains the missing data for a number of samples in this study (Table 3). Increasing the sampled volume could circumvent this issue, however, this may significantly increase sampling efforts depending on its methodology. Furthermore, similar to methane oxidation, preferential oxidation of ethane may also occur in freshwater aquifers (Schout et al., 2017), which can obscure the interpretation of ethane isotopes.

Two exceptions to a mixed or biogenic origin are the samples from wells DGW-10 and DGW-11 (Fig. 6). DGW-10 is located in the east of the Netherlands and screened at a depth of nearly 700 m (Table 3). It has an isotopic composition that is comparable to that of natural gas reservoirs in the Netherlands and therefore presents clear evidence of upward seepage of thermogenic methane up to at least this depth (Fig. 6B). The C₁/(C₂ + C₃) ratio of 5, slightly lower than that of Dutch natural gas, likely indicates that some oxidation has occurred during or after migration. Sample DGW-11 has an enriched δ¹³C-CH₄ composition (−14‰) that is similar to some of the oxidised shallow samples. However, the combination with its δD-CH₄ value of −314‰ places it outside of the traditionally delineated areas for thermogenic and/or biogenic methane, and actually suggests an abiotic methane origin according to classifications in Milkov and Etiope (2018). While such a finding would present a first for the Dutch subsurface (to the authors’ knowledge), additional evidence for this conclusion is presented by the reversed δ¹³C trend of the light alkanes (δ¹³C-C₁ > δ¹³C-C₂ > δ¹³C-C₃, Fig. 7), which is a known characteristic of abiotic gas (Lollar et al., 2006). While no clear evidence is given in the lithological description of the cuttings collected while drilling well DGW-11, intrusions of volcanic rocks are known to be present in the Carboniferous and Permian formations that lie at the base of this well, such as for example the Carboniferous Baarlo Formation (Sissingh, 2004). In the presence of such intrusions, the low concentration of methane in this sample could be the result of catalytic CH₄ generation. In experiments with serpentinised rocks, this process has been shown to occur down to temperatures as low as 20°C (Etiope and Ionescu, 2015), considerably lower than the estimated 37°C in well DGW-11 (assuming a geothermal gradient of 31°C km⁻¹; Verweij et al., 2018). Given its unusual composition compared to the gas from all other groundwater samples, additional evidence would be needed to verify the abiotic origin of dissolved methane in DGW-11.

Dissolved gas in the geothermal wells is distinct from that of Dutch gas wells and has a geochemical signature that falls within the range of values associated with mixing of biogenic and thermogenic gasses (Fig. 6A). However, the depths (1640–2625 m, Table 3) and estimated temperature range of these wells is not typically associated with biogenic methane accumulations. Moreover, wells GTW-01 to GTW-04 are completed in the same formations as (depleted) Lower Cretaceous oil reservoirs, where biodegraded oils are known to be present (Griffioen et al., 2016b). Hence, an oil associated origin is more probable. The dissolved gas isotopic and molecular signature confirms this hypothesis, as it closely matches that of oil associated methane (Milkov and Etiope, 2018). Additionally, the shift in isotopic and molecular composition observed for GTW-02 and GTW-06 matches the shift that is observed for methanogenic biodegradation of oils (Fig. 6). As these are the two shallowest sampled geothermal wells (with estimated temperatures of 60.8 and 72.9°C, respectively), they fall within the range of temperatures for which microbially mediated biodegradation has been observed (<75°C; Milkov, 2010).

Connectivity deep and shallow subsurface

Shallow groundwater samples were collected near both faults and known gas accumulations, where migration of methane from the deep subsurface is most likely. However, naturally occurring thermogenic methane was only observed below the Neogene and Paleogene clay barriers that form the hydrogeological base in the Netherlands. Although no evidence for gas phase accumulations below these layers were found in this study, this suggests that these geological units form an important barrier to upward gas seepage. Their importance in controlling fluid migration is supported by the δ³⁷Cl data for dissolved chloride. Given the non-reactive nature of chloride in aquatic environments, most natural waters tend to have a δ³⁷Cl close to 0‰ (Eggenkamp, 1994). However, since the diffusion coefficient of ³⁵Cl is slightly larger than that of ³⁷Cl, chloride in groundwater with a lighter ³⁷Cl/³⁵Cl composition (more negative δ³⁷Cl value) has been attributed to diffusion-controlled transport (Beekman et al., 2011; Desaulniers et al., 1986).

In this study, a number of samples collected in the south-west of the Netherlands from Paleogene aquifers show very negative δ³⁷Cl values: down to −2.08‰ (Fig. 8). These sandy aquifers are bounded on the top and bottom by low permeable marine clays. They are recharged by fresh rain water infiltrating in neighbouring Belgium, with chloride concentrations increasing from fresh to brackish as the groundwater flows towards the north and west and into the Netherlands (Van Der Kemp et al., 2000). The negative δ³⁷Cl values in these samples therefore point out that transport of solutes is predominantly controlled by diffusion from the under- and overlying clays, rather than advective transport. As chloride concentrations increase further, δ³⁷Cl values increasingly resemble that of the (connate) seawater source in the confining units (0‰), resulting in the observed positive correlation between δ³⁷Cl and Cl.

Figure 8. δ³⁷Cl chloride stable isotope ratio versus chloride concentration (A) and versus percent Modern Carbon (B), for deep and shallow Dutch groundwater samples. Symbols show the age of the stratigraphic unit. Red dashed line is the chloride concentration of modern seawater.

This positive correlation is also observed for the Cretaceous formation waters (Fig. 8). However, chloride concentrations in these aquifers exceed that of seawater. The elevated chloride concentrations result from dissolution of halite, which has δ³⁷Cl values of around 0‰ (Eggenkamp et al., 2019). Subsequent diffusive transport then leads to the observed negative δ³⁷Cl values. Besides the Cretaceous samples, one of the samples from a Paleogene unit in the north eastern corner of the country (DGW-07) also fits this pattern. Here, salinities exceeding that of seawater were indeed shown to result from dissolution of a nearby evaporite diapir (Griffioen et al., 2016b).

Radiocarbon dating shows that groundwater in the Paleogene and older aquifers contains ¹⁴C well above detection limit (0.01–0.1 pMC, Table 3) and a significant fraction of relatively recent groundwater (<40,000 years old) must therefore be present. Significant fractions of modern carbon were also observed in the groundwater samples from the Jurassic, Cretaceous and Carboniferous formations. This confirms that groundwater in these formations is not stagnant, and that some recharge is happening or has happened during the past tens of thousands of years. National-scale groundwater modelling by Valstar and Goorden (2016) indeed pointed out that exchange of groundwater across the Paleogene clay layers is on-going. However, as diffusion appears to be the dominant transport mechanism in these units, significant advective transport of dissolved thermogenic methane from below these layers to the shallow subsurface is unlikely.

Besides the presence of the Paleogene clays, the absence of thermogenic methane in shallow groundwater may also have resulted from a number of other processes that limit its upward migration. First, gas seeping upward through faults that penetrate the sandy Quaternary aquifers that characterise most of the Dutch shallow subsurface would likely become unfocused and hence harder or even impossible to detect (Thielemann et al., 2000). Second, lateral groundwater flow in the aquifers encountered by the migrating gas phase may lead to significant dilution and ultimately stagnation of the upward seeping gas phase, preventing detection in overlying aquifers (Schout et al., 2020). Hence, while no evidence for methane seepage to the shallow groundwater through natural connections with the deep subsurface was encountered, its occurrence cannot be entirely excluded on the basis of this study. Nevertheless, the results of this study do indicate that future observations of thermogenic methane in shallow groundwater are most likely linked to anthropogenic activities.

Conclusions

The origin and distribution of methane in the subsurface of the Netherlands was assessed based on a compiled nationwide dataset of water samples analysed for methane concentrations (n = 12,219) and in conjunction with a dataset of targeted groundwater samples from both shallow aquifers (<500 mbgs) and deep reservoirs that were also analysed for the molecular and isotopic composition of dissolved gasses (n = 40).

The presence of methane is ubiquitous in the shallow Quaternary groundwater system, which is hydraulically separated from deeper groundwater compartments by a succession of thick Neogene and Paleogene marine clays. Concentrations in 7.5% of samples exceed 10 mg L⁻¹ and median concentrations are 0.2 mg L⁻¹, relatively high from an international perspective. The highest methane concentrations occur in the coastal provinces, where Holocene peats and marine clays overly the Quaternary aquifers. There, median concentrations reached 4.0 mg L⁻¹ between 20 to 40 m bgs. However, at greater depths (∼100–160 m bgs) median methane concentrations are lower in the Holocene coastal part of the country than in inland areas, most likely because of the presence of seawater associated dissolved sulphate at these depths preventing microbial methanogenesis in buried Pleistocene aquifers. Besides these shallow Holocene formations, other prominent sources of dissolved methane in Dutch shallow groundwater appear to be the marine Eem Formation and the glacial Peelo Formation.

Isotopic analysis of samples from the shallow groundwater system confirmed that naturally present methane has a biogenic origin, with no perceivable contribution by migration of thermogenic through connections with deeper formations, even though the samples were collected at locations where such connections were deemed most likely (i.e. close to major faults and overlying gas reservoirs). δD-CH₄ of samples with methane concentrations >1 mg L⁻¹ were either around −250‰ or larger, indicating that CO₂ reduction rather than methyl-type fermentation is the dominant methanogenic pathway. A number of samples with low biogenic methane concentrations (<0.1 mg L⁻¹) were severely impacted by anaerobic methane oxidation, resulting in extremely enriched δD-CH₄ values (100–300‰) and δ¹³C-CH₄ and C₁/(C₂ + C₃) values that resemble that of thermogenic gas. In the absence δD-CH₄ analysis, this could therefore lead to incorrect fingerprinting of methane origin.

The composition of dissolved gasses in the majority of samples from deeper groundwater wells (up to 871 m bgs) had a mixed isotopic signature. A clear thermogenic ethane isotopic signature in these samples confirmed that at least trace amounts of thermogenic methane are present at these depths. This shows that some upward seepage of thermogenic methane does occur, or that very early mature ethane is produced more locally. Given that thermogenic methane was not observed in the shallow groundwater above the Paleogene marine clay barriers, they are believed to play a role in constraining further upwards transport of thermogenic gas. Negative δ³⁷Cl values observed in the Paleogene sandy aquifers that are situated in between such clays (up to -2.08‰, vs. SMOC) corroborate this conclusion and indicate that solute transport into these formations is diffusion rather than advection controlled. The results of this study highlight that future observations of thermogenic methane in shallow groundwater of the Netherlands are most likely attributable to migration through anthropogenically induced rather than natural conduits. The presence of such conduits should be considered with the increased use of the deep subsurface for geothermal energy production and CO₂ or hydrogen storage.