Drainage rearrangement as a driver of geomorphological evolution during the Upper Pleistocene in a small tropical basin

Abstract


Introduction
The Ocoa River drainage basin (670 km 2 ) is located in the southern hillside area of the Central Range ("Cordillera Central") of the Dominican Republic, Hispaniola Island (see figure 1-a).The drainage network is a complex system with two main rivers, the southeast-oriented Banilejo River (which has the longer course, 84 km), and the southoriented Ocoa River (see figure 1-b).The two rivers join together near the center of the basin, from where they flow southward as a single stream to meet the Caribbean Sea.
The Ocoa drainage network flows mainly over outcrops of the so-called Peralta Belt (figure 2-a), which is a syntectonic sedimentary basin made up of rocks of the Cenozoic Era (Hernaiz Huerta, 2000;Hernaiz Huerta & Pérez-Estaun, 2002;Hernaiz Huerta et al., 2007;Pérez-Estaún et al., 2010;Mollat, Wagner, Cepek, & Weiss, 2004).The northeastern part of the Ocoa drainage basin consists of igneous and sedimentary materials of late Cretaceous age, which includes part of the volcanic Arroyo Jigüey and Tireo Formations, as well as dark limestone and conglomerate outcrops from El Manaclar Fm. and Los Martínez Fm., respectively (Servicio Geológico Nacional, n.d.).Moreover, the whole complex, especially the materials of the Peralta Belt, has experienced folding, thrusting, overthrusting, and wrenching during the Neogene Period, as a result of the indentation of the Beata Ridge (a submarine ridge that extends south-southwest from Cape Beata into the Caribbean Island Arc of Hispaniola) (Fox, Ruddiman, Ryan, & Heezen, 1970;Mann, McLaughlin, & Cooper, 1991;Hernaiz Huerta & Pérez-Estaun, 2002).
River basins affected by intense tectonic activity, such as the Ocoa Basin, usually harbor excellent examples of drainage rearrangement (Molin & Fubelli, 2005;Ascione, Cinque, Miccadei, Villani, & Berti, 2008;Lavarini, Magalhães Júnior, Oliveira, & Carvalho, 2016;Antón, De Vicente, Muñoz-Martín, & Stokes, 2014;Piacentini & Miccadei, 2014).Bishop (1995) defined drainage rearrangement as a transfer of part or all of a river's flow to another river, a mechanism that may explain river history at different spatial scales.Such changes affect both sediment volume and provenance, which may also have biogeographical implications (Dingle & Hendry, 1984;Loureiro, Duarte, & Zarucki, 2011;Valdovinos et al., 2010).In his thorough review, Bishop defined three forms of drainage rearrangement on bedrock systems: river diversion, beheading, and the traditional concept of river capture.The author summarized the key attributes and processes responsible for the occurrence of all the forms of drainage rearrangement, based on both his own research and contributions from other authors (Davis, 1889;Thompson, 1939;Clark, 1989;Haworth & Ollier, 1992;Small, 1978).Specifically, for river capture (stream piracy), Bishop highlighted the key attributes defining a capturing stream, namely steepness and high energy, as well as an erosive advantage that produces undercutting in soft lithologies (e.g., excessive fluvial incision).Furthermore, tectonic tilting, channel migration, and catastrophic breaching are the processes responsible for river diversion, while scarp retreat is responsible for the beheading rearrangement form.Thus, it is clear that tectonics is an important driver of drainage rearrangement.Because of the influence of tectonics throughout the Ocoa Basin, all forms of rearrangement proposed by Bishop could theoretically exist in it; however, no systematic research has hitherto been conducted to uncover specific cases.
In addition, some authors discuss the role of karst as an understudied driver of drainage rearrangement.Hill and Polyak (2014) and Hill, Eberz, and Buecher (2008) demonstrated that the most suitable explanation as to how the Colorado River crossed the Kaibab uplift was that it occurred by a "karst piracy" process, as first proposed by Hunt involving what he called "piping", a term that is reserved for other uses in karstology (Ford & Williams, 1989).As stated by Hill and Polyak, the key stages of a complete karst piracy process are the following: establishment of a gradient across the topographic divide, followed by the settling of a cave passage that is capable of carrying a significant volume of water, ending with a roof collapse, which creates a sub-aerially exposed gorge that can be widened to form a canyon.In short, karst piracy is a type of drainage rearrangement that involves, under certain conditions, the settling of a cave passage, which is later sub-aerially exposed by roof collapse.Since limestone outcrops occur in the study area, karst piracy is 57 also considered in this research as a potential driver of drainage rearrangement.
The idea for the present study came after a reconnaissance field trip in 2014 to the Parra River Basin, which is part of the Ocoa River Basin.During this trip, a striking variation of alluvial deposits across geomorphological units was detected.Specifically, high terraces consisted mainly of large clasts of sedimentary rock type, while the main channel beds featured small clasts of evenly represented sedimentary and igneous types.Subsequently, supplementary field observations were conducted, and morphometric parameters were computed for the sub-basins of the Parra Basin.These analyses revealed a high variability of both basin topography and drainage network attributes across sub-basins.
Given the current knowledge on the geology and geomorphology of the area, the most plausible explanations considered for such variability were the occurrence of drainage rearrangement; thus, two candidate hypotheses were evaluated.The first hypothesis proposed one or several classical drainage rearrangements affecting specific streams of the network in candidate wind gaps.The second hypothesis proposed a karst piracy model affecting a tributary of the Parra River, along with other minor processes of the classical Bishop approach.The latter provided the most plausible explanation for the evidence collected.Hence, the aim of this paper is to examine data from alluvial deposits, geomorphological features, and basin morphometry parameters to support a suitable explanation for the studied phenomenon of drainage rearrangement.

