Natural Hazards and Earth System Sciences

Abstract. Karymskoye caldera lake is a nearly circular body of water with a diameter of approximately 4 km and a depth of up to 60 m. The sublacustrine, Surtseyan-type eruption in the lake on 2–3 January 1996 included a series of underwater explosions. A field survey conducted the following summer showed signs of tsunami wave runup around the entire coastline of the lake, with a maximum of 29 m runup at the north shore near the source of the eruption, and 2–5 m runup at locations on the east and south shore far away from the source. The tsunami has been simulated using the numerical long wave model COULWAVE, with input from reconstructed realistic pre-eruption bathymetry. The tsunami source was chosen as suggested by Le Mehaute (1971) and Mirchina and Pelinovsky (1988). The initial wave was prescribed by a parabolic shape depression with a radius of R=200 m, and a height of 23 m at the rim of the parabola. Simulations were conducted to show principle directions for wave propagation, wave speed and arrival time for the leading wave group at the shore, and the distribution of wave height throughout the lake. Estimated result for wave runup are of the same order of magnitude as field measurements, except near the source of the eruption and at a few locations where analysis show significant wave breaking.


Introduction
About 25% of all the fatalities directly attributable to volcanoes during the last 250 years have been caused by volcanic tsunamis (Begét, 2000), both by eruptive processes (e.g.1883 Krakatau, 1937 Rabaul), or by collapses of volcanic cones (e.g.Iliwerung 1979, Paluweh 1928, Mount Mayuyama 1792).Volcanic explosions caused 25% of the volcanic tsunamis reported for the four last centuries (Latter, 1981).They usually deposit small volumes of material (< 1 km 3 ) compared to ignimbrites, caldera subsidences or flank collapses (1-100 km 3 ) and travel shorter distances (< 50 km), but may generate deadly local tsunamis, especially in relatively small, shallow bodies of water like bays or lakes.
Such events are known from case studies in Kamchatka (Belousov and Belousova, 2001;Belousov et al., 2000) and in Nicaragua (Freundt et al., 2007), and the possibility of producing an effective warning system has been discussed (Bellotti et al., 2009).The 2-3 January 1996 sublacustrine eruption and subsequent tsunamis that occurred in Karymskoye lake (Kamchatka, Russia) represents a unique opportunity to model tsunamis generated by volcanic explosions.
Karymskoye lake is a nearly circular, bowl-like water body, approximately 4 km in diameter, having steep flanks and a nearly flat bottom with depths of 50-62 m (Figs. 1  and 2).A new, partly submerged crater, with an inner diameter of 500-550 m, was formed by the 1996 eruption near the north shore of the lake (Fig. 3).Belousov and Belousova (2001) described one of the strongest observed explosions as follows: a rapidly rising (∼ 3 s), smooth-surfaced bulbous mass of gas and pyroclasts expands within a shell of water up to 450 m high and then becomes unstable and is pierced by cypressoid jets of pyroclasts with gas and water.Initial velocity of these jets is estimated as 110 m s −1 .Simultaneously, a bore-like elevation of the water surface starts to propagate radially at about 20-40 m s −1 .Within ∼ 15 s from the beginning of the explosion, the jets reach their maximum height of about 1 km and then collapse back toward the lake to produce a   base surge travelling across the lake surface at ∼ 12.5 m s −1 (600 m in 48 s).Other weaker explosions produced a single, near-vertical, column-like jet about 100-150 m high.
The column then collapsed back into the lake and generated a subtle base surge.
Underwater explosions occurred every 4-12 min, with an average interval of 6 min.Considering that the eruption lasted 10-20 h, 100-200 explosions may have occurred in the lake during the eruption (Belousov and Belousova, 2001).Tsunamis periodically affected the shores and forced water from the lake into the only outlet, the Karymskaya River, thus forming lahars.The largest of the tsunamis occurred at the end of the eruption.This may indicate that the explosions that took place near the end of the eruption were stronger than those of the main stage.Belousov and Belousova (2001) estimated that ∼ 2 × 10 5 m 3 of tephra per explosion were emitted, corresponding to kinetic energy of ∼ 2 × 10 12 J.

