On the link between Earth tides and volcanic degassing

Long-term measurements of volcanic gas emissions conducted during the last decade suggest that under certain conditions the magnitude or chemical composition of volcanic emissions exhibits periodic variations with a period of about 2 weeks. A possible cause of such a periodicity can be attributed to the Earth tidal potential. The phenomenology of such a link has been debated for long, but no quantitative model has yet been proposed. The aim of this paper is to elucidate whether a causal link between tidal forcing and variations in volcanic degassing can be traced analytically. We model the response of a simplified magmatic system to the local tidal gravity variations and derive a periodical vertical magma displacement in the conduit with an amplitude of 0.1–1 m, depending on the geometry and physical state of the magmatic system. We find that while the tide-induced vertical magma displacement presumably has no significant direct effect on the volatile solubility, the differential magma flow across the radial conduit profile may result in a significant increase in the bubble coalescence rate at a depth of several kilometres by up to several multiples of 10 %. Because bubble coalescence facilitates separation of gas from magma and thus enhances volatile degassing, we argue that the derived tidal variation may propagate to a manifestation of varying volcanic degassing behaviour. The presented model provides a first basic framework which establishes an analytical understanding of the link between the Earth tides and volcanic degassing.

1 Introduction

Residual gravitational forces of the Moon and the Sun deform the Earth's surface and interior periodically and thus lead to the so-called Earth tides. The tidal potential can be modelled as the result of the interference of an infinite number of sinusoidal tidal harmonics with precisely known frequencies and amplitudes (Darwin, 1883; Doodson, 1921). At the Equator, the tidal potential varies predominantly with a semi-diurnal periodicity. The amplitude of the semi-diurnal cycle is modulated within the so-called spring–neap tide cycle with a periodicity of 14.8 d caused by the interference of the lunar semi-diurnal tide and the solar semi-diurnal tide. The peak-to-peak amplitude of the associated semi-diurnal gravity variations is $a_{\mathrm{astro}}^{\mathrm{st}}=\mathrm{2.4}{\mathrm{\mu }\mathrm{m}\mathrm{s}}^{-\mathrm{2}}$ during spring tide and $a_{\mathrm{astro}}^{\mathrm{nt}}=\mathrm{0.9}{\mathrm{\mu }\mathrm{m}\mathrm{s}}^{-\mathrm{2}}$ during neap tide and is at an intermediate level at other times of the cycle. At midlatitudes, the tidal potential varies predominantly with diurnal periodicity, and at other latitudes both periodicities mix. The spring–neap tide cycle is, however, manifested everywhere and has maximum variability at the Equator (Agnew, 2007). The tidal potential firstly gives rise to a periodical elevation of the Earth's crust with a semi-diurnal peak-to-peak variation of up to about 50 cm (maximum at the Equator), and secondly all crustal compartments exhibit an additional semi-diurnal gravity variation by up to $\mathrm{1.16}\cdot a_{\mathrm{astro}}^{\mathrm{st}}$ (Harrison et al., 1963; Baker, 1984). This gravity variation typically has no effect on the rigid solid crust but can cause fluid movement, e.g. prominently manifested in the form of ocean tides (Ponchaut et al., 2001).

Evidence for tidal impacts on volcanism has been gathered by numerous empirical studies, which detected a temporal proximity between tidal extrema and volcanic eruptions (Johnston and Mauk, 1972; Hamilton, 1973; Dzurisin, 1980) or seismic events (McNutt and Beavan, 1981, 1984; Ide et al., 2016; Petrosino et al., 2018) or found a correlation between the spring–neap tide cycle and variations in volcanic deformation (De Mendoca Dias, 1962; Berrino and Corrado, 1991) or variations in the volcanic gas emissions.

The tide-induced stress variations (∼ 0.1–10 kPa) appear to be negligibly small in comparison to tectonic stresses (∼ 1–100 MPa) or stresses caused by pressure and temperature gradients within a shallow magmatic system (∼1 MPa). The rate of tidal stress change can, however, be around 1 kPa h⁻¹ and thus potentially exceeds stress rates of the other processes by 1 to 2 orders of magnitude (Sparks, 1981; Emter, 1997; Sottili et al., 2007). Furthermore, these subtle stress variations may cause an amplified volcanic reaction, when, for example, the tidal variations cause a widening of tectonic structures (Patanè et al., 1994), a periodic decompression of the host rock (Sottili et al., 2007; Sottili and Palladino, 2012), a variation in the host rock permeability (Bower, 1983; Elkhoury et al., 2006; Manga et al., 2012), self-sealing of hydrothermal fractures (Cigolini et al., 2009), or a mechanical excitation of the uppermost magmatic gas phase (Girona et al., 2018).

First studies on the covariations in tidal patterns and volcanic gas emissions hypothesised a possible tidal impact on the observed sulfur dioxide (SO₂) emission fluxes at Masaya (Stoiber et al., 1986) and Kilauea (Connor et al., 1988). Since the 2000s, automatic scanning networks based on UV spectrometers (e.g. Galle et al., 2010) have provided multi-year time series of volcanic gas emissions of SO₂ and bromine monoxide (BrO). The availability of such data sets enabled extensive investigation of long-term degassing variations. Correlation with the long-term tidal patterns has been reported for the SO₂ emission fluxes of Villarrica and Llaima (Bredemeyer and Hansteen, 2014) and the BrO∕SO₂ molar ratios in the gas plume of Cotopaxi (Dinger et al., 2018). Another possible but less significant correlation has been reported for the SO₂ emission fluxes of Turrialba (with a periodicity somewhere between 9.1 and 16.7 d; Conde et al., 2014). Furthermore, Lopez et al. (2013) reported a periodicity of roughly 16 d in the SO₂ emission fluxes of Redoubt retrieved from the satellite-based Ozone Monitoring Instrument (OMI) – the authors proposed that this periodicity was, however, an artefact of the satellite orbit rather than a tidal signal. In addition, correlation with the long-term tidal patterns have been reported for the diffuse radon degassing of Terceira (Aumento, 2002) and Stromboli (Cigolini et al., 2009).

Cycles in volcanic degassing patterns are not unique to periodicities which match the tidal potential. Many studies reported periodic volcanic degassing patterns with periods of minutes (e.g. Fischer et al., 2002; Boichu et al., 2010; Campion et al., 2012, 2018; Tamburello et al., 2013; Pering et al., 2014; Ilanko et al., 2015; Moussallam et al., 2017; Bani et al., 2017). In contrast, observations of long-term periodicities are rare. Besides the above-mentioned, roughly biweekly periodicities, periodic long-term pattern with periodicities of 50 and 55 d have been observed in the SO₂ emission flux of Soufrière Hills (Nicholson et al., 2013) and Plosky Tolbachik (Telling et al., 2015), respectively.

In view of the growing number of studies revealing similar biweekly patterns in volcanic activity, this paper investigates whether causality between the tidal potential and variations in the volcanic degassing is analytically traceable in a comprehensible way. High-temperature gas emissions of persistently strong, passively degassing volcanic systems are commonly thought to be fed by sustained magma convection reaching the uppermost portions of the volcanic conduit, where volatile-rich low-viscosity magma ascends through essentially degassed magma of higher viscosity, which in turn descends at the outer annulus of the conduit (Kazahaya et al., 1994; Palma et al., 2011; Beckett et al., 2014). Magma ascent rates associated with such convective flow typically vary roughly between 1 and 100 m h⁻¹ (Cassidy et al., 2015, 2018) and thus are orders of magnitudes larger than what we can derive for potentially tide-induced vertical magma displacement rates of at most 0.6 m within 6 h (if not further amplified). A comprehensive model of the tidal impact on the magma motion thus requires a coupling of the convective and the tide-induced transport mechanisms.

Our conceptual model aims to provide the first step by investigating the purely tide-induced transport mechanism acting on the low-viscosity inner magma column, neglecting any interference between the magma ascent and the tidal mechanism, i.e. the model ignores the magma convection in the column. We model the response of such a quasi-static magmatic system (volcanic conduit connected to a laterally more extended deeper magma reservoir) to tide-induced gravity variations analogously to the response of a classical mercury thermometer to temperature variations: the tide drives a periodical expansion of the magma in the reservoir, which leads to a periodical vertical displacement of the low-viscosity magma column in the conduit.

We derive the temporal evolution and amplitude of the vertical magma displacement across the radial conduit profile and examine its impact on the bubble coalescence rate. In order to introduce our novel approach comprehensibly, the modelled processes and conditions are as simplified as suitable; the major simplifications are listed in Appendix A. All findings in this paper are derived analytically. The quantitative model estimates are presented for two exemplary magmatic systems. These examples are intended to match simplified versions of the Villarrica (39.5^∘ S) and Cotopaxi (0.7^∘ S) volcanoes, where covariation between outgassing activity and Earth tidal movements has been observed previously (Bredemeyer and Hansteen, 2014; Dinger et al., 2018). The associated model parameter sets are listed in Table 1. Further, all quantitative estimates are presented for the spring tide, and the consequences of the contrast between spring tide and neap tide are discussed in the last part of this paper.

