Measurement report: Isotopic composition of CH₄ emitted from gas exploration sites in the Transylvanian Basin, Romania

Isotope measurements are increasingly used to constrain the methane (CH₄) budget on various scales, from global to regional. The success of isotope-based source attribution depends to a large degree on the knowledge of the isotope signatures of the various source categories at the point of emission, but this information is in many cases lacking. Here we report the isotopic composition of CH₄ emitted from 48 installations in the gas production region of Transylvania, Romania. The isotopic source signatures are quite homogeneous across the basin with average values of δ¹³C = (−65.6 ± 0.5 ‰) and δD = (−184 ± 1 ‰) confirming the predominantly biogenic origin of the Transylvanian gas, produced by hydrogenotrophic CO₂ reduction. This is similar to values reported previously from natural seeps in Transylvania, to the natural gas exploited in the Dolj region in Southwestern Romania, and to the natural gas in the distribution grid in Cluj-Napoca. However, is more depleted in heavy isotopes than the oil-associated gas emitted in the Southern Romanian Plain, and gas leakages in the city of Bucharest. In addition, we present a step-by-step derivation of the underlying “Keeling plot” mass balance approach that is used to derive isotope source signatures.

1 Introduction

Methane is a strong greenhouse gas and it is important to reduce its emissions to the atmosphere in order to reach the goals of the Paris climate agreement (Nisbet et al., 2019, 2020; Ocko et al., 2021). Emissions from the fossil fuel sector are considered low-hanging fruit in that respect, since a large share of the emissions can be mitigated at little or even no cost (United Nations Environment Programme and Climate and Clean Air Coalition, 2021; Höglund-Isaksson et al., 2020). A prerequisite for emission reduction is knowledge of where the emissions are, which requires direct observations across the value chain. The ROMEO project (ROmanian Methane Emissions from Oil and gas) aimed to identify, attribute and quantify emissions from the oil and gas production infrastructure in Romania, one of the European Union's largest oil and gas production regions. Intensive measurement campaigns with ground-based observations were carried out in 2019 in the South Romanian Plain (Stavropoulou et al., 2023; Delre et al., 2022) and in 2021 in the Transylvanian Basin, manuscript in preparation). In addition, aircraft-borne measurements were carried out to constrain the emissions by both in-situ measurements in 2019 (Maazallahi et al., 2025) and remote sensing in 2021 (Kuhlmann et al., 2025). Those measurements collectively demonstrated that emissions from oil and gas operations in Romania are severely underestimated in National reporting.

The isotopic composition of methane (CH₄) can be used to distinguish CH₄ that is produced via different pathways (Sherwood et al., 2017; Schwietzke et al., 2016; Menoud et al., 2022a; Quay et al., 1999; Brenninkmeijer et al., 2003; Whiticar, 2020; Milkov and Etiope, 2018; Sherwood Lollar et al., 2006; Ojeda et al., 2023). Thermogenic CH₄ is usually associated with relatively high δ¹³C values between −55 ‰ and −30 ‰, and δD generally varies in the range −250 ‰ to −100 ‰. The isotopic composition of fossil reservoirs is additionally influenced by the composition of the fuels (“dry” gas, or “wet” CH₄ in association with oil) and the reservoir maturity (Whiticar, 2020; Menoud et al., 2022b; Milkov and Etiope, 2018). Biogenic CH₄ formed via the hydrogenotrophic pathway is more depleted in ¹³C (δ¹³C between −100 ‰ and −60 ‰) whereas it has relatively similar δD values as thermogenic CH₄ (between −250 ‰ and −150 ‰). Biogenic CH₄ formed via the acetoclastic pathway has δ¹³C values between −70 ‰ and −50 ‰ and is generally depleted in deuterium (δD lower than −250 ‰). Abiotic CH₄ is relatively enriched in ¹³C (δ¹³C > −40 ‰) and can cover a wide δD range between −50 ‰ and −400 ‰. Pyrogenic CH₄ produced mostly during biomass burning is also enriched in both ¹³C and D (δ¹³C between −30 ‰ and −10 ‰; δD > −250 ‰) and thus falls in a similar range as abiotic CH₄ (Sherwood et al., 2017; Menoud et al., 2022b; Whiticar, 2020).

Before the ROMEO campaigns, the existing methane isotopic data from oil and gas fields in Romania were limited to geologic natural emissions related to natural seepage. The measured values (discussed in detail below) may constitute a reference for the current work, as the mentioned sites are representative for the type of gas deposits in the study area.