Study Area
The Parra River Basin, with a surface area of 29.5 km 2 , is located in the eastern half of the Ocoa Basin (see figure 1-b).This basin is made up of igneous and sedimentary outcrops of late Cretaceous age oriented in a northwesterly direction (see figure 2-b), comprising basalt and pyroclastic rocks (Arroyo Jigüey Fm.), limestone, marl, and conglomerate (El Manaclar and Los Martínez Fm.), as well as intermixed volcanic outcrops with limestone and chert layers.
For this research, the outermost polygon of the Parra Basin was conventionally split into three smaller sub-basins, named after their main streams as Naranjal (5.67 km 2 ), Hondo (10 km 2 ), and Parra itself (11.99 km 2 ) (see figure 1-c).
In addition, the remaining small area (1.84 km 2 ) close to the confluence with the Ocoa River (hereafter the Outlet Zone), characterized by a narrow central valley with steep slopes that drains directly to the Parra River, was included in the study.
Rock composition and tectonics vary significantly in the study area.The Naranjal Basin is small in size and elongated in shape, with a wide fluvial valley filled with Quaternary deposits occupying the central area.The area surrounding its central valley is made of Cretaceous materials, predominantly comprising both intermixed volcanic rocks and karstified limestone layers at the outlet of the sub-basin, which are of particular interest for the drainage rearrangement studied in this paper.The western divide of the Naranjal Basin is low and gently sloping, while the eastern side is steeper and higher than the former.The Hondo Basin is poly-lobed and pear-shaped, with a narrow fluvial valley surrounded by steep slopes made of volcanic rocks.Finally, the Parra Basin is elongated with high watershed divides mainly made of sedimentary rocks, and features a V-shaped valley surrounded by steep slopes, which extends also into the Outlet Zone.
Normal and thrust fault lines intersecting the study area are oriented in east/west, northeast/southwest, or northwest/southeast directions (Servicio Geológico Nacional, n.d.;Mollat et al., 2004).A remarkable feature is the ∼7 km-long thrust fault line oriented in a northwest/southeast direction almost parallel to the main channel of the Parra Basin, hereafter the Parra thrust fault (see figure 2 for a simplified geological map).

