Numerical Modeling of Flow Patterns Applied to Analysis of Susceptibility to Movements of the Ground

: Mass movements in deformed areas of natural relief deformed by seismotectonic factors are one of the most destructive and recurrent natural hazards in the Republic of Ecuador, especially during intense rain periods, the El Niño phenomenon, or due to earthquakes such as the one that occurred on 16 April 2016 in the Ecuadorian coastline. This study proposes the application of Hydrological Model D8 and its derived morphometric parameters like slope, orientation of the slope, and curvatures, extracted from the high spatial resolution Digital Elevation Model (DEM), implemented in programs such as Rockworks 7 (gridzo), SURFER (downwards slope), ArcView (ﬂowacc), and SAGA (curvatures) to obtain runoff ﬂow, structural geological lineaments, and superﬁcial deformations of the topographic relief that are the origin of erosion, superﬁcial landslides, lateral propagation, of the rock–soil complex, mass ﬂows, and deep gravitational deformations. This methodology has been validated in three locations with intense deformations: two in Ecuador and one in Spain. The DEM were obtained from the Ecuadorian Spatial Institute (ESI) (spatial resolution of 10 m), the Rural Technological Infrastructure and Information National System (SIGTIERRAS) (spatial resolution of 5 m), and the Council of Andalusia (spatial resolution of 5 m). Conceptualization, M.C. and A.M.; Methodology, M.C. and A.M.; Software, M.C.; Validation, M.C. and A.M.; Formal Analysis, M.C.; Investigation, M.C. and A.M.; Resources, M.C.; Curation, M.C.; Writing-Original Draft Preparation, M.C. and A.M.; Writing-Review & Editing, A.M.; Visualization, Supervision: A.M.; Project Administration, A.M.; Funding Acquisition, A.M.


Introduction
Mass movements in Ecuador are one of the biggest threats that cause recurring disasters, especially during the periods of intense anomalous rains, the El Niño phenomenon, and earthquakes such as the one that occurred on April 16, 2016 in the Ecuadorian coastline.
In 2010, with the creation of the National Secretariat for Risk Management (SNGR), a standardized methodology was developed to elaborate mass movements' susceptibility maps at a municipal level in the entire country using a 1:50,000 scale. This methodology was based on heuristic techniques, taking in consideration a number of variables based on the experience of researchers that have studied these types of risks over a period of several years, and their respective support documents, such as: Generation of Geoinformation for the Valoration of Territory at a National Level, component 3 Geopedology and Geological Risks project developed by the National Secretariat of Planification (SENPLADES); Initial Methodology for the elaboration of preliminary susceptibility maps to mass movements at a 1:50,000 scale, prepared by SNGR-SUN MONTAIN; and Mass Movements in the Andean Region: A Guide for Risk Evaluation, from the Multinational Andean Project (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008).
This cartography identified susceptibility areas as intensity spots in a range going from absence of susceptibility to very high susceptibility and related these susceptible areas as regional geological structures; however, it wasn't able to define structural lineaments, fractures, or topographic relief methodology to extract guidelines or structural geological lineaments and relief deformations has quick processing time, high precision, and it uses free digital elevation models. In the following images, the structural guidelines are shown with the two methodologies applied to the same area (Nicaragua) with the same 20 m spatial resolution digital elevation model. In Figure 1A, the results of the coincidence between the guidelines obtained by remote sensing are shown. Figure 1B presents the structural guidelines obtained by the proposed method of this investigation. Figure 1C shows the coincidences and advantages of the proposed method; in addition to extracting structural guidelines and higher density fractures, relief deformations are also extracted by concentration of these lineaments, as can be evidenced in Figure 1D.

Materials and Methods
In order to extract the flow-path direction and its layout, structural geological guidelines, fractures, and superficial relief deformations, the theory of the basic Hydrological Model 8D was applied with high-resolution digital elevation models obtained from specialized institutions such as Telespazio (radar images), Cosmo SkyMed (10 m spatial resolution, Zone 1), Sigtierras (radar images, 3 m spatial resolution, Zone 2), and the Council of Andalusia (5 m spatial resolution, Zone 3). The used DEM were corrected before analysis. The used coordinate systems were UTM and WGS84-Zone 17S for the Ecuadorian studied areas, and UTM-Zone 30 and ETRS89 in Spain.

