A model of mudflow propagation downstream from the Grohovo landslide near the city of Rijeka ( Croatia )

Mudflows regularly generate significant human and property losses. Analyzing mudflows is important to assess the risks and to delimit vulnerable areas where mitigation measures are required. The smoothed-particle hydrodynamics (SPH) model adopted here considers, in two phases, a granular skeleton with voids filled with either water or mud. The SPH depth-integrated numerical model (Pastor et al., 2009a) used for the present simulations is a 2-D model capable of predicting the runout distance, flow velocity, deposition pattern and the final volume of mudflows. It is based on mathematical and rheological models. In this study, the main characteristics of mudflow processes that have emerged in the past (1908) in the area downstream of the Grohovo landslide are examined, and the more relevant parameters and attributes describing the mudflow are presented. Principal equations that form the basis of the SPH depth-integrated model are reviewed and applied to analyze the Grohovo landslide and the propagation of the mudflow wave downstream of the landslide. Based on the SPH method, the runout distance, quantities of the deposited materials and the velocity of mudflow progression which occurred in the past at the observed area are analyzed and qualitatively compared to the recorded consequences of the actual event. Within the SPH simulation, the Newtonian rheological model in the turbulent flow regime and the Bingham rheological model were adopted and a comparison was made of the application of the Egashira and Hungr erosion law.


