The fluid system within a critically tapering orogenic wedge is governed by complex interactions between topographic drive, thermal gradients, prograde dehydration reactions, internal structure, and regional tectonic compaction. Despite this complexity, topography is widely known to be the primary driving potential responsible for basin-scale fluid migration within the upper 7–10 km of an orogenic wedge. In recent years, investigators have revisited the problem of basin-scale fluid flow with an emphasis on depth-decaying permeability, which is a geologic phenomenon that is seldom accounted for in regional flow models. These recent investigations have shown that depth-dependent permeability at the basin scale strongly influences the relationship between local- and regional-scale flow paths. Here we investigate topography driven fluid flow within an orogenic wedge using a numerical modeling experiment designed to assess first-order fluid system behavior when permeability decreases systematically with depth. Critical taper theory is invoked to define two-dimensional basin geometry, and three subaerially exposed orogenic wedge models are presented with critical taper angles of 2°, 4°, and 10°. To assess the combined influence of topographic slope and depth-dependent permeability, a constant rate infiltration is applied at the wedge surface and a transient simulation is performed within each model for 20 m.y. Our results suggest that (1) depth-dependent permeability severely limits the penetration depth of infiltrating water within broadly tapering orogenic wedge systems, (2) fluid system evolution within a narrowly tapering orogenic wedge (i.e., ≤2°) is governed by local-scale topography superimposed on the regional gradient, (3) the influence of subbasin topography on local-scale fluid circulation is suppressed as the regional topographic gradient increases, and (4) the spatial distribution of groundwater residence time is fundamentally different when topographic slope exceeds 3°.


Investigations into basin-scale fluid system structure can be traced back to the pioneering works of Tóth (1962, 1963) that provided (1) a theoretical basis for topographic controls on regional fluid system structure, and (2) demonstrated that local-scale topographic relief within a drainage basin superimposes local-scale circulation patterns on the regional fluid system. Building upon Tóth’s (1962, 1963) initial contributions, Freeze and Witherspoon (1967) altered the synthetic Tóth basin to include layered heterogeneity and formation anisotropy, the results of which were shown to have a dramatic influence on regional groundwater systems. In the following decades, investigators have expanded greatly on Tóth’s (1962, 1963) contributions by deploying numerical simulation to investigate regional-scale fluid system structure in response to a wide variety of geologic phenomena. For example, Garven and Freeze (1984a, 1984b) used nonisothermal groundwater models to quantify the contribution of gravity driven groundwater flow for the development of Mississippi Valley–type sulfide deposits. Garven (1989) combined the flow equations for heat, groundwater, and oleic phase fluid to demonstrate the influence of regional-scale groundwater flow for the emplacement of oil sand deposits in the Western Canada Sedimentary Basin. Coupled hydromechanical models have also been used to investigate relationships between tectonic compaction and fluid overpressure within foreland basins, the results of which highlight the relationship between fluid overpressure and seismicity along basal detachments (Ge and Screaton, 2005) and provide important constraints on orogen-scale groundwater discharge rates (Ge and Garven 1989, 1994). In addition, crustal-scale, nonisothermal fluid flow models have been used to quantify the distribution of metamorphic fluids at convergent plate boundaries, while showing that thermally produced orogenic fluids strongly influence formation permeability through geologic processes such as mineral precipitation and host-rock alteration (Nunn and Demming, 1991; Cutillo et al., 2003; Lyubetskaya and Ague, 2009).

Investigators have been revisiting generalized models of regional-scale groundwater flow in order to understand the effects of depth-dependent permeability decay, which is a geologic phenomenon that has been convincingly observed (Manning and Ingebritsen, 1999; Saar and Manga, 2004). The principal contributions in this field are investigations that revisit the classic Tóth (1962, 1963) basin, which is a synthetic, two-dimensional (2-D) model domain composed of a ∼1.15° regional topographic slope superimposed with sinusoidal subbasins. Jiang et al. (2009) reconstructed the Tóth basin and showed that depth-decaying permeability inhibits regional-scale groundwater flow patterns, while increasing the penetration depth of subbasin-scale flow systems. On the basis of this result, subsequent work in the Tóth basin has shown that depth-decaying permeability also has a pronounced effect on the distribution of groundwater residence time, where increasing the permeability decay exponent increases groundwater residence times (Jiang et al., 2010; Cardenas and Jiang, 2010; Zlotnik et al., 2011). Jiang et al. (2011) showed that depth-decaying permeability will enhance the compartmentalization of groundwater age distributions, and, if the decay exponent is great enough, and if the basin is relatively shallow, then regional-scale fluid flow can be completely shut off, such that all groundwater flow is due solely to local-scale flow systems.

