Vertical Connectivity Regulates Water Transit Time and Chemical Weathering at the Hillslope Scale

How does hillslope structure (e.g., hillslope shape and permeability variation) regulate its hydro‐geochemical functioning (flow paths, solute export, chemical weathering)? Numerical reactive transport experiments and particle tracking were used to answer this question. Results underscore the first‐order control of permeability variations (with depth) on vertical connectivity (VC), defined as the fraction of water flowing into streams from below the soil zone. Where permeability decreases sharply and VC is low, >95% of water flows through the top 6 m of the subsurface, barely interacting with reactive rock at depth. High VC also elongates mean transit times (MTTs) and weathering rates. VC however is less of an influence under arid climates where long transit times drive weathering to equilibrium. The results lead to three working hypotheses that can be further tested. H1: The permeability variations with depth influence MTTs of stream water more strongly than hillslope shapes; hillslope shapes instead influence the younger fraction of stream water more. H2: High VC arising from high permeability at depths enhances weathering by promoting deeper water penetration and water‐rock interactions; the influence of VC weakens under arid climates and larger hillslopes with longer MTTs. H3: VC regulates chemical contrasts between shallow and deep waters (Cratio) and solute export patterns encapsulated in the power law slope b of concentration‐discharge (CQ) relationships. Higher VC leads to similar shallow versus deep water chemistry (Cratio ∼1) and more chemostatic CQ patterns. Although supporting data already exist, these hypotheses can be further tested with carefully designed, co‐located modeling and measurements of soil, rock, and waters. Broadly, the importance of hillslope subsurface structure (e.g., permeability variation) indicate it is essential in regulating earth surface hydrogeochemical response to changing climate and human activities.

