Geomorphic analysis of transient landscapes in the Sierra Madre de Chiapas and Maya Mountains (northern Central America): implications for the North American – Caribbean – Cocos plate boundary

We use a geomorphic approach in order to unravel the recent evolution of the diffuse triple junction between the North American, Caribbean, and Cocos plates in northern Central America. We intend to characterize and understand the complex tectonic setting that produced an intricate pattern of landscapes using tectonic geomorphology, and available geological and geophysical data. We classify regions with specific relief characteristics and highlight uplifted relict landscapes in 5 northern Central America. We also analyze the drainage network from the Sierra Madre de Chiapas and Maya Mountains in order to extract information about potential vertical displacements. Our results suggest that most of the landscapes of the Sierra Madre de Chiapas and Maya Mountains are in transient stage. Topographic profiles and morphometric maps highlight elevated relict surfaces that are characterized by a low amplitude relief. The river longitudinal profiles display 10 upper reaches witnessing these relict landscapes while lower segments characterized by multiple knickpoints, that adjust to new base-level conditions. These results backed by published GPS and seismotectonic data allow us to refine and extend existing geodynamic models of the triple junction. Relict landscapes are delimited by faults and thus result from a tectonic control. The topography of the Sierra Madre de Chiapas evolved as the result 15 of (1) the inland migration of deformation related to the coupling between the Chiapas Massif and the Cocos fore-arc sliver, and (2) the compression along the northern tip of the Central America Volcanic Arc. Although most of the shortening between the Cocos fore-arc sliver and the North American plate is accommodated within the Sierra de Chiapas and Sierra de los Cuchumatanes, a small part may be still transmitted to the Maya Mountains and the Belize margin through a “rigid” 20


Introduction
The aim of this work is to examine geomorphic features along two key areas (the Sierra Madre de Chiapas and the Maya Mountains, Fig. 1) in order to test existing models of the North American-Caribbean-Cocos plate boundary (e.g., Malfait and Dinkelman, 1972;Burkart, 1983;Guzmán-Speziale 25 et al., 1989;Lyon-Caen et al., 2006;Ratschbacher et al., 2009;Authemayou et al., 2011). The oceanic Cocos plate is subducted beneath the North American and Caribbean plates along the Middle America Trench while the North American-Caribbean plate boundary is a sinistral transform system which accommodates the eastward escape of the Caribbean plate(e.g., Lyon-Caen et al., 2006;Manea and Manea, 2006;Andreani et al., 2008a;Authemayou et al., 2011). 30 The complexity of the plate boundary comes from the fact that both the Caribbean and North American plates are limited to the west by a forearc sliver (Fig. 1) which is coupled to the Cocos slab (e.g., Turner et al., 2007;Phipps-Morgan et al., 2008). This sliver rotates counterclockwise as it is pulled by the eastward escape of the Caribbean plate (e.g., Andreani et al., 2008a;Authemayou et al., 2011). Over the past decades, numerous models regarding the structure and the evolution of the 35 triple junction have been proposed (Authemayou et al., 2011, and references therein). Recent models agree on the fact that the dextral Jalpatagua fault (Fig. 2) coupled with extension along the grabens of Guatemala represents the limit between the forearc sliver and the Caribbean plate (e.g., Lyon-Caen et al., 2006;Andreani et al., 2008a;Authemayou et al., 2011;Franco et al., 2012). However, there is no clear consensus on the boundary between the forearc sliver and the North American plate. This is 40 due to the fact that this boundary is highly diffuse and most of the active deformation is distributed over the Sierra Madre de Chiapas orogenic belt. On top of that, none of the revisited tectonic models of the North American-Caribbean-Cocos plate boundary attempted to include the Maya Mountains ( Fig. 1 and 2). In spite of evidence for recent (Late Neogene or Pliocene) tectonics (e.g., Weidie, 1985;Lara, 1993;Purdy et al., 2003;Bauer-Gottwein et al., 2011), this region was seldom studied Finally, we combine the results from the geomorphic analyses with available geophysical data (GPS and seismicity, Fig. 3) in order to propose a model for the North American-Caribbean-Cocos triple 60 junction.
2 Tectonic setting 2.1 Models of the triple junction Early geodynamic models of the junction between the North American, Caribbean and Cocos plates consisted in a simple transform-trench boundary. For instance, Burkart (1983) proposed a model in which the Polochic fault (and perhaps also the Motagua fault) connected to the Middle America trench. He proposed that the 300 km Neogene motion between the North American and Caribbean plates was accommodated by left-lateral slip along the Polochic and Motagua faults and by extension within the depressions of Honduras and Guatemala. The slip between these two plates was compensated by a left-lateral offset of the trench. 70 Later works challenged the idea of a junction between the Polochic-Motagua fault system and the trench because the Polochic fault trace terminates within the Chiapas Massif and the Motagua fault trace is lost within the Central America volcanic arc (e.g., Guzmán-Speziale et al., 1989;Guzmán-Speziale and Meneses-Rocha, 2000). As a result, more recent models (Fig. 4) argue that the motion between the North American and Caribbean plates does not result in an offset of the trench but is 75 rather absorbed onshore within a complex zone of deformation. However, there is no clear consensus on the geometry of the plate boundary and on how the deformation is distributed inland.
Using seismotectonic data (Fig. 3a), Guzmán-Speziale and Meneses-Rocha (2000) proposed a model in which the Polochic-Motagua fault system is unable to propagate across the Chiapas Massif ( Fig. 4a). As a result, a part of the motion between the North American and Caribbean plates is 80 absorbed by the reverse and strike-slip faults of Chiapas (Fig. 4a). In this model, the reverse faults act as a restraining bend between the strike-slip faults and the Polochic-Motagua fault system. Later on, Guzmán-Speziale (2001) suggested that part of the plate motion is also absorbed by the grabens of Central America, south of the Motagua fault (Fig. 4a).
A second category of models (Fig. 4b) mainly focused on the extensional province located south Motagua fault system and by extension within the grabens of Central America (Fig. 4b). This model was refined by Franco et al. (2012) which proposed a difference in coupling along the subduction 95 interface (Fig. 4b) in order to explain the GPS velocity field in southern Mexico (Fig. 3b). However, these two models did not addressed the tectonic evolution of the Sierra Madre de Chiapas.
The latest models ( Fig. 4c and 4d) attempted to explain both the eastward escape of the Caribbean plate as described by Lyon-Caen et al. (2006) and tectonics of the Sierra Madre de Chiapas. Andreani et al. (2008b) suggested a connection between the strike-slip faults of Chiapas, the Veracruz shear 100 zone which affects the Veracruz basin further north, and the transtension of the Mexican volcanic arc (Fig. 4c). They proposed that these three zones represent the boundary of the so-called 'Southern Mexico block ' Andreani et al. (2008a). In their model (Fig. 4c), the dynamic of the plate boundary is driven by the escape of the Caribbean plate and by the counterclockwise rotation of the Southern Mexico block. However, this model does not solve the interactions between the Southern Mexico 105 block and the Central America forearc sliver. An alternative model was proposed by Authemayou et al. (2011), in which the Central America forearc sliver extends offshore southern Mexico (Fig. 4d).
In their model, the eastward escape of the Caribbean plate induces a counterclockwise rotation of the forearc sliver and a 'zipper' process (i.e., progressive suturing) between the Jalpatagua and Motagua fault as the space between the two faults is left 'empty' by the moving Caribbean plate. According to 110 their model, the Tonalá Shear Zone is a suture resulting from this 'zipper' process. However, both the Motagua and Jalpatagua fault traces are lost west of 91 • W. The western termination of both faults seems to be associated with extension along the Guatemala City graben (Lyon-Caen et al., 2006). Furthermore, the Tonalá Shear Zone clearly connects to the Polochic Fault (Fig. 2, Guzmán-Speziale et al., 1989;Guzmán-Speziale and Meneses-Rocha, 2000). 115 In summary, current views on the plate boundary agree on the fact that the North American and Caribbean plates are limited by a forearc sliver to the west. The counterclockwise rotation of the forearc sliver results in transpression north of the Polochic-Motagua fault system while the eastward escape of the Caribbean plate is accommodated by the strike slip motions of the Motagua and Jalpatagua faults and by extension along the grabens of Ipala and Guatemala city. However, the 120 connections between the western corner of the Caribbean plate and the transpressional tectonics of the Sierra Madre de Chiapas are still unclear. In addition, none of these recent models attempted to integrate the recent tectonics of the Maya Mountains and Yucatán peninsula into the dynamics of the plate boundary.
A part of the motion between the North American and Caribbean plates is also accommodated along at least thirteen N-trending grabens (e.g., Mann and Burke, 1984;Guzmán-Speziale, 2001) located south of the Motagua fault. The main structures form the Guatemala City, Ipala, Esquipulas, 140 Sula and Comayagua depressions (Fig. 2). Rogers et al. (2002) related these grabens to the uplift of the Central America plateau following a slab break-off underneath Central America. The onset of extension has been determined at 11-8 Ma (Gordon and Muehlberger, 1994;Ratschbacher et al., 2009). However, most of the present-day extension (11-12 mm.yr −1 ) seems to be accommodated along the westernmost grabens (Lyon-Caen et al., 2006;Rodriguez et al., 2009;Franco et al., 2012).

