Hydrologic Connectivity and Patch‐To‐Hillslope Scale Relations in Dryland Ecosystems

In drylands, runoff during storms redistributes water and nutrients from bare soil areas to vegetated patches, subsidizing vegetation with additional resources. The extent of this redistribution depends on the interplay between surface roughness and permeability; greater permeability in vegetated patches promotes run‐on to vegetation, but greater surface roughness diverts runoff, producing tortuous flow paths that bypass vegetation. Here, this interplay is examined in virtual experiments using the 2D Saint Venant Equations to measure runoff connectivity. Flowpaths are delineated using tracers advected by the flow. Distances between tracer sources and sinks along flowpaths measure hydrologic connectivity at two lengthscales: connectivity to the hillslope outlet and within‐slope source‐sink connectivity. Differences between these connectivity lengthscales indicate how flow may “by‐pass” vegetated patches within hillslopes. At the hillslope scale, a derived power‐law relation between the runoff coefficient and outlet connectivity describes hillslope water losses, providing a foundation for identifying landscapes likely to shed water.

Connectivity provides a framework for characterizing water fluxes within hillslopes and across scales (Bracken & Croke, 2007;Cammeraat, 2004;Jencso et al., 2009;Turnbull & Wainwright, 2019). The connectivity of runoff-and the associated transfer of resources from patch to hillslope scales-provides an indicator of ecohydrologic function Turnbull & Wainwright, 2019). Changes in hydrologic connectivity can lead to abrupt shifts in ecosystem productivity by shifting the role of runoff from redistributing resources to vegetation, to exporting resources from the landscape (Bracken & Croke, 2007;Cammeraat, 2004;Moreno-de las Heras et al., 2012;Puigdefábregas, 2005;Sidle et al., 2017;Wainwright et al., 2011). For example, under degraded conditions with limited vegetation cover, increasing connectivity of bare soil areas may shift the function of runoff from resource-enhancing to resource-shedding (Okin et al., 2015). Increasing "run-around" behavior is expected to increase the connectivity of runoff, whereas increasing run-on behavior would decrease this connectivity.
This change in the function of runoff should be reflected in the metrics describing both within-hillslope connectivity and the connectivity of transport from the hillslope. If flow is substantially diverted around vegetation, the lengthscales connecting runoff sources to the outlet would exceed those connecting runoff sources to infiltration sinks (vegetation). Therefore, connectivity lengthscales would be expected to reflect how much flow bypasses or infiltrates into vegetated patches, linking the constitutive properties of the hillslope (e.g., permeability, roughness, spatial pattern) to hillslope-scale hydrologic outcomes such as the runoff coefficient, C, defined as the fraction of incident rainfall that leaves the hillslope domain as runoff. Thus, connectivity could act as a conceptual "bridge" between local properties of sites within a hillslope, and the hillslope's bulk hydrological behavior (i.e., C).
The measurements needed to define relations between roughness/permeability, connectivity and runoff coefficients, however, are not generally available. Efforts to quantify runoff connectivity on hillslopes-and its relation to landscape structure, vegetation and rainfall-have focused on quantifying runoff connectivity to the outlet (Antoine et al., 2009), within-slope connectivity (Marchamalo et al., 2016;Turnbull & Wainwright, 2019), or water-balance partitioning at the hillslope scale (Crompton et al., 2019). Generally, not all three components  Dunkerley (2018). Panels b and c show the maximum flow depth h for two rainfall intensities: p = 5 cm/hr (panel b) and p = 1 cm/hr (panel c), with t R = 60 min, F v = 0.3, σ x = 4, K v = 3 cm/hr, and K b = 1 cm/hr. Cumulative infiltration is shown projected onto the horizontal plane. In the case with p > K v (b), vegetative resistance contributes to increased ponding in the vegetation; conversely, with p < K v (c), greater permeability reduces the depth in the vegetation relative to the surrounding bare soil. Panel f shows the depth profiles for a transect in Panels b and c, with gray lines showing the 1-dimensional analytic approximation for steady-state flow (derived by combining conservation of mass conservation and force-balance, S fx = S x ). Panels d and e illustrate the paths of tracers that infiltrate after falling on the hillslope and tracers that escape, respectively, for the case in panel b. The average lengths of the tracer paths define L infl and L esc , with L infl = 52 m and L esc = 84 m for the case shown here.
3 of 11 are simultaneously reported, or reported in conjunction with varying vegetation patch permeability and roughness. Although spatial insight into overland flow is often difficult to obtain through field experiments (but see Dunkerley, 2004), virtual experiments offer a means to predict flow characteristics across a wide variety of conditions that can be controlled (Fatichi et al., 2016;Turnbull & Wainwright, 2019). Passive tracers in virtual experiments reveal flow paths from patch-to hillslope-scales (H. G. Smith et al., 2015) and therefore allow the characterization of source-sink and source-outlet connectivity (Keesstra et al., 2018;Mügler et al., 2011) along with their co-variation with rainfall properties, vegetation roughness, and vegetation pattern.
Here, runoff on patchily vegetated dryland hillslopes is predicted for discrete rainfall events using FullSWOF_2D (Delestre et al., 2014), an open-source two-dimensional Saint Venant Equation (SVE) model that has been previously validated on analytic solutions to the SVE and real rainfall events (Delestre & Esteves, 2010;Delestre et al., 2013). The SVE are used to link overland flow and runoff for a large range of hillslope vegetation and rainfall characteristics, and runoff connectivity is described using massless tracers advected with the local flow.
The simulation results are used to address two research questions: 1. How do source-sink and outlet connectivity lengthscales relate to each other and to the runoff coefficient? 2. Is a tendency for flow to be diverted around vegetation reflected in connectivity lengthscales?

