Kusatsu-Shirane volcano eruption on January 23, 2018, observed using JMA operational weather radars

A phreatic eruption suddenly occurred at Motoshirane (Kusatsu-Shirane volcano, Japan) at 10:02 JST on January 23, 2018. A member of the Japan Self-Defense Force was killed by volcanic blocks during training in Motoshirane, and 11 people were injured by volcanic blocks or fragments of broken glass. According to a field survey, ash fall was confirmed in Minakami, about 40 km east-northeast from Motoshirane. Although the eruption was not captured by a distant camera, the eruption plume/cloud was captured by three of the Japan Meteorological Agency’s operational weather radars. These radars observed the echo propagated to the northeast in the lower troposphere, and to the east in the middle troposphere. This is generally consistent with the observed ash fall distribution. Using the modified probabilistic estimation method, the maximum plume height was estimated to be about 5580 ± 506 m (1σ) above sea level. Estimates of the erupted mass based on the range of plume heights from radar observations and the duration of volcanic tremor during the eruption (about 8 min) do not match that obtained from a field survey (3.0–5.0 × 107 kg). This discrepancy confirms that estimates of erupted mass based on plume heights must account for eruption style parametrically, which can only be constrained by case studies of varied eruption styles.


Introduction
Weather radar is usually used to observe rain, but it is also a useful tool for detecting volcanic eruptions. In fact, many volcanic eruptions have been captured by weather radar (e.g., Sawada 2003;Marzano et al. 2013). However, it is difficult to estimate plume heights by using remote sensing instruments that transmit and receive radio waves, such as weather radar, because of error resulting from the influence of the beam width. To reduce this error, Sato et al. (2018) developed a probabilistic estimation (PE) method and applied it to the 2016 eruption of Mt. Aso. They derived the total amount of ejecta from the radar data and obtained a value consistent with a field survey. Nonetheless, it is important to confirm the validity of their method for other cases and eruption styles, particularly phreatic eruptions. For magmatic eruptions, several empirical formulas have been proposed to estimate the amount of ejecta using the plume height (e.g., Costa et al. 2016), but it is not clear whether they are also valid for phreatic eruptions. Therefore, one of the objectives of this paper is to test the validity of the empirical formulas for phreatic eruptions. Sato Earth, Planets and Space (2021) 73:117 eruptive activity is characterized by phreatic eruptions at Shirane, with most eruptions occurring around Yugama crater (e.g., Terada 2018). Therefore, many geochemical studies have been conducted on Kusatsu-Shirane volcano, especially on the hydrothermal system around Shirane (e.g., Ossaka et al. 1980Ossaka et al. , 1997Sano et al. 1994;Ohba et al. 1994Ohba et al. , 2000Ohba et al. , 2008Ohba et al. , 2019. Nurhassan et al. (2006), Matsunaga et al. (2020) and Tseng et al. (2020) revealed the three-dimensional structure of the Kusatsu-Shirane hydrothermal system using the magnetotelluric (MT) or audio-frequency magnetotelluric method. In particular, Matsunaga et al. (2020) reported the existence of an extensive low resistivity layer at depths of 1-3.5 km beneath the summit of Motoshirane using MT data collected in 2015 and 2016. This layer was interpreted as a hydrothermal fluid reservoir.
In 2018, the first historical eruption at Motoshirane occurred.

Phreatic eruption on January 23, 2018
The Coordinating Committee for Prediction of Volcanic Eruptions (CCPVE 2018) stated that an eruption occurred at Motoshirane around 10:02 JST on January 23, 2018. CCPVE also estimated that the eruption occurred along an east-west row of new craters 500 m long on the north side of the crater of the Kagamiike-kita pyroclastic cone (crater raw A). After that, another crater row (crater raw B) was found near Kagamiike (Japan Meteorological Agency (JMA) 2018). This eruption ejected large volcanic bombs over 1 km from the crater (CCPVE 2018). The eruption was accompanied by about 8 min of volcanic tremor with large amplitude that began at 09:59JST (CCPVE 2018). One Japan Self-Defence Force (SDF) member who was training at Motoshirane was killed during the eruption, and seven other members were injured nearby (Fire and Disaster Management Agency 2018). The eruption caused the ropeway at Kusatsu International Ski Resort to stop, stranding 81 people in the gondolas, four of whom were injured by a windowpane that was broken by the impact of volcanic bombs. According to a field survey conducted by the Joint Research Team for Ash Fall in Kusatsu-Shirane 2018 eruption (2018b), ash fall was also confirmed in Minakami, a town about 40 km east-northeast of Motoshirane.
The total erupted mass was estimated to be 3.0-5.0 × 10 7 kg by the Joint Research Team for Ash Fall in Kusatsu-Shirane 2018 eruption (2018a). As a result of the investigation of the ejecta, the National Institute of Advanced Industrial Science and Technology and the National Research Institute for Earth Science and Disaster Prevention (2018) concluded that the eruption had likely been a phreatic eruption. Yaguchi et al. (2019) also analyzed the volcanic ash of the 2018 eruption of Motoshirane, and showed that the explosion had occurred at a depth reaching the basement rocks (i.e., more than 100 m deep) and in a relatively gas phase-rich condition as compared with the 1982 eruption of Shirane.
Using the volcanic tremor data, Yamada et al. (2021) revealed that the 2018 eruption had been triggered by hydrothermal fluid injection through a different pathway from that of unrest activities at Yugama since 2014. Himematsu et al. (2020) analyzed the crustal deformation associated with the eruption using satellite synthetic aperture radar data. The dislocation plane inferred by them can be considered as a degassing pathway to the surface associated with the eruption.
A JMA camera at Ainomine, located 800 m from the crater raw A, did not observe the eruption, but three JMA weather radars (Nagano, Niigata, and Tokyo; Fig. 1) captured the echo of the eruption plume/cloud. The radar echo was observed to move northeast in the lower troposphere, and east in the middle troposphere (Fig. 2). These multi-station radar observations are consistent with the ash fall distribution reported by the Joint Research Team (2018a; b).

Modified probabilistic estimation method
In order to investigate the height of the eruption plume from its radar echo, the author developed a modified probabilistic estimation (MPE) method, an improved version of the PE method (Sato et al. 2018).
The path of radio waves in the spherically stratified atmosphere is expressed as (e.g., Hartree et al. 1946;Battan 1973): where h is the altitude above the surface of the Earth, s is the horizontal distance, R is the radius of the Earth, and n is the refractive index of the atmosphere.
Assuming (dh/ds) 2 ≪ 1 , n ≈ 1 and h ≪ R , Eq. (1) is simplified to: If the elevation angle θ e of the beam path is small, it can be assumed that dh/ds ≈ θ e , and Eq. (2) becomes: where k e is the equivalent radius ratio (e.g., Doviak and Zrnic 1993; Sato et al. 2018).
Whereas the PE method assumes the 4/3 R equivalent earth radius model, the MPE method incorporates the ellipticity of the Earth by using Euler's (1767) radius of curvature: Here, M ϕ is the radius of curvature of the meridian at latitude ϕ , N ϕ is the radius of curvature of the prime vertical, and α is the azimuth of the target relative to the radar. Using R α ϕ obtained for each radar, Eq. (3) becomes: We obtained dn/dh from the mean field of the lower atmosphere (< 5 km altitude) between the radar and the target by using the initial value of the JMA Mesoscale Model (MSM) at 0:00 UTC on January 23, 2018 as the atmospheric field. Since the MSM does not explicitly (1) predict the refractive index of the atmosphere, it was calculated as: where N is refractivity, p is atmospheric pressure, T is air temperature, and e is vapor pressure (e.g., Bean and Dutton 1968). Saturated water vapor pressure e s was approximated by using the expression used by the World Meteorological Organization (e.g., Mizuno 2000): where T c is the atmospheric temperature in degrees Celsius. In the MPE method, these formulas are used to obtain the equivalent earth radius ratio for each radar. Next, we derived the probabilistic estimation equations used in the MPE method. In the PE method, the radar beam height that captures the plume is calculated by: where h 0 is the antenna height. The JMA weather radars perform 28 scans per 10 min over a range of elevation angles from about 0 to 25 degrees. Furthermore, for the low elevation angles used in this study (about 0-5°), most of the angles are scanned twice, and as a result, the time resolution is about 5 min. Whether the target echo was caused by volcanic plumes or not is judged from the continuity between elevation angles.
One of the improvements from the PE method is a geoid correction. Whereas the PE method does not consider the geoid height for height estimations, the MPE method includes a geoid correction. Specifically, because the antenna height is usually given relative to the geoid, the geoid height at the location of the radar h G,radar is added to convert the antenna height to an altitude relative to the Earth ellipsoid. Then, after deriving the altitude of the target, the geoid height at the location of the target h G,target is subtracted to convert the target height to an altitude relative to the geoid. For the geoid height, GSIGEO2011 by the Geospatial Information Authority of Japan (2016) was used. By this correction, Eq. (8) becomes: The obtained heights are treated as altitudes above sea level (m ASL).
(6) N = (n − 1) × 10 6 = 7.76 × 10 −5 p T + 3.75 × 10 5 e T 2 , h = k e R + h 0 + h G,radar cosθ e cos θ e + s/(k e R) − 1 The probability density function of the heights obtained by each radar is calculated as: Here, h j = j h is the discretized height above sea level (j = 0,1,2…, ∞), and µ i and σ i are the median height and standard deviation of the ith radar, respectively. In addition, we assume that µ i = H c , σ i = β(H u − H l )/2 and β = 1 . H c is the altitude at the center of the beam, H u is the altitude at the top of the beam, and H l is the altitude at the bottom of the beam. Finally, the composite probability density combining all radars is obtained as:

Erupted mass estimation methods
One of the author's goals is to quickly estimate the scale of the eruption, that is, the erupted mass, using the eruption plume height estimated by weather radars. We used four empirical methods in Costa et al. (2016) to estimate the total mass of ejecta by the 2018 Kusatsu-Shirane eruption. For simplicity, each model is reproduced as closely as possible from the literature, and the variables in this section are independent of those used in the preceding section. We briefly describe each method and the references therein which provide more detailed information.
The where ρ a0 = 1.105kg/m 3 is the density of the atmosphere, g ′ = 41.289m/s 2 is the reduced gravity at the plume source, α = 0.1 is the radial entrainment coefficient, N = 0.0134s −1 is the average buoyancy, z l = 2.8 is the maximum non-dimensional height of Morton et al. (1956), β = 0.5 is the wind entrainment coefficient, and (10) (12) ln Ṁ = ln b 1 H n 1 + cWH, 3935m/s is the average wind velocity over the plume height. Mastin et al. (2009; hereinafter 'M09') estimated Ṁ as: where ρ m = 2500kg/m 3 is the density of magma. M09 estimated the volumetric flow rate (dense-rock equivalent) V m 3 /s at first, and converted it to Ṁ with this magma density. Finally, the method of Woodhouse et al. (2013Woodhouse et al. ( , 2016; hereafter 'W16') uses: and, Here, W s is a dimensionless wind strength, α = 0.09 , ρ a0 = 1.104kg/m 3 , θ 0 = 1273K and θ a0 = 268.8K are the source temperature and the atmospheric temperature, respectively, N = 0.014s −1 , γ = 0.007s −1 is the shear rate of the atmospheric wind, n 0 = 0.03 is the mass fraction of gas at the vent, C v = 1810J/ kg K , C s = 1100J/ kg K , and C a = 1000J/ kg K are the specific heat capacity at constant pressure of water vapor, solid pyroclasts, and dry air, respectively. Table 1 shows the range of eruption altitudes obtained by each radar and the final combined results. Using Eq. (5), k e is calculated from R α ϕ and dn/dh derived for each radar. The ranges of eruption altitude were calculated from Eq. (9) using k e . The MPE method assumes a probability density function with a normal distribution (Eq. 10), and derives the center of the beam (median value) and the standard deviation from it. The maximum plume height of the 2018 Kusatsu-Shirane eruption obtained by using Eq. (11) results in a composite estimate of 5580 ± 506 m ASL ( Fig. 3; Table 1). Because there is no large variation among the results obtained from each radar, we consider that anomalous propagation did not occur. Table 2 shows that estimates of the erupted mass by each empirical method, based on the range of the composite plume height in Table 1. All methods except DB12 underestimate the total erupted mass compared to the field survey.

Discussion
One of the reasons for the underestimation in most models (Table 2) seems to be related to the application of the empirical formulas created from magmatic eruptive events and less applicable to phreatic activity. The same tendency was confirmed when we applied the empirical formulas to the 2014 Mt. Ontake eruption, which is also a phreatic eruption. In this calculation, the maximum eruption height is assumed to be 8.0 km ASL (5.3 km above crater) obtained by a photograph (Sato et al. 2016).
It is also assumed that most of the ejecta erupted within about 25 min (Maeno et al. 2016). In other words, the duration of the 2014 Mt. Ontake eruption is assumed to be 25 min.
Only the DB12 method overestimated the erupted mass (Table 2), possibly because it overemphasizes the effect of entrainment by wind (i.e., the second term on the righthand side of Eq. (13). Suzuki and Koyaguchi (2015) estimated the entrainment coefficient β to be 0.1-0.3 from the calculation results of a 3-D model, therefore, there is a possibility that only the coefficient was a problem.
Nonetheless, none of the methods provided a satisfactory estimate of the erupted mass. Each method requires minor modifications according to the eruption style, in this case, a phreatic eruption. Therefore, it is necessary to accumulate far more eruption case studies and clarify the detailed relationships between the various physical quantities and the total erupted mass for the range of phreatic eruptions.

Summary and conclusion
The eruptive plume/cloud echo associated with the Kusatsu-Shirane eruption on January 23, 2018, was captured by three operational JMA weather radars. The echo moved northeast in the lower troposphere, and east in the middle troposphere. These radar observations were consistent with the observed ash deposition patterns. Using the MPE method, the plume height was estimated to be about 5580 ± 506 m ASL (1σ). The amount of the ejecta estimated from this result and the duration of volcanic tremor associated with the eruption (about 8 min) did not match the mass estimated by the field survey (3.0-5.0 × 10 7 kg). More specifically, the estimates of the erupted mass by the three methods (C14, M09 and W16) were less than half of the field survey, while only DB12 overestimated it. One of the reasons seems to be that these empirical formulas are not applicable to phreatic eruptions without change. This implies that estimates of the total erupted mass based on plume height must account for eruption style by selecting reasonable model parameter values, which can only be constrained by further case studies incorporating varied eruption styles.    Woodhouse et al. (2013Woodhouse et al. ( , 2016 1.2 × 10 7 -3.9 × 10 7 Field survey 3.0 × 10 7 -5.0 × 10 7