Sierra Madre de Chiapas
The Sierra Madre de Chiapas is constituted by four main structural domains: the Chiapas Massif, the Central Valley, the Sierra de Chiapas and the frontal fold-and-thrust belt (Guzmán-Speziale and Meneses-Rocha, 2000;Andreani et al., 2008a;Witt et al., 2012b). The Chiapas Massif mostly consists of Permian igneous rocks affected by medium to high grade metamorphism and Middle to late 150 Miocene intrusions related to the subduction of the Cocos plate (Damon and Montesinos, 1978;Weber et al., 2007;Molina-Garza et al., 2015). The western flank of the Chiapas Massif is bounded by the left-lateral Tonalá Shear Zone which connects to the south to the western termination of the Polochic fault (Wawrzyniec et al., 2005;Weber et al., 2005). The displacement along the Tonalá Shear Zone is synchronous with Miocene magmatic intrusions and could have reached 100 km 155 (Molina-Garza et al., 2015). The Central Valley is a ∼170 km long and ∼30 km wide depression corresponding to a NW-trending synclinorium with superimposed smaller folds (Witt et al., 2012b).
The Sierra de Chiapas has a roughly sigmoidal shape and is bounded to the west by the prominent NW-trending left-lateral Tuxtla fault. It is constituted by blocks culminating at heights between 2000 m and 2400 m and delimited by E-trending left-lateral faults defined as the High Sierra Fault 160 System by Witt et al. (2012b). These faults probably connect to the Malpaso fault to the west while strike-slip motion is absorbed at their eastern terminations by thrusting and folding. The Chiapas frontal fold-and-thrust belt is located between the Sierra de Chiapas and the western border of the Yucatán platform. It is constituted by closely spaced faulted folds which are rooted in shallow levels (Witt et al., 2012b). 165 Thermochronological, tectonic and stratigraphic evidence suggest that renewed exhumation and topographic growth occurred along the Chiapas region during the middle Miocene (16-10 Ma) and late Miocene-Pliocene (6-5 Ma), following a phase of rapid exhumation to upper crustal levels at ∼30 Ma Witt et al., 2012a). The Tonalá shear zone may have accommodated significant deformation since 16 Ma while the displacement along the transpressive Tuxtla 170 and Malpaso faults occurred during the last 6-5 Ma and could have reached 50 to 70 km, involving 0.5 to 0.8 cm.yr −1 of left-lateral motion Witt et al., 2012b).

