Monitoring crustal CO₂ flow: methods and their applications to the mofettes in West Bohemia

Monitoring of CO₂ degassing in seismoactive areas allows the study of correlations of gas release and seismic activity. Reliable continuous monitoring of the gas flow rate in rough field conditions requires robust methods capable of measuring gas flow at different types of gas outlets such as wet mofettes, mineral springs, and boreholes. In this paper we focus on the methods and results of the long-term monitoring of CO₂ degassing in the West Bohemia/Vogtland region in central Europe, which is typified by the occurrence of earthquake swarms and discharge of carbon dioxide of magmatic origin. Besides direct flow measurement using flowmeters, we introduce a novel indirect technique based on quantifying the gas bubble contents in a water column, which is capable of functioning in severe environmental conditions. The method calculates the mean bubble fraction in a water–gas mixture from the pressure difference along a fixed depth interval in a water column. Laboratory tests indicate the nonlinear dependence of the bubble fraction on the flow rate, which is confirmed by empirical models found in the chemical and nuclear engineering literature. Application of the method in a pilot borehole shows a high correlation between the bubble fraction and measured gas flow rate. This was specifically the case for two coseismic anomalies in 2008 and 2014, when the flow rate rose during a seismic swarm to a multitude of the preseismic level for several months and was followed by a long-term flow rate decline. However, three more seismic swarms occurring in the same fault zone were not associated with any significant CO₂ flow anomaly. We surmise that this could be related to the slightly farther distance of the hypocenters of these swarms compared to the two ones which caused the coseismic CO₂ flow rise. Further long-term CO₂-flow monitoring is required to verify the mutual influence of CO₂ degassing and seismic activity in the area.

1 Introduction

Long-term monitoring of crustal fluids activity provides a unique opportunity to better understand the relationships among tectonic processes, seismic activity, and migration of fluids in the Earth's crust. Carbon dioxide of deep origin represents a link between deep-seated magmatic sources of CO₂, the fluid migration paths in the crust, which are controlled by the tectonic stress field, and the Earth surface's. The presented study is focused on the monitoring of CO₂ degassing in the West Bohemia/Vogtland area, which is located in the western part of the Bohemian Massif (BM), the largest coherent surface exposure of basement rocks in central Europe. The western BM is hosting a junction of three tectonometamorphic units, Saxothuringian, Teplá–Barrandian, and Moldanubian (Franke, 2000). It is intersected by two regional tectonic structures, the NE–SW-trending Eger Rift (ER) and NNW–SSE-trending Mariánské Lázně Fault (MLF) (Fig. 1).

Figure 1 Relief map of the West Bohemia/Vogtland region with the fault network, granite units, and major fault zones as Mariánské Lázně Fault (MLF) and Eger Rift. Seismic events, the major focal zone of Nový Kostel (NK area), CO₂ monitoring stations, CO₂ mofettes and mineral springs, and Quaternary volcanoes are also indicated (see the legend).

The Tertiary ER is a 300 km long striking structure characterized by elevated heat flow and Cenozoic volcanism, and its formation is thought to be related to Alpine collision (Ziegler, 1992). The Late Variscan MLF was reactivated several times during the geological history up to Cenozoic when it participated in the formation of the Cheb Basin (CB). CB is typified by a blocky structural fabric due to a network of faults. Besides the NNW and NW morphologically expressed marginal faults, faults striking NE, E–W, and N–S were also identified within the basin (Špičáková et al., 2000; Bankwitz et al., 2003).

The present geodynamic activity is manifested by earthquake swarms, massive CO₂ degassing of mantle origin, and Quaternary volcanism (Fischer et al., 2014). Seismic activity in the form of earthquake swarms is concentrated in the area of CB, in particular the Nový Kostel focal zone (Fig. 1), where more than 80 % of seismic energy is released in frames of the whole seismogenic region. Here the hypocenters form a N–S-trending, steeply dipping belt in the depth range from 6 to 10 km; however, no clear fault outcrop has been identified that would match the focal zone geometry. The prevailing focal mechanisms coincide very well with the orientation of the fault zone striking 169^∘ derived from the hypocenter trend. Inversion of focal mechanisms for stress field yields maximum compression direction in the range N135–155^∘ E, which coincides well with the average direction N144^∘ E in western Europe (Fischer et al., 2014). This direction is, however, parallel to the strike of the MLF, which indicates a passive role of the MLF in the present stress field (Vavryčuk, 2011).

The strongest earthquakes usually do not exceed M_L 4.5, as was the case of all the eight major instrumentally recorded swarms between 1985 and 2018.

The concentration of the geodynamic phenomena in this small region is not clearly understood. Some authors relate this seismic activity to intersecting crustal faults (e.g., Bankwitz et al., 2003) or to fluids of mantle origin (e.g., Bräuer et al., 2003), which could originate from active magmatic underplating (Hrubcová et al., 2017).

The degassing occurs in the form of CO₂-rich mineral waters and wet and dry mofettes in several degassing fields. Carbon dioxide is the carrier phase for mantle-derived minor components such as helium, the isotope ratios of which are the best tool to determine whether the fluids are of crustal or mantle-derived origin; high ³He∕⁴He ratios indicate that ascending gases are of mantle origin (Bräuer et al., 2003). The highest portions of mantle-derived helium (up to 6 R_a, where R_a corresponds to the ³He∕⁴He ratio of the atmosphere) were found in the CB; the Karlovy Vary (KV) degassing center has the lowest ³He∕⁴He ratios of 2.5 R_a. Lower He-isotope ratios (e.g., ³He∕⁴He < 6 R_a) probably reflect the gas mixing with crustal-derived He along fluid pathways (Bräuer et al., 2008). Also, the δ¹³C values in the CO₂-rich gas escapes indicate their origin in the upper mantle (Weinlich et al., 1999; Bräuer et al., 2003). CO₂ flow monitoring in West Bohemia has been conducted since the 1990s in a rather discontinuous way. The longest-running observation project is probably the monitoring of radon activity in Bad Brambach (Heinicke and Koch, 2000; Koch et al., 2003), which has been conducted since 1989.