2 Tide-induced magma displacement in the conduit

2.1 Model set-up

Figure 1 Panels (a) and (b): sketch of the model set-up. The model compartments are indicated by white boxes but not depicted to scale. (b) The semi-diurnal tide causes a radial magma displacement profile in the conduit with different amplitudes during spring tide and neap tide. (c) Concept of the tide-enhanced bubble coalescence: two bubbles which are initially close to each other (see “without tide”) exhibit differential vertical tide-induced displacements, which enhances the chance for bubble coalescence (here “at low tide”).

We model the magmatic system analogously to established convection models (Kazahaya et al., 1994; Palma et al., 2011; Beckett et al., 2014), with the exception that the descending high-viscosity magma annulus is assumed to be not affected by the tide-induced dynamics and therefore is considered as an effective part of the host rock, while “conduit” refers in our model exclusively to the ascending low-viscosity magma column. We assume the conduit to be a vertically oriented cylinder with length L_c, radius R_c, and cross-sectional area $A_{\mathrm{c}}=\mathit{\pi }\cdot R_{\mathrm{c}}^{\mathrm{2}}$, which is confined by the penetrated host rock (and high-viscosity magma annulus), connected to a deeper, laterally more extended magma reservoir with volume V_r and centre of mass at a depth D_r, and either exhibiting an open vent or capped by a gas-permeable solid plug (Fig. 1). The magmatic melt in the conduit is modelled as a mixture of a liquid phase and a gas phase having a mean density ρ_melt, which varies with pressure and thus depth, a constant kinematic bulk viscosity ν, and homogeneous local flow properties. The magma compressibility β(ϕ) strongly depends on the gas volume fraction ϕ and lies between the compressibility $\mathit{\beta }\left(\mathrm{0}\right)=\mathrm{2}\cdot {\mathrm{10}}^{-\mathrm{10}}{\mathrm{Pa}}^{-\mathrm{1}}$ of volatile-rich rock and the compressibility ${\mathit{\beta }}_{\mathrm{0}}\left(\mathrm{1}\right)\approx p^{-\mathrm{1}}$ of an ideal gas (see, e.g., Tripoli et al., 2016). The magmatic melt in the reservoir is modelled to be volatile-rich but hosting no gas phase of significant volume and thus having a constant compressibility β_r≈β(0). Further, the quasi-static condition implies a steady-state density stratification within the magma and also with respect to the host rock (no neutral buoyancy; Parfitt et al., 1993). In this equilibrium, we assume a constant hydrostatic pressure gradient (∇p)_vert.

2.2 Response of the host rock on tidal stresses

Magma pathways are often located at intersection points of large-scale fault systems (Nakamura, 1977; Takada, 1994) or in fault transfer zones (e.g. Gibbs, 1990), where the surrounding host rock geometry is relatively sensitive to directional changes in pressure. The vertical and horizontal components of the tidal force exert additive shear tension on the host rock, potentially causing a compression of the host rock (Sottili et al., 2007) or a differential slip between both sides of the fault system (Ide et al., 2016). Both mechanisms can cause an increase in the areal conduit cross section. Connected to the magma reservoir, such an increasing conduit volume is accompanied by decompression and thus causes magma to flow from the reservoir to the conduit, which pushes the initial magma column in the conduit upwards until the initial hydrostatic pressure gradient is re-established. Vice versa a relative decrease in the areal conduit cross section leads to an effective descent of the initial magma column in the conduit. For a given periodic area increase ΔA_c, the amplitude Δz_hr of this additive elevation-descent cycle of the centre of mass of the initial magma column is given by

$$\begin{array}{}\text{(1)}& {\displaystyle }\mathrm{\Delta }{z}_{\mathrm{hr}}={\displaystyle \frac{L_{\mathrm{c}}}{\mathrm{2}}}\cdot {\displaystyle \frac{\mathrm{\Delta }{A}_{\mathrm{c}}}{A_{\mathrm{c}}+\mathrm{\Delta }{A}_{\mathrm{c}}}}\approx {\displaystyle \frac{L_{\mathrm{c}}}{\mathrm{2}}}\cdot {\displaystyle \frac{\mathrm{\Delta }{A}_{\mathrm{c}}}{A_{\mathrm{c}}}}.\end{array}$$

The quantitative scale of tide-induced conduit cross section variations is presumably hardly accessible. The theoretical horizontal components of the tide-induced ground surface displacement are up to about ±7 cm (Baker, 1984). Slip-induced dilation of faults with widths in the sub-centimetre range thus appear to be plausible. For illustration, a conduit radius increase by ΔR_c=1 mm would result in an additive vertical centre of mass displacement by Δz_hr=0.33 m for Villarrica and Δz_hr=0.13 m for Cotopaxi. As a remark, these mechanisms do not require a cylindrical conduit and fault–slip mechanisms would lead to an unidirectional area increase rather than a homogeneous radial increase. Furthermore, the tide could also cause a variation in the host rock permeability (Bower, 1983; Elkhoury et al., 2006; Manga et al., 2012). This mechanism and its possible interference with the concept presented here are ignored in our model.

2.3 Tide-induced magma expansion in the reservoir

The semi-diurnal tide causes a sinusoidal variation in the gravitational acceleration with angular frequency ω_sd = 1.5 ⋅ 10⁻⁴ rad s⁻¹ and amplitude (equals the half peak-to-peak amplitude) $a_{\mathrm{0}}^{\mathrm{st}}$ = 1.4 µm s⁻² during spring tide and $a_{\mathrm{0}}^{\mathrm{nt}}=\mathrm{0.5}{\mathrm{\mu }\mathrm{m}\mathrm{s}}^{-\mathrm{2}}$ during neap tide. Besides those host rock mechanisms triggered by the tidal stresses, these tide-induced gravity variations may also cause a periodical elevation of the magma in the inner conduit.

The compressible magma in the reservoir is pressurised by the hydrostatic load whose weight is proportional to the local gravitational acceleration g. A reduction in the local gravitational acceleration by a₀ leads to a decompression and thus expansion of the magma in the reservoir by $\mathrm{\Delta }{V}_{\mathrm{r}}=\frac{a_{\mathrm{0}}}{g}\cdot (\mathrm{\nabla }p{)}_{\mathrm{vert}}\cdot D_{\mathrm{r}}\cdot {\mathit{\beta }}_{\mathrm{r}}\cdot V_{\mathrm{r}}$. The tidal force can accordingly lead to a periodical magma expansion–shrinkage cycle in the reservoir with a semi-diurnal periodicity and an amplitude modulation within the spring–neap tidal cycle of up to ΔV_r∼𝒪(100–1000 m³).

The realisation of this additional magma volume implies a displacement and thus compression of the host rock at the location of maximum host rock compressibility. This is typically true for the conduit. Assuming that the magma expansion in the reservoir ultimately and exclusively causes an increase in the conduit volume, the volume increase causes an elevation of the centre of mass of the initial magma column in the conduit by

$$\begin{array}{}\text{(2)}& {\displaystyle }\mathrm{\Delta }{z}_{\mathrm{dec}}={\displaystyle \frac{\mathrm{\Delta }{V}_{\mathrm{r}}}{A_{\mathrm{c}}}}={\displaystyle \frac{a_{\mathrm{0}}}{g}}\cdot (\mathrm{\nabla }p{)}_{\mathrm{vert}}\cdot D_{\mathrm{r}}\cdot {\mathit{\beta }}_{\mathrm{r}}\cdot {\displaystyle \frac{V_{\mathrm{r}}}{\mathit{\pi }\cdot R_{\mathrm{c}}^{\mathrm{2}}}}.\end{array}$$

In the general case, the additional volume could be realised by a slight increase in the conduit radius by $\mathrm{\Delta }{R}_{\mathrm{dec}}\approx \frac{R_{\mathrm{c}}}{\mathrm{2}}\cdot \frac{\mathrm{\Delta }{z}_{\mathrm{dec}}}{L_{\mathrm{c}}}\sim \mathcal{O}\left(\mathrm{1}\mathrm{mm}\right)$ caused, for example, by the tidal stresses. If the magmatic system has an open vent, the additional volume can alternatively be realised by an elevation of the lava lake level and thus without a host rock compression.

Analogously, the tide-induced gravity variations result in an expansion of the initial magma column in the conduit. This effect is, however, typically negligible compared to the reservoir effect for sufficiently large reservoirs (volume contrast between reservoir and conduit of more than 1000; see Table 1); thus, for simplicity we neglect the effect of the expansion of the initial magma column in the conduit.