Denser and systematic isotopic investigations all over Romania, combining analysis of samples collected from boreholes and from surface manifestations, are needed in order to geochemically characterize the hydrocarbon deposits, and also to better understand the methane transfer to the atmosphere from the oil and gas industry. During the ROMEO campaigns, air samples were collected in emission plumes to investigate the origin of the emitted CH₄ in more detail using stable isotope analysis. Menoud et al. (2022b) reported the isotopic composition of samples collected at 83 ground locations and 24 samples collected on aircraft flights. They showed a wide range of isotope signatures and confirmed that the gas across the Romanian Plain is mostly associated with oil production and of thermogenic origin with average values of δ¹³C = −50 ‰ and δD = −189 ‰. A few reservoirs of microbial origin were also found. Overall, the isotope composition of gas emitted from oil and gas production sites in Romania was significantly more depleted in ¹³C than commonly used values for the global fossil fuel emissions.

This study aims to provide a better isotopic characterization of CH₄ emissions associated with gas production in the Transylvanian Basin. We report the isotopic composition of air samples collected during phase B of the ROMEO project, conducted in summer 2021, at 48 individual gas production locations across the Transylvanian basin.

Figure 1 Orographic map of the campaign area, highlighting the sub-surface gas fields (pink areas) and the sampling locations for the samples collected in this study. The 12 small dual isotope plots on the left and right indicate the distribution of the individual source signatures derived for the 12 gas fields that were visited. The orange diamond is a common reference point (δ¹³C = −65 ‰, δD = −185 ‰). (Location of the gas fields: modified after Map of the mineral resources, 1:1 000 000, Institute of Geology and Geophysics, Bucharest 1984).

2 Methods

2.1 Campaign region

The Transylvania region is located in the central part of Romania, enclosed between the Apuseni Mountains in the West and the Eastern and Southern Carpathians. With over 100 gas fields scattered throughout the Transylvanian Basin, it remains the foremost gas producer among Central and South-Eastern European countries. The main petroleum system in the Transylvanian Basin corresponds to Neogene deposits, mainly hosting microbial methane (Popescu, 1995; Krézsek et al., 2010). Over the past century, it has yielded an estimated 30 TCF (trillion cubic feet) of gas (Krézsek, 2011), while undiscovered and confirmed reserves amount to approximately 20 TCF (Pawlewicz, 2005) (1 m³ = 35.315 cubic feet). There has been no discovery of commercial oil in the region of concern. Figure 1 shows the distribution of the gas reservoirs in the central-western part of the Transylvanian Basin, and the locations where air samples were collected for this study.

The study area is located in the central part of the Transylvanian Basin, a back-arc basin which is characterized by a substantial accumulation of Middle-Upper Miocene detrital sediments (Badenian to Pannonian). These deposits, formed due to fast subsidence, may exceed 5000 m in certain areas of the basin. A regressive event in the Middle Badenian provided optimal conditions for the accumulation of a substantial salt layer, potentially reaching 300 m in thickness (Krézsek and Filipescu, 2005). The salt tectonics is responsible for the creation of brachyanticlines in the central part of the basin, and diapirs on the margins (Tiliţă et al., 2013). The commercial gas plays are mainly associated with brachyanticlines within the post-salt Badenian–Sarmatian deposits, featuring mild flank dips, typically ranging from 2 to 6°. These multi-layered structures may encompass up to 15 gas-bearing intervals, or even more in particular cases (e.g. Filitelnic – 23 pay levels) (Paraschiv, 1979). The depth of the pay intervals varies significantly, ranging from several hundred meters to over 3000 m.

2.2 Sample collection and isotope measurement

When emission plumes had been identified during ground-based surveys with real-time CH₄ sensors, air samples were collected in the plumes by the ground teams. Air samples were pumped into 2 L volume flex-foil bags using a small pump (KNF Neuberger) via Teflon tubing and a Magnesium Perchlorate dryer. Usually, two air samples were collected in the emission plume downwind of gas production installations, and several more “background samples” in clean air in the respective region, in order to determine the isotopic source signature by a Keeling plot approach (see below). A total of 96 samples from 48 individual production installations and 30 background samples were collected from the study area between 13 June to 4 July, 2021. Figure 1 shows a map of the production regions and the locations where samples were collected.