Maya Mountains and Belize margin
The Maya mountains and the Belize margin are located north to the Polochic-Motagua fault system and east to the Yucatán platform. Due to a dense vegetation cover and scarce roads, this region 175 was seldom studied in comparison with other areas related to the North American-Caribbean plate boundary.
The eastern border of the Yucatán peninsula is crosscut by a series of NNE-trending normal faults (e.g., Weidie, 1982Weidie, , 1985, which are referred as East Yucatán fault zone (EYFZ) in Fig. 2. The fault zone is ∼80 km wide and extends over 500 km between the NE tip of the Yucatán peninsula and the 180 Maya Mountains. The surface expression of this fault zone is seen in the alignment of topographic scarps, hydrological features (cenotes, lakes and rivers) and coastal bays (Weidie, 1985;Lesser and Weidie, 1988;Bauer-Gottwein et al., 2011). The East Yucatán fault zone represents the onshore continuation of an extensive horst and graben system which affects the western margin of the Yucatán Basin (Weidie, 1985;Rosencrantz, 1990).

185
According to Rao and Ramanathan (1988) and Purdy et al. (2003), the Maya Mountains correspond to a roughly NE-trending structural high where the Paleozoic basement is uplifted. Paleozoic rocks constitute the highest elevations of the Maya Mountains. The orogen mainly consists in metamorphosed Late Carboniferous to Middle Permian volcano-sedimentary rocks overlying Late Silurian granites and are bounded by the Northern and Southern Boundary faults (Kesler et al., 1974;190 Bateson and Hall, 1977;Steiner and Walker, 1996). The Maya Mountains are delimited by faultbounded E-to NE-trending depressions: the Corozal Basin located to the north and the Belize Basin that borders the Maya block offshore to the east and onshore to the south. According to Purdy et al. (2003), unroofing of the Cretaceous carbonate cap of the Maya Mountains siliciclastic sediment source did not occur until late Neogene, perhaps no earlier than late Pliocene. Faulting along the 195 Belize margin is not accurately documented. The most recent tectonic event is Pliocene or younger and the resulting structures affected Quaternary carbonate deposition. However, Lara (1993) related this event to transtensional faulting while Purdy et al. (2003) rather interpreted this event as the result of a transpression.

Swath topographic profiles
Swath topographic profiles condense elevation data of a complex region into a single profile (e.g., Isacks, 1992;Masek et al., 1994;Telbisz et al., 2013;Hergarten et al., 2014). Topography is extracted from a rectangular swath using a series of parallel profiles, rather than using a single line as in conventional topographic profiles. Elevation data are then projected onto a vertical plane parallel 205 to the longitudinal axis of the swath rectangle, and statistical parameters (the maximum, minimum and mean elevations) are calculated. The curve for maximum elevation corresponds to the ridgelines and helps to identify topographic features, such as relicts of paleo-surfaces. The curve for minimum elevation corresponds to the valley floors or river beds. A quick estimate of the incision is given by the arithmetic difference between the maximum and minimum elevations in a given longitudinal 210 distance (window).
Swath topographic profiles were extracted from 3 arc-seconds SRTM data from CIAT (Jarvis et al., 2008) using a MATLAB script. The swath width was fixed to 20 km. This width is large enough to contain both elevated surfaces and major rivers and small enough to avoid topographic features that are too oblique with respect to the swath axis. Elevation data were sampled using ∼220 215 parallel profiles separated by 90 m. Elevations along each individual profile were also sampled using a 1-pixel (90 m) interval.

Surface analyses
The combined use different morphometric indices proves to be an efficient way to classify landscapes according to their state of dynamic equilibrium (e.g., Andreani et al., 2014;Domínguez-González 220 et al., 2015). Hypsometric integral (HI) efficiently highlights topographic scarps, surface roughness (SR) substantially increases with incision and relief anomaly (RA) and surface index (SI) highlight elevated low relief landscapes.
The hypsometric integral shows the distribution of landmass volume remaining beneath or above a basal reference plane (Strahler, 1952;Schumm, 1956). This index proved to be efficient in eval-225 uating the response of landscapes to active tectonics (e.g., Pérez-Peña et al., 2009;Mahmood and Gloaguen, 2012;Siddiqui and Soldati, 2014;Andreani et al., 2014;Domínguez-González et al., 2015). According to Pike and Wilson (1971), the hypsometric integral (HI) for a given area can be approximated with Eq. 1: with h mean , h min and h max being the mean, minimum and maximum elevations of the analyzed area.
Surface roughness can be described by several parameters (Smith, 2014, and references therein).
In this work we used the area ratio approach which evaluate the similarities between a topographic surface with a given area and a flat surface with the same geographic extend (e.g., Hobson, 1972;235 Grohmann, 2004;Grohmann et al., 2009;Shahzad and Gloaguen, 2011b). The ratio is close to 1 for flat areas and increases rapidly as the real surface becomes irregular. The method used to obtain the 'real' and flat surfaces is adapted from the GRASS-R algorithm of Grohmann (2004). First, a slope map is produced using the neighborhood algorithm included in ArcGIS (Burrough and Mcdonell, 1998). The 'real' surface S R is then obtained for each pixel using Eq. 2: where res is the DEM resolution in meters and α is the pixel slope in degrees. The flat area S F is defined for each pixel by S F = res × res. The surface roughness is then obtained by summing the pixel values of S R and S F within a moving window and by dividing the sum of S R pixels by the sum of S F pixels.

245
We computed the hypsometric integral and surface roughness from 90 m resolution SRTM data (CIAT, Jarvis et al., 2008) using TecDEM, a MATLAB-based toolbox (Shahzad and Gloaguen, 2011b). Each pixel of the output raster represents the hypsometric integral and surface roughness values for a 15 × 15 km moving window. The choice of the window size depends on the objectives.
The size of the window will affect the amount of smoothing on the relief (low-pass filtering). To 250 get rid of high-frequency signal such as erosion patterns, the window must encompass a sufficient portion of the analyzed landscape and thus be larger than the wavelength of the undesired signal.
Otherwise, computed values will reflect local-scale variations in relief such as hillslopes. Smallest window scales, aside from being noisy due to high-frequency signals, can also be more influenced by errors (e.g, Albani et al., 2004;Tarolli et al., 2012;Sofia et al., 2013Sofia et al., , 2014. In the present case, 255 a 15 km window allows to filter out the effect of the drainage system as the window size represents an area that will include several valleys and ridges. (Andreani et al., 2014).
Finally, we used the relief anomaly (Scotti et al., 2014) and the surface index (Andreani et al., 2014) to highlight elevated and low relief landscapes. The relief anomaly (RA) represents the elevations normalized by the local relief. We computed this index for a given area using Eq. 3: with h mean , h min and h max being the mean, minimum and maximum elevations of the analyzed area. Highest values are obtained for flat and elevated surfaces.
The surface index (SI) combines the elevations from the DEM with the computed maps of hypsometric integral and surface roughness. It allows to discriminate areas with low local relief landscapes 265 from areas with a more rugged topography. To compute this index, rasters of elevations, hypsometric integral and surface roughness are normalized by their respective minimum and maximum values in order to obtain pixels values between 0 and 1. We then combine the newly created rasters using Eq. 4:

Modeling of the Drainage Network
We extracted the drainage network from 30 m resolution SRTM data. The extraction was done using 275 TecDEM (Shahzad and Gloaguen, 2011a) by calculating flow directions and contributing area for each pixels using the D8 algorithm (O'Callaghan and Mark, 1984;Fairfield and Leymarie, 1991;Jones, 2002). Streams were identified using a minimum contributing area of 1 km 2 . This threshold is introduced in order to take into account the transition from colluvial to stream-flow (fluvial) dominated channel, which is usually observed at ∼ 1 km 2 in temperate humid regions (e.g., Montgomery

280
and Foufoula-Georgiou, 1993; Wobus et al., 2006). It is assumed that a theoretical flow accumulation that drains an area larger than 1 km 2 can be associated to a real water channel. Extracted streams were then organized hierarchically using Strahler (1957) order.
A DEM-based procedure allows to easily extract and analyze a regional-scale drainage network.
However, there are several uncertainties related to the DEM and methods. Extracted drainage net-285 work and contributing areas are affected by the quality and sampling of the DEM. Free and commonly used data include the 1 arc-second (ca. 30 m) and the 3 arc-second (ca. 90 m) resolution SRTM. For both datasets the absolute vertical error is reported to be less than 20 m. In some areas, the available 1 arc-second SRTM data presented voids (i.e. pixels with no data) mainly located in ridges. We filled these voids with 3 arc-second SRTM data.

290
The flow-routing method we use (D8 algorithm) may introduce bias in flow path orientation as it discompose flow directions into units of 45 • (Fairfield and Leymarie, 1991). This is especially the case in flat areas. Additional known bias are related to methodological aspects. In our study area the commonly encountered problem concerned nested depressions which are related to DEM imperfections. These pits need to be filled to create flow directions. Other artifacts are mainly found 295 along entrenched rivers and canyons, where bad pixels values 'block' the path of the extracted rivers.
In two cases (Sumidero and La Venta canyons) the DEM filling resulted is a significant deviation of the modeled rivers with respect to the actual rivers. We thus manually carved the bad pixels blocking these two canyons in order to obtain a more accurate flow path. The original topography is also affected by human artifacts such as dammed rivers. This is especially the case for the Grijalva river 300 in southern Mexico. Artificial flats related to DEM filling or dammed rivers introduces errors in extracted river paths. These errors are easily detected in river longitudinal profiles and were taken into account when interpreting the knickzones in extracted rivers.

Analysis of River Longitudinal Profiles
Deviations from the typical concave-up shape of stream longitudinal profiles, such as knickpoints or 305 convex segments, indicate a disequilibrium state resulting from tectonic, base-level or lithological perturbations (Kirby and Whipple, 2001;Chen et al., 2003;Troiani and Della Seta, 2008;Pedrera et al., 2009;Font et al., 2010). The normalized steepness index (k sn ) is widely used to investigate tectonically-induced perturbations in river longitudinal profiles and has been used to propose patterns of uplift (Kirby and Whipple, 2001;Wobus et al., 2006;Whittaker et al., 2008). The relationships 310 between slope and catchment area which define the equilibrium state channel gradient are given by Eq. 5 and 6 (Flint, 1974;Kirby and Whipple, 2001;Wobus et al., 2006): where S is the local channel slope, θ is the channel concavity, k s is the steepness index, A is the upstream drainage area, U is the rock uplift rate and K is the dimensional coefficient of erosion. As suggested by Wobus et al. (2006, and the references therein) a normalized steepness index k sn is used, since k s and θ are strongly correlated.
In some cases, the upstream portions located above prominent knickpoints are associated with an 320 upper-relict landscape. Landscapes tend towards an equilibrium in which rivers are graded to sea level or local base-level. Tectonic or climatically-induced base-level falls modify the equilibrium of the drainage. The result is an erosion wave propagating upstream and the areas not yet affected by the erosion wave form an upper-relict landscape (e.g., Clark et al., 2005;Reinhardt et al., 2007;Pérez-Peña et al., 2015). The reconstruction of the original stream profile downstream of the convexity 325 provides an estimate for the amount of base-level change and subsequent incision (Schoenbohm et al., 2004;Clark et al., 2005;Gallen et al., 2013). To reconstruct the original stream profile, we used the power law between slope and distance defined by Eq. 7 (Hack, 1957): S is the local channel slope and D is the distance from the drainage divide. Parameters i and j are 330 obtained by regressing the upper segment of the stream profile in a logarithmic plot of slope against distance.
We analyzed river longitudinal profiles using TecDEM (Shahzad and Gloaguen, 2011a). Normalized steepness indices were computed from Eq. 5 using regressions of slope versus catchment area in logarithmic plots with a given reference concavity θ ref = 0.45 (Kirby and Whipple, 2001;Wobus 335 et al., 2006). Prominent knickpoints or convex anomalies can be observed directly on river longitudinal profiles. However, logarithmic plots of slope against catchment area allow a more detailed analysis, as minor anomalies in the gradient of rivers can be easily detected. For each longitudinal profile, we selected and regressed several segments delimited by changes in the gradient of the river.
We then plotted regressed segments and their assigned k sn values on a map. The reconstruction of 340 stream profiles was done by regressing the segments located upstream of knickzones in logarithmic plots of the slope against the distance (Eq. 7). We used a Matlab-based script implemented in TecDEM by Andreani et al. (2014).
The integer format of DEMs locally produces multiple flats with zero slope, which cannot be handled in logarithmic plots (e.g., Kirby and Whipple, 2001;Wobus et al., 2006). This issue is com-345 monly solved by smoothing extracted river profiles (Wobus et al., 2006). However, this may induce bias in extracted geomorphic indices, especially for segments located close to major knickpoints.
In this study we smoothed extracted river profiles using a 20 pixels (ca. 600 m) moving window.
We used a bootstrapping approach in order to address uncertainties related to both DEM artifacts and methods (Andreani et al., 2014). Linear regressions used to compute both normalized steepness 350 indices and reconstruction of stream profiles are based on subsamples. Typically, we randomly selected from logarithmic plots 75 % of the points corresponding to the analyzed segment in order to obtain 500 subsets. We then performed a regression on each subset. We obtain the mean values as well as a range for both normalized steepness indices and reconstructed profiles.

355
The spatial distribution of stream height (isobase map) is a useful proxy for investigating geologic or tectonic processes (e.g., Dury, 1952;Filosofov, 1960;Golts and Rosenthal, 1993;Grohmann et al., 2007Grohmann et al., , 2011. Drainage networks are commonly organized according to Strahler (1957) stream order. According to Golts and Rosenthal (1993) streams of similar orders are of similar geological age and are related to similar geological events. Hence the interpolation of isobase lines, which 360 connect streams with a similar order, produce a surface resulting from the same erosional events. As suggested by Grohmann et al. (2011), isobase maps can be seen as a smoothed version of the original topographic surface, from which was removed the 'noise' of the 1st Strahler order streams erosion.
Isobase maps were computed using the drainage network extracted from the 30 m resolution SRTM data (see section 3.3). Using SAGA GIS software we selected streams with a Strahler order ≥ 2 and we interpolated the elevations from extracted rivers using a natural neighbor method.

Swath topographic profiles
We analyzed the topography of the Sierra Madre de Chiapas and Maya Mountains along five swath profiles (location in Fig. 5a). For each profile ( Fig. 6) we extracted the minimum, mean and maximum elevation along the swath. We estimated the local incision by subtracting the maximum and minimum elevation. We also plotted the annual averaged precipitation derived from the Tropical   The Central Valley (profiles 2 and 3, Fig. 6) corresponds to a ∼40 km wide depression. Despite a mean elevation of ∼600 m the local incision appears to be less than 100 m. The topographic scarp associated to the Concordia fault is barely marked. In contrast, La Venta fault and the Tuxtla fault delimitate a ∼1500 m topographic high (referred as "La Venta block", profile 1 in Fig. 6) which acts as a barrier for rivers coming from the Chiapas Massif.

395
In profile 1, the Tuxtla and Malpaso faults represent the boundary between the La Venta block and the Chiapas fold-and-thrust-belt. In profiles 2 and 3 The Tuxtla fault represents the limit between the Central Valley and two topographic highs: the Sierra de Chiapas and the Comitán High. The Sierra de Chiapas forms an impressive bulge which dominates both the Central Valley and the Chiapas fold-and-thrust-belt. Its central area is relatively flat and culminates at ∼2500 m. The topography 400 of the Comitán High is flat with very low incision by the drainage network and the mean elevation is ∼1600 m. This flat surface is limited to the east by a topographic high (referred as "Leyva Velázquez block", profile 3 in Fig. 6) which belongs to the Chiapas fold-and-thrust-belt. While in profile 1 the elevation of the fold-and-thrust-belt is constant, in profiles 2 and 3 the topography of the belt is asymmetric and elevations gently decrease toward the flat Tabasco coastal plain.
Peaks in precipitation indicate orographic effects along the Pacific side of the Chiapas Massif (profiles 1 to 3 in Fig. 6) and along the northern part of the Chiapas fold-and-thrust belt (profiles 1and 2 in Fig. 6). However, the shadowing effect of these ranges appear limited as we observe no significant drop in the amount of precipitation between the inner parts of the Sierra Madre de Chiapas and suroounding areas (Pacific coast and Tabasco coastal plain). However, the peak in precipitation is located along the Ixcán fault. In spite of a SW gradient (from 2500 to 1500 mm/yr), most of the Sierra de los Cuchumatanes receives as much as rain as the Petén is less than 100 m) and appears to be tilted towards the NW.

Morphometric maps
The hypsometric integral values ( Fig. 7a) appear to be essentially controlled by tectonic features. and Yucatán platform. In other areas, the distribution of the low relief landscapes appears more fragmented. This is the case for the areas located between the Ixcán and Motagua faults, where our maps highlight the remnants of the Mayan paleosurface described by Brocard et al. (2011).
We propose to compare the distribution of geomorphic indices ( Fig. 7 and 8), precipitation data 490 (Fig. 5a) and geological contours (Fig. 1b) within the Sierra Madre de Chiapas. Our aim is to detect possible effects of climate or lithology on the location of incised areas and low relief surfaces. At first glance, peaks in precipitation (Fig. 5a) seem to coincide with more rugged topography (along the Pacific coast, the northern Sierra de Chiapas and the eastern part of the Sierra de los Cuchumatanes). However, peaks in precipitation are located at the limits of the orogenic belt while the 495 rugged topography extends inward the Chiapas orogen. On top of that, some areas display a rugged topography in spite of relatively low precipitation (e.g., west of the Sierra de los Cuchumatanes) and vice-versa (e.g., eastern part of the Chiapas fold-and-thrust belt). To compare the distribution of precipitation and geomorphic indices, we randomly selected 20000 points within the Sierra Madre de Chiapas. For each point we extracted the corresponding pixel values from surface roughness, surface 500 index, hypsometric integral, local relief and precipitation maps. We then compared the distribution of the values using scatterplots (Fig. 9). The co-linearity of the scatterplots was estimated using the Pearson correlation coefficient (PCC). PCC values vary between +1 and -1 (1 is total positive correlation, 0 is no correlation, and -1 is total negative correlation). PCC values for plots of precipitation against several geomorphic indices (Fig. 9) suggest no significant correlation between the amount of 505 precipitation and the ruggedness of the landscape.
We also compared the distribution of geomorphic indices for different geological units. We made three sets of points located within the Sierra Madre de Chiapas. Each set of points corresponds to one of the main geological units of the Sierra Madre de Chiapas: the Paleogene terrigenous deposits, the Cretaceous limestones and the Permian granodiorite of the Chiapas Massif (Fig. 2). For each nal profiles (i.e., knickpoints or convex segments) is a suitable approach to explore these phenomena.

Tonalá Shear Zone and volcanic arc
We analyzed 36 rivers draining the western flank of the Chiapas Massif and the Central America volcanic arc (Fig. 11). We focused our analysis on the areas located between the main drainage divide and the alluvial Pacific coastal plain. Longitudinal profiles were extracted from the modeled 525 drainage network ( Fig. 11 and 12). The upper knickzone extends a few kilometers upstream of the fault which makes the junction between the Tonalá Shear Zone and the Polochic fault system. However, the lower knickzone appears 540 more prominent and thus seems to be associated to a larger base-level change.
We extracted two parameters from analyzed longitudinal profiles in order to illustrate the alongstrike changes in the geometry of rivers. The first one is the steepness index (k sn ) extracted from logarithmic plots of slope against area (see examples in Fig. 12). Graded profiles show a single k sn value while non-equilibrated profiles are characterized by changing k sn values associated to several 545 segments. In this latest case we considered the k sn value of the lowermost segment as it is the one along which rivers adjust to new base-level conditions (e.g., Kirby and Whipple, 2001;Wobus et al., 2006;Whittaker et al., 2008). The second parameter is the estimated base-level change. If the profile displayed a well defined upper reach, we projected the segment located above the knickzone using a logarithmic plot of slope against distance (see equation 7). We then estimated the base-level change 550 according to the difference in elevations between the projected and the actual profiles at the outlet (see examples in Fig. 12).
Both k sn values and estimated base-level changes were plotted along a profile passing through the outlets of the analyzed catchments (Fig. 11). The results suggest a substantial alongstrike increase in k sn values, from ∼100 along the northern tip of the Chiapas Massif to >200 along the volcanic 555 arc. Estimated base-level changes seem to be low for the Chiapas Massif (<200 m) but increase significantly along the southern tip of the Tonalá Shear Zone (where it connects to the E-trending Polochic fault) and are higher within the volcanic arc (1000 to 2000 m).

Grijalva river
The Grijalva river (profile 1 in Fig. 13) flows from the northern tip of the volcanic arc towards the dams (e.g., La Angostura, Chicoasén and Malpaso dams) and canyons (e.g., Sumidero), we can infer a major knickzone located between the Tuxtla fault and the Chicoasén dam. We attempted to infer the geometry of the Grijalva river profile using the segments upstream and downstream of the 570 flooded areas. The segments of the river upstream and downstream of La Angostura dammed lake could be fit by the same regression in a logarithmic plot of slope vs area (see profile 1 in Fig. 13).
We are confident that no knickpoints are hidden below this dammed lake otherwise it would have been witnessed by a change in gradient of the river). In contrast, the gradient of the river upstream and downstream of the Chicoasen dammed lake differ and can not be fit by a single regression. The 575 gradient of the Grijalva river is 200 m between the Tuxtla fault and the Chicoasén dam (∼30 km apart) while it is less than 50 m between La Angostura dam and the Tuxtla fault (∼60 km apart).
We thus suspect that the Chicoasen dammed lake hides a major flooded knickzone which marks a ∼350 m base-level fall between the central depression of Chiapas and the Tabasco coastal plain. For the Malpaso and Peñitas dam the gradient of the river became too low and the scatter in slope-area 580 data too high to make such tests.
The Seleguá river (profile 2 in Fig. 13) is a right bank tributary of the Grijalva river which flows from the southern tip of the Sierra de los Cuchumatanes into the central depression of Chiapas. A prominent knickzone is observed upstream of the Polochic fault. This knickzone separates an upper segment with a gentle gradient (k sn < 50) from a lower and much steeper (k sn > 300) segment.

585
The upper segment is associated to the relicts of the Mid-Miocene Mayan paleosurface which is preserved on top of the Sierra de los Cuchumatanes . Downstream of the Polochic fault the profile displays an overall convex shape related to numerous knickzones. The lowermost segment is located within the central depression of Chiapas and shows a gentle gradient.
Longitudinal profiles 3 to 5 (Fig. 13) are from left-bank tributaries located within the Chiapas   Fig. 14) suggest a 500 to 800 m base-level change along the Tuxtla fault scarp.