The investigations by Jiang et al. (2009, 2010, 2011), Cardenas and Jiang (2010), and Zlotnik et al. (2011) show that depth-dependent permeability exerts a first order control on fluid system structure within the Tóth basin. However, these investigations relate solely to the effects of depth-decaying permeability within a model domain of constant geometry, i.e., regional gradient and subbasin topography. In order to further our understanding of depth-decaying permeability in geologic environments, we examine regional-scale groundwater flow within a subaerially exposed orogenic wedge subject to depth-decaying permeability and varying topographic slopes. We use a numerical modeling experiment to evaluate fluid system evolution within three synthetic orogenic wedge models composed of 2°, 4°, and 10° taper angles. Each model domain is parameterized with identical depth-decaying permeability structure, and water infiltrates the wedge surface at a constant rate. Our results highlight the differences between local- and regional-scale fluid system evolution with increasing topographic slope and depth-decaying permeability. The distinction between these two scales of fluid system behavior is particularly important in orogenic belts, which are known to exhibit complex, spatial-temporal variations in fluid system behavior and structure (e.g., Evans, 2010; Cooley et al., 2011; Fitz-Diaz et al., 2011; Vandeginste et al., 2012). Before we can reliably decipher the interconnected diagenetic, thermal, paleomagnetic, deformational, and fluid generation and/or migration history of an orogenic belt, we must first constrain the scale over which these processes operate. As our results demonstrate, these scales are intimately associated with the permeability distribution and regional topographic slope of an orogenic wedge.


Fluid system evolution within a subaerial orogenic wedge is governed by dynamic interactions between topography driven fluid flow, thermal gradients, prograde dehydration reactions, internal geologic structure, and regional tectonic compaction. Despite this complexity, topography is the primary driving potential responsible for basin-scale fluid migration within the upper 7–10 km of an orogenic wedge (e.g., Garven, 1995; Lyubetskaya and Ague, 2009). In these settings, topographic relief may arise in response to continental collisions, erosion, and isostatic rebound (England and Molnar, 1990; Davis et al., 1983). As convergent margins undergo uplift, regional-scale fluid systems develop that can drive fluids tens to hundreds of kilometers away from their source, sometimes resulting in large-scale ore, mineral, and/or hydrocarbon accumulations (Garven and Freeze, 1984a, 1984b; Garven, 1989; Bethke and Marshak, 1990).

The topographic expression of an orogenic wedge is largely dependent on rock strength within the wedge and at the basal detachment. This relationship is described by the critical taper angle, which is the sum of the topographic and basal detachment dip angles. Davis et al. (1983) showed that critical taper angle results from the interaction between basal friction along the detachment horizon and internal wedge strength; increasing basal friction results in a larger critical taper angle, whereas a stronger wedge results in a smaller critical taper angle. An orogenic wedge underlain by strong basal lithologies may achieve a critical taper angle of 8° or more; e.g., the Taiwan fold and thrust belt is stable with a critical taper angle of 9° (Davis and Engelder, 1985). Within these broadly tapering orogenic wedges, the detachment more effectively resists failure, which results in fold-thrust belts, crustal thickening, and a larger critical taper angle. In contrast, orogenic wedge systems underlain by weak basal lithologies (e.g., the salt-detached Jura foreland basin) are characterized by fold-dominant internal structure and critical taper angles <4° and in some cases ≤1° (Davis and Engelder, 1985; Davis et al., 1983).

