Catchment Coevolution and the Geomorphic Origins of Variable Source Area Hydrology

Features of landscape morphology—including slope, curvature, and drainage dissection—are important controls on runoff generation in upland landscapes. Over long timescales, runoff plays an essential role in shaping these same features through surface erosion. This feedback between erosion and runoff generation suggests that modeling long‐term landscape evolution together with dynamic runoff generation could provide insight into hydrological function. Here we examine the emergence of variable source area runoff generation in a new coupled hydro‐geomorphic model that accounts for water balance partitioning between surface flow, subsurface flow, and evapotranspiration as landscapes evolve over millions of years. We derive a minimal set of dimensionless numbers that provide insight into how hydrologic and geomorphic parameters together affect landscapes. Across the parameter space we investigated, model results collapsed to a single inverse relationship between the dimensionless relief and the ratio of catchment quickflow to discharge. Furthermore, we found an inverse relationship between the Hillslope number, which describes topographic relief relative to aquifer thickness, and the proportion of the landscape that was variably saturated. While the model generally produces fluvial topography visually similar to simpler landscape evolution models, certain parameter combinations produce wide valley bottom wetlands and non‐dendritic, trellis‐like drainage networks, which may reflect real conditions in some landscapes where aquifer gradients become decoupled from topography. With these results, we demonstrate the power of hydro‐geomorphic models for generating new insights into hydrological processes, and also suggest that subsurface hydrology may be integral for modeling aspects of long‐term landscape evolution.