Data Sources and Methods
Multiple approaches are used to study river rearrangement, the most common being the analysis of both network and basin morphometry parameters (Molin & Fubelli, 2005;Ascione et al., 2008;Antón et al., 2014).In addition, analysis of rock type composition is a useful tool, because it allows examination of the sediment sources corresponding to the inferred pre-and post-drainage rearrangement of a given basin (Bishop, 1995).Thus, the necessary evidence for discussing explanatory hypotheses of the studied phenomenon came from analysis of the composition of selected alluvial deposits, description of geomorphological features,, and computation of basin morphometry parameters.
The analysis of alluvial deposits was undertaken to reveal the rock type composition of coarse gravel-and cobble-sized particles in selected geomorphological units.To this end, six samples were collected and labeled S1, S2, and so on, up to S6 (see figure 2-b and figure 3).Samples S1 to S4 were aimed at classifying 330 randomly picked clasts into three major rock types: igneous, sedimentary, and metamorphic.In particular, samples S1 and S2 were collected in the Outlet Zone, the former from a low terrace (+2 m above the Parra River level), the latter from a coarse-grained bed point of the mainstream.Samples S3 and S4 were collected from coarse-grained stream beds as well, but at the outlet areas of the Hondo and Naranjal Rivers, respectively.Samples S5 and S6 were taken for the purpose of detecting the presence of igneous types by randomly picking at least 30 clasts each from scarps of high alluvial terraces in the Outlet Zone, a challenging task owing to problems with accessibility.Specifically, sample S5 was collected at the mid-level of a +10/12 m-high terrace close to the outlet of the basin, and S6 was collected almost at the summit of a +75/80 m-high terrace, located near the Parra Village (heights relative to the local base level).It is worth highlighting that, even though the samples collected from alluvial deposits may not adequately represent the entire basin, the results obtained from analyzing them were consistent enough to support the geomorphological evolution proposed.The morphometric parameters computed, which are summarized in table 1, comprised network and basin properties.All the computations were performed using a semiautomatic processing workflow in R (R Core Team, 2018) and GRASS GIS (GRASS Development Team, 2017), with the package rgrass7, which is an interface between the two programming environments (Bivand, 2018).Furthermore, QGIS was used as a graphical interface to display results (QGIS Development Team, 2019).A summary of the processing workflow is described in the following.Horton (1945) Circularity ratio (area of the basin divided by the area of a circle having the same perimeter of the basin) Miller's formula as cited in Gray (1961) Compactness coefficient (perimeter of the basin divided by the diameter of the circle with the same area of the basin) Gravelius's formula as cited in Bendjoudi and Hubert (2002) Profile concavity Snow and Slingerland (1987) Hypsometric integral, numerically integrated area under the hypsometric curve Giandotti's formula as cited in Grimaldi et al. (2012) Transverse topographic symmetry vector and probability of obtaining a greater mean vector magnitude by pure chance combination of random vectors Cox (1994); Curray (1956) Legend: A 1 = numerically integrated area that lies between the profile curve and a straight line connecting the profile endpoints; A 2 = triangular area below that straight line and above a horizontal axis connecting with the profile's downstream endpoint; A c = area of a circle with equal perimeter as the basin; a i = area of the trapezoid i under the dimensionless height-area curve (hypsometric curve); d a = distance from the active meander-belt midline to the basin midline; d d = distance from the basin divide to the basin midline; e = base of the natural logarithm; l i = length of the stream segment i; L m = length of the main channel (in km); h = difference between the mean basin elevation and the outlet elevation (in m); n = number of basin segments measured; N ω number of segments of order ω; P = perimeter of the basin; T = mean vector magnitude of the transverse topographic symmetry T (in percent); ∆z i = elevation difference of segment i; Φ = topological diameter; ω = stream order; † A in km 2 .
First, the source data used as input were from a SRTM 1 Arc-Second digital elevation model (NASA LP DAAC, 2000), hereafter DEM, which was downloaded from the EarthExplorer site (https://earthexplorer.usgs.gov).
The DEM was warped onto a 30 m-resolution grid of the UTM projection (WGS84 Datum) using the GDAL library (GDAL Development Team, 2017).Subsequently, the DEM was pre-treated using the denoising algorithm designed by Sun, Rosin, Martin, and Langbein (2007), configured with the parameters suggested by Stevenson, Sun, and Mitchell (2010).Afterwards, sinks were removed using a hydrological conditioning algorithm based on a minimum impact approach (Lindsay & Creed, 2005).Visual inspection and statistical comparisons with respect to the original DEM were also performed, which suggested that the optimized DEM was suitable for the extraction of parameters.Finally, multiple GRASS GIS add-ons were used to compute the morphometric parameters, e.g.r.watershed, r.stream.*and r.basin (Ehlschlaeger, 1989;Metz, Mitasova, & Harmon, 2011;Jasiewicz & Metz, 2011).
It is worth noting two important configuration details of the algorithms used for computing the stream network.First, the surface flow accumulation and the drainage direction rasters were derived from the optimized DEM using the r.watershed add-on configured in multiple flow direction (MFD) mode (Ehlschlaeger, 1989;Holmgren, 1994).Second, several thresholds of upstream area accumulation for the stream extraction algorithm were tested, from which the value "22 cells" (19800 m 2 ) was considered a suitable option.Although the network generated with such a value showed a higher density of streams than that generated in similar studies (e.g., as in Ozulu & Gökgöz, 2018), it was considered a suitable network for the purpose of this research, since the Parra Basin is relatively small (29.5 km 2 ).In addition, particular care was taken to prevent the generation of artifacts, so orthorectified aerial photos were used to visually confirm that the automatically extracted streams matched those of the actual network (Martínez Batlle, 2018).
A novel approach using the longest flow paths (hereafter LFPs) of the sub-basins draining to the main channel of each basin was implemented to characterize the drainage network.For each LFP, elevation drop, planimetric length (or simply length), and profile concavity (see table 1) were computed using R and GRASS GIS scripts, for which the following procedure was applied.First, each variable was normalized using the formula z i = x i −min(x) max(x)−min(x) , where x = (x 1 , . . ., x n ) and z i is the corresponding normalized value of x i , so that the newly calculated values are laid in the interval [0,1] (Gower, 1971;Pavoine, Vallet, Dufour, Gachet, & Daniel, 2009).The normalization was intended both to improve the performance of further clustering steps and to assign equal weight to each variable.Afterwards, the normalized values were used to compute dissimilarity matrices of euclidean distance across the different LFPs.For testing purposes, several clustering methods were used to generate dendrograms, including the unweighted pair group method with an arithmetic mean (UPGMA), Ward's minimum variance, single linkage, and complete linkage.In order to estimate the degree of dissimilarity information retained by each dendrogram, cophenetic distance matrices were computed and compared with the original dissimilarity matrix using Pearson's r correlation (i.e., cophenetic correlation) (Borcard, Gillet, & Legendre, 2018).Subsequently, classifications with k and k + 1 clusters were generated from the dendrograms, in which k was chosen following the criteria of maximum average silhouette width (Borcard et al., 2018).Finally, the classification generated from the dendrogram with the highest cophenetic correlation, which at the same time produced a meaningful number of clusters (e.g., without single-element clusters), was chosen as optimal for separating and characterizing the LFP of each basin.

Alluvial deposits composition
Composition by rock types in alluvial deposits varied significantly across samples S1 to S4 (see figure 2-b for location map, and figure 3 for context view photos).Both the overrepresentation of sedimentary rock type and tbe scarcity of igneous rock type in S1 were responsible for the significant variation found (see table 2 and figure 4).Furthermore, a closer look at the streambed samples alone revealed a quite homogeneous pattern of rock type composition.Moreover, samples S5 and S6 consisted of coarse gravel and boulders of the sedimentary type; no clasts of the igneous type were found.Therefore, when pooling all the samples together, the terraces (S1, S5, and S6) showed a greater dominance of sedimentary type clasts, in contrast with those of the streambeds (S2 to S4).It is worth noting that increasing the sampling effort in terraces may help in finding clasts of the igneous type, but in any case the number would be small.

Morphometric parameters of Parra Subbasins
The areal traits, as well as the topographic and network parameters of the three basins studied, are summarized in table 3. The surface area and the set of topographic parameters computed (e.g., H and S) show that both the Parra and Hondo Basins are medium-sized, mountainous, and highly dissected, while the Naranjal Basin is small, with a rolling wide valley at the center.Moreover, the Parra and Hondo drainage networks reached a maximum Strahler order of 4, while the Naranjal network reached an order of 3. The Hondo and Parra drainage networks obtained roughly similar bifurcation ratios (∼ 3.85), in contrast with the Naranjal network, which yielded an unexpectedly high bifurcation ratio (> 5).In summary, table 3 shows that the Hondo and Parra networks are very similar in terms of topographical and network parameters, while the Naranjal network is quite different from both.An exception to this generalization is that all the basins have similar values of drainage density.The LFPs draining to the Hondo Basin were classified in three groups using a UPGMA cluster analysis, based on profile concavity, elevation drop, and stream length.The groups are characterized as follows (see figure 5): Group A, short and straight profile LFPs, having a small elevation drop, mostly located to the north of the main channel; Group B, medium-sized and concave profile LFPs with a moderate elevation drop, distributed alternately on both sides of the main channel; Group C, three large slightly concave profile LFPs with high elevation drop, located exclusively at the south-center of the basin.Overall, the LFPs draining to the Hondo Basin show both varied length and elevation drop, as well as moderate profile concavity (C As = 0.10).The LFPs draining to the Naranjal Basin are dominantly concave, with a mean concavity of 0.28 (see figure 6 and table 4).Two groups were distinguished, namely A and B. Group A is composed of large LFPs with a high elevation drop, whereas group B is characterized by short LFPs with a small elevation drop.In both groups, the LFPs showed pronounced concave profiles.Finally, a conspicuous steep slope forms the downstream ending of the main channel (c.40 m of elevation drop).This slope is made of a tufa deposit, with an active waterfall during the humid season, located within the remnants of a sinkhole (figure 6).Seen in plan view, the sinkhole has an elliptical shape of 115 × 95 m, in major and minor axes, respectively, and 0.85 Ha in area.It is floored with large limestone blocks (> 5 m in width) from El Manaclar Formation, which in most cases are covered with tufa concretions (figure 7).
The LFP profiles draining to the Parra main channel exhibit a predominantly rectilinear pattern (see figure 8), with a mean concavity of 0.02 (table 4).Three groups of LFP were clearly distinguished: Group A consists of LFPs with both short length and small elevation drop, as well as straight profile, mainly located in the south portion of the basin and draining to the main channel from both margins; Group B is composed of medium-length LFPs with varied profile concavity, as well as a large elevation drop, and distributed all over the basin along both sides of the main channel; Group C is made up of two large LFPs with concave profile and high elevation drop, one on each margin of the main channel, which meet it close to the center of the basin.
The hypsometric curve and its numeric counterpart, the hypsometric integral, shown in figure 9, summarize the cumulative distribution of elevations across the basins analyzed.Both the Parra and Hondo Basins have slightly concave to moderately concave curves, as well as average values of the hypsometric integral, the Parra Basin having the highest of the ensemble (I = 0.39).Furthermore, the Naranjal Basin showed a well-defined concave hypsometric curve with the lowest integral value (I = 0.24).The Naranjal hypsometric curve also displays a sharp vertical section at the lower end, which corresponds to a short path preceded by a steep elevation drop.The significance of this feature in relation to this research is discussed in detail in the next section.
The transverse topographic symmetry vectors of the studied basins, as well as the corresponding polar plots, presented in figures 10 and 11 respectively, show that the active meander-belts of both the Hondo and Parra Rivers have significantly moved away from their corresponding basin midlines, while the Naranjal River meander-belt shows a migration pattern that is not significant.The greatest magnitude of migration in the ensemble is observed in the Hondo meanderbelt, with a mean symmetry vector of 0.29, 329 • .Also, migration vectors with magnitude 0.20 or greater are distributed throughout the whole Hondo meander-belt.Furthermore, the Parra meander-belt shows a consistent migration pattern, with a mean symmetry vector of 0.24 oriented to the northeast (54 • ).In addition, the Parra sub-basin shows a cluster of vectors of 0.20 or greater magnitude, which are concentrated in the catchment area, at the southwest quarter of the basin.Finally, the Naranjal meander-belt has migrated randomly either to the west-northwest (mainly in the center of the basin) or to the east (mainly in the north portion), with a mean symmetry vector of 0.24, 356 • .

Discussion
This section is divided into two parts.In the first part, morphometric parameters, composition of alluvial deposits, and geomorphological features are analyzed as supportive evidence for an evolution of the Naranjal drainage network, which includes a karst piracy model at the final stage of the sequence (Hill & Polyak, 2014).In the second part,

The Naranjal River Basin
First, it may be noted that the circularity ratio and the compactness coefficient reflect the elongated nature of the Naranjal Basin.Early evolution of the basin was likely to be controlled by a set of faults oriented in a northeastsouthwest direction (see figure 10), which also favored the elongated shape.In addition, the high concavity indexes, the low hypsometric integral of the basin, and the extremely high bifurcation ratio of the network (R b > 5) suggest that the Naranjal Basin experienced a long-term evolution, initially at a high uplift rate, along with fluvial incision, and subsequently a basin uplift at a gradual and steady rate, along with alluvial aggradation (Snow & Slingerland, 1987;Strahler, 1952;Horton, 1945).Also, the wide fluvial valley at the center of the basin, filled with Quaternary deposits, supports the proposed long-term evolution model.Moreover, the co-existence of two well-differentiated groups of stream paths in the network (see figure 6) is interpreted as evidence for the proposed evolution at two consecutive rates of uplift.
At a later stage, the base level might have experienced steady and gradual drawdown with respect to the central fluvial valley, and both a cave passage and a ponor were developed in El Manaclar Fm. limestone (southernmost border of the Naranjal basin, see figure 12-a).This stage may be considered as the beginning of the karst piracy process.As long as the cave passage existed, it was still intricate enough to prevent coarse-grained bedload to pass, so only water and fine-grained sediments could reach the Parra River, in a way very similar to that Hill and Polyak proposed for the Colorado River.This assumption supports the conspicuous absence of igneous clasts from Arroyo Jigüey Fm. in the terrace deposits of the Outlet Zone.However, using the available evidence, it is difficult to state with certainty whether or not the Naranjal River flowed through an alternate sub-aerially exposed channel before the establishment of the cave passage.More research is needed to clarify this issue.
Afterwards, the cave passage might have expanded laterally and the limestone roof collapsed over the underground river, which developed the present sinkhole floored with large limestone blocks from El Manaclar Fm. (figure 7-b).
After the collapse, many of these blocks have been covered with tufa concretions deposited by the newly sub-aerially exposed stream of the Naranjal River (figure 7-c).During a visit in 2015 to the sinkhole, a local reported a welldeveloped cave, namely "Cueva de los indios", and described it as a vertical/sub-vertical entry passage of 6 m drop with a horizontal gallery at its lowest level.The only known entrance to this cave is an extremely narrow breach under a large block settled in the center of the sinkhole (shown in the background on figure 7-b).He also reported that the gallery had no permanent water circulation, even during the moist season, which suggests that the cave became part of the vadose zone of the karst, or maybe got clogged after the collapse.A speleological survey could help to determine whether this cave is a remnant of the original underground river or a passage developed after the collapse as a response to the steady drawdown of the base level.
Given that the limestone is a hard rock, the sinkhole itself and its surrounding area had undergone little fluvial incision since the collapse event and until the present.The same is true for the Naranjal River valley, which has finally evolved as a fluvial valley of fine-grained and gravel-sized deposits with no remarkably fluvial incision.As T obtained a not significant value, it might be interpreted that the migration pattern of the Naranjal meander-belt has been most recently controlled by internal fluvial processes and not by external ones (e.g., tectonic tilting) (Cox, 1994).
In order to include an absolute chronology in the geomorphological evolution, radiocarbon dating was performed on a series of carbon films covering gravel-sized clasts collected in the same terrace where sample S1 was taken (see figure 3-a).A conservative age between 30730 and 30220 calibrated years before the present (cal yr B.P.) was obtained.The stratigraphic features of this +2 m terrace suggests that gravel-sized clasts of sedimentary type (mainly marl, sandstone, and limestone), were deposited in a floodplain pond and subsequently covered with organic matter that later turned into the carbon films dated for this research.The scarcity of igneous rocks suggests that the formation of the terrace is correlative of a condition when the Naranjal River flowed through the cave passage, such that no complementary source of igneous clasts was delivered to the Parra River.Although more research is needed, this evidence ultimately suggests that the roof collapse of the Naranjal River cave passage might have occurred after ∼ 30 cal kyr B.P.
Finally, a negative test result of the evolution of the Naranjal drainage network is reported.A hypothesis about the former courses of the La Vaca and Naranjal Rivers within the same basin was also evaluated (see figure 1-c for location).Such a model was initially deemed possible since the gently sloping western divide of the Naranjal Basin hosts at least The filled circle depicts the approximate location of Parra Village (labeled as "PV").The black curves with arrowhead symbol denote elbows of river diversion.The label "M" points to a sharply bent meander in the Hondo River.The labels "WG" point to wind gaps on the western divide of the Naranjal Basin.The solid grey lines in (b) depict fault-lines extracted from the geological map at 1:50,000 scale, sheet 6171-IV, Servicio Geológico Nacional (n.d.).
three wind gaps (see figure 12-b), which were detected using photogrammetric techniques (Martínez Batlle, 2018).Specifically, the hypothesis proposed that either the La Vaca or the Naranjal River might have flowed to the other through one of the wind gaps, in a context in which both streams had a higher base level.Subsequently, after a hypothetical process of diversion, the Naranjal River might have started to flow through a path similar to the current one.The whole model was finally rejected, since no alluvial deposits were found in the wind gaps.

The Hondo and Parra River Basins
Several morphometric parameters, as well as sediment samples, suggest that both the Hondo and Parra Basins experienced faster uplift than did the Naranjal Basin.A geomorphological evolution of these basins must consider the influence of tectonics, since thrust and/or normal faults may have led to the drainage rearrangement by means of block tilting, lateral migration, or breaching (Bishop, 1995).Thus, apart from proposing a general evolution model for both the Parra and Hondo Basins using morphometric parameters, additional striking morphological features observed in the network are discussed as evidence of drainage rearrangement, for which the classical three forms proposed by Bishop are equally deemed possible.
The highly elongated shape of the Parra Basin, which is denoted by both a circularity ratio and compactness coefficient, seems to be strongly influenced by the Parra thrust fault.Also, the low concavity indexes of the LFP analyzed in the Parra network, as well as the nearly 1000-m elevation drop of the basin and the high hypsometric integral value, support the proposed model of rapid and continued uplift (see figure 8) (Snow & Slingerland, 1987;Strahler, 1952).Moreover, a significant value of T suggests that the main channel has been migrating to the northeast, which is likely to be driven by the Parra thrust fault as well (see figure 10 and figure 11).
Furthermore, the Hondo Basin may have experienced a rapid uplift initially, but at a lower rate than that of the Parra Basin, and it was likely more influenced simultaneously by both the Parra thrust fault and the transverse/oblique normal faults.A remarkable feature of this drainage network is the mixture of concavity traits (see figure 5), which may be attributed to the referred convergence of fault line systems.Such concavity traits, along with the hypsometric integral and topographical parameters (e.g., elevation drop and bifurcation ratio), suggest that the Hondo Basin has experienced a gradual and steady uplift.
Most of the morphometric parameters calculated are consistent with the geomorphological features observed in the eastern part of both the Hondo and Parra Basins; therefore, no drainage rearrangements are proposed in the evolution of their eastern half.However, the existence of two geomorphological striking features in the vicinity of the Outlet Zone suggests that at least one drainage rearrangement form may have occurred in the western half of the study area during the late stages of the evolution.A geomorphological interpretation of these features is detailed as follows.
The first feature analyzed is a sharp meander of the Parra River situated just downriver from the current confluence with the Hondo River.River diversion is the most plausible explanation in this case.The sequence proposed starts when an early mountainous drainage network was already established, and the paleo main channel of the Parra River generated an alluvial fan close to where Parra Village is presently sited (see figure 12-a).Samples taken in high terraces near Parra Village (see figure 3-e and f) seem to be correlative of the formation of this fan, since they consist mainly of coarse gravel and boulders of sedimentary rock type, which is also consistent with the large proportion of sedimentary outcrops observed in geological maps of the Parra sub-basin (Mollat et al., 2004;Servicio Geológico Nacional, n.d.).Topographic profiles generated from photogrammetric techniques (Martínez Batlle, 2018) show that Parra Village is likely situated over the remnants of the alluvial fan and also suggest that the Hondo River formed a small fan, which coalesced from the east with the larger one.
Following the formation of the alluvial fan, the Parra main channel may have started to divert to the northeast, which was likely induced by the Parra thrust fault.Also, fluvial incision may have been initiated by the drawdown of the base level, which affected both the Parra and Hondo Rivers.This stage most likely coincided with the cave passage development in the Naranjal River, suggesting that the diversion of the Parra River may have occurred before the collapse of the cave passage roof.At a final stage, the river diversion and incision led to the formation of the current meander embedded in a narrow valley of steep slopes surrounding the alluvial fan close to Parra Village.The current position of the meander is clearly delimited by two elbows depicted in figure 12-b.
The second feature analyzed is a sharp meander in the Hondo main channel, close to its confluence with the Parra River (see figure 12-b).The most suitable explanation for this feature -among several evaluated, which when tested were rejected -proposes a simple association with either a change in lithology or tectonic tilting.This assumption is consistent with both the change in magnitude of T and the geological set in the area of the meander, and it also fits well with the evolution of the Parra drainage network already described figure 10).