Introduction
In this study, a portion of the Rječina River near the city of Rijeka (Croatia), which was affected by a 1908 mudflow event, was used to investigate and determine the possible flow phenomena of unbounded fine-grained material.What is known from the records is that the mudflow event was initiated by the Grohovo landslide, near the Grohovo village in which several families lived.The mudflow event had a great significance for the Rječina River catchment area where in the early 20th century several washing mills and workshops were built.The mudflow event was caused by heavy rainfall over a short period of time (estimated value was around 220 mm in 7 h), but it was also affected by an earlier rockmass instability near the Grohovo village.According to the historical records (Croatian State Archive in Rijeka, JU 49 -Box 13, JU 51 -Box 45), the mudflow event lost its momentum in the middle of the canyon part of the Rječina River, between the Pašac bridge and the Žakalj village.According to the present terrain configuration of the Rječina watercourse (which has not significantly changed in the meantime), the runout distance of mudflow propagation was estimated to have been between 2300 and 2500 m.
Numerous historical records, images and maps describing the history of landslides in the area surrounding the village of Grohovo in the Rječina River valley were found in the Croatian State Archive in Rijeka (Hungarian Royal Cultural-Engineering Office of 1st District, 1998;Benac et al., 2002Benac et al., , 2005Benac et al., , 2009;;Oštrić et al., 2011;Arbanas et al., 2010;Vivoda et al., 2012;Žic et al., 2014).Sliding was first recorded in 1758 after the appearance of a large number of slips and landslides caused by an earthquake in 1750 with its epicenter in Rijeka.Significant sliding caused by rainfall and flood-Published by Copernicus Publications on behalf of the European Geosciences Union.ing was recorded on both banks of the Rječina River near the village of Grohovo at the end of the 19th century (Žic et al., 2014).A large slide occurred in 1870 on the SW of the hillslope and was again reactivated in 1885 (Fig. 1).On that occasion, a large portion of Grohovo was buried by a rock avalanche.A huge landslide was triggered on the SE slope of the Rječina River in 1853 at the location of the current active landslides (Benac et al., 2002(Benac et al., , 2005)).The channel of the Rječina River shifted approx.50 m to the south.Numerous landslides occurred during the first half of the 20th century without significant consequences.New landslides occurred during the construction of the Valići dam in 1960, when a landslide occurred on the NE slope immediately adjacent to the dam.In the northeastern valley, the largest active landslide along the Croatian Adriatic Sea region was reactivated in December 1996 within the landslide body from 1893.
Comprehensive rehabilitation of that landslide was never implemented, but further extent of the sliding body was significantly reduced.
According to the classification of mass movement types as proposed by Varnes in 1978, which was later modified by Cruden and Varnes (1996) and refined by Hutchinson (1988) and Hungr et al. (2001), flow is one of the basic features of landslide and can be divided into rock flows and soil flows.Soil flows can be classified as debris flows, debris avalanches, earth flows or mudflows.According to a further, more detailed classification of landslide types given by Varnes (1978) and Hungr et al. (2014), flow can be divided into rock flow (rock creep), debris flow (talus flow, debris flow, debris avalanche, solifluction and dry sand flow) and earth flow (dry sand flow, wet sand flow, quick clay flow, earth flow, rapid earth flow and loess flow).
Mudflow is defined as the propagation of fine-grained (silty) material whose composition (silt and/or clay) has greater plasticity and whose liquid index during movement is greater than 0.5 (Hutchinson, 1971;Laigle and Coussot, 1997;Komatina and Ðor dević, 2014).Mudflow represents a very rapid to extremely rapid flow of saturated, plastic, finegrained material in the channel, including significant water content in proportion to the source material (index of plasticity I P > 5 %) (Hungr et al., 2001(Hungr et al., , 2014;;Iverson, 1997).The velocity of mass movement can range from 0.5 to 15 ms −1 , but this limit may be exceeded in some extreme events, with flow reaching a maximum velocity of 25-30 ms −1 .The degree of fluidity was determined by the observed movement velocity or by the distribution and morphology of the sediments formed.Mudflows belong to a gradation series of processes involving water, clay and rock debris (rock fragments) in various proportions.The water content in mudflows can reach 60 %.The degree of water binding, determined by the clay content (particles the size of clay) and the mineralogy of the solid particles (mineral composition of the particles), has a critical effect on the viscosity of the matrix (mixture) and on the flow velocity and morphology (Hungr et al., 2014).
One of the most significant geomorphological features of mudflow is the total travel distance, which is defined as the length of travel path over which the flow of unbound grained materials are in interaction with water (Varnes, 1978;Cruden and Varnes, 1996;Fannin and Wise, 2001).When describing the mudflow, two categories of parameters should be considered: terrain properties and flow properties.Terrain properties are characterized by the ground surface slope and the erodibility of the channel bottom.Flow properties include the sediment concentration, density of particles, amount of water, flow velocity, and parameters that describe the stress and the initial and final (deposited) volume of the mudflow materials (Laigle et al., 2007;Blanc, 2008).In general, the output parameters of the mudflow numerical simulation are flow velocity, flow depth, total deposited volume and runout distance of the muddy deposited material.
The main threat in the Rječina River valley is that landslides could cause possible rearrangement of riverbeds and the creation of a natural lake.Due to large amounts of rainfall, such a lake would fill rapidly, and the accumulated water would then overtop the dam built by the sliding mass.After the collapse of the dam due to overflow, the flood wave would then pass through a narrow canyon of the Rječina River (near the village of Pašac) in its lower section, which could potentially cause the loss of human life and serious damage to buildings in the central part of Rijeka (Oštrić et al., 2011;Žic et al., 2013a).Additional danger lies in the possible occurrence of landslides on the slopes above the Valići accumulation (useful capacity 0.47 million m 3 , located approx.300 m upstream of the Grohovo landslide), if this flysch mass were to slide into the accumulation with significant consequences.Heavy precipitation (> 100 mm) or earthquake events, separately or in combination, might become effective triggers of mudflows.

Geomorphological, geological and hydrological properties of the study area
The dominant tectonic structure in the study area in the Rječina River valley is a portion of a major geomorphological unit that strikes in the direction Rječina River valley-Sušačka Draga valley-Bakar Bay-Vinodol valley (Blašković, 1999;Benac et al., 2002Benac et al., , 2005Benac et al., , 2011)).The Rječina River extends through three distinctive geomorphological units.The first geomorphological unit extends from the karstic spring of the Rječina River in the foothills of the Gorski Kotar mountains to the village of Lukeži; the second from Lukeži to the entrance of a portion of the Rječina River canyon; and the third from that canyon to the alluvial plain at the mouth of the Rječina River in the center of Rijeka.The upstream and central sections of the Rječina River valley are relatively narrow and formed in Paleogene flysch.This portion of the valley also consists of Upper Cretaceous and Paleogene limestone.The downstream section of the watercourse flows through a deep canyon cut into Cretaceous and Paleogene carbonate rocks (Benac et al., 2005(Benac et al., , 2011)).The central section of the watercourse, between the Valići Dam and the Pašac Bridge, is 1.8 km long and 0.8 to 1.1 km wide, as shown in Fig. 2.
The origin of a landslide is preconditioned by the geological structure and morphogenesis of the Rječina River valley.The Rječina River valley is geomorphologically younger than other nearby valleys formed in flysch.Due to its geological and morphological conditions, both slopes in the Rječina River valley between the villages of Drastin and Pašac are on the boundary of a stable equilibrium state.
The flysch bedrock is characterized by its heterogeneity, with frequent vertical and lateral alternations of different lithological sequences.Microscopic petrological analysis of the bedrock showed the presence of silty marl, laminated silt to silty shale and fine-grained sandstone.From the orientation of the sandstone layers, the flysch appears to strike towards the northwest, i.e., downslope.Analysis of the soil indicates that silt is the dominant size fraction, although the clay fraction is also significant, varying between 17 and 38 % (Fig. 3).
To obtain mineralogical, physical and mechanical properties of the soil and rock materials from the Grohovo landslide body, 22 representative samples were selected from the flysch deposit, 18 of which were taken from the drilling cores (1999), while the remaining four were taken from the ground surface in 2006 (Benac et al., 2014).Further analysis of fine-grained fractions (up to 1 mm) were conducted for the mineralogical analysis.The standard geotechnical laboratory tests were conducted on the 13 samples of borehole and four on the surface samples.Grain-size analysis was performed according to the methods of screening and hydrometric following for ASTM standard (IGH, 2000).
Sedimentological analysis of the grain size (Fig. 3a) and geotechnical analysis (Fig. 3b) indicates that in all samples the silt and clay are dominant.It can therefore be concluded that the investigated area is characterized by clayey silt or muddy clay.Figure 3b shows that the average particle size (D 50 ) ranges from 0.004 to 0.042 mm and in the analysis of sediment grain size from 0.0028 to 0.056 mm.Index of plasticity of the tested soil was in the range of I P = 14-22 %, from which it can be concluded that the material has a low to medium plasticity.The liquid limit was in the range of W L = 32-43 %.Quantitative mineralogical analysis of the material composition of the samples has revealed the presence of the following clay minerals: kaolinite, illite, chlorite, mixed-layer clay minerals and -in some samples -vermiculite and smectite (IGH, 2000) (Fig. 4).
Quartzite, calcite and phyllosilicates constitute 86-96 % of the mineral composition across various samples.Laboratory test results using a direct shear test on eight samples have shown measured peak values of the friction angle in the range of 23.7 • < ϕ < 26.1 • and the cohesion within the range of 1 < c < 9.5 kPa (Benac et al., 2014).Based on the laboratory tests results it can be argued with high probability that silty-clay materials are prevalent within the lower part of the colluvial material from the landslide body of landslide.
The section from the spring of the Rječina River to the Grohovo landslide has a meandering shape, low longitudinal slope (approx.5-7 %) and a U-shaped cross section.From the Grohovo landslide to the mouth of the canyon, the Rječina riverbed has a V-shaped cross section and a steep slope (approx.20-30 %).The deposits are large and dragged, and folds are common (the Žakalj folds cause a waterfall).From the Rječina River canyon to the mouth into the sea, the slope of the riverbed decreases approx.4-6 %, and the riverbed was carved into carbonate rock mass.The flow from the total catchment area of Rječina River runoff into the river, which corresponds to the hydrometric profile at Grohovo (194.3 m a.s.l.), includes more than 75 % of the average rainfall for the catchment area of 2250 mm (Ri danović, 1975).
The basin of the Rječina River extends NW-SE.The altitudes in the basin are in the range 0-649 m a.s.l.(above sea level), and the slope generally varies in the range from 0 to 30 • .The Rječina is a typical karstic river originating from a strong karstic spring located at the foot of the Gorski Kotar mountains (325 m a.s.l.).The watercourse is 18.63 km long and has a direct (orographic) catchment area of approx.76 km 2 , but the catchment area of all sources that feed the Rječina and its tributaries is much larger, approx.400 km 2 .The annual average flow of the Rječina spring is 7.76 m 3 s −1 , with maximal flow rates ranging from 0 to over 100 m 3 s −1 (Karleuša et al., 2003).The Rječina River has a few tributaries (Sušica, Mudna Dol, Lužac, Zala, Zahumčica, Golubinka, Ričinica, Borovščica and Duboki jarak) with Sušica the most important tributary (Fig. 2).After the catastrophic flood in 1898, extensive channel regulation was performed in the upper central section of the Rječina watercourse.The majority of the regulation work was completed to reduce flood effects and consisted of transversal structures to prevent deepening of the channel and the formation of landslides (Žic et al., 2014).
Significant, very intensive, short-term rainfall events greatly influence both the surface and groundwater discharge (Fig. 5).The entire area is occasionally subject to very intense rainstorms, which can cause serious damage through flash floods and mass movements.
The natural groundwater flow velocity ranges from 0.2 to 4 cm s −1 , and the hydraulic gradient varies from 0.03 to 0.06 (Biondić, 2000).One indicator of the complexity is the discrepancy between the amount of rainfall and the river network density, which amounts to 0.2 km km −2 in this drainage area (Knežević, 1999).Runoff on the slopes is mostly present in the flysch area in the middle of the basin.The springs at the foot of the landslide remain active even in dry periods.Their capacity is estimated at 2 L s −1 in the dry period and more than 20 L s −1 in the rainy period.A spring with  Depth (m) 5.0 1.0 2.5 4.5 3.0 3.4 7.0 10.3 14.0 1.0 4.5 8.5 0.0  (Benac et al., 2014, modified by Elvis Žic). a capacity of 30 L s −1 was also observed at the foot of the coarse-grained slope deposits after periods of intense precipitation.The groundwater level changed by less than 67 cm in two boreholes (G-5 and G-7) located in the upper part of the landslide but varied by several meters in the boreholes (G-1 and G-3) (Žic et al., 2013b) in the lower part of the landslide, as observed in Figs. 2 and 6.Measurements of groundwater levels have been realized by Mini Diver instruments, used to measure groundwater levels and temperature, connected to a wired ribbon down to the bottom of galvanized steel piezometers of circular shape with a diameter of 10 cm.Installation depth of piezometers G5 and G7 was 8-12 m, that of piezometer G3 was about 9 m and that of piezometer G1 was 6 m (viewed from the ground level at the site of embedded piezometers).
In torrential watercourses such as Rječina, floods are not unusual.Large variations in the discharge, short flood-wave propagation time, high sediment transport and the narrow corridor available for the evacuation of flood waves require a specific approach to flood control problems.One such method is numerical modeling of flood wave propagation, which enables water management professionals to examine various possible flood scenarios and, by varying different parameters directly affecting the occurrence of floods, to select the optimum solution for the protection of the city of Rijeka.

