## Abstract

We investigate the ability of a submarine landslide to generate the tsunami waves observed on the Bulgarian coast of Black Sea on 7 May 2007. In our simulations, a landslide is presented as a quasi-deformable body moving along a curvilinear slope under action of the forces of gravity, buoyancy, water resistance and bottom friction. We employ the fully non-linear weakly dispersive model for tsunami wave simulations. The computations show that the initial landslide position on the real slope is extremely important for its dynamics and the wave generation process. We constructed some model landslides which generated similar waves to those observed. Moreover, these landslides stopped in the same region. Finally, we evaluated the significance of the frequency dispersion effects in the simulations.

Tsunami waves were observed at several points along the Bulgarian Black Sea coast on 7 May 2007. Ranguelov *et al.* (2008) pointed out the following facts about this tsunami:

‘no significant earthquakes in the region were recorded on this same day’;

‘the maximum sea level rise and lowering were respectively +1.2 m and −2.0 m’;

‘the chief period of the oscillations was between 4 and 8 minutes at most places’;

‘turbulence, strong water currents, mud waters and foam in some sites (e.g. in Balchik and Kavarna) were observed’.

Due to the absence of submarine earthquakes, Ranguelov *et al.* (2008) investigated the hypothesis of a submarine-landslide origin, using numerical simulations of the event. Four initial positions of the landslide were considered, located at a depth of 100 m. The volumes of the landslide were assumed as between 30 and 60 million m^{3}, with thickness ranging over 20–40 m. The model landslide was presented as a set of constant-volume blocks (Tinti *et al.* 1997). In the computations of Ranguelov *et al.*, all landslides stopped at a depth of 1000 m with a runout of *c.* 20 km. Velocity peak was *c.* 20 m s^{−1} after 200–300 s from motion initiation. The tsunami wave generation was simulated using the shallow-water approximation of the Navier–Stokes equations. In all cases, a two-beam propagation pattern of the computed tsunami was observed. Comparison of the obtained numerical solution with observational data showed satisfactory agreement at some points, but the highest waves were computed for the southern coast of Bulgaria where no reports of the tsunami were registered. Ranguelov *et al.* (2008) therefore showed that a submarine landslide could have generated waves compatible with the observations, but this hypothesis requires further research.

Later, Vilibic *et al.* (2010) investigated the hypothesis of a meteotsunami. The numerical simulations showed that sudden variations in air pressure, which were captured by some barograms, could have generated tsunami waves similar to those observed. However, the obtained distribution of the maximum waves along the coast was not in a good agreement with the observational data (again, some of the highest computed waves were in Burgas Bay, which is not consistent with the field data).

In the present paper, we continue to develop the hypothesis of the submarine-landslide origin of the observed tsunami. Namely, we focus on the consideration of numerous initial positions of the model landslide at a depth of up to 1500 m. Furthermore, the landslide is presented in our simulations as a quasi-deformable body moving along the real slope (Beizel *et al.* 2012). The wave propagation is computed with a fully non-linear weakly dispersive model (Fedotova & Khakimzyanov 2010). The significance of the frequency dispersion in the computations is discussed. Our simulations show that some landslides could have generated waves which agree perfectly with the observations. In contrast to Ranguelov *et al.* (2008), these model landslides start from depths of >800 m and do not generate high waves on the southern coast.

## Tsunami model