The responses of the overall magmatic system on the tidal stresses and tide-induced gravity variations act simultaneously and in phase with the tidal force. The overall vertical tide-induced magma displacement in the conduit Δz_max can thus be larger then the individual mechanisms; i.e. $\mathit{\{}\mathrm{\Delta }{z}_{\mathrm{hr}},\mathrm{\Delta }{z}_{\mathrm{dec}}\mathit{\}}\le \mathrm{\Delta }{z}_{\max }<\mathrm{\Delta }{z}_{\mathrm{hr}}+\mathrm{\Delta }{z}_{\mathrm{dec}}$. In the following we focus on the reservoir expansion mechanism only in order to keep the derivation of the model parameters strictly analytical. The host rock mechanism is therefore reduced to establishing the required areal conduit cross section increase of ΔR_dec.

Baker (1984)Baker (1984)Zhang et al. (2007)Rust and Cashman (2011)Palma et al. (2011)Palma et al. (2011)Palma et al. (2011)Martel et al. (2018)Baker (1984)Baker (1984)

Table 1 Choice of model parameters, motivated by conditions at (1) Villarrica volcano located at 39.5^∘ S hosting a persistent lava lake of basaltic composition and (2) Cotopaxi volcano located at 0.7^∘ S, which preferentially erupts andesitic magma and intermittently is capped by a solid plug. If not stated otherwise, all numerical values in this paper are calculated with these parameters.

2.4 Radial flow profile in the conduit

The tide-induced vertical magma displacement in the conduit is delayed and extenuated by a viscosity-induced drag force. We access the temporal evolution and amplitude of the tide-induced displacement via the force (per unit mass) balance acting on the centre of mass of the magma column in the conduit:

$$\begin{array}{}\text{(3)}& {\displaystyle }\stackrel{\text{inner force}}{\overbrace{\underset{\text{drag force}}{\underbrace{\mathit{\gamma }\cdot \dot{z}\left(t\right)}}}}=\stackrel{\text{external force}}{\overbrace{\underset{\text{tidal force}}{\underbrace{a_{\mathrm{0}}\cdot \sin ({\mathit{\omega }}_{\mathrm{sd}}\cdot t)}}-\underset{\text{restoring force}}{\underbrace{{\mathit{\omega }}_{\mathrm{0}}^{\mathrm{2}}\cdot z\left(t\right)}}-\underset{\text{inertial force}}{\underbrace{\ddot{z}\left(t\right)}}}},\end{array}$$

where the two model parameters are the bulk damping rate γ and the eigenfrequency ω₀ of the magma column. The restoring force ensures that the centre of mass displacement tends towards the current “equilibrium” displacement associated with the current strength of the tidal force, i.e. $a_{\mathrm{0}}={\mathit{\omega }}_{\mathrm{0}}^{\mathrm{2}}\cdot \mathrm{\Delta }{z}_{\max }$. We further assume a Newtonian bulk drag force proportional to the flow velocity.

The continuity condition implies that the magma flows faster in the conduit centre than close to the boundary between the low-viscosity and high-viscosity magma or host rock. Accordingly, we assume a no-slip condition at the conduit boundary r=R_c and derive the analytical solution of the tide-induced parabolic vertical displacement profile z(r,t) in the conduit:

$$\begin{array}{}\text{(4)}& {\displaystyle }\left\{\begin{array}{ll}z(r,t)& =\mathrm{\Psi }\cdot \left[\mathrm{1}-{\left({\displaystyle \frac{r}{R_{\mathrm{c}}}}\right)}^{\mathrm{2}}\right]\cdot \sin ({\mathit{\omega }}_{\mathrm{sd}}\cdot t-{\mathit{\phi }}_{\mathrm{0}})\\ \mathrm{\Psi }& ={\displaystyle \frac{\mathrm{2}\cdot a_{\mathrm{0}}}{\sqrt{({\mathit{\omega }}_{\mathrm{0}}^{\mathrm{2}}-{\mathit{\omega }}_{\mathrm{sd}}^{\mathrm{2}}{)}^{\mathrm{2}}+(\mathit{\gamma }\cdot {\mathit{\omega }}_{\mathrm{sd}}{)}^{\mathrm{2}}}}}\\ {\mathit{\phi }}_{\mathrm{0}}& =\arctan \left({\displaystyle \frac{\mathit{\gamma }\cdot {\mathit{\omega }}_{\mathrm{sd}}}{{\mathit{\omega }}_{\mathrm{0}}^{\mathrm{2}}-{\mathit{\omega }}_{\mathrm{sd}}^{\mathrm{2}}}}\right)\\ \mathit{\gamma }& ={\displaystyle \frac{\mathrm{8}\cdot \mathit{\nu }}{R_{\mathrm{c}}^{\mathrm{2}}}}\\ {\mathit{\omega }}_{\mathrm{0}}^{\mathrm{2}}& ={\displaystyle \frac{a_{\mathrm{0}}}{\mathrm{\Delta }{z}_{\mathrm{dec}}}}={\displaystyle \frac{g\cdot \mathit{\pi }\cdot R_{\mathrm{c}}^{\mathrm{2}}}{{\mathit{\beta }}_{\mathrm{r}}\cdot V_{\mathrm{r}}\cdot D_{\mathrm{r}}\cdot (\mathrm{\nabla }p{)}_{\mathrm{vert}}}},\end{array}\right.\end{array}$$

with the radial coordinate $\mathrm{0}\le r\le R_{\mathrm{c}}$, the maximum vertical magma displacement amplitude Ψ (which equals twice the centre of mass displacement) and the phase shift φ₀ between tidal force and magma displacement in the conduit (see Appendix C).

For Villarrica, the model implies a tidal displacement amplitude of ${\mathrm{\Psi }}_{\mathrm{vill}}^{\mathrm{st}}=\mathrm{0.45}\mathrm{m}$, which lags behind the tide by ${\mathit{\phi }}_{\mathrm{0},\mathrm{vill}}\cdot {\mathit{\omega }}_{\mathrm{sd}}^{-\mathrm{1}}=\mathrm{2.0}\mathrm{h}$, where the displacement is predominantly limited by drag force. For Cotopaxi, the tidal displacement amplitude is ${\mathrm{\Psi }}_{\mathrm{coto}}^{\mathrm{st}}=\mathrm{0.09}\mathrm{m}$ and lags by ${\mathit{\phi }}_{\mathrm{0},\mathrm{coto}}\cdot {\mathit{\omega }}_{\mathrm{sd}}^{-\mathrm{1}}=\mathrm{0.2}\mathrm{h}$, where the displacement is predominantly limited by the restoring force. In comparison, the direct tide-induced gravity variations leads to a variation in the hydrostatic pressure of 10–100 Pa. In the context of the hydrostatic pressure gradient, this pressure variation has a similar effect as a vertical magma displacement by about 1 mm, thus rendering the direct tidal impact negligible compared to the indirect mechanism derived here.

3 Tide-enhanced bubble coalescence

Integrated over a semi-diurnal cycle, the tides do not result in a net magma displacement. A link from tides to degassing thus requires tide-enhanced mechanisms which irreversibly change the state of the magmatic gas phase. Bubble growth constitutes a predominantly exergonic and thus irreversible mechanism because the bubble surface tension inhibits or at least damps bubble shrinkage and dissolution (Prousevitch et al., 1993). Within a tide-induced radial displacement profile, neighbouring gas bubbles can exhibit differential tide-induced vertical displacements potentially enhancing the bubble coalescence rate (see Fig. 1c and Appendix D). The variation in the bubble coalescence rate leads to bigger bubbles and thus the tide can indeed modify an irreversible mechanism.

In this section, we set up a simplified formalisation of the magmatic gas phase and the typically predominant mechanisms which govern the bubble coalescence rate and estimate the relative tide-induced enhancement of the bubble coalescence by a comparison with these classical mechanisms. We consider a magma layer in the conduit at a particular depth; accordingly, the parameters discussed in the following describe the local conditions within a small volume of magma and should not be confused with the integrated bulk values for the total magma column. The variation in the tide-induced enhancement at different magma depths is discussed in the subsequent section.

3.1 Gas bubbles in magmatic melt

The dominant part of the magmatic volatile content is typically water, followed by carbon dioxide, sulfur compounds, and minor contributions from a large number of trace gases such as halogen compounds (Oppenheimer et al., 2014). For simplicity, we assume that all macroscopic properties of the gas phase are dominated by the degassing of water, in particular that the gas volume fraction ϕ exclusively consists of water vapour. The volatile solubility of magmatic melts is primarily pressure dependent, with secondary dependencies on temperature, melt composition, and volatile speciation (Gonnermann and Manga, 2013). The pressure dependency of the water solubility $C_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ in magmatic melt is given in a first approximation by $C_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}\left(p\right)=\sqrt{K_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}\cdot p}$ with the corresponding solubility coefficients $K_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ (find an empirical formulation in Zhang et al., 2007). For the local gas volume fraction ϕ(p) at a depth associated with the pressure p, we obtain