Gas flow is concentrated in three degassing centers: Cheb Basin, Mariánské Lázně, and Karlovy Vary (KV) (Weinlich et al., 1999; Geissler et al., 2005; Kämpf et al., 2007). They are characterized by a high gas flow with daily discharge of dozens of tons of gas (Nickschick et al., 2015) with CO₂ concentrations of more than 99 vol %. Cheb Basin also has the highest concentration of seismic activity, which makes it ideal for studying the relations between seismicity and gas flow. Interestingly, many studies of the local earthquake swarms show that they may be related to pressurized fluids in the crust and the ascent of gas. This has been pointed out by numerous researchers including Špičák and Horálek (2000), Hainzl and Fischer (2002), Fischer and Horálek (2005), and Hainzl et al. (2016), based on space-time analysis of the seismicity, Horálek et al. (2002), Vavryčuk (2002), and Vavryčuk and Hrubcová (2017), on the basis of moment tensor analysis, and Dahm and Fischer (2014) and Bachura and Fischer (2016), based on V_p∕V_s analysis of the volume of hypocenters. The last two studies show that compressible fluids are required to explain the low velocity ratio observed in the course of seismic activity.

Another long-time monitoring was carried out as part of the “Research of CO₂ pressure field in the area of West Bohemian spas” project funded by the Ministry of the Environment of the Czech Republic from 1996 to 2005. Gas flow in open boreholes and gas pressure in closed boreholes were monitored at 11 gas escape sites in the Cheb Basin and near Mariánské Lázně (Škuthan et al., 2001; Hron and Škuthan, 2006). Monitoring of pressure in a closed well was preferred at many project sites since the functioning of mechanical flowmeters was unreliable due to condensation and freezing. A different type of CO₂ flow monitoring was carried out by Faber et al. (2009), who measured diffuse gas flow by determining CO₂ concentration in soil gas at two stations in the Nový Kostel fault zone. CO₂ flow monitoring was also conducted by Jens Heinicke (personal communication, 2019) in the Bublák mofette from 2008 to 2014 by recording the acoustic noise of bubbles below the water table, a method which is similar to that used by Koch et al. (2003). No convincing observation of seismogenic CO₂ flow anomaly was, however, presented in the abovementioned studies.

Mapping of CO₂ emanations was conducted in the area by Nickschick et al. (2015). They used an infrared gas analyzer and accumulation chambers to measure CO₂ flux and CO₂ soil concentration in the mofette field of Hartoušov and found that the diffuse gas flow in dry vents accounts for a high portion of the mofette field's total gas production.

The measurement of CO₂ flow presented in this paper began in 2009 in the Hartoušov mofette field with the use of a laboratory chamber flowmeter. Despite problems from the condensation of moisture and freezing temperatures, which resulted in time series gaps, we observed a massive post-seismic CO₂ flow increase shortly after the first M_L 3.5 mainshock of the 2014 seismic sequence. A comparison with the fault valve model showed a striking fit, which indicated that the earthquake fracture released gas accumulated in the reservoir beneath hypocenters (Fischer et al., 2017). This gave us a reason to extend the monitoring and test different, more durable, gas flow measurement methods. In this paper, we introduce the principles for our approaches and give a basic comparison of them. We also present the data recorded from the Hartoušov, Bublák, Soos, and Prameny stations (see Fig. 1) and evaluate their response to air pressure and temperature and their possible relation to seismicity.

2 Data and methods

Two types of CO₂ degassing are observed in West Bohemia/Vogtland: (i) diffuse gas flow in soil and (ii) massive gas discharge in mofettes and mineral springs. While gas diffusion in soil is influenced by soil moisture and other local conditions, among other factors, gas flow in massive sources is independent of environmental conditions and should reflect the influence of the gas source in the depth. The deep roots of CO₂ mofettes were also documented by a massive increase in CO₂ flow in the Hartoušov mofette that began about 4 d after the start of seismic activity in 2008 and 2014 (Fischer et al., 2017). This points to the relatively fast speed of gas migration in the upper crust and qualifies mofettes as favorable places to monitor the amount of leaking gas. Since 2015, the current monitoring at Hartoušov has been extended to other places in order to provide robust measurements capable of recording possible future gas anomalies at multiple sites. Because the conditions differed among the monitored sites, different measurement methods were designed. In this study, we distinguish between direct and indirect gas flow measurement methods (Camarda et al., 2006). The direct methods directly record the volume of gas per minute and require that gas flow be captured by a funnel or borehole. The indirect methods either involve deriving gas flow from the bubble fraction in water (pressure probes are placed beneath the water table) or rely on measuring gas overpressure in a closed borehole, or, finally, they calculate CO₂ flux from the concentration of gradients in the soil (Baubron et al., 1990). The dynamic concentration method is based on measuring the CO₂ content in a mixture of soil gas and air obtained by a special probe placed vertically in the soil. The dynamic concentration is proportional to the soil CO₂ flux according to an empirical relationship, which depends on soil permeability (Gurrieri and Valenza, 1988).

Table 1 CarbonNet monitoring network.

2.1 Monitoring network

Five gas escape sites were monitored in the period described: Hartoušov, Bublák, Soos, Dolní Částkov, and Prameny (see Table 1, the map in Fig. 1 and photos in Fig. 2). While the first three are located in mofette fields, the remaining are boreholes which tap mineral spring sources.

Figure 2 Photos of selected CO₂ monitoring fields – Hartoušov, SOOS, Bublák, and Prameny.

The pilot site of Hartoušov is located in a wooden hut above a 28.2 m deep F1 borehole, which taps a CO₂-saturated, pressurized aquifer. The plastic borehole casing, with an inner diameter of 115 mm, is perforated in the depth range of 20–28 m. Water level measurements date back to 2007, and gas flow has been measured here using a drum chamber gas flowmeter since 2009. The sensitivity of this type of instrument to environmental conditions (freezing or evaporation of the working liquid) caused gaps in the recorded time series. Since 2013, there have only been brief gaps thanks to the use of a different type of working liquid, improvements in the condenser separation, and thermal insulation. This direct field gas flow measurement is used as a reference for testing different flow measurement devices prior to their installation at other sites. Additional permanent measurements include water pressure in several depth levels, water temperature, and air temperature and pressure. In 2016, a 108.5 m deep borehole F2 was drilled in the Hartoušov mofette with the aim of studying geo–bio interactions (Bussert et al., 2017). It showed a CO₂ overpressure of 5 bar and was converted to a closed monitoring borehole with continuous measurements of downhole pressure and temperature, and wellhead pressure. A broadband seismometer was installed in a depth of 70 m in the year 2019.