For wave propagation modelling we use the fully non-linear dispersive (FNLD) model on a rotating sphere (Fedotova & Khakimzyanov 2010). The model was derived without making any assumptions about the smallness of wave amplitude, and all the non-linear terms associated with dispersion are stored; it can therefore be used for calculation of the surface wave propagation over an uneven bottom both in deep water and in a coastal zone. In the equations of this model, the bathymetry is presented as a time-dependent function; the FNLD model therefore also allows to simulate the tsunami wave generation by submarine landslides, which extends the range of problems that can be solved within the framework of the known dispersive models on a sphere (Dao & Tkalich 2007; Lovholt *et al.* 2008; Kirby *et al.* 2013). It should be noted that these models can be used in simulations of landslide-generated tsunamis if the other models for the generation stage are employed, such as SAGE (SAIC's Adaptive Grid Eulerian hydrocode) in the study of Lovholt *et al.* (2008) and NHWAVE (Non-Hydrostatic Wave model) (Ma *et al.* 2012) in a more recent study of Schnyder *et al.* (2016). We use only FNLD model for all stages of modelling in this paper.

Moreover, the FNLD model can be written in a quasi-conservative form of mass and momentum balances (Fedotova & Khakimzyanov 2014). It also has a balance equation for the total energy in agreement with the similar equation for the three-dimensional (3D) model, which not only confirms the physical consistency of the model but also allows calculations to be monitored.

Instead of a direct approximation of the derivatives in the equations, we split the original system into a scalar equation of elliptic type and a set of hyperbolic equations (Gusev & Khakimzyanov 2015; Gusev & Beisel 2016; Khakimzyanov *et al.* 2018). The latter contains the first-order derivatives only, while the former is the uniformly elliptic equation for the dispersive component of the integrated pressure, which does not contain time derivatives of the main functions.

Eliminating dispersive components, the non-dispersive non-linear shallow-water (NLSW) model on a rotating sphere is obtained.

We construct the numerical algorithm for the extended system by alternately solving the hyperbolic and the elliptic sub-problems at each time step. This approach allows the use of the most suitable methods for each sub-problem. The hyperbolic system is addressed using a second-order predictor-corrector finite-difference scheme, and the elliptic equation is solved by integro-interpolation and successive over-relaxation methods. For the full solution, we develop the algorithm of second-order approximation in space and time. For a more detailed description of the algorithm and its properties (such as correctness, stability, numerical dispersion and dissipation), see Gusev & Khakimzyanov (2015) and Khakimzyanov *et al.* (2018).

## Landslide model

We use a simplified landslide model due to the lack of information about its typical physical characteristics. In our investigation, the landslide is presented as a quasi-deformable body. Beizel *et al.* (2012) and Dutykh & Kalisch (2013) described this model in detail, so we only provide the main aspects here.

We first assume that the landslide motion is local in the Earth scale, and employ Cartesian coordinates for its calculation using a trivial transformation from the spherical coordinates:
*x*, *y*) and (*λ*, *φ*) are the Cartesian and spherical (geographical) coordinates, respectively; *R* is the Earth's radius; and the computational domain in geographical coordinates is determined as [*λ*_{min}, *λ*_{max}] × [*φ*_{min}, *φ*_{max}].

We then define the bottom function *z* = *−h*(*x*, *y*, *t*) as a sum of the static component *z* = *h*_{bt}(*x*, *y*) and the landslide component *z* = *h*_{sl}(*x*, *y*, *t*). Here and below, *t* denotes time. At each time moment, the position of the landslide is determined by the point *D* of the function *h*_{sl}) moving along the bottom according to the law of non-free motion of a material point. At the initial time *t* = 0 the landslide is at rest, the function x = *D*_{0} specifies its initial form, and the point *x*_{c}(0) = *y*_{c}(0) = *z*_{c}(0) = *t* > 0 the landslide surface is given by the function

The landslide surface will therefore be deformed according to irregularities of the slope, and on a flat slope it will move like a rigid body. Separate points of a modelled landslide, having different horizontal coordinates, are allowed to move relative to one another in the vertical direction only. The landslide thickness is constant during its motion, which may lead to some overestimation of wave heights generated by the landslide in shallow water (as compared to ‘fully’ deformable landslides). However, in many studies of landslide-generated tsunamis (e.g. Grilli & Watts 1999, 2005; Enet & Grilli 2007; Whittaker *et al.* 2015; Gylfadottir *et al.* 2017) the landslide was presented as a rigid body. Moreover, Watts & Grilli (2003, p. 370) noted that the deformation of the landslide decreases its thickness but increases its acceleration, and ‘it is not clear how tsunami amplitude would be affected by deformation’. The approach used in the present paper allows the landslide motion to be simulated on a real bathymetry and the numerical algorithm to be verified (see Gusev *et al.* 2013; Gusev 2014; Khakimzyanov *et al.* 2018) by comparing the obtained solutions with the known experimental data, where the landslide was modelled as a rigid body on a flat slope.

The forces determining the landslide motion are the forces of gravity, buoyancy, friction on the irregular bottom and water resistance (Beizel *et al.* 2012; Dutykh & Kalisch 2013). In the equations of landslide motion the following parameters are significant: the landslide volume *V*; the coefficients of associated mass *C*_{w}, friction *C*_{fr} and hydrodynamic resistance *C*_{d}; *γ* = *ρ*_{sl}/*ρ*_{w} > 1, where *ρ*_{w} is the water density and *ρ*_{sl} is the density of landslide mass; and *C*_{fr} = tan *θ**, where *θ** is the angle of friction between the landslide and the bottom.

In our 2D calculations, we use the following formula for the landslide shape:
*D* = [*x*_{c} *−* *b _{x}*/2,