$$\begin{array}{}\text{(5)}& {\displaystyle }\mathit{\varphi }\left(p\right)={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{melt}}\left(p\right)}{{\mathit{\rho }}_{\mathrm{gas}}\left(p\right)}}\cdot \left(C_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}^{\mathrm{0}}-\sqrt{K_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}\cdot p}\right)\end{array}$$

with the total water weight fraction $C_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}^{\mathrm{0}}$ of the magmatic melt and the mass densities of the gas phase ρ_gas and of the overall melt (liquid + gas) ρ_melt.

The gas phase consists of separated bubbles as long as the gas volume fraction is below the percolation threshold of ${\mathit{\varphi }}_{\mathrm{perc}}=\mathrm{0.3}-\mathrm{0.7}$ (the variation is due to the range of different magmatic conditions;  Rust and Cashman, 2011). Bubbles typically vary in size following a power law (Cashman and Marsh, 1988; Blower et al., 2003) or a mixed power-law exponential distribution (Le Gall and Pichavant, 2016) and in shape from spherical to ellipsoidal (Rust et al., 2003; Moitra et al., 2013). While models based on polydisperse bubble size distributions are available (Sahagian and Proussevitch, 1998; Huber et al., 2013; Mancini et al., 2016), a common starting point to analyse the temporal evolution of the bubbles is nevertheless the assumption of a monodisperse size distribution of spherical bubbles (Prousevitch et al., 1993; Lensky et al., 2004).

We note the bubble size distribution ${\mathit{\delta }}_{\mathrm{b}}^{\mathrm{size}}(f\in {\mathbb{R}}^{+})$ with respect to the bubble radius (rather than the volume); i.e. the bubble radius is given by $r_{\mathrm{b}}=f\cdot R_{\mathrm{b}}$ with the hypothetical bubble radius R_b(p) of a monodisperse bubble size distribution. An estimate of a power-law bubble size distribution would require three parameters: the exponent and the lower and upper truncation cut-off (Lovejoy et al., 2004). An estimate of a mixed power-law exponential bubble size distribution would require at least two further parameters. The following analysis is conducted for an arbitrary bubble size distribution; nevertheless, for a basic quantitative estimate, we mimic a proper polydisperse bubble size distribution by the simpler single-parametric

$$\begin{array}{}\text{(6)}& {\displaystyle }{\stackrel{\mathrm{̃}}{\mathit{\delta }}}_{\mathrm{b}}^{\mathrm{size}}(f;q)=\left\{\begin{array}{lll}\mathrm{1}-q& :& f=\mathrm{1}\\ q& :& f=\sqrt[\mathrm{3}]{\mathrm{2}}\end{array}\right.,\end{array}$$

with $\mathrm{0}\le q<\frac{\mathrm{1}}{\mathrm{2}}$, which represents a monodisperse distribution except for a fraction of q bubbles which emerged from a past coalescence of two bubbles with f=1.

3.2 Bubble motion and bubble coalescence

Diffusion-driven volatile degassing can only take place in the immediate vicinity of a bubble and when the supersaturation pressure is larger than the bubble surface tension (Proussevitch and Sahagian, 2005). The volatile degassing rate is thus controlled by the spatial bubble distribution as well as the bubble size distribution (Lensky et al., 2004). Both distributions change during bubble rise, which is caused by a vertical ascent of the overall magma column or parcel with velocity v_melt and a superimposed bubble buoyancy with a velocity v_buoy which reads for a bubble with radius r_b (Stoke's law):

$$\begin{array}{}\text{(7)}& {\displaystyle }{v}_{\mathrm{buoy}}\left(r_{\mathrm{b}}\right)={\displaystyle \frac{\mathrm{2}\cdot g\cdot r_{\mathrm{b}}^{\mathrm{2}}}{\mathrm{9}\cdot \mathit{\nu }}}\cdot \left(\mathrm{1}-{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{gas}}}{{\mathit{\rho }}_{\mathrm{melt}}}}\right)\approx {\displaystyle \frac{\mathrm{2}\cdot g\cdot r_{\mathrm{b}}^{\mathrm{2}}}{\mathrm{9}\cdot \mathit{\nu }}}.\end{array}$$

If the buoyancy velocity is negligible compared to the magma ascent, the bubble flow is called “dispersed”; if the bubble buoyancy velocity contributes significantly to the overall bubble ascent, the bubble flow is called “separated” (Gonnermann and Manga, 2013). Rising bubbles grow continuously because of (1) decompression and (2) the increasing volatile degassing rate due to the associated decreases in the magmatic volatile solubility and of the bubble surface tension. Bubble coalescence accelerates the bubble growth.

Bubble coalescence requires two bubble walls to touch and ultimately to merge. Once two bubbles are sufficiently close to each other, near-field processes such as capillary and gravitational drainage cause a continuous reduction in the film thickness between the bubble walls until the bubbles merge after drainage times ranging from seconds to hours depending on the magmatic conditions (Herd and Pinkerton, 1997; Castro et al., 2012; Nguyen et al., 2013).

For small gas volume fractions, however, the initial distance between bubbles is large compared to the bubble dimensions and the coalescence rate is dominated by bubble transport mechanisms acting on longer length scales. Because bubble diffusion is typically negligibly small, bubble walls can only approach when a particular mechanism leads to differential bubble rise velocities or by bubble growth. In magmas with a sufficiently separated bubble flow, two neighbouring bubbles of different size can approach each other vertically due to the differential buoyancy velocities (Manga and Stone, 1994; Lovejoy et al., 2004). In magmas with a dispersed bubble flow, in contrast, the relative position of bubble centres remains fixed; thus, bubble coalescence is controlled by the bubble expansion rate caused by the ascent of the overall magma column (or affected magma parcel).

Figure 2 Relative contribution of the tidal mechanism (magnitude given by H_tide) on the bubble coalescence rate for a purely separated bubble flow (magnitude given by H_buoy) depending on the reference bubble radius R_b and the degree of polydispersity q. The reference bubble radius is reciprocally linked to the depth of the particular magma layer.

3.3 Comparison of bubble coalescence mechanism

The proposed tide-induced bubble transport mechanism is compared in the following with the classically predominant bubble transport or approaching mechanisms in order to estimate the relative contribution of the tidal mechanism on the overall coalescence rate. We access the (absolute) strength of a particular transport mechanism by its “collision volume” H_i (see Appendix D). The tidal mechanism is noted by H_tide. For comprehensibility, we focus on a comparison of the tidal mechanisms with the two “end-member” scenarios of a purely separated (H_buoy) and a purely dispersed (H_disp) bubble flow, respectively. A more comprehensive formulation of the classically predominant bubble transport or approaching mechanisms has been proposed, e.g. by Mancini et al. (2016).

For a separated bubble flow, the relative tidal contribution on the bubble coalescence rate depends reciprocally on the reference bubble radius R_b and on the degree of polydispersity q (Fig. 2). For $q=\mathrm{0.1}-\mathrm{0.4}$, the tidal mechanism contributes at least 10 % to the overall bubble coalescence rate for a range of reference bubble radii of R_b= 32–65 µm for Villarrica and R_b= 37–78 µm for Cotopaxi. For comparison, Le Gall and Pichavant (2016) obtained from basalt decompression experiments mean bubble radii of (at most, depending on the volatile content) 23 µm for a pressure of 100 MPa (∼ depth of 3.7 km) and of 80 µm for a pressure of 50 MPa (∼ depth of 1.9 km) and concluded an extensive bubble coalescence rate at depth associated with 50–100 MPa. Similarly, Castro et al. (2012) obtained from rhyolite decompression experiments mean bubble radii of 15 µm for a pressure of 100 MPa (∼ depth of 3.7 km) and of 30 µm for a pressure of 40 MPa (∼ depth of 1.5 km). For andesitic magma, the dependency of the bubble size on the pressure is presumably between the values for the basaltic and the rhyolitic magma. We conclude that the tidal mechanism can significantly contribute to the bubble coalescence rate in magma layers at a depth greater than 1 km, associated with bubble radii of 30–80 µm. In contrast, the tidal contribution becomes negligible at shallow levels once the bubble radii are in the millimetre-range which corresponds to the bubble size range at which bubbles efficiently start to segregate from the surrounding melt.