Review of theoretical background
An exact mathematical description of explosion-generated water waves is not available at the moment, and would in any case be difficult to incorporate as initial condition for a long wave, shallow water model.Le Mehaute (1971) suggested an empirical formula for waves at some distance from the eruption, with initial water displacement η and water velocity u given as respectively (see also Le Méhauté and Wang, 1996).The formula for η includes the maximum surface elevation η 0 , the characteristic radius R of the initial surface deformation, and the radial coordinate r giving the distance from the source point.This formula was found to give the best agreement between linear wave theory and experiments with low energy (E ∼ 2 × 10 6 − 3 × 10 10 J) under-water explosions (Mirchina and Pelinovsky, 1988).
The formula used by Le Mehaute (1971) is not a suitable initial condition for long wave models, because it prescribes a discontinuity of the water level at r = R.We have therefore modified the formula for the free surface displacement where the factor P r decides the rate of decline in surface elevation for r > R. In the simulations for small amplitude initial disturbance we have used P r = 200.0,which results in a profile that is numerically equivalent to the Le Mehaute profile.
Le Mehaute (1971) estimated the maximum wave length, group velocity, and amplitude in shallow water at a distance r from the eruption center as where g = 9.81 m s −2 is the acceleration due to gravity The relationship between the maximum initial surface displacement (in meters) and eruption energy (in Joule) can be estimated as Based on predictions by Belousov and Belousova (2001) that the kinetic energy of the largest explosions were of the order E ≈ 2 × 10 12 J, the corresponding maximum initial surface displacement according to Le Mehaute's formula is η 0 = 23.0 m.The characteristic length scale R = 200.0m is determined by the radius of the caldera created by the eruptions.

Field work
A survey of Karymskoye lake was conducted in September 1996.Wave runup was recorded at 23 locations around the lake, as shown in Fig. 2 (note that TS12 is absent from the records).The runup represents the difference in elevation between the highest lake level during the eruption and the level reached by the inundating wave (which was marked by various material deposited by tsunami).Details can be found in Belousov and Belousova (2001).The runup data has mean value and standard deviation of 7.3 ± 5.9 m, with a maximum of 29 m runup at the north shore near the source of the eruption, and 2-5 m runup at locations on the east and south shore far away from the source.

Reconstruction of pre-eruption model area
Topographic maps of the area surrounding the lake are available, both for the pre-eruption and post-eruption states, but a pre-eruption bathymetry is not available.The digital elevation model of the Karymskoye lake bathymetry and shore region was constructed on a 274 × 334 grid with grid resolution x = y = 21.29 m.A pre-eruption bathymetry was reconstructed based on present day data, combining topographic data from the ASTER satellite with bathymetry data obtained from a survey with an echo-sounder device ("Eagle Seacharter" with integrated GPS).Topographic changes associated with the eruption deposit were removed manually, with guidance from pre-and post-eruption aerial photographs (see Fig. 3).The pre-eruption bathymetry at the location of the new crater was estimated by extrapolating bathymetry data from the surrounding area into the region where the eruption deposits had been removed.