The Bublák station has been located in a natural mofette in a swamp since 2015. To avoid interfering with natural conditions at this site, the equipment is buried underground, which does not allow for direct gas flow measurement. Instead, the differential water level is measured and used as a proxy for the volumetric fraction of free gas in water; see Sect. 2.3. Because of the rising bubbles, the water does not freeze in winter, making this measurement quite stable. The Soos station has been located in a natural mofette field since 2015, and the gas from a single mofette is captured by a funnel allowing for direct gas flow measurement. The small size of the metal box shelter and the need to use battery power, however, do not make it possible to prevent the freezing of the system in winter. The water level and temperature in the mofette and the volumetric fraction of free gas are measured here using an electric resistivity probe. In Dolní Částkov, the gas escapes both through a shallow borehole and the surrounding soil, which makes the flow measurements rather unstable. The Prameny station is located on top of a 100 m deep closed borehole (HJ-3A, drilled 1994) with degassing mineral water. Conditions at this site allow only for the measurement of the water level, temperature and wellhead gas pressure, which have been available since 2009.

2.2 Direct CO₂ flow measurement methods

Long-term gas flow monitoring in the field must meet various requirements. It should provide sufficiently accurate data of gas flow, which may contain dirt particles and moisture in changing field conditions of temperature, humidity, and air pressure. The presence of carbon dioxide further creates a highly corrosive environment, which the sensors should withstand. Commercial flowmeters are usually not designed to meet these demands. We have tested (at SOOS and Dolní Částkov stations) the MEMS (Micro-electro-mechanical systems) flowmeter, which is based on heat convection in moving gas. It works on the principle of a Wheatstone bridge, where changes in the resistivity of the resistor are measured according to the temperature changes caused by the flow of gas through a heater placed in the middle of the sensor (Dmytriw et al., 2007). These low-cost sensors, however, failed in our tests. None of the MEMS flowmeters tested measured for longer than 4 months, despite the installation of filters to capture solid particles and moisture from the gas before entering the sensor. A popular way of measuring gas flow is the Venturi-type flowmeter, which works by measuring the drop-in pressure at a constriction in a tube. Our tests of similar devices failed due to temperature drifts of the sensor and electronics, which were of the same order as the CO₂ flow variations. Direct flow measurement methods also include the acoustic method based on the Doppler effect, which is commonly used for water flow measurement. This, however, does not appear to be suitable for gas, which contains fewer particles acting as diffractors than liquid.

The standard flowmeters with rotating mechanical parts driven directly by the gas flow were also found not suitable due to the corrosive CO₂ environment. Better performance was achieved with a drum-type chamber flowmeter, which contains a revolving measuring drum within a packing liquid (we use low-viscosity oil). The measuring drum compulsorily measures volume by periodically filling and emptying four rigid measuring chambers. This chamber laboratory instrument was found suitable for field measurement, where sufficient space and nonfreezing temperatures can be guaranteed. It has been used as the primary flowmeter at the Hartoušov F1 borehole.

2.3 Indirect CO₂ flow measurement methods

2.3.1 Gas pressure in a closed borehole

In a closed borehole tapping a gas-saturated aquifer, an overpressure builds up whose magnitude has been speculated to reflect the amount of gas entering the aquifer from below (e.g., Hron and Škuthan, 2006). However, a profound discussion of exactly how the deep CO₂ leakage affects the measured overpressure is still absent, to the best of our knowledge. Considering a CO₂ flux $q=q_{\mathrm{1}}+q_{\mathrm{2}}$ summing the flux through the top of the aquifer in the vicinity of the borehole, q₁, and the possible gas leakage through the borehole, q₂, and assuming simply that the borehole overpressure p is proportional to the both, then p follows the equation

$$p={\displaystyle \frac{q}{K_{\mathrm{1}}+K_{\mathrm{2}}}},$$

where K₁ and K₂ are the permeability factors related to the ceiling of the aquifer and to the borehole sealing, respectively. Hence, the measured overpressure is not only proportional to the gas flux controlled by deep processes but also influenced by the permeability of the superficial layer as well as by any possible leaks through the wellhead. In particular, any variation in sealing layer properties, caused, e.g., by the actual weather conditions, is then directly projected onto the pressure measured.

Accordingly, in spite of the easy implementation of the pressure measurements in a closed borehole, we used this method only at the Prameny site, where technical and logistic conditions did not allow the installation and maintenance of a flow measurement. Excessive influence of K₁ on the measured pressure can be suppressed by introducing a controlled leakage in the wellhead, which ensures that K₂ is not small in comparison to K₁ as has been implemented at Prameny station.

2.3.2 Bubble fraction in water

We have used the bubble fraction monitoring method since observing a striking coincidence between the gas flow rate and groundwater level (see later in this section) increase in the Hartoušov F1 borehole during the 2014 seismic sequence (Fischer et al., 2017). Within a few months after the beginning of the sequence, the gas flow rate in the borehole increased 5-fold and the measured water level by more than 1 m. Since then, both quantities have indicated an overall gradual decrease back to their original levels.

Instead of the notion of groundwater level, adopted in Fischer et al. (2017) and other works, we stick to the more strictly defined terms pressure head and hydraulic head in the present paper, which is due a few explanatory comments. Within a steady water column resting in a borehole (or a narrow mofette), the hydraulic head (defined as the sum of the pressure head and elevation) is independent of the elevation and is referred to as the groundwater level, as it coincides with the elevation of the free water surface observed in the borehole. This is why the exact elevation of the actual placement of the pressure probes in the borehole is usually disregarded, and the term groundwater level is used somewhat loosely without risk of any confusion. In the case of a continuous bubbly flow through the borehole, however, hydraulic head is not a depth-independent quantity but rather inevitably increases with elevation. An intuitive explanation is that the mean density of the water–gas bubbles mixture is markedly lower than that of water (this is, however, merely an approximation; see also Sect. 2.5). Following this simple concept, the density of the mixture would be (disregarding the density of the gas CO₂ as negligible)

$$\begin{array}{}\text{(1)}& \mathit{\varrho }\left(z\right)=(\mathrm{1}-\mathit{\varphi }(z\left)\right){\mathit{\varrho }}_{\mathrm{w}},\end{array}$$

where ϕ(z) denotes the volumetric fraction of bubbles in the water column profile at elevation z, and ρ_w stands for the mass density of the water in the well (say, clear water at a given constant temperature). Denoting by ψ(z) the pressure head, related to the actual pressure p(z) through $p\left(z\right)={\mathit{\rho }}_{\mathrm{w}}g\mathit{\psi }\left(z\right)+b$, with g being the gravitational acceleration and b the barometric pressure (Fig. 3a), and denoting by $h\left(z\right)=\mathit{\psi }\left(z\right)+z$ the hydraulic head, we assert that the difference in the measured pressures, pressure head and hydraulic head, equals