For a dispersed bubble flow, the relative tidal contribution on the bubble coalescence rate depends reciprocally on the magma ascent rate, hardly on the gas volume fraction ϕ, but it depends approximately linearly on the volatile content $C_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}^{\mathrm{0}}$ of the magma (Fig. 3). The tidal contribution causes an enhancement of the bubble coalescence rate equivalent to the enhancement caused by an increase in the magma ascent velocity by about 0.5 m h⁻¹ for Cotopaxi and 2.5 m h⁻¹ for Villarrica for the $C_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}^{\mathrm{0}}$ listed in Table 1. For comparison, the magma ascent velocities in passively degassing volcanic systems vary roughly between 1 and 100 m h⁻¹ (Cassidy et al., 2015, 2018). The tidal mechanism can accordingly contribute by at least several percent but potentially up to several multiples of 10% to the overall bubble coalescence rate. For gas volume fractions exceeding the minimum percolation threshold of ϕ_perc≈0.3, the model assumption of independent spherical bubbles increasingly loses its validity.

4 Discussion and conclusions

Our model implies a tide-induced periodical vertical magma displacement in the conduit within every semi-diurnal cycle in the order of 0.1–1 m due to magma expansion in the reservoir. At Villarrica, the modelled vertical magma displacement of 0.45 m implies a periodic variation in the lava lake level (whose areal cross section is about 10 times larger than for the conduit; Goto and Johnson, 2011) of about 5 cm. At Cotopaxi, the modelled vertical magma displacement of 0.09 m may apply additive stress on the solid plug.

Figure 3 Relative contribution of the tidal mechanism (magnitude given by H_tide) on the bubble coalescence rate for a purely dispersed bubble flow (magnitude given by H_disp) depending on the gas volume fraction and the initial water weight fraction of the magmatic melt. The corresponding values for ϕ are calculated with Eq. (5) assuming an ideal gas and magma temperatures of 1200 ^∘C for Villarrica and 1000 ^∘C for Cotopaxi. The relative tidal contribution is displayed as the equivalent to an enhancement of the magma ascent rate, which would have the same effect on the bubble coalescence rate. The model increasingly loses validity above the percolation threshold of ϕ_perc≈0.3.

We linked this magma displacement to bubble coalescence and compared the relative strength of the tide-induced bubble transport mechanism with respect to the classically predominant bubble transport mechanisms in magmas hosting a purely separated or a purely dispersed bubble flow. For both scenarios, we found that the tidal contributions to the overall bubble coalescence rate can be in the order of at least several percent up to several multiples of 10 % at a depth of several kilometres. At shallower depth, the direct tide-induced contribution to the overall bubble coalescence rate is rather negligible because the classical transport mechanisms become more efficient.

The tide-enhanced bubble coalescence rate at greater depth can nevertheless affect the gas phase in the overlying shallower layer because the additionally coalesced bubbles have a larger buoyancy velocity as well as a reduced surface tension and can thus stimulate on the one hand enhanced volatile degassing from the melt phase to the gas phase and on the other hand enhanced bubble coalescence rates in overlying layers (Prousevitch et al., 1993). These enhancements can ultimately cause the percolation of the gas phase at a somewhat greater depth compared to the tide-free scenario. In consequence, the magma becomes gas-permeable at this greater depth potentially causing enhanced volcanic gas emissions (Rust and Cashman, 2011; Gonnermann et al., 2017). The additional contributions from this greater depth to the volcanic gas emissions may also slightly shift the chemical composition of the overall gas emissions towards the chemical composition of the gas phase at this greater depth when compared to the tide-free scenario (Burton et al., 2007).

The quantitative results have been derived for the tidal forcing during spring tide. In contrast, the amplitude of the tide-induced mechanism is smaller by a factor of 3 during neap tide. Accordingly, the amplitude of the additional tide-induced contributions to the coalescence rate varies within a spring–neap tide cycle entailing a periodical signal with a period of about 14.8 d superimposed on the (nevertheless potentially much stronger) tide-independent coalescence rate. For a dispersed bubble flow scenario with rather fast magma ascent, a propagation of this superimposed signal from the enhanced coalescence rate via a variation in the percolation depth to the volcanic gas emissions is comprehensible. For a separated bubble flow scenario, however, the gas bubbles may need much more time than one spring–neap tide cycle to rise from a depth of several kilometres to the percolation depth. Magmatic systems can, however, become permeable already at a depth of 1–3 km (Edmonds and Gerlach, 2007; Burton et al., 2007), i.e. where the derived tidal effects are the strongest. In such a scenario, the tide-enhanced bubble coalescence rate could accordingly cause enhanced degassing without a significant delay.

In a scenario with a shallower percolation depth, the periodic pattern could nevertheless propagate to the degassing signal because several crucial parameters such as the mean bubble radius R_b and the gas volume fraction ϕ typically vary rather monotonously with pressure and thus depth (Gonnermann and Manga, 2013), implying a depth dependency of the relative tidal contributions to the bubble coalescence rate. Convolved along the vertical conduit axis, the tide-enhanced coalescence rate may accordingly preserve an overall periodicity driven by the dominant contributions from those magma layers which are particularly sensitive to the tidal mechanism. Moreover, this pressure dependency implies that gas contributions originating from the particularly tide-sensitive depths are more pronounced in the subsequent volcanic gas emissions during spring tide. Therefore, tide-induced variations in the chemical composition within the volcanic gas plumes may be particularly manifested in the relative molar degassing ratios (e.g. Burton et al., 2007; Bobrowski and Giuffrida, 2012; Balcone-Boissard et al., 2016) associated with these depths.

In conclusion, we traced a possible tidal impact from the tidal potential to a magma expansion in the reservoir, to a vertical magma displacement profile in the conduit, and to an enhanced bubble collision rate (and thus an enhanced coalescence rate), and this ultimately motivated a link between the tide-enhanced bubble coalescence rate and the periodical signal in the observed volcanic gas emissions. Furthermore, illustrative quantitative calculations indicated that the proposed tide-induced mechanism could lead to an enhancement of the bubble coalescence rate by up to several multiples of 10%. If propagated from enhanced bubble coalescence to a variation in the magnitude or chemical composition of the volcanic gas emissions, a periodical spring tide signal would be large enough to explain the observed about 2-weekly variations in volcanic gas emissions.

Nevertheless, our conceptual model only aimed at a proof of concept. Future studies may increase the complexity of the model by, e.g., (1) lifting several of our numerous simplifications (Appendix A), (2) incorporating macroscopic tidal mechanisms affecting the host rock explicitly, (3) adding several further microscopic mechanisms such as a tide-induced loosening of bubbles attached to the conduit walls or the tidal impact on crystal orientation, and (4) investigating possible non-linear interferences between the tide-induced dynamics and the tide-independent magma convection flow.

Appendix A: List of applied mayor simplifications

In our model we applied several simplifications regarding the shape and physical properties of the magma plumbing system. This we did for the sake of clarity and, even more importantly, in order to isolate the tide-induced effect on magma flow and degassing. To achieve this, we (1) modelled the tide-induced magma flow in the conduit neglecting any tide-independent magma dynamics such as magma convection, which implies an initial mechanical and thermodynamic equilibrium between magma and adjacent host rock. The only exception is the discussion of the impact of a constant magma ascent on the bubble coalescence rate. (2) The expansion of the initial conduit magma is neglected. We assume (3) a gas-tight host rock, (4) a cylindrical volcanic conduit, (5) a no-slip condition between conduit wall and magma, and (6) homogeneous magma flow properties. (7) The viscosity of the magma in the conduit is assessed by the effective bulk viscosity. (8) The radial tide-induced magma displacement is neglected. Moreover, (9) bubble coalescence is modelled by bubble collision, neglecting near-field drainage processes, bubble deformation processes, and post-collision coalescence processes. (10) Simple bubble size distributions are chosen, and (11) it is assumed that the volcanic gas phase exclusively consists of water vapour.

Appendix B: Quantitative estimates for the geometrical model parameters

The conduit radius is a crucial model parameter. The uppermost 200 m of Villarrica's conduit have frequently been exposed during the decades prior to the 2015 eruption due to pronounced oscillations of the lava lake level (Moussallam et al., 2016; Johnson et al., 2018b). The cross-sectional area of the conduit has a radius of about 30 m (Goto and Johnson, 2011), which at greater depths, however, narrows down to a mean radius of the order of R_c=6 m as is implied by studies based on gas emission magnitudes (Palma et al., 2011) and seismoacoustic properties (Richardson et al., 2014). The active vent of Cotopaxi was capped by an area of hot material with a diameter of 116–120 m during the eruption in 2015 (Johnson et al., 2018a). Although missing an empirical evidence, it is plausible that the mean conduit radius is somewhat narrower, and therefore we assume a (rather conservative) value of R_c=40 m.

