TerraceM-2: A Matlab® Interface for Mapping and Modeling Marine and Lacustrine Terraces

The morphology of marine and lacustrine terraces has been largely used to measure past sea- and lake-level positions and estimate vertical deformation in a wealth of studies focused on climate and tectonic processes. To obtain accurate morphometric assessments of terrace morphology we present TerraceM-2, an improved version of our Matlab R (cid:13) graphic-user interface that provides new methodologies for morphometric analyses as well as landscape evolution and fault-dislocation modeling. The new version includes novel routines to map the elevation and spatial distribution of terraces, to model their formation and evolution, and to estimate fault-slip rates from terrace deformation patterns. TerraceM-2 has signiﬁcantly improves its processing speed and mapping capabilities, and includes separate functions for developing customized workﬂows beyond the graphic-user interface. We illustrate these new mapping and modeling capabilities with three examples: mapping lacustrine shorelines in the Dead Sea to estimate deformation across the Dead Sea Fault, landscape evolution modeling to estimate a history of uplift rates in southern Peru, and dislocation modeling of deformed marine terraces in California. These examples also illustrate the need to use topographic data of different resolutions. The new modeling and mapping routines of TerraceM-2 highlight the advantages of an integrated joint mapping and modeling approach to improve the efﬁciency and precision of coastal terrace metrics in both marine and lacustrine environments.

The morphology of marine and lacustrine terraces has been largely used to measure past sea-and lake-level positions and estimate vertical deformation in a wealth of studies focused on climate and tectonic processes. To obtain accurate morphometric assessments of terrace morphology we present TerraceM-2, an improved version of our Matlab R graphic-user interface that provides new methodologies for morphometric analyses as well as landscape evolution and fault-dislocation modeling. The new version includes novel routines to map the elevation and spatial distribution of terraces, to model their formation and evolution, and to estimate fault-slip rates from terrace deformation patterns. TerraceM-2 has significantly improves its processing speed and mapping capabilities, and includes separate functions for developing customized workflows beyond the graphic-user interface. We illustrate these new mapping and modeling capabilities with three examples: mapping lacustrine shorelines in the Dead Sea to estimate deformation across the Dead Sea Fault, landscape evolution modeling to estimate a history of uplift rates in southern Peru, and dislocation modeling of deformed marine terraces in California. These examples also illustrate the need to use topographic data of different resolutions. The new modeling and mapping routines of TerraceM-2 highlight the advantages of an integrated joint mapping and modeling approach to improve the efficiency and precision of coastal terrace metrics in both marine and lacustrine environments.

INTRODUCTION
Marine and lacustrine wave-cut terraces are geomorphic markers largely used to estimate past sea-and lake-level positions (e.g., Zeuner, 1952;James et al., 1971), from which surface deformation rates may be assessed (e.g., Bradley and Griggs, 1976;Strecker et al., 1986;Bloom and Yonekura, 1990;Padgett et al., 2019). The increasing availability of highresolution topography has substantially contributed to the study of coastal geomorphic markers since the last decade (e.g., Palamara et al., 2007;Bowles and Cowgill, 2012). However, the identification and precise assessments of geomorphic markers using high-resolution topography requires dedicated morphometric analyses and data visualization techniques, which are in most cases not available in standard GIS platforms. TerraceM has been designed to close this gap by providing open-source tools to interactively map and model the formation of marine and lacustrine wave-cut terraces using high-resolution topography and Matlab R programming. So far, TerraceM has been used for diverse purposes in different coastal settings worldwide (e.g., Melnick, 2016;Benjelloun, 2017;Jara-Muñoz et al., 2017;Melnick et al., 2017Melnick et al., , 2019Pedoja et al., 2018;Bilbao-Lasa et al., 2019;de Gelder et al., 2019;Jara-Munþoz et al., 2019;Racano et al., 2019). These studies include research topics such as coastal dynamics, paleo-geography, paleo-climate, neotectonics, and earthquake geology, highlighting the broad spectrum of applications and the advantages of applying integrated joint mapping and modeling techniques in the Earth sciences.
TerraceM was originally designed as a Matlab R Graphical User Interface (GUI) to assess marine terrace elevations using high-resolution topography. The new version "TerraceM-2" goes further, including new modeling tools and processing functions inside user-friendly GUIs. Furthermore, TerraceM-2 is based on a novel programmatic architecture improving the computational efficiency and processing speed (see Table 1). Topotoolbox functions (Schwanghart and Scherler, 2014) are nested in TerraceM-2 to perform basic analyses, such as swath profile extraction, calculation between topographic parameters, and DEM visualization. The TerraceM-2 architecture is based on independent functions, whose outputs are stored as data structures, enabling memory efficiency, and processing greater number of swath profiles and larger DEMs, at least eleven times faster than the previous release (see speed test in SP1 and Supplementary Figure 1SP). The three new TerraceM-2 modules are devoted to mapping marine and lacustrine terraces, estimating metrics of past sea-level positions, exploring their formation and degradation mechanisms, and using these metrics to estimate deformation processes. The shoreline-angle mapping tool in the previous version (Jara-Muñoz et al., 2016) has been complemented with an interface to compute surface classification models, which allow mapping marine and lacustrine terraces based on the distribution of morphometric parameters extracted from the topography. Additionally, TerraceM-2 includes a new interface to model terraces formation using a Landscape Evolution Model, which simulates the effect of sea-erosion and vertical deformations on the development of wave-cut terraced topography. Furthermore, the spatial distribution of terrace elevations and surface uplift estimated from past sea-or lakelevel metrics can be compared with fault dislocation models to derive fault slip rates. The integration of TerraceM-2 interfaces and routines, in addition to user-friendly GUIs, optimizes the time required for mapping and modeling marine and lacustrine terraces for inexperienced users in Matlab R scripting. However, the TerraceM-2 architecture comprises independent functions allowing experienced users to create customized workflows. Taken together, we envisage that the evolving technological advances and the increased accessibility of high-resolution topography, in addition to the advanced morphometric and modeling methods in TerraceM-2, will promote studies on Quaternary sea-and lake-level changes and neotectonics. In light of the impacts of global sea-level change and seismic activity at plate-margin settings, such studies are crucial, particularly regarding the adequate assessment of hazards associated with seismicity and sea-level change, and a combination thereof.