Simulation framework -smoothed-particle hydrodynamics (SPH) method
In recent decades, modeling of the propagation stage has been largely performed within the framework of continuum mechanics, and a number of new and sophisticated computational models have been developed.Most of the available approaches treat the heterogeneous and multiphase moving mass as a single-phase continuum.Mesh-free methods provide accurate and stable numerical solutions for integral equations or partial differential equations (PDEs) with a variety of possible boundary conditions and a set of arbitrarily distributed nodes (or particles) without using a mesh to provide the connectivity of these nodes or particles (Monaghan and Latanzio, 1985;Monaghan, 1992Monaghan, , 1994;;Monaghan and Kocharyan, 1995;Libersky and Petschek, 1990;Libersky et al., 1993;Liu and Liu, 2003;Liu, 2009).Smoothed-particle hydrodynamics is one of the mesh-free particle methods that was originally proposed for modeling astrophysical phenomena and was later widely extended for applications to problems of continuum solid and fluid mechanics (Lucy, 1977;Gingold and Monaghan, 1977).In the SPH Lagrangian method, the state of a system is represented by a set of particles that possess individual material properties and move according to the governing conservation equations (Liu and Liu, 2003).
The SPH 2-D depth-integrated numerical model is adopted here (code by M. Pastor, 2007 version) (Pastor, 2007).The model is capable of predicting the runout distance of mudflow, flow velocity, composition of the deposition and final volume of mudflow (Pastor et al., 2009a, b;SafeLand project, 2012;Cascini et al., 2014).The basis of the mathematical model is linking the depth-integrated model of the connection between the flow velocity and the pressure using Biot-Zienkiewicz equations.The rheological modeling corresponds to the constitutive equations.
The formulation of SPH is often divided into two key steps.The first step is the integral representation or the socalled kernel approximation of the field functions.The second step is the particle approximation.In the first step, the integration of the multiplication of an arbitrary function and a smoothing kernel function gives the kernel approximation in the form of the integral representation of the function (Gingold and Monaghan, 1982;Oñate and Idelsohn, 1998;Liu, 2009).The integral representation of the function is then approximated by summing up the values of the nearest neighbor particles, which yields the particle approximation of the function at a discrete point or a particle (Vignjević, 2002;Liu and Liu, 2003;Li and Liu, 2004;Hitoshi, 2006).