Proposed timeline and concluding remarks
Using the evidence available, a proposed timeline of the geomorphological evolution, highlighting the role of drainage rearrangement, is discussed in this section (see summary in table 5).Four stages, labeled A to D, are described as follows.
Conventionally, stage A started when the basins experienced rapid uplift.In the Naranjal Basin, this uplift was particularly intense in the eastern half, which led to the development of both steep slopes and incised streams.Furthermore, in the Hondo and Parra Basins, the uplift was particularly controlled by the Parra thrust fault and other fault lines, which favored the development of steep slopes upstream, as well as an alluvial fan at the confluence of these two rivers.
During stage B, the Naranjal Basin may have experienced a steady and gradual uplift, which led to the widening and aggradation of the central valley.At this stage, the development of both a ponor and a cave passage may have been initiated in El Manaclar Fm. limestone, near the confluence with the Parra River.On the other hand, the Hondo and Parra Basins may have experienced further uplift, along with drainage rearrangement processes downstream, of which the most remarkable was the stream diversion of the Parra River surrounding the alluvial fan in the Outlet Zone.
At stage C, the Naranjal River may have expanded the cave passage, followed by the collapse of the cave roof.Furthermore, both the Hondo and Parra Rivers may have incised further within their respective alluvial fans.In the Parra River, the fluvial incision along with the stream diversion may have also sharpened the meander surrounding the alluvial fan.
Finally, at stage D, the Naranjal River may have experienced little or no fluvial incision, since its newly sub-aerially exposed channel began to flow over hard limestone bedrock, and igneous gravel began to pass towards the Parra River.In addition, tufa concretions started to cover large blocks of the collapsed cave roof, which seems to have been occurring until the present time.Both the Parra and Hondo Rivers likely experienced little fluvial incision, with little or no valley widening.To conclude, this research paper presents evidence of drainage rearrangement processes in the geomorphological evolution of the Parra River network, focusing particularly on a karst piracy model as the most suitable explanation for the striking variations in both the composition of alluvial deposits and geomorphological traits.Karst piracy is an understudied process, and since the Dominican Republic hosts many excellent examples of karstified limestone outcrops, new research in this area will contribute to a better understanding of the interaction between rivers and karst systems.