*x*

_{c}+

*b*/2] × [

_{x}*y*

_{c}

*−*

*b*/2,

_{y}*y*

_{c}+

*b*/2],

_{y}*T*is landslide thickness and

*b*and

_{x}*b*are its length in the

_{y}*x*and

*y*directions, respectively. Along these directions, the landslide shape is a cosine arc. In contrast to a Gaussian shape which was used, for example in the studies of Grilli & Watts (2005) and Enet & Grilli (2007), the shape specified by Equation (1) has the finite support

*D*. The semi-elliptic shape used by Grilli & Watts (1999) and Whittaker

*et al.*(2015) causes singularities in simulations on the bounds of its support, while the shape used here does not. Finally, the volume of this landslide can be easily calculated by the formula

Instead of Equation (1) another landslide profile *h*_{sl} can be used, but it should be smooth enough for the computations with non-linear dispersive models.

We emphasize here the connection between the landslide motion and the wave generation process in the model. Once again, the equations of the FNLD model contain both the time and space derivatives of the bathymetry function (Gusev & Khakimzyanov 2015; Khakimzyanov *et al.* 2015, 2018; Gusev & Beisel 2016); any transformation of this function produces the water flow. The numerical algorithm of landslide-tsunami generation is the following. At the initial time, the water is at rest. At each time step, we first compute the values of the bathymetry function through the landslide law of motion. Using these new values, we then compute the wave generation process with the predictor-corrector scheme.

## Simulation concept

Following the study of Ranguelov *et al.* (2008), we compare our computed results with observations made at seven coastal points of the Bulgarian Black Sea (Table 1).

Our previous study (Khakimzyanov *et al.* 2015) showed that the initial position of the landslide on a real slope is very important for its dynamics and the distribution of the maximum wave heights on the coast. For this reason, we focus the present study on searching for initial positions from which a landslide could have generated tsunami waves similar to those observed.

The concept is the following. First, we fix a basic set of the landslide parameters and perform the computation for a particular fixed initial position. In this landslide model, the variation of the thickness of the landslide does not affect its dynamics and trajectory, but it affects the computed heights of the associated tsunami waves. We therefore choose the thickness in such a way that computed wave heights fit best with the field data. To estimate the quality of the fit, we use the sum of squared deviations of the computed maximum and minimum (negative) water levels from the observed data. It should be noted that we do not need to perform the computations for all considered landslide positions with some range of the thickness to choose the ‘best’ thickness (it would take too much CPU time). Instead, we assume here that the computed wave heights will increase (decrease) proportionally along the coast with the variation of the thickness. We therefore perform only one computation for each initial position and normalize the obtained wave heights, dividing these heights on a characteristic height of the computation (e.g. an average positive height of the maximum waves computed in the considered coastal points). Note that the value of this characteristic height will be specific for each computation. The observed heights are also normalized using the characteristic height of the observations.

We then calculate the fit of the normalized heights in each computation to the normalized heights of the observations. This fit indicates the potential of a landslide to generate the observed tsunami, starting its motion from this position. Finally, the landslides on the ‘best’ initial positions are simulated with some variation of the thicknesses to obtain the best fit. We noted that the ratio of the ‘best’ thickness to the basic thickness is a bit less than the ratio of the characteristic wave height in the observation to the characteristic wave height in the basic computation.

The fixed parameters in the landslide law of motion are:
*et al.* (2015), which allowed the trajectory of the model landslide to fit well with the documented landslide (Kazantsev & Kruglyakov 1998). The form of the landslide was set according to Equation (1), with *b _{x}* =

*b*= 5 km,

_{y}*T*= 40 m and

*V*=

As mentioned above, for every initial position the thickness (and consequently the volume) was varied by up to 400 m to reach the best fit with the observations; the maximum volume of the model landslide was therefore 2.5 km^{3}. We have no information on the presence of slide masses near the Bulgarian coast of Black Sea, but the submarine landslide of 40 km^{3} was documented (Kazantsev & Kruglyakov 1998) near the Russian coast (see the study of Khakimzyanov *et al.* 2015 for details and simulation results). The considered volumes do not seem very huge. Previously, we investigated another horizontal extent of the landslide (*b _{x}* =

*b*= 2500 m), but the ‘best’ thicknesses in these simulations were about four times larger than with

_{y}*b*=

_{x}*b*= 5 km, which seems less realistic. The landslide lengths

_{y}*b*and

_{x}*b*can be different from each other, but these cases are not investigated in the present study.

