# Persistent CO₂ emissions and hydrothermal unrest following the 2015 earthquake in Nepal

Frédéric Girault, Lok Bijaya Adhikari, Christian France-Lanord, Pierre Agrinier, Bharat Koirala, Mukunda Bhattarai, Sudhan Mahat, Chiara Groppo, Franco Rolfo, Laurent Bollinger, et al.

<figure><img src="image_1.png" /></figure>

► To cite this version:

Frédéric Girault, Lok Bijaya Adhikari, Christian France-Lanord, Pierre Agrinier, Bharat Koirala, et al. Persistent CO₂ emissions and hydrothermal unrest following the 2015 earthquake in Nepal. Nature Communications, 2018, 9, pp.2956. 10.1038/s41467-018-05138-z. cea-01880809

HAL Id: cea-01880809

https://cea.hal.science/cea-01880809

Submitted on 7 Jun 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

ARTICLE

DOI: 10.1038/s41467-018-05138-z

OPEN

# Persistent CO₂ emissions and hydrothermal unrest following the 2015 earthquake in Nepal

Frédéric Girault ID <fnref n="1" />, Lok Bijaya Adhikari², Christian France-Lanord ID ³, Pierre Agrinier⁴, Bharat P. Koirala², Mukunda Bhattarai², Sudhan S. Mahat⁵, Chiara Groppo ID ⁶, Franco Rolfo ID ⁶, Laurent Bollinger⁷ & Frédéric Perrier<fnref n="1" />

Fluid-earthquake interplay, as evidenced by aftershock distributions or earthquake-induced effects on near-surface aquifers, has suggested that earthquakes dynamically affect permeability of the Earth's crust. The connection between the mid-crust and the surface was further supported by instances of carbon dioxide (CO₂) emissions associated with seismic activity, so far only observed in magmatic context. Here we report spectacular non-volcanic CO₂ emissions and hydrothermal disturbances at the front of the Nepal Himalayas following the deadly 25 April 2015 Gorkha earthquake (moment magnitude M<sub>w</sub> = 7.8). The data show unambiguously the appearance, after the earthquake, sometimes with a delay of several months, of CO₂ emissions at several sites separated by > 10 kilometres, associated with persistent changes in hydrothermal discharges, including a complete cessation. These observations reveal that Himalayan hydrothermal systems are sensitive to co- and post-seismic deformation, leading to non-stationary release of metamorphic CO₂ from active orogens. Possible pre-seismic effects need further confirmation.

<footnote n="1">¹Physics of Natural Sites, Institut de Physique du Globe de Paris, Sorbonne Paris Cité, Université Paris Diderot, CNRS UMR 7154, 1 rue Jussieu, F-75005 Paris, France. ²Department of Mines and Geology, National Seismological Centre, Lainchaur, Kathmandu, Nepal. ³Centre de Recherches Pétrographiques et Géochimiques, Université de Nancy, CNRS UMR 7358, F-54500 Vandoeuvre-lès-Nancy, France. ⁴Stable Isotopes Geochemistry, Institut de Physique du Globe de Paris, Sorbonne Paris Cité, Université Paris Diderot, CNRS UMR 7154, 1 rue Jussieu, F-75005 Paris, France. ⁵Sanjen Jalavidhyut Company Limited, Lazimpat, Kathmandu, Nepal. ⁶Department of Earth Sciences, IGG-CNR, University of Turin, Via Verdi, 8, I-10124 Turin, Italy. ⁷CEA, DAM, DIF, F-91297 Arpajon, France. Correspondence and requests for materials should be addressed to F.G.(email: girault@ipgp.fr)</footnote>
Understanding fluid-earthquake interplay has long received a lot of attention¹. For the last 20 years, there has been a growing number of evidence for fluid-driven earthquakes². Deep fluids have been shown controlling aftershock distribution in the tectonic contexts of rifting³⁻⁵, subduction⁶,⁷, and reverse⁸ and strike-slip faulting⁹. Contemporaneously, numerous observations of earthquake-induced effects on near-surface aquifers have been reported¹⁰, including mainly changes in stream and spring discharge¹¹⁻¹³, groundwater level¹⁴,¹⁵ and temperature¹⁶⁻¹⁹ in various tectonic contexts. These observations have suggested that earthquakes dynamically affect permeability of the Earth's crust²⁰,²¹. Carbon dioxide (CO₂) emissions were observed in association with seismicity in the case of the Matsushiro swarm⁹ in Japan or, recently, at the Lassen volcano²² (Cascades Range, USA) and in the Eger Rift²³ (Czech Republic), which suggests connection between the mid-crust and the ground surface through gas transport. However, such examples so far were only observed in the presence of magmatic activity.

The Himalayan orogen results from the India-Eurasia collision²⁴, the Main Himalayan Thrust (MHT) accommodating at 2 cm year⁻¹ half of the shortening between the two continents²⁵. The largest earthquake in Nepal before 2015, the 1934 Bihar-Nepal earthquake (moment magnitude M<sub>w</sub>~8.2) in Eastern Nepal, which claimed >15,000 lives, ruptured the MHT up to the surface over 150 kilometres along the Main Frontal Thrust²⁶. Inter-seismic deformation is associated with intense background seismicity²⁷, with 4-5 events of local magnitude M<sub>L</sub> > 5 per year,

concentrated between 10 and 25 kilometres depth at the foot of the Himalayan topographic rise²⁸. This region – the Main Central Thrust (MCT) zone²⁹ – also exhibits numerous hydrothermal systems³⁰ (Fig. 1). Following evidence of degassing from chemical and isotopic analysis of hot springs and rivers³¹,³², large CO₂ emissions were discovered near hot springs³⁰,³³, with CO₂ fluxes at places similar to diffusive fluxes from active volcanoes³⁰,³⁴. The seasonal and yearly stability of the soil-gas radon concentration time-series³⁴,³⁵, the invariant radon-CO₂ fluxes relationship³⁰,³⁶, and the results of watering experiments³⁶ at selected sites attest to the remarkable temporal stability of these hydrothermal systems, even during monsoon. These non-volcanic CO₂ emissions are characterised by radiogenic helium, high radon content, and carbon isotopic compositions suggesting metamorphic CO₂ production at >5 kilometres depth³⁰⁻³²,³⁴. In the accepted conceptual model³¹,³², from the decarbonation source at >5 kilometres depth, CO₂ percolates through fracture networks in the MCT zone, where it mixes with meteoric water. Near the water table, a fraction of CO₂ may degas, and water discharges eventually at the surface as a hot spring. Degassed CO₂ may also be transported faster to the surface through a network of interrelated faults without interaction with hydrothermal circulations³⁷.

The M<sub>w</sub> = 7.8 25 April 2015 Gorkha earthquake³⁸ (Fig. 1) caused >9000 deaths, with fatality rates >1% in mountainous areas north of Kathmandu³⁹, reaching 100% at some places⁴⁰. It partly ruptured the MHT along a 120-kilometre-long segment east from the epicentre⁴¹, whose northern limit coincides with the

<figure><img src="image_4.png" /></figure>