Saint Venant Equations (SVE)
Sheet flow during rainfall events typically results in shallow flows (<5 cm), which can locally vary in depth and velocity, as water passes over surfaces of variable permeability and obstructions. Such shallow, unsteady and heterogeneous flows are conventionally represented by the SVE or "shallow water equations" (Brutsaert, 2005). The SVE combine the depth-averaged continuity and conservation of momentum equations and have been used to reliably predict runoff in laboratory and field settings (Cea et al., 2014;Kirstetter et al., 2016;Mügler et al., 2011;Taccone et al., 2020). The two-dimensional form is given as: where U, V denote the depth-averaged flow velocities in the x and y-directions, respectively; h is the depth of flow; S x and S y are the bed-slope in the x and y directions; t is time; and g represents gravitational acceleration. The terms S fx and S fy are the x-and y-components of the friction slope, which represent the effect of bed and other shear stresses on the flow, and are typically represented by empirical friction formulas (M. W. Smith et al., 2007). Rainfall inputs (p) and infiltration losses (i) add or remove mass from the system, but are assumed to have a negligible impact on local momentum transfers.
The imbalance between rainfall p and infiltration i represents a runoff-generation term that can vary with position and time. Rainfall is represented by a spatially uniform and temporally-invariant intensity p for a storm duration t R . While the modeling framework can accommodate variability of the infiltration capacity with antecedent wetness, we simplify the problem by assuming that infiltration capacity is dominated by gravitational drainage. This simplification implies that sorptive effects dominating at early times in dry soils are not accounted for, and thus the infiltration capacity may be underestimated.
The SVE do not form a closed system of equations, so a "closure" model in the form of a resistance formulation must be specified to represent the net effects of bed and other shear stresses (e.g., presence of obstructions) on S f . This closure, also written in one-dimensional form for simplicity, is given by: where | | = √ 2 + 2 and R h is determined by the vegetation and soil properties and consequently differs between vegetated and bare soil areas. To accommodate both surface cover types in a single formulation, the approach by James et al. (2004) is employed, summarized as: where f is the Darcy-Weisbach friction factor for the surface, D is the stem diameter, ϕ v is the stem density (the volume fraction occupied by the vegetation), = 4 ( 2 ) −1 is the number of stems per unit area, and C d is a vegetation drag coefficient, which in turn depends on the vegetation characteristics and the flow velocity |U|. To represent dryland grasses or shrubs, a "cylinder array" formulation for the drag coefficient was employed (Cheng & Nguyen, 2011), described in Text S1 in the Supporting Information S1.
Anticipating low bulk Reynolds numbers (Re = Uh/ν < 2,000) for the simulated p and hillslope dimensions (see Section 2.2), a modified laminar formulation for f was selected, described in Text S1 in the Supporting Information S1. To summarize, bare soil areas are characterized by ϕ v = 0, so that R h = f/(8gh). In vegetation, both stem drag and ground friction can be significant, with stem drag dominating for deeper flows and ground friction dominating for shallower conditions. In all flow simulations, the flow is assumed to be subcritical at all times and locations so that the Froude number Fr = U(gh) −1/2 remains below unity. This was subsequently confirmed in the simulation results.
The choice of resistance formulation is grounded in empirical runoff studies and extensive work characterizing resistance to flow through vegetation (see Tanino & Nepf, 2008, and references therein). While the prediction of resistance to overland flow through vegetation continues to be a challenge (M. W. Smith et al., 2007;Melis et al., 2019), the formulation described above captures the essential feature of spatial resistance contrast between bare soil and vegetation. The results are not anticipated to depend on the choice of resistance formulation (see e.g., the studies of Mügler et al. (2011), Cea et al. (2014, and Crompton et al. (2020), demonstrating that multiple resistance formulations can produce acceptable agreement with runoff measurements through calibration).

SVE Simulations
The SVE model was used to simulate storm runoff across a broad range of storm and landscape conditions, which are summarized in Table 1. The model domain is a planar hillslope with slope S o , length W x = 200 m, and width W y = 50 m, where x increases in the downslope direction from x = 0 at the divide.
The hillslope is assumed to be a binary mosaic of vegetated and bare soil areas, with the presence/absence of vegetation used to prescribe infiltration and surface roughness. This binary treatment of the hillslope is similar to previous modeling efforts, field, and remote sensing studies (e.g., Imeson & Prinsen, 2004;Moreno-de-las-Heras et al., 2020;Rodríguez et al., 2018). An open (Neumann) boundary condition is imposed at the base of the hillslope, and no-flux boundaries are imposed on the sides and top of the hillslope, ensuring that all incident rainfall is lost either by infiltration or runoff at the outlet.
The saturated hydraulic conductivity in vegetated patches K v was selected based on rates reported in field studies on runoff run-on processes in dryland regions (e.g., Cerdà, 1997;Dunkerley, 2002;Vásquez-Méndez et al., 2010). The bare sites ranged from nearly impermeable, with K B = 0.1 cm/hr (representing a limiting case), to K B = 1 cm/ hr. Roughness characteristics are similarly determined by whether the surface is bare or vegetated. To include a range of vegetation types, two values of the stem density ϕ v were included to represent, for example, sparse versus dense grass cover.
To generate many different spatial patterns of vegetation cover, binary patterns were created with varying lengthscales and degrees of anisotropy (Crompton et al., 2019). Patterns were generated from arrays of random numbers between 0 and 1, which were binarized to obtain a given vegetation cover fraction, F v , by adjusting the threshold value used to binarize the array ( Figure S1 in the Supporting Information S1). The characteristic patch lengthscale was adjusted by applying a two-dimensional Gaussian filter prior to binarizing. This filter uses a kernel defined by standard deviations in the x and y directions (σ x along-slope, σ y across-slope), which then set the lengthscale of vegetation patches in each direction. Because dryland vegetation often forms spatially-organized patterns (Noy-Meir, 1979), we adjusted the ratio of σ y /σ x to generate isotropic patterns (σ y = σ x ) and anisotropic patterns that extend along the hillslope contour (σ y /σ x > 1). Figure S2 in the Supporting Information S1 shows sample patterns with varying vegetation fractions, patch lengthscales, and degrees of anisotropy incorporated into the SVE.
In total, 1944 simulations were run. To illustrate the results, Figure 1 shows the maximum flow depth for two simulation cases (panels b and c), and the corresponding depth profile for a sample transect (panel f), showing reasonable agreement with the analytic solution for 1-dimensional, uniform steady-state flow. Additionally, animations are included in the supporting information (Movies S1-S3), with corresponding static visuals in Figure S4 in the Supporting Information S1.

Particle Tracers and Validation
Tracers, representing "tagged" water parcels, were used to describe runoff flowpaths. For each SVE simulation, 1,000 tracers were initialized at random positions on the hillslope and random times ranging from 0 to t R (the duration of the rain). The path of each tracer was determined assuming advection by the computed 2-dimensional velocity vector. At each timestep (dt = 60 s), a given tracer infiltrated with probability i · dt/(h′ + p · dt) where i is the local infiltration rate (K v or K b ) and h′ is the flow depth at the tracer's location.
The tracer computations were evaluated by comparing the tracer runoff coefficient C L -defined as the fraction of tracers that "escape" the domain rather than infiltrating-to the SVE-predicted runoff coefficient C. A small value of the difference, |C L − C|, indicates agreement between the SVE model and tracer computations (i.e., the escape fraction matches the runoff coefficient predicted by mass balance). Releasing 1,000 tracers per simulation was sufficient to achieve errors below 5% (see Figure S3 in the Supporting Information S1, which shows the distribution of errors in C L across the ensemble of simulations). For each simulation, the tracer trajectories were then used to quantify the runoff connectivity, as described next.

Measures of Connectivity and Hydrologic Outcomes
To characterize the connectivity of runoff within and out of the hillslope domain, the trajectories of the tracers were summarized using two indices: the source-sink connectivity L infl , and the outlet connectivity L esc . For a given simulation, L infl is defined as the mean downslope (x-direction) displacement of the subset of tracers that infiltrate within the domain. L esc is similarly defined as the mean displacement of tracers that escape the domain as runoff. While these definitions do not consider tortuosity in the flowpaths, the downslope displacements of the tracers are a reasonable approximation for the path lengths because V/U ≪ 1 (V is the across-slope velocity).
Both L infl and L esc are normalized by the hillslope length (W x = 200 m) so that they range from 0 to 0.5, and their normalized values are labeled with a tilde: ̃ and ̃ . The upper limit reflects the fact that, if every tracer path spans the hillslope length, then the mean distance traveled will follow a uniform distribution with mean 0.5, because the tracers were initially randomly distributed in the domain.
The definitions of ̃ and ̃ are illustrated in Figure 1, which shows the infiltrating and escaping paths (Panels b and c, respectively), for the simulation case with rain intensity p = 5 cm/hr, duration t R = 60 min, vegetation fraction F v = 0.3, and stem density ϕ v = 0.15. The yellow paths in Panel b indicate the subset of tracers that infiltrate, and blue paths in Panel c indicate tracers that escape.
As a measure of distances traveled by flow prior to leaving the domain, ̃ also provides an estimate for the portion of the hillslope contributing runoff to the outlet (i.e., the contributing source length). In analogy to catchment-and watershed-scale studies relating contributing source length to streamflow (Costa et al., 2020;Spence & Mengistu, 2019), we fit a power law to all of the simulations, to describe the relation between the runoff coefficient C and ̃ : where b and k are constant coefficients. In power law relations between contributing area and stream flow, the coefficients are known to vary as a function of the return interval and location (Costa et al., 2020). In this context, it is unknown whether the relation between contributing area and runoff coefficient is governed by a power law relation, and whether the relation varies across storm characteristics. Here, the suitability of a power law to explain the relation between C and ̃ is determined from the simulations.
As noted above, differences between within-slope and outlet connectivities may reflect whether flow is preferentially routed around the vegetation. To measure this effect, a tortuosity metric η was calculated as the mean standard deviation of each tracer's y position over the storm duration and subsequent drainage period: where N = 1,000 is the number of tracers, and std(y i ) is the standard deviation of tracer i's y-positions (perpendicular to the hillslope gradient). This definition reflects the extent to which the flow deviates from a 1-dimensional path parallel to the hillslope gradient. Larger values of η indicate that flow is diverted around vegetated patches.
The flow tortuosity is expected to increase with increasing resistance contrast between vegetation and bare soil, which in turn depends on the surface roughness (e.g., the within-patch stem density ϕ v ) and the runoff intensity. This intensity is described by the bulk Reynolds number Re = Uh/ν-computed as the domain and time-averaged value over the full simulation period. As shown in Text S1 in the Supporting Information S1, the sensitivity of the friction slope S f to ϕ v scales as

Results
The relation between outlet connectivity ̃ and the runoff coefficient C is stationary across rainfall and vegetation characteristics (see Figure 2a, where each marker corresponds to one of the SVE simulations). Upon fitting a power law, the coefficients in Equation 6 are estimated as k = 1.62 ± 0.02 and b = 2.93 ± 0.05 for all storm intensities and durations. It is relevant that k > 1 because increases in ̃ result in even larger increases in C and concomitant water loss from the hillslope.
The relation between source-sink connectivity ̃ and the runoff coefficient C exhibits more variation, in particular, greater sensitivity to the rainfall duration (marker colors in Panel b). ̃ is smaller than ̃ in most simulation cases (points fall below the one-to-one line in Panel c), indicating that the distance traveled by tracers that infiltrate is less than that traveled by those that escape. The tracer paths in Figure 1 suggest that this difference may be attributable to flow diversion around vegetation patches. The paths of escaping tracers display more tortuosity within the domain (Figure 1c), whereas the infiltrating tracers show more linear paths (Figure 1b). Thus, Figure 1 suggests that ̃ includes the tracers that are preferentially diverted around the vegetation, whereas ̃ is biased toward tracers that flow into the vegetation patches and likely infiltrate, thereby shortening their path lengths. To characterize this difference, Figure 3a shows ̃− as a function of log(η · t R ), indicating that the difference between ̃ and ̃ increases with increasing η and t R . The tortuosity increases with increasing vegetation density and Reynolds number (see Figure 3b, showing η as a function of Re   2∕3 ).

Discussion
The study illustrates an approach to relate connectivity lengthscales to fluxes observed at domain boundaries (in this case the hillslope) (Sidle, 2021)-an ongoing challenge for hydrologic modeling-illustrating how such relations can be derived from the physics of the system (or a best approximation of it). Across the range of conditions explored, the relation between outlet connectivity ̃ and hillslope runoff coefficient C (Equation 6) is stationary, illustrating how an emergent power law relation describing the transfer of information from patch to hillslope scales arises from a connectivity framework (Gomi et al., 2008;Keesstra et al., 2018).
Interestingly, the relation between ̃ and C is not stationary and depends on the storm duration t R , with the implication that the relation between ̃ and ̃ depends on t R as well. That the difference between ̃ and ̃ increases with increasing flow tortuosity is expected; however, that it increases with increasing storm duration could not have been foreseen. A plausible explanation is that tracer paths have more opportunity to evolve during longer-duration storms, with the tracer trajectories influenced by both hydrologic and hydraulic controls along the hillslope. This explanation is consistent with the fact that the resistance contrast between vegetated and  bare soil areas increases with increasing Reynolds number Re, producing greater flow tortuosity η (see Figure 3b) and suggesting that connectivity lengthscales do reflect the tendency for flow to be diverted around vegetation.
Comprehensive testing across variation in all relevant variables was not conducted, and Equation 6 may be sensitive to such factors as antecedent soil moisture (e.g., impacting sorptivity effects) or microtopography. For example, soils prone to erosion and rill formation can promote microtopography changes that divert runoff around vegetation patches (Chen et al., 2013), decreasing run-on inputs to vegetation and increasing path lengths. The results are also limited in scope to storm and hillslope scales; extension to the catchment scale would involve including subsurface transport and stream networks in the modeling framework, and is beyond the scope of the study.
The results presented here inform efforts to explore connectivity-mediated ecohydrological feedbacks and regime shifts in drylands. Studies invoke dual concepts of "global" and "local" connectivity to describe the potential for runoff to remove water and resources from the landscape (global connectivity) versus redistribute these through runoff run-on processes (local connectivity) (Mayor et al., 2019;Turnbull & Wainwright, 2019). Taking ̃ as a proxy for local connectivity, and ̃ as a measure of global connectivity, the present study demonstrates the interdependence of these processes. Insofar as local and global ecohydrological feedbacks influence resource redistribution and vegetation dynamics in opposite directions, such feedbacks are not independent. Figure 2b suggests that ̃ and ̃ are closely related, but that the relation between these terms depends on storm duration. Characterizing the relation between ̃ and ̃ thus may inform attempts to explore trade-offs between "local" and "global" connectivity in studies of dryland vegetation pattern formation and resilience (Mayor et al., 2019).
The study results are also relevant to monitoring strategies, as relations between runoff run-on path lengths and water balance partitioning enables the use of hillslope-scale observations to constrain within-slope processes, and vice versa. For example, total runoff observed at the outlet could be used to estimate source-sink path lengths (̃) , as a proxy for the typical bare soil area contributing run-on to vegetation patches (Svoray et al., 2021). For such an application, further work is needed to predict how source-sink path lengths relate to the hillslope runoff coefficient across storm and hillslope conditions.
Applications also include mathematical models characterizing pattern formation in drylands (e.g., HilleRisLambers et al., 2001;Rietkerk et al., 2002, and others), which typically replace stochastic rainfall events by a constant annual rate; however, this idealized representation is known to influence the predicted pattern morphology (Crompton & Thompson, 2021;Gandhi et al., 2020;Siteur et al., 2014). If C = 0 is assumed-optimized rainfall capture by the vegetation, and equivalent to the periodic boundary conditions typically prescribed in these models-the study approach could be used to parameterize and constrain the distance that surface water travels before it infiltrates into the soil (Gandhi et al., 2020), and thereby insert realistic storms and runoff run-on processes into mathematical models describing dryland vegetation patterns.
For experimental and field applications, the approach presented here-in which tracers were delivered with the rainfall-needs to be extended to describe sediment (particles with inertia), dye or seed tracers (Masselink et al., 2017). With an extension to describe sediment connectivity, the approach presented in this study could additionally be applied to estimate the area of landforms involved in sediment cascades from sediment discharge measured at a hillslope or catchment outlet (Buter et al., 2022). Such applications would likely be context-specific, particularly because source-sink distances are expected to vary across the hillslope.

Conclusions
In a series of hillslope runoff simulations using a SVE model-paired with particle tracing to identify flow paths-the relation between outlet connectivity ̃ and runoff coefficient C was stationary across rainfall and vegetation characteristics, suggesting a consistent power law scaling between ̃ and C. By contrast, the relation between C and source-sink connectivity ̃ exhibits more variability, in particular, sensitivity to the rainfall characteristics. For longer durations and more intense rainfall events, ̃<̃ , suggesting that tracers that escape are preferentially those that get diverted around vegetation patches. The results open up a number of opportunities for future study and management tools, including the use of water balance partitioning to constrain runoff connectivity and vice versa. Further research is required to generate methods to estimate ̃ and ̃ from field studies and predict ̃∼̃ across varying storm conditions.