Crisis hazard assessment for snow-related lahars from an unforeseen new vent eruption: the 2018 eruption of Kusatsu-Shirane volcano, Japan

Two-thirds of the 111 active volcanoes in Japan are covered with snow for several months during winter and demonstrate high hazard and risk potentials associated with snow-related lahars during and after eruptions. On 23 January 2018, a sudden phreatic eruption occurred at the ski field on Kusatsu-Shirane (Mt. Motoshirane) volcano, Japan. This new vent eruption from the snow-clad pyroclastic cone required forecasting of future snow-related lahars and crisis hazards zonation of downslope areas including Kusatsu town, a popular tourist site for skiing and hot springs. In order to achieve a prompt hazard assessment for snow-related lahars, a multidisciplinary approach was carried out involving characterization of proximal tephra deposits, snow surveys, and numerical lahar flow simulations using the Titan2D model. To determine the input parameters for the flow model, the consideration of snow water equivalent (SWE) immediately after the eruption (on 29 January) and in the post-eruptive period (on 12 March), was significant. In the case of Kusatsu-Shirane volcano during the winter of 2018, linear relationships between altitude and SWE, obtained at different elevations, were used to estimate the snow volume around the new vents. Several scenarios incorporating snow and snowmelt (water), with or without the occurrence of a new eruption, were simulated for the prediction of future lahars. Three lahar scenarios were simulated, including A) rain-on-snow triggered, B) ice/snow slurry, and C) full snowmelt triggered by a new eruption, and indicated the flow paths (inundation areas) and travel distances. These were useful for lahar hazard zonation and identification of potential high-risk areas. Since the input parameters required for the Titan2D flow model can be relatively easily determined, the model was suitable for the 2018 eruption at Motoshirane where historical and geological lahar records are not available for calibration. The procedure used in the study will enable rapid lahar prediction and hazard zonation at snow-clad volcanoes. Further consideration for simulating a cohesive-type flow, which was predicted by the primary deposits containing large amounts of clay minerals and could not be expressed in the Titan2D flow model, is necessary.


Introduction
Two-thirds of the 111 active volcanoes in Japan are seasonally snow-covered during winter (occasionally up to 6 months) and demonstrate high hazard and risk potentials associated with snow and snowmelt during and after an eruption (Tada and Tsuya 1927;Waythomas 2014;Kataoka et al. 2018). In Japan, a snowmelt lahar, triggered by a hydrothermal eruption and associated flank collapse at Tokachi-dake, Hokkaido in 1924 caused 144 fatalities (Tada and Tsuya 1927;Uesawa 2014). Based on lessons learned from this event and from a snowmelt lahar tragedy at Nevado del Ruiz, Columbia in 1985 (Pierson et al. 1990), most snow-related lahar risk assessments and hazard maps in Japan are based on a single scenario in which snowmelt is directly triggered by a magmatic and pyroclastic eruption, such as that at Nevado del Ruiz. However, the types and triggering mechanisms for snow/ ice-related mass flows vary and include volcanic mixed avalanches triggered by pyroclastic flows sweeping across snow (Pierson and Janda 1994), and ice slurry lahars initiated by phreatic eruptions (Kilgour et al. 2010), as reported worldwide. The trigger, size, and type of mixed flows (including lahars) of volcanic sediments and rocks, snow, and/or water can vary with the timing of the event, whether it occurs during the rainy season, snow season, or snowmelt season. The mechanism by which contrasting lahars occur at seasonally snow-clad volcanoes has been evaluated by a study conducted on lahar deposits and via river monitoring, after the September 2014 phreatic eruption at Ontake volcano, central Japan . At Ontake, a rain-triggered lahar occurred on 5 October 2014, causing the formation of clay-rich debris flow deposits (clay content of 10-20 wt%). In contrast, a lahar flow that occurred in April 2015 under rain-on-snow (hereafter ROS : Kattelmann 1997;Sui and Koehler 2001) conditions during snowmelt season was reportedly more water-rich and erosive, resulting in the formation of a clay-poor hyperconcentrated flow deposit . Such different flow types (e.g., cohesive or non-cohesive; debris flow or hyperconcentrated flow) may exhibit different characteristic travel distances, travel times, inundation areas, and flow transformations. Because lahars at snow-clad volcanoes can have these many variations, hazard assessment and zonation of lahars should be achieved by considering several lahar scenarios for different (i.e., snow, snowmelt, and rainy) seasons.
On 23 January 2018, a sudden phreatic eruption occurred at the ski field on Kusatsu-Shirane (Mt. Motoshirane) volcano in Honshu, Japan (Ogawa et al. 2018). The eruption originated from new vents located near the Kagamiike-kita cone at Motoshirane (Figs. 1 and 2), which had been considered dormant Ishizaki et al. 2020). Therefore, no volcanic hazard maps (including lahar hazards) based on an eruption from Motoshirane had been published prior to the 2018 eruption. Previous geophysical and geochemical studies have intensively focused on another peak (Mt. Shirane) of Kusatsu-Shirane volcano around the Yugama crater in Fig. 1 (e.g., Nurhasan et al. 2006;Ohba et al. 2019;Tseng et al. 2020). Kusatsu, a town located on the downslopes of the volcano, is extremely popular for sightseeing, skiing, and onsen (hot spring bathing); owing to such aspects, people come in close proximity to a volcanic hazard-prone area. Due to an elevated volcanic eruption risk since the January 2018 eruption, the hazards associated with future eruptions and lahars occurring during both snow and snowmelt seasons have garnered attention as major concerns for safety of the town, local residents, and inhabiting communities.
Most lahar hazard maps reported worldwide are background (long-term) maps prepared before eruption-associated unrest. Lahar hazards on glaciated stratovolcanoes in North and South America have been well researched using historical and geological evidence (e.g., Mothes 1992;Pierson 1999;Delgado-Granados et al. 2015). In contrast, mapping of lahar hazards associated with sudden eruptions of seasonally snow-covered volcanoes such as Motoshirane is extremely rare.
Titan2D, a numerical flow simulation tool Pitman and Le 2005), has been used mostly for hindcasting of past lahar events (e.g., Williams et al. 2008;Córdoba et al. 2015). It has also been used in forecasting events prior to the 2007 Crater Lake outburst lahar at Ruapehu, New Zealand (Procter et al. 2012) and during the 2012 Tongariro eruption crisis documented in New Zealand: the latter had significantly less historical or geological data available (Leonard et al. 2014). The present study applied Titan2D to lahar hazard forecasting at Motoshirane, without historical and geological records of past lahars. In this study, snow-related lahar mass flows were simulated by considering the implications of snow conditions and seasonality. To constrain the parameter values for the flow simulations, we characterized primary tephra deposits that can be a source material for future lahars and conducted snow surveys to identify snow conditions and accumulations. A scenario-based hazard zonation considered potential lahar triggers, such as the type of future eruptions, and anticipated rain, snow, and snowmelt conditions. The process described in the present study for determination of scenarios and for estimation of parameter values will be helpful for forecasting snow-related lahars, without substantial reference geological and historical lahar data, during a volcanic crisis caused by an unforeseen new vent eruption on a seasonally snow-clad volcano.