The fluid system within a critically tapering orogenic wedge is defined by fluid type, fluid source, rock type, fluid flow paths, and all mechanical and chemical interactions between them at regional and local scales. Similarly, fluid system structure (i.e., flow paths, fluid fluxes, fluid distribution) is governed primarily by taper angle and internal geologic structure (Fig. 1). During uplift, a regional topographic gradient is established that increases hydraulic potential at the hinterland and results in regional-scale fluid flow from hinterland to foreland (Bethke and Marshak, 1990). This suggests that the magnitude of topographic relief across a wedge will control the magnitude of the hinterland to foreland flow, with thicker wedges having increased topographic drive when compared with thinner wedges. In addition, the presence of kilometer-scale and larger scale faults and folds can act as fluid migration pathways or barriers that will control how fluids move and are distributed within the wedge (Sibson, 1981; Caine et al., 1996; Travé et al., 2000; Fischer et al., 2009; Evans and Fischer, 2012). These large-scale fluid systems are important in societally significant geological processes ranging from hydrocarbon generation and accumulation (Oliver, 1986; Moretti et al., 1996; Garven, 1995; Moore et al., 2004; Dewever et al., 2013) to earthquake generation (Sibson, 1990, 1996), hydraulic fracturing (Agarwal et al., 1979; Rummel and Hansen, 1988), and the formation of hydrothermal and Mississippi Valley–type ore deposits (Bethke and Marshak, 1990; Garven at al., 1993).


The cross-sectional geometry of a contractional orogen is reasonably described as a doubly vergent wedge (Fig. 2A; Koons, 1990; Willett et al., 1993). In keeping with this description, the model domains used for this study represent the cross-sectional geometry of 3 idealized and subaerially exposed, foreland orogenic wedges with critical taper angles of 2°, 4°, and 10°. At the hinterland boundary, each wedge has a maximum vertical dimension of 3 km, a foreland-dipping topographic slope (α) of 1°, 3°, or 9°, and a hinterland-dipping basal detachment (β) of 1° (Figs. 2B–2D). Each model domain has a maximum lateral extent of 100 km; however, the lateral extent for each wedge is determined by the intersection between the topographic slope and basal detachment. The 10° wedge terminates ∼16 km from the hinterland boundary; the 4° wedge terminates ∼40 km from the hinterland boundary; and the 2° wedge terminates ∼80 km from the hinterland boundary (Figs. 2B–2D). Each model comprises a 2-D Cartesian grid with 50 m × 50 m grid cells and 50 m out-of-plane thickness, which results in a series of 2-D models comprising 85,174 grid cells for the 2° wedge; 50,313 for the 4° wedge; and 22,719 for the 10° wedge. The geology of each model is composed of two structurally isotropic subdomains representing either wedge interior or continental crust (Figs. 2B–2D).

To account for depth-dependent permeability decay, k(z), within each subdomain (wedge interior or crystalline basement), we use a scaling coefficient (ζ) to reduce near-surface permeability (ks) as a function of depth (z): 
Here, the formulation for ζ(z) invokes the permeability scaling model proposed by Saar and Manga (2004), which is a piecewise function based in part on the power law scaling model originally proposed by Manning and Ingebritsen (1999). In this model, the depth-dependent permeability scaling coefficient, ζ(z), is given by 
and illustrated graphically in Figure 3A. In this formulation, the decay constant for ζ(z ≤ 800 m) is equivalent to the A parameter described by Jiang et al. (2009) and defined over a range from 0 m–1 (homogeneous k) to 0.001 m–1 (strongly depth-decaying k). To complete the permeability structure for each subdomain, near-surface permeability values (ks) are assigned using generic estimates for fractured sedimentary rocks (wedge interior) and unfractured crystalline rocks (continental crust): the wedge interior is assigned a near-surface (ks) permeability of 5 × 10–14 m2, and the continental crust is assigned a permeability of 5 × 10–17 m2 (Freeze and Cherry, 1979). The resulting 2-D permeability structure of each orogenic wedge is shown in Figures 3B–3D. Similarly, the porosity of each subdomain is assigned using generic estimates such that the wedge interior porosity is 0.15 and the continental crust porosity is 0.10 (Freeze and Cherry, 1979).