Depth and volume of the magma reservoir constitute further crucial model parameters whose empirical estimates come with an even larger uncertainty. Seismic observations conducted at Villarrica imply the existence of a shallow magma reservoir with a lateral diameter of at least 5 km and a vertical extent of about 2.5 km whose centre of mass is located at a depth of around D_r=3 km below the summit (Mora-Stock, 2015), implying a conduit length of about L_c=2 km. Assuming an ellipsoidal magma reservoir, this implies a magma reservoir volume of V_r=35 km³ at Villarrica. The magmatic system of Cotopaxi in contrast seems to be more complex and hosts a rather small magma pocket (2 km³) beneath the SW flank at a depth of about 4 km below the summit (Hickey et al., 2015). Furthermore, seismic observations revealed fluid movements (magma and/or hydrothermal fluids) within a centrally located 85 km³ column spanning 2 to 14 km depth below the summit (Ruiz et al., 1998). This fluid column is assumed to connect the laterally offset shallow pocket with two much larger deeper magma reservoirs, which are situated between 7 and 11 km and somewhere at a depth greater than 16 km below the summit (Arias et al., 2015; Mothes et al., 2017; Martel et al., 2018). For heating 85 km³ of rock, these deep-seated magma reservoirs may be rather large. Missing any accurate volume estimate, we estimate that the upper of the two deep-seated reservoirs hosts a magma volume of V_r=35 km³ with a centre of mass depth of D_r=8 km. The choice of equal reservoir volumes for both Villarrica and Cotopaxi allows for a better comparison of the impact of varying the other volcanic parameters. Further, we assume the small magma pocket as the lower end of the conduit, i.e. with a conduit length of L_c=4 km.

Appendix C: Calculation of tide-induced conduit flow

Oscillating centre of mass displacement. After a negligible settling time, the driven oscillator described by Eq. (3) oscillates with semi-diurnal periodicity, and we obtain the general long-term solution

$$\begin{array}{}\text{(C1)}& {\displaystyle }\left\{\begin{array}{ll}z\left(t\right)& =z_{\mathrm{0}}\cdot \sin ({\mathit{\omega }}_{\mathrm{sd}}\cdot t-{\mathit{\phi }}_{\mathrm{0}})\\ z_{\mathrm{0}}& ={\displaystyle \frac{a_{\mathrm{0}}}{\sqrt{({\mathit{\omega }}_{\mathrm{0}}^{\mathrm{2}}-{\mathit{\omega }}_{\mathrm{sd}}^{\mathrm{2}}{)}^{\mathrm{2}}+(\mathit{\gamma }\cdot {\mathit{\omega }}_{\mathrm{sd}}{)}^{\mathrm{2}}}}}\\ {\mathit{\phi }}_{\mathrm{0}}& =\arctan \left({\displaystyle \frac{\mathit{\gamma }\cdot {\mathit{\omega }}_{\mathrm{sd}}}{{\mathit{\omega }}_{\mathrm{0}}^{\mathrm{2}}-{\mathit{\omega }}_{\mathrm{sd}}^{\mathrm{2}}}}\right)\end{array}.\right.\end{array}$$

Navier–Stokes equation for periodical pipe flow. When exposed to a constant force (per unit mass) $f_{\mathrm{ext}}^{\mathrm{0}}$, a viscous fluid in a cylindrical pipe with radius R_c flows with a parabolic velocity profile v⁰(r), $\mathrm{0}\le r\le R_{\mathrm{c}}$:

$$\begin{array}{}\text{(C2)}& {\displaystyle }{v}^{\mathrm{0}}\left(r\right)={\displaystyle \frac{R_{\mathrm{c}}^{\mathrm{2}}\cdot f_{\mathrm{ext}}^{\mathrm{0}}}{\mathrm{4}\cdot \mathit{\nu }}}\left[\mathrm{1}-{\left({\displaystyle \frac{r}{R_{\mathrm{c}}}}\right)}^{\mathrm{2}}\right].\end{array}$$

When exposed to a periodically varying and thus time-dependent external force $f_{\mathrm{ext}}\left(t\right)=f_{\mathrm{ext}}^{\mathrm{0}}\cdot e^{i\mathit{\omega }t}$, the analytical solution of the flow profile is more complicated (Spurk, 1997):

$$\begin{array}{}\text{(C3)}& {\displaystyle }v(r,t)=\stackrel{\mathrm{‾}}{v^{\mathrm{0}}\left(r\right)}\cdot \mathrm{\Re }\left[-i\cdot {\displaystyle \frac{\mathrm{8}}{N^{\mathrm{2}}}}\cdot e^{i\mathit{\omega }t}\cdot \left(\mathrm{1}-{\displaystyle \frac{J_{\mathrm{0}}\left(\sqrt{-i}N\frac{r}{R_{\mathrm{c}}}\right)}{J_{\mathrm{0}}\left(\sqrt{-i}N\right)}}\right)\right],\end{array}$$

with the centre of mass velocity $\stackrel{\mathrm{‾}}{v^{\mathrm{0}}\left(r\right)}$ of a constant forcing (see Eq. C2), the real part $\mathrm{\Re }[..]$, the imaginary unit i, the Bessel function $J_{\mathrm{0}}(..)$, and the dimensionless parameter $N=\sqrt{\frac{\mathit{\omega }}{\mathit{\nu }}}\cdot R_{\mathrm{c}}$. In the limit N→0, the velocity profile asymptotically adopts the time dependency as well as the magnitude of the external force. For N=1 the exact magnitude is already $\mathrm{0.98}\cdot f_{\mathrm{ext}}^{\mathrm{0}}$, and the radial profile shows hardly any deviation from a parabolic profile. For the chosen model parameters (Table 1) and ω=ω_sd, we obtain N≈0.2 and thus Eq. (C3) reduces in very good approximation to the familiar

$$\begin{array}{}\text{(C4)}& {\displaystyle }v(r,t)\approx {\displaystyle \frac{R_{\mathrm{c}}^{\mathrm{2}}\cdot f_{\mathrm{ext}}\left(t\right)}{\mathrm{4}\cdot \mathit{\nu }}}\left[\mathrm{1}-{\left({\displaystyle \frac{r}{R_{\mathrm{c}}}}\right)}^{\mathrm{2}}\right].\end{array}$$

Derivation of the equation of motion (Eq. 4). The vertical velocity of the centre of mass can be obtained as $\dot{z}\left(t\right)$ = z₀ ⋅ ω_sd ⋅ $\cos ({\mathit{\omega }}_{\mathrm{sd}}\cdot t-{\mathit{\phi }}_{\mathrm{0}})$ from Eq. (C1) and as $v\left(t\right)=(\mathit{\pi }\cdot R_{\mathrm{c}}^{\mathrm{2}}{)}^{-\mathrm{1}}\cdot {\int }_{\mathrm{0}}^{R_{\mathrm{c}}}v(r,t)\cdot \mathrm{2}\mathit{\pi }r\mathrm{d}r=\frac{R_{\mathrm{c}}^{\mathrm{2}}}{\mathrm{8}\cdot \mathit{\nu }}\cdot f_{\mathrm{ext}}(t)$ from Eq. (C4). Further, we know $f_{\mathrm{ext}}\left(t\right)=f_{\mathrm{int}}\left(t\right)=\mathit{\gamma }\cdot \dot{z}\left(t\right)$ from Eq. (3). Applying f_ext(t) to Eq. (C4) reveals $\mathit{\gamma }=\frac{\mathrm{8}\cdot \mathit{\nu }}{R_{\mathrm{c}}^{\mathrm{2}}}$ and ultimately the fully parameterised equation of motion in Eq. (4).

Appendix D: Calculation of the collision volumes

As is common for most coalescence models (including those cited above), we consider spherical bubbles only. Two spherical bubbles with radii f₁⋅R_b and f₂⋅R_b (f₁ and f₂ drawn from ${\mathit{\delta }}_{\mathrm{b}}^{\mathrm{size}}\left(f\right)$) collide as soon as the distance between their bubble centres is $r_{\mathrm{coal}}=(f_{\mathrm{1}}+f_{\mathrm{2}})\cdot R_{\mathrm{b}}$. We introduce the collision volume $H(f_{\mathrm{1}},f_{\mathrm{2}};\mathrm{\Delta }t)$ associated with a bubble with radius f₁⋅R_b as the volume enclosing all possible initial locations of the bubble centre of another bubble with radius f₂⋅R_b such that both bubbles collide (and thus coalesce) at the latest after a time interval Δt. All bubble collision mechanisms are derived as enhancements of the initial static collision volume

$$\begin{array}{}\text{(D1)}& {\displaystyle }{H}_{\mathrm{0}}(f_{\mathrm{1}},f_{\mathrm{2}})={\displaystyle \frac{\mathrm{4}\mathit{\pi }}{\mathrm{3}}}\cdot R_{\mathrm{b}}^{\mathrm{3}}\cdot (f_{\mathrm{1}}+f_{\mathrm{2}}{)}^{\mathrm{3}},\end{array}$$

and we consider only those bubble pairs which have not collided already in the initial state. The absolute enhancement of the collision volume due to a particular bubble collision mechanism divided by Δt thus gives the enhancement of the bubble collision rate contributed by the particular mechanism. Because the tide-induced mechanisms are derived for a semi-diurnal cycle, the relative strengths of all coalescence mechanisms are compared with respect to this time interval Δt_sd.

The collision volumes of the different collision mechanisms are all derived with the same approach: we fix the position of a bubble with arbitrary radius f₁⋅R_b and derive $H(f_{\mathrm{1}},f_{\mathrm{2}};\mathrm{\Delta }t)$ with respect to the relative motion of another bubble with arbitrary radius f₂⋅R_b. In each case the initial collision volume H₀(f₁,f₂) is subtracted either already tacitly in the motivation or explicitly mathematically. Higher-order details such as the influence of a third bubble on the numeric results are ignored.

Tide-enhanced bubble collision volume. We fix the horizontal coordinates $(r,\mathit{\phi }{)}_{\mathrm{bubble}\mathrm{1}}=(r_{\mathrm{0}},\mathrm{0})$, $\mathrm{0}\le r_{\mathrm{0}}\le R_{\mathrm{c}}$, of the first bubble, where the cylindrical symmetry of the conduit allows us to pick the azimuth angle without loss of generality and vary the horizontal coordinates $(r,\mathit{\phi }{)}_{\mathrm{bubble}\mathrm{2}}=(r,\mathit{\phi })$ of a second bubbles. The horizontal distance h between the two bubbles is thus given by $r^{\mathrm{2}}=r_{\mathrm{0}}^{\mathrm{2}}-\mathrm{2}\cdot r_{\mathrm{0}}\cdot h\cdot \cos \left(\mathit{\phi }\right)+h^{\mathrm{2}}$. Within a semi-diurnal cycle, the peak-to-peak differential tide-induced vertical displacement of two bubbles at the radial coordinates r and r₀ is given by Δ z_tide(r,r₀) = $\mathrm{2}\cdot \left|z_{\mathrm{0}}\right(r)-z_{\mathrm{0}}(r_{\mathrm{0}}\left)\right|$ (see Eq. 4). The tide-induced collision volume is then the integral of Δz_tide(r,r₀) integrated over a circle with radius r_coal:

$$\begin{array}{}\text{(D2)}& {\displaystyle }{H}_{\mathrm{tide}}\left(r_{\mathrm{0}}\right)& {\displaystyle }=\underset{\mathrm{0}}{\overset{r_{\mathrm{coal}}}{\int }}\mathrm{d}hh\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\mathrm{d}\mathit{\phi }\mathrm{\Delta }{z}_{\mathrm{tide}}(r,r_{\mathrm{0}})\text{(D3)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{4}\cdot \mathrm{\Psi }\cdot r_{\mathrm{0}}}{R_{\mathrm{c}}^{\mathrm{2}}}}\underset{\mathrm{0}}{\overset{r_{\mathrm{coal}}}{\int }}\mathrm{d}h{h}^{\mathrm{2}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\mathrm{d}\mathit{\phi }\left|\cos \right(\mathit{\phi })-{\displaystyle \frac{h}{\mathrm{2}{r}_{\mathrm{0}}}}|.\end{array}$$

This integral has to be split into two integrals at the angles where the sign of the absolute function changes, which is the case at $\pm {\mathit{\phi }}^{\prime }=\pm \arccos \left(\frac{h}{\mathrm{2}{r}_{\mathrm{0}}}\right)\approx \pm \frac{\mathit{\pi }}{\mathrm{2}}$:

$$\begin{array}{}\text{(D4)}& {\displaystyle }{H}_{\mathrm{tide}}\left(r_{\mathrm{0}}\right)& {\displaystyle }={\displaystyle \frac{\mathrm{16}\cdot \mathrm{\Psi }\cdot r_{\mathrm{0}}}{R_{\mathrm{c}}^{\mathrm{2}}}}\underset{\mathrm{0}}{\overset{r_{\mathrm{coal}}}{\int }}\mathrm{d}h{h}^{\mathrm{2}}\underset{\approx \mathrm{1}\text{for}h\ll r_{\mathrm{0}}}{\underbrace{\left[\sin \left({\mathit{\phi }}^{\prime }\right)-\cos \left({\mathit{\phi }}^{\prime }\right)\cdot {\mathit{\phi }}^{\prime }\right]}}\text{(D5)}& {\displaystyle }& {\displaystyle }\approx {\displaystyle \frac{\mathrm{16}\cdot \mathrm{\Psi }\cdot r_{\mathrm{0}}}{R_{\mathrm{c}}^{\mathrm{2}}}}\cdot {\displaystyle \frac{r_{\mathrm{coal}}^{\mathrm{3}}}{\mathrm{3}}}\text{(D6)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{4}\cdot \mathrm{\Psi }\cdot r_{\mathrm{0}}}{\mathit{\pi }\cdot R_{\mathrm{c}}^{\mathrm{2}}}}\cdot H_{\mathrm{0}}(f_{\mathrm{1}},f_{\mathrm{2}}).\end{array}$$