Materials and Methods
In order to extract the flow-path direction and its layout, structural geological guidelines, fractures, and superficial relief deformations, the theory of the basic Hydrological Model 8D was applied with high-resolution digital elevation models obtained from specialized institutions such as Telespazio (radar images), Cosmo SkyMed (10 m spatial resolution, Zone 1), Sigtierras (radar images, 3 m spatial resolution, Zone 2), and the Council of Andalusia (5 m spatial resolution, Zone 3). The used DEM were corrected before analysis. The used coordinate systems were UTM and WGS84-Zone 17S for the Ecuadorian studied areas, and UTM-Zone 30 and ETRS89 in Spain.
The used DEM were processed with public-use programs such as SAGA [18], ArcView [19] and commercial software like Rockworks 7 [20] and Surfer [21] to obtain the superficial water-runoff flow paths, structural geological lineaments, morphologic rupture lines, and deformations of the topographic relief. The realized fieldwork (300 control points) allowed us to verify the response of this methodology and collect data in susceptible sectors as well as types of associated mass movements. The operation of the basic Hydrological Model D8 and the determination of the direction and tracing of the surface runoff flows from precipitation are explained in the process diagram of Figure 2.
Geosciences 2018, 8, x; doi: FOR PEER REVIEW www.mdpi.com/journal/geosciences paths, structural geological lineaments, morphologic rupture lines, and deformations of the topographic relief. The realized fieldwork (300 control points) allowed us to verify the response of this methodology and collect data in susceptible sectors as well as types of associated mass movements. The operation of the basic Hydrological Model D8 and the determination of the direction and tracing of the surface runoff flows from precipitation are explained in the process diagram of Figure 2. The first derived products of the application of the D8 Model in the DEM were the direction and tracing of the runoff flow paths from the rain surface, which are the base of the information needed to obtain the following products: structural geological lineaments, fractures, and surface relief deformations of erosive and seismotectonic origin.
The calculation of flow directions and line tracing represents that each cell of the digital elevation model is the starting point to find hydrological characteristics in the topographic relief. From the topographic surface, a raster was generated that showed the flow directions and the exit direction of each line of each cell is drawn. Each valid exit direction is related with each adjacent cell where the flow can travel. This approach is commonly known as the eight-way flow (D8) [1][2][3].
Flow direction is determined by the steepest-descent direction or the maximum fall from each cell, and it is calculated as follows (Equation (1)): The first derived products of the application of the D8 Model in the DEM were the direction and tracing of the runoff flow paths from the rain surface, which are the base of the information needed to obtain the following products: structural geological lineaments, fractures, and surface relief deformations of erosive and seismotectonic origin.
The calculation of flow directions and line tracing represents that each cell of the digital elevation model is the starting point to find hydrological characteristics in the topographic relief. From the topographic surface, a raster was generated that showed the flow directions and the exit direction of each line of each cell is drawn. Each valid exit direction is related with each adjacent cell where the flow can travel. This approach is commonly known as the eight-way flow (D8) [1][2][3].
Flow direction is determined by the steepest-descent direction or the maximum fall from each cell, and it is calculated as follows (Equation (1)): where P: maximum fall, Z0−Z1: change of elevation value, and d: distance between each cell center. As distance is calculated between cell centers, if the cell size is 1, the distance between two orthogonal cells is 1 (creating an isosceles triangle). The distance between these two diagonal cells if √ 2 = 1.414. It may occur in this process that, when the maximum descent of various cells is the same, nearby sectors grow until the steepest descent and the flow direction of the cell are codified with the value that represents the mentioned direction. If all next elements are higher than the processing cell, it is considered a sink that fills with the lower value of the following elements and has a flow direction towards the mentioned cell. However, if a cell is next to the physical edge of the DEM or has at least a NoData cell as a neighbor, it doesn't fill or doesn't have value because the neighbor's data are insufficient to be considered a true cell and a true flow and cell receptor. In order to be a receptor, the information around the cell has to be present to continue the process [22].
This procedure was applied to each cell of the DEM, Cn (I,j), with the exception of the borders and their 8 neighbors. Thus, water that enters a cell can exit through one of the 8 mentioned directions.
Workflow starts from a C0 cell (painted in yellow) with height starting point Z0. For each cell from C1 to C8 in a mesh, of designated heights by the DEM, slope P is calculated for each one of its neighbors to the north, northeast, east, southeast, south, southwest, west and northwest, as it is shown in Figure 3A. The objective of this process is to calculate the tangent of the angle formed by the horizontal with the line that joins the point under study with the neighbor that is the direction of the flow ( Figure 3). The process is completed using Equation (2) ( Figure 3B): Geosciences 2018, 8, x FOR PEER REVIEW 6 of 18 In Figure 3C, you can see how the maximum-slope line is drawn (maximum fall) inside the DEM. In a 10 m spatial resolution DEM, represented by the grid in Figure 3D, the lowest point is 750 masl. The construction of a morphologic flow and a rupture line from an initial C0 cell is carried In Figure 3C, you can see how the maximum-slope line is drawn (maximum fall) inside the DEM. In a 10 m spatial resolution DEM, represented by the grid in Figure 3D, the lowest point is 750 masl. The construction of a morphologic flow and a rupture line from an initial C0 cell is carried out in 4 phases: In Phase 1, the initial point of the line is established in a time t = 0, with the initial state C0 = 1, and assigned height by the DEM = 1068.44 m. Figure 3E,F shows the configuration of the data and their spatial position in the E, SE, NE, N, NW, W, or SW directions from the starting cell C0 (i,j).
In Phase 2, the slopes are calculated towards the 8 nearest neighbors, from the C0 cell towards the C1, C2, C3, C4, C5, C6, C7 and C8 cells, according to the local slope obtained from the slope map calculated with Equation (1), where d is the horizontal distance between the cells, which is 10 m in this study. This calculation was completed in the base of the previous equations. Orthogonal distance d from orthogonal cell centers is 1, and the distance between two diagonal cells is √ 2 = 1.414 (Equation (3)): The spatial distribution of the slopes can be seen in Figure 3F. Three cases are possible in the drawing of flow lines, flow patterns, or flow ( Figure 4A-C). In Case 1, cell C1 is surrounded by 8 cells higher than it. In this case, Cell C0 is a "pit" in point B ( Figure 4A). In Case 2, cell C1 is part of a group of cells with the same dimension that represents flat areas. In the event that the flat area is surrounded by higher cells, it is known as a sink ( Figure 4B). In Case 3, more than one cell are adjacent to C1 with a slope equal to the maximum. This problem can be called an indeterminacy, and it occurs when C1 is on a coastal edge or DEM boundary. However, if the DEM has abrupt discontinuities or inflection zones, these "indeterminate edges" within the DEM can determine the presence of regional tectonic structural systems such as faults, fractures, and terrain deformations ( Figure 4C).  Figure 4A). In Case 2, cell C1 is part of a group of cells with the same dimension that represents flat areas. In the event that the flat area is surrounded by higher cells, it is known as a sink ( Figure 4B). In Case 3, more than one cell are adjacent to C1 with a slope equal to the maximum. This problem can be called an indeterminacy, and it occurs when C1 is on a coastal edge or DEM boundary. However, if the DEM has abrupt discontinuities or inflection zones, these "indeterminate edges" within the DEM can determine the presence of regional tectonic structural systems such as faults, fractures, and terrain deformations ( Figure 4C).  In Phase 3, the propagation rule is applied. The cell with the maximum negative slope defines the cell that will change state. In our example, the cell is C8 = −89.28; therefore, C8 is the cell that will light up (yellow in Figure 3F), so that the path of the flow line grows on the matrix, and its length depends on the iterations necessary to reach some of the completion conditions, thus creating the system of flow paths. The process is repeated from the new cell state with a negative maximum slope in an iterative process throughout the study grid ( Figure 5A1,A2).
The chosen point (maximum slope) was incorporated in the downstream flow line and used as a basis to return to Phase 2, with negative and positive lateral-flow transport zones, flow accumulation zones, and flat areas. In Figure 5B1, the center of cell C0 is the start of the flow lines. In Figure 5B2, the drawing of the flow "downstream" line is displayed; and in Figure 5B3, the final plot of the flow line upstream and downstream from the starting cell C0 is depicted. Figure 5C1,C2 demonstrate the concentration downstream of the runoff flows that form channels and drainage networks that collect the waters of the entire basin and pour them into a single river at the mouth of the basin.  Figures 5C1 and 5C2 demonstrate the concentration downstream of the runoff flows that form channels and drainage networks that collect the waters of the entire basin and pour them into a single river at the mouth of the basin. In Phase 4, the plane and profile curvatures in the matrix of the slopes of Figure 3A are calculated with Equations (4)-(7) below [23]: Profile Curvature : PLC = 2 DG 2 + EH 2 + FGH In Phase 4, the flat and profile curvatures are calculated, which are one of the factors that contributes to relief deformation and instability. These curvatures were extracted directly from the DEM, apart from aspect and slope. The layouts of the flat curvature (concavity degree/convexity along an isoline perpendicular to the slope profile) and the profile curvature (concavity degree/convexity of the slope profile) were selected due to their influence over the water content that is stored in the soil and is responsible for mass movements when intense precipitations occur. These traced curvatures, such as continuous isolines, brake slopes and orientation originating natural fractures processes [24].
For the flat curvature, negative values show that the surface is concave upwards and positive values show that the surface is convex upwards. For the profile curvature, negative values show that the surface is convex upwards and positive values show that the surface if convex upwards.
A cero-curvature value shows that the surface is flat for both types of curvature. Additionally, for each DEM cell, from the C0 cell, if the line that represents flow direction or curvatures is infinitely grown until the DEM edge is located, the structural lineaments (parallel lines to the curvatures) can be found, and deformation zones of the topographic relief in the intersections of these lines with the runoff flow paths from superficial water can be also found. The intensity of these concentrations values the deformation and susceptibility degree to mass movements in the study area. The curvature values are calculated in all cells of the DEM and intermediate meshes of slope change are drawn, obtaining the curvature vector (inflection that indicates concavity or convexity from node to node in degrees. Through this second operation, continuous facets of similar orientation within the DEM are connected via a direction vector from the cell centers. Vectors are where surface runoff water enters from the upper parts and erodes the deformed and weathered geoforms (morphological break line) ( Figure 6A).
Within the digital elevation model, when the runoff and rupture lines were concentrated, erosion channels and areas of subsidence and breakage of the geoform were generated. Climate and topographic relief influence the pattern of a network, but the underlying geological structure is usually the most influential factor for the generation of this type of line and erosive processes. Hydrographic patterns are closely related to geology, which can help identify faults and structural lineaments [23,[25][26][27].
The final active surface deformation of a region (for example, any DEM) with internal mass flows, can be spatially determined by the intersection of a set of morphological break lines that define critical points and regions of movement susceptibility ( Figure 6B1-6B3). By applying a density filter to these lines, we obtained Morphological Sliding Plates ( Figure 6C1-6C3), which relate morphology and geological structures, within which can be found fault escarpments, triangular facets of erosion, closing hills, pressure hills, deformation of alluvial cones, inversions of river beds, and other morphotectonic aspects associated with types of mass movements, such as rotational, translational, lateral expansion, falling blocks, flows in general, solifluction, and creep.
This methodology was validated in two areas of Ecuador ( Figure 6D1) and one in Spain ( Figure 6D2). Zone 1 is located between the towns of El Cadial, Chilcaloma, and Angas in the province of Bolívar, Ecuador. Zone 2 is located in the town of Pedernales, province of Manabí [28], Ecuador, and Zone 3 is in Arenas Del Rey in Granada, Andalusia, Spain [29]. Bolívar, Ecuador. Zone 2 is located in the town of Pedernales, province of Manabí [28], Ecuador, and Zone 3 is in Arenas Del Rey in Granada, Andalusia, Spain [29].

Results
Zone 1: in all of Zone 1, structural guidelines, fractures, surface deformations and their type of associated mass movements, verified with field work, were found in three great plates of intense deformation, framed in red-colored polygons, ( Figure 7A). Zone 1 in El Cadial-La Chorrera-La Chorrera de Arriba (340 ha); Zone 2 in Chiriyacu (262 ha); and Zone 3 in Agas-Saltuco-Chilcaloma-El Salto (370 ha). Besides that, other smaller deformation plates were found in Las Juntas and Santa Lucia. The lower zones of the river valleys, such as in Balzapamba, with an intense deformation with blue tendencies, represent mass-flow areas that generally affect these zones during intense precipitation times. For this zone, the digital elevation model of 10 m of spatial resolution was enough to obtain in detail the surface-relief deformations, structural lineaments and their relation with mass flow such as debris flow in ravines, eroded hillsides, and translational landslides ( Figure 7B).

Results
Zone 1: in all of Zone 1, structural guidelines, fractures, surface deformations and their type of associated mass movements, verified with field work, were found in three great plates of intense deformation, framed in red-colored polygons, ( Figure 7A). Zone 1 in El Cadial-La Chorrera-La Chorrera de Arriba (340 ha); Zone 2 in Chiriyacu (262 ha); and Zone 3 in Agas-Saltuco-Chilcaloma-El Salto (370 ha). Besides that, other smaller deformation plates were found in Las Juntas and Santa Lucia. The lower zones of the river valleys, such as in Balzapamba, with an intense deformation with blue tendencies, represent mass-flow areas that generally affect these zones during intense precipitation times. For this zone, the digital elevation model of 10 m of spatial resolution was enough to obtain in detail the surface-relief deformations, structural lineaments and their relation with mass flow such as debris flow in ravines, eroded hillsides, and translational landslides ( Figure 7B). In Figure 7C-F, real photographs of the deformed hillsides were compared, without the application of the methodology proposed in Figure 7C,E, and with the application of the methodology proposed in Figure 7D,F, where intense deformations were obtained in these hillsides and their mass-movement systems associated with translational landslides, falling blocks, and mass flows. It also allowed for better comprehension of the relation between deformation of the topographic relief, mass movements, and the vulnerability of the population and its agroproductive systems to this kind of threat. An important characteristic of the obtained result with 10 m of spatial resolution is that the hidden geomorphological characteristics can be discovered as circular structures that can apparently be collapsed volcanoes, ancient megalandslides, and surface fractures caused by the location of intrusive locks (intrusive of Balsapamba). This circular structure, to the east of the Angas, is shown in a series of mass flows ( Figure 7G).
Zone 2: This zone is located on the town of Pedernales, province of Manabi, in the Ecuadorean coast. Pedernales' population was the most affected by being near the epicenter of the earthquake of April 16, 2016, with a 7.8 score on the Richter Scale ( Figure 8A). In Figure 8B we can observe the intensity of the concentration of structural geological lineaments and fractures (red-colored stripes) in the Zero Zone of Pedernales, obtained with a 3 m spatial resolution DEM before the earthquake. In the Zero Zone, the population suffered major destruction and building loss, as well as the loss of human lives (more than 183). In Figure 8C, we can observe destroyed buildings by the earthquakes, delimited by red-colored polygons. In Figure 8D, a detail of the structural lineaments that possibly amplified the superficial Rayleigh waves that destroyed houses and civil infrastructure is shown. In Figure 8E, we can observe a detail of the destroyed or severely damaged houses, marked with a red-colored cross and revised in the field after the earthquake. This infrastructure was raised (tendency to the yellowish-red color) or lowered (tendency to a violet-red color) as it is shown in the topographic relief image, obtained with an HSV saturation filter (saturation matrix). These houses were found near or on top of structural lineaments, surface fractures, and deformed zones obtained with the presented methodology. The houses with green colors show no destruction or cracks.
Zone 3: The third zone that theoretically validated this proposed methodology without field-work verification was in Arenas del Rey, Granada, Andalusia, España ( Figure 9A). According to the chroniclers of the time, it was moved by an earthquake occurring on 25 December 1884, with its epicenter in Arenas del Rey, and a magnitude between 6.2 and 6.5 points in the Richter Scale, lasting for approximately 20 seconds, with a hypocenter between 40 and 50 km underground, causing between 1050 and 1200 deaths and approximately double that number in injuries. A total of 90% of civil structures were destroyed. Figure 9B and 9C shows surface structural geological lineaments that exist today in Arenas del Rey obtained using only the proposed methodology. Figure 9C shows the limits of the city (solid red-colored lines) and its location on top or near geological structural lineaments, with NS, NE-SW, and NW-SE directions, and surface fractures that possibly amplified the seismic waves of the 1884 earthquake and destroyed the city. This situation puts the population at high risk. Figure 9D shows the zoning of the seismic vulnerability of city plots in function of the actual terrain deformations that could raise or lower civil structures with a similar earthquake to the 1884 one. Figure 9E shows the concentrations (tendency to red color) of the lineaments and structural fractures in the Federico Garcia Lorca church, Alfonso XII plaza, and Jaen and Infanta Eulalia streets, which are one unified deformation train. These old deformations were possibly activated with the 1884 earthquake and destroyed the church, plazas, and other nearby buildings (90% of the population was killed). In the attached images ( Figure 9F) of the document "The Earthquake of Alhama of Granada of 1884 and its Impact", the destruction of the church and the houses near the Alfonso XII plaza are represented. This evidence corroborates the actual deformation of those plots that place these buildings in danger.

Discussion
This standardized methodology by the SNGR, using heuristic and variable weighing to elaborate susceptibility to mass-movement maps, doesn't identify geological structural lineaments, morphological ruptures lines, fractures, and surface-relief deformations. It also does not identify different types of mass movements, subsidence, and specific subsidence in small areas. However, these maps are a base reference to incorporate in them the geomorphologic traces mentioned, obtained with the proposed methodology that uses flow directions and its vector line with magnitude and the direction it represents. It also traces morphologic rupture lines using slope values, orientation of the slope, and flat and longitudinal curvatures, in public-license software such as SAGA and Arcview, and in commercial software such as Surfer and Rockworks 7 used in this study. This geomorphological traces were obtained in a quick fashion, unlike with the heuristic methodology, weighted variables, and analog methodologies such as the Mora Varhson image analysis that need larger thematic coverage and a larger amount of time for their processing.
Digital elevation models of 10 m, 3 m, and 5 m of spatial resolution used were enough to obtain relief-deformation zones based on structural lineaments and relief fractures. These deformation zones possibly contain weak geomechanic characteristics, can weather easily, and generate mass movements by rain such as in Zone 1, or by earthquakes with more than 6 points in the Richter Scale, such as the ones that shook the Ecuadorean coastline (2016) and Arenas del Rey, Andalucia, Spain (1884).
This methodology can be used with DEM of different spatial resolutions with ranges from 1 m (Lidar) up to big resolutions obtained from SRTM, AsterGdem, ALOS-PALSAR, Open Topography, and GTOPO 30, among other servers.
The ease to obtain digital elevation models with their public application allowed for quick obtaining of results to delimit borders, instability borders, and safe areas that can be used to locate meeting points, space routes, and safe areas to mobilize the population in the event of this type of threat. The application of this methodology can strengthen the operation of early warning systems to this type of threat.

Conclusions
The zoning of susceptibility by mass movements based on flow directions with the basic D8 Hydrological Method uses digital elevation models processed in free and commercial software, facilitating the identification and quick localization of structural geological lineaments, fractures, relief-deformation areas, and different types of mass movements.
A special characteristic of the presented methodology to obtain structural lineaments, fractures, and surface-relief deformations is that, depending on the digital elevation model resolution, it is possible to extract hidden geomorphologic characteristics such as extinct volcanic calderas, relief deformation by intrusive rocks and fractured pressure mounts, and deformed hillsides, among other geomorphologic traits.
The location of deformed areas by concentration of structural lineaments and relief fractures obtained from this methodology, with the field verification done in Ecuador (300 field control points) and in the base of chroniclers in Spain, have corroborated what happened in the areas of study.
Finally we conclude that the proposed methodology, unlike the standardized and analog methods, is economic, with a high precision degree, and can be applied to any geographic region before a disaster occurs caused by natural risks, such as intense rains or earthquakes with a ranking higher than 6 on the Richter scale. Thus, this methodology can be used as a forecasting mechanism through which maps of structural lineaments, fractures, and quaternary deformations with their own associated mass-movement types can be created, so populations located on or near these geomorphologic traits can be evaluated by their vulnerability degree to subsidence, structural collapse, landslides, and other types of mass movements. This will allow the elaboration of prevision and reduction maps to this type of risk in very high susceptibility zones and will strengthen municipal planning of the territory as well as its risk-management units, improving response and reduction protocols.