The isotopic composition of CH₄ (δ¹³C and δD) in the air samples was analysed at Utrecht University using a continuous-flow isotope ratio mass spectrometry system (Brass and Röckmann, 2010; Menoud et al., 2022b). First, CH₄ is separated from ambient air samples and purified using temperature-controlled traps and gas chromatography. The isotopic composition of the purified CH₄ is then determined in an isotope ratio mass spectrometer. Individual measurements have a precision better than 0.1 ‰ for δ¹³C and 2.0 ‰ for δD (Brass and Röckmann, 2010; Röckmann et al., 2016; Menoud et al., 2020). The system has been carefully calibrated and participated in inter-laboratory comparisons (Brass and Röckmann, 2010; Umezawa et al., 2018). It has been used in numerous previous projects to characterize CH₄ isotopic composition (Röckmann et al., 2011, 2016; Maazallahi et al., 2020; Menoud et al., 2020, 2021, 2022a; Lu et al., 2021; Fernandez et al., 2022; Fiehn et al., 2023).

2.3 Determination of isotopic source signatures

When emissions of CH₄ into the atmosphere lead to a clearly measurable enhancement in the CH₄ mole fraction, the observed mole fraction (mf_obs) is the sum of a background component (mf_bg) and a source (mf_src) component according to the mathematical equation

$$\begin{array}{}\text{(1)}& {\text{mf}}_{\text{obs}}={\text{mf}}_{\text{bg}}+{\text{mf}}_{\text{src}}.\end{array}$$

A similar equation is valid for each individual isotopologue, e.g. for the mole fraction of the ¹³C-substituted CH₄ ¹³mf:

$$\begin{array}{}\text{(2)}& {}^{\mathrm{13}}{\text{mf}}_{\text{obs}}={}^{\mathrm{13}}{\text{mf}}_{\text{bg}}+{}^{\mathrm{13}}{\text{mf}}_{\text{src}}.\end{array}$$

The equation for Deuterium-substituted CH₄ is exactly analogous and not shown. These two mass conservation equations can be combined to determine the isotope signature of the source that is responsible for the observed mole fraction enhancement. The approach first used in (Keeling, 1961) assumes that the background component (both mole fraction and isotopic composition) remain constant over the course of the measurement. This is valid for measurements carried out over a short time close to strong emitters like the ones presented below.

Equation (2) can be written as

$$\begin{array}{}\text{(3)}& \begin{array}{rl}& {\displaystyle \frac{{}^{\mathrm{13}}{\text{mf}}_{\text{obs}}}{{}^{\mathrm{12}}{\text{mf}}_{\text{obs}}}}{}^{\mathrm{12}}{\text{mf}}_{\text{obs}}={\displaystyle \frac{{}^{\mathrm{13}}{\text{mf}}_{\text{bg}}}{{}^{\mathrm{12}}{\text{mf}}_{\text{bg}}}}{}^{\mathrm{12}}{\text{mf}}_{\text{bg}}+{\displaystyle \frac{{}^{\mathrm{13}}{\text{mf}}_{\text{src}}}{{}^{\mathrm{12}}{\text{mf}}_{\text{src}}}}{}^{\mathrm{12}}{\text{mf}}_{\text{src}}\\ & {{}^{\mathrm{13}}R}_{\text{obs}}{}^{\mathrm{12}}{\text{mf}}_{\text{obs}}={{}^{\mathrm{13}}R}_{\text{bg}}{}^{\mathrm{12}}{\text{mf}}_{\text{bg}}+{{}^{\mathrm{13}}R}_{\text{src}}{}^{\mathrm{12}}{\text{mf}}_{\text{src}}.\end{array}\end{array}$$

where R is the heavy-to light isotope ratio (in this example ¹³R = ${}^{\mathrm{13}}\mathrm{C}/{}^{\mathrm{12}}\mathrm{C}$). When the heavy isotope has a much lower abundance than the light isotope, the approximation ¹²mf∼mf is valid and Eq. (3) can be approximated as

$$\begin{array}{}\text{(4)}& {{}^{\mathrm{13}}R}_{\text{obs}}{\text{mf}}_{\text{obs}}={{}^{\mathrm{13}}R}_{\text{bg}}{\text{mf}}_{\text{bg}}+{{}^{\mathrm{13}}R}_{\text{src}}{\text{mf}}_{\text{src}}.\end{array}$$