Motivation
Landscape geomorphology is inextricably connected to runoff generation.Topographic slope is often a strong predictor of hydraulic gradient (Haitjema & Mitchell-Bruker, 2005), whereas topographic curvature affects how water is concentrated or dispersed as it moves downslope (Lapides et al., 2020;Prancevic & Kirchner, 2019;Troch et al., 2003), all of which affects the likelihood of surface runoff.Subsurface porosity and permeability further affect these quantities, as they affect the subsurface capacity to infiltrate water and transmit it laterally toward streams (Horton, 1933;O'Loughlin, 1981).At the same time, runoff can alter geomorphic properties of landscapes because it drives erosion and ultimately the incision of river channels, which then affect the morphology of adjacent hillslopes (Callahan et al., 2019;Dietrich et al., 2003;Roering et al., 2001).
Landscape evolution models (LEMs) are essential tools for understanding topographic change over long timescales (e.g., reviews by Bishop (2007), Chen et al. (2014), Pelletier (2013), Tucker and Hancock (2010), Valters (2016), andWillgoose (2005)) and thus could be useful for understanding relationships between topography and runoff.However, such models usually simplify hydrology such that feedbacks between topography and runoff dynamics cannot be examined.Recent studies have made progress in representing hydrologic processes more explicitly in LEMs, and show that drainage density scales linearly (Luijendijk, 2022) or non-linearly (Litwin et al., 2021) with transmissivity when runoff is generated by saturation excess overland flow.Although these studies broke new ground by revealing how runoff generation affects topography, it is still unclear how this coevolution affects hydrological function.Understanding how current hydrological function is affected by landscape history has the potential to transform how we understand Earth's critical zone, and how we make hydrological predictions (Harman & Troch, 2014;Singha & Navarre-Sitchler, 2022;Troch et al., 2015).
Relationships between various stores and fluxes are used to describe hydrological function.One of the most fundamental is the catchment water balance, which describes the long-term partitioning of water into storage, evapotranspiration, runoff, and deep recharge.On shorter timescales, storage-discharge relationships aid interpretations of rainfall-runoff response and catchment recession (McMillan, 2020).Hydrological function also has a spatial component, including the dynamics of saturated areas and flowing stream length (Jensen et al., 2017;Latron & Gallart, 2007;Prancevic & Kirchner, 2019;Warix et al., 2021).When indicators of hydrologic function are compared across many sites it may be possible to map hydrological function to catchment attributes, which can improve hydrological predictions where historical data sets are short or not available (Wagener et al., 2007).Understanding how hydrological function coevolves with catchment attributes may provide a deeper understanding of these mappings.Dunne (1978) provided a succinct framework for understanding the relationship between climate, landscape morphology, and runoff generation mechanisms.In humid climates with minimal anthropogenic disturbance, thick soils in steep landscapes produce primarily subsurface variable source areas (Hewlett & Hibbert, 1967), with minimal surface expression of changes in aquifer storage.By contrast, in humid environments with shallow soils and more gentle topography, subsurface lateral flow capacity may be exceeded, and surface water runoff may be generated by groundwater exfiltration and precipitation on emergent saturated areas (Dunne & Black, 1970).A recent re-examination of the relationships between topographic properties and runoff generation mechanisms (Wu et al., 2021) further subdivided runoff generation mechanisms and controls, but found that the fundamental relationships identified by Dunne (1978) still emerged.

Runoff Generation and Saturated Areas
However, research so far has not provided sufficient explanation for why certain combinations of topographic properties and runoff generation mechanisms appear (Li et al., 2014).This question is not particularly new.In an early distributed hydrological modeling study, Freeze (1980) noted: The simulations carried out in this study have placed the author in some awe of the delicate hydrologic balance on a hillslope.If one fixes the mean hydraulic conductivity of a hillslope, then there is only a very narrow range of topographic slopes that can lead to runoff generated by the Dunne mechanism.If one fixes the topographic slope of a hillslope, then there is only a very narrow range of hydraulic conductivities that will lead to a water table that is high enough to allow the Dunne mechanism to be operative in a given climatic regime.The fact that the Dunne mechanism is so common in nature in spite of these theoretical limitations on its occurrence infers a very close relationship between climate, hydraulic conductivity, and the development of geomorphic landforms.
What Freeze (1980) observed in simulations indicates that some sort of catchment coevolution (Troch et al., 2015) might be needed to explain a tendency toward saturation excess variable source area runoff generation (the "Dunne mechanism").The literature exploring the evolution of climate-morphology-runoff generation relationships is minimal.Here our goal is to provide a broad picture of what kinds of landscapes and hydrological behavior emerge as we allow climatic, hydrologic, and geologic properties to vary, assuming that runoff generation is driven by saturation excess from shallow groundwater flow and subsurface stormflow.We provide a synthesis of some of our results in the context of relationships identified in field data, including those of Dunne (1978), which indicate that at least certain features of that relationship are emergent products of catchment coevolution.

Dynamic Hydrology in LEMs
Precipitation variability has long been recognized as an important factor in landscape evolution.Ijjász-Vásquez et al. (1992) considered steady-state subsurface flow and an exponential distribution of rainfall depths, and showed that the statistical distribution of resulting erosion rates effectively smoothed hillslope-valley transitions.Tucker and Bras (1998) found similar smoothing of hillslope-valley transitions with a comparable model that used a steady state partitioning of flow between surface and subsurface for storm events drawn from exponential distributions of depth, duration, and interstorm duration.Similar approaches with stochastic precipitation but steady-state hydrological models are still widely in use (e.g., Barnhart et al., 2018).However, these models are limited in that antecedent conditions are not considered; the runoff and sediment transport rate during each event is independent from previous events.
Other studies have taken different approaches to capture some of the effects of memory and event sequence on runoff and erosion.Lague et al. (2005) examined the effects of discharge variability on channel long profile evolution by using a power law distribution of runoff rates, forgoing an explicit model of the processes that convert rainfall to runoff.Deal et al. (2018) further advanced understanding of how runoff distributions affect channel long profiles using the stochastic hydrological model developed by Botter et al. (2007) to generate runoff from a coupled, spatially lumped soil moisture and single reservoir groundwater model.This approach accounts for the important effects of antecedent water storage and evapotranspiration on runoff generation, which ultimately affects fluvial erosion (Rossi et al., 2016).Because of these features, the model developed by Deal et al. (2018) shows promise for understanding how climate translates into erosion events and long-term evolution.However, models of channel long profile evolution cannot quantify spatially distributed hydrological features of interest, such as variably saturated areas.Moreover, it remains unclear exactly how the hydrological parameters (e.g., the reservoir coefficient, or reservoir size) needed for the model developed by Deal et al. (2018) are best linked to the evolving channel profile or surrounding hillslopes.
Although a few previous studies have used LEMs that resolve spatially distributed hydrological features, none have investigated the hydrological function that emerges at geomorphic dynamic equilibrium.Huang and Niemann (2006) used a coupled groundwater-LEM to examine how topographic evolution changed runoff generation at a well studied site, but evolved the landscape for far less time than needed to achieve dynamic equilibrium.Huang and Niemann (2008) investigated long-term evolution with a coupled groundwater-LEM, but examined the sensitivity of modeled topography to hydrologic parameters by prescribing changes onto the slope-area relationship rather than directly simulating the evolution to dynamic equilibrium, which makes evaluating the role of coevolution between runoff generation and topography challenging.Lastly, Zhang et al. (2016) presented a highly detailed coupled hydrological and LEM, but the model has only been used as a proof of concept.

Approach
Here, we focus on the coevolution of topography and runoff generated by groundwater return flow and precipitation on saturated areas.To do this, we use the streampower-diffusion LEM called DupuitLEM that was developed by Litwin et al. (2021), in which runoff produces the shear stress for detachment limited erosion, and topography sets the boundary conditions for the groundwater system.To capture time-varying runoff generation, we include stochastic storm generation and a simplified representation of vadose zone dynamics.We evolve the coupled model toward geomorphic dynamic equilibrium where the denudation rate is approximately equal to the uplift rate.At this point the hydrological function of the landscape is in some sense in equilibrium with topography-what exactly this equilibrium is and how it emerges are central to our results and discussion.
Hydrologic function of emergent landscapes likely depends on geomorphic, hydraulic, and climatic parameters in the model.However, this parameter space is large, and combinations of parameters do not necessarily result in unique model outputs.Dimensional analysis of the model allows us to approach both of these problems.We begin with the nondimensionalization developed by Litwin et al. (2021), and expand it to include the effects of the vadose zone and time-varying climatic forcing.We produce a minimal set of dimensionless groups that both uniquely determine the model output and provide insight into the competing processes that affect evolved morphology and hydrologic function.
Here we focus on how dimensionless groups that express climate and subsurface hydraulics affect (a) topography and drainage dissection, (b) water balance and flux partitioning, (c) spatial patterns of saturation, and (d) temporal relationships between saturated area, discharge, and storage.Based on the results, we present a perceptual model of the emergence and persistence of variable source area hydrology and show that the relationships between geomorphology and runoff generation in humid landscapes that were identified by Dunne (1978) can be obtained through coevolution.Finally we discuss the potential for using landscape history to understand present hydrological function.

Topographic Evolution
The LEM used here considers the evolution of topographic elevation z(x, y, t) by water erosion E f (x, y, t), erosion resulting from the divergence of hillslope regolith transport E h (x, y, t), and uplift or baselevel change U.
Litwin et al. ( 2021) derived the water erosion term from excess shear stress, arriving at a form that is similar to the detachment-limited streampower law, but using the area per contour width a instead of upslope area A. Bonetti et al. (2018) define a(x, y) as the scalar field satisfying ∇ ⋅ (a ∇z |∇z| ) = 1, which is an elegant analytical definition of the concept usually defined as a = A/v 0 in the limit of small contour width v 0 .The erosion law also scales linearly with the dimensionless discharge Q* = Q/(pA), where Q is the volumetric discharge and p is the mean precipitation rate, which we derived from the hydrological model as discussed below.The rate of water erosion is: where K is the streampower incision coefficient, and E 0 is a threshold below which no water erosion occurs.
Although erosion thresholds can have important effects on morphology (e.g., Tucker, 2004), here we only present results for E 0 = 0, as we found the threshold to have little effect on the hydrological behavior of interest in this study.
The term E h describes gravity-driven movement of sediment via processes such as frost heave, animal burrowing, and tree throw.We used linear diffusion law in Litwin et al. (2021), E h ∼ ∇ 2 z.This assumption produces unrealistically steep toe slopes (|∇z| > 1) as hillslopes become long.Landscapes with high relief and long hillslopes generally have a form better described by nonlinear sediment flux laws, where flux increases super-linearly with slope.Near ridges and when relief is low, the law produces near-parabolic topography (like the linear diffusion law), but as the slope gradient increases it produces increasingly planar hillslopes.This replicates a shift from short-range transport to longer-distance transport processes such as dry ravel and shallow mass failures (e.g., Doane et al., 2018;Gabet, 2003;Roering et al., 1999;Tucker & Bradley, 2010).Data compiled by Godard and Tucker (2021) showed that most documented field case studies of hillslope morphology, transport efficiency, and erosion rate fall within the nonlinear transport regime.We chose the hillslope transport model described by Ganti et al. (2012), which is a Taylor expansion of the critical slope model (3) used by Andrews and Bucknam (1987) and Roering et al. (1999), but is more computationally tractable for landscape evolution simulations.
where D is the transport coefficient and S c is a critical slope.This expression represents the first two terms of the Taylor expansion, which Ganti et al. (2012) showed to be a close approximation of the original partial differential equation.
Combining the water-and gravity-driven erosion terms with a constant rate of baselevel change U, we arrive at the governing equations of the LEM: where the angled brackets 〈⋅〉 indicate the time-averaged value of the quantity.At geomorphic equilibrium, the erosion terms must be just sufficient to remove sediment at the rate U everywhere, as is the case in most LEMs.However, our model differs from most in that the fluvial erosion term depends on the presence of surface water that we model explicitly.Surface water is only present where and when the unconfined aquifer interacts with topography, as we will discuss in Section 2.4.

Subsurface Model
As in Litwin et al. (2021), we consider only a homogeneous, surface-parallel layer of permeable material with thickness b above impermeable bedrock.We will refer to this as the "permeable thickness" as we do not distinguish between mobile regolith and weathered or fractured bedrock.Although deeper groundwater flow can be important for runoff generation, here the permeable thickness sets the lower boundary for groundwater circulation.We will sometimes refer to "regolith" when discussing hillslope sediment transport, although ultimately the geomorphic model is agnostic to the composition of the subsurface.That is, we assume there is always enough regolith to meet the slope-based hillslope flux law.Lastly, the term "soil" is used when referencing or drawing comparison with field studies that use this term, in which case the analogous term in our model is permeable thickness.
How exactly the subsurface of real landscapes evolves to keep pace with surface evolution is an active subject of research that so far has no consensus (Riebe et al., 2017).Surface-parallel permeability structure is sometimes (but not always) observed in the field St. Clair et al. (2015), and also emerges at geomorphic steady state with the widely used exponential production model (e.g., Rosenbloom & Anderson, 1994;Tucker & Slingerland, 1997).
Fixing the permeable thickness rather than tracking its evolution does have limitations, as discussed in Section 6.3, but ultimately we decided to keep the subsurface representation simple to focus on the dynamics of topographic and hydrologic evolution.

Hydroclimatological Model
Given the long timescales of landscape evolution relative to runoff generation, it was necessary to make compromises between process representation and computational efficiency.Our goal was to develop a minimally complex model that captured the emergence of catchment and hillslope scale hydrological function (sensu Wagener et al., 2007) including water balance partitioning, and the presence of surface and subsurface variable source areas.We therefore aimed to construct a model that incorporated the following elements: • rainfall, and therefore recharge, varies in time, • rainfall is partitioned between quickflow and storage, • storage is partitioned between ET and baseflow, • ET is limited by energy in humid climates, and by water availability in dry climates.
To address the first element, we generated stochastic storm depth d s , duration t r , and interstorm duration t b using exponential distributions, following Eagleson (1978), and many papers in the hydrology (e.g., Botter et al., 2007; Water Resources Research 10.1029/2023WR034647 LITWIN ET AL.Rodriguez-Iturbe et al., 1999) and landscape evolution literature (e.g., Barnhart et al., 2018;Tucker & Bras, 2000).Previously we introduced the mean precipitation rate, p, which is related to the above parameters as p = 〈d s 〉/(〈t r 〉 + 〈t b 〉).The distributions for storm depth, duration, and interstorm duration are: To address elements 2 and 3, we needed to account for storage in the unsaturated zone as well as the saturated zone.A thorough treatment of coupled saturated-unsaturated zone dynamics would be computationally prohibitive for landscape evolution simulations.We opted instead for a one-way coupling between a simple unsaturated zone model and the Dupuit-Forcheimer groundwater model, which is capable of capturing important features.Schenk (2008) presented a simple model (called here the Schenk model) for vadose zone dynamics in a 1-dimensional profile that serves our purpose well.
The model is based on the assumption that plants extract water from the shallowest depth where water is available, and use all available water at that depth before extracting water from deeper in the profile.Conversely, precipitation fills available storage at the ground surface first, and displaces water already present deeper into the profile.Schenk (2008) showed that the distribution of depths from which water is extracted in this model mimics the plant rooting depth distributions in a wide range of climates.This is useful to our study because the depths of root water uptake emerge as a result of the climate and subsurface hydraulic properties selected rather than requiring an additional parameter.We then calculate recharge at each location on the grid as the amount of water that has infiltrated below the water table.A complete description of the Schenk model can be found in Section A.

Groundwater Flow and Runoff Generation
Runoff is generated by exfiltrating subsurface lateral flow and from precipitation (i.e., recharge from the unsaturated zone model) on saturated areas.We use a quasi 3-dimensional shallow unconfined aquifer model based on the Dupuit-Forcheimer approximations (Childs, 1971) for groundwater flow above a sloping impermeable boundary.We solve for the (saturated) aquifer thickness h(x, y, t) based upon the lateral groundwater flux q(x, y, t), local runoff production q s (x, y, t), and recharge r(x, y, t).Surface water discharge Q(x, y, t) is calculated by instantaneously routing q s over the area upslope from a given location.A conceptual cross section of this model is shown in Figure 1.The governing equations for the hydrological model are: Figure 1.Conceptual hillslope with labeled variables for the groundwater model.The topographic elevation is z, and on hilltops the curvature is approximately h g /ℓ 2 g (see Section 4 for the introduction of these characteristic scales).Permeable thickness b and aquifer thickness h have lower bounds at the aquifer base, which has elevation z b .The angle between z b and the horizontal is θ.The recharge rate is r and the lateral groundwater specific discharge is q.When the water table approaches the surface, local surface runoff q s is produced, which includes both excess of lateral groundwater flow and recharge.The inset plot shows a planview of the square model domain, with three zero flux boundaries and one fixed value boundary.where n e is the drainable porosity, which we assume to have a constant value, k s is the saturated hydraulic conductivity, and θ(x, y, t) is the slope of the aquifer base.The regularization function G( ⋅ ) is equal to zero when the argument is less than 1, and approaches 1 as the argument approaches 1.In this case the argument h/b represents the portion of the total permeable thickness b that is occupied by the aquifer with thickness h.The ramp function R( ⋅ ) is zero when the argument is less than zero and is equal to the argument when it is greater than zero.Thus, Equation 12 says that runoff will occur when the ground is saturated to near the surface and the recharge exceeds the divergence of the groundwater flux.

Modeling Platform
The governing equations were solved on a 125 × 125 square raster grid.The grid cell size is best considered within the framework of the nondimensionalization we use, which will be discussed in Section 4. The top, right, and left boundaries are zero-flux boundaries, while the bottom boundary is a fixed value (Dirichlet) boundary, where the land surface and water table elevation are coincident (see inset plot in Figure 1).The initial condition is a near flat, roughened surface at baselevel.The domain can be considered in a moving reference frame, where the bottom boundary is an adjacent lateral stream (albeit with zero slope) incising at a rate U.The vadose profile was discretized such that each depth increment is equal in size and has a maximum unsaturated storage ≤1% of the mean storm depth.
The model is implemented as a Python package called DupuitLEM that is built on process components from the Earth surface modeling platform Landlab (Barnhart et al., 2020;Hobley et al., 2017).The LEM is solved using existing process components in a loosely coupled scheme, where diffusion is solved with a forward Euler finite volume method and the streampower erosion module is solved with an implicit method based on Braun and Willett (2013).The groundwater model (Litwin et al., 2020) is solved with an approach that combines explicit calculation of lateral groundwater flow and an analytical solution for groundwater rise and exfiltration based on the regularization presented by Marçais et al. (2017).

Upscaling Discharge to Geomorphic Time
A scaling scheme that relates the timescales of geomorphic and hydrologic processes is needed to make the coupled model computationally feasible.There are two primary approaches to this temporal scaling: online updating and offline updating.An online updating approach matches one hydrological time step to one geomorphic time step, and simply scales the effect of each event to represent some longer duration of time.An offline updating approach simulates hydrology on a fixed landscape over some duration from which meaningful average quantities can be derived.The average quantities, including discharge, are then used to evolve the landscape over some longer duration.We used an offline updating approach, as it allows for greater continuity of the hydrological state.We ran the hydrological model for 25 storms before using the results to evolve the landscape forward in time.The scaling factor between hydrologic and geomorphic time varies by the simulation from 250 to 64,000, depending on the duration of the mean storm-interstorm period.Because Equation 5 is linear in the average dimensionless discharge 〈Q*〉, the average is an ordinary time-weighted mean.Further considerations would be needed to account for an incision threshold or different exponents in the incision model.

Scaling Analysis
We used dimensional analysis to identify groups of parameters that affect the solutions to the governing equations in related ways.Litwin et al. (2021) nondimensionalized a simpler version of the model presented here with linear hillslope diffusion and uniform recharge.The nondimensionalization applied the concept of symmetry groups (Barenblatt, 1996), minimal sets of parameters that, when scaled by a constant factor, leave the governing equations unchanged.We nondimensionalized the governing equations systematically by choosing constant factors for each symmetry group and introducing definitions of equivalent dimensionless variables.
Characteristic drainage timescale The approach used characteristic scales derived from the model parameters to isolate dimensions of the model.These are listed in Table 1.Applying the symmetry group approach to the continuum equations of the model used here, we identified the following dimensionless governing equations.
where prime indicates a dimensionless equivalent variable (defined in Section B).Five dimensionless groups, plus the critical gradient S c , remain as parameters.These are listed in Table 2.
Of the dimensionless groups, we expect β and γ to be the most important controls on emergent hydrological behavior, as they affect critical aspects of hydrological function.The aquifer relief index β describes the geomorphic height scale relative to aquifer thickness, which was called the hillslope number in Litwin et al. (2021) because its form is analogous to Brutsaert (2005, their Equation 10.139) used to understand shallow groundwater dynamics.However, we found that it was harder to interpret β as a hillslope number in this study due to the combination of evolving landscape form, time-variable recharge and evapotranspiration.We will return to the discussion of the hillslope number and the role of β in Section 6.5.The drainage capacity γ describes the permeable thickness relative to characteristic aquifer thickness, or equivalently the ratio of a characteristic Darcy flux to precipitation on a hillslope with length ℓ g .The characteristic gradient α is the ratio of geomorphic height and length scales, which we will keep fixed in this paper, but was explored in Litwin et al. (2021).The timescale factor δ is the ratio of hydrologic to geomorphic timescales, which we expect to be small in all cases given the large difference between hydrologic and geomorphic process rates.Lastly, λ is the domain scale factor, where l is the domain side length.λ is large in all cases considered here, and is not expected to affect our results (Anand et al., 2022;Bonetti et al., 2020;Litwin et al., 2022).
By nondimensionalizing the landscape evolution and hydrological equations together to arrive at these groups, we reveal how hydrological function is a product of hydrological and geomorphic parameters together.For instance, U appears prominently in many of the dimensionless groups, honoring the wellestablished importance of uplift rate for the development of relief.Relief is in turn an important control on the development of aquifers.As such, the dimensionless groups allow us to effectively explore emergent hydrological behavior in the parameter space without neglecting the importance of the geologic and geomorphic setting.
The stochastic forcing introduced three additional parameters: mean storm duration 〈t r 〉, mean interstorm duration 〈t b 〉, and mean storm depth 〈d s 〉, while adding vadose zone dynamics and evapotranspiration introduced two additional parameters: evapotranspiration rate pet and plant-available water content n a .We found that four dimensionless groups were needed to represent the additional parameters, shown in Table 3.These groups were chosen to provide intuition into competing processes.The dimensionless forms of the Schenk model are given in Appendix C, Equations C8 and C15.
Additional dimensionless parameters needed to characterize the vadose model and stochastic climate are shown in Table 3.The event storage index σ describes the competition between the maximum possible saturated zone storage and mean storm depth, where smaller values indicate that local saturated zone storage is more easily exceeded by storms.Ai is the duration-corrected rate of potential evapotranspiration relative to rainfall, which has a critical effect on how much water becomes recharge.The precipitation steadiness index ρ is the proportion of time in which rainfall is occurring, which in the limit of ρ → 1 is the steady case considered by Litwin et al. (2021).Lastly, the moisture content index ϕ describes the vadose zone plant-available water content relative to the saturated zone drainable porosity.

Results
We conducted simulations to explore the effects of subsurface properties and climate on morphology and runoff generation by varying the dimensionless groups identified in Section 4.Although the simulations do not cover the entire parameter space, they are sufficient to show a range of hydrologic behaviors that emerge from the coupled model.The results are presented in three sets of simulations.First, we varied σ and Ai, holding other dimensionless groups constant (except the timescale factor δ, which we assume to have negligible effect).Second, we ran the same combination of σ and Ai but decreased β by a factor of 10 to examine a case where aquifers are thicker relative to relief.Finally, we fixed Ai to a humid value of 0.5 and varied γ and σ to explore the interaction of storage and drainage in the subsurface.
All our results focus on the climatic end-member where storm durations are short relative to the time between storms (ρ = 0.03).We used the average value of α investigated by Litwin et al. (2021) (α = 0.15), the domain scale factor was fixed at λ = 250, and others were chosen to be physically reasonable, including a critical slope S c = 0.5 and moisture content index ϕ = 1.5.Hydrological fluxes are averaged over 2,000(〈tr〉 + 〈t b 〉).Further details on the hydrological analysis can be found in Appendix D.

Effects of Climate on Topography
Before examining the hydrological function of the evolved landscapes, we will begin by examining their topography.Aridity index Ai and event storage index σ play important roles in determining the partitioning of precipitation into evapotranspiration, surface flow, and subsurface flow, which ultimately affects the amount of water available to shape topography through erosion.The hillshades in Figure 2a show the development of characteristic ridge-valley topography when Ai < 1, where drainage dissection decreases with increasing aridity.When Ai ≥ 1 drainage networks are minimal or nonexistent (given the domain size, boundary conditions and other parameters used).Relief also increases with increasing aridity, as relief increases with decreasing drainage dissection while α and S c are held constant.
The event storage index σ modulates the relationship between topography and aridity.When σ is small, the subsurface has a small capacity to store water relative to the average storm depth, and consequently surface runoff is produced more frequently across more of the landscape, increasing dissection and lowering relief.In the results  presented, we decreased σ by reducing rainfall frequency and increasing intensity.Figure 2 shows that drainage networks can form under higher aridity climates when σ is small, as large infrequent storm events have more potential to generate surface water runoff than if the same annual precipitation were spread amongst more frequent storms.
Cross sections through the subsurface (Figures 2b-2e) show differences in relief and mean water table position between selected model runs.Humid landscapes with small σ have the least relief and maintain water tables close to the surface (Figure 2d).Arid landscapes with small σ have the highest relief, and the aquifer is absent except near stream channels (Figure 2e).When σ is large (small storms relative to storage), aridity has a highly non-linear effect on topography, in which effectively no stream channels emerge in cases where Ai > 1.

Water Balance Partitioning
We examined the water balance at two levels, first partitioning precipitation into actual evapotranspiration (AET ) and total runoff and interpreting the results using the Budyko framework (Budyko, 1974), and then partitioning total runoff into baseflow and quickflow and interpreting the results using the L'vovich framework (L'vovich, 1979).This approach was applied to all three sets of model runs (varying Ai and σ, varying Ai and σ with low β, and varying σ and γ), which have corresponding topographies in Figures 2, Figures S5, and S6 in Supporting Information S1, respectively.
The Budyko plots in Figures 3a and 3d, show how precipitation is partitioned to evapotranspiration (rather than discharge) as a function of aridity, and the constraints that energy and mass balance place on this partitioning (dashed lines) for high and low β cases.In energy-limited environments, the maximum ratio 〈AET〉/〈P〉 is 〈PET〉/ 〈P〉 ≈ Ai, whereas in water-limited environments, the maximum ratio is one.Model results in Figure 3a closely follow respective energy and water limitations at each aridity, indicating actual ET is occurring at close to the potential ET rate.In contrast, Figure 3d shows that when the aquifer relief index β and event storage index σ are small (i.e., thinner aquifers relative to relief and smaller storage capacity relative to storm depth) and the climate is humid, substantially less precipitation becomes evapotranspiration (and more becomes discharge) than in the previous case.Figure 3g shows how this partitioning is affected by drainage capacity γ for a constant aridity Ai = 0.5.Here poorly drained landscapes (low γ) appear to produce less actual ET relative to precipitation, although the effect is smaller than that of Ai.In the particular stochastic simulations we ran, 〈PET〉/〈P〉 > Ai, so we place the horizontal line at 〈PET〉/〈P〉 to show that actual ET still does not exceed potential ET.
The L'Vovich framework allows us to more deeply understand the catchment water balance by decomposing discharge into quickflow Q f that leaves the watershed rapidly during storms, and baseflow Q b that is released more slowly.We examine (a) how precipitation is partitioned into quickflow and storage, and then (b) how storage is partitioned into ET and baseflow.
We first consider the fraction of precipitation that becomes quickflow, shown in Figures 3b, 3e, and 3h.These show that quickflow fraction is sensitive to all dimensionless groups considered (γ, β, σ, and Ai).In Figure 3b, the quickflow fraction decreases rapidly with increasing aridity, until almost no quickflow is generated when Ai ≥ 1.In contrast, Figure 3e shows that when the aquifer relief index β is small, quickflow is more sensitive to event storage index σ, and for σ = 8 the quickflow is greater that 50% even when Ai = 1.Quickflow fraction declines nonlinearly with increasing drainage capacity γ (Figure 3h), similar to the effect of aridity shown in Figure 3b.
Second, we consider how the remaining precipitation (that has become storage rather than quickflow) is partitioned into evapotranspiration and baseflow (Figures 3c,3f,and 3i).In Figure 3c, we see that the baseflow fraction declines linearly with aridity in humid climates, and is small when Ai > 1.This is also true when the aquifer relief index β is small (Figure 3f), provided the event storage index σ is large.However, when σ is small and the climate is humid, partitioning to baseflow is less sensitive to aridity, similar to the sensitivity seen in Figures 3d and 3e.Although quickflow fraction decreases with drainage capacity γ and aridity (Figures 3b and  3h), baseflow fraction generally increases with γ (Figure 3i), the opposite of the trend between baseflow fraction and aridity (Figure 3c).The baseflow fraction increases with γ until it levels out at a constant value as actual ET approaches potential ET (Figure 3g).