_{y}In the numerical simulations, we used a geographic grid with the angular resolution of 3.75′′. Due to this resolution, the landslide profile is described by 61 × 44 grid nodes. The bathymetric data were obtained by applying bilinear interpolation to data retrieved from the GEBCO One Minute Grid (published in 2008). The main fragment of the computational domain is presented in Figure 1. The simulation time was three hours, the grid size was 2881 × 1921 nodes and the time step was about 0.6 s, which led to a CPU time for one single case of about 31 hours (with the FNLD model). Note that in the computations the shoreline is approximated by a polygon with sides parallel to (spherical) coordinate axes. We did not employ run-up algorithms and mounted a vertical wall near the shoreline for simplification. This simplification may decrease the computed wave heights on the coast, but we believe that it does not seriously affect the pattern of the distribution of the maximum and minimum wave heights recorded during the whole simulation time.

Firstly, we chose 40 different initial locations along the Bulgarian coastline, located mainly at depths of 200, 1000 and 1500 m. These locations are depicted with black rectangles in Figure 1. We could then delimit the area where the hypothetical landslide could have taken place. In a second iteration, this area was refined with an additional 171 initial landslide locations marked with crosses in Figure 1. The simulations allowed us to select six trajectories (*L*_{1}–*L*_{6}) and thicknesses for landslides which best fit the observations. The initial positions, thicknesses and initial depths of these landslides are listed in Table 2.

The landslides starting at the other positions or with different thicknesses either had a trajectory similar to *L*_{1}–*L*_{6} landslides or generated tsunami waves with a poor fit to the field data.

## Dynamics of the landslides

To show the dynamics of landslides *L*_{1}–*L*_{6} we present the graphs of absolute velocities of the respective mass centre. In the landslide model, these velocities are equal to the absolute velocities of the point *L*_{2} and *L*_{3} stopped abruptly at the end of their trajectories. Generally, the maximum velocities are somewhat larger than those in the simulations of Ranguelov *et al.* (2008) (*c.* 20 m s^{−1}). In contrast to that study, all of our model landslides did not stop at a depth of 1000 m, but rather at a depth of 1850 m. The interesting result is that landslides *L*_{1}–*L*_{6} stopped at the same region, which can be explained by the fact that the trajectories of landslides *L*_{1}–*L*_{6} lie in an amphitheatre-like area on the slope. By reducing the coefficient of friction, it is possible to advance the model landsides a little further. However, we can recommend this region of convergence for underwater explorations.

The maximum of local Froude numbers Fr = *et al.* (2015). However, it should be noted that in our computations the absolute acceleration ^{−2} (mostly when an abrupt stop occurred), while Whittaker *et al.* (2015) considered accelerations of 1 and 1.5 m s^{−2}. In a future paper, we will present the comparisons of our computations with the Whittaker *et al.* (2015) experiments and discuss the capability of different hydrodynamic models to reproduce the measured tsunami waves.

## Comparison with field data

Comparisons of the computed minimum and maximum wave heights with the observed data (Fig. 3) demonstrate a best fit with landslide *L*_{1}. In this case, there is a perfect agreement of the computed wave heights with those observed for all points.

For the considered landslides, the distribution of the maximum wave heights is quite different. Similar to the study of Khakimzyanov *et al.* (2015), we see the great importance of initial position of the landslide. In contrast to Ranguelov *et al.* (2008), the distribution only sometimes had a two-beam pattern. In these cases, the maximum and minimum wave heights are observed only at two coastal regions placed some distances from the closest coastal point from the landslide initial position. In simulations with the ‘best’ landslides *L*_{1}–*L*_{6}, this pattern is not observed. For example, Figure 4 shows the maximum height distribution for the landslides *L*_{1} and *L*_{5}. In general, it is very difficult to summarize when the two-beam pattern appears. It may depend on the landslide features (such as shape, initial depth and trajectory) and on the bathymetry. We will investigate this question in the future study.

The main challenge for the simulations is to reproduce the high waves at Balchik (point 5) and the low waves at Galata (point 7). We emphasize that in contrast to the simulations of Ranguelov *et al.* (2008) and Vilibic *et al.* (2010), landslides *L*_{1}–*L*_{6} generated no high waves along the southern coast of Bulgaria.

The synthetic marigrams for the seven coastal points of the simulation with landslide *L*_{1} are presented in Figure 5. It can be seen that the period of the oscillations was mainly 3–12 min, which is in satisfactory agreement with the observations. Note that at most points the first wave was positive. For comparison, we also present the marigrams for landslide *L*_{6}, where the first wave was negative at most points.