615
Longitudinal profiles 4 to 6 are located within the northern parts of the Sierra de Chiapas. They are complex, as rivers flow through the E-trending left-lateral faults and then through the NW to NNW-trending structures of the Chiapas fold-and-thrust belt. As a result, these profiles display several knickzones (e.g., profile 4 in Fig. 14). However, the upper knickzones are more prominent as they mark a sharp transition between the uppermost segments characterized by gentle gradients and low 620 k sn values (<100) and the lower segments which show steeper gradients and higher k sn values (between 150 and 250). Reconstructions of upper reaches suggest a 700 to 1100 m base-level change with respect to the Grijalva river (profiles 5 and 6 in Fig. 14) or the Tabasco coastal plain (profile 4 in Fig. 14).

Steepness index map 625
The analysis of river longitudinal profiles shows that major base-level changes (i.e., higher than 100 m) are associated to prominent knickzones and to a substantial increase in steepness index (k sn ) values. We propose thus to map the variations of k sn values within the Sierra Madre de Chiapas in order to assess the locations of main drainage network perturbations. We computed k sn values from extracted river profiles using a 5 km moving window. We then interpolated obtained values using a 630 natural neighbor (NN) method. The result is shown in Figure 15.
The overall distribution of k sn appears to be tectonically controlled, as highest values are found close to major tectonic features. In Chiapas highest k sn values (>150 in Fig. 15) are mainly dis-