Before simulations, initial conditions were computed to establish the hydrostatic pressure distribution for each model domain. The initial conditions represent fully saturated conditions ranging from atmospheric pressure (0.101 MPa) at the model surface to ∼29 MPa at 3 km depth, and follow a linear surface-parallel gradient (Fig. 4A). Groundwater flow is simulated in each model for 20 m.y., which provides sufficient time to achieve steady flow within each wedge scenario. The code selection for this work is TOUGH2-MP (Zhang et al., 2008) compiled with EOS1 (Pruess et al., 1999), the equation of state module for simulating nonisothermal, multiphase conditions, including a two waters formulation (i.e., water and water with tracer). Dirichlet boundaries are imposed at the wedge surface to maintain atmospheric pressure, and along the foreland edge of the domain to maintain the hydrostatic pressure gradient in the far field. In keeping with other basin-scale flow models, adiabatic (no flux) boundaries are imposed across the bottom and hinterland edge of each model domain (Fig. 4B). In a physical sense, the hinterland boundary condition represents a permeability contrast between sedimentary rocks within the orogenic wedge and comparatively much lower permeability rocks within the interior of the orogen. In order to isolate for the combined effects of depth-decaying permeability and topographic slope, all simulations are isothermal at 20 °C and temperature effects are not taken into account. Although thermal effects will likely influence the deep fluid system at the orogen scale, preliminary simulations showed that differences in fluid system architecture between each wedge scenario are most pronounced within 1 km of the wedge surface where thermal effects are likely a second order process in comparison to topographic driving potential and permeability decay of 10–2.

In order to simulate the effects of infiltrating meteoric water, a source term is used in the row of grid cells immediately below the surface boundary (Fig. 4, blue line). Within each infiltration grid cell, water enters the system at a constant mass rate () of 7.7 × 10–2 kg/s, which is equivalent to a recharge flux (R) of 974 mm/yr on the basis that each grid cell surface area (A) is 50 m × 50 m in horizontal extent, and water density (ρw) at 20 °C is 998 kg/m3, e.g., R = /(ρwA). This recharge flux is based on simulations that indicate groundwater recharge rates of 300–1000 mm/yr across a number of orogenic belts (Wada et al., 2010, fig. 1A therein). In order to prevent infiltration from artificially migrating upward into the constant pressure boundary, the model mesh is augmented to change the connection distance between the infiltration cell and its downward interface from 25 m to an infinitesimally small value. This effectively places the infiltrating meteoric water at the interface between the source cell and its underlying neighbor, which is at a depth of 75 m below the wedge surface and suitable for assessing kilometer-scale fluid system evolution over geologic time. In addition, this invocation of source blocks obviates the complexity associated with boundary conditions of the third kind (pressure and flux). During simulations, the two waters function of TOUGH2-MP/EOS1 (Zhang et al., 2008; Pruess et al., 1999) is invoked so that saturations for the infiltrating water and initial water mass are separately tracked through space and time within each model domain. In the two waters formulation, individual mass balances are solved for each water, while keeping identical thermophysical water properties within each grid cell. Due to the high Péclet number (≥2 × 105) for the simulations presented here, molecular diffusion is not taken into account and mixing between infiltrating water and the initial water mass is purely a function of advection.


The model results presented here assess the combined influence of varying topographic slope and depth-decaying permeability on fluid system structure within a subaerially exposed orogenic wedge. Figure 5 presents the overall fluid system architecture for each wedge scenario (α = 1°, 3°, and 9°) after 16 k.y., by which time distinct differences in each wedge scenario are apparent. The temporal evolution of meteoric water through each wedge scenario is presented in Figure 6, and the distribution of groundwater residence time at 20 m.y. is presented in Figure 7. These model results reveal that depth-decaying permeability and increasing regional topographic slope act in opposition to each other, and in the following we explore the implications of these effects in terms of (1) local- and regional-scale controls on fluid system structure, and (2) the temporal evolution of this structure, as identified by meteoric water distribution and flow in each system.

Orogenic Fluid System Structure: Local or Regional Control?