## Frequency dispersion effects

In numerical simulations of tsunamis, the significance of some physical effects (such as non-linearity and dispersion) are often evaluated. In particular, the non-linearity effects are known to be crucial in a coastal zone where wave heights become comparable with water depth. The significance of the frequency dispersion effects decreases with wave length and increases with water depth and propagation time (Glimsdal *et al.* 2013). Dispersion generally reduces the height of a leading wave and contributes to the appearance of a train of shorter waves which propagate with lower velocities. Landslide tsunamis are known (Lovholt *et al.* 2008; Tappin *et al.* 2008; Glimsdal *et al.* 2013) to be more dispersive than seismically generated tsunamis. For example, Gylfadottir *et al.* (2017) simulated the sub-aerial landslide tsunami in Askja Lake 2014 and concluded that the observed wave heights may be reproduced much better using the dispersive model.

The significance of the frequency dispersion effects in simulations of tsunami waves generated by an underwater landslide is still under debate, but obviously depends on the parameters of the problem such as shape of the landslide, its dynamics, water depth and propagation distance (time). Moreover, irregularities in the bottom relief can also influence the extent of dispersion.

In Figure 3, we also present results computed using a non-dispersive NLSW model. Note that computations with the numerical algorithm for the NLSW model take about 10–20% of CPU time needed for computations using the FNLD model. The computed wave heights occasionally differ significantly from the computations with the FNLD model. It can be concluded that dispersion is not the crucial factor in this particular problem. We believe that the main reason for such a weak display of dispersion is the consideration of the nearest coast only. In the previous study (Khakimzyanov *et al.* 2015), we noted that in simulations with large submarine landslides the dispersive effects became visible only after propagation of the tsunami through deep water. Again, the main differences in the computations between dispersive and non-dispersive models in Gylfadottir *et al.* (2017) were observed on the opposite coast with respect to the initial landslide position. Even in our computational domain, the synthetic marigrams recorded with gauges *M*_{1} and *M*_{2} for deep-water sites show significant differences in the FNLD and NLSW computations for the landslide *L*_{1} (Fig. 6). It can be seen that in the records of *M*_{1} the dispersion changed the form of the wave, but the maximum and minimum free-surface levels are similar in the computations. In contrast, the marigrams for *M*_{2} show significant differences both in the form of the waves and in the heights.

We should note here that the model landslides with *b _{x}* =

*b*= 2.5 km (i.e. two times ‘shorter’) generate more dispersive waves, but the thickness of the landslides becomes less realistic. Apparently, landslides with larger horizontal extensions and lower thickness will generate less dispersive waves.

_{y}## Conclusion

In this paper, we studied the ability of a submarine landslide to generate tsunami waves similar to those observed on the Bulgarian coast on 7 May 2007. In total, 211 initial positions of the landslide were considered in our simulations. We selected six model landslides which generated wave heights that fitted well with the field data. The remarkable fact is that all these landslides stopped at the same region. We recommend this region for underwater explorations.

Initial position of a landslide is crucial for its dynamics. In our computations, the absolute velocity maximum was usually 20–40 m s^{−1}, which led to a Froude number of 0.15–0.3. The landslide absolute acceleration was seldom greater than 0.25 m s^{−2}.

Considering the nearest coast with respect to the landslide position, it can be concluded that employment of dispersive models is unnecessary in such cases. However, the dispersion effects were visible even in the local computational domain. Note that these effects depend on many factors such as landslide form, its dynamics and bathymetry configuration; all have to be considered for every specific problem.

This work is an extension of the study of Ranguelov *et al.* (2008). It shows that a landslide could have generated the observed tsunami waves but, in contrast to Ranguelov *et al.* (2008), such a landslide could have started from the deep water. Moreover, it may have generated significant wave heights at the northern coast of Bulgaria and small wave heights at the southern coast. The peculiarity of the previous studies (Ranguelov *et al.* 2008; Vilibic *et al.* 2010) is the highest simulated waves at the southern coast where no reports of the tsunami were registered.

## Acknowledgements

We thank Richard N. Hiscott and an anonymous reviewer for their comments.

## Funding

The work was supported by the Russian Science Foundation project No. 14-17-00219 awarded to L. C. and the Council for Grants of the President of the Russian Federation for State Support Leading Scientific Schools (grant No. NSh-7214.2016.9) awarded to Y.I. Shokin.

- © 2018 The Author(s). Published by The Geological Society of London. All rights reserved