By dividing Eq. (4) by the isotope ratio of the international standard (index ST) and subtracting Eq. (1), Eq. (4) can be formulated in terms of δ values as follows:

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& {\displaystyle \frac{{{}^{\mathrm{13}}R}_{\text{obs}}}{{{}^{\mathrm{13}}R}_{\text{ST}}}}{\text{mf}}_{\text{obs}}={\displaystyle \frac{{{}^{\mathrm{13}}R}_{\text{bg}}}{{{}^{\mathrm{13}}R}_{\text{ST}}}}{\text{mf}}_{\text{bg}}+{\displaystyle \frac{{{}^{\mathrm{13}}R}_{\text{src}}}{{{}^{\mathrm{13}}R}_{\text{ST}}}}{\text{mf}}_{\text{src}}\\ & {\displaystyle \frac{{{}^{\mathrm{13}}R}_{\text{obs}}}{{{}^{\mathrm{13}}R}_{\text{ST}}}}{\text{mf}}_{\text{obs}}-{\text{mf}}_{\text{obs}}={\displaystyle \frac{{{}^{\mathrm{13}}R}_{\text{bg}}}{{{}^{\mathrm{13}}R}_{\text{ST}}}}{\text{mf}}_{\text{bg}}-{\text{mf}}_{\text{bg}}\\ & +{\displaystyle \frac{{{}^{\mathrm{13}}R}_{\text{src}}}{{{}^{\mathrm{13}}R}_{\text{ST}}}}{\text{mf}}_{\text{src}}-{\text{mf}}_{\text{src}}\\ & \left({\displaystyle \frac{{{}^{\mathrm{13}}R}_{\text{obs}}}{{{}^{\mathrm{13}}R}_{\text{ST}}}}-\mathrm{1}\right){\text{mf}}_{\text{obs}}=\left({\displaystyle \frac{{{}^{\mathrm{13}}R}_{\text{bg}}}{{{}^{\mathrm{13}}R}_{\text{ST}}}}-\mathrm{1}\right){\text{mf}}_{\text{bg}}\\ & +\left({\displaystyle \frac{{{}^{\mathrm{13}}R}_{\text{src}}}{{{}^{\mathrm{13}}R}_{\text{ST}}}}-\mathrm{1}\right){\text{mf}}_{\text{src}}\\ & {\text{mf}}_{\text{obs}}\times \mathit{\delta }{}^{\mathrm{13}}{C}_{\text{obs}}={\text{mf}}_{\text{bg}}\times \mathit{\delta }{}^{\mathrm{13}}{C}_{\text{bg}}+{\text{mf}}_{\text{src}}\times {\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}}_{\text{src}}.\end{array}\end{array}$$

Noting that mf_bg, ${\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}}_{\text{bg}}$ are assumed to be constant over the period of the measurement, Eq. (5) can be expressed as linear equation of the type $y=mx+a$ where $y={\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}}_{\text{obs}}$, $m={\text{mf}}_{\text{bg}}\times ({\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}}_{\text{bg}}-{\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}}_{\text{src}})$ and $x=\frac{\mathrm{1}}{{\text{mf}}_{\text{obs}}}$.

$$\begin{array}{}\text{(6)}& {\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}}_{\text{obs}}={\displaystyle \frac{{\text{mf}}_{\text{bg}}\times \left({\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}}_{\text{bg}}-{\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}}_{\text{src}}\right)}{{\text{mf}}_{\text{obs}}}}+{\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}}_{\text{src}}\end{array}$$

The Keeling plot approach is a graphical approach where a linear fit is applied to a correlation plot of ${\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}}_{\text{obs}}$ versus $\frac{\mathrm{1}}{{\text{mf}}_{\text{obs}}}$, and the y-axis intercept of the linear fit equation then returns the isotopic signature of the source, ${\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}}_{\text{src}}$.

Figure 2 Dual isotope plot (δD versus δ¹³C) of all individual isotope source signatures derived for the different sampling locations across the Transylvanian Basin. The different symbol colours represent different gas fields.

We note that in cases where mf_bg and ${\mathit{\delta }{}^{\mathrm{13}}\mathrm{C}}_{\text{bg}}$ are not constant, but can be specified (e.g. for analysis of longer time series), the mass conservation equations Eqs. (1) and (2) can be rewritten differently in the so-called Miller–Tans approach (Miller and Tans, 2003) to determine isotope source signatures. The differences between the two approaches have been investigated in detail recently (Defratyka et al., 2025). That study also investigated the effect of different mathematical methods to apply linear regression analysis to Eq. (6). In our study we use the orthogonal distance regression method (Boggs et al., 1988).

3 Results