responsive to changing climate and human activities; the deeper subsurface is not as sensitive to surface variations but can offer much-needed "" capacity that buffers surface perturbations.
Subsurface properties vary significantly with depth. Shallow soils are weathered materials that are depleted in geogenic solutes such as Ca and Mg (Brantley, Holleran, et al., 2013;Jin et al., 2010) but are organic-rich with plant roots, leaves, and microbiota (Billings et al., 2018;Herndon et al., 2015;Jobbágy & Jackson, 2000). The deeper subsurface includes parent and partially weathered bedrocks that are enriched with geogenic elements. Shallow soils generally have orders of magnitude higher permeability compared to the deeper, unweathered bedrock (Welch & Allen, 2014). As the incipient rainwater enters the subsurface, it tends to flow through the very thin surficial soil, leading to water ages of a few months (Jasechko et al., 2016;Kirchner, 2016b). In contrast, deeper zones harbor much older water, with ages ranging from hundreds of years to millennia (Jasechko et al., 2017).
These distinct shallow and deep characteristics have strong implications for connecting stream chemistry, reaction rates, flow paths, and solute depth profiles (Benettin et al., 2020;Godsey et al., 2009b;Li, Bao, et al., 2017;Torres et al., 2015). Concentrations of geogenic solutes often increase with depth and their concentration-discharge (CQ) relationship exhibit dilution patterns, often become more pronounced with the presence of sharp contrasts between shallow and deep subsurface (Ameli et al., 2017;Zhi et al., 2019). Aged dissolved organic carbon (DOC) in arid basins and base flow often associate with deeper flows paths (Barnes et al., 2018), whereas concentrations of weathering products (e.g., Si) typically decrease with increasing stream discharge (Clow & Drever, 1996;Johnson, 1969). Concentrations of silica have been observed to positively correlate with groundwater age (Burns et al., 1998;Burns et al., 2003) but not with the age of younger stream water estimated using surface water dating techniques (e.g., Frisbee, Wilson, et al., 2013), raising the question of whether surface water dating have largely missed the long tail of TTDs. Along similar lines of ideas, the end-member mixing analysis (EMMA) have long been used by the hydrology community to differentiate the contributions from different source waters to streamflow (Hooper et al., 1990;Pinder & Jones, 1969) (Klaus & McDonnell, 2013). Recent meta analysis further support the shallow and deep hypothesis that vertical chemical contrasts in shallow and deep waters regulateThis idea however has emerged from studies of individual catchments much earlier (Boyer et al., 1997;Creed et al., 1996).
Reactive transport modeling of homogeneous systems have emphasized the importance of the Damköhler number, the relative time scales of reactions versus water flow in regulating weathering and solute chemistry (e.g., (Lebedeva et al., 2015;Maher, 2010;Salehikhoo & Li, 2015). The water transit time here reflects the time that water interacts with the omnipresent reacting minerals in homogeneous systems. In heterogeneous systems where reacting minerals are present only in specific locations, additional reactive transit times or exposure times are needed in the upscaled rate laws in order to account for the water-mineral interactions (Wen & Li, 2017, 2018. Causal and quantitative linkages among hillslope geomorphic and subsurface structure, transit time, and chemical weathering however have not been established (Benettin et al., 2020;Harman & Cosans, 2019;Sanz-Prat, Lu, Amos, et al., 2016;van der Velde et al., 2010).
Here we ask the questions: how and to what extent does hillslope structure (hillslope shape and subsurface permeability) regulate its hydrogeochemical functioning (water transit times, VC, solute export, and weathering)? Which characteristics exert a stronger control for what function? We focus on weathering systems that involve the dissolution of primary minerals (silicates) and precipitation of secondary minerals (kaolinite). In contrast to recent hillslope studies that use a single component (Ameli et al., 2017;Brantley & Lebedeva, 2021;Harman & Cosans, 2019), this work uses a multi-component model. Numerical experiments were carried out for hillslopes of the same homogenous mineralogy but different shape and subsurface permeability. The mean discharge conditions range from 0.02 to 2.00 mm/d, representing arid (annual mean discharge of 7 mm/yr) to wet (700 mm/yr) climate. Comparison among cases can identify the most influential factors for hydrogeochemical functioning.

The Hillslope Model
In the hillslope model illustrated in Figure 1, infiltrating water (=annual rainfall and snow fall − annual evapotranspiration, or P − ET) can flow laterally through the shallow subsurface to the stream, or vertical-ly to the deeper subsurface, eventually leaving the hillslope from deeper depths as streamflow. Along the flow paths, water interacts with soils and rocks, dissolving primary minerals and precipitating secondary minerals. The model ignores variable water saturation in the shallow subsurface and does not simulate the dynamics of storm events and seasonality. Permeability variations in the horizontal direction are not included either. The hillslope model therefore represents long-term averaged behaviors, and does not represent switching dominant flow paths across seasons and events.
CrunchFlow, a reactive transport code, was used to solve for the flow field and reactive transport equations (Steefel et al., 2015;Steefel & Lasaga, 1994). In setting up simulation domains and properties, observations from the Shale Hills catchment were used for baseline values (Brantley et al., 2018). The relief (elevation change from top of hillslope to bottom of hill) is ∼40 m, and the horizontal projection of hillslope length is ∼100 m. The domain has 50 × 40 grids with a resolution of 2 m (horizontal) × 1 m (vertical) × 1 m (thickness). The hillslope domain includes the "rock" grids where flow and reactive transport occur; the "air" grids were included in the simulation domain because of the code requirement for regularly shaped domains. They were set as inactive in all simulations.
The flow field was calculated using Darcy's law based on permeability distribution and boundary flow conditions. The permeability decreases with depth following a widely used exponential function, κ(z) = κ 0 e −αz (Harr, 1977;Weyman, 1973). Here κ is permeability, α represents how fast permeability decreases with depth, and z is depth in meters (Table 1). In all simulations, the minimum permeability values were set at 10 −16 m 2 . Streamflow was a mixture of water fluxes from different depths at grids within 10 m from the stream outlet.

Numerical Experiments
To tease apart the influence of permeability distribution and hillslope shape, we designed 12 cases with different combinations of these two variables (Table 1). All cases have homogeneous minerals. They have linear or non-linear hillslope shapes and variable permeability distribution (α from 0 to 0.6) (Table 1). In natural systems, permeability typically decreases with depth (Welch & Allen, 2014). Here we used κ 0 = 10 −12 m 2 based on observations at Shale Hills (Kuntz et al., 2011). Values of κ quickly drop to XIAO ET AL.  Water flow paths and velocities in two convex hillslopes with homogeneously (left, high vertical connectivity) and heterogeneously (right, low vertical connectivity) distributed permeability. Incoming water enters the hillslope and flows into the stream via the shallow and deep flow paths, dissolving soils and rocks along the way. The flow paths (light gray lines) and velocity fields (colors) are model outputs from two cases with permeability distribution following the equation κ(z) = κ 0 e −αz , where κ is permeability, κ 0 at the ground surface, z is depth (mbls, meters below land surface), α quantifies how fast permeability decreases with depth (see Table 1 in Section 2). The orange line illustrates the position of soil-weathered rock interface. In the homogeneous case with high vertical connectivity (α = 0), the flow lines are relatively evenly distributed with depth and > 55% stream water is from depths deeper than 6 m); in the heterogeneous case where permeability decreases sharply (α = 0.6), vertical connectivity is low and the flow lines mostly concentrate in the shallow soil and less 5% of stream water is from depths deeper than 6 m. High permeability at depth (high VC) enhances deeper water penetration and increases flow velocities at depth, therefore promoting weathering at depths and more dissolved solutes in stream water. 10 −14 m 2 at 6 m depth when α = 0.6, but remain at κ 0 in the homogeneous case (Figure 1b). The heterogeneous cases (with non-zero α values and permeability decreasing with depth) represent characteristics of natural systems with low permeability at depth. In the homogeneous case, permeability at depth is the same as in shallow soil, representing conditions where high permeability can develop in deeper zones, for example, in karst formations (Hartmann et al., 2014), or systems with dense fractures at depth (Somers & McKenzie, 2020) or deepening roots (Nippert et al., 2012). In all cases, model results show that >95% of water flow comes out from the top 6 m (as observed at Shale Hills (Sullivan et al., 2016)) and the remaining water is from depths below 6 m. Five numerical experiments were run with discharge from 0.02 to 2.00 mm/d, representing very dry (annual mean discharge about 7 mm/ yr) to wet (700 mm/yr) conditions.

Weathering Reactions
The model considered a generic set of weathering reactions, which includes plagioclase feldspar as the primary dissolving mineral, kaolinite as the precipitating mineral, and quartz as an inert mineral ( The reaction consumes Al and CO 2 (aq) and releases Na, Ca, and Si. Chloride (Cl) is conservative and originates from rainfall and snowfall. Because of the rate difference between plagioclase dissolution and kaolinite precipitation, Al produced from plagioclase dissolution is not consumed right away. Soil organic carbon (OC) was used as a "proxy" acidity source for the production of CO 2 through soil respiration and the presence of organic acid (Brantley, Lebedeva, & Bazilevskaya, 2013;Li, Maher, et al., 2017). Previous reactive XIAO ET AL.
10.1029/2020WR029207 4 of 21 a Each case was run with five values of discharge (0.02, 0.2, 1.0, 1.5, 2.0 mm/d) and with the same initial and boundary conditions (supporting information, Table S1).  transport modeling has shown that weathering rates without soil CO 2 can be orders of magnitude lower and cannot reproduce soil water chemistry observations (Heidari et al., 2017;Moore et al., 2012). The OC(s) was set at a volume fraction of 5% in the top 2 m and 1% below 2 m with much lower specific surface area (10 −3 times the number for the top 2 m), as organic matter is typically more abundant and reactive in shallow soils. We assumed isovolumetric weathering that ignores contraction or compaction processes. The model outputs of concentrations are within the concentration ranges of soil water, groundwater, and streams in literature (Heidari et al., 2017;Jin et al., 2014;Zhi et al., 2019).

Reactive Transport Model
CrunchFlow has been used previously for one-dimensional systems with constant flow velocities (Heidari et al., 2017;Hubbard et al., 2014;Maher, 2010;Navarre-Sitchler et al., 2011;Navarre-Sitchler et al., 2009), and two-dimensional domains (Bao et al., 2014;Wen et al., 2016). To our knowledge, this is the first time CrunchFlow is used for hillslope simulation. The code solves mass conservation equations, considering advection, diffusion/dispersion, and weathering reactions. A representative equation for solute j in grid i is as follows: where t is time (s), ϕ is porosity, C j,i is the concentration of primary species j in grid i (mol  m −3 water), D is the hydrodynamic dispersion tensor (m 2 porous media  s −1 ), u is the Darcy flow velocity (m 3 water  m −2 porous media  s −1 ), r j,i is the total rate of species j (mol  m −3 porous media  s −1 ) in multiple reactions, and N is the total number of primary species.
Weathering rates follow Transition state theory (Lasaga, 1998) where k 1 , k 2 , k 3 are the kinetic rate constants (mol m 2 s −1 ),  H a and  OH a represent activities (product of activity coefficient and concentrations) of the hydrogen and hydroxyl ions, n1 and n2 describe the rate dependence on activities, A is the surface area per unit volume of porous media (m 2 m −3 ), IAP is the ion activity product, and K eq is the equilibrium constant. The term IAP/K eq quantifies the distance from equilibrium. The model outputs include concentrations of individual solutes and reaction rates in time and space.

Particle Tracking
We followed a previously published particle tacking approach (Ameli et al., (2017); Ameli et al., (2016b)) to track flow paths and to quantify water age. Particles enter the subsurface at time t 0 and follow the flow field, eventually leaving the hillslope at the outlet at time t out (Figure 1). A total of 250 particles was injected into the hillslope at the ground surface. For each particle p, its coordinates at any time step k depend on its coordinates at the previous step k−1, flow velocities in the x and z direction v xr and v zr , and the size of the time step (set as one hour). The coordinates of particle p at the time step k (x p k and z p k ) are: With the random movement of particles, v xr and v zr are random numbers sampled from an exponential distribution with the means of v x and v z being the flow velocities at x p k−1 and z p The t out of each particle at the stream outlet quantifies its total travel time across the hillslope. The distribution of t out for the 250 particles was considered to follow a Gamma distribution that has been extensively used for TTD (Ameli, Amvrosiadi, et al., 2016;Harman, 2015;McGuire & McDonnell, 2006): Here p(t) is the probability density function of Gamma distribution, A is the shape parameter. Note that here A is often referred to as α in the Gamma distribution equation. We use A instead of α to avoid confusion XIAO ET AL.
10.1029/2020WR029207 6 of 21 Figure 2. Spatial distribution of flow lines, weathering rates, and solute concentrations in hillslopes with different permeability distributions, slope shapes, and discharge rates. All cases have homogeneously distributed reacting minerals. The models were run to steady state when solute concentrations no longer change with time. Plagioclase is weathered to kaolinite during exposure to water flowing through the hillslope. Except in homogeneous hillslopes (α = 0), the majority of the water (>95%) flows through the top 6 meters of the subsurface. Water chemistry in the top 6 m is therefore most influential in determining stream chemistry. The concentrations of Na, Ca, and Si are lower in shallow zones (top 2 m), and Al concentrations are higher in shallow soils. Chloride concentrations are relatively similar within the top 6 meters.
with α in the permeability distribution equation. The β is the rate parameter (the reciprocal of the scale parameter), t is the transit time, and     is the Gamma function. The product of A and β is the statistic mean of t , that is, the mean transit time (MTT). The Gamma distribution is positive and skewed. The shape parameter A describes the weight in the tails of the distribution. The smaller A, the larger variability in transit times with higher proportions of young and old waters.

Transit Time, Young and Old Water Fractions, and Water Age
The transit time is defined as the difference between the time when a particle enters and the time when a particle leaves the hillslope and enters the stream (McGuire & McDonnell, 2006;Rodhe et al., 1996). The TTD is the probability density function of the transit times of all particles (McGuire & McDonnell, 2006). The MTT is the product of the two parameters (μ = A × β) in Equation 5. The Water age is the average transit time of particles in particular locations (e.g., stream). The young water fraction (F yw<1yr ) is defined as the fraction of particles (out of all 250 particles) in the stream with transit times shorter than 1 year. The corresponding old water fraction is F ow>1yr = 1 − F yw<1yr . Note that this differs from the typical young water fraction defined as less than 2.3 months, or 0.2 years (Kirchner, 2016a;Kirchner, 2016b). As discussed later, this one year cutoff for young and old water fractions is somewhat arbitrary but roughly matches the cutoff between fast-and slow-flow through Shale Hills (Sullivan et al., 2016). The threshold of 0.2 and 5 years were also used, as shown in Figure S2 in the supporting information and discussed in the Results section.

Vertical Connectivity
Vertical connectivity is a measure of the connection between shallow soil and deeper groundwater. Their values were quantified as the water fluxes below the soil-rock interface (deeper than the soil zone, averaged across the entire hillslope) divided by the total water fluxes at all depths (into the stream). Higher VC means higher water fluxes into (and out of) the deeper zone. The top few meters of the subsurface were considered as soil. Three depths (1 m, 2 m, and 3 m) were used as the depths of the soil-rock interface for analysis. As will be shown later, different soil depths lead to similar VCs. We therefore show the analysis for the 2 m soil depth in the Results section. The results for 1 and 3 m are in supporting information, Figure S1.

Concentration-Discharge Relationships
Stream concentrations (C stream ) were calculated as the flux-weighted average over the entire depth. The grids within 10 m from the outlet were used to avoid boundary effects ( Figure 1). Stream CQ relationships describe concentration variations responding to changes in mean streamflow (Godsey et al., 2009a;Herndon et al., 2018;Moatar et al., 2017). The simple power law C = aQ b was used to classify CQ patterns of chemostatic (b ∼ 0), dilution (b < 0), and flushing (b > 0) patterns.

Solute Concentrations in Shallow and Deep Zones
The model output includes solute concentrations in individual grids. Shallow soil water concentrations (C sw ) were calculated as the volume-weighted average concentrations from the top 2 m of the soil; groundwater concentrations (C gw ) were computed as the volume-weighted average concentrations from below 2 m. The concentration contrast C ratio was calculated as C sw /C gw . Analysis was also done using 1 and 3 m as the soil depth but did not yield significantly different MTT and C ratio values (supporting information, Figure S2).

Hillslope Weathering Rates
The summation of solute fluxes (or loads, = concentration (C) × discharge (Q)) at the stream outlet quantifies weathering fluxes at the hillslope scale. Some may argue that these quantities are solute fluxes instead of weathering rates. In natural systems where other geochemical processes (e.g., sorption or ion exchange) occur, the two quantities can differ. They may also differ at the event or daily scales with transient flow fields. Here as we simulate weathering without other reactions and under steady state flow conditions, the solute fluxes and weathering rates are expected to be the same. This is further confirmed by the resultant same values arising from C × Q at the outlet and the sum of local rates at individual grids.

Effects of Permeability Distribution
Here we discuss two linear hillslopes with different permeability distributions. In each hillslope, water particles that enter the valley generally have shorter flow paths compared to those that enter the ridge top. In the homogeneous case (α = 0), permeability remains high at depth; in the heterogeneous case (α = 0.6), permeability decreases by two orders of magnitude within the top 5 meters. The flow lines in the heterogeneous case are concentrated in the shallow, high-permeability subsurface (Figure 2). The flow lines in the homogeneous case, especially those starting at the ridge top, penetrate much deeper. These flow paths are much longer and water travels through a much larger volume of the subsurface than those in the heterogeneous case.
In the homogeneous hillslope, water fluxes are homogeneous with depth; in the heterogeneous hillslope, water fluxes concentrate mostly within the top 5 m (Figure 3). Water fluxes decrease and water ages increase with depth. Much more water exits the hillslope at shallow depths in the heterogeneous case. The water fractions exiting the top 2 m are 11% and 64%, respectively, in the cases of α = 0 versus α = 0.6. The corresponding water fractions for the top 6 m are 38% and 95%, respectively. The water leaving via shallow flow paths is considerably younger than that via deeper paths. The TTDs differ significantly in different cases (Figures 3c and 3d). The heterogeneous hillslopes have short transit times with much larger fractions of young water compared to the homogeneous hillslopes.

Effects of Hillslope Shape
The hillslope shape does not have as much influence as permeability variations. The Concave case has the highest fraction of younger water (F yw<1yr ), and the Convex case has the lowest F yw<1yr (Figures 3g and 3h). The differences among different hillslope shapes, however, are much smaller than the differences among different permeability distributions with the same hillslope shape (Figures 3c and 3d). The hillslope shapes generally have much stronger influences in the homogeneous case, and in younger water end of the TTD (below 1 year). The older water end of the TTD almost overlaps among different hillslopes. Figure 4a shows that as α values vary from 0.6 to 0.0, VC increases. This is expected because more water penetrates to greater depths at low α values. The VC values however are almost identical in hillslopes of different shapes but identical α values, indicating negligible impacts of hillslope shape. Values of MTT (Equation 5), increase from ∼1 year to ∼10 years as α decreases. In the homogeneous cases (α = 0) with high VC, the hillslope shape has large impact on MTT: The Concave case has an MTT value of 4.1 years compared to 9.5 years for the Convex case. At α = 0.6, however, MTT values in these hillslopes overlap, indicating negligible impact of hillslope shape (Figure 4b). The younger water fraction F yw<1yr generally is much more influenced by hillslope shape than old water. Different hillslope shapes lead to distinct young water fractions, with the Convex and Concave cases resulting in the highest and lowest values, respectively. Supporting information, Figure S1 show similar trends for MTT versus VC at soil depths of 1 m, 2 m, and 3 m. The young water fractions with different threshold values (0.2, 1, and 5 years) also show similar trends with respect to dependence on VC (supporting information, Figure S2). Figure 2 shows negligible effects of hillslope shape: weathering rates and solute concentrations are similar among the linear, concave and convex cases (Figures 2b-2d). The permeability distribution however strongly regulates flow paths, weathering rates, and solute distribution. The flow lines penetrate as deep as 20 meters in the homogeneous hillslope, but only to ∼5 meters in the heterogeneous case (Figure 2, row 1). Plagioclase generally dissolves faster in shallow compared to deeper zones, because newly infiltrated water in the shallow zones are farther away from equilibrium conditions (row 2, plagioclase rates). With the faster release of weathering products in shallow zones, kaolinite precipitation also occurs at faster rates (row 3). The shallow zones are therefore where weathering rates are highest; reaction rates plummet by orders of magnitude beneath the shallow zones. With deeper water penetration, weathering can occur at greater depths in the homogeneous hillslope.

Hillslope Distribution of Weathering Rates and Solutes
Weathering transforms plagioclase to kaolinite (Equation 1). Plagioclase dissolution releases Na, Si, Al, and Ca, but kaolinite precipitation consumes Si and Al. Plagioclase releases more Si and less Al than the amounts consumed in kaolinite precipitation, as indicated in the Si:Al stoichiometry in Equation 1. The net weathering reaction therefore consumes Al but releases Si, Na, and Ca. Transforming one mole of plagioclase to one mole of kaolinite consumes 0.387 mole of Al and generates 0.387 moles of Na, Si, and Ca. As such, concentrations of Na, Ca, and Si increase along flow paths, leading to higher concentrations at depth in the top 6 meters (Figures 5a and 5b). Si decreases with depth beyond 6 m in some of the most heterogeneous hillslopes, indicating plagioclase dissolution rates are lower than kaolinite precipitation rates under these conditions with very low flow. In contrast, Al is consumed as kaolinite precipitates along flow paths such that its concentrations decrease with depth ( Figure 5c). The non-reactive Cl is homogeneously distributed over all depths in the homogeneous cases and in the top 6 meters of the heterogeneous cases (Figure 5d

Effects of Discharge
As discharge decreases from 1.0 to 0.2 mm/day (Figures 2 and 5), weathering rates decrease but solute concentrations increase. In addition, at lower discharge, the effects of permeability on depth profiles lessen. Figure 6a shows that plagioclase dissolution rates increase as discharge increases. Kaolinite precipitation rates mirror those of plagioclase dissolution, although at a slightly lower magnitude. The net weathering rates are positive and increase with discharge, because shorter transit times under high discharge maintain the system away from equilibrium. The homogeneous cases generally yield higher weathering rates and higher concentrations of Na, Si, and Ca compared to heterogeneous cases. Solute concentrations vary less as discharge varies. Concentrations of Al increase with discharge in the heterogeneous cases. The concentrations are at levels commonly observed in natural systems. For example, Al can reach as high as tens of μmol/L that is, within the high end of observed concentration ranges (Herndon et al., 2015;. The concentrations of Na, Si, and Ca exhibit dilution CQ patterns with higher concentrations at low discharge ( Figure 6). The Al concentrations demonstrate flushing patterns with high concentrations at high discharge, because the mobilized solutes are flushed out rapidly without as much kaolinite precipitation under high discharge. The concentrations of Cl exhibit a chemostatic pattern with little variation, ultimately converging to the infiltrating water concentrations (1.2 × 10 −5 mol/L) at high discharge. Figure 5a suggests that total weathering fluxes (C × Q, the same as export rates) are higher at higher discharge, although stream concentrations of Na, Si, and Ca are lower in these cases.
The CQ patterns do not change under different permeability distributions or hillslope shape, although the extent of dilution or flushing (the value of the slope b) may differ (Figure 7). For example, Na exhibits a dilution pattern in all cases but with much steeper CQ at higher α values (lower VC), indicating more negative b values (Figure 7a). This may be due to increasing water fluxes through the top zone soil at high Q. Similarly, Al exhibits flushing patterns in all cases but exhibit more positive b values in the more heterogeneous cases. The heterogeneous cases generally have XIAO ET AL.    Figure 2). The CQ relationships of Na, Si, and Ca suggest dilution (decreasing C with increasing Q); Al exhibits a flushing pattern (increasing C with increasing Q); Cl exhibits a chemostatic pattern (invariable C with Q). The dilution and flushing patterns are more pronounced (larger absolute b values) in the heterogeneous cases. more pronounced, steeper CQ relationships. Compared to the effect of variations in the distribution of permeability, variation in hillslope shape results in negligible changes in CQ patterns. The Concave hillslope notably has a more significant impact than other hillslope shape.
The CQ power law slope b values (in C = aQ b ) vary with different measures of transit time and VC (Figure 8). In all cases, Na, Al, and Cl exhibit dilution, flushing, and chemostatic patterns, respectively, regardless of the permeability distributions. The b values are negative at low C sw /C gw ratios (<1.0), are close to zero at C sw /C gw ratios roughly equal to 1, and are positive at high C sw /C gw ratios (>1.0) ( Figure 8a). As VC increases, the shallow and deep contrasts decrease such that the b values approach 0 ( Figure 8b). Similarly, when MTT and F ow increase, b values approach 0. In other words, high VC indicates more connected shallow and deeper zones, reduced shallow-versus-deep contrast, and more chemostatic CQ relationships. Figure 9 shows that weathering rates generally increase with VC and discharge. The magnitude of the increase depends on discharge. Under dry conditions (0.02 mm/d), weathering rates vary negligibly (Figures 9b  and 9c). As conditions get wetter, the dependence on VC becomes more significant. At each Q, weathering rates increase with F ow (Figure 9b) and with MTT (Figure 9c). The stream concentrations of Na and Si correlate with MTTs under non-arid conditions (Q > 0.2 mm/day) via a power law relationship C stream = C 0 (MTT) ν (Figure 10). For Na, C 0 = 10 −3.4 mol/L and ν = 0.58 (R 2 of 0.90); for Si, C 0 = 10 −3.1 mol/L and ν = 0.53 (R 2 = 0.88). These correlation equations do not include concentrations under arid conditions (Q ≤ 0.02 mm/day). Under these arid conditions, concentrations become relatively constant at 10 −2.6 , 10 −2.0 , 10 −6.9 , and 10 −1.8 mol/L.

Weathering Over Geological Time
To explore the effects of permeability distributions over geological time, the models were run at 1.00 mm/d for 10,000 years for the Convex, homogeneous (α = 0) and heterogeneous (α = 1.0) cases. Note that here we simulate a "stationary" hillslopes where supply from below or uplift equals soil erosion such that these two competing processes "cancel" each other and are not considered in the model. Figure 11 shows the τ values (mass transfer coefficient) for plagioclase dissolution, where τ = 0 indicates the same plagioclase abundance as the parent bedrock and τ = −1.0 represents 100% depletion of plagioclase compared to the parent rock (isovolumetric weathering). In the homogeneous hillslope, water penetrates deeply and has longer water-rock contact times, therefore generating a deeper reaction front where τ value change from −1 to 0. The depths of the weathering front at 10,000 years are 11 m, 10 m, and 4 m, for the ridge top, mid slope, and valley positions, respectively. In contrast, in the heterogeneous case, most water flows through a thin soil layer of high permeability. After 10,000 years, weathering fronts are situated at 4 m, 5 m, and 7 m, for ridge top, mid slope, and valley floor, respectively. In other words, the weathering fronts are shallower in the heterogeneous hillslopes. Overall, the total weathered mass for the entire hillslope is 3.70 × 10 5 mol in the homogeneous case, compared to 1.99 × 10 5 mol in the heterogeneous case. In other words, weathering rates almost doubled in the homogeneous case. In addition, the homogeneous Hillslope develops the deepest weathering front at the ridge top, as the water penetrates deepest there; the heterogeneous case develops the deepest front in the valley, where all water is eventually channeled. XIAO ET AL.

Discussion
This work used numerical thought experiments to understand the relative importance of different hillslope structure (permeability variation with depth and hillslope shape) in regulating hydrogeochemical functioning. The experiments made a few simplification assumptions. The hillslopes were modeled using relatively simplified, albeit multicomponent weathering reaction with dissolving primary mineral (plagioclase) and precipitating secondary mineral (kaolinite). The minerals were homogeneously distributed with depth in order to focus on effects of physical structure. The experiments used steady state flows that represent average hydrodynamic conditions and ignored details of unsaturated and saturated flows. The experiments therefore do not represent event-scale dynamics such as transient horizontal connectivity, antecedent conditions, and variable storm sizes (Hopp & McDonnell, 2009;Tromp-van Meerveld & McDonnell, 2006;Zimmer & McGlynn, 2017). Instead, they represent more of hillslopes of same homogeneous distributed mineralogy but different shape and permeability distribution under different hydroclimate conditions from arid to wet climates. These experiments linked hillslope horizontal and vertical structure with their hydrological (flow paths and transit times) and geochemical functioning (weathering and solute export). As we discuss below, the model predictions are consistent with many field observations in literature, indicating the model represents essential processes that captured first-order dynamics at the hillslope scale. Here we distill the results into three hypotheses, discuss existing supporting evidence, and lay out future work that can test these hypotheses.

Hypothesis 1
The variations of permeability with depth more strongly influence the mean transit times of stream water than hillslope shape does; hillslope shapes instead influence more of the younger water fractions of stream water.

Existing Evidence and Implications
In hillslopes with high permeability at depth (α = 0), more water penetrates deeper, leading to more evenly distributed water fluxes from different depths, high VC, high MTT, and low F yw (Figures 3 and 4). About 40% of infiltrated water enter the stream via the top 6% and and 12% via top 2 meters of the subsurface, respectively. Where permeability decreases sharply (α = 0.6), ∼95% of infiltrated water enters the stream via the top XIAO ET AL.
10.1029/2020WR029207 12 of 21 Figure 8. The concentration-discharge (CQ) power law slope b as a function of (a) C ratio (=C sw /C gw ), the ratio of concentrations in shallow soil water (C sw ) and those in groundwater concentration (C gw ), (b) VC, (c) mean transit time (MTT), and (d) old water fraction F ow>1yr . The CQ power law slope b was calculated based on C stream at Q from 0.02 to 2.0 mm/d. The values of b generally increase with C sw /C gw ; they gradually approach zero (chemostasis) as vertical connectivity, MTT, and F ow>1yr increase. 6 meters and ∼65% via the top 2 meters (Figure 4). These fractions of deeper water, or the complementing percentages that add up to one, are within the range of base flow indices reported in literature (Somers & Mc-Kenzie, 2020). The results underscores the essential roles of permeability variation as regulators of shallow and deep water partitioning.
MTTs have been shown to reflect the long tails of skewed TTDs and are not good predictors of seasonal cycles of non-reactive tracers such as Cl (Kirchner, 2016a;Kirchner, 2016b). Instead, F yw<0.2 yrs has been shown to be a better measure of seasonal variations and topographic characteristics (Kirchner, 2016a;Kirchner, 2019;von Freyberg et al., 2018). Results here underscore that F yw and MTTs in fact reflect different metrics of TTDs and distinct flow paths. Values of F yw reflect more of the shallow, fast flows shaped by hillslope shape whereas MTTs signal more of deep flow paths regulated by subsurface permeability variation. These observations corroborate existing literature on interacting controls of topography and permeability on transit times and flow paths. For example, steep catchments have been observed to have lower transit times; such topographic control weakens in flat and subdued terrains, where transit times depend more on subsurface permeability (Tetzlaff et al., 2009).
The correlations between MTTs and permeability distribution indicate that MTTs can be a good measure of the properties of hydrologically active soil and bedrock (Asano & Uchida, 2012). These correlations from modeling output echo observed correlations between MTTs, mean length of flow paths (McGuire & McDonnell, 2010), and Si concentrations that reflect water-rock contact time (Asano & Uchida, 2012;. Low-permeability granite hillslopes have been observed to have a higher fraction of young, event water in the stream, as most water flows through the shallow layer. In contrast, hillslopes overlaying fractured, high-permeability shale had more older water due to the longer flow paths via deep bedrock (Onda et al., 2006). Similarly, longer MTTs have been observed in catchments with permeable, fractured sandstone bedrock compared to those with tight, volcanic bedrock (Hale & McDonnell, 2016;Pfister et al., 2017).
Over geological time, under similar tectonic and soil erosion conditions, thick soils can develop under high VC conditions ( Figure 11). Deeper water penetration also connects a larger proportion of the subsurface to the stream and induces larger dynamic water storage (Berghuijs and Kirchner, 2017;Birkel et al., 2015;Buttle, 2016;McGuire & McDonnell, 2010;McNamara et al., 2011;Pfister et al., 2017;Remondi et al., 2018). Although not explicitly examined here, the relationship between subsurface permeability and dynamic storage can have important implications for storage-discharge relationships, streamflow generation, and the capacity to "buffer" external perturbations.

Hypothesis 2
High vertical connectivity (VC) arising from high permeability at depth enhances weathering because it promotes deeper water penetration and water-rock interactions in wet climates; the influence of VC weakens under arid conditions and in large hills with longer mean transit times.

Existing Evidence and Implications
Under wet conditions, weathering rates increase with both old water fractions F ow and with MTTs (Figure 9). Solute concentrations correlate to MTTs via a power law relationship ( Figure 10). VC however is less of an influence under droughts and arid climates where slow flows and long transit times drive weathering to equilibrium. This hypothesis echoes the Wen and Li upscaled rate law that emphasizes two influential factors of weathering (Equation 14 and Figure 11 in Wen and Li (2018)). One is the total water flow, or the mean discharge in this work, that determines transit time and the Damköhler number Maher & Chamberlain, 2014;Maher & Druhan, 2014;Salehikhoo & Li, 2015;Salehikhoo et al., 2013). As mean discharge is regulated by water input from precipitation (rainfall and snowfall) and water demand by evaporation and transpiration (discharge = precipitation-evapotranspiration) (Anderson et al., 2019;Wlostowski et al., 2021;Zapata-Rios et al., 2015), this in fact reflects the first order control of hydroclimatic and vegetation conditions. The other factor in the Wen and Li rate law is the water fraction that interacts with reacting minerals (Wen & Li, 2017;Wen & Li, 2018). This water fraction or water partitioning is determined by the heterogeneity factor, that is, permeability contrast between reactive and non-reactive zones (Wen & Li, 2017;Wen et al., 2021). In this work, the lateral versus vertical partitioning is quantified by VC values and is regulated by permeability variation with depth. High permeability at depth and VC enable the water contact with a larger volume of reacting minerals and reactive surface area, and therefore high weathering rates. Where more water can penetrate deep, stream water is older, resulting in higher concentrations of geogenic solutes ( Figure 10). This concentration-water age relationship in fact echo many field observations that link water age dating to stream concentrations (Benettin et al., 2015;Frisbee, Wilson, et al., 2013;Rademacher et al., 2001). Natural conditions that can enhance water-conducting capabilities at depth include deepening plant roots (Beven, 2018;Hasenmueller et al., 2017;Wen et al., 2021), dense fractures (Brantley & Lebedeva, 2021;Brunet et al., 2016;Gu et al., 2020;Pandey & Rajaram, 2016;Wen et al., 2016), and large conduits in karst formation (Hartmann et al., 2014;Husic et al., 2019;. The Wen and Li rate law also predicts that the influence of the heterogeneity factor (or VC) weakens when water transit times are high. Model results here show that permeability distribution does not matter much under arid condition (long MTTs) (Figure 8), which echoes the prediction from the rate law (Wen & Li, 2018). In large hillslopes with long flow paths, long MTTs can also drive weathering to equilibrium such that the control of VC diminishes. The Damköhler coefficient that explains CQ relationships of geogenic solutes (Ibarra et al., 2016;Maher & Chamberlain, 2014;Torres & Baronas, 2021;Wymore et al., 2017) indicate that long transit times lead to chemical equilibrium and chemostatic patterns. It implies that the dependence of weathering rates on VC and MTT will weaken as spatial scales increase. Some threshold spatial scales may XIAO ET AL. exist, beyond which VC and subsurface permeability do not matter. Such threshold lengths will likely vary with minerals of distinct reaction thermodynamics and kinetics and specific conditions.

Hypothesis 3
Vertical connectivity (VC) regulates chemical contrast between shallow and deep waters (C ratio ) and solute export patterns encapsulated in the concentration-discharge (CQ) power law slope b. Higher VC leads to similar shallow and deep water chemistry (C ratio ∼ 1), leading to more chemostatic CQ patterns.

Existing Evidence and Implications
The solute export patterns, or CQ power law slope b values, are generally consistent with solute depth profiles (Figures 2 and 5). Concentrations of Na, Si, and Ca increase with depth as weathering occurs along flow paths; Al concentrations decrease with depth, as Al continues to be depleted along flow paths by kaolinite precipitation. We therefore observe dilution patterns for Na, Si, and Ca and flushing patterns for Al. For Cl, its non-reactive nature and atmospheric source lead to relatively similar shallow and deep water chemistry and chemostatic patterns (Figure 7).
These results support the shallow and deep hypothesis arising from watershed reactive transport modeling and continental scale metaanalysis Zhi et al., 2019). Namely, the shallow-versus-deep concentration contrasts determine solute export patterns, and CQ power law slope b values can be predicted based on C ratio . In general, solutes enriched in shallow zones demonstrate flushing patterns, whereas those abundant in deep zones show dilution patterns. This hypothesis echoes the use of EMMA and distinct geochemical signatures for hydrograph separation and source water differentiation (Hooper et al., 1990;Klaus & McDonnell, 2013;Pinder & Jones, 1969), although these approaches do not explicitly emphasize the depths of source waters and the deep and shallow water chemistry contrasts. The shallow and deep hypothesis is supported by extensive CQ data from individual sites and from across scale analysis (Creed et al., 1996;Ebeling et al., 2021;Moatar et al., 2017;Musolff et al., 2017). A recent global data analysis also shows the depth of solute generation, not climate drivers, determines CQ slope b (Botter et al., 2020).
A natural follow-up question is then what determines the shallow-versus-deep concentration contrast (C ratio ). This work underscores the predominant role of physical heterogeneity (permeability distribution and VC), even in chemically homogeneous hillslopes. Low VC stratifies the flows via different depths, leading to distinct chemistry in older waters in deeper zones and in younger waters in shallow soils (Jin et al., 2014;Richardson et al., 2020;Sullivan et al., 2016). Large shallow and deep water contrast and non-chemostatic (dilution and flushing) patterns in these hillslopes suggest that chemical heterogeneity is required for non-chemostatic patterns, which is similarly observed in other studies (Torres & Baronas, 2021). This may appear surprising but align with reactive transport principles: Concentration gradients do develop even in idealized, homogeneous systems, as water continues to interact with minerals and dissolve out solutes along its flow paths, although mineral heterogeneity may intensify concentration gradients. In natural hillslopes where parent rocks are often at depth, CQ relationships may exhibit more pronounced dilution patterns with more negative b values. Solutes that are enriched in shallow soils, such as DOC, almost always show flushing CQ patterns (Zarnetske et al., 2018). Metals that form complexes with DOC often exhibit flushing patterns (Herndon et al., 2015;Kiewiet et al., 2020;McIntosh et al., 2017).
The hillslope model here shows that long-term (e.g., annual) average CQ patterns of most geogenic solutes exhibit dilution patterns, similar to observations in recent global data . Earlier CQ data have shown chemostatic patterns for geogenic solutes (Godsey et al., 2009b) and for nutrients (Basu et al., 2010). It is possible that chemostatic patterns from earlier data result from infrequent (often bimonthly) data mostly in low and intermediate flow regimes. Recent, more frequently sampled stream chemistry data rarely show chemostatic patterns (Knapp et al., 2020;Torres et al., 2015;. At minute and hourly time scales, highly resolved temporal solute data have demonstrated highly dynamic, contrasting CQ patterns in events of different size and antecedent conditions (Burns et al., 2019;Duncan et al., 2017). High frequency stream chemistry data may reflect distinct source water chemistry at different depths that relate to rising water table during the events. In other words, different storm events "sample" source waters with distinct chemistry from different depths of water table under different hydrological conditions. Broadly, it is possible that stream chemistry may "mirror" the subsurface water chemistry, and stream chemistry at different time scales (events, seasonal, annual) may reflect subsurface water chemistry at different spatial resolutions.

Further Hypothesis Testing
Further testing of these hypotheses requires co-located measurements in individual sites with collated hydrogeochemical functioning data (stream chemistry and water age) and structure characterization (permeability distribution and hillslope shape). Water chemistry measurements will need to include soil water and groundwater at different depths as well as stream chemistry data at sufficient frequency covering the full spectrum of discharge regimes. Geophysical characterization, soil and rock chemistry, and physical measurements can help us "see" the subsurface, and they often require time-consuming and labor-intensive sampling campaigns (St. Clair et al., 2015). Multiple measures of TTD can be quantified using different tracers, with stable water isotopes targeting water fractions at the time scales of months to years, and groundwater dating techniques (e.g., radioactive carbon) and geochemical tracers for older water at ages ranging from tens to thousands of years (Frisbee, Wilson, et al., 2013). As indicated in the relationships between MTTs, stream chemistry, and VC, benchmark relationships can potentially be established between these data, such that stream chemistry can be used to mirror sparsely measured subsurface properties and weathering rates. The juxtaposition of different datasets can offer a multi-faceted view and reveal structure-function relationships in intensively measured sites.
The collation of data across different sites can then shed light on the first-order controls across distinct hillslope characteristics and climate conditions. Testing these hypotheses across sites requires carefully chosen catchments that differ in hillslope shape or permeability distribution but have similar conditions for other variables. Such studies have been done previously but are often challenging, as natural systems rarely differ only in one variable. Other variables, such as drainage area and soil properties, inevitably differ and mask the signature of compared variables (Hale & McDonnell, 2016;Li et al., 2018). Under such conditions, numerical thought experiments using process-based models, paired with carefully collected data, can help differentiate the effects of confounding variables, as has been shown before (Berghuijs & Kirchner, 2017;Remondi et al., 2018;Xiao et al., 2019).

Conclusions
We used a 2D hillslope reactive transport model and particle tracking to answer the question: How does hillslope structure (hillslope shape and permeability variation) influence hillslope hydrogeochemical functioning? Results show that permeability distribution is the first-order control over VC, flow paths, transit time, and MTT. In hillslopes with low permeability at depth, >65% of water is channeled through the top 2 m and >95% is channeled through the top 6 m of the subsurface, respectively (Figure 3). Where high permeability is sustained at depth, water penetrates deeper and interacts with parent rocks, leading to high VC, MTT, and larger fractions of old water (F ow>1yr > 50%). In contrast, the hillslope shape, a proxy to topography, more strongly influences the TTD of younger water that enters streams (Figure 4). These results imply that MTT and old water fractions, often considered poor metrics of watershed characteristics such as topography and hydroclimatic conditions, may encode essential information about subsurface structure and VC.
Model results also show the decisive role of permeability variation in influencing geochemical functioning. High permeability at depth and VC enhance weathering, especially under wet climate conditions. Higher VC reduces shallow versus deep solute chemistry contrast and leads to more chemostatic export patterns (Figures 5-8). The results support the shallow and deep hypothesis that vertical water chemistry contrasts regulate concentration discharge (CQ) patterns. Concentrations of Cl vary negligibly with depth, resulting in chemostatic CQ patterns. Concentrations of Na, Si, and Ca increase with depth, leading to dilution CQ patterns; Al concentrations decrease with depth, leading to flushing patterns. Higher VC can reduce the shallow and deeper water chemistry contrast, leading to CQ patterns toward chemostasis. These results underscore the important linkage betweèn hillslope structure (permeability distribution) in regulating hydrological (flow paths and TTD), and geochemical functioning (weathering and solute export). They underscore the potential of using streams as mirrors to infer subsurface biogeochemical reactions and scarcely measured subsurface physical and chemical properties. Establishing such linkage is essential for understanding and projecting the response of earth surface processes to future climate and human activities.

Data Availability Statement
Simulation data from this work will be archived in the Consortium of Universities for the Advancement of Hydrologic Science, Inc. (CUAHSI) data website (http://www.hydroshare.org/ resource/8362807be163491dbf3df933c37a16ef).

Acknowledgments
This work was financially supported by funds from the National Science Foundation for the Susquehanna Shale Hillslopes Critical Zone Observatory (SSHCZO) to S. L. Brantley (EAR 13-31726). D. Xiao was additionally supported by funds from the Penn State Earth and Environmental Systems Institute. This research was conducted in Penn State's Stone Valley Forest, which is funded by the Penn State College of Agriculture Sciences and the Department of Ecosystem Science and Management and managed by the staff of the Forestlands Management Office. We acknowledge discussions with members in the SSHCZO and Brantley and Li Reactive Water research group. We are deeply grateful for the help from Brandon Forsythe for data organization and management. We are thankful for the careful text editing by Bryn Stewart from the Li Reactive Water Group. Constructive comments from two anonymous reviewers and the associate editor has significantly improved the manuscript. We particularly acknowledge Dr. Daniel Ibarra for his suggestion of framing the discussion into hypotheses, which has significantly sharpened the paper.