Maya mountains
We analyzed the drainage network of the Maya Mountains using a map of isobases from rivers and several longitudinal profiles (Fig. 16). Longitudinal profiles 1 to 5 (Fig. 16) are part of the Belize river watershed, which drain the northern and western side of the Maya Mountains. All these profiles display prominent knickpoints which mark the transition between the upper reaches which 650 are usually gently concave and the lower segments which are either concave with a steeper gradient or convex. The knickpoints are located upstream of the ENE-to E-trending Northern Boundary fault for profiles 1 to 3 (Fig. 16) and upstream of a NE-trending fault zone in profiles 4 and 5 (Fig. 16).
The reconstructed profiles made using the upper reaches suggest that the drainage network faced a 300 to 500 m base level change.

655
Profiles 6 to 9 (Fig. 16) are located in the southern and eastern parts of the Maya Mountains and flow directly into the Caribbean sea. Only southernmost profiles display prominent knickpoints which delimitate a well defined upper reach (e.g., profile 6 in Fig. 16). These knickpoints are located upstream of the NE-trending Southern Boundary fault and a ENE-trending fault zone. By contrast, most of profiles located along the eastern side of the Maya Mountains display an almost smooth and 660 concave shape such as profile 7 (Fig. 16). In particular, there is no significant change in the gradient of the rivers along the Southern Boundary fault. However, in some cases a knickpoint associated to a relict upper reach is found close to the drainage divide (profiles 8 and 9 in Fig. 16). The reconstructed upper reaches in profiles 5 and 8 (Fig. 16) suggest a 500 m base level change.
Previous works indicate the occurrence of such relict landscapes along the Polochic-Motagua fault system and within the Chortis block (western corner of the Caribbean plate). Rogers et al. the Ixcán and Motagua faults (Fig. 17) and was subsequently uplifted and deformed. Our results

