Low variability , snowmelt runoff inhibits coupling of climate , tectonics , and 1 topography in the Greater Caucasus 2 3

This manuscript has been submitted for publication in Earth and Planetary Science Letters. This version has not undergone peer review and subsequent versions of this manuscript may have slightly different content. If accepted, the final version of this manuscript will be accessible via the “Peer-reviewed Publication DOI” link on the righthand side of this webpage. Please feel free to contact the authors, we welcome feedback.

Daily records of discharge, which we convert to runoff by dividing by drainage area (Fig.  204 2), for the Caucasus region comes from the Global Runoff Data Centre (GRDC) and 205 was also originally presented elsewhere (Forte et al., 2016). We reprocess the runoff 206 data here to remove basins whose variability may be artificially low due to dams and to 207 describe variability as a power-law fit of the right tail of the distribution, which we 208 describe in a later section. 209 Prior analysis of the GC runoff data found extremely low variability, which was 210 speculatively linked to the dominance of snowmelt runoff (Forte et al., 2016). To better 211 understand the cause of low daily runoff variability, we partitioned daily flows into 212 annual, seasonal, and event components in each basin (Table S1)  To develop a climatology of daily flows, we also calculate mean daily runoff as a 221 function of day of year and apply a 31-day moving mean to smooth over the influence of 222 historic events. Similar analyses on mean daily rainfall from TRMM are only used to determine the timing of peak rainfall in the main text, though full time series are shown 224 in Fig. S23. 225 226