Regional setting of Kusatsu-Shirane volcano and the 2018 phreatic eruption
Kusatsu-Shirane volcano (2160 m above sea level (a.s.l.) at Mt. Shirane, and 2171 m a.s.l. at Mt. Motoshirane) is located on the border of Nagano and Gunma prefectures on Honshu Island. The summit of the volcano is predominantly composed of two clusters of young pyroclastic cones and craters, the Shirane group and the Motoshirane group Ishizaki et al. 2020;Figs. 1 and 2), some of which are occupied by crater lakes. Greater active volcanism has been observed in the Shirane group with at least 10 phreatic eruptions reported during the twentieth century. Motoshirane was dormant until the phreatic eruption in 2018. Geological records of past volcanic activity at Motoshirane suggest that the magmatic eruptions that led to the formation of the Kagamiike-kita cone occurred approximately 1500 years ago, while the Kagamiike cone was formed approximately 4800 cal yr BP (Ishizaki et al. 2020  On 23 January 2018 (10:02 a.m. Japanese Standard Time), a phreatic eruption suddenly occurred at vents near the Kagamiike-kita cone (ca. 2050 m a.s.l.), part of Mt. Motoshirane. Tremors and tilting were recognized immediately before and after the eruption (Japan Meteorological Agency 2018; Terada et al. 2021;Yamada et al. 2021). Several new vents were opened, ejecting tephra and ballistic projectiles estimated to total between 30,000 and 50,000 metric tons (The Joint Research Team for ash fall in Kusatsu-Shirane 2018 eruption 2018; Kametani et al. 2021). Fumaroles were observed until the end of February at vents located north of the Kagamiike-kita crater.
Proximal tephras, including ballistic projectiles, were reported to be the direct cause of one fatality and injured 11 people who were skiing or were in cable cars at the Kusatsu-Kokusai (now renamed as Kusatsu-Onsen) ski field (Fire and Disaster Management Agency 2018). The medial to distal tephra distribution was investigated by several research groups immediately after the eruption (The Joint Research Team for ash fall in Kusatsu-Shirane 2018 eruption 2018; Kametani et al. 2021). Airborne and ground surveys indicated that the tephra was distributed to the northeast location of the vents (Fig. 2).
Creeks and rivers where potential lahars could have occurred after the 2018 eruption were identified in this study ( Figs. 1 and 2). The head of Furiko-zawa Creek is situated on the northern slope of the Kagamiike-kita cone. This creek runs through a valley formed at the edge of a lobe of Furikozawa lava (Hayakawa 1983;Uto et al. 1983) and was the main component of a ski course. Shimizu-sawa Creek, a part of which was also used as a ski course, spans from 500 m east of the Kagamiike-kita The topographic map is based on open source data provided by the Geospatial Information Authority of Japan (GSI). Distribution of tephra and isopleths and dark gray-colored and limit of gray-colored areas indicating ash covered surfaces are based on a preliminary report by the Geological Survey of Japan (released on January 26, 2018) and the Joint Research Team for ash fall in Kusatsu-Shirane 2018 eruption (released on 6 February 2018). New vent locations after GSI data (https:// www. gsi. go. jp/ BOUSAI/ kusat sushi rane-index. html). b Schematic longitudinal profiles of the Furiko-zawa and Shimizu-sawa creeks and the Yazawa-gawa River cone, flows in the ENE direction, and joins Furiko-zawa Creek in the middle of the ski field. After the confluence, Furiko-zawa Creek flows down the lower ski slope and finally joins the Yazawa-gawa River, the upper catchment of which is located on the eastern slope of the Yugama crater of Mt. Shirane ( Figs. 1 and 2). A lahar hazard map is based on the assumption of lahar generation from Yugama occurring as an event with a probability of 1-in-200 years. This map was prepared before the 2018 eruption, when there was no expectation that lahars could originate from an eruption at Motoshirane.
Glaciers are absent, but the Kusatsu-Shirane volcanic region is seasonally snow-clad from November to May. The climatological normals for 1981-2010 over three winter months (December, January, and February) at the Kusatsu weather station (1223 m a.s.l., in center of Kusatsu town) operated by Japan Meteorological Agency showed a mean air temperature of − 3.1 °C and precipitation of 183.5 mm. For the three spring months (March, April, and May) mean temperature is 5.5 °C and rainfall is 346.3 mm. During the rainy season from mid-June to mid-July, this station receives, on average, 370.8 mm of rainfall. Before the 2018 eruption, the maximum daily precipitation recorded at the Kusatsu weather station was 222 mm in 1982; however, it was surpassed by the record of 250 mm in 2019.

Procedure and methods for lahar hazard assessment
A crisis (short-term) hazard map for potential lahar hazards and risks after the 2018 eruption was generated using a stepwise procedure (Fig. 3). As the eruption occurred in the middle of the snow season, consideration of various conditions of snowfall and snowmelt, which can affect lahar generation and flow behavior, was essential (Kilgour et al. 2010;Córdoba et al. 2015;Kataoka et al. 2018). The approach (Fig. 3) consists of (1) characterization of eruption deposits (petrography, grain size, clay mineralogy, and physical properties); (2) meteorological observations and snow surveys, and (3) flow simulations using actual values for parameters obtained via (1) and (2). Additionally, (4) real-time river monitoring was conducted for the hazard-potential Furiko-zawa Creek and Yazawa-gawa River (Fig. 2) to obtain data on flow features and to estimate the travel times in the event when a lahar is generated during the post-eruption period (Appendix 1).

Characteristics of the 2018 eruption deposits Materials
Proximal tephra deposits are considered the most important information source for prediction of lahar flow type. The primary tephra deposits were sampled on 24 January 2018 (one day after the 23 January eruption) from no. 29 cable car at the Kusatsu-Kokusai ski field. The cable car was operational and was located ca. 440 m NNE of the new vents during the eruption (Location A; Fig. 2); it was damaged by ballistic ejecta (Fig. 4a). Owing to the proximity of the car to the vents, it could be ascertained whether the transport process of tephra occurred by fall, pyroclastic density current, or both. The tephra deposits that covered the seats in the car and clung to its window rims ( Fig. 4b) were collected, and the following analyses were performed.
Numerical lahar flow simulations using Titan2D (of single and two-phase models)

Petrography, grain size, clay mineralogy, and physical property of the eruption deposits
For petrography, the tephra samples were observed under binocular and petrographic microscopes. The samples were then subjected to washing over a sieve with a 63-µm mesh size to remove the mud fractions. Grain-size analysis was performed as follows: the tephra samples were dried, weighed, soaked in distilled water, and then wetsieved. Grains larger than 63 μm (sand and gravel) were sieved on a phi scale, whereas grains smaller than 63 μm (i.e., silt and clay) were analyzed using a laser grain-size analyzer (Malvern Mastersizer 3000) at Niigata University. For the particle analyzer, samples were ultrasonically stirred for 2 min, and five randomly collected portions of the sample were measured five times. The duration of each measurement was set at 15 s. The mean values of all measurement were considered for grain-size distribution. The gray-colored eruption deposits observed in no. 29 cable car were 7-9 mm thick and composed of armored lapilli-like aggregates (i.e., individual grains coated with mud) up to a few millimeters in size (Fig. 5a). Pebble-to cobble-sized volcanic rocks were scattered on the floor of the car. The deposits mainly consist of white hydrothermally altered rock fragments, quartz, feldspar, and extremely fine pyrite grains (Fig. 5b). Glassy rock fragments are rarely found. Visually, although the deposits seemed to be moderately sorted, grain-size distributions ( Fig. 6) indicate poor sorting; however, this is probably attributable to the mud content. The mud population (< 0.063 mm) occupies 42-48 wt% of the deposits  and 45-50 wt% of finer fractions (< 2 mm). Clay-sized (< 4 µm) particles, mostly derived from hydrothermally altered materials are ~ 10 wt%.
X-ray diffraction (XRD) analysis (Rigaku Ultima IV at Niigata University) was performed to understand the mineral composition and to determine the clay mineralogy. After separation by hydraulic settling and centrifugation-based separation, bulk and fine fraction (clay mineral) analyses were conducted. For clay minerals, analysis was performed as follows: (1) orientation, (2) treatment with ethylene glycol, (3) heating at 400 °C, and (4) heating at 550 °C for identification of the minerals with overlapping of specific cell parameters. The measurement conditions included exposure to CuKα radiation, 40 kV and 40 mA, and a scanning speed of 2°/min.
For numerical flow simulations (described in the latter sections), the internal friction angle of the flows is necessary. For lahar flow simulations, the angle of repose of deposited materials was used as a surrogate for internal friction angle (e.g., Procter et al. 2010). In this study, the angle of repose (mean of four measurements) of tephra deposits, snow grains (1-2 mm in size), and admixture of tephra with snow was measured using a protractor in a cold room (under 0 °C) at the Research Institute for Natural Hazards, Niigata University. The angles of repose for dried tephra and snow grains were 33.6° and 39.3°, respectively. The angle of repose for a mixture of ca. 20 wt% of snow and 80 wt% of tephra was 36.3°.

Snow survey
The snow water equivalent (SWE: amount of the liquid water that would be released upon complete melting of the snowpack) around the new vents was estimated, based on the relationship between altitude and SWE at various sites on the slopes of Kusatsu-Shirane volcano, because in situ measurements of SWE around the vents immediately after the eruption were not possible. Field measurements of SWE were conducted at four sites (from 1059 to 1512 m a.s.l.) on the eastern side and at four sites (from 1150 to 1977 m a.s.l.) on the western side of the volcano on January 29, and again at five selected sites among the above-mentioned sites on March 12. A snow sampler (a metal tube with a cutter) was used to obtain the snow core samples, and the samples were then weighed on an electric balance. SWE assessment can help provide values for the input parameter of thickness (i.e., volume) of snow and snowmelt for numerical simulations. In the case of Kusatsu-Shirane volcano during the winter of 2018 (Fig. 8), the relationships between altitude H (m) and SWE (mm) on January 29 and March 12 can be described as follows:  During the snowmelt season, snow pit observations of snow stratigraphy, grain shape and size, density, hardness, and liquid-water content were performed on 10-11 April 2018 along the National Highway 292. At Location B (Fig. 2: 1838 m a.s.l.), on the proximal slope, NE of Motoshirane, a snowpack of ~ 0.6 m thickness above the 2018 tephra layer and a total snow thickness of 2.75 m were observed (Fig. 9a). A 10-cm-thick ice layer was observed immediately above the 2018 eruption deposits (depicted in dark gray, 2 cm thick), which, in turn, was overlain by a water-saturated snow layer (Fig. 9b). The tephra and ice layers were impermeable, thereby diminishing the infiltration capacity and enhancing the anticipated surface runoff from the snow layer, above the impermeable basal boundary, when it was saturated with water from melting snow and heavy rain. These observation occurred after the initial lahar hazard zonation because of the limitation imposed on visiting high-risk areas. However, the result confirmed the parameter setting in the scenario models for snowrelated lahars.

Flow models of Titan2D
Titan2D is a numerical model ) that is used for studying dry granular flows and is based on the Savage-Hutter model (Savage and Hutter 1989). One of the assumptions of the Savage-Hutter model is where V x and V y represent the flow velocities in the x and y directions, g x and g y represent gravity accelerations in the x and y directions, and φ int and φ bed represent the internal and bed (basal) friction angles. Titan2D provides the option of using a two-phase flow model developed by Pitman and Le (2005). This model is also based on the Savage-Hutter model, but with the consideration of solid and fluid phases. In practice, the two-phase model requires the inclusion Fig. 9 Snowpack with intercalation of the 2018 tephra fall deposits on the northern and proximal slope nearby Motoshirane (1838 m a.s.l.) on April 10, 2018 in the snowmelt season. a A 2.75 m depth snowpack with patchy tephra exposure on the snow surface, and b A-10-cm thick ice layer developing above the 2018 eruption deposits (dark gray, 2 cm thick), overlain by water-saturated snow layer. Partly dark-colored part below the ashfall deposits is due to the wet snow (not ash infiltration); it became darker when water flowed from the above layer during digging the pit of the volume fraction of the solid phase as an input parameter. The two fluid system equations are based on those reported by Anderson and Jackson (1967), and the mass conservation equations for the two phases are as follows: where ρ s and ρ f represent the specific densities of the granular material (solid) and interstitial fluid, respectively, and ϕ is the solid volume fraction. u and v represent the velocities of the fluid and solid phases, respectively. The momentum conservation equations are expressed as: where T s and T f represent solid and fluid stress tensor, respectively. f represents the interaction force which is expressed as: The phenomenological function β is based on the experimental results of Richardson and Zaki (1954) expressed as: where v T is the terminal velocity of a typical solid particle falling in the fluid under gravity, g is the magnitude of the gravitational force, and m represents the Richardson and Zaki exponent related to the Reynolds number of the flow. The solid is assumed to be a Mohr-Coulomb incompressible granular fluid, as modeled by Savage and Hutter (1989). The model is then based on the assumption of a frictional boundary condition on the smooth basal surface, b(x, y) . The fluid is assumed to be inviscid; however, it exerts a drag at the basal surface as a Navier slip at b, which accounts for fluid dissipation. After averaging the values for depth and simplifying certain assumptions, the entire set of equations with model parameters and their notations are listed in Appendices 2 and 3. A simulated flow using Titan2D is commenced spontaneously with an initial pile that collapses because of its morphology and instability as a granular pile, and the simulation ends when it reaches the maximum number of iterations (10,000 iterations were performed in the initial phase of this study). The maximum number of iterations was set to 20,000-110,000 as this rendered the simulation time close to the estimated flow duration (Table 1). In this study, we used a cylindrical pile because it is easy to calculate an initial radius for the cylinder from the initial estimated volume. The input parameters of the Titan2D simulations include the internal friction angle, bed friction angle, initial pile height, initial pile radii, and volume fraction of the solid phase. A digital elevation model (DEM) is also necessary to calculate the granular surface flow.
The DEM was provided by the Geospatial Information Authority of Japan (GSI). The original DEM was used with a grid at 1 m; however, it was resampled to a 5-m

Scenarios and input parameters for flow simulations
In this study, field and laboratory data were used to determine the values of the input parameters for the simulations (Fig. 3). Historical and geological records of past lahars at Motoshirane as well as data on the distributions, depositional characteristics (thickness, grain size, cohesiveness, and so on), and triggering mechanisms necessary for the calibration of the simulations were not available. Therefore, several scenarios were considered, and each involved the triggering of a lahar due to snow and snowmelt interactions (Fig. 10, Tables 1 and 2).
The two-phase model used in Titan2D simulations requires the initial pile height and volume (i.e., area), the solid fraction, internal friction angle, and bed friction angle. In the simulated flows, the initial pile height and volume are equivalent to the entrained materials at the starting point of the flow, including tephra, snow, and water (rain). In this study, snow was considered as either a solid (particles) or a liquid/fluid (after melting) state of matter depending on the conditions and lahar triggers (discussed below). Tephra volumes in the proximal areas, which can be source material for a lahar, were considered based on a preliminary investigation reported on February 6, released by The Joint Research Team for ash fall in the Kusatsu-Shirane eruption (2018) consisting of groups belonging to the Geological Survey of Japan (GSJ), Earthquake Research Institute, the University of Tokyo (ERI), and National Research Institute for Earth Science and Disaster Resilience (NIED). The tephra distribution and volume estimation were referred from individual isopachs/areas of 1.5 m thick (GSJ) or 0.51 m thick (ERI) tephra deposition over an area of 100 m × 100 m, and of > 1 m thickness over an area of 70 m × 40 m (NIED) near vents. For snow accumulation, the snow survey helped determine volumes above the 2018 tephra layer (Figs. 8 and 9).
During the snowmelt season, a rain-triggered lahar can occur; therefore, the rain condition for the ROS lahar after the 2014 eruption at Ontake (3067 m a.s.l.), which recorded a total input ~ 400 mm (300 mm rain + 100 mm snowmelt) in 44 h (~ 2 days) at the Tanohara weather station (2195 m a.s.l.; Kataoka et al. 2018), was referred to set a maximum input value. In addition, the recorded maximum daily precipitation was 222 mm in Kusatsu town (1223 m) before the 2018 eruption. In near-vent areas (2050 m), more rainfall can be expected than in the lower elevated areas.
Internal friction values were based on laboratory observations of the angles of repose of the tephra and tephra and snow-mixed materials. Procter et al. (2010) reported that internal friction angles might affect the lateral movement of the flowing mass, and their results (between 30° and 35°) accorded with the representative angle of repose of debris flow deposits. In addition, the actual values of the internal friction angles in the Titan2D model are not sensitive for the simulation results and may subtly affect the lateral confinement of the flow Williams et al. 2008;Procter et al. 2010). The bed friction angle was variable and subjected to adjustments through repeated simulation runs, but was also based on previous lahar and snow avalanche simulations using Titan2D (Williams et al. 2008;Procter et al. 2010Procter et al. , 2012Takeuchi et al. 2018).
Three different scenarios were selected for the lahar simulations (Fig. 10, Table 1). Scenario A is based on the assumption of a non-eruption trigger (i.e., post-2018 eruption), whereas Scenario C is based on the flow triggered by a future eruption event. The intermediate setting, scenario B, is based on the assumption of both post-2018 and future (new) eruptive events.  Scenario A: A rain-triggered situation, at the end of the snow season and without further eruption, was considered. It was based on the assumption of an ROS condition with an input of approximately 400 mm per event. The SWE (Fig. 8) suggested that ca. 0.21 m of meltwater, above the tephra layer, was generated at 2050 m a.s.l. Considering that flow generation could occur widely, the representative value was set as 0.15 m. As a result, the solid fraction contained only tephra. The initial pile was set on the northern slope of the Kagamiike-kita cone overlapping the upper catchment of Furiko-zawa Creek, where proximal tephra deposits of approximately 1 m thickness were thought to be concentrated.
Scenario B: This scenario, for the middle to the end of the snow season, was based on the assumption of an ice/snow slurry triggered by melting of snow under warm temperatures and/or a small phreatic/hydrothermal eruption (e.g., Kilgour et al. 2010). The hypothesized slides were generated from the boundary above the tephra layer, and this assumption was validated later by snow pit observation (Fig. 9b). These slides can occur over a relatively wide area on the northern slope near the vents, resulting in the formation of the largest initial pile among the three scenarios. During the snowmelt season, patchy exposure of tephra and/or snow-covered slopes occurs (Fig. 9a) in areas proximal to the vents. Surface albedo modification by dark-colored ash cover may enhance snow and ice melt (Driedger 1981;Manville et al. 2000;Richardson and Brook 2010); however, thicker tephra layers may reduce ablation (Mattson et al. 1993;Richardson and Brook 2010) as the insulation reduces melting (Mattson et al. 1993;Brock et al. 2007). Therefore, the thickness of tephra was set as 0.5 m for this scenario. Water input was estimated to be ca. 0.21 m based on the infiltration of snowmelt water from the surroundings and upper catchment above the tephra layer (calculated by SWE at 2050 m a.s.l.; Fig. 8) and was augmented by a 100-mm input (assuming the intensive rainfall). In this scenario, snow can behave as a solid in the ice/snow slurry. The snow proportion was estimated by considering the same SWE, which was converted to snow depth (multiplied by 2.8, when snow density is ~ 360 kg/m 3 ).
Scenario C: This scenario was based on the assumption of a "full" snowmelt lahar generated by a minor magmatic, pyroclastic eruption during the middle to the end of the snow season; considering an event with similar volumes of proximal tephra fall deposits to the 2018 phreatic eruption. The lahar was sourced at the eruption center; thus, a narrow area was adopted for the initial pile. The tephra thickness contributing to the total pile height was considered to be 1.5 m (the maximum thickness suggested by GSJ). All snowpacks, including the tephra layer near the vents, resulted in 0.73 m of snowmelt water at 2050 m a.s.l. In this scenario, there was no additional rainfall input, and only melted snow contributed to the fluid phase of the simulated flows.
The solid fraction for each scenario was > 0.6 (Table 1), implying that the lahars of debris flow type following a definition based on sediment concentration (Scott et al. 1995;Pierson 2005). As a consequence of pile height estimation, a solid fraction of approximately 0.8 was used to simulate the snow slurry of Scenario B. However, the plausibility of this value should be assessed further, because empirically, a non-volcanic slushflow (Hestnes 1998) exhibits more watery characteristics, whereas the hydrothermal eruption-triggered ice/snow slurry lahar at Ruapehu demonstrated a volume of snow that was at least 60 times greater than the expelled water from the Crater Lake (Lube et al. 2009).
For comparison, the present study was conducted using the single-phase Coulomb model, which fundamentally applies to single-phase granular flows such as snow avalanches and debris avalanches (not wet) (Manville et al. 2013), although it has been previously applied for lahar hazard assessments (Procter et al. 2004(Procter et al. , 2012. For simplification, the sum of the solid and fluid volumes in individual scenarios was considered as the initial volume for the single-phase model used in this study. Additionally, the peak flow discharge and the duration for the hypothesized lahars in individual scenarios (Table 1) were estimated using the empirical methods described by Pierson (1998). The calculated durations were used to determine the maximum number of iterations for the simulations.

Comparison of the results
The simulated flow distributions are illustrated in Figs. 11, 12 and 13. In terms of the bed friction angle, the conduction of several two-phase flow simulations suggested that a value of 25° was representative of the present situation. Additionally, a simulated flow using the single-phase Coulomb model, with the same bed friction angle, covered a markedly shorter distance (Fig. 11b). Simulations using the two-phase model with lower friction angles of around 10° of previous lahar simulations (Williams et al. 2008;Procter et al. 2010Procter et al. , 2012 have been conducted; however, this resulted in numerical instability in the simulated flows in this study. Other cases were calculated using the single-phase Coulomb method, with minimum bed friction angles of 5° for Scenario A and 8° for Scenarios B and C. The simulations were based on the adoption of minimum bed friction angles for estimation of the maximum runout distances of the lahars. The minimum angles were defined by decreasing the values by 1° steps until the flow results showed instability, including flows progressing towards highly elevated zones. This study then considered application of the identified Shimizu-sawa C.

Kusatsu Town
a Two-phase, bed 25 c Single-phase, bed 5

Fig. 11
Simulated flow distribution of the Scenario A (rain-on-snow), a two-phase with the bed friction angle of 25°, b single-phase with the bed friction angle of 25°, and c single-phase with the bed friction angle of 5° minimum bed friction angle plus 1° for the conduction of single-phase simulations. The Scenario A (an ROS condition) lahars demonstrated a downward flow from the Furiko-zawa Creek until the cable car terminal, with a maximum flow thickness of ~ 0.5 m, consequently reaching the Yazawagawa River (Fig. 11a). A lahar flow simulated using the single-phase model (Fig. 11c) was found to be more valley-confined and demonstrated a maximum flow depth of ~ 1 m, even in the Yazawa-gawa. In Scenario B (snow-slurry), lahars bifurcated and reached the Furikozawa Creek and into a valley south of the Furiko-zawa. The two branches converged in the middle of the ski field and finally flowed into the Yazawa-gawa River (Fig. 12). The maximum flow thickness was almost steady and was approximately 0.5-1 m from upstream to downstream.
Around the cable car terminal, the topography split the lahar into two; one flowed further eastward, reaching Mononuguno-ike pond and the other flowed through a narrow valley to the Yazawa-gawa River. The flow paths of Scenario C (magmatic, small-scale eruption-triggered snowmelt lahar) under conditions of application of the two-phase model were similar to those of Scenario B with a maximum flow thickness of approximately 0.5 m, whereas the flow simulated using the single-phase model was similar to that observed in Scenario A (Fig. 13). Furiko-zawa C.
Kusatsu Town NH292 b Single-phase, bed 8  (Figs. 1, 2, 11a, c, 12 and 13). These locations are in ski areas during winter and are near the major transportation route during other seasons. Additionally, all numerical simulations conducted using variable input parameters suggest that there is a low probability of the occurrence of lahars directly reaching Kusatsu town and/or spilling from the Yazawa-gawa River (flowing in a deep gorge), and finally inundating the town. The graphs of travel distance versus time (Fig. 14) depict breakpoints in the lines, these indicate changes in the flow velocity at certain times. Therefore, the flow velocity was calculated separately in the segments before and after the breakpoint (Table 3). The breakpoints occur ca. 2500 m downstream, reflecting the influence of topography such as the valley becomes widen and the slope gradient changes around the location.

Snow-related lahar hazard assessment without the availability of geological and historical data
The 2018 eruption occurred at new unforeseen vents on the snow-clad Motoshirane. Geological and historical lahar data for probabilistic assessment of snow-related lahar hazards were unavailable. Therefore, a deterministic approach, using the characteristics of proximal tephra Furiko-zawa C.
Kusatsu Town NH292 Fig. 13 Simulated flow distribution of the Scenario C (full snowmelt): a two-phase with the bed friction angle of 25° and b single-phase with the bed friction angle of 8° deposits and snow conditions, was adopted and flow was simulated using a numerical model.

Flow type of the predicted lahars after the 2018 phreatic eruption
Prediction of the flow type (cohesive, non-cohesive, debris flow, hyperconcentrated flow, and so on) of potential future lahars, even at qualitative levels, is important because it affects flow energy, flow transformation, travel distance, inundation areas, and so on, and these affect the outcomes of hazard mitigation strategies implemented, including evacuation planning. Syn-and post-eruption lahars tend to exhibit the characteristics of eruption deposits. For example, a lahar derived from magmatic activity can contain fresh magmatic materials such as pumice, volcanic glass shards and mineral grains. Lahar deposits associated with magmatic pyroclastic eruptions include gravel and sand-sized particles (e.g., Pierson et al. 1996) and can vary from matrix-rich (not clay-rich but silt-and sand-rich) to clast-rich, depending on the eruption style, size, and explosiveness (e.g., Mothes and Vallance 2015). In contrast, phreatic and magmatic hydrothermal activity through extant hydrothermal systems can produce clay-rich ejecta and hence clay-rich, syn-and post-eruptive lahars Kataoka and Nagahashi 2019). Also, gravitational collapse of hydrothermally altered parts of the volcanic edifice can produce clay-rich lahars (Vallance and Scott 1997;Vallance 2005). Therefore, immediately after a phreatic/hydrothermal eruption, clay mineral assemblages and clay content (a cohesive lahar can be defined by the lahar deposit matrix containing > 3-5 wt% of clay: Scott et al. 1995) should be evaluated to predict the cohesiveness and its effect on the flow type of possible future lahars. The 2018 eruption deposits at Motoshirane contain ~ 10 wt% clay-sized particles, together with a mineral assemblage (Fig. 7) that suggests they were sourced from hydrothermally altered zones. The clay content in these primary ashfall deposits is lower than that in the deposits of the 2014 Ontake phreatic eruption (> 30 wt%), which produced cohesive lahars containing 10-20% matrix clay early in the eruption aftermath . The primary tephra characteristics of the 2018 phreatic eruption at Motoshirane suggested that an extremely cohesive flow following the eruption was unlikely. However, the presence of kaolin group minerals and clay content in the 2018 eruption deposits (Figs. 6 and 7) indicates that future lahars, derived from the tephra deposits, may be cohesive as such clay minerals can increase the yield strength of the flow (Hampton 1975;Pierson 2005). It is also predicted that the winnowing of fine particles from primary

Constraining snow and snowmelt volumes for snow-related lahar simulation
This study considered the SWE in near-vent areas to determine the values of input parameters for Titan2D flow simulations (Table 1). In Japan, no governmental agency performs either manual or automatic SWE measurements. At weather stations operated by the Japan Meteorological Agency, snow depth is automatically measured in snowy regions, in addition to other meteorological variables such as air temperature and precipitation. However, as in many other countries, most stations are installed at low altitudes (the Kusatsu weather station at 1223 m a.s.l. in Kusatsu town is the third high-altitude located station in Japan where snow depth is measured (Suzuki 2018), even though the area of new vents at Motoshirane is located at ~ 2050 m a.s.l.). No other reliable data on SWE and/ or snow depth are available, and thus, it is impossible to directly assess SWE and snow depth in the near-vent area. Snow-related reports from ski resorts located close to the volcanoes can be obtained; however, the information on snow depth is confined to middle altitudes and reflects the prevailing local conditions. In this study, using the snow survey data, we estimated the SWE in the near-vent area based on the relationship between SWE and altitude in areas of different elevations around Kusatsu-Shirane volcano (Fig. 8). It is widely known that a linear relationship between SWE and altitude can be observed in many mountain areas in Japan, including active volcanoes such as Daisetsu (Yamada et al. 1979), Ontake (Iyobe et al. 2016), and Adatara and Azuma (Matsumoto et al. 2019). These studies also highlight the fact that these relationships vary with season and slope orientation (to the prevailing wind direction in winter). Thus, it is important to understand the characteristics of the spatial distribution of SWE on individual snow-clad volcanoes for snow-related assessment. The SWE can also be measured during a background period, i.e., before volcanic unrest, and can be considered to constrain the values of input parameters.

Recommendations and limitations of the Titan2D flow simulation tool
Lahar flow models have been widely applied using a variety of methodologies (Manville et al. 2013). A popular empirical method such as LAHARZ (Iverson et al. 1998), used to predict the inundation area of lahar flow, was not suitable for assessment of the 2018 eruption at Kusatsu-Shirane (Motoshirane) because of the poor definition of the spatial distributions of previous lahar deposits that might have been considered for comparison. In contrast, application of the Titan2D (Pitman et al. 2003;Pitman and Le 2005) simulation tool depends on input parameters, such as the initial volume of flow, the solid fraction in flow, the internal friction angle, and the bed friction angle. These parameters can be easily estimated for the Motoshirane case through field survey and laboratory work. Furthermore, for application of the simulation tool, a normally specified desktop computer (Procter et al. 2012) is necessary, and the tool can be used to compute the results in a relatively short time. This enabled the authors to conduct numerous simulations under different conditions (such as a wide range of bed friction angles) and under several scenario settings. Therefore, the Titan2D tool can be widely applied in forecasting snow-related lahars, even without the availability of historical and geological data for calibration. The estimated flow velocities (Table 3) are higher than the reported or calculated values (e.g., Pierson 1986;Pierson and Costa 1987) for moderate to large lahars. Williams et al. (2008) reported velocities of 3.1-7.0 m/s for observed and estimated small-scale lahars, at Vulcán Tungurahua with volumes (54,000 m 3 ) of the same order as the hypothesized lahars for the 2018 eruption at Motoshirane (Table 1). These overestimated flow velocities are possibly attributable to the steep topography in the upper reach, which enhanced the simulated flow acceleration. Procter et al. (2012) also reported several difficulties in simulating flow travel times in proximal sections of a catchment, using the Titan2D model, owing to topography, although flow paths were appropriately defined by using the model. The single-phase and twophase models used in Titan2D lahar simulations indicate different results (Figs. 11,12 and 13). Simulated lahars using a two-phase model travelled further than those which used a single-phase model under conditions of the same bed friction angle (Fig. 11a, b). The two-phase flow simulations with lower bed friction angles occasionally depict numerical instabilities, indicating caution when assessing lahar hazards and risks even with plausible input values. In contrast, one-phase simulations showed stability in flows progressing through the valleys. Using lower values of bed friction in the single-phase model, as indicated in the present study, can aid the simulation of more fluidal and far-reaching flows. For wider applications, calibration of bed friction angles using past lahar deposits on other volcanoes will be helpful in constraining the models. The cohesiveness of the flow cannot be precisely expressed in the Titan2D model, although it is an important parameter for clay-rich lahar flow. As most previous studies on lahar flow simulations using the Titan2D flow model have reported the involvement of sandy and gravelly lahars (i.e., less cohesive) for hindcasting/forecasting (Williams et al. 2008;Córdoba et al. 2015;Procter et al. 2010Procter et al. , 2012, verification of the model for cohesive flows, such as that in the Motoshirane case and as determined by estimation of the clay content and mineralogy of the lahar source material, will be necessary in future studies.

The consequences of the 2018 eruption and suggestions for future crisis hazard assessment
A prompt hazard assessment, for the benefit of residents, communities, and emergency planners, is required in terms of information on flow paths and travel distances (inundation areas), rather than accurate simulation of physics within a model. Data on lahar hazard zonation using Titan2D flow models can be derived from the flow paths, areas, and distances of lahars under several scenarios using both single-phase and two-phase models (Figs. 11, 12 and 13). The information on the identification of key high-risk locations in the ski areas and near National Highway 292, as indicated by the simulations of lahar flows, was shared with the emergency planning board at the end of March 2018. At the time, it was recommended to close the gates on the highway when either further volcanic unrest or rapid changes in weather with heavy rain and/or snowmelt, were detected.
The monitoring system in Furiko-zawa Creek and the Yazawa-gawa River ( Fig. 2 and Appendix 1) did not detect any major lahar events (debris flows and hyperconcentrated flows) from early March 2018 to mid-June 2019, although stage and turbidity changes in the rivers were recorded when heavy rains occurred. The data obtained by the monitoring systems were shared among scientists and emergency planners in Kusatsu town (Fig. 3). An internet-based remote camera, which acquires images every 10 min (and videos when certain movements are detected) transmits data to an internet cloud server. This is useful for understanding river and snow conditions. Consequently, the absence of a lahar after the 2018 eruption was probably because of the formation of lower volumes of initial tephra deposits on the upper slopes (30,000-50,000 tons: The Joint Research Team for ash fall in Kusatsu-Shirane eruption 2018) and the occurrence of less torrential rains during the early stage of the posteruption period.
The findings of this study do not imply that Kusatsu-Shirane (Motoshirane) volcano and its river systems have not experienced lahars in the past. It is also important to study and find the past lahar deposits sourced from both Mts. Motoshirane and Shirane in the drainage system. Future lahars and resultant hazards at Motoshirane depend on the eruption magnitude, type, and interaction with snow and rain events. The crisis hazard maps for snow-related lahars created in this study (Figs. 11,12 and 13) will be helpful as a background hazard map for conducting hazard zonation for the assessment of future Motoshirane eruptions. The revised parameter values applicable to the tephra and snow conditions at the time of future eruptions, together with the bed friction angles considered in this study, will result in effective crisis hazard assessment for any eruption. Even if future eruptions occur from new vents, without the availability of historical and geological lahar data, the procedure used in the study (Fig. 3) will enable rapid lahar prediction and hazard zonation.

Conclusions
Several previous studies, based on the application of numerical simulations for lahar hazard zonation, have hindcast past lahar events or forecast future lahars by calibrating historical and geological lahar records (e.g., Williams et al. 2008) prior to the next volcanic unrest. In contrast, the present study demonstrated the use of a deterministic approach to forecast snow-related lahars during the 2018 eruption crisis that originated at unforeseen vents at Motoshirane. Even without the availability of historical and geological lahar records for the volcano, characterization of proximal tephra deposits, snow surveys, and numerical lahar flow simulations using the Titan2D model enabled the prediction of flow type and flow paths, and the conduction of hazard assessment and zonation in a relatively short timeframe.
The characterization of primary tephra deposited proximal to the vents indicated that a non-juvenile eruption originated from hydrothermally altered zones. The tephra deposits contain approximately 10 wt% clay and clay minerals, which contributed to the cohesiveness of the predicted lahar flows.
Several scenarios for numerical flow simulations were applied. This included snow and snowmelt being incorporated to predict the future lahars, with or without a new eruption. To determine the input parameters for the flow model, the consideration of SWE, immediately after the eruption and in the post-eruptive period, effectively helped to estimate snow volume at the high-altitude regions near the vent area.
The present study reported the application of not only the two-phase model of Titan2D, commonly used for lahar assessment, but also the single-phase Coulomb model, which is relatively more stable numerically. Aspects of the results, such as apparent overestimation of flow velocities, highlight the limitations of the models. Nevertheless, the flow paths (inundation areas) and travel distances, calculated in repetitive runs with variable scenario settings, are useful for lahar hazard assessment and for identification of potential high-risk areas intended to be monitored closely whenever further volcanic unrest or abrupt changes in weather conditions are detected. The procedures used in the present study provide insights for forecasting of lahars sourced from new vent eruptions, as well as for indicating an appropriate use of latest eruption-and snow-related information.