Spatial Structure of Saturated Areas
The location and extent of saturated areas vary in response to changing recharge, water storage, and topographic states.To examine spatial and temporal patterns of saturated area, we defined a metric of saturation occurrence and classified the landscape into zones that are wet, variably saturated, or dry.We classify surface saturation as where the water table is within 0.025h g of the ground surface.This metric approximates the "squishy boots" test used to identify variable source areas (e.g., Dunne et al., 1975).Areas that are saturated at the end of more than 95% of storms and interstorms are classified as wet, whereas locations that are saturated after less than 5% of storms and interstorms are classified as dry.Variably saturated areas are all others not in either of the previous classes.For our purposes, the classification is relatively insensitive to the choice of threshold values; for details, see Figure S4 in Supporting Information S1.
Figure 4 shows that the model produces widespread variably saturated areas organized around the interface between the channel network and adjacent hillslopes.In humid landscapes where the event storage index σ is small, channel networks are permanently saturated, and hillslopes can occasionally become fully saturated (frequency >0.05), when storm depths approach or exceed local saturated zone storage capacity.With increasing σ, variable source areas retreat to localized areas of topographic convergence.With increasing aridity, the water table tends to interact with the surface less frequently, leading to intermittent channel networks when Ai ≥ 1.
The transition to intermittent saturation in valley bottoms is also affected by the drainage capacity γ, due to its influence on the partitioning of water between surface and subsurface flow (Figure S10 in Supporting Information S1).When γ is large and the event storage index σ is small, storms are large relative to storage, but subsurface drainage is efficient.Consequently, ridges remain dry, but saturation in valley bottoms is more variable than cases with lower γ or higher σ (see Figure S10 in Supporting Information S1, simulation 4).
We found discontinuous wet sections of the channel network emerge from the model without any introduced heterogeneity or spatial variation in permeable thickness.This appears in Figure 4 subplots 3, 4, 10, and 11, which have high relief relative to aquifer thickness (large β), intermediate aridity Ai, and large storms relative to storage capacity (small σ).They also appear when drainage capacity γ is large (Figure S10 in Supporting Information S1), but are largely absent when β is smaller (Figure 5).These patterns are driven by differences in the local convergence and downslope conveyance capacity associated with topographic curvature and slope, and could be indicative of discontinuous stream channels with perennial and ephemeral reaches.However, our saturation  2, where the aquifer relief index β = 0.5).Surface saturation is determined on the basis of timevariable water table proximity to the surface.Locations are classified as dry if they experience surface saturation at <5% of the ends of storms and interstorms, and are classified as wet if they are saturated at the end of >95% of storms and interstorms.Variably saturated areas are everywhere that does not meet either of these criteria.The results show the extent of variably saturated areas is greatest when σ is small.Non-permanent streams emerge in some cases as aridity increases, including simulations 4 and 5, in which the channel network contains discontinuous reaches that are always wet.
metric describes only the proximity of the water table to the surface, and does not include the presence of water routed from upslope, which in our model cannot re-infiltrate.
Other unusual features emerge as the water table relief approaches topographic relief (when β is small).First, we observed that particular combinations of parameters produce trellis-like drainage networks that are nearly divergent (Figure 5).This is highly unusual in LEMs, particularly those with single-direction flow routing and fluvial incision like this one.We say "close to" divergent because subtle drainage divides do exist such that surface flow is only routed in one direction at any particular topographic state-that is, saddle-points.However, saturation can extend up to these saddle-point divides, and flow directions near them may change frequently with evolving topography.
Second, some model results show the presence of persistently saturated valley bottoms with widths greater than one pixel (e.g., Figure 5 subplots 10, 11, 17, 23, and 28).This is also uncharacteristic of the type of LEM formulation used here, which will generally incise valleys only one grid cell wide (Tucker & Hancock, 2010).This illustrates that erosion by runoff on saturated areas near the toes of hillslopes can help account for the formation of valleys that are substantially wider than the channels they contain.Landscape evolution models have generally only achieved valleys wider than one pixel by explicitly representing valley widening by lateral channel migration (Langston & Tucker, 2018).These permanent lowland wetland features and trellis-like drainage networks are evidence of the strong influence that the aquifer structure can exert on surface drainage organization, even in a relatively simple model.