Modern convergence rates 227
To compare the erosion rates to modern convergence rates, we follow prior 228 efforts which divided GPS velocities into either a Greater Caucasus or Lesser Caucasus 229 domain (Avdeev and Niemi, 2011;Forte et al., 2014) and calculated average velocities 230 along-strike using a sliding 50-km moving window (Fig. S1). Convergence between the 231 Lesser and Greater Caucasus is the difference between these velocities along-strike. 232 Our results are similar to prior estimates (Forte et al., 2014), but incorporate updated 233 GPS velocities (Sokhadze et al., 2018

Cosmogenic Erosion Rates from Alluvial 10 Be Inventories 250
Prior to field sampling, we vetted basins that appear to be in local topographic 251 steady-state (i.e., lacking major knickpoints; outside the influence of LGM glaciation) so 252 that basin-averaged 10 Be erosion rates could be reliably related to ksn (Fig. S3). This 253 analysis motivated the sampling of 76 basins across the southern range front of the Greater Caucasus. From these, a subset of 47 were processed for erosion rates (Table  255 S2). Low abundance of quartz and difficulty in processing some samples due to 256 lithology (see discussion in Supplement) resulted in usable amounts of quartz for 34 257 samples. For each sample, we selected the 0.25-1 mm size fraction and used a 258 combination of traditional HF and HNO3 leaches and the 'hot phosphoric acid' technique 259 (Mifsud et al., 2013) to isolate and purify quartz. Additional details for quartz purification 260 are described in the Supplement. Samples were spiked with either commercial or 261 custom low-background 9 Be carrier, Be was extracted through liquid chromatography, 262 and BeO was analyzed by accelerator mass spectrometry at PRIME Lab, Purdue 263 University. To convert 10 Be concentrations into erosion rates, we calculated effective 264 latitude and elevations to determine basin-averaged 10 Be production rates (Portenga 265 and Bierman, 2011), and used these in v3.0 of the online calculator formerly known as 266 the CRONUS calculator (Balco et al., 2008). Erosion rates are reported for a time 267 independent scaling scheme (Stone, 2000). Additional details with regards to site 268 selection, sample processing, and erosion rate calculations are provided in the 269 Supplement. All relevant parameters needed to reproduce erosion rates are provided in 270   Table S3. 271 Due to low quartz yields, we also carefully examined the bedrock geology for 272 each basin (Fig. S24-S57, Table S4) to assess the influence of variable quartz sourcing 273 on derived erosion rates. By recalculating topographic metrics and erosion rates after 274 removing portions of basins with lithologies unlikely to provide quartz, we found no 275 meaningful difference in the E-ksn patterns (Fig. S4, Table S3). We also considered the 276 end-member scenarios where we assume that quartz is entirely sourced from the upper 277 or lower 50% of each basin and recalculated topographic metrics and erosion rates (

Parameterization of SPIM 283
To assess which SPIM parameters best characterize the relationship between 284 channel steepness and 10 Be erosion rates, we fit eq. 3 to the measured E and ksn data.
To do this, we linearize eq. 3 using a log-transform and then fit the data using the 286 orthogonal distance regression (ODR) algorithm in SciPy. To estimate ranges of 287 acceptable fits, we tested both a bootstrap and Monte Carlo method. The latter is similar 288 to the method used by Adams et al. (2020). While results are comparable, the bootstrap 289 approach results in wider estimates of uncertainty. As such, we report uncertainties 290 using the bootstrap fits as more conservative estimates. Additional details of fitting are 291 laid out in the Supplement. In fitting the data, we exclude data from one basin whose 292 uncertainty exceeds its mean value (Fig. S11). We also test the sensitivity of fits to the 293 two highest erosion rates. While removal of these two rates suggest a lower n, the 294 range of uncertainties inclusive and exclusive of these data substantially overlap (Fig.  295   S11). Given the lack of any meaningful reason to exclude these high erosion rate data, 296 all reported fits include these high erosion rate basins.  Because none of the 10 Be basins are gauged, we generalize stochastic 308 parameters in gauged GRDC basins for attribution to ungauged 10 Be basins (Figs. 3-4). 309 To estimate ̅ in ungauged basins, we use the relationship between ̅ (Fig. 3) and 310 TRMM mean daily rainfall (MDP) in gauged basins. To do this, we first converted mean 311 daily discharge to mean daily runoff assuming a linear relationship between drainage 312 area and discharge (Fig. 2). Next, we used ground based precipitation stations (Klein 313 Tank et al., 2002) to bias-correct the TRMM 3B42 data and derive MDP for each 314 gauged basin (Fig. 3). The linear fit between MDP and runoff ratio ( ̅ /MDP) for gauged 315 basins is used to estimate ̅ in ungauged basins (Fig. 3). We use a linear fit to runoff ratio as opposed between MDP and ̅ , because the former is better approximated with 317 a linear fit than the latter. 318 Runoff variability is characterized using the shape parameter ( ) of the daily 319 distribution, which is estimated by fitting a power-law to the upper 1% of flows in each 320 gauged basin (Fig. 2). These results are generalized to ungauged basins using linear 321 regressions to the maximum elevation of the basin and the standard deviation of 322 monthly mean snow cover (Fig. 4, Table S3). We found these two metrics be the best 323 proxies for k after testing a variety of topographic and climatic metrics. Given that the 324 two metrics make slightly different predictions for individual basins and lack a clear 325 basis to choose one, we averaged the two estimates to derive values for the 326 ungauged basins. 327 The scaling between channel width and discharge ( ) is an important, and hard 328 to constrain, hydraulic geometry relationship that strongly controls the shape of the E-  Given the lack of direct constraints on either and the infeasibility of leaving both as free 345 parameters, one needs to be fixed to calibrate the model. Although we primarily report 346 solutions where is free and is fixed, we do test the alternative case (Fig. 5). While 347 values of or differ between optimizations, the E-ksn pattern is unchanged. Our goal 348 is find a single, best-fit value of that can be used as representative of the entire 349 erosion rate dataset. To arrive at this, we first treat each basin independently and 350 estimate . Following DiBiase and Whipple (2011) and fixing at 45 Pa, we used 351 STIM to find for each basin that most closely reproduces measured E for the known 352 value of ksn and estimated and ̅ (e.g., Fig. 3 Alternatively, we also independently estimate values for ̅ , , and by treating these 367 as free parameters and fit STIM to the measured ksn and E values using an ODR fit. The 368 associated algorithm for performing this fit is archived in the GitHub repository. systematically increases from the southern flanks of the range towards the core, 377 reaching a peak south of the topographic crest (Fig. 7B). Despite the wide range of E, 378 all data lie on a single, highly nonlinear relationship between ksn and E (Fig. 8A). Similar 379 relationships exist between E and mean basin slope due to the quasi-linear relationship 380 between ksn and slope in this setting (Fig. 8C, S9). Remarkably, over erosion rates from 381 ~300 to >5000 m Myr -1 , channel steepness is essentially invariant, ranging between 382 ~300-500 m (Fig. 8). While there is substantial scatter in these high E and ksn basins, 383 this is not unusual for these kinds of datasets. Moreover, detailed interrogation of 384 potential confounding factors reveals no meaningful way to subdivide these data into 385 different physically interpretable populations (Fig. S8). 386 387

River Incision Modeling 388
Fitting data with SPIM (eq. 3) suggests an n of 3.1 to 4 with a median value of 389 3.5 (Fig. 8, Fig. S11). This is in the range of n found elsewhere, but well above the 390 global mean value of ~2.5 (Harel et al., 2016;Lague, 2014). To see if direct 391 relationships exist between mean climate and either E and ksn, Figure 8B shows the 392 relationship between mean rainfall and erosion rate color-coded by channel steepness. 393 While E does not systematically vary with MDP ( Fig. 8B) or ̅ , both E and ksn do 394 increase where variability decreases (i.e., increasing ) (Fig. 8A). Given this outcome, 395 we turn to STIM which explicitly accounts for daily runoff variability, to see how well it 396 explains the strong nonlinearity in the empirical E-ksn relationship. 397 We use measured E and ksn and estimated ̅ and to estimate for each 398 sampled basin. Figure 6 shows these results and suggests that optimized varies over 399 six orders of magnitude and quasi-linearly varies with ̅ (Fig. 6). We do not think this values. We suspect channel width scaling to be a key source of inter-basin 407 heterogeneity. However, because we do not observe any clear relationship between channel width and either ̅ or E (Fig. S13), we do not think this important source of 409 uncertainty is systematically biasing STIM predictions. 410 To further explore variations in the parameterization of , we compare results to 411 an alternative optimization where is the free parameter and is fixed (Fig. 7C-D). 412 This exercise also produces an apparent relationship between the free parameter ( in 413 this case) and ̅ , albeit over a relatively narrower range of values. Whether optimizing 414 for or , the basic result is that STIM predictions for each basin retain a runoff 415 dependence that cannot be resolved with our data. Interestingly, this cross-correlation 416 does not appear for runoff variability, where no relationship emerges between and 417 optimized or ( Fig. 6A or C). Furthermore, neither optimized or systematically 418 vary with erosion rate or topography of the basins (Fig. S15). The wide range of 419 optimized or and their correlation with ̅ may reflect dynamics not included in 420 STIM like rules for sediment flux and bed cover (Sklar and Dietrich, 2006) (Fig. 8). The alternative approach of using an ODR fit to the measured ksn and E 427 to estimate ̅ , , and had very little impact on model results, independently 428 suggesting similar values for these three parameters as those found from the median 429 value approach (Fig. 8). In comparison to SPIM, both applications of STIM performs 430 similarly in goodness of fit metrics (Fig. 8, S16). In detail, the different models deviate 431 from measured values in different ways (e.g., SPIM shows better correspondence to 432 lower E and ksn data than STIM and vice versa; Fig. 8A). Despite comparable goodness 433 of fit, we favor STIM results because it enables a data-driven interpretation to the highly 434 nonlinear relationships observed. shortening rates applied to north-dipping structures with reasonable dips (Fig. 6), exists, likely due to local structural complexity, this result strongly contrasts with the 453 poor correlation between E and mean rainfall or estimated runoff (Fig. 8). As such, a 454 simple climatic control on E in this setting is unsupported and thus requires more careful 455 consideration of hydro-climatic controls on bedrock river incision itself, the focus of the 456 rest our discussion. 457 458

Strengths and Limitations of STIM 459
The ability of STIM to reproduce observed E-ksn relationships (Fig. 8, S16) 460 suggests that the shape of this relationship in the GC is aided by considering the local 461 hydro-climatology, namely the systematic decrease in runoff variability with elevation 462 (e.g., Fig. 4B). We relate these orographic patterns in variability, and thus the extreme Supplement, but this analysis suggests our data is best explained by three clusters 492 (e.g., Fig. 9A, S17). Cluster 1 has moderate ̅ (< 4 mm/day) with very low variability 493 ( >4). Cluster 2 has moderate ̅ with low variability (2> <4). Cluster 3 has high ̅ (>4 494 mm/day). Clusters were used both to evaluate model fits (Fig. 9) and to aid 495 interpretation of the underlying driver for differences in the runoff between basins (Fig.  496 10), which are discussed in turn below. 497 First, we consider model performance when clusters are considered separately. 498 Using the median values of ̅ and for each cluster but keeping fixed to the 499 population median, we model each cluster using STIM (Fig. 9D). This approach explains channel steepness patterns in low runoff basins (Clusters 1 and 2), but not in 501 high runoff basins (Cluster 3, Fig. 9D). Spatial patterns in E for Cluster 3 are consistent 502 with the tectonic forcing (e.g., Fig. S20) suggesting that ksn is anomalously high for at 503 least three of the fours basins in this cluster (i.e., those in Fig. 9D that lie substantially 504 above the modeled STIM relationship). Lithological differences do not explain 505 anomalously steep basins (e.g., Fig. S8) indicating that other model parameters must 506 differ for these basins and/or vary systematically vary with runoff, a finding supported by 507 where this cluster falls in the relationship between ̅ and (Fig. 9C). Nevertheless, 508 Clusters 1 and 2 represent the bulk of the data and show that subsampling the 509 population by differences in runoff variability is comparable to or improves upon STIM 510 predictions derived from the population as whole (Fig. 10, S21). 511 Next, we relate the clusters to hydro-climatological controls. Figure 10A shows warranted where rainfall also peaks in the summer. Cluster 1 basins show a strong 516 seasonal signal that is systematically offset from peak precipitation. All these basins 517 have very low baseflow in the fall and winter. Cluster 2 basins exhibit muted to non-518 existent seasonality. These basins show less systematic relations to the timing of peak 519 precipitation, though the overall phase lag is lower than in the other two clusters (Fig.  520 10A, S23). Cluster 3 basins all show strong seasonality and a systematic offset with the 521 timing of precipitation, like observed in Cluster 1. However, baseflow in the fall and 522 winter is typically very high compared to Cluster 1 and some lower elevation basins 523 show a second, lower peak in runoff in the winter. Regardless of cluster, higher 524 elevation basins typically show summer seasonality, reinforcing our interpretation that 525 snowmelt is the dominant driver of seasonal flows throughout the Caucasus region. 526 Figure 10A does not fully characterize the regularity of flows because data were 527 smoothed to develop a seasonal climatology. To probe whether and how well 528 streamflow seasonality explains the runoff variability parameter, k, we partitioned time 529 series data into three components: event, seasonal, and annual fractions which together 530 sum to the total water fluxing through each river. For gauged basins, the seasonal component was the strongest and only correlate to daily runoff variability, especially for 532 basins in the GC proper (Fig. 10B). Given that we attribute the seasonal component to 533 spring/summer snowmelt with modest contributions from seasonal rainfall in some 534 basins, we interpret patterns in runoff variability to be principally driven by the 535 contribution of snowmelt to runoff. Thus we interpret that our application of STIM is 536 accounting for orographic patterns in runoff variability that embed the long hydrologic 537 response times associated with snowmelt runoff (Deal et al., 2018). topographically sensitive to either climatic or tectonic changes. These areas: (1) have 558 higher runoff variability due to a lesser influence of snowmelt (Fig. 8-10), and (2) are in 559 the quasi-linear portion of the E-ksn relationship (Fig. 8). Conventional approaches 560 toward accounting for orographic precipitation in landscape evolution have focused on 561 elevation-dependent fluxes of mean annual rainfall (Bookhagen and Burbank, 2006) or 562 snowfall (Anders et al., 2008). This work highlights the critical role of the transition from 563 rainfall-to snowmelt-driven hydrology in mediating runoff variability itself (Rossi et al., 564 2020), an important complexity rarely considered in landscape evolution studies. 565 Transitioning from rainfall-to snowmelt-driven hydrology is dictated by the elevation 566 distribution within a mountain range and presents a possible direct relation between 567 climate and erosion rates in orogenic systems, albeit not in the traditional sense where 568 there is a positive correlation between erosion and precipitation or runoff rates (Ferrier 569 et al., 2013) . Importantly, a snowmelt control on runoff variability may be relevant to 570 many mountain ranges where the growth of topographic relief has undermined the 571 erosive ability of higher mean annual precipitation via distributing flows over longer 572 duration snowmelt events. 573 574

Conclusions 575
We present a large suite of new basin-averaged 10 Be erosion rates from the 576 Greater Caucasus that are consistent with longer term exhumation and shorter-term 577 decadal scale rates. Erosion systematically varies with convergence rates between the 578 Greater Caucasus and Lesser Caucasus and is uncorrelated to mean annual rainfall, 579 favoring a tectonic control on erosion rates. The relationship between erosion and 580 channel steepness is extremely nonlinear in this setting. However, careful consideration 581 of regional hydro-climatology incorporated into a stochastic threshold incision model of 582 river incision reveals that extremely low variability, snowmelt runoff is driving this 583 nonlinearity thus explaining why prior efforts failed to recognize a clear climatic imprint 584 on topography in the mountain range. 585 Our results also highlight the importance of both: (1) considering regionally 586 constrained relationships between topography and erosion when assessing potential 587 climate-tectonic interactions, and (2) understanding the underlying mechanism(s) 588 setting that form. In the Greater Caucasus, significant climate-tectonic interactions are 589 precluded because topography becomes insensitive to changes in forcing at uplift rates 590 exceeding 300-500 m Myr -1 . This is in contrast to other settings where relationships 591 between erosion and topography may be more linear. We emphasize that the observed 592 nonlinearity between erosion rates and channel steepness in the GC is not a global 593 solution to an apparent lack of coupling between climate and tectonics. Rather, the wide 594 range of such relationships around the world likely reflects important landscape specific, 595 hydro-climatic details that must be considered when applying erosion models. Our 596 results also show that spatial and temporal patterns in precipitation phase that alter 597 flood frequency may be an underappreciated governor on the degree of climate-tectonic 598 coupling possible in mid-latitude mountain ranges not heavily influenced by glacial 599 erosion. Relationship between optimized and estimated ̅ . The conversion between and 685 the threshold parameter (using the fixed of 45 Pa) is shown with the right-hand y-686 axis. (C) Same as 6A but colored by optimized . (D) Relationship between optimized 687 and estimated . Equation S19 is used to convert between and D50 and assumes 688 a shields parameter of 0.3. As in 6B, conversion between and the threshold 689 parameter (using the fixed of 2.24e-10) is shown with the right-hand y-axis. Cosmogenic 10 Be erosion rates vs distance from the center of the swath (colored by 698 mean elevation of sampled basins). Pearson's correlation coefficient (r) is shown 699 comparing erosion rates and distance from the swath center, along with respective p 700 value. (C) Cosmogenic 10 Be erosion rates vs distance along the swath (colored by 701 distance from the swath center). Grey regions indicate estimated vertical component of 702 uplift from the GC-LC convergence (Fig 1C) assuming convergence on a 2° or 45° 703 dipping thrust. Correlation coefficient between E and GC-LC convergence (Fig 1C) is 704 shown (see also between erosion rate and rainfall along with p-value is shown, note that this suggests 719 non-statistically significant correlation between these variables. (C) Mean basin ksn 720 compared to mean hillslope gradient, colored by E. Note that the linear relationship 721 between ksn and gradient reflects that both ksn and gradient become insensitive to 722 increases in erosion rate at ~500 m Myr -1 (Fig. S9).