695
suggest that relict landscapes are also preserved within the Sierra Madre de Chiapas and Maya Mountains. Topographic profiles (Fig. 5 and 6) and both relief anomaly and surface index maps Our results suggest a limited effect of climate. Within the Sierra Madre de Chiapas we found no clear correlation between peaks in precipitation and the distribution of rugged topography. There is no scaling between the amount of precipitation and the ruggedness of the topography or the amount of incision (Fig. 9). We observe a clear orographic effect on precipitation (swath profiles in Fig 6).

710
This effect is confined to the borders of the orogen and may locally be superimposed to a rugged topography. However the incision by the drainage network extends more inward where precipitation are no significantly different from the surrounding coastal plains. The erosion is related to the uplift of the Sierra Madre de Chiapas rather than to a change in climatic conditions. Indeed, recent studies suggest a dominance of tectonics over climate in the denudation of active mountain ranges (e.g., 715 Godard Fuchs et al., 2015). There is a possible effects of lithology on the location of incised areas and low relief surfaces. The distribution of morphometric indices suggests slight changes related to lithology (Fig. 10) We attempted to map the distribution of elevated relict landscapes from northern Central America (Fig. 17). The surface index and relief anomaly maps (Fig. 8) highlights well areas where elevated and low relief surfaces are extensively preserved and provide thus a good basis for mapping the extend of relict landscapes. In some areas where elevated surfaces are incised by the drainage network and for the Chiapas Massif where surfaces are nested between ridges, we refined our map using slopes and topographic contours derived from the SRTM data.
Relict landscapes from the Sierra Madre de Chiapas can be divided in four main domains. The 730 westernmost domain consist in a monadnock-type landscape which developed over the Chiapas Massif (Fig. 17). In morphometric maps, this area is characterized by low hypsometric integral and surface roughness values (Fig. 7). In fact, these geomorphic indices respond to an extensive system of large and relatively flat valleys. The drainage system appears in steady-state as rivers display smooth concave profiles and low k sn values (upper segments of profiles 4 and 5 in Fig. 13). The isobase sur-  However, this isobase surface is now located ∼400 m above the present-day sea level. We suspect thus that the topography of the Chiapas Massif was much lower than it is today and that it has been subsequently uplifted.
The large valleys of the Chiapas Massif connects to the east to two morphological domain: the Central Valley of Chiapas and La Venta plateau (Fig. 17). The eastern boundary of these two domains The last domain encompasses the elevated plateaus east of the Tuxtla fault (Comitán High and Sierra de Chiapas, Fig. 17). The northwest part of this plateau is heavily incised as shown by high surface roughness values (Fig. 7). By contrast, the areas south of the Tectapan fault appear to be less 760 dissected. In spite of elevations up to 2400 m, surface roughness values (Fig. 7) remains low. Indeed, the top of the plateau is associated to a low relief topography which is disrupted by incisions along the E-trending faults. These differences are probably related to several factors, including lithology (the southern areas consist mainly in Mesozoic carbonates with subterranean drainage limiting surface erosion), climate (the northern areas may receive more precipitation from the Gulf of Mexico), and 765 local base-level changes. The analysis of river longitudinal profiles (Fig. 13) suggests a ∼500-700 m base level change along the Tuxtla fault and ∼700-1100 m in the northern part of the Sierra de Chiapas.
The age of the relict landscapes in Sierra Madre de Chiapas is unknown. Thermochronological and stratigraphic evidence suggest two main periods of topographic growth. The first one occurred dur- Indeed, Witt et al. (2012a) show that the middle to late Miocene period is associated to marked increase in erosion of the Chiapas Massif and sedimentation within the basins of the Tabasco region.
The distribution of morphometric indices (Fig. 7) and the location of knickzones upstream of major faults such as the Tuxtla fault (profiles 1 to 3 in Fig. 14) indicate that this late Miocene erosional 780 topography was subsequently uplifted and fragmented by tectonics. This probably happened during the latest period of deformation (∼5 Ma to present). A middle Miocene age for the erosional topography of the Sierra Madre de Chiapas would be consistent with the timing of similar topographic features, such as the so-called "Mayan paleosurface" which covered most of the Polochic-Motagua sliver Brocard et al., 2011). According to Brocard et al. (2011), this 785 widespread planation surface formed at low elevation during the middle Miocene (between 13 and 7 Ma).
Our results show that most of the Maya Mountains consist in an elevated relict landscape (Fig. 17).
Geomorphic maps ( Fig. 7 and 8) as well as isobases from rivers (Fig. 16) indicate that low relief surfaces have been extensively preserved in the central and western parts of the Maya Mountains, 790 while the eastern part of the range has been more eroded. Most of the analyzed river longitudinal profiles (Fig. 16) display an upper reach associated to these low relief surfaces. Here again, the age of this relict landscape is unknown. According to Purdy et al. (2003), uplift in the Maya Mountains did not occur until late Neogene and perhaps no earlier than late Pliocene. The low relief surfaces we observe in the Maya Mountains may have formed during the Middle Miocene as the Mayan 795 paleosurface. In the southern part of the range, the low relief surface of the Maya Mountains gently dip towards the west and seems to connect to the flat topography of the Petén basin.

