Persistent CO2 emissions and hydrothermal unrest following the 2015 earthquake in Nepal

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 (CO2) emissions associated with seismic activity, so far only observed in magmatic context. Here we report spectacular non-volcanic CO2 emissions and hydrothermal disturbances at the front of the Nepal Himalayas following the deadly 25 April 2015 Gorkha earthquake (moment magnitude Mw = 7.8). The data show unambiguously the appearance, after the earthquake, sometimes with a delay of several months, of CO2 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 CO2 from active orogens. Possible pre-seismic effects need further confirmation.

U nderstanding fluid-earthquake interplay has long received a lot of attention 1 . For the last 20 years, there has been a growing number of evidence for fluid-driven earthquakes 2 . Deep fluids have been shown controlling aftershock distribution in the tectonic contexts of rifting [3][4][5] , subduction 6,7 , and reverse 8 and strike-slip faulting 9 . Contemporaneously, numerous observations of earthquake-induced effects on nearsurface aquifers have been reported 10 , including mainly changes in stream and spring discharge [11][12][13] , groundwater level 14,15 and temperature [16][17][18][19] in various tectonic contexts. These observations have suggested that earthquakes dynamically affect permeability of the Earth's crust 20,21 . Carbon dioxide (CO 2 ) emissions were observed in association with seismicity in the case of the Matsushiro swarm 9 in Japan or, recently, at the Lassen volcano 22 (Cascades Range, USA) and in the Eger Rift 23 (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 24 , the Main Himalayan Thrust (MHT) accommodating at 2 cm year −1 half of the shortening between the two continents 25 . The largest earthquake in Nepal before 2015, the 1934 Bihar-Nepal earthquake (moment magnitude M w~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 26 . Inter-seismic deformation is associated with intense background seismicity 27 , with 4-5 events of local magnitude M L > 5 per year, concentrated between 10 and 25 kilometres depth at the foot of the Himalayan topographic rise 28 . This regionthe Main Central Thrust (MCT) zone 29 also exhibits numerous hydrothermal systems 30 (Fig. 1). Following evidence of degassing from chemical and isotopic analysis of hot springs and rivers 31,32 , large CO 2 emissions were discovered near hot springs 30,33 , with CO 2 fluxes at places similar to diffusive fluxes from active volcanoes 30,34 . The seasonal and yearly stability of the soil-gas radon concentration time-series 34,35 , the invariant radon-CO 2 fluxes relationship 30,36 , and the results of watering experiments 36 at selected sites attest to the remarkable temporal stability of these hydrothermal systems, even during monsoon. These non-volcanic CO 2 emissions are characterised by radiogenic helium, high radon content, and carbon isotopic compositions suggesting metamorphic CO 2 production at >5 kilometres depth [30][31][32]34 . In the accepted conceptual model 31,32 , from the decarbonation source at >5 kilometres depth, CO 2 percolates through fracture networks in the MCT zone, where it mixes with meteoric water. Near the water table, a fraction of CO 2 may degas, and water discharges eventually at the surface as a hot spring. Degassed CO 2 may also be transported faster to the surface through a network of interrelated faults without interaction with hydrothermal circulations 37 .
The M w = 7.8 25 April 2015 Gorkha earthquake 38 (Fig. 1) caused >9000 deaths, with fatality rates >1% in mountainous areas north of Kathmandu 39 , reaching 100% at some places 40 . It partly ruptured the MHT along a 120-kilometre-long segment east from the epicentre 41  MCT zone (Fig. 1), releasing only partially the stored elastic energy sufficient to produce a possible M w~9 mega-earthquake 42 . The mainshock was followed by an intense aftershock sequence 43  In this paper, we report spectacular outbursts of CO 2 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 2 from active orogens.

