Access the full text.
Sign up today, get DeepDyve free for 14 days.
J. Vandemeulebrouck, D. Stemmelen, T. Hurst, J. Grangeon (2005)
Analogue modeling of instabilities in crater lake hydrothermal systemsJournal of Geophysical Research, 110
A. Hurst, R. Dibble (1981)
Bathymetry, heat output and convection in Ruapehu crater lake, New ZealandJournal of Volcanology and Geothermal Research, 9
E. Sartori (2000)
A critical review on equations employed for the calculation of the evaporation rate from free water surfacesSolar Energy, 68
S. Fitzgerald, A. Woods (2003)
THE STABILITY OF TWO-PHASE GEOTHERMAL RESERVOIRS
Rudolph Merwe, E. Wan (2004)
Sigma-point kalman filters for probabilistic inference in dynamic state-space models
B. Scott (1994)
Cyclic activity in the crater lakes of Waimangu hydrothermal system, New ZealandGeothermics, 23
D. Huang, T. Schneider, A. Stuart (2021)
Iterated Kalman methodology for inverse problemsJ. Comput. Phys., 463
T. Hurst, B. Christenson, J. Cole-Baker (2012)
Use of a weather buoy to derive improved heat and mass balance parameters for Ruapehu Crater LakeJournal of Volcanology and Geothermal Research, 235
G. Brown, H. Rymer, David Stevenson (1991)
Volcano monitoring by microgravity and energy budget analysisJournal of the Geological Society, 148
(1966)
New Zealand Volcanology, Central Volcanic Region, 1st edn. Department of Scientific and Industrial Research
C. Torrence, G. Compo (1998)
A Practical Guide to Wavelet Analysis.Bulletin of the American Meteorological Society, 79
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations
S. Särkkä (2013)
Bayesian Filtering and Smoothing, 3
S Särkkä (2013)
Bayesian Filtering and SmoothingCambridge University Press, Cambridge.
D. Stevenson (1992)
Heat transfer in active volcanoes: models of crater lake systems
J. White (2018)
A model-independent iterative ensemble smoother for efficient history-matching and uncertainty quantification in very high dimensionsEnviron. Model. Softw., 109
A. Hurst, H. Bibby, B. Scott, M. Mcguinness (1991)
The heat source of Ruapehu crater lake; deductions from the energy and mass balancesJournal of Volcanology and Geothermal Research, 46
B. Christenson, A. Reyes, R. Young, A. Moebis, S. Sherburn, J. Cole-Baker, K. Britten (2010)
Cyclic processes and factors leading to phreatic eruption events: Insights from the 25 September 2007 eruption through Ruapehu Crater Lake, New ZealandJournal of Volcanology and Geothermal Research, 191
N. Fournier, F. Witham, M. Moreau-Fournier, L. Bardou (2009)
Boiling Lake of Dominica, West Indies: High-temperature volcanic crater lake dynamicsJournal of Geophysical Research, 114
G. Brown, H. Rymer, J. Dowden, P. Kapadia, D. Stevenson, J. Barquero, L. Morales (1989)
Energy budget analysis for Poás crater lake: implications for predicting volcanic activityNature, 339
Friedlander B (1898) Some notes on the Volcanoes of the Taupo District
H. Rauch, Fang-Wu Tung, C. Striebel (1965)
On the maximum likelihood estimates for linear dynamic systems
RR Dibble (1966)
93
E. Adams, Douglas Cosler, K. Helfrich (1990)
Evaporation from heated water bodies: Predicting combined forced plus free convectionWater Resources Research, 26
(1978)
Thermodynamic and Transport Properties of Fluids : SI Units, 2nd edn. Basil Blackwell
T. Hurst, T. Hashimoto, A. Terada (2015)
Crater Lake Energy and Mass Balance
(2002)
A New Approach to Linear Filtering and Prediction Problems
H. Togt (2003)
Publisher's NoteJ. Netw. Comput. Appl., 26
S. McNutt (2005)
VOLCANIC SEISMOLOGY
Volcanic lakes often capture a significant amount of volcanic heat emission and thus provide a unique opportunity to monitor changes inside the volcano. We present a Bayesian inversion method to automatically infer changes in volcanic heat emission over time at the base of a volcanic lake from lake monitoring data using a non-linear Kalman Smoother. Our method accounts for the, sometimes large, uncertainties in observations and the underlying physics- based model to generate probabilistic estimates of heat emission. We verify our results using a synthetic test case and then estimate the daily heat input rate into Mt. Ruapehu’s Crater Lake, New Zealand, between 2016 and 2022. Time- frequency analysis of the heat input rate shows dominant periods of heating cycles ranged between 100 - 250 days. The period between 2017 and 2020 was dominated by shorter cycles and greater-than-average heat input rate which points to changes in the magmatic heat supply and the hydrothermal system during this time. Keywords Crater lakes, Non-linear Kalman Smoother, Volcanic heat emission Monitoring changes in temperature, water mass, and Introduction ion concentration in these lakes can, therefore, serve The amount and rate of heat emitted by a volcano over as proxies for changes in gas- and steam-release by the time provides important clues about volcanic processes hydrothermal system and, ultimately, the melt zone at depth. For example, a sudden rise in the heat output beneath the volcano. rate combined with increased microgravity signals may Ruapehu Crater Lake (Te Wai ā-moe in Te Reo Maori result from the emplacement of fresh magma at shallow and referred to as RCL from hereon) on Mt. Ruapehu, depth (Brown et al. 1991). A slow decrease in the heat an active andesitic stratovolcano in the center of New emission rate together with a high carbon-sulfur ratio of Zealand’s North Island (Fig. 1), is one such example with the emitted gas may point to the formation of a hydro- records of observations dating back to the 19th century thermal seal (Christenson et al. 2010). (Friedlander 1898). Mt. Ruapehu showed eruptive activ- Volcanic lakes help with estimating heat emissions. ity on at least 603 days since 1830 with many phreatic Often, the majority of a volcano’s heat and gas output and phreato-magmatic eruptions and two major mag- passes through the crater lake. Crater lakes act thus like a matic episodes (Scott , 2013; Historic Eruptive Activity filter, integrating heat and gas input from the vents enter - Mt. Ruapehu , 2022). Although the latest magmatic event ing the lake. between June 1995 to November 1997 and a subsequent dam collapse changed the shape and volume of the lake, *Correspondence: it has been largely unchanged since then and the lake’s Yannik Behr y.behr@gns.cri.nz bathymetry is well known from times when it was empty. GNS Science, 114 Karetoto Rd, 3384 Taupo, New Zealand © The Author(s) 2023. 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 licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence 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 licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/. The Creative Commons Public Domain Dedication waiver (http:// creat iveco mmons. org/ publi cdoma in/ zero/1. 0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data. Behr et al. Journal of Applied Volcanology (2023) 12:1 Page 2 of 10 Fig. 1 Image of Mt. Ruapehu’s crater lake from 2008 (©Visual Media Library, GNS Science, ID 6277, photographer Karen Britten). The red circle marks the location of the lake outlet. The location of Mt. Ruapehu is shown as a red triangle on the inset map First attempts to quantify the amount of heat entering One way to handle these difficulties is to make edu - the lake date back to 1966 (Dibble, 1966). Since 2016 cated guesses about unconstrained parameters, inter- water temperature and lake level have been monitored polate between irregularly sampled data, and ignore continuously and water samples for ion concentration uncertainties. We will call this the deterministic measurements have been taken on an approximately method. While it can lead to reasonable results, they monthly basis. This enables us to continuously estimate have the following problems: heat input into the lake and, in principle, automatically detect changes as soon as they occur. 1 It is not possible to quantify our confidence in the We typically infer the heat input from a mass and accuracy of the results. energy balance calculation: the total mass and energy 2 Comparing results to other models becomes very dif- flow into and out of the lake has to account for the ficult. observed changes in lake temperature, water level and ion concentration (Hurst and Dibble 1981; Hurst et al. To illustrate the second point let us imagine that we 1991; Stevenson 1992; Fournier et al. 2009; Scott 1994). make an estimate of intruded magma volume based on At RCL, and probably at most other volcanic lakes, this measured gas flux. Estimating the heat flux based on is an underdetermined problem, that is the number of this magma volume using, for example, the methods unknown parameters in the mass and energy balance described in Stevenson (1992), will most likely not match model exceeds the number of independent observa- the heat flux estimated with the deterministic method. tions. Hence, it is necessary to constrain some model The reason for that is that gas flux, the assumptions on parameters using prior knowledge. Further, some amount of exsolved gas from a magma body, and the observations and input parameters come with large mass and energy balance model all carry large uncertain- uncertainties and are sampled irregularly. The mass and ties. Instead of trying to find a single estimate of magma energy balance model itself also comes with consider- volume that fits all observations perfectly we better try to able uncertainty as it is only a rough approximation of find a range of magma volumes that are compatible with the physical processes within the crater lake. our models and observations within their uncertainties. B ehr et al. Journal of Applied Volcanology (2023) 12:1 Page 3 of 10 The result would be a probabilistic model that can pro - the reader to previous publications for more detailed dis- vide lower and upper bounds on magma volume. cussions (Stevenson 1992; Hurst et al. 2015). There is no In the following, we propose a probabilistic (Bayesian) consensus on the best quantitative model of evaporation inference approach of the time-varying heat input rate from warm lakes and, in fact, the most suitable equation which naturally handles these difficulties and transpar - may depend on the location of and available observations ently propagates uncertainties of observations and prior from the lake. Appendix A describes the equations used assumptions into posterior uncertainties. in this study in detail which include the recent develop- ment on forced (i.e. wind-induced) convection (Sartori 2000). It is important to note here that Q is a function of Probabilistic inference method e T and the wind speed above the lake. Physical model dT Following from Eq. 3 the temperature change can The mass and energy balance model can be written in dt then be written as: terms of three coupled ordinary differential equations (Hurst et al. 1991; Stevenson 1992). The rate of lake water dT 1 T ̇ ̇ ̇ ̇ ̇ ̇ = Q + Q − c TM − Q − c TM − M i r w o e w s (4) mass change ( M ) is the difference between the rate of dt c M M steam input through the vent system beneath the lake Finally, the change in ion concentration expressed in total ( M ), the water inflow rate through melt and precipita - 2+ ˙ ˙ ion amount (X) of, for example, Mg can be defined as: tion ( M ), the outflow rate through overflow ( M ), and s o the rate of mass loss through evaporation ( M ): dX X =−M (5) dM dt M ˙ ˙ ˙ ˙ = M + M − M − M (1) i s o e dt This equation is only true under the assumption that X only changes due to dilution by the inflow of meltwater The energy content of the lake can be described with the and steam and there is no re-supply between eruptive notion of enthalpy, H, which is defined relative to a ref - episodes. Equation 5 gives an estimate of the average out- erence enthalpy. In the following we take the enthalpy flow over periods on the order of months but does not of the lake at temperature T = 0 C = T as reference capture shorter term variations in the outflow. In the H = H (T ) . At constant pressure, the enthalpy of the 0 0 following we will also only account for water outflow lake can then be written as a function of T, T , M and the at the outlet and ignore seepage through the lake floor. specific heat of water, c : This approximation is supported by the observation that amongst the various rivers at Mt Ruapehu only the one H (T ) = H + c Mτdτ (2) 0 w fed directly by RCL’s outlet bears RCL’s chemical signa- ture (B. Scott, 2022, GNS Science, pers. commun.). For all following equations we take c = constant , such The main challenge in solving the mass and energy bal - that the enthalpy of lake water at temperature T can be ance equations for Q arises from the non-linear relation- written as H (T ) = H + c MT . 0 w ship between Q and T as shown in Fig. 2. The rate of enthalpy change of the lake water can then be defined as a function of the rate of energy input from Inference volcanic sources (Q ), the rate of energy gain due to solar The objective of the inference, or inverse, problem is irradiation from short-wavelength radiation (Q ), the rate to estimate the rate at which heat is entering the lake of energy loss through outflow/seepage ( c TM ), the w o ˙ through the hydrothermal system below (Q ) from the rate of surface energy losses due to evaporation (forced given observations and the mass and energy balance and free convection), sensible heat, and long-wavelength model described in Section 2.1. Let y be the time series 1:t radiation (Q ), and the rate of energy loss from heating e ˙ of t observations of T, M, X, and M and x be the n incoming surface water ( c TM ): w s dH d(c MT ) ˙ ˙ ˙ ˙ ˙ ˙ ˙ (3) = = c TM + c MT = Q + Q − c TM − Q − c TM w w i r w o e w s dt dt Because RCL is at an altitude of ∼ 2500 m and big parts model parameters described in Section 2.1 at time step of its catchment are covered by snow and ice all year k. Note that x also includes T, M, X and M . We dis- k o round we assume the temperature of incoming surface tinguish here between the state of the lake and observa- water to be at 0 C . Several empirical and theoretical tions of this state (T, M, X, M ) which are to some degree equations exist to describe the surface losses and we refer uncertain. Solving the inverse problem probabilistically Behr et al. Journal of Applied Volcanology (2023) 12:1 Page 4 of 10 Fig. 2 Surface energy loss (Q ) with respect to lake water temperature (T). The model parameters that have been held fixed are shown in the title means that we are looking for a probability distribu- p(x x )p(x y ) k+1 k k+1 1∶t p(x y ) = p(x y ) dx (9) k 1∶t k 1∶k k+1 tion that can best explain our observations taking into p(x y ) k+1 1∶k account the mass and energy balance model and all rel- with evant uncertainties. More concisely we are looking for: x ∼ p(x |y ) k k (6) 1:t p(x |y ) = p(x |x )p(x |y )dx (10) k+1 k+1 k k k 1:k 1:k which means that the variable x has the probability dis- tribution p(x |y ) . The latter is the conditional prob - k 1:t p(x |y ) = p(y |x )p(x |y ) (11) ability of x given all observations until time t. k 1:k k k k 1:k−1 Solving Eq. 6 using standard strategies like Markov Chain Monte Carlo is computationally very expensive as the number of variables grows exponentially with the z = p(y |x )p(x |y )dx (12) k k k k k 1:k−1 length of the time series and the number of model param- eters. Alternatives are Bayesian Filtering and Smoothing algorithms of which the Kalman Filter and Smoother are p(x |y ) = p(x |x )p(x |y )dx k k k−1 k−1 k−1 1:k−1 1:k−1 the closed form solution for linear models (Kalman 1960; Rauch et al. 1965; Särkkä 2013). Note that to infer x at (13) time step k we use all available observations, that is t ≤ k where we made use of the Chapman-Kolmogorov and t > k . This is called a Bayesian or Kalman Smoother. equation: For the Kalman Filter only t ≤ k are used. We assume that the states of the lake form a Markov p(x |x ) = p(x |x )p(x |x )dx 3 1 3 2 2 1 2 chain, meaning that the state of the lake at time step k, x , only depends on the state at adjacent time steps: Equations 9 to 13 describe a recursive way to compute p(x |x , y ) = p(x |x ) Eq. 6 in which we only need to know x , p(x |x ) , k 1:k−1 k k−1 (7) 0 k k−1 1:k−1 and p(y |x ) . Equation 9 and 11 are also known as the and Bayesian Smoothing and Bayesian Filtering equation, respectively (see Appendix B for more details on the deri- p(x |x , y ) = p(x |x ) k−1 k:t k−1 k (8) k:t vation). If p(x |x ) and p(y |x ) are linear transforma- k k−1 k tions with normally distributed errors, the Kalman filter With this assumption and using Bayes’ rule, Eq. 6 can be (Kalman 1960) and Kalman smoother (Rauch et al. 1965) rewritten (Särkkä 2013) as: B ehr et al. Journal of Applied Volcanology (2023) 12:1 Page 5 of 10 represent closed form solutions. In our case, these trans- (see Section 3.2). As with real data, average and standard formations are: deviation of all values within a day are calculated. The standard deviation is taken as a proxy for the uncertainty p(x |x ) =f (x ) + q k k−1 k−1 (14) k−1 in observations ( R ). To simulate the disparity between sampling intervals of different data streams we randomly p(y |x ) =x [1 : m]+ r k k k (15) remove most of the observations of total ion amount (X), ˙ ˙ water outflow rate ( M ), and water inflow rate ( M ). This o s where f (x ) is the non-linear physical model described k−1 dataset is then inverted using the Unscented Kalman in Section 2.1 propagated forward in time using a fourth- Smoother to recover the original heat input rate. order Runge-Kutta scheme. If we assume the errors, Figure 3 shows the result of the probabilistic infer- r and q , to be normally distributed ( q ∼ N (0, Q ) , k k k k ence of heat input compared to the true input of the r ∼ N (0, R ) ), we can compute an approximate solu- k k synthetic model. The observations fit well within their tion to Eq. 9 using the Unscented Kalman Smoother (Van error bounds, but we see some discrepancies between der Merwe , 2004; Särkkä , 2013), a non-linear Kalman the original and the recovered heat input rates, especially Smoother. for very steep changes in the input rate. This is the result We estimate R from the daily variations in the obser- of preferring a smooth solution as already mentioned in vations, but we do not know x and Q . S etting x to an 0 0 Section 2.2. arbitrary value within physically reasonable bounds, but It is noteworthy that the mass outflux rates are recov - large uncertainty works well in practice and has very lit- ered quite accurately despite providing very few synthetic tle influence on the results. While Q can be formally observations. optimized we treat it here as a design parameter that we modify such that the observations are matched within their error bounds and the inferred heat input rate var- Real data ies smoothly. The latter is motivated by the fact that the Temperature at RCL is measured by an integrated-circuit physical model described in Section 2.1 is not able to temperature device (Texas Instruments LM35) located model short-term transient changes in lake volume, tem- roughly 1.9 m below the lake’s current overflow level. perature, and ion concentration. These can be caused by The water level sensor is a hydrostatic pressure sensor rain or increased meltwater influx, partial collapse of the (First Sensor KTE/KTU/KTW6000...CS Series) and it lake walls, or the influx of ion-enriched fluids. The under - is co-located with the temperature sensor. Both sensors lying assumption of a perfectly mixed lake typically does have a sampling interval of 15 minutes and transmit data not apply in these situations. several times a day via satellite to GNS Science, which The ability to account for these epistemic uncertainties monitors New Zealand’s geohazards through its GeoNet explicitly is a particular strength of the Bayesian Filtering program (GeoNet Home 2022). Due to the type of sat- and Smoothing equations. It allows us to gain valuable ellite connection only about 30% of the temperature and insights even from very simple models. level measurements arrive at GeoNet’s datacenter; the rest of the data is lost. Water samples are taken manu- Results ally 8-12 times a year from which subsequently, amongst ++ − Synthetic test many other analyses, concentrations of Mg and Cl To evaluate the inference scheme we design a synthetic are determined at GNS Science’s Geothermal Analytical ++ test, consisting of a very large cylinder, the same volume Laboratory. We will focus here on Mg concentration as RCL. We keep the wind speed and surface water inflow as it shows less deviation from the assumption of con- constant at 4.5 m/s and 10 kt/day, respectively. To model stant dilution than Cl . Because of the sparse sampling, ++ the outflow we assume the outlet to be a pipe 0.2 m in new information on Mg concentration is available cross-section and then use Bernoulli’s equation to model approximately once per month. the water outflow based on the water level in the cylinder To estimate uncertainties of temperature and water above the pipe. We assume that the incoming heat (Q ) level observations, we compute their daily averages and ++ follows a combination of a third-order polynomial and a standard deviations. Uncertainties of Mg concentra- sinusoidal function and that the enthalpy of the incoming tions are estimated directly when two samples were taken steam is 2.8 MJ/kg which corresponds to the enthalpy of on the same day. Otherwise, their uncertainty is taken to saturated steam at hydrostatic pressures similar to those be the median of directly estimated uncertainties from at the bottom of the lake (Mayhew and Rogers 1978). To when multiple samples were taken. the resulting synthetic observations we add uncertain- Lake outflow ( M ) is measured once or twice per year, ties similar to those estimated from real observations and we further know from field observations that the Behr et al. Journal of Applied Volcanology (2023) 12:1 Page 6 of 10 Fig. 3 Inversion result for the synthetic test described in Section 3.1. Turquoise lines and markers represent synthetic observations; red dashed lines are the inversion results and black dashed lines are the true input to the synthetic test. Error bars and shaded areas represent ± 3 standard deviations lake is not overflowing beneath a certain water level. We over the more common windowed Fourier transform for fit a sigmoid function to these observations assuming a its ability to detect non-stationary signals (Torrence and maximum outflow rate of 250 l/s. We assume the uncer - Compo 1998). Figure 4c shows the global wavelet spec- tainty in M as the 95 percent confidence interval of the trum which is an average of the wavelet analysis in Fig. 4b curve-fitting. over the whole analysis period. The global wavelet spec - As there is no permanent weather station close to RCL trum can be interpreted similarly to the Power Spectral we assume the wind speed to be 5 m/s. This was the aver - Density. Figure 4c shows that periodicity is dominated age wind speed derived from a temporary deployment of by periods between 100 and 250 days. Which period is a weather buoy on RCL and adjacent permanent weather dominant at any one time changes significantly, and we stations (Hurst et al. 2012). will expand on this further in the following section. Figure 4a shows the heat input rate (Q ) for RCL for lake measurements between 4 March 2016 and 1 Febru- Discussion and conclusion ˙ ˙ ary 2022. During this time Q reached a maximum of 753 Inferring the heat input rate (Q ) into a crater lake is an i i MW and was on average 165 MW. The average standard important part of estimating the total energy budget of deviation was 34 MW and generally increased for lower a volcanic system (Brown et al. 1989). Similar to a volca- values of Q . To put this in context, currently installed no’s gas budget, changes in the heat input rate can either geothermal power plants in New Zealand have a com- indicate changes in the magmatic heat source or in the bined output of ∼ 1000MW. hydrothermal system. To also investigate changes in periodicity we computed Experiments on analog models of hydrothermal sys- the continuous wavelet transform of Q using a complex tems suggest that changes in the heat source and the Morlet wavelet (Fig. 4b). We chose the wavelet transform hydrothermal system should also change the length of B ehr et al. Journal of Applied Volcanology (2023) 12:1 Page 7 of 10 Fig. 4 a) Inference results for the heat input rate Q at RCL between 4 March 2016 and 1 February 2022. The solid line shows the marginal probability and the uncertainty is displayed as shaded region; b) Continuous wavelet transform of the time series shown in a; c) Power spectral density of the heat input rate shown as a solid line with uncertainty displayed as shaded region heating cycles (Vandemeulebrouck et al. 2005; Fitzger- in 2017 a period of higher than mean values starts towards ald and Woods 1997). Figure 4b shows the heating cycles the end of 2017 and ends early 2020. before mid-2017 and after mid-2019 were dominated by To further explore hypotheses like this, the next logical longer periods. From mid-2017 to mid-2019 cycles were step is to include numerical models of magmatic gas and mostly of shorter periods which may indicate a change in heat release as well as hydrothermal heat and fluid flow. the volcanic system during this time. Our inversion method is well suited for such more complex This is further supported by the cumulative difference models as its runtime scales linearly with the number of to the mean for the heat input rate as shown in Fig. 5. This parameters. Several studies have demonstrated the value of type of graph helps to highlight periods when the heat non-linear Kalman Filters for inverse problems with non- input rate was above or below the long-term mean. Before linear, high-dimensional numerical models (e.g. White 2017 and after 2019 the difference appears to fluctuate 2018; Huang et al. 2022). more or less around the mean value. After a dip beginning Fig. 5 Cumulative difference to the mean for the heat input rate into RCL (solid line) and its uncertainty (shaded region) calculated from the time-series in Fig. 4a Behr et al. Journal of Applied Volcanology (2023) 12:1 Page 8 of 10 Previous studies of RCL’s mass and energy balance found where that they required very high values of enthalpy for the incoming steam to avoid violating physical boundary con- T , T are the virtual surface and air temperatures sv av ditions such as negative values for the inflow rate of mete - respectively. This correction accounts for the oric water ( M ) (Hurst et al. 1991, 2015). Their explanation extra buoyancy of water vapor compared to air was that as rising hot fluids enter the lake, relatively cool (see Eq. 18 for details) lake water flows back into the vent. This would lead to an e , e are the saturation vapor pressure and the s a overall lower addition of mass to the lake and resulting in ambient air vapor pressure respectively an enthalpy value that is double the expected for saturated L is the latent heat of evaporation steam at hydrostatic pressures likely to be found at the bot- u is the wind speed tom of the lake. By including uncertainties in our observa- X is the fetch, or characteristic length, of the lake tions, this explanation is still a possibility but no longer a P is the atmospheric pressure requirement to satisfy the physical constraints. Numerical modelling of hydrothermal processes may also shed further light on the heating processes and whether the lake and the vent form a single or two separate convective systems. T T = for x = s or a xv (18) In conclusion, the combination of the mass and energy (1 − 0.378e /P ) x a model, developed in previous studies (Hurst et al. 1991, 2015), with the probabilistic inference method we devel- A.3 Sensible heat oped here provides a powerful method to continuously This describes the heat loss due to conduction and con - infer the heat input rate into RCL. It allows us to combine vection above a hot lake and according to Stevenson disparate data streams, invert them for probabilistic esti- (1992) the ratio between sensible heat, Q , and heat loss mates of the heat input rate and thereby keep a finger on sh due to evaporation, Q can be written as: Mt. Ruapehu’s energy budget. evap 1/3 Q ρc (T − T ) (T − T ) sh a s a s a = R = sh Q L (e − e ) (T − T ) evap s a sv av Appendix A: Surface losses A.1 Long‑wavelength radiation (19) The net loss through long-wavelength radiation can be writ - ten as: where 4 4 Q = Aσ (ǫ T − ǫ T ) (16) rad w a s a ρ i s the air density c is the specific heat of air where 2 We decided to ignore term A as it tends to be very close A is the surface of the lake (in m ) −8 W to 1. σ is the Stefan-Boltzmann coefficient ( 5.67e ) 2 4 m K ǫ is the lake emissivity A.4 Putting it all together T is the water surface temperature (here assumed to The final equation to compute surface heat losses is then: be 1 C less than the measured temperature) ǫ is the effective atmospheric emissivity a Q = Q + Q (1 + R ) e evap rad sh (20) T is the effective air temperature (here assumed to be and the mass-loss due to evaporation is: 0.9 C) M = Q /L e evap (21) A.2 Evaporation As proposed by Hurst et al. (2015) we use the equation of Adams et al. (1990) but replace the term for forced convec- tion by that of Sartori (2000): 1/3 2 0.8 −0.2 2 1/2 (17) Q =[(2.2(T − T ) (e − e )) + (L(0.00407u X )(e − e )/P ) ] evap sv av s a s a a 0 B ehr et al. Journal of Applied Volcanology (2023) 12:1 Page 9 of 10 Appendix B: Derivation of the Bayesian Filtering p(X , X |Y ) = p(X |X , Y )p(X |Y ) K K +1 1:t K K +1 1:t K +1 1:t and Smoothing equations p x , x |y = p x |x , y p x |y k k+1 k k+1 k+1 1:t 1:t 1:t In the following we will reiterate the derivation of the Bayesian Filtering and Smoothing equations presented Combining the joint and the conditional distributions by Särkkä (2013). and using again the Chapman-Kolmogorov equation we get the Smoothing Equation (Eq. 9) that solves Eq. 6: B.1 Bayesian Filtering p(x x )p(x y ) k+1 k k+1 1∶t For the Bayesian Filtering equation we only consider p(x y ) = p(x y ) dx k 1∶t k 1∶k k+1 p(x y ) k+1 1∶k observations up to time step k. Note that Eq. 9 is again computed recursively where p(x |y ) = p(y |x , y )p(x |y ) k k k (22) 1:k k 1:k−1 1:k−1 p(x |y ) is the solution to the previous smoothing step. k 1:t The recursion starts with the final step of the Bayesian Filtering recursion. = p(y |x )p(x |y ) k k (23) k 1:k−1 Abbreviations RCL Ruapehu Crater Lake where: M Lake water mass change M Steam input z = p(y |x )p(x |y )dx M Melt water and precipitation inflow k k k k k 1:k−1 M Lake overflow H Enthalpy and, using the Chapman-Kolmogorov equation: T Lake temperature c Specific heat of wat er Q Energy input from volcanic sources p(x |y ) = p(x |x )p(x |y )dx ˙ k 1:k−1 k k−1 k−1 1:k−1 k−1 Q Solar irradiation Q Sur face energy loss X Total ion amount in the lake Where we made use of Bayes’ rule for Eq. 22 and the Markov property (Eq. 7) for Eq. 23. Acknowledgements The authors would like to thank Nico Fournier for providing feedback and This is a recursive scheme where p(x |y ) depends 1:k ideas throughout the project, Sophie Pearson-Grant for her constructive and on the previous step’s solution p(x |y ) . To start k−1 1:k−1 insightful feedback on this manuscript, and Yuri Taran and one anonymous the recursion we have to guess p(x ). reviewer for their helpful comments. GNS Science’s volcano monitoring group has been an invaluable team to test and improve this work. Graphs have been produced using the open-source graphing libraries Plotly (https:// plotly. com/ B.2 Bayesian Smoothing python/), matplotlib (https:// matpl otlib. org/) and cartopy (https:// scito ols. Bayesian Smoothing takes into account all observations org. uk/ carto py/ docs/ latest/). Kalman Smoothing was implemented using the FilterPy library (https://github.com/rlabbe/filterpy). until time step t. This is sometimes referred to in the literature as history matching. Authors’ contributions Because the states x form a Markov chain, they are Y.B. wrote the main manuscript text, produced graphs, and implemented the Kalman Smoothing. S.S. worked on data collection and interpretation. T.H. independent of future measurements given x : k+1 developed and implemented the mass and energy balance model. All authors reviewed the manuscript. p(x |x , y ) = p(x |x , y ) k k+1 1:t k k+1 1:k Authors’ information Combining the Markov property with Bayes’ rule we can All authors are affiliated with GNS Science. TH is Emeritus Volcano Geophysi- cist, SS is Senior Volcano Geophysicist and YB is Volcano Science Operations express the probability distribution of x given x and k k+1 Specialist. y as follows: 1:t Funding p(x |x , y ) =p(x |x , y ) k k+1 k k+1 1:t 1:k This project was funded in parts by the Strategic Science Investment Funding p(x , x |y ) to GNS Science within the Data Science program and by the GeoNet program. k k+1 1:k p(x |y ) k+1 1:k Availability of data and materials The source code for this work can be accessed at https:// github. com/ yanni p(x |x , y )p(x |y ) k+1 k k 1:k 1:k kbehr/ pumahu. All data used in this study is provided free of charge by the p(x |y ) k+1 1:k GeoNet program (https:// www. geonet. org. nz). p(x |x )p(x |y ) k+1 k k 1:k p(x |y ) Declarations k+1 1:k Ethics approval and consent to participate The joint distribution of x and x given y can be k k+1 1:t Not applicable. expressed as: Behr et al. Journal of Applied Volcanology (2023) 12:1 Page 10 of 10 Consent for publication Scott BJ (1994) Cyclic Activity in the Crater Lakes of Waimangu Hydrothermal Not applicable. System. New Zealand. Geothermics 23(5–6):555–572. https:// doi. org/ 10. 1016/ 0375- 6505(94) 90019-1 Competing interests Scott BJ (2013) A Revised Catalogue of Ruapehu Volcano Eruptive Activity: The authors declare that they have no competing interests. 1830-2012. Technical Report 45. GNS Science, New Zealand Stevenson DS (1992) Heat Transfer in Active Volcanoes: Models of Crater Lake Systems. PhD thesis, The Open University Received: 19 April 2022 Accepted: 28 October 2022 Torrence C, Compo GP (1998) A Practical Guide to Wavelet Analysis. Bull Am Meteorol Soc 79(1):61–78. https:// doi. org/ 10. 1175/ 1520- 0477(1998) 079< 0061: APGTW A>2. 0. CO;2 Van der Merwe R (2004) Sigma-Point Kalman Filters for Probabilistic Inference in Dynamic State-Space Models. PhD thesis, Oregon Health & Science University References Vandemeulebrouck J, Stemmelen D, Hurst T, Grangeon J (2005) Analogue Adams EE, Cosler DJ, Helfrich KR (1990) Evaporation from Heated Water Bodies: Modeling of Instabilities in Crater Lake Hydrothermal Systems. J Geophys Predicting Combined Forced Plus Free Convection. Water Resour Res Res Solid Earth 110(B2). https:// doi. org/ 10. 1029/ 2003J B0027 94 26(3):425–435. https:// doi. org/ 10. 1029/ WR026 i003p 00425 White JT (2018) A Model-Independent Iterative Ensemble Smoother for Brown G, Rymer H, Dowden J, Kapadia P, Stevenson DS, Barquero J, Morales Efficient History-Matching and Uncertainty Quantification in Very High LD (1989) Energy Budget Analysis for Poás Crater Lake: Implications for Dimensions. Environ Model Softw 109(May):191–201. https:// doi. org/ 10. Predicting Volcanic Activity. Nature 339(6223):370–373. https:// doi. org/ 1016/j. envso ft. 2018. 06. 009 10. 1038/ 33937 0a0 Brown GC, Rymer H, Stevenson D (1991) Volcano Monitoring by Microgravity and Energy Budget Analysis. J Geol Soc 148(3):585–593. https:// doi. org/ Publisher’s Note 10. 1144/ gsjgs. 148.3. 0585 Springer Nature remains neutral with regard to jurisdictional claims in pub- Christenson BW, Reyes AG, Young R, Moebis A, Sherburn S, Cole-Baker J, Brit- lished maps and institutional affiliations. ten K (2010) Cyclic processes and factors leading to phreatic eruption events: Insights from the 25 September 2007 eruption through Ruapehu Crater Lake. New Zealand. J Volcanol Geotherm Res 191(1):15–32. https:// doi. org/ 10. 1016/j. jvolg eores. 2010. 01. 008 Dibble RR (1966) Volcanic Seismology. In: Thompson BN, Kermode LO, Ewart A (eds) New Zealand Volcanology, Central Volcanic Region, 1st edn. Depart- ment of Scientific and Industrial Research, New Zealand, pp 93–98 Fitzgerald SD, Woods AW (1997) The Stability of Two-Phase Geothermal Res- ervoirs. In: Proceedings of the 22nd Workshop on Geothermal Reservoir Engineering. Stanford, p 3 Fournier N, Witham F, Moreau-Fournier M, Bardou L (2009) Boiling Lake of Dominica, West Indies: High-temperature Volcanic Crater Lake Dynamics. J Geophys Res Solid Earth 114(2):1–17. https:// doi. org/ 10. 1029/ 2008J B0057 73 Friedlander B (1898) Some notes on the Volcanoes of the Taupo District. Trans. New Zealand Institute, 31:498–510 GeoNet Home. https:// www. geonet. org. nz/. Accessed 2022-02-01 Historic Eruptive Activity Mt. Ruapehu. (2022) GeoNet. https:// github. com/ GeoNet/ data/ tree/ main/ histo ric- volca nic- activ ity. Accessed 09 Mar 2022 Huang DZ, Schneider T, Stuart AM (2022) Iterated Kalman methodology for inverse problems. J Comput Phys 463:111262. https:// doi. org/ 10. 1016/j. jcp. 2022. 111262 Hurst AW, Bibby HM, Scott BJ, McGuinness MJ (1991) The Heat Source of Ruapehu Crater Lake; Deductions from the Energy and Mass Balances. J Volcanol Geotherm Res 46(1–2):1–20. https:// doi. org/ 10. 1016/ 0377- 0273(91) 90072-8 Hurst AW, Dibble RR (1981) Bathymetry, Heat Output and Convection in Ruap- ehu Crater Lake. New Zealand. J Volcanol Geotherm Res 9(2–3):215–236. https:// doi. org/ 10. 1016/ 0377- 0273(81) 90005-6 Hurst T, Christenson B, Cole-Baker J (2012) Use of a Weather Buoy to Derive Improved Heat and Mass Balance Parameters for Ruapehu Crater Lake. J Volcanol Geotherm Res 235–236:23–28. https:// doi. org/ 10. 1016/j. jvolg eores. 2012. 05. 004 Re Read ady y to to submit y submit your our re researc search h ? Choose BMC and benefit fr ? Choose BMC and benefit from om: : Hurst AW, Hashimoto T, Terada A (2015) Crater Lake Energy and Mass Balance. In: Lakes Volcanic (ed) Rouwet D, Christenson B, Tassi F, Vandemeuleb- fast, convenient online submission rouck J. Springer, Berlin, Heidelberg, pp 307–321 thorough peer review by experienced researchers in your field Kalman RE (1960) A New Approach to Linear Filtering and Prediction Problems. J Basic Eng 82(1):35. https:// doi. org/ 10. 1115/1. 36625 52 rapid publication on acceptance Mayhew YR, Rogers GFC (1978) Thermodynamic and Transport Properties of support for research data, including large and complex data types Fluids : SI Units, 2nd edn. Basil Blackwell, Oxford • gold Open Access which fosters wider collaboration and increased citations Rauch HE, Striebel CT, Tung F (1965) Maximum Likelihood Estimates of Linear maximum visibility for your research: over 100M website views per year Dynamic Systems. AIAA J 3(8):1445–1450 • Särkkä S (2013) Bayesian Filtering and Smoothing. Cambridge University Press, Cambridge. https:// doi. org/ 10. 1017/ CBO97 81139 344203 At BMC, research is always in progress. Sartori E (2000) A Critical Review on Equations Employed for the Calculation Learn more biomedcentral.com/submissions of the Evaporation Rate from Free Water Surfaces. Sol Energy 68(1):77–89. https:// doi. org/ 10. 1016/ S0038- 092X(99) 00054-7
Journal of Applied Volcanology – Springer Journals
Published: Feb 7, 2023
Keywords: Crater lakes; Non-linear Kalman Smoother; Volcanic heat emission
You can share this free article with as many people as you like with the url below! We hope you enjoy this feature!
Read and print from thousands of top scholarly journals.
Already have an account? Log in
Bookmark this article. You can see your Bookmarks on your DeepDyve Library.
To save an article, log in first, or sign up for a DeepDyve account if you don’t already have one.
Copy and paste the desired citation format or use the link below to download a file formatted for EndNote
Access the full text.
Sign up today, get DeepDyve free for 14 days.
All DeepDyve websites use cookies to improve your online experience. They were placed on your computer when you launched this website. You can change your cookie settings through your browser.