Mathematical model
This section is largely based on the work of Pastor (2007) and Pastor et al. (2009a) and is included here for complete-ness.Soils are geomaterials with pores that can be filled with water, air and other liquids.They are, therefore, multiphase materials with a mechanical behavior that is regulated by all phases.When the soil is considered a mixture, the continuity equation, momentum balance equations and the constitutive equations can be formulated for each phase.Darcy's relative velocity (ω α ), which represents the velocity of the liquid phase with respect to the velocity of the solid phase, connects the liquid phase velocity (v α ) with the solid phase velocity (v s ).The total Cauchy stress, σ , within the mixture can be separated into solid phase stress, σ (s) , and pore liquid phase stress, σ (w) .Simultaneously, the pore air phase stress, σ (a) , is usually separated in continuum mechanics into hydrostatic and deviatoric components.Generally, all three phases (solid, liquid and air) are present in the soil mixture; hence the total Cauchy stress can be represented by three partial stresses: where the incides s, w, a refer to the partial stresses in the solid, water and air phases; p is the average pressure; σ is the effective stress; n represents porosity; s α = dev(σ α ) and S α stand for the deviatoric stress component and the degree of saturation, respectively, for the liquid and air phase (labeled "nphases"); and I represents the identity tensor of the second order.For more details, readers are also referred to Blanc (2008).The general model consists of the following equations: 1. the mass-balance equations for the solid and liquid phases: (3) 2. the momentum-balance equations for the solid and liquid phases: where b is external force and k α is the permeability (leakage) of phase α; 3. the kinetic equations that connect the velocity to the strain rate tensor: where D is the rate of deformation tensor.Assuming that the relative velocities between the fluid phase and its acceleration are small, model v-p w can be formulated as a function of the velocity of a solid skeleton and the relative velocity of the fluid within the skeleton (Blanc, 2008;Pastor 2009a, b;Blanc et al., 2011).
Rapid flow includes two physical phenomena: the consolidation and dissipation of the pore pressure and the propagation.Axes x 1 and x 2 are on a slope near the plane, or horizontal axes, whereas axis x 3 is normal (perpendicular) to the plane (Fig. 7).
Following Pastor et al. (2009a), it is assumed that the velocity can be separated as v = v0 + v1 and the pore pressure is decomposed as pw = pw o + pw 1 .In this way, v 1 can be identified as the velocity corresponding to the 1-D consolidation, and v 0 is the velocity of propagation (Blanc, 2008;Haddad et al., 2010).The propagation-consolidation model consists of a set of partial differential equations.Equations are integrated along the normal direction of the surface using the Leibnitz and Reynolds theorem (Pastor et al., 2009a).
The erosion is considered by introducing the erosion rate, e r = − ∂z ∂t , which yields ∂ ∂t (h + z) = ∂h ∂t − e r and must be integrated into the mass-balance equation.Therefore, the depthintegrated mass-balance equation is The linear balance momentum equation is integrated over the depth and yields assuming that the stress on the surface equals 0, and the stress at the bottom of the channel is τ B = −ρ g h grad Z − τ b .The model considers the existence of saturated layers of the height, h s , at the bottom of the flow (Hungr, 1995).Therefore, a reduction of the pore pressure is caused by the vertical consolidation of this layer.Finally, the depth-integrated consolidation equation has the following form: where c v = 0.000006 m 2 s −1 is accepted as the coefficient of consolidation (Sridharan and Rao, 1976;Olson, 1986;Robinson and Allam, 1998) 1-D consolidation equation.The resulting mass-balance, momentum-balance and pore pressure dissipation equations are ordinary differential equations, which can be integrated in time using a scheme such as leap-frog or Runge-Kutta (2nd or 4th order).The results depend on the rheological model chosen, from which it is possible to obtain the basal friction and the depth-integrated stress tensor.Further details may be found in Pastor et al. (2009a) and Blanc et al. (2011).