Numerical simulation
Numerical simulations have been performed using the COULWAVE model, which has been developed to simulate the evolution of fully nonlinear (wave amplitude to water depth ratio of a/ h ∼ 1) long waves over a variable bathymetry (Lynett and Liu, 2002;Lynett et al., 2002).The model is based on a set of weakly dispersive Boussinesqtype equations, which is a suitable framework for studying the propagation and transformation of nonlinear long waves on fairly large computational domains.This particular model has previously been used in studies of tsunami waves from various sources (Lynett and Liu, 2002;Korycansky and Lynett, 2007;Geist et al., 2009), as well as for smaller long wave phenomena such as ship waves (Torsvik et al., 2009a,b).The original model has been modified by the authors to include the initial conditions for an under-water explosion.
T. Torsvik et al.: Numerical simulation of tsunami from volcanic eruption Simulations were performed on a x = y = 10 m grid, where data for the bathymetry was further refined by interpolation.The total simulation time was 5 min, which is sufficient time for the leading waves to propagate across the lake.A flow relaxation scheme was used on the boundaries of the numerical domain to absorb any outgoing waves, but this was of little practical use due to the absence of open boundaries.Waves reaching the coastline were reflected back into the lake.
The numerical model became unstable for large values of η 0 , and we were not able to run the model with η 0 = 23 m.Simulations have been performed with initial wave heights of η 0 = 0.23 m, 0.46 m, and 2.3 m.For the two first cases, the leading wave group is expected to behave like linear waves, whereas weakly nonlinear effects can be expected in the last case.The equivalence source function Eq. ( 1) and formulas (3) for maximum wave length, amplitude and velocity are based on linear wave theory.
The model of Le Mehaute as well as the COULWAVE model can not be used for description of the real, strongly nonlinear wave field at the center of the eruption, but the numerical and theoretical results are expected to converge towards the realistic wave field with increasing distance from the eruption center.That is why we have applied an equivalence source function that is good only for wave description in the far field.In the case of Karymskoe lake, where the model area is bounded within a fairly small area, it is difficult to make a clear distinction between the near field and far field region.Therefore we also need to include the near field region in the simulations.In order to perform the calculations without inducing too large numerical errors, we first decrease the real wave amplitude by 100 times and run the simulations for almost linear waves of small amplitude, and thereafter the same factor of 100 is multiplied back into the wave amplitude before the estimated runup heights are calculated.
The procedure described above has obviously some limitation with respect to accuracy.We test the accuracy of the procedure in two ways.Firstly, we make a comparative study between runup measurements and predictions based on the simulations and a theoretical runup formula.The results can be expected to deviate greatly for coastal regions close to the eruption center, but the deviations should decrease with increasing distance from the eruption center.Secondly, we make a comparative study between simulated results using different reduction factors (10, 50, and 100).The purpose of this study is to ensure that the simulated results are not highly sensitive to the scaling factor.

Linear waves
A simulation with η 0 = 0.23 m was used for the linear wave case.Figure 4 shows the maximum water elevation for each grid point throughout the lake.Only the maximum of the 3 leading waves have been included in the record, in order to avoid noise due to waves reflected from the shore, and in order not to take into account short crested waves that are not accurately represented by the model equations.Note that color scale in Fig. 4 is logarithmic.A circular feature indicates the extension of the initial profile, and the maximum at the center of the circle is produced when the initial profile collapses.The wave amplitude is reduced rapidly as the waves propagate away from the source area, but increase again near shore due to wave shoaling.However, the runup height is only of the same order of magnitude or less than the wave amplitude at the source, and there is generally no wave amplification in the basin.Wave amplitude near the shore is shown in more detail in Fig. 5, which shows maximum water elevation at the 5 m depth contour along the shore line.The strongest wave impact is seen on the north and west shores.
Figure 2 shows the location of the 91 numerical wave gauges used in the simulations, where the wave amplitude is recorded at each time step throughout the simulation.Some gauges are located along lines radiating out from the location of the eruption, while others are placed offshore measurement sites for wave runup.The panels in Fig. 6ae show the depth along different lines, maximum wave amplitude of the 3 leading waves at each wave gauge, the maximum wave amplitude predicted by Le Mehaute's theory Eqs.(3), and the mean wave speed between wave gauges.The wave speed was estimated by detecting the passing of the first wave trough at each wave gauge, as this measure is less influenced by wave-wave interaction than the wave peaks.The comparison between theory and simulations seem to be satisfactory along all radial directions, except in the case of westward wave propagation (Fig. 6a).
Wave propagation along the westward line in Fig. 6a differs from the other records since this line is located in shallow water close to the north shore.Note that the scale for the wave amplitude in Fig. 6a is larger than for Fig. 6be, because the first point is located nearly at the edge of the initial source.Waves propagating westward are highly influenced by wave refraction due to the close proximity to the north shore, and it is therefore not expected that the leading wave group propagates along the straight line formed by the wave gauges.This is also indicated by Fig. 6f, which show the Froude numbers along each line of propagation, calculated from the estimated wave speeds V between wave gauges, the acceleration of gravity g, and the local depth h.The large values in Froude number for westward propagating waves indicate that the leading wave group does not follow the line of wave gauges, but instead originate from further offshore.The large jumps at the tails of the other Froude number records may also indicate that the wave trough measure for wave speed is inappropriate in shallow water when waves are shoaling.
Comparison between the theoretical prediction of wave amplitude and the simulated result shows that there is a similar trend in amplitude reduction with distance, but there is only a good quantitative match for waves propagating to the east (Fig. 6e).The results do not match near the eruption source and in shallow water where waves are shoaling, but this is expected since the theoretical prediction assume that waves have propagated over a long distance in water of constant depth.However, the simulated and theoretical results in Fig. 6b-d are converging in the deep part of the lake as the distance increases, indicating that a better match between theory and simulations could be obtained in a larger water basin.
The wave speed of the leading wave group is up to 24 m s −1 in the deepest part of the lake, and the first wave trough reaches the south-east coast approximately 140 s after the eruption.After an initial acceleration, the waves are seen to propagate at Fr h = 0.95 − 0.97 in the deep part of the lake.By Le Mehaute's theory we should have λ max ≈ 1.7R = 340 m, and the Froude number for linear waves in water of depth h = 60.0 m can be estimated as which is significantly slower than the simulated results.A possible reason for the deviation between the theoretical and simulated results may be that the source length scale R is not very large compared to the water depth h, and therefore linear shallow water theory may not be well suited for this problem.The difference in initial wave amplitude does not have any significant influence on characteristic features of the leading wave group.Both the travel time and the group structure is preserved, but the wave amplitude for η 0 = 2.3 m deviates slightly from the two other results.Although the nonlinear effects become larger with increasing wave amplitude, nonlinear effects are likely to remain weak in the long wave part of the wave spectrum, including specifically the waves in the leading wave group.Wave breaking is also likely to be an important factor, both near the wave generating source and in the near shore region.However, the source function used in the simulations is applicable to linear waves far from the eruption center, and consequently also far from the region near the eruption center where wave transformation due to wave breaking is most likely to occur.Since the source function is partly empirical, this means that the effect of wave breaking is already parameterized within the source function formulation.As a consequence, wave predictions near the eruption site, where wave transformation due to nonlinear effects is more pronounced, are not as reliable as results for the far field.