$$\begin{array}{}\text{(2)}& \begin{array}{rl}p\left(z_{\mathrm{1}}\right)-p\left(z_{\mathrm{2}}\right)& =p_{\mathrm{1}}-p_{\mathrm{2}}=\underset{z_{\mathrm{1}}}{\overset{z_{\mathrm{2}}}{\int }}\mathit{\rho }\left(z\right)g\mathrm{d}z\\ & \Rightarrow {\mathit{\psi }}_{\mathrm{1}}-{\mathit{\psi }}_{\mathrm{2}}=\underset{z_{\mathrm{1}}}{\overset{z_{\mathrm{2}}}{\int }}\left(\mathrm{1}-\mathit{\varphi }\left(z\right)\right)\mathrm{d}z\\ & \Rightarrow h_{\mathrm{2}}-h_{\mathrm{1}}=\underset{z_{\mathrm{1}}}{\overset{z_{\mathrm{2}}}{\int }}\mathit{\varphi }\left(z\right)\mathrm{d}z.\end{array}\end{array}$$

That is, the hydraulic head h(z) measured in the borehole increases with elevation by a factor equal to the volumetric fraction of the bubbles in the borehole profile. Here, for the sake of brevity, we abstract from the time dependence of all quantities.

Figure 3 Measurement of pressure in the borehole with ascending bubbles: the measured pressure p and the related pressure head ψ at two different depths d₁ and d₂ within the bubble column using the (a) differential method; pressure within the bubble column p_m and beneath p_e used to determine the mean bubble fraction using the (b) integral method. Note that the difference in altitudes $z_{\mathrm{2}}-z_{\mathrm{1}}=d_{\mathrm{e}}-d_{\mathrm{m}}$.

The mean bubble fraction within the measured section of the water column can thus be defined as

$$\begin{array}{}\text{(3)}& {\mathit{\varphi }}_{\mathrm{12}}={\displaystyle \frac{h_{\mathrm{2}}-h_{\mathrm{1}}}{z_{\mathrm{2}}-z_{\mathrm{1}}}}=\mathrm{1}-{\displaystyle \frac{{\mathit{\psi }}_{\mathrm{1}}-{\mathit{\psi }}_{\mathrm{2}}}{z_{\mathrm{2}}-z_{\mathrm{1}}}}.\end{array}$$

As the ascending gas bubbles expand due to the decreasing pressure, both the volumetric flux of the gas and the bubble fraction ϕ(z) increase correspondingly with elevation. In order to obtain a quantity independent of the depths of the pressure probes, a further correction needs to be applied. A reasonable approximation can be obtained based on the following simplification. We assume that the gas expands isothermally, so that its volumetric flux is inversely proportional to pressure, and that the bubble fraction is approximately proportional to the volumetric flux, so that we can write

$$\mathit{\varphi }\left(z\right)={\mathit{\varphi }}_{\mathrm{0}}{\displaystyle \frac{p_{\mathrm{0}}}{p\left(z\right)}},$$

where ϕ₀ represents the bubble fraction at the reference pressure p₀ (which we later set as 100 kPa). Further, approximating the pressure profile between the two pressure probes by a linear function

$$p\left(z\right)=p_{\mathrm{1}}+{\displaystyle \frac{z-z_{\mathrm{1}}}{z_{\mathrm{2}}-z_{\mathrm{1}}}}(p_{\mathrm{2}}-p_{\mathrm{1}})$$

we obtain (by substituting ϕ(z) into Eq. (2) and integrating) the formula for ϕ₀; let us call this the projected bubble fraction,

$$\begin{array}{}\text{(4)}& {\mathit{\varphi }}_{\mathrm{0}}={\mathit{\varphi }}_{\mathrm{12}}{\displaystyle \frac{p_{\mathrm{2}}-p_{\mathrm{1}}}{p_{\mathrm{0}}\ln \left(\frac{p_{\mathrm{2}}}{p_{\mathrm{1}}}\right)}}.\end{array}$$

One should note that the quantity obtained here is subject to some uncertainty due to a number of simplifications and that it only gives the approximated volumetric fraction and not the gas flow rate itself (see Sect. 2.5).

In the Hartoušov F1 borehole, the pressure head had been measured until September 2018 by one pressure gauge in the depth of 8 m, well above the bubble entry point. In Fischer et al. (2017), the corresponding hydraulic head was referred to as the groundwater level. As proposed in the paper, we split its time variation into two parts: the variation (a) in the hydraulic head h_e(t) at the bottom of the bubble flow column and (b) in its increase through the column due to the gas bubbles. The optimal solution to obtain data for (a), which was implemented in Hartoušov F1 in late 2018, is to measure the pressure head in a depth beneath the bubble entry point directly by a dedicated pressure probe (Fig. 3b). While direct measurement was unavailable, it was supposed that (a) is given only by the surrounding hydrogeologic situation and is unaffected by the gas flow. A single measurement of the pressure at the bubble occurrence depth by Fischer et al. (2017), corrected by a continuous pressure head record from a nearby observation well in Hrzín 8 km away, which is not affected by the CO₂ gas flow, was used as h_e(t). Note that h_e(t) also describes the hydraulic head in any depth beneath the occurrence of bubbles. While Fischer et al. (2017) considered the possibility that the gas exsolution depth varies with time, we argue here (see Sect. 2.4) that the gas bubbles have to appear at the penetrated section of the Hartoušov F1 borehole. This allows us to determine the mean volumetric fraction of the bubbles using Eq. (3) with h₁(t)=h_m(t) being the hydraulic head measured at the depth d_m=4 m and h₂(t)=h_e(t) being the hydraulic head measured at the bubble entry depth (or anywhere below), which we suppose to be at the upper part of the penetrated section at d_e=20.5 m (Fig. 3).

In Fig. 4 the record of h_e(t) and h_m(t) and the resulting projected bubble fraction ϕ₀(t) defined by Eq. (4) are shown for the whole period studied in Fischer et al. (2017). We refer to this method as the integral method.