Figure 1 :
Figure 1: (a) The Ocoa River Basin, showing its location in the Dominican Republic.(b) Map of the Ocoa Basin, showing its location on the southern hillside of the Central Range.The solid blue lines in (b) depict both main courses of the Banilejo River and the Ocoa River itself.(c) Detailed map of the three main sub-basins of the Parra River Basin, namely Naranjal, Hondo, and Parra, as well as the area conventionally referred to as the Outlet Zone.The background is a relief view (red-white is highland, green is lowland), based on a 30-m Shuttle Radar Topography Mission -digital elevation model (SRTM DEM) (NASA LP DAAC, 2000.SRTM 1 Arc-Second Global.https://earthexplorer.usgs.gov/.Published September 2014) .

Figure 3 :
Figure 3: Context view photos of the sample locations.(a) Terrace scarp close to the Parra River (sample S1), November 2015; (b) Stream bed of the Parra River during the dry season (sample S2), December 2014; (c) The Hondo River, flowing from right to left in the image (sample S3), October 2018; (d) Stream bed of the Naranjal River (sample S4), December 2018; (e) and (f) Terrace scarps in the Outlet Zone (samples S5 and S6), November 2015 and December 2018, respectively .

Horton ( 1945 )
Total length of streams L = Σl i Horton (1945) Drainage density, or average length of streams per unit of area D = L/A

*
Legend: R c = circularity ratio; K = compactness coefficient; T c = time of concentration in decimal hours and in time format (mm:ss); C Am = profile concavity (C A ) of the main channel; C As = average C A of the LFPs of the sub-basins draining into the main channel; I = hypsometric integral; T = mean vector magnitude of T and number of basin segments sampled (n, enclosed in parentheses); p = probability of obtaining a greater mean vector magnitude by pure chance combination of random vectors; sign.: significance of T at α = 0.05, where 'n.s.' stands for non-significant result and the symbol * stands for significant result.See table 1 for formulas.Furthermore, basin shape and symmetry parameters are summarized in table 4. Both the circularity ratio and the compactness coefficient show that all sub-basins are quite elongated (see figure1-c), the Hondo sub-basin being the least elongated of the ensemble and having the shortest time of concentration.The profile concavity indexes (C Am ) show that the main channels of the sub-basins are essentially concave, that of the Hondo sub-basin being the most concave of the ensemble.However, on average, the LFP profiles of the Naranjal Basin are far more concave than those of the other two basins.The morphometric patterns of the LFPs by basin are thoroughly described in the following section.

Figure 5 :
Figure 5: Elevation profiles (a) and plan view sketch (b) of both the Hondo Basin main channel (thick black line labeled 1) and the LFPs of the sub-basins draining into it (colored lines labeled 2, 3, . . ., n, where n is the number of LFPs).The color scheme represents the groups classified from a UPGMA cluster analysis; the scheme is consistent across all panels.Panel (c) shows the concavity, elevation drop, and stream length of the LFP for classified groups.See text for details.

Figure 6 :
Figure 6: Elevation profiles (a) and plan view sketch (b) of both the Naranjal Basin main channel and the LFPs draining into it.Panel (c) shows concavity, elevation drop, and stream length of LFPs for classified groups.Tthe numbered key and the symbology follow the same format as figure 5. See text for details.

Figure 7 :
Figure 7: (a) Partial view of the sinkhole (0.85 Ha in area and ∼ 40 m in depth) located at the outlet area of the Naranjal River.Three levels of detail are shown as follows: (b) Two limestone blocks of El Manaclar Formation, the first of c. 5 m in width (foreground, indicated by the black arrow), the second of 10 m in height (background).(c) Enlarged view of the 5 m limestone block (hammer for scale).Tufa deposits cover a significant part of the block's surface.(d) Hand-sized sample of the block, which shows the dark limestone of El Manaclar Formation.The yellow ellipses show a 1.60 m-high person for scale.

Figure 8 :
Figure 8: Elevation profiles (a) and plan view sketch (b) of both the Parra Basin main channel and LFPs of the sub-basins draining into it.Panel (c) shows concavity, elevation drop, and stream length of LFPs for classified groups.Both the numbered key and the symbology follow the same format as figure 5. See text for details.

Figure 10 :
Figure 10: Sample vectors of the transverse topographic symmetry ( T ) of the Hondo, Naranjal, and Parra main channels.Samples of T with magnitude 0.2 or greater are indicated by the red arrows (maximum = 0.6), whereas open circles and black dots indicate the position of 0.1 ≤ T < 0.2 and T < 0.1, respectively.The symbol * denotes significance of T at α = 0.05.Both normal and thrust fault lines intersecting the study area were extracted from the geological map at 1:50,000 scale (sheet 6171-IV), Servicio Geológico Nacional (n.d.), http://repo.sgn.gob.do/mg50/MG6171-IV LaCienaga.pdf

Figure 11 :
Figure 11: Polar plots of the transverse topographic symmetry vectors ( T ) of the Hondo, Naranjal, and Parra Basins.Angular units in degrees from north in clockwise direction.At the center, magnitude of T = 0; at margin, magnitude of T = 1.Average magnitude and bearing denoted by a red triangle.The symbol * denotes significance of T at α = 0.05

Figure 12 :
Figure 12: (a) Planform of the hypothetical paleogeography during the pre-rearrangement stage of the studied drainage network, depicting the main channels and relevant geomorphological features.(b) Planform of the present main channels, geomorphological features, and fault-lines intersecting the study area.See text for details.Legend:The grey solid and dashed lines in (a) represent the present main channels and watershed divides, respectively.The circle filled white in (a) near the Naranjal River outlet depicts the possible location of a ponor (labeled as "P").The grey shaded polygons overlaid by dashed lines depict alluvial fans.The dot pattern in the Naranjal Basin represents fluvial deposits.The solid blue lines represent, respectively, estimated and present main channel paths in (a) and in (b).The black semicircular hashed symbol in (b) depicts the location of the Naranjal River sinkhole (labeled as "S").The filled circle depicts the approximate location of Parra Village (labeled as "PV").The black curves with arrowhead symbol denote elbows of river diversion.The label "M" points to a sharply bent meander in the Hondo River.The labels "WG" point to wind gaps on the western divide of the Naranjal Basin.The solid grey lines in (b) depict fault-lines extracted from the geological map at 1:50,000 scale, sheet 6171-IV,Servicio Geológico Nacional (n.d.).

Table 1 :
Summary of the morphometric parameters computed from the DEM, formulas used and relevant references

Table 2 :
Observed frequencies of clasts by major groups of rock types in samples S1 to S4 = 96.89,p ≪ 0.05 for the entire sample; χ 2 = 8.58, p = 0.07 for samples S2 to S4.The p-values were computed by Monte Carlo simulation with 2000 replicates.The frequencies of undifferentiated clasts were not included in the tests.

Table 3 :
Areal, topographic, and drainage network parameters of the Parra, Hondo, and Naranjal Basins = area of the basin (in km 2 ); H min , H max , H, ∆H = minimum, maximum, mean, and drop values of the basin elevation, respectively (in m); S, S m = mean slope of the basin and main channel in %, respectively; N = number of streams; U = maximum Strahler stream order; M = number of first-order streams (basin magnitude); F 1 = number of first-order streams per km 2 ; R b = bifurcation ratio; L = total length of streams in km; D = drainage density in km/km 2 .See table 1 for formulas.

Table 4 :
Basin shape and symmetry parameters of the Parra, Hondo and Naranjal Basins

Table 5 :
Proposed timeline of the geomorphological evolution in relation to drainage rearrangement