Comparing model results from each wedge scenario (α = 1°, 3°, and 9°) after 16 k.y. reveals the intuitive findings that (1) greater topographic slope results in higher flow rates, and (2) groundwater generally flows from hinterland to foreland in response to the regional hydraulic head gradients (Figs. 5A–5C). Within the narrowly tapering (α = 1°) wedge model, however, the regional hinterland to foreland flow pattern is superimposed with local-scale groundwater circulation patterns that are controlled by local-scale topography along the wedge surface (Fig. 5A). In this modeling experiment, local-scale topography is a consequence of grid discretization along the wedge surface, which creates periodic 50 m elevation drops. Because the topographic slope (1°) in the narrowly tapering wedge is approximately the same as in the Tóth basin (1.15°), the periodic 50 m elevation drops create corresponding 50 m head drops, and the resulting local-scale groundwater circulation patterns (Fig. 5A, inset) are the same phenomenon predicted by Tóth (1962, 1963). Although these discretization-induced head drops are a consequence of grid discretization, the impact of these head drops influences the fluid system in the same manner as local-scale topography, thus providing an opportunity to assess the relative importance of local versus regional topography and depth-dependent permeability on fluid system structure.

The influence of depth-decaying permeability on regional- and local-scale groundwater flow paths has been assessed in recent investigations within the synthetic Tóth basin; both Jiang et al. (2009) and Cardenas and Jiang (2010) found that permeability decay with depth inhibits regional-scale groundwater flow paths and increases the penetration depth of local-scale flow patterns in the Tóth basin. The combination of sloping terrain and depth-decaying permeability results in decreasing permeability toward the highlands along any horizontal transect within the basin (Fig. 3A, inset). In terms of the orogenic wedge geometry presented here, water within the thicker portion of the wedge occurs within lower permeability rocks, and lateral hinterland to foreland flow paths are suppressed, which allows for enhanced development of local-scale groundwater circulation patterns. The narrowly tapering (α = 1°) wedge model has approximately the same topographic slope as the Tóth basin (∼1.15°), and as a result, the penetration depth of local circulation patterns is ∼57% deeper after 500 k.y. than an identical model comprised of homogeneous and isotropic internal wedge with k = ks (Van Dusen, 2014).

Although the works by Jiang et al. (2009) and Cardenas and Jiang (2010) are sufficient to explain fluid system structure within the narrowly tapering wedge model, the effects of depth-decaying permeability have not been investigated for basin-scale fluid systems with topographic slopes in excess of ∼1.15°. Within the moderately tapering (α = 3°) wedge model, the penetration depth of local-scale flow is much shallower than within the narrowly tapering wedge (Fig. 5B, inset), and local-scale groundwater circulation is almost entirely suppressed within the broadly tapering (α = 9°) wedge (Fig. 5C, inset). The inability of the 3° and 9° wedge models to develop well-defined local circulation patterns results from steeper topographic slopes generating much stronger hinterland to foreland head gradients (0.03 and 0.15 m m–1, respectively), which overwhelm the effects of local-scale topography across the wedge surface (Figs. 5B, 5C, inset). In addition, the steeper surface slope results in rapidly increasing permeability along horizontal transects toward the foreland (Fig. 3, inset). All three orogenic wedge models have identical permeability structure at the hinterland boundary (Figs. 3B–3D); however, the steep topographic slopes within the moderately and broadly tapering wedge models result in faster permeability gains over shorter lateral distances (Fig. 3A, inset). Instead of well-developed local-scale circulation patterns, fluid system structure within the moderately (α = 3°) and broadly (α = 9°) tapering wedge models is characterized by a regional-scale hinterland to foreland flow direction with a stratified flux magnitude that is surface parallel and decreases with depth (Figs. 5B, 5C). This stratification of flux magnitudes is a direct result of the depth-decaying permeability structure within the wedge interior for these models (α = 3° and 9°), and is not seen in identically parameterized wedge flow models with internally isotropic and homogeneous permeability structure (Van Dusen, 2014).

Time Variation of Fluid System Structure

To assess the temporal evolution of fluid system structure within each model scenario, the spatial distribution of infiltrating water is assessed at discrete time intervals (Fig. 6). Within the narrowly tapering (α = 1°) wedge model, the influence of local-scale groundwater circulation is apparent in the U-shaped patterns of high infiltrating water saturation (Fig. 6). As time progresses from 1 k.y. to 500 k.y., these U-shaped circulation patterns become progressively deeper, resulting in the transport of infiltrating water to greater depths within the model domain until the wedge is fully saturated with infiltrating water by 20 m.y. This suggests that transport of meteoric waters within the narrowly tapering (α = 1°) orogenic wedge is primarily governed by local-scale topographic variability across the wedge surface, and, as shown by Jiang et al. (2009), the effects of depth-decaying permeability enhance this phenomenon by increasing the penetration depth of local-scale circulation patterns. In addition, this result implies that meteoric waters can fully saturate a narrowly tapering orogenic wedge under long-term, geologically stable conditions.