Co-Variant Dynamics of Spatially-Averaged Saturation, Storage, and Discharge
The relationship between saturated area and baseflow discharge is a useful indicator of the relationship between landscape morphology, subsurface properties, and runoff generation (Latron & Gallart, 2007).We examined baseflow rather than total flow because the total flow has greater dependence on storm intensity, whereas baseflow is primarily exfiltrating groundwater, which should vary more systematically with aquifer properties and topography.The classes are the same as in Figure 4.In contrast to the higher β case, here we see variably saturated areas from valleys to ridges in much of this parameter space.In the transition zone between widespread variable saturation and the zone without any channels, we see unusual channel forms, including trellis-like drainages (10,16,22,23,28), and extensive valley bottom wetland zones (11,17,23,24,28).
Figure 6 shows the dimensionless baseflow discharge where Q b is the total baseflow discharge for the model domain, A sat is the total saturated area, and A tot is the total domain area.Saturated areas are calculated with the same criterion as in the spatially distributed figures.For reference, light gray points were added to indicate the total dimensionless discharge . Baseflow points are colored by the dimensionless saturated storage S* = S/(bn e A tot ).
As expected from the spatial patterns of saturation in Figure 4, the range of the dimensionless saturated area decreases with increasing storage capacity relative to storm depth (σ).When σ is large, approximately 10% of the watershed area is saturated in the most humid case, which decreases with increasing aridity.When σ is small, increasing aridity does not prevent the landscape from reaching near full saturation (A * sat = 1) , but does lower the minimum saturated area, increasing the observed range of A * sat .Model runs with the same aridity tend to have similar minima of A * sat , but with decreasing σ, the maxima generally increase.
Despite differences in the saturation-baseflow discharge relationship with the parameters shown, there are underlying patterns that may reveal coevolution.Primarily, we notice that the relationship between baseflow and saturated area has a concave up form in most cases, where the rate of change of saturated area increases with baseflow discharge.Furthermore, the simulations with the largest range in (log-transformed) A * sat (e.g., in Figure 4) appear to have a sigmoidal relationship, which can be divided into three regimes: rapid increase in saturated area with low baseflows, moderate increases in saturated area with moderate baseflows, and again rapid increases in saturated area with the highest baseflows.Several reference lines are included for comparison with these regimes: , which are indicative of the rates of change in the upper and middle regimes, respectively.sat space back to their underlying saturation patterns.We chose to examine the results in Figure 6 subplot 4 (σ = 8.0, Ai = 1.0) in more detail, as it displays the three-regime form well. Figure 7 shows this mapping, where subplots B-E are hillshades colored blue where the ground is saturated.Corresponding points are labeled in subplot A. Subplot (B) shows saturated areas cover only the second order and higher stream channels under low-flow conditions.In (C), saturated areas extend up through first order channels, and some channel-adjacent areas.In (D), saturated areas emerge in many unchannelized concave regions, while by (E), saturation is widespread on all concave and planar regions, extending toward ridges.Critically, we can see that the inflection point near (C) represents the threshold above which saturated areas emerge outside of the channel network.
The geomorphic transition between in-channel and out-of-channel saturated areas at the transition between the middle and high baseflow discharge regimes translates to cases that do not display all three regimes.Figure S13 in Supporting Information S1 shows how saturation patterns are related to points in the baseflow saturated area relationship for the case shown in Figure 6 subplot 14 (σ = 8.0, Ai = 0.25).The point of maximum curvature is still associated with increasingly widespread saturation outside of the incised channel network.This supports the idea that the saturation-baseflow relationship embeds information about landscape morphology (at least hillslopechannel transitions).However, as the cases with large event storage index σ demonstrate (Figure 6 subplots 28-31), the extent and variability of saturated areas affect how much of the morphology is visible in this relationship.Varying the drainage capacity γ affects the shape of the relationship between baseflow discharge and saturated area.The slope of the middle regime decreases with increasing γ, and the transition between the middle and upper regimes sharpens.Topographically, high γ cases also have greater relief, lower drainage density, and sharper transitions between channels and hillslopes (Litwin et al., 2021).In contrast, when the aquifer relief index β is small, the relationship between saturated area and baseflow discharge weakens, as shown in Figure S11 in Supporting Information S1.We will return to additional synthesis of these relationships in the discussion.