Morpho-tectonic interpretations
Landscapes from the Sierra de Chiapas and the Maya Mountains were fragmented and locally rejuvenated by tectonics. We propose to combine our topographic profiles ( Fig. 5 and 6), geomorphic 800 maps ( Fig. 7 and 8) and analyzed river profiles (Fig. 11 to 16) with published GPS and seismotectonic data in order to locate main active deformation zones. We attempted to define coherent tectonic blocks (Fig. 18) in order to produce a comprehensive map of the plate boundary.
Under constant climate and lithologies, the steepness index (k sn ) correlates with uplift rates (e.g., Kirby and Whipple, 2001;Wobus et al., 2006;Kirby and Whipple, 2012). We suspect thus higher 820 uplift rates along the southern border of the Chiapas Massif (north of the Polochic fault and south of the Necta fault). The Necta fault separates the northern domain, where the topography is almost in steady-state, from the southern domain, where landscapes are being rejuvenated due to an increase in uplift rates. We propose to define a rigid block (referred as 'Chiapas Massif sliver" in Fig. 18b) delimited by the Tonalá shear zone to the west, by the La Venta and Tuxtla faults to the east and by 825 the Necta fault to the south. The Tonalá shear zone seems to accommodate a 2.5 mm.yr −1 extension while the Tuxtla fault instead accommodates a 2.5 mm.yr −1 shortening (Fig. 18b).
The northern tip of the Central America volcanic arc forms an uplifted plateau with a mean elevation of ∼3000 m. The drainage network developed prominent knickzones as it responded to high-amplitude base level changes ( Fig. 11 and 12 and k sn map in Fig. 15). The differential mo-830 tion between GPS stations TPCH and SOL (Fig. 18) indicates a ∼5.8 mm.yr −1 NNE shortening between the Pacific coastal plain and the Sierra de los Cuchumatanes. We also noted that the motion of GPS stations located along the volcanic arc (MAZ, CHL and SSIA in Fig. 18) is intermediate between the Cocos forearc sliver and the Caribbean plate. From the motion of these GPS stations we infer a dextral fault zone bounding the southern edge of the volcanic arc. We suspect that the 835 volcanic arc acts as a tectonic sliver which accommodates a significant part of the motion between the Cocos forearc sliver and the overriding plates. East of the Graben of Guatemala City, the sliver is possibly delimited by two right-lateral fault zones which accommodates each half of the motion (∼13 mm.yr −1 ) between the Cocos forearc sliver and the Caribbean plate. West of the Graben of Guatemala City, the volcanic arc sliver is pinned between the Cocos forearc sliver to the SW, the 840 Polochic sliver to the NE and the Chiapas Massif to the north (Fig. 18b). Though a 6.3 mm.yr −1 shortening is expected between GPS stations MAZ and JOY, we could not find a clear limit between this volcanic arc sliver and the Motagua sliver.
Within the Sierra Madre de Chiapas orogenic belt, most of the present day seismicity (Fig. 3b) is located along the strike-slip faults of Chiapas, which delimitate the La Venta block ( (García-Acosta and Suárez-Reynoso, 1996;Peraldo and Montero, 1999;Guzmán-Speziale, 2010) and the 23 September 1902 (Böse, 1903;Figueroa, 1973) earthquakes. Historical earthquakes of this 850 area were attributed either to the Concordia (Guzmán-Speziale, 2010) or to the Tuxtla (Andreani et al., 2008a) faults. The differential motion between GPS stations CONC and SOLE suggests a ∼2.5 mm.yr −1 NE-trending compression which is perpendicular to the azimuth of mapped faults.
Our geomorphic data suggest that the Tuxtla fault is more likely to have produced these earthquakes.
The Concordia fault trace is not associated to any major fault scarps and rivers display a rather 855 low degree of perturbation. By contrast the Tuxtla fault scarp displays a very young and prominent morphology (high hypsometry and low surface roughness) and the drainage network responded to a significant base-level fall (>500 m, profiles 1 to 3 in Fig. 14). Moreover, major topographic uplift took place east of the Tuxtla fault (Comitán High and the Sierra de Chiapas) where the elevation of relict landscapes is between 1600 and 2400 m.

Proposed model of the plate boundary 870
The proposed model for the North American-Caribbean-Cocos plate boundary is displayed in Fig. 19.
This model incorporate main ideas from previous works. East of 91 • W the dynamic of the plate boundary is guided by the eastward escape of the Caribbean plate. This motion results in sinistral shear along the Polochic and Motagua slivers, dextral shear along the Jalpatagua fault and extension mainly concentrated along the grabens of Guatemala City and Ipala ( Fig. 18b and 19; e.g., Lyon-Caen et al., 2006;Authemayou et al., 2011). Recent models of the triple junction also agree on the fact that both North American and Caribbean plates are limited by a rigid forearc sliver to the west (Turner et al., 2007;Phipps-Morgan et al., 2008;Authemayou et al., 2011). The eastward escape of the Caribbean plate results in a counterclockwise rotation of the rigid forearc sliver (Phipps-Morgan et al., 2008;Authemayou et al., 2011).

880
The interactions between the forearc sliver and the North American plate were summarized as distributed compression north of the Polochic fault (Guzmán-Speziale and Meneses-Rocha, 2000; Andreani et al., 2008a;Witt et al., 2012b). Our interpretation, based on the characterization of elevated relict landscapes and published geophysical data, is that the compression between the forearc sliver and the North American is mainly absorbed at the boundary of small tectonic blocks which 885 act as 'buffers' (Fig. 18b). As a result, significant topographic uplift occurred along two narrow domains. The first one encompasses the Sierra de Chiapas and the western tip of the Polochic sliver (Sierra de los Cuchumatanes) and the second one correspond to the northern tip of the volcanic arc ( Fig. 19).
Within the Sierra Madre de Chiapas, the distribution of highest relict landscapes (>1500 m) in-  (Manea and Manea, 2006;Manea et al., 2013). This flattening would have resulted in higher coupling along the subduction interface, as suggested by Franco et al. (2012). As a result, the areas west of the Tuxtla fault (which is located along most of the present day volcanic centers) became partially incorporated into the forearc domain. This would explain why significant Plio-Quaternary topographic growth (> 905 1000 m) mainly occurred east of the Tuxtla fault rather than along the Tonalá Shear Zone.
The compression between the rotating forearc sliver and the North American plate resulted in topographic uplift along the Sierra de los Cuchumatanes and the northern tip of the Central America volcanic arc (Authemayou et al., 2011, this study). In existing models, the westernmost corner of the Caribbean plate is stretched along NE-trending normal faults and the counterclockwise rotation 910 of the forearc sliver is associated to a zipping (i.e. a suture zone) between the inferred prolongation of the Jalpatagua and Motagua faults Franco et al., 2012). However, a suture zone between the Jalpatagua and Motagua faults may not fully explain why the volcanic arc is uplifted, mainly because most of the elevated Volcanic arc plateau is located west of the area where both faults would be zipped. To explain the atypical motions of GPS stations along the volcanic arc 915 and the topographic elevation of the plateau (MAZ, CHL and SSIA in Fig. 18a) we infer a fault zone bounding the western edge of the volcanic arc ( Fig. 18 and 19) and Mountains. This residual motion may have reactivated NE-trending faults inherited from the Eocene 930 opening of the Yucatán basin (Rosencrantz, 1990;Leroy et al., 2000).

Conclusions
We demonstrate that a geomorphic analysis allows to classify zones of similar relief patterns that we assume witnessed different tectonic and erosive histories. Using DEM-based geomorphic indices, we examined the topography of the Sierra Madre de Chiapas and the Maya Mountains. We 935 used topographic profiles and morphometric maps in order to understand the spatial distribution of landscapes and we also analyzed in detail the disequilibrium of drainage network in order to map potential vertical displacements. Finally, we combined our results with existing GPS and seismological data in order to better understand the interactions between tectonics and landscapes within the highly diffuse North American-Caribbean-Cocos plate boundary.  Fig. 18b). Within the Sierra de Chiapas and Comitán High, relict landscapes are found at elevations between 1600 and 2400 m. The areas between the Tuxtla and La Venta faults to the east and the Tonalá shear zone to the west were also 960 uplifted but to a lesser extend (400 to 600 m). Major uplift is also observed along the western tip of the Polochic sliver (Sierra de los Cuchumatanes) and the northern tip of the volcanic arc. Using GPS data and morphotectonic analyses, we propose that the northern tip of the Central America volcanic arc is part of a tectonic sliver which has been pinned between the Cocos forearc sliver and the North American plate (Fig. 18b)     . Comparison between precipitation(derived from TRMM data) and morphometric indices within the Sierra Madre de Chiapas. Points represents pixel values extracted from the precipitation (Fig. 5b) and morphometric maps ( Fig. 7 and 8). The linear correlation between the variables was estimated using the Pearson correlation coefficient (PCC).         (Franco et al., 2012). Colors emphasize vectors with similar azimuths and velocities.
The diagram of northern and eastern components (bottom left corner) allows to deduce motions between paired GPS stations. B: Main tectonic blocks defined using known major faults (Fig. 2), our geomorphic maps ( Fig. 7 and 8) and published GPS data (Fig. 18a) and seismicity (Fig. 3a). Black arrows show the motion along blocks boundaries deduced from paired GPS stations (values are in mm.yr −1 ). See Fig. 2 for faults abbreviations.