Within the moderately tapering (α = 3°) wedge model, local-scale groundwater circulation patterns govern the distribution of infiltrating water in early time (<4 k.y.), after which the fate of this water is governed primarily by regional-scale hinterland to foreland flow (Fig. 6). Although this model suggests a theoretical capacity for infiltration to fully penetrate the moderately tapering (α = 3°) orogenic wedge, complete saturation from infiltrating waters occurs only within the upper 750 m after 20 m.y. while saturation at depth is limited to 0.8 (Fig. 6). This result arises due to the combined effects of regional topographic slope and depth-decaying permeability. The regional head gradient (∼0.03 m m–1) that develops with a 3° topographic slope is sufficient to lessen the penetration depth of local-scale fluid circulation patterns; this in turn decreases the penetration depth of infiltrating water along the wedge interior. However, this head gradient is not strong enough to prevent downward penetration along the hinterland boundary, the result of which allows infiltration to fully penetrate along this boundary. As meteoric water penetrates downward along the hinterland boundary, the depth-decaying permeability structure encourages lateral flow because of rapid permeability gains from hinterland to foreland as wedge thickness decreases (Fig. 3A, inset). For applications to real-world scenarios, this moderately tapering (α = 3°) wedge model suggests that geochemical interactions between meteoric waters and connate fluids are theoretically possible at any depth within the moderately tapering orogenic wedge because a non-negligible fraction of connate fluids will remain in the system over geologic time scales.

The temporal evolution of infiltrating water within the broadly tapering (α = 9°) wedge model is characterized by shallow penetration depth, comparatively low saturation, and a rapid return to steady flow (Fig. 6). The shallow penetration of infiltrating water results from the comparatively high regional hydraulic head gradient (0.15 m m–1), which overwhelms local-scale fluid circulation patterns with high horizontal flow rates. To complement this effect, depth-decaying permeability inhibits the diffusivity of infiltrating waters at depths >∼1000 m, where the permeability drop is approximately two orders of magnitude (Fig. 3). The combined effects of depth-decaying permeability and steep topographic slope encourage meteoric waters to accumulate at shallow depths (Fig. 6) and travel downslope under the regional head gradient (Fig. 5C). As a result, infiltration cannot penetrate a significant portion of the broadly tapering (α = 9°) orogenic wedge. This result is in stark contrast to results from an identically parameterized model with internally homogeneous permeability (no depth dependence), in which infiltrating water can fully penetrate and fully saturate each of the wedge scenarios discussed here (Fig. 6). These results may have important consequences for the development and spatial distribution of diagenetic and ore-forming processes in geologic settings with similar geometric configurations, as well as for the development of models that seek to understand such settings.

Groundwater residence time is an important criterion for understanding the relationship between fluid system structure and a number of regional-scale geologic processes, including ore formation, diagenesis, and oil sand emplacement (Garven, 1995). In order to quantify the spatial distribution of groundwater residence time (τ) within the model scenarios presented here, τ is computed within each grid cell using the traditional formulation: τ = Vw/Q, where Vw [L3] is the volume of water within each grid cell and Q [L3 T–1] is the net volumetric flow rate through the grid cell (L is length, and T is time). Using TOUGH2-MP/EOS1 (Zhang et al., 2008; Pruess et al., 1999) model output at 20 m.y. of simulation time, Q is computed for each grid cell by dividing mass flow rate () [M T–1] by water density (ρw) [M L–3]: Q = w (M is mass). Within the three model scenarios presented here, groundwater residence time ranges between 100–102 yr near the wedge surface to 105–106 yr at 2–3 km depth along the hinterland boundary (Figs. 7A–7C).