Figure 4 Daily averages of the reference pressure head h_e (grey) unaffected by gas flow and in the monitored Hartoušov F1 well, h_m (green); the volumetric fraction of bubbles (black) determined by Eq. (4) and CO₂ flow (blue) in the Hartoušov F1 well for the period 2008–2016. Because no reference pressure measurement below the bubble entry point was available for the studied period, a pressure head in nearby Hrzín borehole, which was free of gas flow, was used as h_e.

The method presented above is applicable only in boreholes and narrow tube-like mofettes. The borehole should tap the underground water, and there should exist a continuous column of gas bubble flow from a certain depth to the surface. Also, independent measurement of the hydraulic head in the aquifer/reservoir beneath the bubble flow column should be possible, either in the same well or, at least, in a nearby well free of gas flow. These conditions are not fulfilled in natural mofettes, which are usually only less than 2 m deep and communicate with the surface water significantly. In such cases, the difference in pressure heads along a fixed depth interval within the bubble flow column can be measured and used to define the mean bubble fraction ϕ₁₂(t) and the projected bubble fraction ϕ₀(t).

This differential method has been tested in the Hartoušov F1 well and Bublák mofette stations since 2015 using two analog water-level sensors attached at a 1 m distance on a metal rod. The obvious disadvantage is that both measurements (instead of just one) are subject to fluctuation due to the bubbly flow and that the noise in the resulting bubble fraction data is inversely proportional to the distance between the probes. To suppress the noise, an RC circuit with a 100 s time constant is applied.

An alternative way of determining the bubble fraction is based on the electric resistivity measurement of the water–bubble mixture. Unlike the pressure difference method, this method does not need to be focused on the vertical chain of bubbles, but it can assess the fraction of bubbles in a 3D volume defined by the geometry of electrodes. For this purpose, two water resistivities are measured by a special probe in the mofette: the reference resistivity of the water free of bubbles R_R and the resistivity of the water–bubble mixture R_M. The bubble fraction is then derived as

$$\mathit{\varphi }\left(t\right)=\mathrm{1}-c{\displaystyle \frac{R_{\mathrm{R}}\left(t\right)}{R_{\mathrm{M}}\left(t\right)}},$$

where c is the geometric calibration constant. This type of measurement has been tested in the Soos mofette since 2015.

2.4 Depth of gas bubbles appearance

It is possible to speculate that the exsolution of the gas bubbles from the water with dissolved CO₂ takes place at a certain depth in the borehole, while below that depth the pressure is sufficient to contain the CO₂ in the dissolved phase. In this view, the exsolution depth d_e could vary in time, as considered by Fischer et al. (2017), following variations in h_e(t) and in the CO₂ supply from the reservoir. Let us note, however, that such a scenario is only possible for gas fluxes much lower than those observed in the Hartoušov F1 borehole or in the presence of significant water discharge, such as in mineral springs.