Figure 2 shows a dual isotope plot of all source signatures determined at the individual sampling locations visited in this study. The locations are color-coded by gas field. Table 1 provides the numerical values, and the individual δ¹³C values are also shown as color-coding in the map of Fig. 1. It is evident that the isotopic composition of the gas produced in the investigated part of the Transylvanian basin is quite homogeneous. Most of the δ¹³C source signatures at individual sampling locations fall in a range between −70 ‰ to −60 ‰. δD values at most individual locations fall within a narrow range between −200 ‰ and −280 ‰. This characterizes the Transylvanian gas as microbial, produced by the hydrogenotrophic pathway. All of the outliers have high uncertainties in the determination of the source signatures, implicating large scatter of individual air samples around the linear fits to the Keeling plots. Such large scatter usually indicates that the assumption of the mass balance model (Eqs. 1 and 2) may not be met, and in many cases this is because of other interfering sources. For example, the “high” outlier at the Dumbravioara gas field may be caused by an interference from combustion emissions.

Figure 3 Isotope source signatures of each individual gas field visited during the ROMEO-B campaign for δD (top) and δ¹³C (bottom). The circles are the individual measurements, the boxes represents the 25–75 percentile of values, the mean is indicated as horizontal line, and the whiskers show the 5 %–95 % percentiles. Note that the boxes are more an indication than a statistically robust evaluation due to the limited number of measurements. The horizontal red line shows the source signature derived from a single Keeling plot analysis combining all samples. The final category shows the distribution of all individual source signatures estimates as box plot.

In order to investigate possible differences between gas fields, Fig. 3 shows the derived Keeling plot intercepts for the different gas fields, i.e., all samples from a certain gas field were combined. Figure 3 also includes the average source signature that is derived when all of the samples are combined in one single Keeling plot analysis.

A Keeling analysis of all samples collected across the Transylvanian basin returns a y-axis intercept of δ¹³C = (−65.6 ± 0.5 ‰) and δD = (−184 ± 1 ‰), where the uncertainties state the 1σ uncertainty of the intercept. This is the average isotopic composition of the CH₄ emitted from the investigated gas production installations across the Transylvanian basin. The complete dataset is available at https://doi.org/10.18160/4SJW-ST8W (Röckmann et al., 2025).

4 Discussion and conclusion

Methane from several surface gas manifestations in the Transylvanian Basin has been isotopically characterized previously, showing a distinctive microbial footprint in the case of seeps for the central part of the basin, with δ¹³C in the range of −60.3 ‰ to −67.4 ‰ and δD between −189.5 ‰ and −192.2 ‰ (Etiope, 2009; Baciu et al., 2018). The isotopic source signatures obtained from gas production infrastructure presented above are very similar to the natural gas seeps in central Transylvania, suggesting that gas that escapes via natural seepage in Romania originates from similar underground reservoirs as the produced gas characterized in the present study. The formation pathway that is typically associated with these source signatures is hydrogenotrophic CO₂ reduction (Whiticar, 2020; Milkov and Etiope, 2018; Menoud et al., 2022a). The isotopic composition is also similar to what is found for identified gas leakages from the gas distribution network in the city of Cluj-Napoca (van Es et al., 2024), confirming that the gas in the distribution grid originates from Transylvanian reservoirs. The present data allows to better define the Transylvanian petroleum system in terms of CH₄ isotopic composition with respect to the global CH₄ isotopic dataset.

Reported isotope values show that CH₄ from surface seepage becomes more thermogenic and enriched in N₂, CO₂ and He towards the eastern margin of the Transylvanian Basin, due to the thermal influence of the volcanic range of the Eastern Carpathians (Etiope et al., 2011). By contrast, the few available CH₄ isotopic analyses from the exterior of the Carpathian arcuated range have shown a dominant thermogenic origin of the gases. To this category belong the mud volcanoes and everlasting fires from the Carpathian Fysch and Foredeep, as Pâclele, Fierbători, Beciu, Andreiaşu (Etiope, 2009), Răiuţi, Lopătari, Lepşa (Baciu et al., 2018), or the seep from Bacău–Moldavian Platform (Baciu et al., 2008).

The emitted CH₄ associated with oil production in the southern part of Romania has a very similar δD signature (−189 ± 38 ‰) but a very different δ¹³C signature (−49.7 ± 6.4 ‰) (Menoud et al., 2022b), confirming the different (thermogenic) sub-surface formation pathway. During the ROMEO-A city campaign in the city of Bucharest, located in the southern part of the country, Fernandez et al. (2022) measured methane isotopic values of δ¹³C = −50 ‰ and δD = −196 ‰ from leakages in the natural gas distribution system in Bucharest, similar to the associated gas in Southern Romania.

In summary, CH₄ emitted from 48 gas production sites in the Transylvanian basin exhibits a homogeneous isotopic composition of δ¹³C = (−65.6 ± 0.5 ‰) and δD = (−184 ± 1 ‰), confirming the biogenic origin of the gas in the central-western region of the basin.