CONCEPTUAL FRAMEWORK AND NEW MODULES IN TERRACEM-2
Here we present several concepts to define the morphogenetic characteristics of marine and lacustrine terraces, which are essential to understand the functionalities of TerraceM-2. For instance, the term coastal terraces refer to marine and lacustrine wave-cut terraces, generated by the effects wave erosion. In contrast to wave-built marine or lacustrine terraces formed by the accumulation of beach deposits (e.g., Dietz, 1963;Jara-Muñoz and Melnick, 2015), these terraces are purely erosional and may be covered by a thin veneer of sediments. The term coastal-terrace sequence denotes several coastal terraces that occur at successively higher positions above a reference datum, i.e., sea-level or lake level. The terraces form during rather high water levels and may form staircase-like sequences if the climatically controlled water-level changes are superposed to vertical crustal movements (e.g., Lajoie, 1986;Pedoja et al., 2011Pedoja et al., , 2014. There are several methods and markers that have been used to estimate the elevation of marine terraces; for instance, by identifying slope changes along topographic profiles (e.g., Armijo et al., 1996;Saillard et al., 2009;Regard et al., 2010;Roberts et al., 2013), or more sophisticated methods like semi-automated mapping of the inner edge (e.g., Palamara et al., 2007). However, using inadequate markers of terrace elevation may result in disparate errors. Therefore, the accurate mapping of marine terraces is crucial to obtain consistent elevation estimates and to derive unequivocal uplift rate estimates. The morphology of coastal terraces usually comprises a sub-horizontal paleo-platform bounded by a paleocliff inlands, the intersection between the paleo-cliff and the paleo-platform defines the shoreline angle, a geomorphic marker that represents the maximum reach of the sea or lake level during the formation of the terrace (Lajoie, 1986). Thus, the shoreline angle represents a geodetic marker that can be used to estimate relative water-level changes and vertical deformation (e.g., Lajoie, 1986;Burbank and Anderson, 2011;Jara-Muñoz et al., 2016). TerraceM-2 uses these geomorphic relationships and provides quantitative methodologies to analyze the metrics of coastal terraces by assessing their morphological elements and by modeling the topographic arrangement of coastal terrace sequences. TerraceM-2 comprises three modules, (a) Maptools (MAP) (b) Landscape evolution model (LEM), and (c) Dislocation Modeling (DIM), designed to be used both as independent and complementary GUI's ( Figure 1). Various ancillary tools are available for data analysis and interpretation as well, each module include an information interface that displays verbose operations to guide the user through the GUI's workflow, the World Catalog of Coastal Sequences (WCCS, Figure 1) comprising, wave-cut and wave-built marine terraces, coral reefs, and strandlines from the last high-stands (MIS 5, MIS 7, and MIS 9). This catalog is an update of the compilation of Pedoja et al. (2014Pedoja et al. ( , 2011, which currently includes 1010 sites based on 1232 references. The WCCS interface displays the location of marine terrace studies within a defined distance from the TerraceM-2 mapping area, allowing the user to rapidly integrate previous estimates of marine terrace elevations, terrace ages, dating methods, uplift rates, data quality index, and references. Another ancillary tool of TerraceM-2 is the Uplift Calculator that allows quick determination of uplift rates using both the typical method correlating terrace levels with sea-level high-stands (e.g., Lajoie, 1986), and the synchronous correlation method (Roberts et al., 2013) that considers the elevations of all terrace levels of a sequence simultaneously, by comparing observed and predicted marine terrace elevations based on the assumption of constant uplift rate along time. TerraceM-2 is designed to run in Mac R , Linux R , and Windows R computers using Matlab R 2017 or later (see Table 1), minimal RAM memory required is 2GB and more than 5 GB free disk space, which are also the minimal requirements to run Matlab R . The inputs and outputs of each TerraceM-2 module are based on data structures (Supplementary Figures 2SP, 3SP) allowing flexible organization and exchange between the different modules. TerraceM-2 inputs are: topography in raster format (geotiff or Arcgrid ascii), and polygon(s) in shapefile format that may be created in any GIS platform (for example QGIS or ArcGIS) containing rectangular boxes defining the location of swath profiles. Inputs are loaded and processed in the MAP module, and these results may be used as inputs for the LEM and DIM modules facilitating data exchange between the different modules. The outputs of each module may be exported to shapefiles and Matlab R structures (.mat format). This latter file also stores all the data that defines a particular project, allowing for subsequent editing and visualization. Furthermore, each module comprises independent functions, which can be used to

TerraceM-2 Maptools
TerraceM-2 MAP is a module designed for mapping coastal terraces (Figure 2). This module includes two complementary methods to measure the elevation and characterize the spatial distribution of coastal terraces.

Swath Profile Mapping
The swath profile mapping interface is an improved GUI that allows mapping and visualizing shoreline angles (Figure 3). TerraceM-2 uses the maximum elevations of swath profiles, which resemble the original coastal terrace morphology without the effects of post-abandonment modification, principally resulting from fluvial incision. As in the previous version of TerraceM, linear regressions are fitted upon user-defined segments of the paleo-cliff and paleo-platform of coastal terraces (Figure 3); the shoreline angle is then located by intersecting the extrapolations of linear regressions. Errors are determined by intersecting the 2σ confidence intervals of both regression lines. The swath profile mapping interface includes different methods (see TerraceM in Jara-Muñoz et al. (2016) for further details): (a) Staircase analysis, which is based on the intersection between paleo-platform and cliff, (b) the Fixed cliff analysis, which includes the option of adding field measurements to calibrate the paleo-cliff slope, and (c) the Stack analysis, to determine the shoreline angle in rough coastal areas comprising sea stacks, and (d) the Cliff diffusion analysis, determines the geometry of the colluvium wedge to reconstruct the original terrace profile as well as determining the geomorphic age of the terrace scarp (e.g., Hanks et al., 1984). Swath mapping can be also called from command line as an independent function (see SP2.1 and Supplementary Figures 4SP, 5SP). TerraceM-2 MAP includes ancillary visualization tools, such as Profile and Project Plotter, to create publication-quality plots of the mapping results.

Surface Classification Model
The second mapping method of TerraceM-2 MAP is based on Surface Classification Models (SCM) (Figure 4), designed to semi-automatically detect low-relief and gently inclined smooth  areas, often associated with coastal terrace paleo-platforms and steeper areas that may represent terrace paleo-cliffs (Eqs. 1 and 2). The SCM is based on the method developed by Bowles & Cowgill (2012) that uses a combination of topographic slope (SLP) and roughness (RG), the former being defined as the standard deviation of the slope (Frankel and Dolan, 2007).
The distribution of slope and roughness are truncated by userdefined cut-offs that characterize the morphometric attributes of paleo-platforms or paleo-cliffs (SLPr and RGr). The SCM is then calculated by normalizing and combining slope and roughness ranges, represented by patches masking marine terrace platforms or cliffs depending of the selected analysis (SCM cliff or SCM plat , Eqs. 1 and 2). False positive classifications may usually occur at valley bottoms and slopes; TerraceM-2 allows including a drainage network map and user-defined distance buffers to filter out such erroneous classifications (Figure 4). The resulting SCM patches are intersected with the topography to extract the elevations. Finally, elevation histograms are used to classify the SCM and to differentiate distinct terrace levels. The SCM algorithms are included in a GUI that allows for a quick adjustment and visualization of the morphometric parameters and interactive mapping of coastal terrace surfaces (Figure 4). The results of the SCM and swath profile mapping can be readily exported in shapefile, geotiff raster, or Googlearth R formats.

Architecture and Inputs
The application of LEMs in wave-dominated coastal systems have allowed studying diverse processes shaping the landscape, such as relative sea-level variations, vertical deformation, coral reef growth rates, sediment transport and deposition, cliff retreat and diffusion, eolian transport, and river erosion among others (e.g., Hanks et al., 1984;Anderson et al., 1999;Storms and Swift, 2003;Nakamura and Nakamori, 2007;Refice et al., 2012;Thébaudeau et al., 2013;Kline et al., 2014;Shikakura, 2014;Melnick, 2016;Husson et al., 2018;Limber and Barnard, 2018). The application of LEMs for modeling the development of sequences of coastal terraces sequences can provide valuable chronological information in the absence of terrace ages (e.g., Jara-Muñoz et al., 2017;Bilbao-Lasa et al., 2019). Furthermore, LEMs allow simulating a diversity of possible geomorphic scenarios that may be difficult to evaluate using other approaches; for instance, complex histories of marine terrace reoccupation (e.g., Melnick, 2016), or marine terraces affected by variable uplift rates along time (Racano et al., 2019).
The TerraceM-2 LEM has been adapted from the wave erosion and energy dissipation model initially formulated by Sunamura (1992) and further developed by Anderson et al. (1999). This model has been already used to simulate the development of drowned shorelines in southern Chile (Jara-Muñoz et al., 2017), rasa surfaces in northern Chile (Melnick, 2016) as well as staircase marine terrace sequences in Turkey (Racano et al., 2019). The model simulates the morphology of wave-cut terraces formed under oscillating sea level conditions and variable or constant uplift rates (Figure 5). The erosive power of wave abrasion is defined as a linear function of the rate of wave energy dissipation (E1, Eq. 3) that follows an exponential increase as depth decreases landwards (Sunamura, 1992).
where Sl is the sea-level position at time t and Zn is the local water depth, (dz/dt) represents the vertical seabed erosion rate, which depends of wave base depth (Wb, Eq. 4).
beyond Wb water particle motions are too low to erode the sea floor (e.g., Bradley, 1958;Sunamura, 1992), Wb is related to the breaking wave height (HB) by the dimensionless constant B (Eq. 4). This constant has been assessed in several studies (e.g., Hallermeier, 1978;Thornton and Guza, 1983;Kline et al., 2014) ranging between 0.1 and 0.5. The wave energy dissipates as the wave approach to the coast and depends of the water depth and the near shore platform morphology (Bradley, 1958;Komar, 1998). E0 is the initial energy in the far field, which resembles the erosive power of waves, and is assumed to be constant. When dissipation is lower than E0 the remaining energy in the system (Ec, Eq. 5) can drive cliff retreat leading to an expansion of the wave-cut platform (Sunamura, 1992;Anderson et al., 1999).
The initial model geometry is a linear synthetic slope on which the terraces are carved; the slope angle (φ) as well as the size and resolution of the modeling domain can be customized. The simulated topography is adjusted incrementally depending of the uplift rate, which can be constant or variable over time. Likewise, slope diffusion (Hanks et al., 1984) may be included to account for the effect of local climate and surface-process dynamics on topography, to reproduce realistic marine terrace morphologies. To account for sea-level variations, we included a compilation of 10 published sea-level curves (Imbrie et al., 1984;Hemleben et al., 1996;Rohling et al., 1998Rohling et al., , 2009Lambeck et al., 2002;Lea et al., 2002;Waelbroeck et al., 2002;Siddall et al., 2003;Bintanja et al., 2005;Spratt and Lisiecki, 2016). The amplitude of lowstand periods of the sea-level curves can be modified to reproduce the influence of local bathymetry, as for instance the sill in the Gulf of Corinth (e.g., de Gelder et al., 2019). In addition, the effects of the Holocene sea-level rise can be removed to preserve low-stand marine terraces. TerraceM-2 LEM provides various outputs including the evolution of the coastline position, time series of cliff erosion rate (Figures 6A,B), synthetic topographic profiles, synthetic shoreline angle elevations for each modeled terrace level, and chronology of the marine terrace sequences ( Figure 6C). The LEM outputs can be exported as image (.pdf), video (.mp4), and GIS files (.shp).

Comparing LEM Results With Marine Terrace Topography
The synthetic topography generated by TerraceM-2 LEM can be compared with measured swath profiles to determine ages of undated coastal terrace levels and estimate the uplift rate history. This analysis is recommended for areas of well-preserved coastal terraces, which include at least three terrace levels per profile. The horizontal scaling of modeled LEM topography and mapped TerraceM-2 MAP swath profiles is adjusted using a single shoreline angle of known age as piercing point. After modeled and measured profiles are scaled, the slope, diffusion, initial erosion rate, uplift rate, as well as other LEM parameters, can be adjusted to fit the width and height of the paleo-platforms and cliffs to reproduce the entire sequence of coastal terraces.
To search for the best-fitting modeled to measured topography, we developed two statistical approaches: (a) iteration of LEM parameters based on defined ranges and increments; and (b) a random optimization method based on random sampling Monte Carlo -based scheme. The first method is useful to search for discrete time-independent parameters. The second method can be applied when higher degrees of freedom are involved in the computation, such as for time-dependent scenarios involving variable uplift or erosion rates. Best-fitting models are selected based on Root Mean Square (RMS, Eq. 6) and Normalized Root Mean Square Error (NRMSE, Eq. 7) statistics, where Zobs is the maximum distribution of elevations of swath topography and ZLEM is the LEM elevations. TerraceM-2 LEM also includes the option to find the best fit by minimizing the area between model and observations. To perform the analyses, we developed the "LEM iteration and fit interface, " which allows an interactive selection of the modeling parameters and quick analysis of profiles ( Figure 5). The best-fit models can be exported as graphic (.pdf) and Matlab R data (.mat) files. The LEM model and LEM fit functions can be called independently from command line to develop custom workflows (see SP2.2 and Supplementary Figure 6SP).

TerraceM-2 Dislocation Modeling
Active faults in coastal areas may produce local relative sealevel changes, surface ruptures, and increasing the magnitude and arrival times of tsunami waves (e.g., Wendt et al., 2009;Melnick et al., 2012;Jara-Muñoz et al., 2017). Therefore, identifying seismically active structures and estimating their slip rate are crucial to assess the spatial and temporal characteristics of seismic hazards in coastal areas. TerraceM-2 includes a novel routine to model fault dislocations and surface deformation based on the spatial distribution of coastal terrace elevation or uplift rates (Figure 7). TerraceM DIM is based on the fault-dislocation theory of Okada (1985) that provides an analytical solution for surface deformation resulting from slip along a fault in a homogeneous elastic half-space. Fault parameters such as fault-slip rate and length are crucial to infer potential earthquake scenarios as well as recurrence times and earthquake magnitude (e.g., Anderson and Menking, 1994;Wells and Coppersmith, 1994;Jara-Muñoz et al., 2017).
Frontiers in Earth Science | www.frontiersin.org 7 October 2019 | Volume 7 | Article 255 TerraceM-2 DIM is designed as a workflow including three steps: (a) loading inputs, (b) sensitivity tests to define the ranges of fault parameters, and (c) iterating fault parameters to search for a best-fitting model. Fault parameters include strike, length, slip, rake, down-dip, up-dip, and dip angle; fault geometries can be imported from GIS platforms as shapefile (.shp). Furthermore, the TerraceM-2 DIM functions can be executed from command line to create custom workflows (see SP2.3 and Supplementary Figure 7SP). The TerraceM-2 DIM GUI is designed to perform quick sensitivity tests to determine the approximate solution and set the parameter ranges (Figure 7). Fault parameters are then iterated using random optimization for an unlimited number of faults. The Weighted Root Mean Square Error (WRMSE) is calculated by comparing modeled uplift or elevations with measured shoreline angles, and minimized to determine the best-fitting model using the WRMSE (Eq. 8), where Zobs are the elevations of shoreline angles, ZobsE are the associated errors, and Zmodel are the surface displacements of the dislocation model at the location of the shoreline angles. The best-fit model can be exported as geotiff raster and added as an

EXAMPLES
Three examples are provided to illustrate the capabilities of TerraceM-2 modules and also the use of topographic data with different spatial resolution. The spatial resolution ranges from very-high (0.5 m/px) for drone surveys in the Dead Sea, high (2 m/px) for airborne Light Detection and Ranging (LiDAR) surveys in California, and moderate (12-m) for radar topography such as TanDEM-X along the coast of Peru. These need for moderate, high, or very-high resolution data will depend on the spatial scale of the geomorphic features of interest and the analysis method, which will be illustrated in the following case studies.

Mapping Lacustrine Shorelines in the Dead Sea
The Dead Sea is the deepest continental basin on the Earth, located along the left-lateral transform boundary Dead Sea transform fault (DSF, Figure 8A), between the Arabian and Sinai plates (Garfunkel, 1981;Garfunkel and Ben-Avraham, 1996), characterized by ongoing tectonic activity. During the Late Pleistocene the Dead Sea basin experienced a dramatic lakelevel drop of ∼200 m recorded in the sedimentary sequences of the Lisan Formation and as sequences of regressive lacustrine shorelines that straddle the coast of the Dead Sea (Figures 8B,C) (Bowman, 1971;Bartov et al., 2002;Ghazleh and Kempe, 2009;Torfstein et al., 2013). The outstanding exposure and excellent preservation of dozens of fossil lacustrine shorelines represent a perfect site to test the mapping capabilities of TerraceM-2 and to track permanent deformation accumulated by the DSF. Lacustrine shorelines at Ein Gedi site, Israel, comprise a sequence of 21 levels between -185 and -355 m asl Frontiers in Earth Science | www.frontiersin.org 9 October 2019 | Volume 7 | Article 255  (Bowman, 1971(Bowman, , 1997. The shorelines morphology is characterized by paleo-platforms of 5-40 m width partly covered by clays and silt; in contrast, the shoreline cliffs are characterized by accumulations of boulders and slight changes in slope ( Figure 8C). The site in Jordan is located in front of the Lisan Peninsula, near the Potash company area, and comprises a sequence of 23 shorelines formed between -170 and -355 m asl and almost continuously exposed along the coast in this area (Ghazleh and Kempe, 2009). The morphology of the shoreline cliffs and platforms is similar to Ein Gedi. The overall morphology of both sites resembles a bar-code formed by light and dark fringes corresponding to the paleo -platforms and -cliffs ( Figure 8B). The ages of the lacustrine shorelines at both sides of the Dead Sea range between 20 and 30 ka (Ghazleh and Kempe, 2009;Jara-Munþoz et al., 2019), and record the progressive drop of the lake during the Late Pleistocene. We used the TerraceM-2 MAP module to estimate the amount of permanent deformation accumulated by the DSF, comparing the distribution of lacustrine shorelines at both flanks of the Dead Sea. Lacustrine shorelines along the Dead Sea are characterized by a smooth topographic expression that complicates the mapping of shoreline angles using topographic swath profiles. Therefore, instead of using swath profiles we studied the metrics of lacustrine shorelines using the SCM that can be adjusted to enable the distinction between slight changes of slope, differentiating shoreline paleo-platforms of different levels (Figures 9A,B). We performed this analysis using highresolution topography, and complementing our mapping with field geomorphic observations. We surveyed high-resolution topography in Israel and Jordan using a photogrammetric drone, the images were processed using structure from motion software and vertically georeferenced with ground control points, obtaining a 0.5 m/px digital terrain models (Figures 9A,B).
We adjusted the roughness and slope thresholds of the SCM between 0.1 and 0.3 and between 8 • and 12 • respectively. We converted the shoreline elevations to probability density functions, using the mean value of the distribution of the elevations of each shoreline level to perform the comparisons between sites (Figure 9C). Due to the lack of chronological constraints of individual shoreline levels we attempt to correlate shorelines between Potash and Ein Gedi comparing first levels of lower elevations considering that accumulated less deformation and thus are preserved at similar elevation and by assuming that the sequence of shorelines at both sites is continuous. Our results display a linearly incremental pattern of vertical offsets between shorelines at both sites (Figure 9D), ranging from 0 to 4 m for the lower shorelines (levels 13 to 22) and between 5 and 14 m for the higher shorelines (levels 1 to 12). The shoreline offsets represent relative vertical displacements between both sites; therefore, we cannot identify which of the sites has been uplifting or subsiding.
However, it is possible to distinguish that lacustrine shorelines at the Potash site are higher than at Ein Gedi ( Figure 9C). If we consider that the maximum age of these shorelines is ∼30 ka (Ghazleh and Kempe, 2009), then the mean differential vertical velocity between both sites would be ∼0.5 m/ka. Furthermore, both sites are located on the footwalls of the main fault strands Frontiers in Earth Science | www.frontiersin.org 11 October 2019 | Volume 7 | Article 255 of the DSF (Figure 8A), located offshore along the Israeli and the Jordan coasts and at similar distance from the sites (∼5 km). When compared to the DSF strand in Israel, our results suggest a faster slip rate of the DSF along the fault strand near the coast of Jordan. Despite the limited age constraints and the possible uncertainties in our correlation method, this example highlights the applicability of TerraceM-2 morphometric analyses in mapping lacustrine shorelines and in efforts to detect and measure subtle manifestations of surface deformation. These new approach open new possibilities for potential application in other terminal lakes exhibiting emerged lacustrine shorelines, such as at Uyuni salt lake in Bolivia and lake Bonneville in US among others (Chen and Maloof, 2017;Zeilinger et al., 2018).

Modeling Fault Slip Rates Using Marine Terraces at Santa Cruz, California
We compare fault slip rates derived from marine terrace deformation patterns obtained by previous studies and by TerraceM-2 DIM in the Monterrey bay, California. This area is characterized by continuous exposure of several levels of extensive marine terraces that can be followed over distances of hundreds of kilometers along the coast (Alexander, 1953;Bradley and Griggs, 1976). The marine terraces are preserved at elevations between ∼30 and ∼200 m asl and displaying a notorious warping apparently related to deformation of the plate boundary strike slip San Andreas fault (SAF, Figure 10A). Anderson and Menking (1994) associated marine terrace deformation patterns to coseismic fault slip during earthquakes similar to the Loma Prieta earthquake (Mw 6.9) occurred in 1989, which epicenter was located 15 km east of the bay (Lisowski et al., 1990;Wald et al., 1991). Using elastic dislocation models, the deformation patterns of the whole sequence of marine terraces, and the rupture parameters of the Loma Prieta earthquake, Anderson and Menking (1994) determined a long term vertical fault slip rate between 2.4 and 4.1 m/ka.
In this example we focus on the Highway terrace, the lowest and best-preserved terrace level in the area bounded by a prominent paleo-cliff. We use TerraceM-2 MAP and LiDAR data (2-m/px) available from the National Oceanic and Atmospheric Administration (NOAA 1 ) to determine shoreline angles from 46 swath profiles ( Figure 10B). We obtained a continuous record of shoreline angles along 25 km of coast, ranging between ∼30 and ∼60 m asl and increasing in elevation eastwards. Anderson and Menking (1994) tentatively correlated the Highway terrace with MIS 5; however, due to the lack of chronological constraint at that time the authors were unable to distinguish between the MIS5e (∼125 ka) and MIS 5c (∼105 ka) levels. Therefore, for this example we compared fault slip rates considering both possible ages. We use TerraceM-2 DIM to generate a dislocation model based on the fault parameters of the Loma Prieta earthquake; we consider a single reverse fault oriented N50 • W and with a 70 • SW dip (Wald et al., 1991;Anderson and Menking, 1994), 50 km length (Smith and Sandwell, 2006), and up-and downdip depths of 1 and 17 km, respectively (Lisowski et al., 1990;Anderson 1 www.noaa.gov and Menking, 1994). To simplify our model, we neglect lateral displacements using a rake of 90 • . We iterate 1000 models ( Figure 10D) using random slip values within the range of previous estimates, determining the best fit for a vertical dip and fault slip of 255.2 m (Figure 10C), equivalent to a fault slip rate of between 2.1 and 2.6 m/ka for MIS 5e and MIS 5c ages, respectively, which is similar to the previous estimates of Anderson and Menking (1994). This example and highlights the advantage of integrated mapping and modeling routines of TerraceM-2 and also the quick and easy application of modeling methods, which would usually require advanced knowledge in computer programming.

Modeling the Evolution of Marine Terraces in Southern Peru
We compared the uplift rates estimated with TerraceM-2 LEM and previous estimates at Cerro el Huevo, southern Peru. This area is characterized by remarkably well-preserved and laterally continuous marine terrace sequences that form a staircase morphology. The Cerro el Huevo site is located adjacent to San Juan de Marcona bay, above the south-eastern flank of the Nazca Ridge, which is subducted beneath the South American plate (Figure 11A). Previous studies in this area (Hsu, 1988(Hsu, , 1992Ortlieb and Macharé, 1990;Saillard et al., 2011) have described staircase sequences of uplifted wave-cut marine terraces comprising more than 10 levels reaching ∼250 m asl and etched on the western slope of Cerro el Huevo ( Figure 11B). The chronology of these terraces has been partly determined by amino acid racemization spanning the period between the marine isotope stages (MIS) 3 and 11 (∼60 to ∼400 ka) (Hsu, 1988(Hsu, , 1992; additionally, Ortlieb and Macharé (1990) determined U/Th ages for the period between MIS 5 and 7 (∼125 to ∼230 ka). However, most of these ages were reinterpreted in the context of cosmogenic radionuclide dating that provided ages between MIS 7 and 11 (∼230 to ∼400 ka) (Saillard et al., 2011). Saillard et al. (2011) compared marine terraces ages to sea-level high-stands using the compilation of sea-level variations of Siddall et al. (2007), obtaining uplift rate estimates between 0.4 and 0.9 m/ka. However, according to Saillard et al. (2011), these rates vary along time from 0.6 m/ka between 400 and 100 ka to 0.8 m/ka after 100 ka, associated with the progressive southward migration of the Nazca ridge. Considering the complex uplift history of this area, and the availability of marine terrace ages and TanDEM-X topography data (12 m/px), Cerro el Huevo represents a perfect site to test the modeling and mapping capabilities of TerraceM-2.
To reproduce the Cerro el Huevo marine terrace sequence, we apply the LEM with the sea-level curve of Rohling et al. (2009), which is based in part in Siddall et al. (2007), but provides a better resolution for high-stand periods, allowing the generation of detailed sequences of synthetic LEM terraces. The initial slope of the LEM was fixed to 11 • based on the angle between the highest elevation point and the coastline along the swath profile. The optimal E0 was determined after sensitivity tests, obtaining an initial value of 0.05 m/yr that better reproduce the average width of the marine terrace platforms. To simplify our model, we neglected cliff diffusion. We measure the elevation of marine terraces using 12 m/px TanDEM-X topography, provided by the Deutschen Zentrums für Luft-und Raumfahrt (DLR), extracting a swath profile of 2.5 × 0.3 km that included the whole sequence of marine terraces ( Figure 11B). The modeled and measured marine terraces topography were tied using the shoreline angle of the MIS 5e terrace level ( Figure 12A) located at ∼100 m asl (Saillard et al., 2011). To find the best-fitting model that reproduces both the topography and ages, we used an optimization technique, generating random patterns of variable uplift rates within a range that included the previous uplift rate estimates of Saillard et al. (2011). We iterated 150 uplift rate histories and determined the best fit model by minimizing the NRMSE; this model comprises a history of uplift variability, including rates between 0.1 and 0.2 m/ka during the period between 500 and 350 ka, 0.3 and 0.5 m/ka between 350 and 150 ka, and between 0.9 and 1 m/ka since 150 ka ( Figure 12B). This difference between our LEM-based uplift rate history and the previous estimates of Saillard et al. (2011) is interesting as it points to a common flaw in the methodology to determine uplift rates using shoreline angles. The LEM minimization method identifies the modeled topography that optimally fits the maximum swath topography; therefore, it may be argued that the discrepancy between our results and those of Saillard et al. (2011) could be related to an incorrect fit between terrace levels of different ages. However, as may be observed in Figure 12A, the modeled ages and the terrace ages from Saillard et al. (2011) are in good agreement. In fact, Saillard et al. (2011) estimated uplift rates by comparing the presentday shoreline angle elevation with the sea-level position and the age of the associated high-stand. They missed out on taking into account the fact that surface uplift is incremental, and therefore temporal variations in uplift rates will affect cumulative vertical displacements. Thus, uplift rate determinations that exclusively consider the present-day elevation of marine terraces may include systematic biases. For instance, if uplift rate increase over time, older marine terraces will accumulate the sum of vertical displacements, and thus uplift rates estimates based on the present-day elevation of shoreline angles will be overestimated. We suggest that this was the case for the measurements of Saillard et al. (2011), which overestimated uplift rates of the older marine terraces (MIS 7 to 11, Figure 12B). Our uplift rate history for the Cerro El Huevo in southern Peru suggests a continuous increase of uplift rates over the past ∼500 ka, which agrees with reconstructions of the southward-migrating collision between the Nazca ridge and the South American margin (Hampel, 2002). This example highlights the advantage of LEMs to resolve vertical displacements under complex uplift rate histories, which would be difficult using conventional marine terrace analysis methods.

SUMMARY
Coastal terraces have been widely used to estimate past water level positions and vertical deformation rates in marine and lacustrine landscapes. Since its release in 2016, TerraceM has been used in different coastal settings around the world to map terraces in marine and lacustrine realms quantitatively. The upgraded version of TerraceM-2 improves processing efficiency and includes new modules and ancillary functions to map, process, analyze, model, export, and display morphometric and spatial characteristics of marine and lacustrine terraces. Next, we describe the advantages and limitations of TerraceM-2.
Advantages of TerraceM-2: • Interface: TerraceM-2 includes GUIs that facilitate the analysis of terraces for inexperienced users, and also includes independent functions designed for custom command line programming for advanced users.
• User experience: The intuitive TerraceM-2 interface design includes interactive message windows to guide the user through the different modules. • Integration: TerraceM-2 modules are designed to exchange data between them allowing to easily integrate the mapping and modeling capabilities. • Efficiency: TerraceM-2 improved its processing speed at least eleven times in comparison with the previous release of TerraceM.
Limitations of TerraceM-2: • Objective: TerraceM-2 is designed for mapping only marine and lacustrine terraces. However, the SCM interface may be used to map any planar landform, such as peneplains or fluvial terraces associated with both aggradation or erosion. • Inputs: TerraceM-2 is designed to work with highresolution topography. However, the resolution will • depend on the spatial scale of the targeted geomorphic features. We present three examples that rely on data at different resolutions to illustrate the use of TerraceM-2 to analyze geomorphic features with different spatial scales.
• Software platform: TerraceM-2 is open source but runs under Matlab R 2017 or later, which is a paid software. However, MathWorks R provides free Matlab R licenses to students and low-cost licenses for educational and research institutions. We favor Matlab R above other open-source programming software because of its versatility with GIS data. • Efficiency: TerraceM-2 can be used in a normal desktop or notebook computer running Windows R , Linux R , or Mac R operational systems. However, its efficiency depends of the technical characteristics of each computer (e.g., amount of memory or processor speed).
We present three examples to highlight the new features in TerraceM-2: lacustrine shoreline mapping using very-highresolution drone topography (0.5 m/px) in the Dead Sea, estimating fault slip rates using the deformation pattern of marine terraces in California using high-resolution LiDAR data (2 m/px), and estimating temporal variability of coastal uplift rates associated with ridge subduction in southern Peru using TanDEM-X data (12 m/px). These examples provide a guide on the resolution of the topographic data needed to map, model, and analyze different wavecut geomorphic features quantitatively and highlight the novel methodological improvements of TerraceM-2, which will help the research community studying marine and lacustrine terraces.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the manuscript/Supplementary Files.

AUTHOR CONTRIBUTIONS
JJ-M developed TerraceM and wrote the manuscript. DM participated in the conceptual and technical development of the TerraceM-2 functions and prepared the manuscript. KP compiled a worldwide database of marine terraces included in this version of TerraceM-2 and actively participated in the development of the manuscript and in the discussion of marine terrace concepts. MS actively assisted in the development of this manuscript, the discussion of marine terrace concepts, and quantitative mapping methods.

FUNDING
This study was supported by the Deutsche Forschungsgemeinschaft (DFG), the Project "LIFE: Linking Lake-level variations, surface deformation and seismogenesis in the Dead Sea" Grant JA 2860/1-1 funded by DFG, Grant STR 373/30-1, project ME 3157/4-1 funded by DFG, and the Millennium Science Initiative (ICM) of the Chilean government through Grant NC160025 "Millennium Nucleus the Seismic Cycle along Subduction Zones, Cyclo, " National Fund for Scientific and Technological Development (FONDECYT) grants 1150321, 1181479, and 1190258. MS was supported by the German Israeli Foundation (GIF) (Grant No. I-1280-301.8/2014.