How Do Topography and Hydrology Coevolve in DupuitLEM?
The purpose of the model developed in this paper is to help us better understand how real topography and hydrologic dynamics coevolve.A clear conceptual understanding would make it far easier to comprehend the sensitivity of the results to variations in the parameters presented above, and to ascertain where the model may provide insight and where it is deceptive.A visualization of our conceptual understanding is illustrated in Figure 8.The "perennial aquifer" corresponds to 95% exceedance probability for aquifer thickness, and the "variable water table" corresponds to 5% exceedance probability (Lower panel) selected model results, where key dimensionless parameters have been varied from a base case, which is shown in the middle.The bold headings should be read as the difference between the model results at the tip and tail of the gray arrows.For example, the upper left result is the same as the base case, except for increased transmissivity (associated with increasing γ).Likewise, the lower right panel is the same as the upper right panel, except for increased storm frequency and decreased storm size (increased σ).

Balance Between Water and Sediment Supply
In DupuitLEM, hillslope morphology evolves to simultaneously shed water and regolith at the rates they are supplied.Water is supplied by rainfall and is lost to runoff (when it falls on saturated ground), subsurface drainage, and ET.Regolith is supplied by uplift or baselevel change U, may be redistributed by diffusive hillslope transport E h , and removed by water erosion E f .Note that diffusive hillslope transport (unlike water erosion) tends to spread regolith out, rather than removing it from the model domain (except at the boundaries).Therefore, the simulated landscape morphology must therefore evolve toward a condition where the production of runoff is just sufficient to remove regolith by water erosion at (areal-averaged) rate U.
The key to achieving this balance is the perennial aquifer that forms in areas of topographic convergence.For the visualizations in Figure 8 the perennial aquifer has been defined based on 95% exceedance probability of aquifer thickness, and therefore corresponds to the "wet" saturation zones in Figures 4 and 5, Figure S9 in Supporting Information S1.The perennial aquifer appears as dark blue in Figure 8.When storms are large and infrequent (i.e., small σ) and transmissivity is moderate (i.e., γ > 1), the perennial water table determines which areas experience variable and perennial surface saturation and fluvial erosion.Perennial saturation occurs where the aquifer reaches the surface.Variable saturation occurs where the perennial water table is shallow enough that storms can raise the water table to the surface (the variable water table shown in Figure 8 is defined based on 5% exceedance probability of aquifer thickness).Therefore, the landscape morphology must evolve such that the spatial extent of the perennial aquifer is large enough to ensure sufficient surface runoff production.This is accomplished by balancing supply and demand.The aquifers are continually draining, and so to remain at or above their minimum level, the aquifers must receive a continual supply of lateral subsurface inflows from adjacent hillslopes.Both the supply rate and drainage rate are controlled by the topography.Therefore, by controlling these rates (and therefore the extent of the perennial aquifer) a topography can emerge that ensures sufficient runoff production to remove regolith at the supplied rate U.
The rate lateral flow is supplied to the perennial aquifer is determined by the size of the accumulated area upgradient, and the recharge rate in that area.That recharge is exported downgradient toward convergent areas as subsurface flow.When the drainage capacity γ is sufficiently large, the subsurface flow is sufficiently efficient that uplands never experience surface saturation (these are the beige areas in Figure 8).Consequently, upland regolith must be exported to convergent areas via diffusive hillslope transport.
The rates of regolith and water export from the uplands to the convergent areas must strike a delicate balance.The regolith export must be small enough that it does not overwhelm the capacity of water erosion in the convergent areas to remove it.The water export must be large enough that it can sustain the perennial aquifer that makes that water erosion possible.However both the regolith export and water export rates will depend on the accumulated area at the transition from uplands to convergent areas.Therefore, the drainage density must adjust until these demands are in balance.If the drainage density is too small, excess lateral flow from the uplands will expand the perennial aquifer, leading to increased surface saturation and water erosion.If the drainage density is too large, lateral flow will be insufficient to maintain the perennial aquifer and promote water erosion, and so diffusive regolith flux will gradually fill the convergent areas.

