Rock diversity drives river network reorganizations that boost continental biodiversity

Supporting Information for 1 2 Rock diversity drives river network reorganizations that boost 3 continental biodiversity 4 Pedro Val1*, Daniel Peifer2 5 1Queens College, City University of New York, NY, USA 6 2University of Tübingen, Tübingen, Germany 7 *Corresponding author: [email protected] 8 9 This Supporting Information contains the following information: 10 11 Supporting texts S1-S5 12 Text S1 – Geologic maps and Lithology-Topography relationships 13 Text S2 – Across-divide relief asymmetry histograms 14 Text S3 – Best-fit concavity for each study area 15 Text S4 - c - c ' comparison 16 Text S5 - Systematics and sensitivity of the landscape response to lithologic change 17 18 Supporting figures S1-S28 19 Fig. S1-S5 – Geologic and Topographic maps 20 Fig. S6-S9 – Topographic and river profile analysis of study areas 21 Fig. S10-S13 – D H distributions for each drainage divide 22 Fig. S14-S15 – Topographic analysis of numerical models 23 Fig. S16 – Numerical modeling results showing basin elevation and area change. 24 Fig. S17 – Numerical modeling results and the D HDc ’ plot. 25 Fig. S18 – River captures disturb the D HDc ’ plot. 26 Fig. S19-S22 – Best-fit concavity distributions for each study area 27 Fig. S23 – Pulses in sediment flux resulting from river captures due to lithologic change. 28 Fig. S24 – River capture sensitivity to background rock uplift rates. 29 Fig. S25-S28 – Comparison between D HDc ’ and D HDc relationships 30 31 Supporting Tables S1-S4 32 Table S1 – Divide asymmetry data for study areas 33 Table S2 – Model setups 34 Table S3 – Slope of the D HDc ’ relationship for different values of m/n 35 Table S4 - Regression coefficients of D HDc ’ and D HDc relationships 36 37 Supporting Video S1-S7 38 Video S1 – Numerical landscape evolution model revealing the topographic and river 39 network response to the exhumation of a resistant lithology in the first 50 My of model 40 run. Also shown are the responses of river profiles in a growing (blue) and shrinking 41 (pink) basin. 42 Video S2 – Numerical landscape evolution model revealing the evolution of drainage 43 divides and the temporal evolution of the DH-Dc ’ plot for model BLU4d, m/n = 0.45. 44 Video S3 – Numerical landscape evolution model revealing the evolution of drainage 45 divides and the temporal evolution of the DH-Dc ’ plot for model BLU4d, m/n = 0.3. 46 Video S4 – Experiments with variable outcrop length (BLU5a and BLU5b) revealing the 47 dependence of the positive area-loss feedback on outcrop length. 48 Video S5 – Erosion vs Uplift evolution of model BLU4c and sediment flux through time. 49 Video S6 – Same a Video S3, but for model BLU4c, m/n = 0.45, n = 2, and an order of 50 magnitude contrast in lithologic erodibility. 51 Video S7 – Same a Video S7, but two orders of magnitude contrast in lithologic 52 erodibility. 53 Supporting videos available at: https://doi.org/10.5281/zenodo.13328316 54 55 Supporting dataset available at for peer review: 56 https://zenodo.org/records/13328316?preview=1&token=eyJhbGciOiJIUzUxMiJ9.eyJpZCI6IjI157 ZjBmMmQ2LTNlZjUtNGRmNi1hYzcyLTgzYjdhYmNmYzhhZSIsImRhdGEiOnt9LCJyYW5kb20iO58 iJmZWQyMmE1NjgzYzRiYjExMGRkMzg2YzFjZTRkMjFmYyJ9.8oyYQGo2LfyRP5z8Zs2IqAZeJ59 WG3l55E-fnFpBoQPBJCohhYCXmdOe-hH2v2CGg8RDxZdBbTSJM_CgXQqfwcog 60 61 Text S1: Geologic maps and relationships with the topography of the study areas 62 In the study areas, the lithologic transition often marks a distinct escarpment 63 perpendicular to the main regional flow direction (Fig. S1): (A) between the Guiana shield and 64 the sedimentary rocks of the Amazon sedimentary basin (Fig. S2); (B) between the Serra Geral 65 basalts and silicified sedimentary sequences as part of the mesozoic-cenozoic megasequences 66 (escarpment and to due west; Menegazzo et al., 2016) and the non-resistant sandstones and 67 siltstones that form the paleozoic-mesozoic megasequences to the east (Almeida, 1949) (Fig. 68 S3); (C) between the Uinta formation and the less resistant rocks of the Green River Formation 69 (east) (Fig. S4); (D) between the Glen Canyon Group demarcating the edges of the San Rafael 70 Swell and the Moenkopi Formation within the swell (Fig. S5). Here, the Navajo Fm. is the most 71 resistant unit, but the Kayenta Fm. also forms an escarpment. In all of the study areas, note that 72 the observed divide asymmetries are not necessarily coincident with cross-divide lithologic 73 contrasts. 74 75 Text S2: Across-divide relief asymmetry histograms 76 Across-divide relief asymmetry is determined as the divide-to-channel head relief 77 (vertdist2stream.m from TopoToolbox; Schwanghart and Scherler 2014). For each main drainage 78 divide (Fig. S2-S9), we show the distributions of asymmetries as kernel densities (Figures S10-79 S13). From these data we determine the average asymmetry using a 95% bootstrap confidence 80 interval on the mean (Matlab function bootci.m). Note that the peaks of the distributions are non-81 zero, demonstrating that we accurately capture drainage divide asymmetries. 82 83 Text S3: Best-fit concavity 84 The histograms below show the distribution of best-fit concavity computed from hundreds of 85 drainage basins with 108 m2 drainage area. The concavity dictates the slope of the delta-plots 86 shown in Fig. 3 of the Main text (Table S3). 87 88 Text S4: Systematics and sensitivity of the landscape response to lithologic change 89 To guide our analysis and interpretation of the transient landscape features observed in 90 the study areas, we use numerical models of landscape evolution based on the stream power 91 incision law (see Methods). Once the hard rock is exhumed, the divide asymmetry ensues along 92 the entire drainage divide and initiates divide migration which, in turn, triggers in-sequence 93 transient geomorphic responses: (a) base level rise everywhere upstream of the exhumed hard 94 rock (Fig. 2) (Fig. S16); (b) larger drainage basins reach an equilibrium profile faster, creating 95 relief next to smaller basins (Fig. S16); (c) across-divide asymmetry ensues first in the areas near 96 the lithologic contact, but quickly propagate in the upstream direction (Supporting Video 1); (d) 97 divide migration and river captures then progress from the headwaters (far from the lithologic 98 contact) down to the lithologic contact; (d) over time, the differential relief grows due to an area-99 loss positive feedback (c.f. Yang et al, 2015) until the losing drainage basin is completely 100 consumed by the expanding basin. River captures lead to breaks in the slope of the chi-profiles 101 just as those observed in the real landscapes. Until all basins in the model equilibrate, divide 102 migration and river captures continue for up to hundreds of millions of years (Fig. 4 in main 103 text). Once basins reach an equilibrium between rock uplift and erosion at the lithologic 104 transition, drainage area exchange ceases (Fig. 4c, d in main text). This process forms low-relief 105 surfaces and plateau-top relief variability at various spatial scales such as in the Guiana Shield 106 (Fig. 1), especially if influenced by geodynamic or lithospheric reactivated uplift (Gernon et al., 107 2024; Ruetenik et al., 2024). 108 The lithologic influence on divide asymmetry produces continuous drainage area 109 exchange and divide migration. As the shrinking basin loses drainage area, the average elevation 110 increases while the expanding basins grows and decreases its average elevation (Figure S14). 111 This dynamic exchange creates relief between basins and can be detected in a delta-plot which 112 accounts for both the lithologic change as well as the drainage area loss (Figure S15). The delta-113 chi-prime proxy (Dc’) predicts the across-divide asymmetry (DH) while accounting for lithologic 114 disturbance and differences in the drainage area between basin pairs. The slope of the Dc’- DH 115 relationship depends on channel concavity (Table S3). Every river capture produces a transient 116 deviation from the expected relationship (Fig. S17, S18; Supporting Videos 2, 3) while large 117 river captures produce a Z-shape in the delta-plot (Fig. 3b). We make this observation based on 118 numerical models computed using a concavity index (m/n) of 0.3 which reveal larger river 119 captures due to the longer response times. We use these model observations to produce the 120 simplified conceptual plots in Fig. 3 of the Main Text. Importantly, the scatter in the Guiana 121 Shield data is consistent with these numerical models and supports the interpretation that the 122 scatter in the delta-plot (Main Text Fig. 3) is at least partly related to large river captures (Fadul 123 et al., 2022). 124 For a linear dependence of erosion on channel slope, the divide migration rates (therefore 125 the response time) are independent of rock uplift rate (Fig. 4c, d in main text), consistent with 126 theoretical predictions (Whipple et al., 2017). For models with nonlinear dependence on channel 127 slope (i.e. slope-exponent n > 1), higher erodibility contrasts are needed to drive the positive 128 feedback (Supp. Video 6, 7). As erodibility contrasts increase beyond an order of magnitude, 129 river basins may become internally drained and cause deviations from the relationships shown in 130 Figure 4 and 5 in the main text. These lithologic controls persist across a range of rock uplift 131 rates from 1 to 200 m/Ma (Fig. S24). 132 Finally, the river captures that ensue in response to the lithologic control of base level 133 produce transient pulses of erosion (Supp. Video S5). When converted to an erosive flux and for 134 a background uplift rate of 20 m/My, these erosive pulses produce sedimentary fluxes on the 135 order of 1-5 x 103 km3/My. These values are consistent with the magnitudes and temporal 136 patterns of sediment fluctuations observed in passive margins (Pazzaglia and Brandon, 1996; 137 Contreras et al., 2010). 138 139 Text S5: c - c ' comparison 140 We check our chi-prime metric against the regular chi measurement (c.f. Perron and Royden, 141 2013). Plots (Fig. S25-S28) and regression data (Table S2) show that the regular chi data 142 captures sthe same trends against the across-divide relief. However, not accounting for the 143 lithologic differences reduces the confidence on the relationship. 144 145 146 Figure S1: Location of select transient topography upstream of hard rocks in continent interiors. a) Eastern Guiana 147 Shield (EGS) in northern Pará, BR; b) Oblique view of the plateau-top landscape transients and shrinking basin; c) 148 Eastern Paraná Basin (EPB) in Sao Paulo, BR (Depressão Periférica Paulista); d) San Rafael Swell (SRS) and 149 Tavaputs Plateau, Utah, USA. The chosen areas are characterized by low lithologic complexity (EPB, TP), large 150 river captures (EGS), and high lithologic complexity (SRS) (Fig. 1). Topographic data obtained from Copernicus 151 Digital Elevation Models with 30m and 90m resolution. Location of hard rock units obtained from the Brazilian and 152 Utah geological surveys. 153 154 155 Figure S2: Geology and topography of the Eastern Guiana Shield. Top: 1:1,000,000 geologic map of northern 156 Brazil (Bahia et al., 2004; Faria et al., 2004a, 2004b; Faraco et al., 2004). Vector data provided in Supporting 157 Dataset. Bottom: Digital elevation model (COP 90 m) showing topographic escarpments and coincidence with 158 resistant lithologic units of the Trombetas Group (sandstones). 159 160 161 Figure S3: Geology and topography of the Depressão Periférica Paulista, Eastern Paraná Basin. Top: 1:1,000,000 162 geologic map of southeast Brazil (Heineck et al., 2003; Leite et al., 2004; Lopes et al., 2004; Valente et al., 2004). 163 Bottom: Digital elevation model (COP 90 m) showing topographic escarpments and coincidence with resistant 164 lithologic units of the Serra Geral Fm. (basalts) and Botucatu Fm., Marília Fm. (silicified sandstones and 165 conglomerates). The legend is simplified to the sedimentary rock units where the divide asymmetry measurements 166 were made. Basement rocks are largely not involved in the analysis. 167 168 169 Figure S4: Geology and topography of the southern limits of the Tavaputs Plateau. Top: Geological map of the 170 eastern Tavaputs Plateau (Sprinkel, 2009); Bottom: Digital elevation model (COP 90 m) showing topographic 171 escarpments and coincidence with resistant lithologic units of the Uinta Formation. 172 223 224 Figure S12: Delta-H for the Tavaputs Plateau divides. 225 226 227 228 Figure S13: Delta-H for the San Rafael Swell divides. 229 230 Figure SD14: Output from one numerical landscape evolution model reproducing the empirical data. This model 231 was computed with a background uplift rate of 40 m/Ma and outcrop length of 10 km (BLU5b model, see Supp. 232 Table S2). a) Topographic output with basin outlines (black), rivers (white), and selected rivers (colored) for c-233 profile analysis. Also shown are traces of swath profiles; b) Swath profiles highlighting perched basins with 234 asymmetric drainage divides on their flanks mimicking the real landscapes (Fig. S5-9); d) c plot showcasing a river 235 capture of the shrinking basin’s (pink) headwaters. 236 237 238 Figure SD15: Output from one numerical landscape evolution model. This model was computed with a background 239 uplift rate of 20 m/Ma, outcrop length of 20 km, with slope-exponent n = 2 (m/n = 0.45), and a hard rock with K = 5 240 x 10-8 m0.1/y (BLU4b_n2_kc4 model, see Supp. Table S2). a) Model snapshot after 9.8 My; b) Model snapshot at 9.9 241 My; c) Model snapshot at 10 My; d-e) Zoom-in of river captures and the formation of wind-gaps (c.f. Fig. S5-S9). 242 Throughout these 300 ky, the model highlights the formation of systematically shrinking basins (a), river captures, 243 drainage reversals, and wind-gaps (b-c) reproducing the geomorphic signals of the real landscapes. 244 245 246 Figure S16: Numerical modeling results for model BLU4c revealing the evolution of the drainage area and basin 247 average elevation for a pair of basins undergoing drainage area exchange (same pair shown in Fig. 2 of the main 248 text). The growing difference in average elevation is consistent with the growing divide asymmetry (i.e. across 249 divide relief) shown in Fig. S17. 250 251 252 Figure S17: Numerical modeling results for model BLU5d revealing the evolution of the delta plot. The graph 253 shows the progression of the DH-Dc’ trend over model time. This plot uses a background uplift rate of 40 m/Ma, 254 rock outcrop length of 30 km, and concavity index of 0.45. Supporting Videos 2-3, 6-7 show the evolution of the plot 255 over time. 256 257 258 Figure SD18: Numerical model example of large river captures with n = 1, m/n = 0.3. a-c) Snapshots of model 259 outputs highlighting shrinking basins and migrating drainage divides; d-f) Delta-plot showing how river captures 260 disrupt the DH - Dc’ relationship in 280 My and 320 My (c.f. Fig. 3a). 261 262 263 264 Figure S19: Best-fit concavities for basins in the Guiana Shield study area (n = 1513). The average for this plot is 265 0.48. 266 267 268 269 Figure S20: Best-fit concavities for the eastern Paraná basin study area (n = 262). The average for this plot is 0.5. 270 271 272 273 274 275 Figure S21: Best-fit concavities for the Tavaputs Plateau. Basins in this range from 107 m2 to 109 m2, therefore, we 276 computed best-fit concavities for basins with a minimum 107 m2 drainage area (n = 81) and capturing river channels 277 up to a headwater area of 1 x 104 m2. The average best-fit m/n is 0.35. 278 279 280 Figure S22: Best-fit concavities for the San Rafael Swell study area. Given the small size of some drainage basins 281 (i.e. 106 m2), concavities were computed for basins with a 5 x 106 m2 and capturing river channels up to a headwater 282 area of 1 x 104 m2 (n = 431). The average m/n is 0.26. 283 284 285 Figure S23: Sediment fluxes associated with model BLU4c (See Table S2; Supp. Video S5). Large fluctuations in 286 sedimentary flux are associated with river captures. Note that pulses of sediment flux ensue 10s to 100s of millions 287 of years after the exhumation of the hard rock (initiated in the first 100 ky of model run). In a real landscape, these 288 erosive fluxes are likely accompanied by increases in nutrient availability due to erosion of the soil and rock 289 columns. 290 Heineck, C. A.; Vieira, V. S.; Drummond J. B. V.; Leite C.A.L.; Lacerda Filho, J. V.; Valente, 364 C. R.; Souza, J. D.; Lopes, R. C.; Malouf, R. F.; Oliveira, I. W. B.; Oliveira, C. C.; Sachs, 365 L. L. B.; Paes, V. J. C.; Junqueira, P. A.; 2003. Folha SE.23 Belo Horizonte. In: 366 Schobbenhaus, C.; Gonçalves, J. H.; Santos, J. O. S.; Abram, M. B.; Leão Neto, R.; 367 Matos, G.M.M; Vidotti, R. M.; Ramos, M. A. B.; Jesus J. D. A. de; (eds). Carta 368 Geológica do Brasil ao milionésimo, Sistema de informações geográficas, Programa 369 Geologia do brasil. CPRM, Brasília. 370 Leite, C. A. S; Perrotta, M. M.; Silva, L. C.; Silva M. A.; Heineck C. A.; Salvador, E. C.; Vieira 371 V. S.; Lopes, R. C.; Silva, M. G. M.; Drumond, J. B. V.; Malouf, R. F.; Lacerda Filho, J. 372 V.; Valente, C. R.; Gomes, S. D.; Sachs, L. L. B.; Oliveira, I. W. B.; Ramgab, G. E.; 373 Netto, C.; Junqueira, P. A.; Paes, J. V. C., 2004. Folha SF 23 Rio de Janeiro. In: 374 Schobbenhaus, C.; Gonçalves, J. H.; Santos, J. O. S.; Abram, M. B.; Leão Neto, R.; 375 Matos, G.M.M; Vidotti, R. M.; Ramos, M. A. B.; Jesus J. D. A. de; (eds). Carta 376 Geológica do Brasil ao milionésimo, Sistema de informações geográficas, Programa 377 Geologia do brasil. CPRM, Brasília. 378 Lopes, R. C.; Peruffo, N.; Sachs, L. L. B.; Silva V. A.; Batista I. H.; 2004. Folha SF 22 379 Paranapanema. In: Schobbenhaus, C.; Gonçalves, J. H.; Santos, J. O. S.; Abram, M. B.; 380 Matos, G.M.M; Vidotti, R. M.; Ramos, M. A. B.; Jesus J. D. A. de; (eds). Carta 381 Geológica do Brasil ao milionésimo, Sistema de informações geográficas, Programa 382 Geologia do brasil. CPRM, Brasília. 383 Menegazzo, M. C., Catuneanu, O. & Chang, H. K. The South American retroarc foreland 384 system: The development of the Bauru Basin in the back-bulge province. Marine and 385 Petroleum Geology 73, 131–156 (2016). 386 Pazzaglia, F. & Brandon, M. Macrogeomorphic evolution of the post-Triassic Appalachian 387 mountains determined by deconvolution of the offshore basin sedimentary record. Basin 388 Research 8, 255–278 (1996). 389 Perron, J. T., & Royden, L. (2013). An integral approach to bedrock river profile analysis. Earth 390 surface processes and landforms, 38(6), 570-576. 391 Ruetenik, G. A., Jansen, J. D. & Sandiford, M. Surface Uplift Due To Time-Varying Elastic 392 Thickness in Continental Interiors. Geophysical Research Letters (2024). 393 Schwanghart, W. & Scherler, D. Short Communication: TopoToolbox 2 – MATLAB-based 394 software for topographic analysis and modeling in Earth surface sciences. Earth Surface 395 Dynamics 2, 1–7 (2014). 396 Sprinkel, D. A; 2009 Interim Geologic Map of the Seep Ridge 30' x 60' Quadrangle, Uintah, 397 Duchesne, And Carbon Counties, Utah, And Rio Blanco and Garfield Counties, 398 Colorado. Utah Geological Survey Map 549DM. 399 Valente, C. R.; Lacerda Filho, J. F.; Rizzotto, G. J.; Lopes, R. C.; Romanini, S. J.; Oliveira, I. W. 400 B.; Sachs L. L. B.; Silva, V. A.; Batista I. H.; 2004. Folha SE 22 Goiâ. In: Schobbenhaus, 401 C.; Gonçalves, J. H.; Santos, J. O. S.; Abram, M. B.; Matos, G.M.M; Vidotti, R. M.; 402 Ramos, M. A. B.; Jesus J. D. A. de; (eds). Carta Geológica do Brasil ao milionésimo, 403 Sistema de informações geográficas, Programa Geologia do brasil. CPRM, Brasília. 404 Whipple, K. X., Forte, A. M., DiBiase, R. A., Gasparini, N. M. & Ouimet, W. B. Timescales of 405 landscape response to divide migration and drainage capture: Implications for the role of 406 divide mobility in landscape evolution. Journal of Geophysical Research: Earth Surface 407 122, 248–273 (2017). 408 Yang, R., Willett, S. D. & Goren, L. In situ low-relief landscape formation as a result of river 409 network disruption. Nature 520, 526–529 (2015). 410