Results
Outbursts of CO 2 triggered by the Gorkha earthquake. New CO 2 emissions were observed after the Gorkha earthquake, sometimes spatially associated with substantial changes from the pre-existing CO 2 emissions (Table 1 and Supplementary Table 1). Because a large number of CO 2 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 2 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 2 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 2 discharge in the area now reaches 1010 ± 110 mmol s −1 (3.8 ± 0.4 ton d −1 ), 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 2 fluxes ( >10 5 g m −2 d −1 ) showed one-order-of-magnitude smaller CO 2 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 2 emissions and changes in pre-existing CO 2 emissions were also observed at other locations in the same valley ( Supplementary Fig. 5). In Timure, nine kilometres to the north, CO 2 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 2 emission increase of a factor of 8 ± 3 compared with preearthquake values. By contrast, at the Chilime hot spring site, seven kilometres to the west of Timure (Fig. 1), while CO 2 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 2 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 2 emissions appeared at several locations in a tunnel being excavated nearby (Supplementary Fig. 7 and Movie 1). The CO 2 emission in the tunnel was still present in January 2018. The CO 2 concentration in the air of the tunnel then ranged from 4 to Table 1 Post-seismic effects on the carbon dioxide emissions at five hydrothermal sites in Central Nepal. Data of CO 2 flux, total CO 2 discharge, and carbon isotopic ratio of CO 2 emissions are compiled before and after the Gorkha earthquake Site Location Before/after Gorkha earthquake 5 vol%, thus creating a major health hazard, and the total CO 2 discharge was estimated in the tunnel to 580 ± 150 mmol s −1 (2.2 ± 0.6 ton d −1 ), 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 2 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 2 emission was still present in January 2018, with a total discharge higher than 1240 ± 330 mmol s −1 (4.7 ± 1.3 ton d −1 ), similar to the whole CO 2 discharge observed in Syabru-Bensi. In addition to the flux from the bank, an innovative method was used to measure the CO 2 flux from and through the hot water pond (see Methods and Supplementary Fig. 8). This is the most spectacular new CO 2 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 postseismic CO 2 outburst is unusual. For instance, small outbursts following the 2008 Wenchuan earthquake in China lasted only a few months 46 . In the Marsyandi valley, further west ( Fig. 1), by contrast, no comparable phenomenon was observed and peak CO 2 fluxes (5300 to 28,700 g m −2 d −1 ), while significant, are smaller than in the Upper Trisuli valley.
All these CO 2 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 2 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 2 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.  62 . Numbers are carbon isotopic ratios of the gaseous CO 2 (δ 13 C, relative to V-PDB). Total CO 2 discharge amounts to 480 ± 50 mmol s −1 (1.8 ± 0.2 ton d −1 ) before the earthquake and 1010 ± 110 mmol s −1 (3.8 ± 0.4 ton d −1 ) after, giving a post-seismic overall increase of a factor of 2.1 ± 0.3. The surface area of high CO 2 emission ( > 500 g m −2 d −1 ) increased by a factor of about 1.9, with the barycentre shifting eastward by 15 ± 3 metres. Example of CO 2 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 2 -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. Carbon isotopic ratio of the CO 2 emissions (δ 13 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 δ 13 C values range from -6.9 ‰ to -0.1 ‰, and appear similar at the various sites of a given valley. CO 2 emissions in Sanjen and Syabru-Bensi, although separated by eight kilometres, have comparable CO 2 concentration (96-98%) and δ 13 C signature (from -0.9 ‰ to -0.7 ‰), precluding shallow sources, but instead suggesting similar CO 2 source and transport from a crustal-scale reservoir. In addition, the δ 13 C values remain relatively stable after the earthquake ( Supplementary Fig. 4), indicating that the earthquake revealed a pre-existing CO 2 reservoir.
Hydrothermal unrest triggered by the Gorkha earthquake. The CO 2 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 preearthquake 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 preseismic 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  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 SBP0 are compiled from our database and the literature 30 . Timing of the Gorkha earthquake, the Chilime hot spring cessation and the Sanjen CO 2 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 2 flux values obtained in 2006 (ref. 33  sometimes pre-seismic changes in hot spring temperature have been reported in the literature 16,19 , 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 −1 ; 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 2 outburst in Sanjen.
Dissolved Inorganic Carbon (DIC) concentration and isotopic ratio (C DIC and δ 13 C DIC ; 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 DIC 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 δ 13 C DIC 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 2 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 preearthquake conditions, while those more dependent on superficial effects (e.g., SBB5) did not yet. First non-volcanic earthquake-induced gaseous changes. Hydrologic responses to seismic stimulation have been observed in numerous instances 10 , 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 w = 3.85 at 5.7 kilometres distance, corresponding to SED∼0.2 J m −3 and PGV∼0.2 cm s −1 , caused an outburst of geothermal fluids, explained by a two-fold permeability increase 22 . In the case of the Matsushiro earthquake swarm in Japan, modelling indicates the necessity of a 2-order-ofmagnitude permeability increase with a small overpressure of only a few megapascals 9 . In Syabru-Bensi, the Gorkha earthquake produced SED∼63 J m −3 and PGV∼26 cm s −1 (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 station 44,45 (Fig. 1), giving PGV∼49 cm s −1 (Supplementary Fig. 13). While the first clearly documented instances in the Himalayas, our observed near-field CO 2 outbursts and hydrothermal unrest remain compatible with previously compiled earthquake-induced changes ( Supplementary Fig. 14). Our observed changes in the CO 2 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 2 is degassed from hydrothermal waters, needs some modifications at least in part. Indeed, previous work in Syabru-Bensi 34 identified large CO 2 emissions without the presence of an important hot spring in the vicinity. Now, the evidence for an extended reservoir of crustal CO 2 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 2 source can therefore be the hydrothermal circulations, with near surface degassing and subsurface accumulation in the whole fault zone. Possible additional CO 2 can also be released in the footwall directly from the production source. The earthquake caused (Fig. 4b) an increasing water discharge, with increased CO 2 degassing, or better communication to the surface of previously accumulated CO 2 . In the case of the Chilime valley (Chilime, Sanjen and Brapche) (Fig. 4c, d), however, the pervading presence of CO 2 , revealed by the earthquake, kilometres away from hot springs, is better explained if CO 2 is directly accumulated from below, possibly through the water table, without advection by hydrothermal circulations. Strikingly, CO 2 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 2 in the Himalayan crust, suggesting that metamorphic CO 2 produced at depth is huge, as independently shown by petrological estimates 47,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 2 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 2 emissions appear consistent with hydrothermal outburst effects, as observed in Lassen 22 , suggesting that apparent two-fold permeability increases and/or changes in hydrothermal pathways for pre-existing CO 2 emissions could have been major effects of the Gorkha  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 2 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 2 can also move upwards independently, reaching the surface without hot spring (Sanjen, Syabru-Bensi GZ3) and accumulate in the subsurface reservoir revealed by CO 2 emissions following the earthquake. The main effect of the earthquake is an increase of vertical permeability with transient or permanent changes NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-05138-z ARTICLE NATURE COMMUNICATIONS | (2018) 9:2956 | DOI: 10.1038/s41467-018-05138-z | www.nature.com/naturecommunications earthquake. First-order modelling indeed confirms that at least a doubling of permeability could explain CO 2 flux increases (Supplementary Fig. 15). In addition, first-order modelling supports that a lasting permeability increase 18 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 45 . The global pre-seismic water-cooling could be evidence for a preseismic dilatancy effect 49 of the impending rupture zone.
Compared with reported earthquake-induced hydrological effects 10,16,19,46 , some of our observations are unusual: first, the six-month delay of the Sanjen CO 2 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 2 emissions and hydrothermal unrest. In Lassen 22 , the permeability increase was inferred to be maintained for 50-60 days. In the Napa valley 13 , a vertical permeability change, following diffusion of a pressure pulse, was proposed to explain postseismic stream discharge, suggesting basin-scale fluid diffusivity of 0.1 m 2 s −1 . 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 2 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 −1 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 2 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 Octoberbeginning of November), less than 15 kilometres from Syabru-Bensi, producing PGV > 1 cm s −1 (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 2 may in turn play a role in the installation of the next inter-seismic regime. Alternatively, the still evolving CO 2 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 42 , long-term monitoring of CO 2 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 2 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 2 emissions observed in the Narayani basin suggest non-stationary metamorphic CO 2 production, and that its current estimate 30,32 (>1.3 × 10 10 mol year −1 with >(1.0 ± 0.2) × 10 8 mol year −1 from direct gaseous CO 2 emissions), independently confirmed by petrological studies 48 , may be enhanced during a significant fraction of the seismic cycle. Metamorphic CO 2 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 50 was used to measure surface CO 2 flux and to quantify the total CO 2 discharge of a given site and the associated uncertainties, whose assessment is based on numerous systematic tests and our 10-year experience 34 . The method is robust, even in remote locations 51 and during the monsoon 36 , and allows the measurement of CO 2 fluxes over more than five orders of magnitude ( Supplementary Fig. 2). The increase in CO 2 concentration in the chamber is measured using various portable infrared CO 2 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 2 flux is expressed in grams per squared metre per day (g m −2 d −1 ). The total CO 2 discharge, expressed in mol s −1 (or ton d −1 ), is estimated using the CO 2 flux data-set by kriging and interpolation procedures 34 . CO 2 fluxes (n = 1720) and total CO 2 discharges obtained before the Gorkha earthquake (from 2006 to 2011) are published elsewhere 30,33,34,51 . Here we present for the first time CO 2 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 2 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 2 flux from or through water was measured at the new degassing site of Machhakhola (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 2 discharge observed in Machhakhola (Supplementary Movie 2).
Carbon isotopic composition of the gas phase. We sampled gas in the field using evacuated glass tubes. CO 2 fraction of the gas sample was determined manometrically. The δ 13 C of CO 2 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) 32,34 . External repeatability of a given sample was ±0.1 ‰. The twenty δ 13 C values measured before the Gorkha earthquake are published elsewhere 30,33,34,52 . 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 −3°C ) and high-sensitivity (10 −4°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; HI98107 and HI98130 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 screw cap vials of 12 millilitres volume each. DIC concentration (C DIC = [H 2 CO 3 ] + [HCO 3 -] + [CO 3 2-]) and its isotopic composition (δ 13 C DIC ), expressed in mmol L −1 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 53 . The relative experimental uncertainty of C DIC was 1-2% and the experimental uncertainty of a given δ 13 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 δ 13 C DIC values measured before the Gorkha earthquake are published elsewhere [30][31][32][33] . 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 54 . After air-water equilibrium is reached by manual shaking, radon concentration is inferred from scintillation flask sampling and photomultiplier counting, as described elsewhere 51 . Radon concentration in water is expressed in Bq L −1 .
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 55 . Expressed in mBq L −1 , 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 2 degassing areas. The CO 2 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 30,51 . In remote places, the detection of previously unknown thermal springs relied on the use of hand-held thermal infrared cameras 34,35,51 (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 2 emissions and hot springs, we calculated SED 56 and vertical PGV 57 produced by the Gorkha earthquake and its main aftershocks. Empirical equations are used to estimate SED 56 and PGV 58 . Peak Dynamic Stress (PDS) is estimated using: PGV × shear modulus (3 × 10 10 Pa) / shear wave velocity (3500 km s −1 ). 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 2 flux before and after the earthquake, we used a 2-D model, described elsewhere 37 , in which a vertical fault (f) surrounded by two media (a and b) transports the gaseous CO 2 to the surface by advection. Each medium is characterised by connected porosity (ε) and permeability (k). A pressure source is fixed at depth. Darcy's law is defined according to a pressure distribution (p) that follows Δp 2 = 0 (ref. 59 ). 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 37 , 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 37 . 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 60 (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 ρ and specific heat C P . Temperature at depth h is defined by T(h) = βh + T(0), where β is the thermal gradient, which is poorly constrained in this hydrothermal region 61 and taken equal to 55-75°C km −1 . We consider a quasi-static state having characteristic length λ T and we define θ = T-T(0) and α = ρC P Qλ T /(pK T h). Using the initial conditions, the differential equation αh(dθ/dz) + θ = -βz has the following solution: θ = αhβ(1-e -1/α ). 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).