The Importance of Transmissivity
The perennial aquifer is maintained when the rate of lateral inflow equals the downslope drainage rate.The drainage rate is controlled by the transmissivity k s b (which is fixed), the local slope of the aquifer base (which in DupuitLEM is determined by the topography because the permeable thickness b is constant in space), and potentially also by the level of the water table farther downgradient.The latter is only important when the aquifer is thick relative to topographic relief (i.e., β is small).
The roles of the recharge, transmissivity, and local slope at the transition from uplands to convergent areas are captured by the dimensionless parameter γ.The slope of the convex uplands will tend to vary downslope in proportion with distance from the ridge and with ridge curvature.More precisely, at distance x we would expect the slope to be approximately xξ(x)h g /ℓ 2 g .ξ(x) is less than 1 and captures the effect of the nonlinearity in the hillslope diffusion law (in fact ξ(x) = tanh xh g /ℓ 2 g / S c )/ xh g /ℓ 2 g / S c ) for the exact form of the nonlinear diffusion law (Equation 3).The maximum subsurface flow per unit width at that point is therefore xk s bh g ξ( x contour width upslope from that point is a(x) and the recharge is r(x), it follows from the definition of γ that in the vicinity of the transition from uplands to where surface saturation and water erosion becomes important the following is true: The first term on the left is the area per contour width divided by distance from the ridge.This quantity is a measure of the degree of topographic convergence.It will be ≈1 for straight slopes, <1 for divergent areas, and >1 for convergent areas.The second term measures the fraction of precipitation that becomes recharge, and is therefore influenced by the aridity Ai.Therefore, γ sets the degree of upland contributing area convergence needed to produce surface saturation and water erosion, modulated by the effect of water balance on recharge and nonlinear slope processes.This makes it clear why γ is an important control on the drainage density of the coevolved landscapes (see Figure S9 in Supporting Information S1), and why aridity also plays a role (see Figure 4).Both effects are illustrated in Figure 8.

The Effect of Decoupling Topographic and Water Table Gradients
Note that the drainage capacity γ depends on the transmissivity k s b, rather than on the hydraulic conductivity k s alone.That means it is possible to vary the permeable thickness while keeping γ constant by also varying the hydraulic conductivity inversely.Doing so amounts to varying the relief index β-a small β corresponds with a large permeable thickness.This was explored in Litwin et al. (2021), where β was referred to as the hillslope number Hi.This is perhaps regrettable, because although β is closely related to the hillslope number (as we shall see in Section 6.5), they are not the same thing.As with the hillslope number, when β is small, the aquifer thickness becomes large relative to the relief, making it possible for water table gradients to substantially differ from topographic gradients.
As a consequence, when the aquifer relief index β is small, the drainage rate of the perennial aquifer is more dependent on the landscape morphology downgradient.Because the slopes downgradient are gentler in lowland areas, the drainage rate is slower relative to the case with large β.Consequently the rate of lateral inflows can be smaller, and a smaller upslope area is needed to supply those inflows.This explains the larger areas of perennial and variable surface saturation when β is small (Figure 5), compared to when it is larger (Figure 4).Some remarkable effects emerge when β is small, as shown in Figure 5.Under the right circumstances, the large permeable thickness allows the perennial aquifer to be connected across surface topographic divides, resulting in connected loops surrounding isolated "islands" of uplands.Broad, low-gradient areas of perennial and variable saturation emerge (Figure 8), particularly around confluences.The pattern of isolated hills above irregular, trellislike drainage networks is reminiscent of cockpit or cone karst landscapes found in many parts of the world with humid climates and limestone bedrock (Lyew-Ayee et al., 2007;Waltham, 2008).While our model does not capture chemical weathering that is an important component of karst landscape evolution, our results suggest that deep, permeable aquifers in humid climates may be primed for such topography regardless.

The Role of Climate
The conceptual explanation of the model results presented here can also help us understand the variations with event storage and aridity (σ and Ai) that appear in earlier figures.When σ is small, storms are large and infrequent relative to the subsurface storage capacity.Longer times between storms allow for greater drainage of the aquifer, and larger storm depths promote greater water table rise during storm events, both of which increase the variably saturated area (VSA) (Figures 4 and 5, and S8 in Supporting Information S1).The large extent of areas contributing runoff when σ is small allows for widespread water erosion to remove regolith from the landscape.This increases drainage density and lowers hillslope relief (Figure 2) in order to maintain the balance between uplift and hillslope denudation.When σ is large, the perennial aquifer is more extensive because there is not time to contract before the next event, but the smaller storms mean transient surface saturation is less likely.This has the opposite effect on drainage dissection and hillslope relief.
The vadose zone is an important link between climate aridity and topography.In our results, there is a transition between channelized and unchannelized topography between Ai = 0.71 and Ai = 1.41 in Figure 2, and the transition is more abrupt for large σ.This can be attributed to the reduced likelihood of recharge when Ai > 1 as captured in the Schenk model of the vadose zone.Recall that ET creates a storage deficit in the vadose zone that must be satisfied before rainfall from an individual storm event can produce recharge at depth.When Ai < 1, the potential ET between storms tends to be smaller than the rainfall that typically falls in each storm.Consequently the vadose zone tends to be wet at depth, and the effects of ET are limited to generating deficits close to the surface.However, when Ai > 1, the storage deficit that can accrue between storms is larger than the depth of rain that typically falls in each storm.Consequently the vadose zone is dry at depth, and recharge will only occur from a storm large enough to fill the profile, or when storms are clustered together.This becomes increasingly unlikely when σ is large.Given less frequent recharge when Ai > 1, larger hillslopes are needed to supply the lateral flow necessary to sustain a perennial aquifer in areas of topographic convergence.This makes areas of permanent or variable saturation less extensive.

How Realistic Is the Hydrology in DupuitLEM?
As illustrated in the results above, DupuitLEM produces landscapes that not only have the appearance of realistic topography, they function hydrologically as one would expect of a realistic landscape-up to a point.The model results deviate from what we might expect in a real landscape in several ways that are worth highlighting.
In the arid cases, as shown in Figure 3a, almost no runoff is produced by the model, and the resulting landscape is unchannelized.Many real arid landscapes will still produce substantial runoff at an aridity index of 2 (the maximum value we considered) (e.g., Wang & Wu, 2013), and exhibit widespread channelization.However, runoff in these settings is more often produced by infiltration excess rather than interaction of the water table with the ground surface (Wu et al., 2021).This mechanism is not included in the analysis here, but could be easily added to DupuitLEM.
We also observed unexpected deviations from energy and water limitations in the Budyko plot when the aquifer is thick relative to topographic relief (small β) and storms are large relative to storage (small σ), as seen in Figure 3d.Because quickflow is primarily derived from precipitation on saturated areas, the deviations suggests these cases have more extensive saturated areas than their high β counterparts in Figure 3b.The contrast in perennially saturated areas between high and low β cases (Figures 4 and 5) indicates this is the case.A larger fraction of quickflow is more likely when σ is small because in these cases storms more easily overwhelm available storage and produce additional saturated areas.
The other large deviation of 〈AET〉/〈P〉 from 〈PET〉/〈P〉 occurs when the drainage capacity γ is small.This is somewhat paradoxical; Troch et al. (2013) found that landscapes with the longest drainage timescales tend to have the highest ET relative to precipitation because water that stays in the landscape longer is more likely to become ET.This is partially the case in our results.Figure 3i shows that the baseflow fraction of remaining water is smallest when γ is small, indicating that more water is becoming ET.However, this is not controlling the overall water balance behavior.Instead, quickflow behavior controls the decrease in 〈AET〉/〈P〉 with decreasing γ, as the proportion of precipitation that becomes quickflow declines precipitously with increasing γ, as shown in Figure 3h.Increasing this quickflow fraction decreases the water that remains available to become ET.Poorly drained landscapes also would be expected to make more water available for ET, as Troch et al. (2013) showed, but because our model ET cannot access water in the saturated zone, we are not able to reproduce this observation.
The one-way coupling of saturated and unsaturated flow also has implications for runoff generation.For example, with little or no additional recharge, the water table can rise rapidly into the capillary fringe (e.g., Crosbie et al., 2005;Gillham, 1984;Weeks, 2002).If the capillary fringe extends to the surface, as it may in wetter areas like concave hillslopes and valley bottoms, saturated areas could expand rapidly during storm events.Saturation of the soil profile due to wetting front propagation (e.g., Ogden et al., 2017) also could enhance the rapid emergence of saturated areas.On the other hand, ET from the saturated zone where it is near the surface could substantially reduce saturated areas during interstorm periods.Because we do not capture these features, we may substantially underestimate the variability of saturated areas and, depending on their relative importance, we may overestimate or underestimate runoff generation from saturation excess.
We also considered only a one-way coupling of groundwater and surface water, such that reinfiltration was not possible.This limits our ability to interpret discontinuous saturation patterns, as seen in Figure 4. Nevertheless, the emergence of this discontinuous network of saturated areas indicates that the morphology of the landscape, rather than just variability in subsurface properties, may provide a structural control on heterogeneous patterns of surface flow in valley bottoms.This feature is likely to persist in a model that allows re-infiltration, although instead of variably saturated valley bottoms, some areas of the parameter space may instead produce reaches that gain and lose water, again as a function of adjacent landscape morphology.
The form of our groundwater model may also affect features that we observe across our parameter space.In order to have tractable solutions for the LEM, the groundwater flow model we use relies on the Dupuit-Forcheimer approximations, which are valid where the component of flow normal to an impermeable lower boundary is small.This usually occurs when saturated thickness is small relative to hillslope or seepage face length (Bresciani et al., 2014), which may not be valid everywhere in our model parameter space.Field studies have also shown that deeper flow paths not captured by our model are important components of stream runoff, especially during baseflow conditions.Accounting for these deeper flow paths could increase baseflow discharge, changing the L'vovich water balance partitioning shown in Figure 3.

What Processes Were Left Out of DupuitLEM?
In the interest of creating a tractable model, we have left out several key climatic, hydrologic, and geomorphic processes that may affect the coevolution of runoff generation and topography.In addition to the limitations discussed in the previous section our representation of climate is simplified, as we neglected seasonality of precipitation or potential ET, which are important controls on the water balance and on the extent of saturated areas (Latron & Gallart, 2007;Yokoo et al., 2008).Considering only single direction flow routing with no depression storage also limits the development of valley bottoms and wetlands that can be important zones for saturation excess overland flow.
The style of water erosion is also limited in this model, as we consider only detachment-limited fluvial erosion, neglecting fluvial sediment deposition and factors such as groundwater sapping (Abrams et al., 2009;Laity & Malin, 1985) and pore-pressure driven landslides (Montgomery & Dietrich, 1994), which could be the subject of separate studies of coevolution between topography and groundwater systems.
We have also limited our study to understanding the evolution of topography, while the progressive weathering of rock and development of regolith are simultaneous components of critical zone evolution (Anderson et al., 2013(Anderson et al., , 2019;;Brantley, Eissenstat, et al., 2017;Brantley, Lebedeva, et al., 2017;Harman & Cosans, 2019).Furthermore, we considered only cases where the subsurface porosity and hydraulic conductivity are constant in a zone that uniformly parallel to topography, while spatially variable subsurface structure may have important implications for runoff generation.Here we have laid the foundation for future modeling that considers both surface and subsurface features and processes.
Lastly, we exclusively examined results at geomorphic dynamic equilibrium.While this approach is valuable for understanding how hydrologic function is adjusted to morphology, it necessarily limits the landscapes where our results are directly applicable.Further work on the transient evolution of landscapes with hydrology could help clarify the role of time in catchment coevolution (Troch et al., 2015), and could provide useful insights into drainage network initiation (Cullen et al., 2022).

Does DupuitLEM Match Field Evidence for the Relationship Between Saturated Area and Baseflow?
As we have shown, saturated area-discharge relationships contain information about runoff generation.However, variation in saturated area through time has not been widely reported in the literature as the measurements are labor intensive and can be sensitive to the judgment of the observer.Latron and Gallart (2007) compiled many published relationships into a single plot (reproduced in Figure 9) that shows a range of forms the relationship can take.We compared our results with those in this plot by choosing a set of dimensioned parameters to re-project some of the dimensionless results in Figure 6 into the dimensioned world.The caveat to this approach is that the position of the results, especially along the x-axis, is subject to the particular dimensioned parameters we chose.
In humid climates, our results show strong resemblance to the concave up form observed by Dunne (1978) at Sleepers River, VT, USA.Saturated area and baseflow in this relationship were measured in Sleepers River Watershed W-2, which has gentle topography and relatively low permeability soils (Dunne et al., 1975), consistent with our low σ cases.Dunne et al. (1975)  not shown), consistent with our high σ cases.Field relationships by Ambroise (1986), Latron (1990), andMyrabø (1986) shown in Figure 9 have convex forms, where baseflow increases faster than saturated area in logspace.Some of our model results (e.g., Figure 7) also have convex forms for lower baseflow and saturated areas, but these relationships seem to have a different origin.In our case, the low baseflow regime was associated with channel network ephemerality, whereas the field studies are still primarily describing variable source areas in valley bottoms and adjacent hillslopes.In fact, the studies here with the lowest saturated extents appear to have linear or slightly concave relationships.However, the linear relationship shown by Latron and Gallart (2007) is from a terraced landscape with fragmented saturated areas, which obscure the link between topography and baseflow.
Reasons why observed relationships could be different from our model predictions are numerous.Our model has only considered a limited number of runoff generation and landscape evolution processes, and lacks the heterogeneities and complexities of real watersheds.However, it is encouraging that our results agree with field surveys from Dunne et al. (1975), given that ours are emergent features of coevolution.

How Does DupuitLEM
Compare to the Dunne Diagram?Dunne (1978) presented a synthesis of how runoff generation mechanisms are related to topography, subsurface properties, and climate, often called the "Dunne Diagram" (Figure 10a).On the humid half of the diagram, Dunne associated saturation excess overland flow (i.e., Dunne overland flow) with gentle topography, and moderate to poorly drained soils.Subsurface storm flow was considered the opposite end member, and was associated with deeper, more permeable soils and steeper straight to convex topography.
We mapped this concept into a quantitative relationship between hydrological and geomorphic metrics for comparison with our results.We selected topographic variance Z as the geomorphic metric to differentiate between gentle and steep topography in terms of hillslope relief.We calculated the variance of topographic elevation for each horizontal slice of the domain (parallel to the open boundary), found the mean of the slice variances, and took the square root to obtain Z with units of length.We normalized Z with the characteristic height scale h g for consistency with our dimensionless framework.For the hydrological metric, we selected the fraction of quickflow relative to total flow 〈Q f 〉/〈Q〉, because in our model, quickflow is generated primarily by Dunne overland flow.
We found a clear mapping between topographic variance and quickflow fraction (Figure 10b), in which model runs that have gentle topography (low Z) generate the most runoff via Dunne overland flow (high 〈Q f 〉/〈Q〉).
Figure 10c shows the same results, but colored to show the parameters.All cases are humid (Ai = 0.5), but have a range of values for the other parameters discussed in this paper.The most consistent pattern is that cases with low drainage capacity γ tend to produce more quickflow.The storage index σ is a secondary control, with generally gentle topography and more Dunne overland flow produced when σ is small, provided that γ is not too large.Low aquifer relief index β cases are the most likely to break expectations, which is expected given their tendency to evolve unusual drainage networks and wide valley bottoms.
In the Dunne Diagram framework, Dunne overland flow is associated with moderate to poorly drained soils and low relief, gentle topography.In our model, poorly drained soils (relative to climate forcing) produce the associated gentle topography through water erosion to maintain the balance discussed in Section 6.1.Likewise, well-drained soils produce steeper (higher variance) topography by expanding the zone where overland flow and water erosion do not occur.Consequently, they develop hydrological responses dominated by subsurface flow (baseflow).

What Controls Variably Saturated Extent? Role of the Hillslope Number
The results in Section 6.5 convey some information about the VSA (shown in colors); however, close inspection of Figure 10b shows that the relationship between quickflow fraction and VSA is not monotonically increasing.
Figure 9. relationship between saturated area and baseflow discharge for a several well-studied sites, reproduced from Latron and Gallart (2007), along with several of our model runs, 0, 4, 16, and 28, with different aridity index Ai and event storage index σ, from Figure 6.The simulation results have been re-dimensionalized for the sake of obtaining baseflow discharge in the appropriate units.Our results are clipped to the extent originally presented in Latron and Gallart (2007).
Water  This is expected, because permanently saturated areas also can generate Dunne overland flow.However, variably saturated areas are distinct expressions of the transition zone between areas of recharge and discharge, between diffusive transport and perennial water erosion.Could there be unique controls on variably saturated extent that are not captured in Figure 10?
The answer to this question may lie in connection to the hillslope number.In Section 4, we discussed how Litwin et al. (2021) called β = h g /h a the hillslope number, which is defined as hillslope relief divided by the aquifer thickness.However the relief and mean aquifer thickness are emergent products of our coevolving system that we cannot specify ahead of time.We found that the actual emergent hillslope number has some bearing on the extent of variably saturated areas.
Before plotting the hillslope number, we first plotted the proportion of the domain classified as variably saturated against the dimensionless mean topographic variance (Figure 11a).The pattern is similar to that shown in Figure 10b, but with more scatter.The scatter shows greater difference in topographic variance between model runs that have the same VSA but different event storage index σ (among other factors).
Dividing mean topographic variance by the actual mean aquifer thickness 〈h〉 rather than the characteristic height scale h g gives an estimate of the emergent hillslope number, Z/〈h〉, on the y-axis.Figure 11b shows that this produces three tight relationships, separated by differences in aquifer relief index β.
The hillslope numbers that we observe for the high β case are within the range described by Lyon and Troch (2007), who calculated hillslope numbers in the range of 18-96 for several real sites, although this will be sensitive to exactly how the relief is defined.
The importance of β in Figure 11b is expected given its role as a type of characteristic hillslope number based on model parameters.By normalizing the hillslope number with β we obtain a relationship (Figure 11c) that is tighter than the original between VSA and topographic variance (Figure 11a).This indicates that there is a trade-off between the hillslope number and the proportion that is variably saturated: larger normalized hillslope numbers, which are associated with thin aquifers relative to relief, emerge with smaller variably saturated areas; thicker aquifers relative to relief emerge with greater variably saturated areas.The hillslope number has proved to be a useful concept to understand hydrologic response (Lyon & Troch, 2007), and here reveals a connection with emergent landscape features that to our knowledge, has not been shown before.We will not attempt to explain why this relationship exists here, but it certainly demonstrates that there are rich and largely unexplored avenues of research in emergent hydrogeomorphic dynamics.

Does Coevolution Explain Freeze's Observation About the Prevalence of Dunne Overland Flow?
Freeze (1980) observed a "delicate hydrologic balance on a hillslope" where only narrow combinations of parameters produced Dunne overland flow, despite its prevalence in nature and the wide plausible ranges of parameter values.This led Freeze to hypothesize that there is a "very close relationship between climate, hydraulic conductivity, and the development of geomorphic landforms," in nature that leads to preference for Dunne overland flow.We also suggested that there is delicate balance in landscapes, but we did not see a tendency toward Dunne overland flow.Instead at geomorphic dynamic equilibrium water is partitioned to maintain balance between the size of recharge and discharge areas (Section 6.1).Consequently, coevolution reinforces the associations described by the Dunne diagram: places with thick, highly permeable soils evolve steep topography and subsurface-dominated runoff generation, whereas places with thinner less permeable soils evolve gentler topography and more Dunne overland flow (Section 6.5).
However, our results are not conclusive evidence against Freeze's hypothesis.Although we explored the role of coevolving topography ("geomorphic landforms"), we did not account for the coevolution of climate and subsurface properties.Li et al. (2014) found evidence that the coevolution between subsurface properties and climate is important for runoff generation.They conducted a comprehensive study of the prevalence of runoff generation mechanisms using synthetic watersheds where climate, subsurface parameters, and topographic relief were varied independently.Model runs that conformed to the discharge-ET partitioning described by the Budyko hypothesis were also more likely to conform to the topography-runoff generation relationship in the Dunne diagram.
Although the tendency of watersheds to conform to the Budyko hypothesis is not fully understood, it has been associated with coevolution between soil, vegetation, and climate (Troch et al., 2013).
Our study suggests that hydrogeomorphic constraints could complement the Budyko water balance constraint.Li et al. (2014) varied topography in their synthetic watersheds by stretching the vertical dimension of a digital elevation model, effectively decoupling catchment morphology from other attributes, as they intended.Perhaps something like Equation 20 that relates transmissivity and climate to source area size and convergence, or the relationship between relief and quickflow fraction in Figure 10 could be used to provide a hydrogeomorphic behavioral constraint to complement the water balance constraint.Such a constraint would need to be grounded in field evidence, but could use relationships derived from our simulations as plausible hypotheses.This could also be useful for constraining parameters in large-scale predictive models, where the appropriate values of subsurface parameters are unknown and difficult to measure.

Conclusions
Landscape evolution models are powerful tools for understanding the surface processes, acting as testing grounds for theories about how tectonics, climate, and lithology affect geomorphic features we observe today.Hydrology is often the glue that links these forcings and features together, as water is a powerful and ubiquitous agent for transporting solid and dissolved material from headwaters to depocenters.Here we have shown that LEMs have the potential to provide insights into the emergence of hydrological processes as well, provided the mechanisms underlying those processes are resolved in sufficient detail.
We have shown just one potential avenue for using an LEM to answer hydrological questions, in which runoff from shallow groundwater and precipitation on saturated areas provides the shear stress for detachment-limited erosion.Within this scope, we have revealed complex interactions between topography, aquifer properties, and hydrologic function, including water balance partitioning, patterns of recharge and saturated areas, and the emergence of variable source area runoff generation.Most importantly, we found that.
1.At dynamic equilibrium, the size of diffusion-dominated uplands that saturate very infrequently evolves to supply enough recharge to lowlands that they can remain near-saturated, and so experience erosion by surface water sufficient to remove the colluvium supplied by diffusive transport from upslope. 2. As a consequence, drainage dissection increases not only with decreasing drainage capacity, as shown by Litwin et al. (2021), but also with hydroclimatic properties including the subsurface storage capacity relative to storm depth and the aridity index.3. When aquifers are thick relative to relief, aquifer gradients can become decoupled from topographic gradients and unexpected features like wide valley bottom wetlands and trellis-like drainage networks can emerge, which are not possible in simpler models that do not simulate aquifer evolution.4. All of the model simulations we conducted collapse on the same log-linear inverse relationship between dimensionless local relief and the ratio of average catchment quickflow to discharge.Because quickflow in our model is primarily generated by saturation excess overland flow, this result suggests that the relationship between hydrology and geomorphology on the humid side of the Dunne Diagram, in which landscapes with deep soil and steep topography are associated with subsurface flow, and gentle topography and poorly drained soils are associated with saturation excess overland flow, can emerge as a result of coevolution.
5. There is a nonlinear inverse relationship between the extent of variably saturated areas and the Hillslope number, which describes the local relief relative to the aquifer thickness.Further work is needed to understand this and the previous relationships and examine whether they are supported by empirical data.
This study lays the foundations for future work in which LEMs can be used to ask hydrological questions, and dynamic hydrological processes are given more consideration in spatially resolved LEMs.A forthcoming contribution will build on this foundation by parameterizing our model for several field sites and comparing emergent topographic and hydrological metrics.

Notation
Variable definitions are below, with dimensions length L, time T, and mass M. Prime always indicates the dimensionless equivalent, where dimensionless equivalents are defined in the text.where the depth to the water table is b h(x, y, t), the permeable thickness minus the aquifer thickness.We set the maximum profile depth equal to the permeable thickness b, such that d ≤ b, which ensured continuity between saturated and unsaturated zone models.Note that the recharge rate in A3 is equal to the precipitation rate i when the water table is at the surface (b h = 0).The groundwater model discussed in Section 2.4 then determines how this recharge will be partitioned between overland flow and saturated subsurface flow.A full sample calculation of the Schenk model is shown in Figure S2

Appendix D: Details for Hydrological Analysis
We ran simulations for 2000t g , by which time most simulations had reached an equilibrium where mean relief was no longer increasing (Figure S1 in Supporting Information S1).Cases not reaching equilibrium tend to be arid and have poorly developed drainage networks.Once the landscape evolution simulation had completed, we ran the hydrological component of the model (without changing the topography) for 2,000(〈tr〉 + 〈t b 〉) using the final water table as an initial condition to collect more detailed information on hydrological state and fluxes.Spatially distributed output (saturated area, recharge) were recorded at the storm-interstorm timescale, while spatiallylumped data (water balance components, total saturated area, total storage) were recorded at intervals corresponding to 1% of maximum timestep for groundwater model stability (Litwin et al., 2020).
In our analysis, we further divide discharge into fast and slow responding components (quickflow and baseflow).Discharge during interstorm periods is defined as entirely baseflow, whereas baseflow during storm events is estimated by linear interpolation between the pre-storm-event discharge and the post-storm-event discharge.This approach works for our model because all runoff generated during storm events is instantaneously routed to the outlet (Equation 13), leaving only the slowly varying exfiltration to leave as runoff during interstorm periods; see Figure S3 in Supporting Information S1 for an example.Quickflow is then some combination of exfiltration and precipitation on saturated areas.Although the regularization function in Equation 12 makes it difficult to isolate their respective contributions, precipitation on saturated areas is usually the dominant contribution to quickflow as the model lacks mechanisms that would rapidly increase exfiltration during storm events.
Zenodo (Litwin et al., 2023b).Landlab v2.0 (Barnhart et al., 2020) is a core dependency of DupuitLEM.The complete list of input parameter values can be found in Table S1, and in the model output archive.

Figure 2 .
Figure 2. (a) Hillshades of modeled topography with varying aridity index Ai and event storage index σ, showing strong declines in dissection with increasing aridity, and a weaker positive relationship between dissection and σ.Here γ = 4.0, β = 0.5, ρ = 0.03, ϕ = 1.5, α = 0.15, S c = 0.5, λ = 250, and δ = 2.0e 5.The simulation numbers are in the upper left hand corner.(b)-(e) Lateral transects through topography along the red dashed lines for model runs corresponding to the small numbers on subplots in (a).Gray areas are impermeable bedrock, and brown areas are regolith, which is shown behind the mean aquifer thickness in light blue.Note differences in the vertical scales of (b)-(e).

Figure 4 .
Figure 4. Classification of surface saturation for simulations with different aridity index Ai and event storage index σ (the same simulations as in Figure2, where the aquifer relief index β = 0.5).Surface saturation is determined on the basis of timevariable water table proximity to the surface.Locations are classified as dry if they experience surface saturation at <5% of the ends of storms and interstorms, and are classified as wet if they are saturated at the end of >95% of storms and interstorms.Variably saturated areas are everywhere that does not meet either of these criteria.The results show the extent of variably saturated areas is greatest when σ is small.Non-permanent streams emerge in some cases as aridity increases, including simulations 4 and 5, in which the channel network contains discontinuous reaches that are always wet.

Figure 5 .
Figure 5. Classification of surface saturation for simulations with different aridity index Ai and event storage index σ (the same simulations as in Figure S5 of Supporting Information S1 and Figure 3d-3f), where the aquifer relief index β = 0.05.The classes are the same as in Figure4.In contrast to the higher β case, here we see variably saturated areas from valleys to ridges in much of this parameter space.In the transition zone between widespread variable saturation and the zone without any channels, we see unusual channel forms, including trellis-like drainages(10, 16, 22, 23, 28), and extensive valley bottom wetland zones(11, 17, 23, 24, 28).

Figure 6 .
Figure 6.Dimensionless discharge Q* and baseflow Q * b versus dimensionless saturated area A * sat for simulations with different aridity index Ai and event storage index σ (the same that appear in Figure 2).A * sat is calculated using the same saturation criteria as all other figures, and has a maximum value of 1 when all cells are saturated.Each point depicts a model timestep, recorded at intervals corresponding to 1% of maximum timestep for groundwater model stability.Dimensionless discharge Q* is depicted in gray.Dimensionless baseflow Q * b is colored by the dimensionless storage S*, which varies from 1 when aquifer thickness is 0 everywhere to 1 when aquifer thickness is equal to permeable thickness everywhere.All quantities are totals of the model domain, and normalized by total area (we have left off total flux subscripts for simplicity of notation).Gray lines are provided as references for scaling between Q * b and A * sat .Data are absent for simulations 13, 19, 20, 26, 27, and 32-34 because surface runoff was not produced.

Figure 7 .
Figure 7. Detailed view of the dimensionless baseflow versus saturated area plot presented in Figure 6 subplot 4. (b)-(e) The spatial show the spatial distribution of saturation (in blue) at model timesteps corresponding to the locations (b)-(e) in panel (a).Panel (b) shows saturation in second-order channels.(c) Falls right at the inflection point of the saturation-baseflow relationship, and corresponds to saturation just beginning to extend beyond the first-order channel network.Panel (d) shows more extensive in unchanneled concave areas, while (e) shows widespread saturation on concave and planar slopes.

Figure 8 .
Figure 8. (Upper panel) conceptualization of how hillslope morphology and hydrology interact in DupuitLEM.The "perennial aquifer" corresponds to 95% exceedance probability for aquifer thickness, and the "variable water table" corresponds to 5% exceedance probability (Lower panel) selected model results, where key dimensionless parameters have been varied from a base case, which is shown in the middle.The bold headings should be read as the difference between the model results at the tip and tail of the gray arrows.For example, the upper left result is the same as the base case, except for increased transmissivity (associated with increasing γ).Likewise, the lower right panel is the same as the upper right panel, except for increased storm frequency and decreased storm size (increased σ).

Figure 10 .
Figure 10.(a) The Dunne Diagram, reproduced from Li et al. (2014), highlighting the humid environments (red dashed box).(b) The relationship between the quickflow fraction and the mean relief, normalized by the characteristic height scale h g for three sets of parameters, varying σ, γ, and β.All model runs are for humid climates (Ai = 0.5).Other parameters are the same as those used previously (ρ = 0.03, ϕ = 1.5, α = 0.15, λ = 250, S c = 0.5, δ varies from 1e 4 to 4e 6 with β).Colors indicate the fraction of the landscape classified as variably saturated under the definition used in previous sections.(c) The same plot as (b), but colored to show the particular combination of the three parameters varied.Dot size scales with σ, color lightness with γ, and base color with β.

Figure 11 .
Figure 11.(a) Variation in dimensionless relief Z/ h g versus the variably saturated area as a fraction of the total area based on the definition introduced in Section 5.3 using the same model runs shown in Figure 10, showing substantially more scatter than Figure 10a.(b) The hillslope number Z/〈h〉 plotted against the proportion variably saturated, showing parallel but distinct relationships for each value of β.(c) Normalizing the vertical axis in (b) by βcollapses all the relationships, and shows that the β-normalized hillslope number is maximized when the fraction of the watershed that is variably saturated is small, and minimized when the variably saturated fraction is large.
of Supporting Information S1.Additional description of the model implementation can be found in Text S1 of Supporting Information S1.In order to uniquely determine the mean storm duration and interstorm duration Equations 8 and 9, we introduce one final parameter:ρ = 〈t r 〉 〈t r 〉 + 〈t b 〉 Precipitation steadiness index (C12)Likewise, the simplified expression for storm recharge R Equation A2 is:R d = I min(dn a S d ,I) (C13)and the equivalent dimensionless form is:R′ d = I′ min( d′bn a p(〈t r 〉 + 〈t r 〉) S′ d ,I′) (C14)where R d = R′ d p(〈t r 〉 + 〈t r 〉) .When the dimensionless parameter definitions are substituted, this becomes:R′ d = I′ min(d′σϕ S′ d ,I′).(C15)

Table 1
Characteristic Scales That Were Derived to Isolate the Dimensions (Length, Height, and Time) of the Model

Table 2
The Dimensionless Groups That Appear in the Dimensionless Equations 14-19

Table 3
The Additional Dimensionless Groups Needed to Describe the Climatic and Vadose Models LITWIN ET AL.
also observed lower variability in baseflow discharge and saturated areas in a steeper watershed with deeper and more permeable soils (Sleepers River Watershed WC-4, LITWIN ET AL.
above depth d over the time interval Δt is the lesser of the available vadose storage above d and the depth of rainfall during the interval, minus the lesser of the water in vadose storage above d and the evapotranspiration during the interval.We assumed that the recharge R d received by a water table at depth d is the amount of water that has infiltrated below d in the vadose profile: LITWIN ET AL. stored