**Fig. 1** Post-seismic responses of carbon dioxide emissions and hot springs following the 25 April 2015 Gorkha earthquake in Central Nepal. Post-seismic responses are shown as coloured symbols: positive effect (red), negative effect (green), absence of any significant change (blue), or when it is unknown (grey). Aftershock epicentres (M<sub>L</sub> ≥ 3) from 25 April 2015 to 31 December 2017 (Nepal National Seismological Centre, NSC) are plotted as yellow dots⁴³. The Main Central Thrust (MCT) shear zone (in dark grey), which branches at depth on the MHT, separates low-grade Lesser Himalayan Sequence (LHS) and high-grade Greater Himalayan Sequence (GHS) to the north. The red dashed contour shows the area of slip larger than two metres⁴¹. The bottom left inset shows the location of the main study area. The top right inset depicts an enlargement of the Upper Trisuli valley. Seismic stations belong to the network operated by NSC (Nepal) and DASE (France): Gorkha (GKN), Daman (DMN), Kakani (KKN), Phulchoki (PKI), Gumba (GUN) and Jiri (JIRN). Global Positioning System (GPS) stations belong to the Nepal Geodetic network operated by Caltech (USA), DMG (Nepal) and DASE (France): Lamjung (LMJG), Daman (DMN), Chilime (CHLM), Kakani (KKN) and Gumba (GUN). Details of the main CO₂ emission sites and of the hot spring sites as well as their respective responses to the earthquake are given in Table 1 and Supplementary Table 1 and in Table 2 and Supplementary Table 2, respectively. Figure performed using Generic Mapping Tool (http://gmt.soest.hawaii.edu)
<table><caption>Table 1 Post-seismic effects on the carbon dioxide emissions at five hydrothermal sites in Central Nepal. Data of CO₂ flux, total CO₂ discharge, and carbon isotopic ratio of CO₂ emissions are compiled before and after the Gorkha earthquake</caption><thead><tr><th rowspan="2">Site</th><th rowspan="2">Location</th><th rowspan="2">Before/after Gorkha earthquake</th><th colspan="3">CO₂ flux (g m<sup>-2</sup> d<sup>-1</sup>)</th><th rowspan="2">CO₂ discharge (10<sup>-3</sup> mol s<sup>-1</sup>)</th><th colspan="2">δ<sup>13</sup>C of CO₂ (V-PDB) (%)</th><th rowspan="2">Post-seismic effect</th></tr><tr><th>N<sub>meas</sub> (n<sub>point</sub>)</th><th>Min/max range</th><th>Geometric mean</th><th>N<sub>meas</sub> (n<sub>point</sub>)</th><th>Geometric mean</th></tr></thead><tbody><tr><td colspan="10"><em>Budhi Gandaki Valley</em></td></tr><tr><td>Machhakhola</td><td>Near Hot Spring</td><td>After</td><td>63(43)</td><td>21.6/81,800</td><td>514 ± 16</td><td>>1240 ± 330</td><td>4(4)</td><td>-3.2 ± 0.2</td><td>Increase</td></tr><tr><td colspan="10"><em>Upper Trisuli Valley</em></td></tr><tr><td rowspan="3">Sanjen</td><td>Tunnel</td><td>After</td><td>36(33)</td><td>13.2/94,700</td><td>640 ± 40</td><td>>580 ± 150</td><td>4(4)</td><td>-0.1 ± 0.2</td><td>Increase</td></tr><tr><td>Piezometer DH1</td><td>After</td><td>3(1)</td><td></td><td>92,900 ± 5200</td><td>n.e.</td><td>3(1)</td><td>-0.94 ± 0.01</td><td>Increase</td></tr><tr><td>Piezometer DH2</td><td>After</td><td>3(1)</td><td></td><td>15,400 ± 2900</td><td>n.e.</td><td>4(1)</td><td>-0.70 ± 0.01</td><td>Increase</td></tr><tr><td rowspan="3">Chilime</td><td rowspan="3">Gas Zone</td><td>Before<sup>a,c</sup></td><td>192 (92)</td><td>3.4/123,000</td><td>608 ± 12</td><td>470 ± 120</td><td>2(2)</td><td>-1.5 ± 0.1</td><td>Decrease</td></tr><tr><td>After (01/2016)</td><td>138 (131)</td><td>2.6/10,100</td><td>136 ± 2</td><td>91 ± 19</td><td>3(3)</td><td>-1.50 ± 0.06</td><td></td></tr><tr><td>After (01/2018)</td><td>99(93)</td><td>0.3/2620</td><td>21.3 ± 2.0</td><td>19 ± 4</td><td>0</td><td>n.m.</td><td></td></tr><tr><td rowspan="3">Syabru-Bensi</td><td>GZ1-2</td><td>Before<sup>a,d,e</sup></td><td>652 (333)</td><td>2.5/236,000</td><td>196 ± 2</td><td>480 ± 50</td><td>7(3)</td><td>-0.88 ± 0.07</td><td>Increase</td></tr><tr><td>Terrace</td><td>(2006-2011)</td><td>577 (378)</td><td>2.8/226,000</td><td>236 ± 3</td><td>1010 ± 110</td><td>14(4)</td><td>-0.74 ± 0.02</td><td></td></tr><tr><td>After (2015-2018)</td><td>0</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="3">Timure</td><td>Hot Spring FF2</td><td>After</td><td>69(44)</td><td>73/5600</td><td>767 ± 22</td><td>91 ± 23</td><td>1(1)</td><td>-1.6 ± 0.1</td><td>Increase</td></tr><tr><td rowspan="2">Northern Profile</td><td>Before<sup>c</sup></td><td>75(69)</td><td>110/17,400</td><td>1850 ± 50</td><td>270 ± 70</td><td>2(2)</td><td>-2.11 ± 0.02</td><td>Increase</td></tr><tr><td>After (01/2016)</td><td>117(89)</td><td>6.1/175,000</td><td>2170 ± 50</td><td>740 ± 200</td><td>0</td><td>n.m.</td><td></td></tr></tbody></table>

n.e.: not estimated, n.m.: not measured, N<sub>meas</sub>(n<sub>point</sub>): total number of measurements (total number of measurement points)

Data from: a ref.30, b ref.52, c ref.51, d ref.34, e ref.33 and 62

MCT zone (Fig. 1), releasing only partially the stored elastic energy sufficient to produce a possible M<sub>w</sub>~9 mega-earthquake<sup>42</sup>. The mainshock was followed by an intense aftershock sequence<sup>43</sup>, with >90 (M<sub>L</sub> > 5) events over a year. From 25 April 2015 to 31 December 2017, >10,200 (M<sub>L</sub> > 3), >1100 (M<sub>L</sub> > 4) and >135 (M<sub>L</sub> > 5) earthquakes occurred in Central Nepal (Fig. 1), leading to the most important seismic crisis in Nepal since 1934. The Gorkha earthquake is the first large Himalayan earthquake with seismic<sup>43</sup> and Global Positioning System<sup>44,45</sup> (GPS) data, as well as prior data on hydrothermal systems and CO₂ emissions<sup>30</sup>.

In this paper, we report spectacular outbursts of CO₂ and hydrothermal disturbances at several sites in Central Nepal separated by >10 kilometres that have been triggered by the 2015 Gorkha earthquake. These observations, in particular the first earthquake-induced gaseous changes in the absence of magmatic activity, reveal high sensitivity of Himalayan hydrothermal systems to co-, post- and possibly pre-seismic deformation, and non-stationary release of metamorphic CO₂ from active orogens.

## Results

**Outbursts of CO₂ triggered by the Gorkha earthquake.** New CO₂ emissions were observed after the Gorkha earthquake, sometimes spatially associated with substantial changes from the pre-existing CO₂ emissions (Table 1 and Supplementary Table 1). Because a large number of CO₂ flux measurements (see Methods) were made before (n = 1718; December 2006 – January 2011) and after the earthquake (n = 1053; November 2015 – January 2018), a quantitative comparison of CO₂ emissions before and after the earthquake can be undertaken in the Upper Trisuli valley. Spectacular effects were observed first in Syabru-Bensi, 56 kilometres east from the epicentre (Fig. 1), a site that had been extensively studied from 2006 to 2011. The main gas emission zone, located above the hot springs (Fig. 2; Supplementary Fig. 1), on the west bank of the Trisuli river, changed substantially. Significant CO₂ fluxes and discharges appeared near the dwellings above the previously existing gas zone, in a cultivated area where flux level

was close to local background levels before 2009. Hydrogen sulphide is now frequently smelled inside the houses. Total CO₂ discharge in the area now reaches 1010 ± 110 mmol s⁻¹ (3.8 ± 0.4 ton d⁻¹), corresponding to an increase factor of 2.1 ± 0.3 (Fig. 2). Fluxes at two reference locations (K+6 and K+12; Fig. 2), regularly monitored since November 2015 (Fig. 3), peaked in November 2016 and, in January 2018, continued to remain 2-order of magnitude higher than their pre-seismic values. By contrast, the cavity characterised by the largest pre-earthquake CO₂ fluxes (>10⁵ g m⁻² d⁻¹) showed one-order-of-magnitude smaller CO₂ emission after the earthquake (Supplementary Fig. 4). The cavity fluxes returned to pre-earthquake values in January 2018, more than 2.7 years after the mainshock.

New CO₂ emissions and changes in pre-existing CO₂ emissions were also observed at other locations in the same valley (Supplementary Fig. 5). In Timure, nine kilometres to the north, CO₂ emissions tripled after the earthquake (in November 2015 and January 2016) along a profile which had been measured precisely at different times before 2011 (Supplementary Fig. 6). Repeated measurements in September 2017 yielded a post-seismic CO₂ emission increase of a factor of 8 ± 3 compared with pre-earthquake values. By contrast, at the Chilime hot spring site, seven kilometres to the west of Timure (Fig. 1), while CO₂ fluxes remain significant, the total discharge is reduced by a factor of 5.2 ± 1.7 in January 2016, and by a factor of 25 ± 8 in January 2018. The CO₂ emissions in Timure and Chilime thus depict opposite and unstabilized post-seismic responses > 2.7 years after the mainshock.

Two kilometres to the west of Chilime, at the Sanjen hydropower construction site (Fig. 1), spectacular bubbling was suddenly noticed, beginning of November 2015, in two 40-metre-deep piezometers whose water level was being monitored. When construction resumed during 2016, CO₂ emissions appeared at several locations in a tunnel being excavated nearby (Supplementary Fig. 7 and Movie 1). The CO₂ emission in the tunnel was still present in January 2018. The CO₂ concentration in the air of the tunnel then ranged from 4 to
<figure><img src="image_5.png" /></figure>

**Fig. 2** Outburst of carbon dioxide emissions following the 2015 Gorkha earthquake in Syabru-Bensi, Upper Trisuli valley. Interpolated maps of surface CO₂ fluxes above the main hot springs are shown: **a** before the Gorkha earthquake (652 measurements³³,³⁴ from 2006 to 2011); **b** after the earthquake (577 measurements from November 2015 to January 2018). CO₂ fluxes are expressed in logarithmic scale. Topography was determined in 2007 by geodetic levelling⁶². Numbers are carbon isotopic ratios of the gaseous CO₂ (δ¹³C, relative to V-PDB). Total CO₂ discharge amounts to 480 ± 50 mmol s⁻¹ (1.8 ± 0.2 ton d⁻¹) before the earthquake and 1010 ± 110 mmol s⁻¹ (3.8 ± 0.4 ton d⁻¹) after, giving a post-seismic overall increase of a factor of 2.1 ± 0.3. The surface area of high CO₂ emission ( > 500 g m⁻² d⁻¹) increased by a factor of about 1.9, with the barycentre shifting eastward by 15 ± 3 metres. Example of CO₂ fluxes measured along the K profile in the new high emission zone is shown in Supplementary Fig. 2. Radon-222 fluxes were also measured on the terrace and were compared with data obtained before the earthquake. No significant change was found in the CO₂-radon fluxes relationship, suggesting similar gas transport mechanism, source, and travel time before and after the earthquake (Supplementary Fig. 3). While no measurement was performed before and during the earthquake, locals reported the unusual simultaneous death of hen and chicken possibly weeks before the earthquake near house H. Figure performed using PV-WAVE® software (Rogue Wave)

5 vol%, thus creating a major health hazard, and the total CO₂ discharge was estimated in the tunnel to 580 ± 150 mmol s⁻¹ (2.2 ± 0.6 ton d⁻¹), hence of the same order as the main Syabru-Bensi discharge.

In the Budhi Gandaki valley, 16.5 kilometres to the epicentre (Fig. 1), the highest CO₂ emissions ever reported in Central Nepal were observed in January 2017 on the partly flooded riverbank, where no phenomenon had been known to the locals before the earthquake (Supplementary Movie 2). The CO₂ emission was still present in January 2018, with a total discharge higher than 1240 ± 330 mmol s⁻¹ (4.7 ± 1.3 ton d⁻¹), similar to the whole CO₂ discharge observed in Syabru-Bensi. In addition to the flux from the bank, an innovative method was used to measure the CO₂ flux from and through the hot water pond (see Methods and Supplementary Fig. 8). This is the most spectacular new CO₂ emission observed so far in the Himalayas and elsewhere in the absence of volcanic activity. Together with the Upper Trisuli

valley (Syabru-Bensi, Sanjen and Timure), such a persistent post-seismic CO₂ outburst is unusual. For instance, small outbursts following the 2008 Wenchuan earthquake in China lasted only a few months⁴⁶. In the Marsyandi valley, further west (Fig. 1), by contrast, no comparable phenomenon was observed and peak CO₂ fluxes (5300 to 28,700 g m⁻² d⁻¹), while significant, are smaller than in the Upper Trisuli valley.

All these CO₂ emissions were quantified in the field during the dry season or in absence of rain, as shown by the 2015–2017 rainfall data in Dhunche (seven kilometres southwest to Syabru-Bensi) and Timure (Fig. 3). At a given site, the CO₂ flux data were obtained before and after the earthquake at about the same periods under similar meteorological conditions. Besides, the observed changes in the CO₂ emission were done contemporary at several sites separated by more than 10 kilometres. All these care and observations preclude to these changes any environmental, meteorological or shallow origin.
<figure><img src="image_6.png" /><figcaption>**Fig. 3** Temporal changes in selected parameters in Syabru-Bensi, Upper Trisuli valley. a Corrected northward position difference of GPS stations CHLM minus KKN4 (red), cumulated number of aftershocks with vertical Peak Ground Velocity (PGV) higher than 1 cm s⁻¹ in Syabru-Bensi (blue) (Supplementary Fig. 12), Seismic Moment Release Rate (SMRR) (orange) and local magnitude (green and purple) of earthquakes ≤ 25 kilometres from Syabru-Bensi. The CHLM minus KKN4 difference of the vertical position helps cancelling most of the seasonal variations. Aftershock epicentres from 25 April 2015 to 31 December 2017 (NSC) are taken into account. b Water temperature (red) and dissolved inorganic carbon (DIC) concentration (orange) of the main hot spring of Syabru-Bensi (SBPO) from 2001 to 2018, surface CO₂ fluxes (blue) at two locations on the alluvial terrace above the hot spring, K+6 and K+12 (Fig. 2 and Supplementary Fig. 2), and the 2015–2017 rainfall data in Dhunche, seven kilometres southwest to Syabru-Bensi (green), and Timure, nine kilometres to the north (purple) (courtesy of Nepal Department of Meteorology, Kathmandu). Pre-earthquake water temperature and DIC concentration of SBPO are compiled from our database and the literature³⁰. Timing of the Gorkha earthquake, the Chilime hot spring cessation and the Sanjen CO₂ outburst are displayed as vertical grey lines. Possible pre-seismic water-cooling of 0.6 ± 0.2 °C is observed 2 weeks before the earthquake, followed by a progressive warming reaching 4.4 ± 0.2 °C about 1 year after the earthquake, which persisted >2 years after the mainshock. Water temperature returned to pre-earthquake values in January 2018, more than 2.7 years after. Water-cooling was observed a few weeks before the earthquake at three other hot springs of the Upper Trisuli valley (Supplementary Fig. 9). The DIC concentration increased by 29 ± 2% after the mainshock and returned to pre-earthquake values in September 2017. Compared with CO₂ flux values obtained in 2006 (ref. ³³), close to local background level, the post-seismic values were significantly larger in November 2015 and January 2016 (increase factors of 150 ± 30 and 160 ± 20 for K+6 and K+12, respectively). The CO₂ fluxes peaked in November 2016 and continued to show higher values in January 2018 compared with 2006</figcaption></figure>

Carbon isotopic ratio of the CO₂ emissions (δ¹³C; see Methods) was systematically measured at all occurrences (Table 1 and Supplementary Table 1), giving the most comprehensive data set (n = 77) so far in a seismically active area without magmatic activity. Along the rupture zone, average δ¹³C values range from -6.9 ‰ to -0.1 ‰, and appear similar at the various sites of a given valley. CO₂ emissions in Sanjen and Syabru-Bensi, although separated by eight kilometres, have comparable CO₂ concentration (96–98%) and δ¹³C signature (from -0.9 ‰ to -0.7 ‰), precluding shallow sources, but instead suggesting similar CO₂ source and transport from a crustal-scale reservoir. In addition, the δ¹³C values remain relatively stable after the earthquake (Supplementary Fig. 4), indicating that the earthquake revealed a pre-existing CO₂ reservoir.

**Hydrothermal unrest triggered by the Gorkha earthquake.** The CO₂ emissions were associated with changes of hydrothermal activity (Table 2 and Supplementary Table 2). In Syabru-Bensi,

new hot springs appeared after the earthquake, with a few persisting in January 2018, > 2.7 years after the mainshock. The temperature of the main Syabru-Bensi hot spring (SBP0), which had been stable at 60.7 ± 0.1 °C for more than 12 years before the earthquake (Fig. 3), increased to 64.1 ± 0.3 °C after the earthquake (November 2015 – January 2016), along with a significant flow rate increase of 16 ± 1%. Water temperature peaked about 1 year after the mainshock. Water warming of SBP0 persisted more than 2 years after the mainshock, before starting to decrease to pre-earthquake values in January 2018, more than 2.7 years after. Other springs in Syabru-Bensi (SBB5) also showed warming, persisting in January 2018 (Supplementary Fig. 9). A slight pre-seismic water-cooling was possibly detected in Syabru-Bensi (Fig. 3). More significantly, three other hot springs of the same valley (Timure, Langtang and Chilime) also showed cooling a few weeks before the earthquake (Supplementary Fig. 9). In Tatopani (Budhi Gandaki valley), locals reported hot spring temperature decrease a few weeks before the earthquake. While co-, post- and
<table><caption>Table 2 Post-seismic effects on the hot springs at nine hydrothermal sites in Central Nepal. Data of temperature, flow rate, and dissolved inorganic carbon concentration and isotopic ratio of hot springs are compiled before and after the Gorkha earthquake</caption><thead><tr><th rowspan="2">Site</th><th rowspan="2">Location</th><th rowspan="2">Name</th><th rowspan="2">Type</th><th rowspan="2">Before/after Gorkha earthquake</th><th rowspan="2">Spring temperature (°C)</th><th rowspan="2">Spring flow rate (L s<sup>-1</sup>)</th><th colspan="3">Dissolved inorganic carbon (DIC)</th><th rowspan="2">Post-seismic effect</th></tr><tr><th>N<sub>meas</sub></th><th>C<sub>DIC</sub> (mmol L<sup>-1</sup>)</th><th>δ<sup>13</sup>C<sub>DIC</sub> (%)<br>(V-PDB)</th></tr></thead><tbody><tr><td colspan="12" style="text-align: left; font-weight: bold;">Budhi Gandaki Valley</td></tr><tr><td>Machhakhola</td><td>Eastern Bank</td><td>BUDO</td><td>HS<sup>a</sup>,<br>BB,<br>DDS</td><td>After</td><td>59.8 ± 0.2</td><td>n.m.</td><td>1</td><td>26.4 ± 0.3</td><td>-0.5 ± 0.2</td><td>Increase (New)</td></tr><tr><td rowspan="3">Tatopani</td><td rowspan="3">Southern<br>Secondary<br>Springs<br>Main Springs</td><td>BUD3B</td><td>HS<sup>a</sup></td><td>After</td><td>43.2 ± 0.1</td><td>n.m.</td><td>1</td><td>25.6 ± 0.9</td><td>0.9 ± 0.3</td><td>Increase (New)</td></tr><tr><td>BUD4B</td><td>HS</td><td>Before<sup>b</sup><br>After</td><td>30 ± 1<br>22.2 ± 1.8</td><td>n.m.<br>0.061 ± 0.002</td><td>1<br>1</td><td>12.2 ± 1.2<br>15.8 ± 0.3</td><td>n.m.<br>3.7 ± 0.3</td><td>Decrease (Temp.)</td></tr><tr><td>BUD4C</td><td>HS</td><td>Before<sup>b</sup><br>After</td><td>50 ± 1<br>48.5 ± 1.3</td><td>n.m.<br>0.19 ± 0.04</td><td>1<br>1</td><td>19.1 ± 1.9<br>22.0 ± 0.2</td><td>3.3 ± 0.3<br>1.6 ± 0.3</td><td>Decrease (Temp.)</td></tr><tr><td colspan="12" style="text-align: left; font-weight: bold;">Upper Trisuli Valley</td></tr><tr><td rowspan="3">Sanjen</td><td>Tunnel</td><td>TSJ3</td><td>HS<sup>a</sup>,<br>DDS</td><td>After</td><td>20.3 ± 0.1</td><td>n.m.</td><td>1</td><td>&gt;55</td><td>-9.1 ± 0.7</td><td>Increase (New)</td></tr><tr><td>Piezometer</td><td>DH1</td><td>AQ,<br>DDS</td><td>After</td><td>18.9 ± 0.6</td><td>n.m.</td><td>2</td><td>43.8 ± 0.1</td><td>1.6 ± 0.1</td><td>Unknown</td></tr><tr><td>Piezometer</td><td>DH2</td><td>AQ,<br>DDS</td><td>After</td><td>19.1 ± 0.8</td><td>n.m.</td><td>1</td><td>43.3 ± 0.1</td><td>2.1 ± 0.1</td><td>Unknown</td></tr><tr><td rowspan="3">Chilime</td><td rowspan="3">Hot Spring</td><td>CHI</td><td>HS,<br>DDS</td><td>Before<sup>b,c,d,e,f,g</sup><br>After</td><td>48.9 ± 0.4<br>No spring</td><td>5.0 ± 0.2</td><td>3</td><td>13.8 ± 1.3</td><td>8.3 ± 0.4</td><td>Decrease<br>(Cessation)</td></tr><tr><td>GZ3</td><td>HS,<br>DDS</td><td>Before<br>After</td><td>20.3 ± 1.2<br>23.4 ± 0.1</td><td>n.m.<br>n.m.</td><td>1<br>1</td><td>1.9 ± 0.1<br>3.2 ± 0.1</td><td>n.m.<br>3.7 ± 0.1</td><td>Increase (Temp.)</td></tr><tr><td>FF2</td><td>HS<sup>a</sup>,<br>BB</td><td>After</td><td>In Trisuli River</td><td>n.m.</td><td>0</td><td>0</td><td>n.m.</td><td>Increase (New)</td></tr><tr><td rowspan="6">Syabru-Bensi</td><td rowspan="3">Western Bank<br>Secondary<br>Springs</td><td>SBPO</td><td>HS,<br>DDS</td><td>Before<sup>b,c,d,g</sup><br>After</td><td>60.7 ± 0.1<br>64.1 ± 0.3</td><td>0.087 ±<br>0.004</td><td>4</td><td>25.8 ± 1.6</td><td>4.7 ± 0.7</td><td>Increase (Temp.,<br>flow, C<sub>DIC</sub>, δ<sup>13</sup>C<sub>DIC</sub>)</td></tr><tr><td>GZ3</td><td>HS,<br>DDS</td><td>Before<sup>g</sup><br>After</td><td>31.8 ± 0.3<br>34.9 ± 0.2</td><td>0.104 ± 0.007<br>0.282 ± 0.009</td><td>5</td><td>34.3 ± 1.5<br>25.8 ± 3.6</td><td>1.0 ± 0.1<br>0.9 ± 0.1</td><td>Increase (Temp.,<br>flow)</td></tr><tr><td>FF2</td><td>HS<sup>a</sup>,<br>BB</td><td>Before<sup>g</sup><br>After</td><td>50.1 ± 1.9<br>40.2 ± 0.5</td><td>n.m.<br>n.m.</td><td>2<br>1</td><td>29.6 ± 0.9<br>17.8 ± 1.1</td><td>0.2 ± 0.3<br>2.6 ± 0.1</td><td>Decrease (Temp.)</td></tr><tr><td rowspan="3">Main Springs</td><td>SBM</td><td>HS<sup>a</sup>,<br>DDS</td><td>After</td><td>34.9 ± 0.2</td><td>n.m.</td><td>1</td><td>22.3 ± 0.1<br>30.3 ± 0.1</td><td>2.1 ± 0.1<br>-0.6 ± 0.1</td><td>Increase (New)</td></tr><tr><td>SBN</td><td>HS<sup>a</sup>,<br>DDS</td><td>After</td><td>37.3 ± 0.8</td><td>n.m.</td><td>1</td><td>27.4 ± 0.1</td><td>2.3 ± 0.1</td><td>Increase (New)</td></tr><tr><td>SBN2</td><td>HS<sup>a</sup>,<br>DDS</td><td>After</td><td>39.2 ± 0.6</td><td>0.50 ± 0.02</td><td>3</td><td>17.3 ± 1.1</td><td>2.0 ± 0.2</td><td>Increase (New)</td></tr><tr><td rowspan="3">North Syabru</td><td rowspan="3">Western Bank</td><td>TT1</td><td>HS,<br>BB</td><td>Before<sup>b,c</sup><br>After</td><td>24.38 ± 0.07<br>25.3 ± 0.1</td><td>n.m.<br>n.m.</td><td>4<br>2</td><td>35.8 ± 2.9<br>38.3 ± 2.3</td><td>12.3 ± 0.7<br>-0.6 ± 0.1</td><td>Increase (Temp.)</td></tr><tr><td>TIM</td><td>HS,<br>DDS</td><td>Before<sup>b,c,f</sup><br>After</td><td>61.3 ± 2.2<br>70.7 ± 1.0</td><td>0.20 ± 0.02<br>n.m.</td><td>5<br>2</td><td>17.3 ± 1.0<br>14.2 ± 0.2</td><td>3.3 ± 0.1<br>2.1 ± 0.1</td><td>Increase (Temp.)</td></tr><tr><td>LPAH</td><td>HS</td><td>Before<sup>c,e</sup><br>After</td><td>41.0 ± 0.1<br>37.3 ± 0.5</td><td>n.m.<br>n.m.</td><td>2<br>1</td><td>5.1 ± 0.2<br>4.8 ± 0.2</td><td>-4.8 ± 0.1<br>-3.2 ± 0.2</td><td>Decrease (Temp.)</td></tr><tr><td colspan="12" style="text-align: left; font-weight: bold;">Bhote Koshi Valley</td></tr><tr><td rowspan="2">Kodari</td><td rowspan="2"></td><td>KOD</td><td>HS</td><td>Before<sup>b,g</sup><br>After</td><td>44.6 ± 1.2<br>50.7 ± 0.1</td><td>3.0 ± 0.2<br>3.00 ± 0.04</td><td>1<br>0</td><td>8.4 ± 0.3<br>n.m.</td><td>-8.7 ± 0.1<br>n.m.</td><td>Increase (Temp.)</td></tr><tr><td>KOD2</td><td>HS<sup>a</sup></td><td>After</td><td>48.1 ± 0.1</td><td>n.m.</td><td>0</td><td>n.m.</td><td>n.m.</td><td>Increase (New)</td></tr></tbody></table>

n.m.: not measured, HS: hot spring, BB: bubbles, DDS: diffuse degassing structure, AQ: aquifer degassing

<sup>a</sup>New spring that appeared after the Gorkha earthquake

Compilation of new original data and of data compiled in ref.³⁰, and in particular from: <sup>b</sup>ref.³²; <sup>c</sup>refs.³¹ and ⁵²; <sup>d</sup>ref.⁶³; <sup>e</sup>ref.⁶⁴; <sup>f</sup>ref.⁵¹; <sup>g</sup>ref.³³

sometimes pre-seismic changes in hot spring temperature have been reported in the literature¹⁶,¹⁹, a two-year-long or longer warming is unusual.

The most spectacular change in hydrothermal activity, however, is the complete cessation of the Chilime hot spring at the end of October 2015, after periods of unusual intermittence between April and June 2015 (Supplementary Fig. 10). Before the earthquake, this spring had a stable flow rate of > 5 L s⁻¹; it was the pillar of local economy, being the most important in Central Nepal after the Kodari hot spring (Fig. 1). The village elders had previously reported spring intermittence at the time of the 1934 earthquake, but no cessation. Particularly impressive is the fact that it happened at about the time of the CO₂ outburst in Sanjen.

Dissolved Inorganic Carbon (DIC) concentration and isotopic ratio (C<sub>DIC</sub> and δ<sup>13</sup>C<sub>DIC</sub>; see Methods) were systematically measured for all hot spring waters (Table 2 and Supplementary Table 2). In Syabru-Bensi, where the number of measurements

is significant (n = 8), C<sub>DIC</sub> of the SBP0 hot spring increased by 29 ± 2% after the mainshock and, in September 2017, returned to values measured before the mainshock (Fig. 3), while Ca and Na concentrations remained similar within 5%. The δ<sup>13</sup>C<sub>DIC</sub> decreased also significantly (Supplementary Fig. 11). These observations suggest a larger amount of dissolved carbon in the SBP0 hot spring after the earthquake, which is compatible with the aforementioned increase in gaseous CO₂ emissions. Available pH values (Supplementary Fig. 9) also show a slow return to pre-earthquake values, with anomalously high values before the earthquake, a fact to be interpreted with caution given the lack of additional geochemical information. In Syabru-Bensi, dissolved radon and radium concentrations in SBP0 and SBB5 hot springs changed after the earthquake (Supplementary Fig. 9). In January 2018, >2.7 years after the mainshock, the warmest spring in Syabru-Bensi (SBP0) returned to pre-earthquake conditions, while those more dependent on superficial effects (e.g., SBB5) did not yet.
<figure><img src="image_7.png" /><figcaption>**Fig. 4** Conceptual models of the carbon dioxide and hydrothermal transport before and after the Gorkha earthquake. Pre- and post- seismic models are shown separately for **a-b** the Bhote Kosi valley (Syabru-Bensi and Timure) and **c-d** the Chilime valley (Chilime, Sanjen and Brapche). The CO₂ is produced at depth and a significant fraction is transported to the surface by hydrothermal circulations, which efficiently take advantage of fault network; some CO₂ can also move upwards independently, reaching the surface without hot spring (Sanjen, Syabru-Bensi GZ3) and accumulate in the subsurface reservoir revealed by CO₂ emissions following the earthquake. The main effect of the earthquake is an increase of vertical permeability with transient or permanent changes</figcaption></figure>

**First non-volcanic earthquake-induced gaseous changes.** Hydrologic responses to seismic stimulation have been observed in numerous instances¹⁰, and correlated with Seismic Energy Density (SED) or Peak Ground Velocity (PGV) values. In Lassen (Cascades Range, USA) for example, a 2014 volcano-seismic swarm peaking at *M*<sub>w</sub> = 3.85 at 5.7 kilometres distance, corresponding to SED~0.2 J m<sup>-3</sup> and PGV~0.2 cm s<sup>-1</sup>, caused an outburst of geothermal fluids, explained by a two-fold permeability increase²². In the case of the Matsushiro earthquake swarm in Japan, modelling indicates the necessity of a 2-order-of-magnitude permeability increase with a small overpressure of only a few megapascals⁹. In Syabru-Bensi, the Gorkha earthquake produced SED~63 J m<sup>-3</sup> and PGV~26 cm s<sup>-1</sup> (see Methods; Supplementary Table 3), hence strong enough to affect the hydrothermal system, with several aftershocks maintaining high ground motion during the following 31 months (Supplementary Fig. 12). These estimates are confirmed by GPS time-series of the Chilime station44,45 (Fig. 1), giving PGV~49 cm s<sup>-1</sup> (Supplementary Fig. 13). While the first clearly documented instances in the Himalayas, our observed near-field CO<sub>2</sub> outbursts and hydrothermal unrest remain compatible with previously compiled earthquake-induced changes (Supplementary Fig. 14). Our observed changes in the CO<sub>2</sub> emissions are however the first earthquake-induced gaseous changes in a non-volcanic region.

**Discussion**

The observations following the Gorkha earthquake give some indications that the standard model, where CO₂ is degassed from hydrothermal waters, needs some modifications at least in part. Indeed, previous work in Syabru-Bensi³⁴ identified large CO₂ emissions without the presence of an important hot spring in the vicinity. Now, the evidence for an extended reservoir of crustal

CO₂ is overwhelming, which we accommodate with the following conceptual models (Fig. 4). In the case of the Trisuli valley (Syabru-Bensi and Timure) (Fig. 4a, b), hydrothermal circulations are important and increased after the earthquake. The main CO₂ source can therefore be the hydrothermal circulations, with near surface degassing and subsurface accumulation in the whole fault zone. Possible additional CO₂ can also be released in the footwall directly from the production source. The earthquake caused (Fig. 4b) an increasing water discharge, with increased CO₂ degassing, or better communication to the surface of previously accumulated CO₂. In the case of the Chilime valley (Chilime, Sanjen and Brapche) (Fig. 4c, d), however, the pervading presence of CO₂, revealed by the earthquake, kilometres away from hot springs, is better explained if CO₂ is directly accumulated from below, possibly through the water table, without advection by hydrothermal circulations.

Strikingly, CO₂ emissions with similar isotopic anomaly were observed over the whole region, whenever faults, boreholes or a tunnel gave the opportunity. These observations also attest the presence of a large, relatively shallow, reservoir of CO₂ in the Himalayan crust, suggesting that metamorphic CO₂ produced at depth is huge, as independently shown by petrological estimates47,48, and unlikely sequestered. Nevertheless, as hot water could also be found when deep boreholes are available, the overwhelming presence of hot water could match the evidence of an extended CO₂ reservoir, and the debate cannot be considered as closed.

The effects of the earthquake on our hydrothermal systems are not unexpected. Indeed, our observed CO₂ emissions appear consistent with hydrothermal outburst effects, as observed in Lassen²², suggesting that apparent two-fold permeability increases and/or changes in hydrothermal pathways for pre-existing CO₂ emissions could have been major effects of the Gorkha
earthquake. First-order modelling indeed confirms that at least a doubling of permeability could explain CO₂ flux increases (Supplementary Fig. 15). In addition, first-order modelling supports that a lasting permeability increase¹⁸ of 10–20% could accommodate the Syabru-Bensi hot spring warming and increased flow (Supplementary Figs. 16 and 17). These vertical permeability increases appear compatible with post-seismic afterslip⁴⁵. The global pre-seismic water-cooling could be evidence for a pre-seismic dilatancy effect⁴⁹ of the impending rupture zone.

Compared with reported earthquake-induced hydrological effects¹⁰,¹⁶,¹⁹,⁴⁶, some of our observations are unusual: first, the six-month delay of the Sanjen CO₂ outburst and the shut down of the Chilime hot spring, and second, the persistence (and rising at some places) > 2.7 years after in January 2018 of all the increased CO₂ emissions and hydrothermal unrest. In Lassen²², the permeability increase was inferred to be maintained for 50–60 days. In the Napa valley¹³, a vertical permeability change, following diffusion of a pressure pulse, was proposed to explain post-seismic stream discharge, suggesting basin-scale fluid diffusivity of 0.1 m² s⁻¹. Taking this conservative value, for one kilometre spatial scale, the timescale for fluid displacement is ~ 120 days, compatible with our six-month delay time. This suggests that the Sanjen CO₂ outburst and the Chilime cessation likely result from hydrogeological fluid diffusion.

Other explanations are, however, possible. In Central Nepal, post-seismic deformation (afterslip on the MHT) and associated aftershocks remained active several months after the mainshock. About 6 months after the mainshock, contemporaneous with the Sanjen outburst and the Chilime cessation (Fig. 3a), the number of seismic events producing PGV>1 cm s⁻¹ in Syabru-Bensi increased, and a linear trend in the GPS time-series, continuing to October 2017 about 2.5 years after the mainshock, was initiated. Thus, some post-seismic CO₂ emissions and hydrothermal unrest may be related to changes in the state of post-seismic relaxation. Incidentally, the time of the Sanjen outburst and the Chilime cessation coincides with renewal of aftershock activity near Syabru-Bensi in the fall of 2015 (Fig. 3a), including an $M_L$ = 5.3 event a few weeks before. In addition, two $M_L$ > 4 events occurred at about the same time (end of October – beginning of November), less than 15 kilometres from Syabru-Bensi, producing PGV > 1 cm s⁻¹ (Supplementary Table 3).

The Himalayan hydrothermal systems appear highly sensitive to small deformation rates and are therefore in near-critical condition. This suggests that post-seismic relaxation of co-seismic stress may result from pore pressure changes, and that metamorphic CO₂ may in turn play a role in the installation of the next inter-seismic regime. Alternatively, the still evolving CO₂ emissions in Syabru-Bensi, Timure and Chilime may also indicate a currently unstabilized system, able to diverge unpredictably. Given the issue of a pending mega-earthquake in the region⁴², long-term monitoring of CO₂ emissions should seriously be considered as a chance to capture possible pre-earthquake signals.

In this paper, we have presented the first assessment of CO₂ emissions triggered by a major earthquake, demonstrating the coupling between mechanical deformation and fluid transport properties at the crustal scale, highlighting that crustal deformation dynamically affects permeability during the seismic cycle. The large post-seismic CO₂ emissions observed in the Narayani basin suggest non-stationary metamorphic CO₂ production, and that its current estimate³⁰,³² (>1.3 × 10¹⁰ mol year⁻¹ with >(1.0 ± 0.2) × 10⁸ mol year⁻¹ from direct gaseous CO₂ emissions), independently confirmed by petrological studies⁴⁸, may be enhanced during a significant fraction of the seismic cycle. Metamorphic CO₂ and its transport therefore emerge as an essential component of mountain build-up and the associated dynamics of large earthquakes.

## Methods

**Carbon dioxide flux measurement and mapping.** The accumulation chamber method⁵⁰ was used to measure surface CO₂ flux and to quantify the total CO₂ discharge of a given site and the associated uncertainties, whose assessment is based on numerous systematic tests and our 10-year experience³⁴. The method is robust, even in remote locations³¹ and during the monsoon³⁶, and allows the measurement of CO₂ fluxes over more than five orders of magnitude (Supplementary Fig. 2). The increase in CO₂ concentration in the chamber is measured using various portable infrared CO₂ sensors (Testo™ 535, Testo AG, Germany; Airwatch™ PM 1500, Geotechnical Instruments Ltd., UK; Vaisala™ CARBOCAP® Hand-Held GM70, Finland), that are regularly inter-calibrated in the laboratory. The CO₂ flux is expressed in grams per squared metre per day (g m⁻² d⁻¹). The total CO₂ discharge, expressed in mol s⁻¹ (or ton d⁻¹), is estimated using the CO₂ flux data-set by kriging and interpolation procedures³⁴. CO₂ fluxes (n = 1720) and total CO₂ discharges obtained before the Gorkha earthquake (from 2006 to 2011) are published elsewhere³⁰,³³,³⁴,⁵¹. Here we present for the first time CO₂ fluxes (n = 1668) and total discharges obtained after the earthquake in the Marsyandi, Budhi Gandaki and Upper Trisuli valleys relying on seven measurement campaigns carried out in November 2015, in January, May and November 2016, in January and September 2017, and in January 2018. Pre-earthquake and post-earthquake campaigns were performed outside the monsoon periods to reduce the meteorological effects on CO₂ flux data. Every uncertainty is given around one-sigma standard deviation (68% confidence level) and averages are arithmetic means except otherwise stated. Data are summarised in Table 1 and Supplementary Table 1.

**Carbon dioxide flux measurement through a water layer.** The bubbling CO₂ flux from or through water was measured at the new degassing site of Machkhakola (Budhi Gandaki valley) in the following manner. A collecting container was installed upside down on the water with a pipe leading to an accumulation chamber installed on the ground nearby (Supplementary Fig. 8). Then, the flux was measured in this accumulation chamber as described above. This method yielded a minimum value to the large CO₂ discharge observed in Machkhakola (Supplementary Movie 2).

**Carbon isotopic composition of the gas phase.** We sampled gas in the field using evacuated glass tubes. CO₂ fraction of the gas sample was determined micrometrically. The δ¹³C of CO₂ of the gas, expressed in per mil relative to the standard values of Vienna Pee Dee Belemnite (V-PDB), was measured after off-line purification by mass spectrometry on a Finnigan™ MAT-253 mass spectrometer (Thermo Electron Corp., Germany) in CRPG (Nancy, France)³²,³⁴. External repeatability of a given sample was ±0.1 ‰. The twenty δ¹³C values measured before the Gorkha earthquake are published elsewhere³⁰,³³,³⁴,⁵². Fifty-seven measurements were carried out after the earthquake from January 2016 to January 2018. Data are summarised in Table 1 and Supplementary Table 1.

**Water temperature, pH and flow rate measurements.** To measure water temperature of springs, we used thermometers (Generic TP101 Digital Thermometer, China) regularly inter-calibrated in the laboratory with a reference thermometer not currently used in the field (Digital Thermometer model 4400 Ertco Eutechnics, USA), and compared with high-precision (10⁻³ °C) and high-sensitivity (10⁻⁴ °C) thermometers (Seabird™ 39plus, Sea-Bird Scientific, USA). Due to the time response of the instruments, several minutes are needed to measure temperature. Experimental uncertainty of a given measurement was ±0.1 °C. The pH of the water in thermal springs was measured with various pH metres (H170 Portable pH metre, Hach, Germany; H198107 and H198130 pH metres, Hanna Instruments, USA), systematically recalibrated in the field using buffer solutions. Experimental uncertainty of a given measurement was ±0.1. The flow rate of thermal springs was determined using stopwatch and measured cylinders or buckets and was repeated at least three times. Data are summarised in Table 2 and Supplementary Table 2.

**Carbon isotopic composition and DIC concentration in water.** We sampled every water in the field using two to three glass cap vials of 12 millilitres volume each. DIC concentration ($C_{DIC}$ = $[H_2CO_3] + [HCO_3^-] + [CO_3^{2-}]$) and its isotopic composition (δ¹³C$_{DIC}$), expressed in mmol L⁻¹ and in per mil relative to V-PDB, respectively, were determined using a gas chromatograph coupled to an isotope ratio mass spectrometer (GCIRMS, GV 2003, GV Instruments, UK) in IPGP (Paris, France). The whole procedure is described elsewhere⁵³. The relative experimental uncertainty of $C_{DIC}$ was 1–2% and the experimental uncertainty of a given δ¹³C$_{DIC}$ measurement was ±0.1‰. For a given sample, final values correspond to weighted arithmetic averages of two to three measurements. The $C_{DIC}$ and δ¹³C$_{DIC}$ values measured before the Gorkha earthquake are published elsewhere³⁰⁻³³. A total of 71 measurements were carried out after the earthquake from January 2016 to January 2018. Data are summarised in Table 2 and Supplementary Table 2.

**Water radon-222 and radium-226 concentration measurements.** Dissolved radon concentration in water was measured in the field by emanometry in air⁵⁴. After air-water equilibrium is reached by manual shaking, radon concentration is inferred from scintillation flask sampling and photomultiplier counting, as described elsewhere⁵¹. Radon concentration in water is expressed in Bq L⁻¹.
Experimental uncertainty ranged from 5 to 30%. The radium concentration in water was similarly measured in the IPGP laboratory after keeping the bottle closed for at least 50–80 days⁵⁵. Expressed in mBq L⁻¹, the experimental uncertainty was the same as for radon concentration in water. Here we present only the dissolved radon and radium concentrations in the main Syabru-Bensi hot springs (Supplementary Fig. 9).

**Detection of thermal springs and CO₂ degassing areas.** The CO₂ degassing areas were detected based on pervasive hydrogen sulphide odour, measurement of high radon-222 flux (radioactive gas of half-life 3.8 days), surface temperature anomalies, occurrence of water bubbles, presence of cavities, occurrence of inactive or active travertine deposits, and discussion with local people, or a combination of the above³⁰,⁵¹. In remote places, the detection of previously unknown thermal springs relied on the use of hand-held thermal infrared cameras³⁴,³⁵,⁵¹ (model 880-V3 before 2015 and model 875-1i after, Testo™ AG, Germany).

**Determination of SED, PGV and PDS.** At the sites which experienced the most significant post-seismic changes in CO₂ emissions and hot springs, we calculated SED⁵⁶ and vertical PGV⁵⁷ produced by the Gorkha earthquake and its main aftershocks. Empirical equations are used to estimate SED⁵⁶ and PGV⁵⁸. Peak Dynamic Stress (PDS) is estimated using: PGV × shear modulus (3 × 10¹⁰ Pa) / shear wave velocity (3500 km s⁻¹). For the mainshock and six main aftershocks, values of SED, PGV and PDS at eight hydrothermal sites are gathered in Supplementary Table 3. In addition, SED, PGV and PDS produced at sites of the Upper Trisuli valley by three aftershocks located near Syabru-Bensi are also given.

**Determination of Seismic Moment Release Rate.** We calculated the Seismic Moment Release Rate (SMRR) produced by the aftershocks of the Gorkha earthquake using the NSC earthquake catalogue from 25 April 2015 to 31 December 2017. We converted local magnitude ($M_L$) given by the catalogue into moment magnitude ($M_w$), using the following linear calibration: $M_w = 1.109 M_L - 0.626$. The seismic moment was obtained using the classical relationship: $log(M_0) = 1.5M_w + 9.1$. The SMRR, expressed in N.m per 20 days, was calculated for all seismic events around Syabru-Bensi and is shown in Fig. 3.

**Modelling of earthquake-induced changes.** To study the effect of permeability changes on CO₂ flux before and after the earthquake, we used a 2-D model, described elsewhere³⁷, in which a vertical fault (f) surrounded by two media (a and b) transports the gaseous CO₂ to the surface by advection. Each medium is characterised by connected porosity ($\varepsilon$) and permeability ($k$). A pressure source is fixed at depth. Darcy's law is defined according to a pressure distribution ($p$) that follows $\Delta p^2 = 0$ (ref.⁵⁹). An example calculation is shown in Supplementary Fig. 15.

To study the effect of permeability changes on spring water flow rate at the surface, we adapted the 2-D fault model described above³⁷, to the calculation of water flow rate. We consider a vertical fault (f) driving the water flow to the surface, surrounded by two media (a and b). Each medium is characterised by permeability ($k$). Pressure source is fixed at depth. Solutions of the pressure distribution are expressed as a sum of exponential terms on horizontal axis, modulated by a sinusoidal signal along vertical axis³⁷. An example of calculation is shown in Supplementary Fig. 16.

To study the effect of water flow rate changes on hot spring exit temperature, we relied on an analytical first-order model⁶⁰ (Supplementary Fig. 17a). We consider a vertical conduit, characterised by perimeter $p$ and thermal conductivity $K_T$, that drives hot water from depth ($z=h$) to the surface ($z=0$), with flow rate $Q$, density $\rho$ and specific heat $C_p$. Temperature at depth $h$ is defined by $T(h) = \beta h + T(0)$, where $\beta$ is the thermal gradient, which is poorly constrained in this hydrothermal region⁶¹ and taken equal to 55–75 °C km⁻¹. We consider a quasi-static state having characteristic length $\lambda_T$ and we define $\theta = T-T(0)$ and $\alpha = \rho C_p Q \lambda_T / (p K_T h)$. Using the initial conditions, the differential equation $a h(d\theta/dz) + \theta = -\beta z$ has the following solution: $\theta = a h \beta (1-e^{-1/\alpha})$. Examples of calculation of the water temperature as a function of water flow rate, for two values of the thermal gradient, are shown in Supplementary Fig. 17b.

**Data availability.** The data that support the findings of this study are available in the article, in Supplementary Information, and from the corresponding author upon request (girault@ipgp.fr).

Received: 15 October 2017 Accepted: 11 June 2018

Published online: 27 July 2018

**References**

1. Irwin, W. P. & Barnes, I. Tectonic relations of carbon dioxide discharges and earthquakes. *J. Geophys. Res.* **85**, 3115–3121 (1980).

2. Miller, S. A. Chapter 1 - The role of fluids in tectonic and earthquake processes. *Adv. Geophys.* **54**, 1–46 (2013).

3. Noir, J. et al. Fluid flow triggered migration of events in the 1989 Dobi earthquake sequence of Central Afar. *Geophys. Res. Lett.* **24**, 2335–2338 (1997).

4. Miller, S. A. et al. Aftershocks driven by a high-pressure CO₂ source at depth. *Nature* **427**, 724–727 (2004).

5. Terakawa, T., Zoporowski, A., Galvan, B. & Miller, S. A. High-pressure fluid at hypocentral depths in the L'Aquila region inferred from earthquake focal mechanisms. *Geology* **38**, 995–998 (2010).

6. Waldhauser, F., Schaff, D. P., Diehl, T. & Engdahl, E. R. Splay faults imaged by fluid-driven aftershocks of the 2004 M<sub>w</sub> 9.2 Sumatra-Andaman earthquake. *Geology* **40**, 243–246 (2012).

7. Wang, Z., Huang, W., Zhao, D. & Pei, S. Mapping the Tohoku forearc: Implications for the mechanism of the 2011 East Japan earthquake (M<sub>w</sub> 9.0). *Tectonophysics* **524–525**, 147–154 (2012).

8. Sibson, R. H. An episode of fault-valve behaviour during compressional inversion? — The 2004 M<sub>L</sub>6.8 Mid-Niigata Prefecture, Japan, earthquake sequence. *Earth Planet. Sci. Lett.* **257**, 188–199 (2007).

9. Cappa, F., Rutqvist, J. & Yamamoto, K. Modeling crustal deformation and rupture processes related to upwelling of deep CO₂-rich fluids during the 1965–1967 Matsushiro earthquake swarm in Japan. *J. Geophys. Res.* **114**, B10304 (2009).

10. Wang, C.-Y. & Manga, M. Earthquakes and water. *Lect. Notes Earth Sci.* **114**, 235 (2010).

11. Muir-Wood, R. & King, G. C. P. Hydrological signatures of earthquake strain. *J. Geophys. Res.* **98**, 22035–22068 (1993).

12. Sato, T., Sakai, R., Furuya, K. & Kodama, T. Coseismic spring flow changes associated with the 1995 Kobe earthquake. *Geophys. Res. Lett.* **27**, 1219–1222 (2000).

13. Wang, C.-Y. & Manga, M. New streams and springs after the 2014 M<sub>w</sub>6.0 South Napa earthquake. *Nat. Commun.* **6**, 7597 (2015).

14. Roeloffs, E. A. Persistent water level changes in a well near Parkfield, California, due to local and distant earthquakes. *J. Geophys. Res.* **103**, 869–889 (1998).

15. Brodsky, E. E., Roeloffs, E., Woodcock, D., Gall, I. & Manga, M. A mechanism for sustained groundwater pressure changes induced by distant earthquakes. *J. Geophys. Res.* **108**, 2390 (2003).

16. Mogi, K., Mochizuki, H. & Kurokawa, Y. Temperature changes in an artesian spring at Usami in the Izu Peninsula (Japan) and their relation to earthquakes. *Tectonophysics* **159**, 95–108 (1989).

17. Ma, Z. et al. *Earthquake Prediction: Nine Major Earthquakes in China* (1966–1976) (Springer-Verlag, Berlin Heidelberg, 1990).

18. Wang, C.-Y., Manga, M., Wang, C.-H. & Chen, C.-H. Transient change in groundwater temperature after earthquakes. *Geology* **40**, 119–122 (2012).

19. Cox, S. C. et al. Changes in hot spring temperature and hydrogeology of the Alpine Fault hanging wall, New Zealand, induced by distal South Island earthquakes. *Geofluids* **15**, 216–239 (2015).

20. Manning, C. E. & Ingebritsen, S. E. Permeability of the continental crust: Implications of geothermal data and metamorphic systems. *Rev. Geophys.* **37**, 127–150 (1999).

21. Manga, M. et al. Changes in permeability caused by transient stresses: Field observations, experiments, and mechanisms. *Rev. Geophys.* **50**, RG2004 (2012).

22. Ingebritsen, S. E. et al. Hydrothermal response to a volcano-tectonic earthquake swarm, Lassen, California. *Geophys. Res. Lett.* **42**, 9223–9230 (2015).

23. Fischer, T., Matyska, C. & Heinicke, J. Earthquake-enhanced permeability – evidence from carbon dioxide release following the M<sub>L</sub> 3.5 earthquake in West Bohemia. *Earth Planet. Sci. Lett.* **460**, 60–67 (2017).

24. Avouac, J.-P. Mountain building, erosion, and the seismic cycle in the Nepal Himalaya. *Adv. Geophys.* **46**, 1–80 (2003).

25. Ader, T. et al. Convergence rate across the Nepal Himalaya and interseismic coupling on the Main Himalayan Thrust: Implications for seismic hazard. *J. Geophys. Res.* **117**, B04403 (2012).

26. Sapkota, S. N. et al. Primary surface ruptures of the great Himalayan earthquakes in 1934 and 1255. *Nat. Geosci.* **6**, 71–76 (2013).

27. Pandey, M. R., Tandukar, R. P., Avouac, J.-P., Lavfi, J. & Massot, J. P. Interseismic strain accumulation on the Himalayan crustal ramp (Nepal). *Geophys. Res. Lett.* **22**, 751–754 (1995).

28. Bollinger, L., Avouac, J.-P., Cattin, R. & Pandey, M. R. Stress buildup in the Himalaya. *J. Geophys. Res.* **109**, B11405 (2004).

29. Upreti, B. N. An overview of the stratigraphy and tectonics of the Nepal Himalaya. *J. Asian Earth Sci.* **17**, 577–606 (1999).

30. Girault, F. et al. Large-scale organization of carbon dioxide discharge in the Nepal Himalayas. *Geophys. Res. Lett.* **41**, 6358–6366 (2014).

31. Becker, J. A., Bickle, M. J., Galy, A. & Holland, T. J. B. Himalayan metamorphic CO₂ fluxes: Quantitative constraints from hydrothermal springs. *Earth Planet. Sci. Lett.* **265**, 616–629 (2008).
32. Evans, M. J., Derry, L. A. & France-Lanord, C. Degassing of metamorphic carbon dioxide from the Nepal Himalaya. *Geochem. Geophys. Geosyst.* **9**, Q04021 (2008).

33. Perrier, F. et al. A direct evidence for high carbon dioxide and radon-222 discharge in Central Nepal. *Earth Planet. Sci. Lett.* **278**, 198–207 (2009).

34. Girault, F. et al. The Syabru-Bensi hydrothermal system in central Nepal: 1. Characterization of carbon dioxide and radon fluxes. *J. Geophys. Res. Solid Earth* **119**, 4017–4055 (2014).

35. Richon, P. et al. Temporal signatures of advective versus diffusive radon transport at a geothermal zone in Central Nepal. *J. Environ. Radioact.* **102**, 88–102 (2011).

36. Girault, F., Koirala, B. P., Perrier, F., Richon, P. & Rajaure, S. Persistence of radon-222 flux during monsoon at a geothermal zone in Nepal. *J. Environ. Radioact.* **100**, 955–964 (2009).

37. Girault, F. & Perrier, F. The Syabru-Bensi hydrothermal system in central Nepal: 2. Modeling and significance of the radon signature. *J. Geophys. Res. Solid Earth* **119**, 4056–4089 (2014).

38. Avouac, J.-P., Meng, L., Wei, S., Wang, T. & Ampuero, J.-P. Lower edge of locked Main Himalayan Thrust unzipped by the 2015 Gorkha earthquake. *Nat. Geosci.* **8**, 708–711 (2015).

39. Sapkota, S. N., Bollinger, L. & Perrier, F. Fatality rates of the M<sub>w</sub>~8.2, 1934, Bihar–Nepal earthquake and comparison with the April 2015 Gorkha earthquake. *Earth Planets Space* **68**, 40 (2016).

40. Martin, S. S., Hough, S. E. & Hung, C. Ground motions from the 2015 M<sub>w</sub> 7.8 Gorkha, Nepal, earthquake constrained by a detailed assessment of macroseismic data. *Seismol. Res. Lett.* **86**, 1524–1532 (2015).

41. Grandin, R. et al. Rupture process of the M<sub>w</sub> = 7.9 2015 Gorkha earthquake (Nepal): Insights into Himalayan megathrust segmentation. *Geophys. Res. Lett.* **42**, 8373–8382 (2016).

42. Stevens, V. L. & Avouac, J.-P. Millenary M<sub>w</sub> > 9.0 earthquakes required by geodetic strain in the Himalaya. *Geophys. Res. Lett.* **43**, 1118–1123 (2016).

43. Adhikari, L. B. et al. The aftershock sequence of the April 25th 2015 Gorkha-Nepal earthquake. *Geophys. J. Int.* **203**, 2119–2124 (2015).

44. Galetzka, J. et al. Slip pulse and resonance of Kathmandu basin during the 2015 M<sub>w</sub> 7.8 Gorkha earthquake, Nepal imaged with geodesy. *Science* **349**, 1091–1095 (2015).

45. Gualandi, A. et al. Pre- and post-seismic deformation related to the 2015, M<sub>w</sub>7.8 Gorkha earthquake, Nepal. *Tectonophysics* **714–715**, 90–106 (2017).

46. Zheng, G., Xu, S., Liang, S., Shi, P. & Zhao, J. Gas emission from the Qingzhu River after the 2008 Wenchuan Earthquake, Southwest China. *Chem. Geol.* **339**, 187–193 (2013).

47. Groppo, C., Rolfo, F., Castelli, D. & Mosca, P. Metamorphic CO₂ production in collisional orogens: Petrological constraints from phase diagram modeling of Himalayan, scapolite-bearing, calc-silicate rocks in the NKC(F)MAS(T)-HC system. *J. Petrol.* **58**, 53–83 (2017).

48. Rapa, G. et al. Titanite-bearing calc-silicate rocks constrain timing, duration and magnitude of metamorphic CO₂ degassing in the Himalayan belt. *Lithos* **292–293**, 364–378 (2017).

49. Scholz, C. H., Sykes, L. R. & Aggarwal, Y. P. Earthquake prediction – Physical basis. *Science* **181**, 803–810 (1973).

50. Chiodini, G., Cioni, R., Guidi, M., Raco, B. & Marini, L. Soil CO₂ flux measurements in volcanic and geothermal areas. *Appl. Geochem.* **13**, 543–552 (1998).

51. Girault, F., Koirala, B. P., Bhattarai, M. & Perrier, F. Radon and carbon dioxide around remote Himalayan thermal springs. *Geol. Soc. Lond. S. P.* **451**, 155–181 (2016).

52. Becker, J. A. Quantification of Himalayan Metamorphic CO₂ Fluxes: Impact on Global Carbon Budgets. PhD Thesis, University of Cambridge (2005).

53. Assayag, N., Rivé, K., Ader, M., Jézéquel, D. & Agrinier, P. Improved method for isotopic and quantitative analysis of dissolved inorganic carbon in natural water samples. *Rapid Commun. Mass Spectrom.* **20**, 2243–2251 (2006).

54. Girault, F., Perrier, F. & Przylibski, T. A. Radon-222 and radium-226 occurrence in water: A review. *Geol. Soc. Lond. Spec. Publ.* **451**, 131–154 (2016).

55. Perrier, F., Aupiais, J., Girault, F., Przylibski, T. A. & Bouquerel, H. Optimized measurement of radium-226 concentration in liquid samples with radon-222 emanation. *J. Environ. Radioact.* **157**, 52–59 (2016).

56. Wang, C.-Y. & Manga, M. Hydrologic responses to earthquakes and a general metric. *Geofluids* **10**, 206–216 (2010).

57. Elkhoury, J. E., Brodsky, E. E. & Agnew, D. C. Seismic waves increase permeability. *Nature* **441**, 1135–1138 (2006).

58. Akkar, S. & Bommer, J. J. Empirical equations for the prediction of PGA, PGV, and spectral accelerations in Europe, the Mediterranean Region, and the Middle East. *Seismol. Res. Lett.* **81**, 195–206 (2010).

59. Mosley, R. B., Menetrez, M. Y., Snoddy, R. & Brubaker Jr., S. A. The influences of diffusion and advective flow on the distribution of radon activity within USEPA's soil chamber. *Environ. Int.* **22**, S521–S533 (1996).

60. Vasseur, G., Michard, G. & Fouillac, C. Contraintes sur la structure profonde et le fonctionnement du système hydrothermal de Chaudes-Aigues (France). *Bull. Hydrogéol.* **4**, 3–17 (1997).

61. Derry, L. A., Evans, M. J., Darling, R. & France-Lanord, C. Hydrothermal heat flow near the Main Central Thrust, central Nepal Himalaya. *Earth Planet. Sci. Lett.* **286**, 101–109 (2009).

62. Byrdina, S. et al. Dipolar self-potential anomaly associated with carbon dioxide and radon flux at Syabru-Bensi hot springs in central Nepal. *J. Geophys. Res.* **114**, B10101 (2009).

63. Le Fort, P. Himalayas: Collided range, present knowledge of continental arc. *Am. J. Sci.* **A275**, 1–44 (1975).

64. Kotarba, M. Hydrogeological investigations in Seti Khola and Trisuli thermal springs areas (Nepal Himalayas). *Geologia* **12**, 37–51 (1986).

## Acknowledgements

We warmly thank Soma Nath Sapkota and the Department of Mines and Geology for enthusiastic support over the years. We are indebted to Buddha Lama, Niraj Jairu and Ohm Tharum for assistance in the field, and to Nabin Tamang and Dawa Tamang for giving access to their property in Syabru-Bensi. Svetlana Byrdina is thanked for preparing the CO₂ flux data obtained in 2006–2007. We are thankful to Jean-Philippe Avouac, Steve Ingebritsen, Yves Gaudemer, Damien Jougnot, Christian Fouillac, Pascal Bernard, Hélène Lyon-Caen and Marc Chaussidon for fruitful discussions. Thomas Rigaudier is thanked for the measurement of carbon isotopic ratios in CRPG. Carine Chaduteau is thanked for the measurement of carbon isotopic ratios in IPGP. Narayan P. Malla is thanked for giving access to the Rasuwagadhi Hydroelectric Project. Hélène Bouquerel, François Prevot, Marie Bouih & Roxane Ferry are thanked for their participation in the fieldwork. Financial support for this work was essentially provided by public and private funding of the IPGP team "Physics of Natural Sites", partially complemented in 2015 by a Tellus-SYSTER INSU CNRS funding to F.G. and in 2016 by a CNRS attribution to C.F.-L. This is IPGP contribution number 3942.

## Author contributions

L.B.A. provided aftershock catalogue and processed Fig. 1. C.F.-L. performed carbon isotopic measurements of the gas phase in the laboratory. P.A. performed carbon isotopic and dissolved inorganic concentration measurements of the water phase in the laboratory. B.P.K. and M.B. participated to field CO₂ flux measurements. S.S.M. participated in the discussions and managed the access to the Sanjen Hydroelectric Project sites. C.G. and F.R. performed water measurements a few weeks before the Gorkha earthquake at four hot springs. L.B. managed the assistance to NSC, seismic data processing, and provided detailed DEM and shared experience. F.G. performed water measurements, CO₂ and radon flux measurements, and collected all samples. F.P. performed radon flux measurements and processed Supplementary Fig. 5. F.G. and F.P. wrote the paper.

## Additional information

Supplementary Information accompanies this paper at https://doi.org/10.1038/s41467-018-05138-z.

Competing interests: The authors declare no competing interests.

Reprints and permission information is available online at http://npg.nature.com/reprintsandpermissions/

Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_8.png" /></figure>

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2018