We integrate H_tide(r₀) over the local spatial bubble distribution in the conduit in order to obtain the average effect. We parameterise the (isotropic) spatial bubble distribution by the depth-independent ${\mathit{\delta }}_{\mathrm{b}}^{\mathrm{spatial}}\left(r_{\mathrm{0}}\right)=(\mathrm{1}+\mathit{\alpha })\cdot \frac{\mathrm{1}}{R}\cdot (\frac{r_{\mathrm{0}}}{R}{)}^{\mathit{\alpha }}$, which is a homogeneous distribution for α=1 but with all bubbles at the conduit wall if α→∞, respectively. For the averaged tide-induced collision volume, we obtain

$$\begin{array}{}\text{(D7)}& {\displaystyle }{H}_{\mathrm{tide}}& {\displaystyle }=\underset{\mathrm{0}}{\overset{R}{\int }}{\mathit{\sigma }}_{\mathrm{tide}}\left(r_{\mathrm{0}}\right)\cdot {\mathit{\delta }}_{\mathrm{b}}^{\mathrm{spatial}}\left(r_{\mathrm{0}}\right)\cdot \mathrm{d}{r}_{\mathrm{0}}\text{(D8)}& {\displaystyle }& {\displaystyle }=\underset{\text{distribution}}{\underbrace{\left[{\displaystyle \frac{\mathrm{1}+\mathit{\alpha }}{\mathrm{2}+\mathit{\alpha }}}\right]}}\cdot \underset{\text{tidal}}{\underbrace{\left[{\displaystyle \frac{\mathrm{4}\cdot \mathrm{\Psi }}{\mathit{\pi }\cdot R_{\mathrm{c}}}}\right]}}\cdot \underset{\text{scale}}{\underbrace{H_{\mathrm{0}}(f_{\mathrm{1}},f_{\mathrm{2}})}}.\end{array}$$

The “distribution term” is $\frac{\mathrm{2}}{\mathrm{3}}$ for an isotropic bubble distribution and approaches unity if all bubbles are close to the host rock. Arguably, the conditions for crystal nucleation and thus bubble nucleation are better close to the host rock where the magma is cooler and more crystals and thus nucleation possibilities are available. Following this reasoning but also because we want to examine the maximum possible tidal impact, we set the distribution term to unity. The “tidal term” contains the information on the scale of the effective tide-induced impact. The “scale term” contains the information on the actual bubble size distribution, highlighting that the relative tidal enhancement is identical for any bubble size distribution, at least in our simple model.

Buoyancy-induced bubble collision volume. Two bubbles with radii $f_{\mathrm{1}}\cdot R_{\mathrm{b}}\ne f_{\mathrm{2}}\cdot R_{\mathrm{b}}$ have a differential rise velocity $\mathrm{\Delta }{v}_{\mathrm{buoy}}=|f_{\mathrm{2}}^{\mathrm{2}}-f_{\mathrm{1}}^{\mathrm{2}}|\cdot v_{\mathrm{buoy}}\left(R_{\mathrm{b}}\right)$, and thus their relative distance changes during the rise. The two bubbles will collide if the larger bubble is below the smaller and if the horizontal distance between their bubble centres is at most r_coal. Accordingly, the buoyancy-induced collision volume H_buoy is a cylindrical volume with base area $\mathit{\pi }\cdot r_{\mathrm{coal}}^{\mathrm{2}}$ and cylinder length Δv_buoy⋅Δt_sd:

$$\begin{array}{}\text{(D9)}& {\displaystyle }{H}_{\mathrm{buoy}}& {\displaystyle }(f_{\mathrm{1}},f_{\mathrm{2}})=\mathit{\pi }\cdot r_{\mathrm{coal}}^{\mathrm{2}}\cdot |f_{\mathrm{2}}^{\mathrm{2}}-f_{\mathrm{1}}^{\mathrm{2}}|\cdot v_{\mathrm{buoy}}\left(R_{\mathrm{b}}\right)\cdot \mathrm{\Delta }{t}_{\mathrm{sd}}\text{(D10)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{3}\cdot |f_{\mathrm{2}}-f_{\mathrm{1}}|}{\mathrm{4}\cdot R_{\mathrm{b}}}}\cdot v_{\mathrm{buoy}}\left(R_{\mathrm{b}}\right)\cdot \mathrm{\Delta }{t}_{\mathrm{sd}}\cdot H_{\mathrm{0}}(f_{\mathrm{1}},f_{\mathrm{2}}).\end{array}$$

For a given pair of bubbles with radii $f_{\mathrm{1}}\cdot R_{\mathrm{b}}\ne f_{\mathrm{2}}\cdot R_{\mathrm{b}}$, f₁ and f₂ drawn from ${\mathit{\delta }}_{\mathrm{b}}^{\mathrm{size}}\left(f\right)$, the ratio of the contribution from the tide-induced and the buoyancy-induced collision mechanisms is

$$\begin{array}{}\text{(D11)}& {\displaystyle \frac{H_{\mathrm{tide}}}{H_{\mathrm{buoy}}}}={\displaystyle \frac{\mathrm{24}\cdot \mathrm{\Psi }\cdot \mathit{\nu }}{\mathit{\pi }\cdot R_{\mathrm{c}}\cdot |f_{\mathrm{1}}-f_{\mathrm{2}}|\cdot g\cdot R_{\mathrm{b}}\cdot \mathrm{\Delta }{t}_{\mathrm{sd}}}}.\end{array}$$

The bulk ratio (with respect to the local magma layer) can be obtained by a previous and separate integration of H_tide and H_buoy over f₁ and f₂ with respect to the actual bubble size distribution ${\mathit{\delta }}_{\mathrm{b}}^{\mathrm{size}}\left(f\right)$ (rather than integrating Eq. D11). For the explicit bubble size distribution ${\stackrel{\mathrm{̃}}{\mathit{\delta }}}_b^{\mathrm{size}}$ from Eq. (6), we obtain the bulk collision volumes ${\stackrel{\mathrm{̃}}{H}}_{\mathrm{tide}}$ and ${\stackrel{\mathrm{̃}}{H}}_{\mathrm{buoy}}$,

$$\begin{array}{}\text{(D12)}& {\displaystyle }& {\displaystyle \frac{{\stackrel{\mathrm{̃}}{H}}_{\mathrm{tide}}\left(q\right)}{H_{\mathrm{0}}(\mathrm{1},\mathrm{1})}}=(\mathrm{1}+\mathrm{0.89}\cdot q+\mathrm{0.11}\cdot q^{\mathrm{2}})\cdot {\displaystyle \frac{\mathrm{4}\cdot \mathrm{\Psi }}{\mathit{\pi }\cdot R_{\mathrm{c}}}},\text{(D13)}& {\displaystyle }& {\displaystyle \frac{{\stackrel{\mathrm{̃}}{H}}_{\mathrm{buoy}}\left(q\right)}{H_{\mathrm{0}}(\mathrm{1},\mathrm{1})}}=(q-q^{\mathrm{2}})\cdot {\displaystyle \frac{\mathrm{9}}{\mathrm{16}\cdot R_{\mathrm{b}}}}\cdot v_{\mathrm{buoy}}\left(R_{\mathrm{b}}\right)\cdot \mathrm{\Delta }{t}_{\mathrm{sd}},\end{array}$$

and thus the bulk ratio (used for the calculation of Fig. 2):

$$\begin{array}{}\text{(D14)}& {\displaystyle \frac{{\stackrel{\mathrm{̃}}{H}}_{\mathrm{tide}}}{{\stackrel{\mathrm{̃}}{H}}_{\mathrm{buoy}}}}{\displaystyle }=\mathrm{60}\cdot (\mathrm{0.9}+{\displaystyle \frac{\mathrm{1}+q^{\mathrm{2}}}{q-q^{\mathrm{2}}}})\cdot {\displaystyle \frac{\mathit{\nu }\left[{\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}\right]\cdot \mathrm{\Psi }\left[\mathrm{m}\right]}{R_{\mathrm{c}}\left[\mathrm{m}\right]\cdot R_{\mathrm{b}}\left[\mathrm{\mu }\mathrm{m}\right]}}.\end{array}$$

Growth-induced bubble collision volume. In magma with a dispersed bubble flow (v_buoy≪v_melt), a rising bubble exhibits a pressure decrease rate of

$$\begin{array}{}\text{(D15)}& {\displaystyle \frac{\mathrm{\Delta }p}{\mathrm{\Delta }t}}=v_{\mathrm{melt}}\cdot (\mathrm{\nabla }p{)}_{\mathrm{vert}}.\end{array}$$

Ignoring accompanying changes in secondary parameters such as melt temperature and magma composition and assuming for simplicity a monodisperse bubble size distribution (thus $R_{\mathrm{b}}^{\mathrm{3}}\propto \mathit{\varphi }$), for the enhancement of the collision volume due to a rise-driven pressure decrease by Δp≪p₀ (apply Eq. 5 on Eq. D1), we obtain

$$\begin{array}{ll}{\displaystyle }{H}_{\mathrm{disp}}(\mathrm{\Delta }p;p_{\mathrm{0}})& {\displaystyle }=H_{\mathrm{0}}\left(R_{\mathrm{b}}\right(p_{\mathrm{0}}-\mathrm{\Delta }p\left)\right)-H_{\mathrm{0}}\left(R_{\mathrm{b}}\right(p_{\mathrm{0}}\left)\right)\\ {\displaystyle }& {\displaystyle }=H_{\mathrm{0}}(\mathrm{1},\mathrm{1})\cdot {\displaystyle \frac{C_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}^{\mathrm{0}}-\frac{\mathrm{1}}{\mathrm{2}}\sqrt{K_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}\cdot p_{\mathrm{0}}}}{C_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}^{\mathrm{0}}-\sqrt{K_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}\cdot p_{\mathrm{0}}}}}\cdot {\displaystyle \frac{\mathrm{\Delta }p}{p_{\mathrm{0}}}}\\ \text{(D16)}& {\displaystyle }& {\displaystyle }+\mathcal{O}{\left[({\displaystyle \frac{\mathrm{\Delta }p}{p_{\mathrm{0}}}}\right)}^{\mathrm{2}}],\end{array}$$

where we assume that ρ_melt is constant and ρ_gas follows the ideal gas law. Inserting Eq. (D15) in Eq. (D16), we obtain

$$\begin{array}{ll}{\displaystyle \frac{H_{\mathrm{disp}}\left(p_{\mathrm{0}}\right)}{H_{\mathrm{0}}(\mathrm{1},\mathrm{1})}}& {\displaystyle }={\displaystyle \frac{C_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}^{\mathrm{0}}-\frac{\mathrm{1}}{\mathrm{2}}\sqrt{K_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}\cdot p_{\mathrm{0}}}}{C_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}^{\mathrm{0}}-\sqrt{K_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}\cdot p_{\mathrm{0}}}}}\\ \text{(D17)}& {\displaystyle }& {\displaystyle }\cdot v_{\mathrm{melt}}\cdot \mathrm{\Delta }{t}_{\mathrm{sd}}\cdot {\displaystyle \frac{(\mathrm{\nabla }p{)}_{\mathrm{vert}}}{p_{\mathrm{0}}}}.\end{array}$$

The ratio of the contribution from the tide-induced and the growth-induced collision mechanism (used for the calculation of Fig. 3) is

$$\begin{array}{}\text{(D18)}& {\displaystyle \frac{H_{\mathrm{tide}}}{H_{\mathrm{disp}}}}{\displaystyle }=\underset{\approx \mathrm{0.25}-\mathrm{0.5}}{\underbrace{{\displaystyle \frac{C_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}^{\mathrm{0}}-\sqrt{K_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}\cdot p_{\mathrm{0}}}}{C_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}^{\mathrm{0}}-\frac{\mathrm{1}}{\mathrm{2}}\sqrt{K_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}\cdot p_{\mathrm{0}}}}}}}\cdot {\displaystyle \frac{\mathrm{4}\cdot \mathrm{\Psi }\left[\mathrm{m}\right]\cdot p_{\mathrm{0}}\left[\mathrm{MPa}\right]}{R_{\mathrm{c}}\left[\mathrm{m}\right]\cdot v_{\mathrm{melt}}\left[{\mathrm{m}\mathrm{h}}^{-\mathrm{1}}\right]}}.\end{array}$$