Simulated wave amplitude, modeled and observed wave runup
Table 1 includes the observed runup heights at different locations around the coast of the lake, and data for the first wave arriving at numerical wave gauges at offshore locations close to where runup was measured.No suitable offshore point was found for TS19, since the waves most likely propagated along the shore to this point.A source function with initial maximum wave height of η 0 = 0.23 m was used in the simulations.A comparison between measured runup data and simulated maximum offshore wave amplitude shows that the former is larger than the latter by a factor of approximately 300.An amplification factor of 100 is expected by the use of an equally reduced initial wave height at the source.These results indicate that there is an additional amplification factor of approximately 3, which can be considered as a typical factor for the ratio between wave runup and offshore wave amplitude.In the following analysis, an estimate of wave runup based on the simulated results is obtained by using the theoretical formulas for wave runup.Results at locations of the largest wave runup (TS17, TS18, and TS22) deviate from the rest of the data points.These locations are on the north-west and north part of the shore, near the source of the eruption.Highly nonlinear waves moving at an oblique angle along a steep coastline may attain unusually large amplitudes at the point of reflection at the shore.This effect has been studied extensively for solitary waves, where it has been found that Mach stem reflection, due to nonlinear interaction between the incident and reflected wave, can lead to amplitude amplification of up to 4 times the incident wave amplitude (Miles, 1977;Peterson et al., 2003).A similar type of process may occur at the coast near the eruption center, which would account for the large runup heights found at the locations TS17, TS18, TS19 and TS22.
A rough estimate for wave runup can be calculated using analytical formulas derived from nonlinear shallow water theory for sine wave runup on a beach of constant slope (Pelinovsky, 1982;Pelinovsky and Mazova, 1992;  locations, three are in the immediate vicinity of the source (TS17, TS18, and TS22), where the breaking most likely occurs directly at the shore.Wave breaking is likely to be a significant factor for the remaining four locations (TS1-TS3, and TS24), which may explain the large discrepancies between measured and modeled results at these points.Figure 10 show beach profiles included in the runup analysis, indicating both measured and modeled runup heights, for locations where the analysis indicate that wave breaking was not a significant factor.

Conclusions
Waves generated by the sublacustrine, Surtseyan-type eruption in Karymskoye caldera lake on 2-3 January 1996 have been examined by numerical simulations using a long wave Boussinesq-type model.The simulated maximum wave amplitude distribution for the leading waves (Fig. 4) is generally consistent with measured runup data (Fig. 10) along the coast of the lake.Due to numerical stability considerations, simulations were performed with reduced wave amplitude (η 0 = 0.23 m), making the simulated waves essentially linear.The wave speed of the leading wave group is limited by the local water depth, and reach a maximum of 24 m s −1 in the deepest part of the lake.Important features such as the wave speed, period, and group structure of the leading waves are not significantly influenced by the amplitude scaling, as demonstrated by the comparison between η 0 = 0.23 m, η 0 = 0.46 m, and η 0 = 2.3 m results.Both simulated and theoretical results predict wave reduction as waves propagate away from the eruption center, but the simulated waves are generally larger than the theoretical predictions, except for eastward propagation where there is good agreement between theoretical and simulated results.However, the theoretical and simulated results seem to converge also along the other radial directions, except for westward propagation where the waves propagate along the coastline.The theoretical predictions are expected to be valid at a large distance from the wave source, and this situation is not achievable due to the limited extent of the lake.Moreover, the theoretical formula was found to give the best agreement between linear wave theory and experiments for low energy explosions (E ∼ 2 × 10 6 − 3 × 10 10 J), whereas the kinetic energy of the simulated volcanic explosion is estimated to be of the order E ∼ 2 × 10 12 J, which is possibly beyond the applicable range of the theory.It has also been indicated by Le Méhauté and Wang (1996) that the theoretical formulas need to be modified in shallow water, but these modified formulas are not available in any publications known to the authors.Considering these limitations, it is not surprising that there are some discrepancies between the theoretical and simulated results.
Comparison between runup measurements and simulated wave amplitude offshore measurement sites show a clear correlation for locations far from the eruption center.The present model setup does not capture nonlinear effects believed to be particularly important in the coastal region near the eruption site.Estimated results for wave runup are of the same order of magnitude as measurements, except near the source of the eruption and at a few locations where the analysis shows significant wave breaking.
Considering all these parameters, such simulations are suitable for assessing tsunami hazards related to volcanic explosions.Many lakes and sea bays over the world are of volcanic origin or partly occupied by recent eruptive vents, and their shores are often densely populated.Surprisingly, tsunamis are commonly neglected in volcanic hazards scenarii and the present contributions may help in introducing tsunamis in future studies on volcanic hazard assessment in lacustrine and sea bay environments.

Fig. 1 .
Fig. 1.Aerial photograph of Karymskoye caldera lake, taken from the North.The horse-shoe shaped structure in the northern part of the lake is a partly submerged crater formed by the 1996 eruption.Photograph by A. Belousov, September 1996.

Fig. 2 .
Fig. 2. Bathymetry of Karymskoye lake, and topographic contours of the lake shore each 5 m.Markers indicate location of runup measurement sites (cyan, yellow, magenta, black) and of 91 numerical wave gauges (green).

Fig. 3 .
Fig. 3. Vertical aerial photographs of the 1996 eruption site in the northern part of Karymskoye lake.The upper photo shows the preeruption coastline (1984), and the lower photo shows the coastline after the eruption (August 1996).

Fig. 4 .
Fig. 4. Simulated maximum wave amplitude in Karymskoye lake during the 1996 eruption.The bathymetry is represented by dashed lines, which are drawn at 10 m depth intervals.

Fig. 3 .
Fig. 3. Vertical aerial photographs of the 1996 eruption site in the northern part of Karymskoye lake.The upper photo shows the preeruption coastline (1984), and the lower photo shows the coastline after the eruption (August 1996).Images courtesy of Institute of Volcanology and Seismology, Russia.

Fig. 4 .
Fig. 4. Simulated maximum wave amplitude in Karymskoye lake during the 1996 eruption.The bathymetry is represented by dashed lines, which are drawn at 10 m depth intervals.

Fig. 4 .
Fig. 4. Simulated maximum wave amplitude in Karymskoye lake during the 1996 eruption.The bathymetry is represented by dashed lines, which are drawn at 10 m depth intervals.

Fig. 5 .
Fig. 5. Maximum wave amplitude at the 5 m isobath along the coast of Karymskoye lake.The location of the center of the source (crater) is indicated by a black, dashed line.

Fig. 7 .
Fig. 7. Comparison between time series of surface displacement for different values of η 0 .Wave gauge location: (a) WGA and (b) WGB.

Fig. 9 .
Fig.9.Measured and modeled runup, and wave breaking analysis.Points where wave breaking is significant have the largest discrepancy between the measured and modeled wave runup (see Fig.8).