Rheological models
For a full simulation framework, mathematical models need to be completed by defining constitutive or rheological models.The best-known model is the Bingham viscoplastic model (Bingham and Green, 1919;Cantelli, 2009;SafeLand project, 2012;Calvo et al., 2014), which is used for mudflow modeling.In the case of Bingham fluids, the shear stress on the bottom as a function of the averaged velocity cannot be directly obtained.The expression relating the averaged velocity to the basal friction for the infinite mudflow problem is given as where µ is the coefficient of dynamic viscosity, τ Y is the yield stress and τ B is the shear stress on the bottom (Pastor et al., 2009a(Pastor et al., , 2015;;Blanc et al., 2011;SafeLand project, 2012;Calvo et al., 2014).Most depth-integration models use simple rheological laws because of the difficulty of their implementation.The friction model is one such simple model.It follows from the model by Cheng and Ling by neglecting cohesion and viscous terms (Cheng and Ling, 1996), with the vertical distribution of the shear stress τ (z) = ρ g(h − z) sin θ and the Mohr-Coulomb strain s(z) = ρ d g(h − z) cos θ tan ϕ.The symbol h is the depth of flow, z is the elevation, θ represents the slope angle and ρ d is the submerged particle density, equal to ρ s -ρ w .With respect to the base friction, the pore pressure is included as Based on the latter equation, the pore pressure can be concluded to have an effect similar to the reduction of the friction angle.

Erosion
Consideration of erosion activity requires a rheological or constitutive behavior of the interface and it depends on the variables such as the flow structure, density, particle size and on how close the effective stresses at the surface of the terrain are to failure (Iverson, 1997;SafeLand project, 2012;Cuomo et al., 2014).In this study, the erosion laws of Hungr and Egashira were adopted.
The Hungr law employs the erosion rate, which increases in proportion to the depth of flow, resulting in proportional distribution of the depth of the input material and the exponential growth of the mudflow with displacement.Changes in the stress conditions lead to a collapse of the bottom of the flow route and an engagement of material proportional to the change in the total normal stresses on the channel bottom (Hungr, 1990(Hungr, , 1995;;Hungr and Evans, 1997).The empirical law was based on the erosion rate of displacement E s , the so-called "growth rate" (Blanc, 2008).This parameter represents the normal depth of the eroded bottom per unit of flow and displacement.The Hungr law consists of the relationship between the erosion rate, e r , and the rate of growth, E s (Blanc, 2008;SafeLand project, 2012).
The Egashira erosion law ( 2001) is based on the tests of the inlet channel as well as on the numerical and dimensional analysis.Egashira assumed that the slope of the channel bottom is always aligned when mudflow is traveling through the erodible bottom of the channel (Egashira et al., 2001).The Egashira erosion law appears in the form e r = c * v tan(θ − θ e ), where c * is the concentration of sediment volumes of the sediment bottom (of the stationary layer), θ is the slope of the channel bottom and θ e represents the balance slope of the bottom.

Application of the SPH method on the Grohovo landslide
In this study, a Newtonian fluid model for turbulent regimes (Pastor, 2007;SafeLand project, 2012) and the Real Bingham fluid model (Pastor et al., 2004(Pastor et al., , 2007;;Calvo et al., 2014) are used to simulate the mudflow propagation.The choice between the Hungr and Egashira erosion laws for modeling erosion processes within the SPH method was considered.Because of the causes of the instability of the slopes in the Rječina River basin, the topography is provided by a digital elevation model (DEM) that was created on the Geographic Information System platform (ArcGIS 10.1 version) with the equidistant mesh grids of 2, 5 and 10 m (Fig. 8).The digital elevation model of the terrain is used to create the simulation of unbound fine-grained material propagation using the SPH algorithm (Pastor Code -version from 2007).
As stated earlier, the objective of the simulation is to gain a clearer picture of the mudflow which occurred in 1908 in the area downstream of the Grohovo landslide and its propagation to the urban part of Rijeka.With the help of the obtained visualization of the simulation, the volume of deposited finegrained material of the mudflow, the wave velocity of propagation, the depths of the deposited materials and the scope of mudflow in the analyzed area were quantified.The present analysis allows the quantification of the individual input parameters that initiate the formation of a mudflow.Parameters defined by well-established relationships should enable a correlation between the geomorphological and hydrogeological conditions and the identification of the specific field conditions with soil characteristic parameters that may lead to the formation of a mudflow.Defining the critical geomorphological and hydrogeological parameters of the soil that encourage the emergence of a mudflow on a flysch area will allow the assessment of hazards and mitigation measures.
Installation of the monitoring equipments on the Grohovo landslide (part of Croatian-Japanese bilateral scientific research project "Risk identification and Land-Use Planning for Disaster Mitigation of Landslides and Floods in Croatia") started in May 2011.The necessary measurement and research equipment, systems and equipment for meteorological, hydrological and geotechnical observations include meteorological stations and meteorological radar, Mini Diver instruments for measuring water and groundwater levels and instruments for geodetic and geotechnical monitoring.Geodetic monitoring includes geodetic surveys with a robotic total station -measuring 25 geodetic benchmarks (prisms) and GPS master unit with 9 GPS receivers (rovers) -while geotechnical monitoring includes vertical inclinometers, long-and short-span wire extensometers, pore pressure gauges and a weather station.A more detailed description of the monitoring system is given in Arbanas et al. (2014).For the portion of the research activities, a complex, integrated, real-time monitoring system was installed on the Grohovo landslide (Mihalić and Arbanas, 2013).
Soil parameters used in the computational simulation are presented in Table 1 and were determined from the undrained cyclic loading ring shear test and some from older laboratory testing (Benac et al., 2005).The significant geological explorations and measurements at the Grohovo landslide only started at the end of the 20th century.For this reason, authors used the data from Table 1 as the only hitherto relevant geotechnical data for the purpose of creating the numerical model.Although the coefficients k and B ss obtained from the undrained cyclic loading ring shear test were not used in the simulation, the relevant parameter used in the simulation is the excess pore pressure r u .
Long-term rainfall events and the consequent rise of the groundwater level have been the primary triggering factors for landslide occurrences in the Rječina River valley in the past (Vivoda et al., 2012;Žic et al., 2014).This increase in the groundwater level in the model was expressed by the pore pressure ratio values greater than r u = 0.60; the value r u = 0.60 corresponds to a groundwater level at the terrain surface.

Analyses and results
The depth-integrated numerical model SPH has already been benchmarked with problems where the analytical solution exists, such as a depth-integrated solution of a dam collapse across wet or dry channel bottoms.For the rheological model, comparisons can be made only using simple fluids whose rheological properties were obtained in the laboratory.A common solution to validate rheological models is to use numerical models (here the SPH method), implemented as an approximate mathematical model (here the depth-integrated model) and a rheological model, and recalculate observations from past events.
In this case study, the density of the mixture used was 2100 kg m −3 .The rheological models used to simulate this mudflow are the Newton fluid in turbulent regime model (Simulation 1) and the Real Bingham fluid model (Simulation 2) (Pastor, 2007).The parameters found to best fit the reconstructed event from 1908 were the turbulence coefficient value of 200-500 m s −2 , the friction angle of approx.27 • (tan ϕ = 0.466) and zero cohesion.Cohesion at the sliding surface during motion (c m ) 0.0 kPa Benac et al. (2005) Excess pore pressure (r u ) 0.0-0.6Assumption based on laboratory tests of soil samples taken in the area of the Grohovo landslide.By varying the values of the turbulence coefficient within the range of 200-500 m s −2 numerical simulations in the SPH model, no significant change was noted in terms of the mudflow reach, velocity and height of a mudflow, whereas an increased value of the turbulence coefficient 1000 m s −2 led to significant changes of the velocity and depth of mudflow.Subsequent analyses were performed assuming a rheological model with properties ranging within the values given in Table 1.All of the results below were obtained using this set of parameters and several preliminary simulations were executed assuming fully saturated soil.The erosion processes are modeled using the Egashira (Simulation 1) and Hungr (Simulation 2) laws with the following parameters: the sediment concentration of the flow, c = 0.64; the bed sediment concentration, c * = 0.7; and the empirical constant, K = 0.012.As expected, the results demonstrate that erosion processes seem to be strongly dependent on the channel slope.
The first SPH simulation (Simulation 1) comprised Newton's model of the turbulent flow regime with the effect of erosion activity (Figs. 9 and 10).The mudflow propagation of the Rječina River was realized using the shallow-water module (SW module) of the SPH code (Pastor, 2007).Most input parameters for the simulations are presented in Table 1.The spatial domain was discretized with an equidistant mesh with a size of 5 m, resulting in 128 453 nodes.The initiating mass of the naturally formed dam was created with 132 nodes.Each node was given an initial height of the ma-terial.The acceleration of gravity g is taken with the value of 9.81 m s −2 , fluid density ρ with 2100 kg m −3 , the Manning coefficient of roughness is n = 0.04 m −1/3 s, the friction angle during motion is 26 • and the minimum thickness of the layer under shear stress due to flow is assumed to be 0.001 m.Within the SPH code the control parameter for the pore pressure ic pw = 1 was set to 1.0, to account for the reduction of the pore water pressure.Parameter pw prel = 0.6 was adopted as the ratio between the pore pressure P and the liquefaction pressure (pw prel = P /P licuef ).The time increment in the calculation was taken as 1 s.The intention was to provide a simulation of the mudflow propagation along the Rječina River resulting from the formation of muddy deposited materials downstream of the Grohovo landslide and its gradual saturation with groundwater at a level corresponding to the maximum elevation of the deposited materials (fully saturated materials).The overall runout distance of mudflow propagation for this simulation is approx.1745 m, which was reached after 236 s of the initial flow formation.The maximum flow velocity recorded in the simulation is about 20 m s −1 (72 km h −1 ), and the maximum affected area due to mudflow is 4.15 ha.The initial volume of muddy materials is 132 450 m 3 , whereas the final total volume of mudflow propagation is 427 550 m 3 .The total volume of mudflow propagation along the Rječina River is approx.295 100 m 3 .The maximum depth of mudflow-deposited materials is 30.7 m (in a canyon of the Rječina River, near the Pašac Bridge), whereas the minimum depth of deposited material is 10.9 m.In Figs. 10 and 12, the height variability values of the de-   (Pastor et al, 2009a;Blanc et al., 2011).As with the first SPH simulation, the runout distance of deposited materials, their flow velocity, the depth of the deposited materials and the size of the area affected by the mudflow propagation are recorded.In this model, the total mudflow propagation obtained from the simulation has a duration of approx.236 s.The maximum runout distance of the mudflow is 1992 m, which was reached after 221 s.The maximum flow velocity of mudflow propagation is about 21 m s −1 (approx.76 km h −1 ), whereas the maximum affected area due to mudflow propagation is around 4.53 ha.The initial volume of muddy materials is 132 450 m 3 , whereas the final total deposited volume is approx.462 122 m 3 .The difference between the above two volumes yields the total mudflow volume generated within the Rječina River due to the mudflow propagation of 329 672 m 3 .The maximum depth of the mudflow that occurs during its propagation is slightly less than 33 m (in the canyon of the Rječina River, directly upstream of the Pašac Bridge).

Discussion
It can be concluded that the simulations using the Hungr erosion law gave similar results for the deposition pattern, mud volume and the flow velocity as the simulations adopting the Egashira erosion law.The differences in results for the erosion processes are shown in Figs. 13 and 14.The volume of mudflow increases faster using the Egashira erosion law than using the Hungr erosion law, but the final volume using the Hungr erosion law is slightly higher.The Egashira law seems to be better suited to this case study than the Hungr law based on descriptive evidence from old historical documents.In contrast, in the simulation with the Hungr law, the linear erosion rate has a quite high value, which explains why the volume increases along the entire flow path.
The analysis of the erosion processes has shown quite significant oscillations (variations) in the erosion activity along the Rječina River.Indeed, the Egashira erosion law improves some characteristics of the mudflow: the flow velocity and mudflow deposition pattern (height of mud lobes) (Fig. 15).However, the results for the erosion rate and the increased volume are quite similar to those using the Hungr erosion law (Fig. 16).
The above analysis allows for a comparison of the effects due to the two erosion laws that are not based on the same parameters, as the Hungr erosion law is based on the flow velocity and the flow depth, whereas the Egashira law is based on the current velocity and the slope of the terrain.Both of these laws allow the initial volume of the mudflow to increase along the travel path to reach the same final mudflow volume as it happened in the actual event.However, the volume does not evolve in time in the same manner for the two laws.Using the Egashira law, the volume tends to vary more similarly than the real mudflow behavior, which is very roughly described in historical records found in the Croatian State Archive in Rijeka (Benac et al., 2006;Žic et al., 2014).Therefore, it can be argued that the Egashira law results seem to be more realistic than those using the Hungr erosion law.
The distinctive features of mudflow are strictly related to the mechanical and rheological properties of the involved materials, which are responsible for their long travel distances and the high velocities that they may attain.The numerical simulation is very sensitive to the choice of these pa- rameters.Runout predictions are affected by the initial mass and the rheology selected.Good estimates of the initial distribution of the pore pressure and pore pressure dissipation are required.Despite these uncertainties, the prediction of the runout distances and velocities through mathematical modeling of the propagation stage can notably reduce losses due to these phenomena by providing a means for defining hazardous areas, estimating the intensity of the hazard and identifying and designing appropriate protective measures.
Regarding discretization effects, the mudflow mass is discretized using a series of nodes (material points).The accuracy of the simulation greatly depends on the number of nodes.It is possible to perform simplified analyses with a reduced number of nodes.The results of the analyses showed that using a smaller number of material points had a greater effect on the velocity rather than on the flow path; therefore, a smaller number of material points could be used for providing estimates.
Reliable forecast of susceptible propagation areas and the velocities of mudflows is a crucial issue for risk analysis, and the numerical modeling of the propagation stage is a valuable tool to predict these quantities in engineering analyses.However, the irregular topography of natural slopes considerably affects the motion of propagating materials, and accurate DEMs are paramount for realistic simulations and assessments (Delinger i Iverson, 2004;Cuomo et al., 2013Cuomo et al., , 2014)).
Several simulations were created with different spatial domain discretizations (equidistant 2 m × 2 m, 5 m × 5 m or 10 m × 10 m mesh grids) (Table 2).The simulation view of the mudflow propagation in the 10 m × 10 m case was very different from the 5 m × 5 m and 2 m × 2 m cases.In the 10 m × 10 m case, the flow occurred in multiple directions on the terrain and, in the end, the model was seen as too crude to provide a reliable mudflow simulation.
The velocity of the mudflow, its path and the runout distance depend greatly on the terrain topography.For SPH models, structured topographic meshes are more suitable because it is immediately possible to determine the cell to which a given point belongs.Therefore, a first indicator of the precision of the mesh is the product of the second-order derivative of the basal surface height by the square of the mesh size, but this is not sufficient.Based on our experience, it is suggested that at least 10 points should be used to discretize canyons and gullies channeling flow.
In addition to the DEM cell size, there are elements with characteristic sizes smaller than the DEM grid spacing that can affect the propagation path, such as cascades, bridges and large stone blocks that can divert the flow.Proper modeling of these features requires the inclusion of special elements in the analysis as these features may artificially divert the flow.To consider them, special barriers have been included in this study, composed of a series of nodes that interact with those of the flowing material whenever the distance between them is less than a given tolerance, which was here adopted as half the topographic grid size.One of the major practical issues in setting up the simulation was the choice of a particular rheological model and its parameters.Cohesive fluid models, such as Bingham, are recommended for modeling mudflows.Mudflows are usually generated in very loose metastable materials, where the pore pressures generated in the triggering process have largely contributed to the failure, closely associated with the ground-water level in the soil.High groundwater levels (significantly saturated soil) cause sudden launches of muddy materials, resulting in significant propagation velocity at the start and propagation of larger amounts of material downstream.Additionally, the grain size and density of the material and the ratio of the lateral pressure have a great effect on the sen-  sitivity of the numerical model and the propagation of the mudflow (Fig. 17).

Conclusions
Computational simulation using a coupled, SPH depthintegrated model capable of considering pore water pressure dissipation in the mudflow mass was presented.The propagation of the catastrophic mudflow that occurred in the Rječina River valley (Croatia) in 1908 was simulated.The validity of the proposed approach was assessed using two rheological models and two erosional laws.In the first simulation, Newton's model was applied to the turbulent regime, whereas the second simulation considered the propagation of mudflow based on the Real Bingham fluid model.The obtained results highlight the capability of the SPH framework to simulate the propagation stage of such complex phenomena and the relevant role played by the rheological properties in an adequate simulation of the runout distance, velocity, affected area and height of the propagating masses.From the results of these simulations, it can be concluded that the Real Bingham fluid model is better suited to modeling real mudflow propagation from the given input hydrogeological parameters.The objective of this study was to apply and validate the SPH 2-D integrated model on a real terrain configuration and on a real event from the past in order to facilitate simulations that can be used in engineering practice, including the Hungr and Egashira erosion laws.The study suggests that the use of the Egashira erosion law yields better predictions for the velocity and the deposition samples than the use of the Hungr erosion law.However, both of these erosion laws give a good estimate of the final volume.
Due to the very scarce data about the mudflow occurrence that occurred in 1908 in the area near the Grohovo village, the verification of the described model has been limited.It should be noted that a part of the numerical simulation was qualitatively verified on the basis of old historical images, from which the height of mudflow in some places within the Rječina watercourse was reconstructed.The historical pictures of events are in black and white, which complicated the verification.The mudflow occurred very rapidly and no actual measurements were recorded.The heavy precipitation that occurred during and after the event have further hampered any chance of thickness measurements of the suspended sediment, as the fine-grained material was easily flushed away.From the technical records in the old documents it is suggested that the mudflow propagation did not reach the mills in Žakalj village (see Figs. 8,9 and 11), which on that occasion was not damaged.Compared to the citations and statements within the Hungarian project of the river reg-ulation of the central part of the Rječina River catchment area (Žic et al., 2014) it can be the concluded that the presented simulation of mudflow propagation represents a reasonable reconstruction of the actual event.
The considered erosion laws should be further examined in a hydraulics laboratory using the hydraulic flume.The adopted simulation can be applied to other mudflow events from the past to create a database necessary for the calibration and loading to a valuable database of specific parameters.
Based on the presented computational simulations, it can be concluded that the potential mudflow propagation is unlikely to threaten the urban part of the city of Rijeka and that it is unlikely to cause substantial effects on the environment or lead to loss of human lives.

Figure 1 .
Figure 1.Map of landslides on the wider area around the Grohovo village.

Figure 7 .
Figure 7. Reference system and notation used in the numerical modeling.
. The above equation represents the quasi-Lagrangian form of the vertically integrated www.nat-hazards-earth-syst-sci.net/15etal.: A model of mudflow propagation downstream from the Grohovo landslide near Rijeka

Figure 8 .
Figure 8. Digital elevation model of the Rječina River valley with an SPH mesh, showing natural dam materials and the Valići accumulation.
Several numerical simulations were conducted to describe the mudflow propagation based on the Bingham rheological model (Simulation 1) in which the turbulence coefficient varied in the range of 200-1000 m s −2 , while the angle of internal friction was www.nat-hazards-earth-syst-sci.net/15/293/2015/Nat.Hazards Earth Syst.Sci., 15, 293-313, 2015

Figure 9 .Figure 10 .
Figure 9.The simulation view of mudflow propagation on the Rječina River (based on the equidistant mesh grids of 5 m × 5 m), with the propagation of materials from a natural dam formed on the Grohovo landslide, Simulation 1.

Figure 12 .
Figure 12.The results of the mudflow propagation on the Rječina River (based on the equidistant mesh grids of 5 m × 5 m), with propagation of materials from a natural dam formed on the Grohovo landslide, Simulation 2.

Figure 13 .
Figure 13.The simulation view of erosion activity on the Rječina watercourse with the application of the Egashira erosion law, Simulation 1.

Figure 14 .
Figure 14.The simulation view of erosion activity on the Rječina watercourse with the application of the Hungr erosion law, Simulation 2.

Figure 15 .
Figure 15.Visualization of dependence of individual output parameters with respect to time and runout distance of mudflow propagation: Simulation 1 with the application of the Egashira erosion law.

Figure 16 .
Figure 16.Visualization of dependence of individual output parameters with respect to time and runout distance of mudflow propagation: Simulation 2 with the application of the Hungr erosion law.
Figure 17.(a) Comparison of the linear erosion rate and the relative traveled distance and (b) comparison of the volume increase rate and the relative traveled distance for the 1908 Grohovo mudflow event.

Table 1 .
Soil parameters used in the SPH computer simulation.

Table 2 .
The impact of spatial domain discretization on output parameters of mudflow propagation with the application of the Egashira erosion law.