Assuming a steady flow of the dissolved CO₂ up through a borehole section with no penetration below d_e, two transport mechanisms can be considered, convection or molecular diffusion. As for convection, no driving force to induce a flow in a water column in a borehole, in particular no significant temperature variations, has been observed in the Hartoušov F1 borehole. The mass flux due to molecular diffusion, on the other hand, can be estimated as follows, and it appears to be very limited. Assuming that the concentration of the dissolved CO₂ in the resting water column increases with increasing depth as much as allowed by the increasing hydrostatic pressure (with the Henry's law constant being on the order of 10⁻⁵ kg m⁻³ Pa⁻¹; see Sander, 2015), then the corresponding diffusive flow rate (with the diffusivity being on the order of 10⁻⁹ m² s⁻¹) through the borehole of the given cross-section area (say, 10⁻² m²) would be no more than on the order 10⁻¹² kg s⁻¹. This is 8 orders smaller than the flow rate directly observed in the borehole. We thus infer that the gas bubbles enter the Hartoušov F1 borehole in its penetrated part, as we assumed in the previous section, while they originate somewhere in the surrounding media.

2.5 Tests of bubble fraction method

The methodology for indirect gas flow measurement using pressure difference (the differential method) was first tested in the laboratory. The experimental setup consisted of an air pump connected with a valve for controlling the airflow, which was led to the bottom of a plastic tube with an inner diameter of 10.5 cm and a height of 2.5 m simulating the borehole. Two water pressure probes were installed at a fixed distance of 0.5 m on a vertical rod inside the tube, and all the air from this tube was led to the chamber gas flowmeter. The air inflow was increased stepwise, and at each level the data were recorded for a period of 15 min using a 1 Hz sampling rate. The gas flow ranged from 0 to 30 L min⁻¹, which corresponds to the volumetric flux ranging from 0 to 0.06 m s⁻¹. The observed mean bubble fraction appears to increase nonlinearly with the gas flow (Fig. 5). Note that the modification using Eq. (4) is insignificant here, due to the fact that both pressure probes are at a depth of less than 1 m. The bubble fraction values are more scattered than the gas flow rate measured by the flowmeter. The resulting noise was partially suppressed by low-pass filtering of the pressure data using an RC circuit with a time constant of 30 s and additionally by 1 min data sampling to smooth the water level values.

Figure 5 Comparison of gas flow velocity and volumetric fraction of bubbles for the field (Hartoušov mofette) and laboratory measurement. Laboratory measurements are smoothed by RC circuit of 30 s time constant (red points) and additionally by running average of 1 min length; the field measurements from 2016 (grey points) are based on the differential method with sensor distance of 1.0 m; the field measurements from 2018 to 2019 (blue points) are based on the integral method. Blue line is the best fit based on Eq. (5).

It is worth noting that the dynamics of bubbly flow in a borehole is quite a complex issue, which, however, appears to have been studied rather intensively in the chemical and nuclear engineering literature (see, e.g., Ghiaasiaan, 2007; Montoya et al., 2016). The simple considerations introduced in the previous text (Eqs. 3 and 4) correspond to the drift-flux model for a vertical borehole, provided that the water flux through the borehole is negligible. In particular, any momentum exchange with the walls is ignored. While this approach is well justified for the bubbly flow regime observed with small gas fluxes, with increasing volumetric gas one observes different flow regimes of greater complexity, such as the slug flow and churn flow. As the bubbles ascend, they increase in volume, join each other merrily, or even sadly split apart, their brief life being eventually cut short by obstacles such as the pressure probes dangling in the well; these are, however, out of the scope of this paper. Even in the bubbly flow regime, the relation between the bubble fraction ϕ and the volumetric flux of the gas bubbles u (m s⁻¹) has been described, e.g., by the following well-established empirical relation (Zuber and Findlay, 1965).

$$\begin{array}{}\text{(5)}& \mathit{\varphi }={\displaystyle \frac{u}{C_{\mathrm{0}}u+V_{gj}}},\end{array}$$

where C₀=1.13 and, assuming that the density of the gas is negligible when compared to that of water,

$$V_{gj}=\mathrm{1.41}{\left({\displaystyle \frac{\mathit{\sigma }g}{{\mathit{\rho }}_{\mathrm{w}}}}\right)}^{\mathrm{1}/\mathrm{4}}.$$

The curve in Fig. 5 is obtained by taking the surface tension for water and air at the laboratory temperature as σ=0.07 N m⁻¹. It appears that the mean bubble fractions derived from our pressure probe data overestimate the void fractions given by the Zuber–Findlay model, in particular for low flow rates.

For comparison, we also show in Fig. 5 the projected bubble fraction data ϕ₀ plotted against the corresponding gas flow rates measured on the Hartoušov F1 borehole. The comparison to the laboratory data and to the Zuber–Findlay model reveals a discrepancy that indicates some need for further analysis, which is beyond the scope of the present study. Let us briefly note that the difference cannot be explained by the mere parametric differences from the laboratory setting such as the temperature, gas, and water composition.

In Fig. 6 we compare the bubble fractions obtained at the Hartoušov F1 borehole using integral and differential methods. While it appears that the projected bubble fraction data from the Hartoušov F1 site cannot be directly inverted to obtain a reasonable gas flow rate estimate, it is important that they provide a fair correlation (see also Fischer et al., 2017, and Fig. 4 in their paper) and can thus provide a valuable gas flow rate proxy. As expected, the integral method data seem to perform better than the differential data (see Sect. 2.3). Note particularly the high noise of the latter and its lower correlation to the gas flow rate measurement (Fig. 6). Accordingly, using the pressure sensors at a larger distance, and, if possible, placing one of them below the bubble entry depth, seem preferable for indirect gas flow measurement.

Figure 6 (a) Comparison of volumetric bubble fraction derived using differential (top) and integral (bottom) methods with the gas flow at the Hartoušov F1 borehole. Depth of pressure probes for the differential method are 4.45 and 5.45 m below the surface. For the integral method it is 5.45 and 27.2 m. In (b) the cross-correlation functions between the bubble fraction and gas flow rate are shown.

2.6 Environmental effects

The measurements of CO₂ flow, CO₂ pressure, and pressure head are influenced by environmental effects – mainly variations in temperature (diurnal, seasonal), changes in barometric pressure, and tidal effects. Temperature and barometric changes are the most significant, since their influence can be local and can vary even among the stations of the network. Barometric and tidal loading of aquifers has been studied in detail (e.g., Jacob, 1939; Rojstaczer and Riley, 1990; Roeloffs, 1996). Here, we address the basic principles that are relevant to the pressure and production of the up-streaming gas. Both the confined and unconfined response of pressure head are characterized by the barometric efficiency E_B, which expresses the ratio of the change in the hydraulic head Δh caused by the barometric pressure change Δb,

$$\begin{array}{}\text{(6)}& E_{\mathrm{B}}={\mathit{\rho }}_{\mathrm{W}}g{\displaystyle \frac{\mathrm{\Delta }h}{\mathrm{\Delta }b}}\end{array}$$

The net response is always a decrease in the hydraulic head with an increase in barometric pressure. The barometric pressure variations act directly on the open water level in the well and also on the formation composed of the mineral matrix and pore space filled by the water. As a result, the direct effect on the water level in the borehole is partially suppressed by the fraction of the external load borne by the formation water. Hence, the barometric efficiency can also be written as

$$\begin{array}{}\text{(7)}& E_{\mathrm{B}}={\displaystyle \frac{\mathit{\theta }\mathit{\beta }}{\mathit{\theta }\mathit{\beta }+\mathit{\alpha }}},\end{array}$$

where α and β are the compressibilities of the rock matrix and water, respectively, and θ represents porosity within the aquifer. Thus, barometric efficiency can be described as the fraction of specific storage derived from the compressibility of water or, equally, as the fraction of external load change borne by the formation either as compaction or expansion. Accounting for the range of fractured rock compressibilities, E_B of confined aquifers usually ranges between 0.2 and 0.7 (Todd, 1980) and may reach 1.0 for granite with a very low compressibility of the rock matrix (Roeloffs, 1996; Acworth and Brain, 2008). Note, however, that large values of E_B may also correspond to large values of β; this fact is not addressed in the literature for the simple reason that it is usually the rock that varies from site to site and not the water. In this concern, the possible effect of the presence of the compressible CO₂ bubbles within the aquifer surrounding the borehole on the barometric efficiency is a question that has not been addressed in the literature, to the best of our knowledge.

Similarly, the barometric effect to the CO₂ discharge from an aquifer through an open well has not been studied either. One can expect that an increase in the pore pressure due to an increase in the barometric pressure allows for larger amounts of CO₂ to be dissolved in water, which in turn decreases the volume of CO₂ leaking into the well. Similarly, a decrease in barometric pressure may induce increased degassing. Note that there exist many unknowns in this regard, such as the flow paths of the gas ascending through the aquifer, the amount of the mobile and immobile gas bubbles in the porous space, etc. In the Hartoušov F1 borehole, a strong anticorrelation between the gas flow rate and the barometric pressure has been observed.

We correct the measured quantity f (pressure head or gas flow) for demeaned barometric pressure variations b using the equation $f_{\mathrm{c}}=f-E_{\mathrm{B}}$ (b−〈b〉). Barometric efficiency E_B is determined with the target of minimum cross-correlation of f and f_c. To account for the possible time variation in E_B a sliding window of 1 d is used; see Fig. 7a for original and corrected records of pressure head and gas flow in Hartoušov F1. Figure 7b shows the cross-correlation functions between barometric pressure and original and corrected records. The success of barometric correction is indicated both by removing the anticorrelation with air pressure and by minimizing short-period variations in the corrected records. The mean barometric efficiency was 0.76 for the pressure head and 0.46 L min⁻¹ kPa⁻¹ for the gas flow.

Figure 7 The barometric effect to the pressure head and CO₂ flow in the Hartoušov F1 borehole for the period from October 2018 to April 2019. (a) Original measurements are indicated in blue, and those corrected for barometric pressure are in red. The upper panel shows pressure head at the depth below the bubble formation, and the lower panel shows gas flow measured by flowmeter. The success of barometric correction is illustrated in (b) showing the decrease in barometric anticorrelation after correcting.

Other external effects like diurnal temperature variations and Earth tides were found to be much weaker than the influence of barometric pressure. The volumetric fraction of bubbles is not affected by air temperature, since the sensors are placed in groundwater with an almost constant temperature. In addition, the periods of diurnal temperature variations and significant Earth tide components are significantly shorter than the expected durations of anomalies of deep-generated gas flow. Accordingly, we do not apply corrections for temperature variations and Earth tides.

3 Results of CO₂ flow monitoring in West Bohemia

The time series of gas production at all monitored stations, along with seismicity plot, are shown in Fig. 8. The record at Hartoušov F1 for the period from late 2007 to 2019 shows a long-term decrease interrupted by several abrupt massive increases in gas discharge. The maximum flow, reaching 50 L min⁻¹, followed the 2014 seismic sequence in late summer/autumn 2014; the minimum values, below 10 L min⁻¹, were observed prior to the 2014 seismic sequence and at the present time. The fast coseismic increase and long-term, post-seismic decrease are visible both in the gas flow and integral bubble fraction data determined using Eqs. (3) and (4) and are consistent with Sibson's fault valve model (Fischer et al., 2017). Note particularly the abrupt rise in gas flow and CO₂ bubble fraction during the M_L 4.4 seismic sequence of May–August 2014 and in bubble fraction during the October 2008 M_L 3.8 swarm. Next to these striking coincidences of seismic activity and CO₂ release we also find cases of strong seismic activity, which was not accompanied by a significant gas flow anomaly (see the M_L 3.4 swarm of 2011 and the most recent M_L 3.8 swarm of 2018). On the other hand, the CO₂ flow record shows a few positive pulses which are not related to significant seismic activity (Fig. 8b). The most striking one is the gas production increase in the period from the beginning of May till the end of July 2016, which is visible both in the gas flow and bubble fraction records. This is, however, undoubtedly of anthropogenic origin caused by drilling of the nearby 108.5 m deep F2 borehole at a distance of 40 m from the monitored F1 borehole; drilling started on 30 March (Bussert et al., 2017). The drilling reached the top of a CO₂ pressured horizon at a depth of 80 m on 21 April and created a shortcut to the shallower aquifer, which was tapped by the monitored borehole. The 3-month-long gas increase thus represents a delayed response to a nearby drill. Another, less pronounced, positive pulse in the period from mid-September to late November 2016 is of unknown origin. A number of negative pulses and oscillations are found on the bubble fraction record alone, which lower the correspondence between the gas flow rate and the bubble fraction data and indicate a more complex relation between gas flow in a borehole and volume fraction of ascending bubbles, as already noted in Sect. 2.5.

Figure 8 Comparison of seismic activity and CO₂ production at individual monitoring sites in West Bohemia. For Hartoušov F1 the CO₂ flow and gas bubble fraction are shown; for Bublák and Soos the CO₂ bubble fraction is shown, and for Prameny the gas pressure in a closed borehole is plotted; (a) period 2007–2019; (b) detail for the period 2014–2019. Gas bubble fraction was determined using the pressure difference method in the integral version for Hartoušov (blue), in differential version for Bublák (red), and by the resistivity method for Soos (light brown).

The records of gas differential bubble fraction data in Bublák and resistivity-based bubble fraction in Soos indicate in the monitored period since autumn 2015 a steady gas release with only a few bumps, which are most probably of local origin and related to the shallow character of the mofettes. Gas at these sites passes through approximately cylindrical vents of ∼0.5 m diameter and ∼1 m depth filled by surface water. The similarity of bubble fraction increase at Soos and gas flow increase at Hartoušov F1 in summer 2016 is most probably merely accidental, considering the anthropogenic origin of the rise at Hartoušov and the large mutual distance of about 5 km of these sites.

As mentioned in Sect. 2.3, the integral method of bubble fraction measurement provides better results than the differential method. The latter suffers particularly from high noise caused by the placement of both pressure probes in a water column with flowing bubbles as shown in Fig. 6. One can also notice a better correlation of the integral method compared to the differential one. Unfortunately, due to technical problems, we were not able to perform this comparison for the same time window – so time windows of the same length (3 months) free of any technical issues were selected.

4 Discussion

The barometric efficiency E_B of the groundwater pressure head of 0.76, which we obtained, is relatively high. The high values of E_B are generally considered an indication of the small compressibility of the rock matrix that is typical for unweathered granite (Acworth and Brain, 2008). The target aquifer is formed by sedimentary formations of the Cheb Basin composed of sandstones and conglomerates with varying clay contents underlain by mica-schist basement (Bussert et al., 2017). The compressibility of these types of rocks is, however, 3 to 6 times greater than of granite (their bulk moduli range from 10 to 20 GPa compared to 50 GPa for granite). Using Eq. (4), porosity of 30 %, and bulk moduli ratio of matrix and pore fluid equal to 5, one gets E_B=0.5. The level E_B = 0.76 is reached for bulk moduli ratio of 15. Assuming the bulk modulus of aquifer rocks about 10 GPa, one obtains a bulk modulus of the fluid of only about 0.7 GPa, which corresponds to 3 times larger compressibility than for water. This could be explained by the presence of carbon dioxide in the groundwaters in gaseous phase and is worth further research.

The gas flow trend in Hartoušov after the 2014 seismic sequence shows signatures similar to those in the period before 2014, which followed the 2008 swarm. A similar, long-term overall decrease is followed by steady-state behavior with an almost constant flow rate of about 10 L min⁻¹. In terms of the Sibson's fault valve model, this corresponds to the self-sealing phase of the fault due to mineral precipitation (Sibson, 1992) when pressure builds up and in combined action with tectonic loading results in increasing instability of the fault. This inevitably leads to later recurrence of fault failure in the form of seismic activity and regeneration of fault permeability. As indicated above, the coincidence of a massive rise in CO₂ flow and seismic activity has not been observed since the 2014 seismic sequence. Indeed, none of the earthquake swarms since 2014 have been accompanied by a distinct CO₂ degassing anomaly (Fig. 8b). All in all, in the whole period of CO₂ flow monitoring in Hartoušov since 2007, five earthquake swarms with magnitude M_L larger than 3.0 occurred (2008, 2011, 2014, 2017, 2018), and only two of them (2008 and 2014) were accompanied by a strong and long-lasting coseismic increase in CO₂ degassing. This is not surprising in general, because the fault valve mechanism might act only under certain circumstances. And even if a fluid pulse is released during every stronger seismic sequence its volume might not be sufficient to reach the Earth's surface with a detectable amplitude. This is also directly related to the pressure buildup in the fluid reservoir beneath the sealed fault, which is a long-lasting process, and thus earthquakes that occur soon after releasing the accumulated fluid pressure are likely to not be accompanied by a significant fluid release.

In this context it is also of interest to consider the hypocenter cluster geometry in 3D and its relation to the presence of permeable channels in a shallow crust allowing crustal fluids to reach the surface. In Fig. 9 hypocenters of individual earthquake sequences are indicated in a vertically oriented cross section and show a hat-like structure in depths from 6.0 to 10.5 km extending about 10 km north–south. The Hartoušov mofette field is located about 10 km south from the center of the main cluster, which corresponds to 6 km distance from its southern tip. A pronounced segmentation of the fault plane is apparent with the 2008 and 2014 segments in the southern branch of the cluster and the 2011, 2017, and 2018 segments clusters in its northern branch. It is worth noting that the 2014 mainshocks showed unfavorable oriented focal mechanisms and occurred on a fault jog activated by stress concentration resulting from previous swarm activity (Hainzl et al., 2016; Jakoubková et al., 2018). This structural and possibly impermeable boundary within the fault zone was broken by the M_L 3.5 mainshock of the 2014 sequence – the first earthquake of this sequence, which was followed by the massive CO₂ flow rise in Hartoušov.

Recently, since the summer of 2019, the CO₂ flow rate in Hartoušov F1 has decreased below 10 L min⁻¹, which could be a sign of the approaching occurrence of a new seismic swarm, according to the Sibson's fault valve model. However, one should also take into account that the flow rate decrease in 2019 could have been caused by the drought period in the summer. The reduced groundwater pressure in the whole area would lead to the rise of the diffuse component of gas flow reducing the gas discharge in the borehole. Comparing the records of groundwater level in nearby mofettes with the gas flow rate in the F1 borehole, however, gives unequivocal results. While some correlation between gas flow and water level was found for 2018 and 2019, the gas flow rate in 2017 was found independent of the water level in mofettes.

The clear coseismic CO₂ flow rate increase during the 2008 and 2014 seismic sequences indicates the presence of a permeable channel between the southern cluster and the Hartoušov mofette field (Fig. 9). The absence of CO₂ flow anomalies coinciding with the seismic activity in northern clusters could be interpreted to show that the hydraulic connection between these fault patches and the Hartoušov mofette is missing, which could be related to the aforementioned fault jog. Besides, it is also of interest that epicenter distribution and CO₂ degassing occurrence is typically separated in the area (Weinlich et al., 2006; Babuška and Plomerová, 2008); most earthquakes occur in the northern, CO₂-free part of the Cheb Basin.

Figure 9 Vertical, north–south-oriented section through the Nový Kostel fault zone with hypocenters of the earthquake swarms occurred between 2008 and 2018. The position of the Hartoušov mofette showing coseismic CO₂ flow rate increase is indicated on the surface. Red ovals show the position of first events of the 2008 and 2014 seismic sequences.

Other monitored sites such as Bublák and Soos show, similar to Hartoušov, almost constant CO₂ discharge since early 2017. As these stations were not in operation during the 2008 and 2014 seismic sequences showing coseismic CO₂ increase in Hartoušov, no inferences about their correlation to the seismic activity can be drawn. Bubble fractions derived from resistivity measurements of the water–bubble system are found quite stable. For the measurement period of 4 years no maintenance of the probe was required, and, compared to the pressure difference method, the system is less sensitive to the depth of the probe below the water table. However, no seismogenic anomaly of as flow rate has occurred yet that could be used to calibrate the system. Continuous monitoring of CO₂ degassing is required to determine whether future seismic activity in the southern cluster will generate an increase in degassing in either of the monitored sites and enable the verification of the hypothesis that only earthquakes in the southern cluster are capable of generating a CO₂ pulse which reaches the surface.

5 Conclusions

The present study is focused on the long-term monitoring of CO₂ degassing in the form of mofettes and gaseous mineral springs targeted on the West Bohemia/Vogtland region in central Europe, which is typified by the occurrence of earthquake swarms and discharge of carbon dioxide of magmatic origin. The gas flow measurement is applied to two types of sources: natural wet mofettes with gas outflow through surface water pools and boreholes tapping shallow CO₂-saturated aquifers. The different local conditions of the five monitored sites call for different methods of gas capture and flow rate measurement. Besides the direct flow measurement using a drum chamber gas flowmeter, electronic MEMS flowmeters, and Venturi-based probes we introduce a novel, indirect method based on quantifying the gas bubble contents in a water column, which is capable of functioning in severe environmental conditions. The method is based on measuring the pressure difference along a fixed depth interval in a water column, which is proportional to the mean bubble fraction within the measured section. We analyze the dependence of the bubble fraction on depth and project it to the atmospheric pressure to make it directly comparable to the gas flow rate. Laboratory tests indicate the nonlinear dependence of the bubble fraction on the flow rate, which is confirmed by empirical models found in the chemical and nuclear engineering literature. Flow rates and bubble fractions observed in a pilot borehole F1 in the Hartoušov mofette show a high mutual correlation; however, some discrepancy is found between the measured flow rate and that predicted by the empirical models. This discrepancy calls for further analysis.

We also analyzed the long-term monitoring of gas flow and bubble fraction in the pilot borehole for the period 2008–2019. We found a quite strong barometric influence on the hydraulic head of the confined aquifer corresponding to a barometric efficiency of 0.76, which can be attributed to the compressibility of the pore fluids including the gaseous phase of carbon dioxide.

The record of gas flow rate and bubble fraction in Hartoušov F1 shows two high-amplitude coseismic rises coinciding with the occurrence of earthquake swarms in 2008 and 2014. The flow rate increased to a multitude of the preseismic level for several months and was followed by a long-term decay. However, another three seismic swarms occurring in the same fault zone were not associated with any significant CO₂ flow anomaly. We surmise that this may be related to the slightly farther location of hypocenters of these swarms in comparison with the two which caused the coseismic CO₂ flow rise. Further long-term CO₂-flow monitoring is required to verify the mutual influence of CO₂ degassing and seismic activity in the area.