Mathematical and Information Technologies, MIT-2016 — Mathematical modeling Numerical Simulation of Surface Waves Arising from Underwater Landslide Movement Yury Zakharov and Anton Zimin Kemerovo State University, Krasnaya Street, 6, 650043 Kemerovo, Russia zaxarovyn@rambler.ru,sliii@mail.ru Abstract. The main aim of this paper is to construct a model of simul- taneous movement of the landslide, internal currents and surface waves that can come ashore. The idea of a multicomponent fluid movement is used in this paper. We consider soil, liquid and gas as components of non- homogeneous fluid. Movement of such fluid is described by the Navier- Stokes equations with variable density and viscosity and the convection- diffusion equations. Special ratios are used to calculate the density and the viscosity of the medium. The results of test calculations for two- dimensional problems of the wave generation are presented. Keywords: Navier-Stokes equations, surface wave propagation, land- slide movement, inhomogeneous fluid, multicomponent fluid 1 Introduction Tsunami waves generated by underwater and above-water landslides can be dan- gerous for buildings located on the shore and settled lands. Under natural condi- tions the underwater landslide is a movement of some soil mass along the slope of the bottom. Large volumes of moving mass generate surface waves that are close in their characteristics to the waves generated by tsunamigenic earthquake. The overview of historical landslides and tsunamis that they generated can be found in [1,2]. Construction of tsunami wave model generated by landslide movement can be divided into two tasks: construction of a model of wave propagation on a free surface and a model of landslide movement on the bottom of water basin. Free surface of fluid means the border between fluid and gas that is above it. Due to the fact that fluid density is several times greater than gas density, influence of gas on the movement of fluid is often neglected, and it is considered that fluid moves independently of gas movement or, in other words, “freely”. Models of wave fluid dynamics are the examples of such approach. They include shallow water theory equations, ideal fluid movement equations, etc [3]. In case of complex wave movements, big splashes and active interaction of two phases, multicomponent non-homogeneous medium is considered, where fluid and gas are separate components with their own values of density and viscosity [4]. 535 Mathematical and Information Technologies, MIT-2016 — Mathematical modeling According to the manner of discretization, mathematical models of multicom- ponent fluid movement can be divided into Lagrangian and Eulerian. Lagrangian methods are based on recording the equations of the medium movements in La- grangian coordinates connected with particles of the moving medium. A set of nodal points moving together with the medium can be used in order to get the discrete analogues of such equations. It can be grid points (grid Lagrangian methods [5,6]) or point particles that are not connected with each other by grid lines (meshless Lagrangian methods [7-9]). In this approach the position of the interphase boundary is tracked automatically. For discretization of medium movement equations in Euler approaches a fixed computational grid is usually used. At this, the interphase boundary moves on the grid, and special methods are applied in order to determine its position. These methods include MAC [10,11], VOF [12,13], Level Set [14,15]. In laboratory studies the underwater landslide can be imitated either by movement along the slope of a fully submerged solid body [16], or by some granulated soil slipping down the basin [17]. Several approaches are identified for computational modelling of landslide movement. It can be a model of movement of an absolutely solid body [18] or an ensemble of such bodies [19]; or a model of fluid flow that has different density, viscosity etc. [20]. Until recently, the majority of models applied for modelling of the tsunami of landslide type, relied on nonlinear shallow water wave theory [21,22]. In many cases the Boussinesq equations are applied [23]. And also the attempts are made to apply three-dimensional non-homogeneous models based on Navier-Stokes equations and the concentration transfer equations [24]. The aim of this paper is to apply the model of three-component viscous in- compressible fluid with variable viscosity and density and with the presence of mass diffusion between the components for the problems of occurrence of sur- face waves generated by landslide movement. Previously two-component model has been used in problems of substance diffusion in the branched channel [25], cohesive soil erosion [26] and surface wave propagation [27]. 2 Mathematical model The movement of the medium consisting of three incompressible interfusing fluids with constant density ρ1 , ρ2 , ρ3 and viscosities µ1 , µ2 , µ3 is considered. Let the dx particle of mixture be the solution x = x(t) of the Cauchy problem = V (x, t), dt x(0) = x0 , where V (x, t) = (v1 , v2 , v3 ) is a velocity vector of the mixture in point x = (x1 , x2 , x3 ) and t is a time point. Let C1 (x, t), C2 (x, t), C3 (x, t), µ and ρ be volume concentrations of the components, dynamic viscosity and mixture density correspondingly. Components concentrations are interconnected in the following way: C1 + C2 + C3 = 1, 0 ≤ Ci ≤ 1. (1) 536 Mathematical and Information Technologies, MIT-2016 — Mathematical modeling In order to find out the values of viscosity and density of the mixture we have the following dependence on the concentrations of the components:  µ1 µ2 µ3 µ = , C1 µ2 µ3 + C2 µ1 µ3 + C3 µ1 µ2 (2) ρ=C ρ +C ρ +C ρ . 1 1 2 2 3 3 Mass diffusion occurs between the particles of the mixture according to the law: ∂ρ qn = −D , (3) ∂n where D is the diffusion coefficient. We consider the mixture components to possess the incompressibility prop- erty and its interfusion to form an incompressible medium, which density de- pends only on the concentrations. Let ωt be a control moving volume of such medium. Value of ωt remains constant due to the incompressibility: ∫ dτ = const. ωt From known formula for the time differentiation of the integral taken over the moving volume [28] ∫ ∫ [ ] d df f dτ = + f divV dτ (4) dt dt ωt ωt we obtain the condition of incompressibility: divV = 0, (5) d where dt is total time derivative. The equations of mass balance for fluid volume ωt considering the presence of mass diffusion, take the following form: ∫ ∫ d ρ dτ = − qn dσ, (6) dt ωt ∂ωt where qn is defined in (3). From (6) taking into account (3) and (5) we get convection- diffusion equation for density: dρ = D∆ρ, (7) dt where ∆ is Laplace operator. (5) and (7) together provide the condition of mass balance conservation in the medium for three-component incompressible fluid. 537 Mathematical and Information Technologies, MIT-2016 — Mathematical modeling For our objectives we consider three-layered fluid, i.e. we suppose that the first and the third component do not directly interact in the solution: { C3 = 0, C1 ̸= 0, C1 = 0, C3 ̸= 0. Then the density equation (2) will be the following: { C1 ρ1 + C2 ρ2 , C3 = 0, ρ= C2 ρ2 + C3 ρ3 , C1 = 0. Or considering (1) { C1 ρ1 + (1 − C1 ) ρ2 , C3 = 0, ρ= (8) C3 ρ3 + (1 − C3 ) ρ2 , C1 = 0. The diffusion coefficient D can be also expressed as: { D = D1 , C3 = 0, (9) D = D3 , C1 = 0. Then (7) can be transformed into the following equations for the component concentrations:   dC  1 = D1 ∆C1 ,  C3 = 0,   dt dC3 (10)   = D3 ∆C3 , C1 = 0,   dt  C2 = 1 − C1 − C3 . From the integral momentum equation ∫ ∫ ∫ d ρV dx = Pn dσ + ρF dx (11) dt ωt ∂ωt ωt considering (4) and (5) we get: ( ) d ρV = divP + ρF , (12) dt where P is the stress tensor in the mixture, F = (f1 , f2 , f3 ) is the vector of mass forces. Then, considering variable viscosity, we obtain a system of equations for the motion of the mixture of three viscous incompressible interfusing fluids: 538 Mathematical and Information Technologies, MIT-2016 — Mathematical modeling  ( )   d ρV   = −∇p + divI + ρF ,   dt     divV = 0,      dC1    = D1 ∆C1 ,  dt dC3 (13)   = D3 ∆C3 ,   dt     C2 = 1 − C1 − C3 ,     µ1 µ2 µ3     µ= ,   C µ µ 1 2 3 + C 2 µ1 µ3 + C 3 µ1 µ2  ρ = ρ1 C1 + ρ2 C2 + ρ3 C3 . where p is pressure in the mixture, I = µ∆V + (∇µ · ∇) V + (∇µ · JV ) is viscous part of stress tensor, JV is Jacobian matrix, arranged as follows:  ∂v ∂v ∂v  1 1 1 ∂x ∂x ∂x  ∂v1 ∂v2 ∂v3  JV =  2   ∂x1 ∂x2 ∂x3  . 2 2 ∂v3 ∂v3 ∂v3 ∂x1 ∂x2 ∂x3 Thus, the given model consists of the convection-diffusion equations for con- centration of the components, correlations for the determination of density and viscosity, and also hydrodynamic Navier-Stokes equations for incompressible vis- cous fluid. We use a no-slip condition on the solid wall and boundary conditions of the second kind for concentration equations. Some initial distribution for concentra- tions is also given. 3 Solution scheme For discretization of the system (13) in the spacial variables is used a finite- difference method on a rectangular uniform grid with staggered arrangement of nodes [29]: pressure, velocity divergence and component concentration are determined in the centers of cells; velocity vector components are determined on the borders of cells. Application of the staggered grid allows to link speed and pressure values in the adjacent nodes and avoid the appearance of oscillations in solution, which arise when using central differences on a combined grid. Also, a staggered ar- rangement of the nodes automatically allows to satisfy the discrete representa- tion of the continuity equation. Time motion algorithm consists of the following stages: 1. By taking into account the known velocity vectors and concentration distri- bution (and thus density and viscosity values), a time step for the Navier- Stokes equations system is made. 539 Mathematical and Information Technologies, MIT-2016 — Mathematical modeling 2. Using the received values of velocity component, a time step for the “soil” convection-diffusion equation is made. 3. Using the received values of velocity component, a time step for the “air” convection-diffusion equation is made 4. Knowing the distribution for “air” and “soil” concentrations, a value of “wa- ter” concentration is calculated according to the formula (1). 5. Recalculation of density and viscosity values in the medium is carried out according to the formulae (2). After that follows a transition to the first step of the next time iteration. To solve the Navier-Stokes equations system there is used a splitting scheme on physical factors [30] with regard to variable density. It consists of three steps. At the first step a momentum transfer is carried out due to convection and dif- fusion; intermediate velocity field is calculated according to an implicit scheme: Ṽ − V n = − (V n · ∇) Ṽ + ∆t (14) 1( ) + −V n D∆ρ + µ∆Ṽ + (∇µ · ∇) Ṽ + (∇µ · JṼ ) + F , ρ In order to solve the system (14) a prediction-correction method is used [31]. The obtained system of algebraic equations is solved by sweep method. At this, despite the fact that the obtained intermediate velocity field Ṽ does not satisfy continuity equation, it has a physical meaning because it preserves vortex characteristics in internal points. At the second step, with regard to (5) and variable density, the pressure field is calculated according to the obtained intermediate velocity field Ṽ : ∑ ∂ ( 1 ∂pn+1 ) ∇Ṽ = . (15) i ∂xi ρ ∂xi ∆t Solution of the system of equations obtained as a result of discretization of the equation in order to find pressure in (15) is one of the most important and dominant aspects of computational procedure from the viewpoint of machine resources expenses. Operator of this system can have a complex structure, which complicates the task significantly. To solve this stage of computational process a gradient iterative method BiCGStab [32] is used. At the third step the transfer of momentum is carried out only due to pressure gradient: V n+1 − Ṽ 1 = − ∇pn+1 (16) ∆t ρ The equation (15) obtained by taking the divergence of both sides of equation (16) with regard to ∇V n+1 = 0. 540 Mathematical and Information Technologies, MIT-2016 — Mathematical modeling To solve the convection-diffusion equations (10) a prediction-correction scheme with approximation of convective components against the stream is used [31]. The obtained system of algebraic equations is solved by sweep method. The numerical scheme has first-order temporal and spatial approximations. 4 Results 4.1 Collapse of viscous soil There was considered a problem of collapse of viscous and stiff soil on the bot- tom of reservoir that generates waves on the surface of fluid. Here one of the components (more stiff and viscous) models the behavior of bottom soil, another one – liquid, and the third one – aerial environment. Fig. 1 shows the geometry of area and the scheme of component arrangement. Fig. 1. Geometry and initial distribution of components. At the initial time the half circle of the wet soil is located at the center of the area. As time passes, it caves under the influence of gravity F = (0, −9.8) sm2 and causes the movement of the entire medium. Fig. 2 shows the results of the calculation for various time points. Parts of the bottom soil spread in the opposite directions, then reflected from the side walls and connected again in the center of the area. The liquid phase surface followed the bottom soil movements. Fig. 2. Picture of medium motion for various time points t = 0.9 s, 1.9 s, 2.7 s, 5.4 s. 541 Mathematical and Information Technologies, MIT-2016 — Mathematical modeling The following hydrodynamic parameters were chosen: µ1 = 10 P a · s, µ2 = kg kg kg 10−3 P a · s, µ3 = 10−5 P a · s and ρ1 = 3000 m 3 , ρ2 = 1000 m3 , ρ3 = 1 m3 for the soil, liquid and gas phases. The following grid parameters and time step were used: hx = 5 · 10−2 m, hy = 5 · 10−2 m, τ = 10−2 s. All the area walls are solid except the upper one. The atmosphere pressure Patm = 101325 P a is indicated at the top. We consider the boundary between components to take place at C = 0.4. The calculation demonstrates possibility of the model to simulate direct in- teraction between the bottom soil and the surface waves without distinguishing characteristics at the phases boundaries. 4.2 Soil movement on the inclined bottom A problem of the soil slip movement on the inclined bottom that generates waves on the surface of fluid was considered. The scheme of area was taken from [33] (see Fig. 3). The result of calculation was compared with one of the numerical model presented in [33] and with laboratory experiment carried out in [34]. Fig. 3. Geometry and initial condition. At the initial time the wet soil is located on the left side of the inclined bottom. As time passes, it rolls down by gravity F = (0, −9.8) sm2 and causes the movement of the entire structure, simulating the soil slip movement. The kg hydrodynamic parameters were used the same as in [33]: µ1 = 10 m 3 , µ2 = −3 kg −5 kg kg kg kg 10 m3 , µ3 = 10 m3 and ρ1 = 1950 m3 , ρ2 = 1000 m3 , ρ3 = 1 m3 for the soil, liquid and gas phases. The following grid parameters and time step were used: hx = 5 · 10−2 m, hy = 5 · 10−2 m, τ = 10−2 s. All the area walls are solid except the upper one. The atmosphere pressure Patm = 101325 P a is indicated at the top. We consider the boundary between components to take place at C = 0.4. Fig. 4 shows the results of the calculation for various time points. Fig. 5 shows the comparison of water surface elevations with results obtained in [33] and [34]. Calculation demonstrates good agreement with MM3 [33]. This approach produces similar waveforms and slide deformation geometries. However, there 542 Mathematical and Information Technologies, MIT-2016 — Mathematical modeling Fig. 4. Picture of medium motion for various time points t = 0.0 s, 0.4 s, 0.8 s, 1.0 s. 0.10 0.05 Free surface height (m) 0.00 -0.05 -0.10 Model presented here MM3 (Smith et al., 2016) -0.15 Experiment (Assier- Rzadkiewicz et al., 1997) -0.20 -0.5 0 0.5 1 1.5 2 2.5 Distance (m) 0.10 0.05 Free surface height (m) 0.00 -0.05 -0.10 Model presented here MM3 (Smith et al., 2016) -0.15 Experiment (Assier- Rzadkiewicz et al., 1997) -0.20 -0.5 0 0.5 1 1.5 2 2.5 Distance (m) Fig. 5. A comparison of water surface elevations for various time points t = 0.4 s, 0.8 s. Given model (solid red), MM3 (dashed blue, [33]) and experiment (dotted green, [34]). 543 Mathematical and Information Technologies, MIT-2016 — Mathematical modeling is a discrepancy in the fields of high gradients, what can be explained by the form of equation (2) for µ. Differences do not exceed 12% on the whole surface for t = 0.4 s; 20% for t = 0.8 s. The surface shape in numerical simulation is qualitatively similar to that observed in the laboratory experiment. Also it has a slightly better agreement with experiment surface then MM3. 5 Conclusion Presented model of three-component viscous incompressible fluid was used for modeling simultaneous movement of the landslide, internal currents and surface waves. The advantage of this approach is that it allows one to simulate the complex phenomenon of the interaction of waves and bottom soil using a uniform numerical algorithm for a number of different problems without distinguishing characteristics at the phases boundaries. Test calculations for two-dimensional problems of the wave emergence and propagation on the free surface were carried out. The results obtained show the possibility of modeling such a phenomenon. Agreement with other model and experiment was demonstrated. Acknowledgments. The work was carried out with support of state task of Ministry of Science and Education, project number 1.630.2014/K. References 1. Langford, P.S.: Modeling of tsunami generated by submarine landslides. PhD Thesis, University of Canterbury, New Zealand (2007) 2. Papadopoulos, G. A., Kortekaas, S.: Characteristics of Landslide Generated Tsunamis From Observational Data. Submarine Mass Movements and Their Conse- quences Advances in Natural and Technological Hazards Research, 367–374 (2003) 3. Gusev, O. I., Shokina, N. Y., Kutergin, V. A., Khakimzyanov, G. S.: Modelirovanie poverkhnostnykh voln generiruemykh podvodnym opolznem v vodokhranilishche [Numerical modelling of surface waves generated by underwater landslide in a reser- voir]. Vychislitelnyye tekhnologii Vol. 18, No. 5, 74–90 (2013) 4. Khrabry, A. I.: Chislennoe modelirovanie nestatsionarnykh turbulentnykh techeniy zhidkosti so svobodnoy poverkhnost’yu [Numerical simulation of unsteady turbulent flow of liquid with a free surface]. PhD Thesis, Sankt-Peterburg, 154 (2014) 5. Butler, T.D.: LINC method extensions. Proceedings of the Second International Conference on Numerical Methods in Fluid Dynamics Lecture Notes in Physics Vol. 8, 435–440 (1971) 6. Hirt, C.W.: An arbitrary lagrangian-eulerian computing method for all speeds. Jour- nal of Computational Physics Vol. 14, 227–253 (1974) 7. Gingold, R.A., Monaghan, J.J.: Smoothed Particle Hydrodynamics: Theory and Application to Non-Spherical Stars. Monthly Notices of the Royal Astronomical Society Vol. 181, 375–389 (1977) 8. Lucy, L.B.: A Numerical Approach to the Testing of Fusion Process. The Astro- nomical Journal Vol. 82, No. 12, 1013–1024 (1977) 544 Mathematical and Information Technologies, MIT-2016 — Mathematical modeling 9. Afanas’ev, K.E., Makarchuk, R.S.: Calculation of hydrodynamic loads at solid boundaries of the computation domain by the ISPH method in problems with free boundaries. Russian Journal of Numerical Analysis and Mathematical Modelling Vol. 26, No. 5, 447–464 (2011) 10. Harlow, F.H., Welch, J.E.: Numerical Calculation of Time-dependent Viscous In- compressible Flow of Fluid with Free Surface. Phys. Fluids (American Institute of Physics) Vol. 8, No. 12, 2182–2189 (1965) 11. McKee, S., Tome, M.F., Ferreira, V.G., Cuminato, J.A., Castelo, A., Sousa, F.S., Mangiavacchi, N.: The MAC method. Computers Fluids Vol. 37, 907–930 (2008) 12. Hirt, C.W.: Volume of fluid (VOF) method for the dynamics of free boundaries. Journal of Computational Physics Vol. 39, 201–226 (1981) 13. Yakovenko, S.N., Chan, K.S.: Approksimatsiya potoka ob”emnoi fraktsii v techenii dvukh zhidkostei [Approximation of volume fraction stream in flow of two liquids]. Teplofizika i aeromekhanika Vol. 15, No. 2, 181–199 (2008) 14. Tonkov, L.E.: Chislennoe modelirovanie dinamiki kapli vyazkoi zhidkosti metodom funktsii urovnya [Numerical simulation of the dynamics of viscous liqued drop by level set method]. Vestnik Udmurtskogo universiteta No. 3, 134–140 (2010) 15. Nikitin, K.: Realistichnoe modelirovanie svobodnoi vodnoi poverkhnosti na adap- tivnykh setkakh tipa vos’merichnoe derevo [Realistic simulation of the water surface on adaptive grids of octal tree type]. Nauchno-tekhnicheskii vestnik SPbGU ITMO Vol. 70, No. 6, 60–64 (2010) 16. Enet, F., Grilli, S. T.: Experimental study of tsunami generation by threedimen- sional rigid underwater landslides. Waterway, Port, Coastal and Ocean Engineering Vol. 133, 442–454 (2007) 17. Mohammed, F., Frits, H. M.: Experiments on tsunamis generated by 3D Granu- lar Landslides. Submarine Mass Movements and Their Consequences, Advances in Natural and Technological Hazards Research Vol. 28, 705–718 (2010) 18. Watts, P., Imamura, F., Grilli, S. T.: Comparing model simulations of three bench- mark tsunami generation cases. Science of Tsunami Hazards Vol. 18, No. 2, 107–123 (2000) 19. Tinti, S., Bortolucci, E., Vannini, C. A.: Block-based theoretical model suited to gravitational sliding. Natural Hazards Vol. 16, 1–28 (1997) 20. Smith, R. C., Hilla, J., Collinsa, G. S., Piggotta, M. D., Kramera, S. C., Parkinsona, S. D., Wilsond, C.: Comparing approaches for numerical modelling of tsunami gen- eration by deformable submarine slides. Ocean Modelling Vol. 100, 125–140 (2016) 21. Beyzel’, S. A., Khakimzyanov, G. S., Chubarov, L. B.: Modelirovanie poverkhnos- tnykh voln, porozhdaemykh podvodnym opolznem, dvizhushchimsya po pros- transtvenno neodnorodnomu sklonu [Simulation of surface waves generated by an underwater landslide moving along a spatially irregular slope]. Vychislitelnyye tekhnologii Vol. 15, No. 3, 39–51 (2010) 22. Petrukhin, N. S., Pelinovsky, E. N.: Riemann waves in the dynamics of the land- slides above plane slope. Sovremennye problemy nauki i obrazovaniya, No. 6 (2011) 23. Watts, P., Grilli, S. T., Kirby, J. T., Fryer, G. J., Tappin, D. R.: Landslide tsunami case studies using a Boussinesq model and a fully nonlinear tsunami generation model. Natural Hazards and Earth System Sciences Vol. 3, 391–402 (2003) 24. Kozelkov, A. S.: Modelirovanie voln tsunami kosmogennogo i opolznevogo proiskhozhdeniya na osnove uravneniy Nav’ye-Stoksa [Modeling tsunami waves of landslide and cosmogenic origin based on the Navier-Stokes equations]. PhD Thesis, Nizhniy Novgorod, 401 (2016) 545 Mathematical and Information Technologies, MIT-2016 — Mathematical modeling 25. Gummel, E.E., Milosevic, H., Ragulin, V.V., Zakharov, Yu.N., Zimin, A.I.: Motion of viscous inhomogeneous incompressible fluid of variable viscosity. Zbornik radova konferencije MIT 2013, 267–274 (2014) 26. Zakharov, Yu., Zimin, A., Nudner, I., Ragulin, V.: Two-Component Incompress- ible Fluid Model for Simulating the Cohesive Soil Erosion. Applied Mechanics and Materials Vols. 725–726, 361–368 (2015) 27. Zakharov, Y., Zimin, A., Ragulin, V.: Two-Component Incompressible Fluid Model for Simulating Surface Wave Propagation. Mathematical Modeling of Technological Processes Vol. 549, 201–210 (2015) 28. Sedov, L. I.: Mekhanika sploshnoj sredy [Continuum mechanics]. Nauka, Vol. 1 (1970) 29. Patankar, S.: Numerical heat transfer and fluid flow. Hemisphere Publishing Cor- poration (1980) 30. Belotserkovskiy, O.M.: Chislennoye modelirovaniye v mekhanike sploshnykh sred [Numerical modeling in continuum mechanics]. Fizmatlit, Moscow (1994) 31. Yanenko, N.N.: Metod drobnykh shagov resheniya mnogomernykh zadach matem- aticheskoy fiziki [Method of fractional steps for solving multidimensional problems of mathematical physics]. Nauka, Novosibirsk (1967) 32. van der Vorst, H. A.: Bi-CGStab: a fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems. SIAM Journal on Scientific and Statistical Computing Vol. 13, 631–644 (1992) 33. Smith, R., Hill, J., Collins, G., Piggott, M., Kramer, S., Parkinson, S., Wilson, C.: Comparing approaches for numerical modelling of tsunami generation by deformable submarine slides. Ocean Modelling Vol. 100, 125–140 (2016) 34. Assier-Rzadkiewicz, S., Mariotti, C., Heinrich, P.: Numerical simulation of subma- rine landslides and their hydraulic effects. J. Waterw. Port Coast. Ocean Eng. Vol. 123, 149–157 (1997) 546