Within the narrowly tapering (α = 1°) wedge, the spatial distribution of groundwater residence time is primarily controlled by the local-scale fluid circulation patterns that characterize the fluid system structure (Fig. 7A). In particular, the subbasins defined by discretization-induced head drops comprise recharge and discharge cells (Fig. 5A, inset) that circulate comparatively younger waters (100–102 yr). These recharge and discharge cells are separated by zones of stagnation where groundwater residence time is >103 yr within 100 m of the wedge surface. This highly compartmentalized fluid system bears much resemblance to the work of Jiang et al. (2010), who showed that increasing the decay constant within the Tóth basin results in increasing compartmentalization, and a strongly heterogeneous groundwater age distribution.

Within the moderately and broadly tapering wedge models (α = 3° and 9°, respectively), groundwater residence time is characterized by a layered stratification of increasing age with depth (Figs. 7B, 7C); however, a slight influence of local-scale fluid circulation is apparent in the moderately tapering wedge. The groundwater residence time distribution in models with steeper regional topography is primarily a function of the flow rate within each grid cell, which is controlled by the combined effects of the topographic slope and depth-decaying permeability. As previously discussed, the comparatively strong regional head gradients within steeper wedges inhibit local-scale fluid circulation, while encouraging high horizontal flow rates within the upper kilometer of the wedge surface (Figs. 5B, 5C). However, the flow rates decrease rapidly with depth as a function of permeability decay, so that the spatial distribution of groundwater residence is similar in nature to the spatial distribution of permeability, i.e., surface parallel. This is apparent by comparing the permeability distributions for the α = 3° and 9° scenarios in Figures 3C and 3D, respectively, with the comparable groundwater residence time distributions in Figures 7B and 7C. Groundwater flow rates are lower in the moderately tapering wedge than in the broadly tapering wedge, and as a result groundwater residence time increases more rapidly in the moderately tapering wedge, although the general pattern of residence time is the same in both model scenarios.


In this paper a numerical modeling experiment is used to refine our conceptual understanding of orogen-scale groundwater flow within three synthetic model domains comprising geometry based on critical taper theory and subject to depth-dependent permeability. Although basin-scale fluid flow has a rich history in the literature, investigators (Jiang et al., 2009, 2010, 2011; Zlotnik et al., 2011) have revisited the classic Tóth basin and found that depth-decaying permeability plays an important role in the development and structure of regional- and local-scale fluid circulation patterns. The models presented here build on these recent investigations by illustrating how increasing topographic drive influences fluid system evolution under depth-decaying permeability structure.

In particular, we find that fluid system structure within the narrowly tapering orogenic wedge (α = 1°) is similar in nature to the Tóth basin; however, increasing the topographic slope (α) to 3° or more fundamentally changes the relationship between local- and regional-scale flow paths. As shown by Jiang et al. (2009), gently sloping terrain (∼1°) combined with depth-decaying permeability results in a fluid system structure with inhibited lateral flow paths and enhanced local-scale fluid circulation patterns. This is caused by a gradually decreasing permeability toward the highlands along any horizontal transect within the domain that effectively inhibits the development of regional-scale flow paths. As a result, subbasin topographic variability and local-scale flow paths govern fluid systems with low topographic gradients.

This work also shows that increasing topographic slope results in a rapid lateral permeability drop through the domain, which allows regional flow paths to govern fluid system evolution when topographic slope is >3°. In concert with this substantial topographic control, the influence of depth-decaying permeability plays an important role in the penetration depth of meteoric waters and the spatial distribution of groundwater residence time. In particular, broadly tapering orogenic wedge systems subject to depth-decaying permeability structure are characterized by (1) a reduction in the penetration depth of meteoric waters and (2) a decrease in groundwater residence time. Both of these phenomena are more pronounced with increasing topographic slope. Although these models are generic in nature, they provide a quantitative framework that enhances our understanding of topographic driving potential on fluid system evolution within critically tapering orogenic wedge settings subject to depth-decaying permeability.

This project received support from the Northern Illinois University College of Liberal Arts and Sciences, the Division of Research and Graduate Studies, and the Graduate Council Committee for Research and Artistry. Pollyea is grateful to Mark Frank for insightful conversations regarding the implications of this work for large-scale geologic processes. We thank Associate Editor Mike Williams and three anonymous reviewers for generous contributions of time and suggestions, the results of which have greatly enhanced the quality of this manuscript.