We applied multiple geochemical tracers (87Sr/86Sr, [Sr], δ13C, and δ18O) to waters and carbonates of the lower Colorado River system to evaluate its paleohydrology over the past 12 Ma. Modern springs in Grand Canyon reflect mixing of deeply derived (endogenic) fluids with meteoric (epigenic) recharge. Travertine (<1 Ma) and speleothems (2–4 Ma) yield 87Sr/86Sr and δ13C and δ18O values that overlap with associated water values, providing justification for use of carbonates as a proxy for the waters from which they were deposited. The Hualapai Limestone (12–6 Ma) and Bouse Formation (5.6–4.8 Ma) record paleohydrology immediately prior to and during integration of the Colorado River. The Hualapai Limestone was deposited from 12 Ma (new ash age) to 6 Ma; carbonates thicken eastward to ∼210 m toward the Grand Wash fault, suggesting that deposition was synchronous with fault slip. A fanning-dip geometry is suggested by correlation of ashes between subbasins using tephrochronology. New detrital-zircon ages are consistent with the “Muddy Creek constraint,” which posits that Grand Wash Trough was internally drained prior to 6 Ma, with limited or no Colorado Plateau detritus, and that Grand Wash basin was sedimentologically distinct from Gregg and Temple basins until after 6 Ma. New isotopic data from Hualapai Limestone of Grand Wash basin show values and ranges of 87Sr/86Sr, δ13C, and δ18O that are similar to Grand Canyon springs and travertines, suggesting a long-lived spring-fed lake/marsh system sourced from western Colorado Plateau groundwater. Progressive up-section decrease in 87Sr/86Sr and δ13C and increase in δ18O in the uppermost 50 m of the Hualapai Limestone indicate an increase in meteoric water relative to endogenic inputs, which we interpret to record progressively increased input of high-elevation Colorado Plateau groundwater from ca. 8 to 6 Ma. Grand Wash, Hualapai, Gregg, and Temple basins, although potentially connected by groundwater, were hydrochemically distinct basins before ca. 6 Ma. The 87Sr/86Sr, δ13C, and δ18O chemostratigraphic trends are compatible with a model for downward integration of Hualapai basins by groundwater sapping and lake spillover.

The Bouse Limestone (5.6–4.8 Ma) was also deposited in several hydrochemically distinct basins separated by bedrock divides. Northern Bouse basins (Cottonwood, Mojave, Havasu) have carbonate chemistry that is nonmarine. The 87Sr/86Sr data suggest that water in these basins was derived from mixing of high-87Sr/86Sr Lake Hualapai waters with lower-87Sr/86Sr, first-arriving “Colorado River” waters. Covariation trends of δ13C and δ18O suggest that newly integrated Grand Wash, Gregg, and Temple basin waters were integrated downward to the Cottonwood and Mojave basins at ca. 5–6 Ma. Southern, potentially younger Bouse basins are distinct hydrochemically from each other, which suggests incomplete mixing during continued downward integration of internally drained basins. Bouse carbonates display a southward trend toward less radiogenic 87Sr/86Sr values, higher [Sr], and heavier δ18O that we attribute to an increased proportion of Colorado River water through time plus increased evaporation from north to south. The δ13C and δ18O trends suggest alternating closed and open systems in progressively lower (southern) basins. We interpret existing data to permit the interpretation that the southernmost Blythe basin may have had intermittent mixing with marine water based on δ13C and δ18O covariation trends, sedimentology, and paleontology. [Sr] versus 87Sr/86Sr modeling suggests that southern Blythe basin 87Sr/86Sr values of ∼0.710–0.711 could be produced by 25%–75% seawater mixed with river water (depending on [Sr] assumptions) in a delta–marine estuary system.

We suggest several refinements to the “lake fill-and-spill” downward integration model for the Colorado River: (1) Lake Hualapai was fed by western Colorado Plateau groundwater from 12 to 8 Ma; (2) high-elevation Colorado Plateau groundwater was progressively introduced to Lake Hualapai from ca. 8 to 6 Ma; (3) Colorado River water arrived at ca. 5–6 Ma; and (4) the combined inputs led to downward integration by a combination of groundwater sapping and sequential lake spillover that first delivered Colorado Plateau water and detritus to the Salton Trough at ca. 5.3 Ma. We propose that the groundwater sapping mechanism strongly influenced lake evolution of the Hualapai and Bouse Limestones and that groundwater flow from the Colorado Plateau to Grand Wash Trough led to Colorado River integration.

Late Miocene to modern carbonate deposits along the lower Colorado River (Fig. 1) provide a rich source of information about the paleohydrology of spring, lake, and river waters in this region over the past 12 Ma. Recent models for the Bouse Formation have inferred a series of events that took place during a short period of time (bracketed between 5.6 and 4.8 Ma), when first-arriving Colorado River water caused a chain of lakes and internally drained basins to progressively fill and spill in a downward-propagating sequence (Blackwelder, 1934; Spencer and Patchett, 1997; House et al., 2005, 2008; House and Pearthree, 2014; Roskowski et al., 2010; Spencer et al., 2008, 2013). The Bouse Formation provides a record of these events, which some studies conclude took place in less than 50–100 ka (Spencer et al., 2008). The inferred “event stratigraphy” of the Bouse Formation challenges notions of Walther’s law, the resolution of available geochronology, and time scales needed for establishment of continental-scale river systems. Earlier models for deposition of all of the Bouse and Hualapai Limestones in a marine estuary (e.g., Metzger, 1968; Smith, 1970; Buising, 1988, 1990) have been more recently refuted based on isotopic and stratigraphic data (Spencer and Patchett, 1997; Roskowski et al., 2010; Spencer et al., 2013; Pearthree and House, 2014). However, a marine influence in the southernmost Bouse basins is still supported by some data sets (McDougall, 2008, 2010; McDougall and Miranda Martinez, 2014). Resolution of the debate about a possible marine connection for the southern Blythe basin is needed in order to test competing models for the magnitude of post–6 Ma uplift of the Colorado Plateau relative to the sea level (e.g., Lucchitta, 1979; Karlstrom et al., 2007, 2008).

This paper uses new and previously published data on the geochemistry of carbonates to test and refine models for initiation of the Colorado River by downward integration, and explore the possibility of a marine influence in southernmost exposures of the Bouse Formation. We applied multiple geochemical tracers—87Sr/86Sr; [Sr], δ13C, and δ18O—to characterize the geochemistry of waters and carbonates in the southwestern Colorado Plateau region over the past 12 Ma We first characterized modern springs and surface waters of the Grand Canyon region to develop a snapshot of the modern hydrologic system. We then extended this snapshot back in time to about ca. 1 Ma, by showing that Quaternary travertines associated with modern carbonic spring systems have an isotopic character similar to associated waters, and thus these carbonates provide a reasonable approximation of ancient spring-water chemistry. We extended this reasoning further back in time and compared the geochemistry of 2–4 Ma U-Pb–dated groundwater speleothems in Grand Canyon. Finally, we present new and previously published data for the geochemistry of carbonates from the late Miocene to early Pliocene Hualapai Limestone and Bouse Formation, including new measured sections and chemostratigraphy for the Hualapai Limestone. Our organization of topics “follows the water” in a general progression from upstream to downstream.

It is widely agreed that the Hualapai Limestone at the western mouth of Grand Canyon was deposited between 12 and 6 Ma, prior to the arrival of the Colorado River at ca. 5–6 Ma (Spencer et al., 2013). There is ongoing debate, however, regarding the exact age and early history of the Colorado River (Dorsey et al., 2007, 2011; House et al., 2008; Aslan et al., 2010; Spencer et al., 2013; Miller et al., 2014; Pearthree and House, 2014; McDougall and Miranda Martinez, 2014; Harvey, 2014). In this paper, we adopt the view that the first Colorado River waters reached the western edge of the Colorado Plateau after 6 Ma (e.g., the age of upper Hualapai Limestone; Spencer et al., 2001) and before the ca. 4.4 Ma river gravels at the western mouth of Grand Canyon (Faulds et al., 2001). The goals of this paper are to present a new model for the importance of groundwater during Colorado River integration, reinforce published models for downward integration of the river through a series of lakes, and depart from some recent models by showing that existing data permit an intermittent marine influence on the southernmost basin of the Bouse Formation.

Twenty-three new Sr isotope analyses were conducted on a Neptune multicollector–inductively coupled plasma–mass spectrometer (MC-ICP-MS) at the Radiogenic Isotope Geochemistry Laboratory in the Department of Earth and Planetary Sciences at the University of New Mexico. Samples of 30–40 g of calcite were dissolved in 3 N HNO3, run through a 200 µl Teflon column with Eichrom Sr-spec resin, dried, and then redissolved in 3% HNO3 acid. They were introduced into a Cetac Aridus II desolvating nebulizer and run against National Bureau of Standards Sr standard NBS-987, which has an 87Sr/86Sr value of 0.71025. Values are reported with 1σ (%) standard error. These data were combined with a comprehensive compilation of ∼200 87Sr/86Sr analyses on waters and ∼250 on carbonates from the region (summarized in Supplemental Table 11). Our simple two-component Sr isotope mixing models using the 87Sr/86Sr ratio and [Sr] concentrations follow the methods of Ingram and DePaolo (1993) and Crossey et al. (2006).

In total, 230 samples were run for stable oxygen and carbon isotopic analysis; these were drilled, and powders were flushed with He gas and then reacted for 24 h with H3PO4 at 50 °C (Spotl and Vennemann, 2003). Evolved CO2 was measured by continuous-flow–isotope ratio–mass spectrometry using a Gasbench device coupled to a Finnigan Mat Delta Plus Isotope Ratio Mass Spectrometer at the Stable Isotope Laboratory at the University of New Mexico. Oxygen isotope ratios in waters were determined using the CO2 equilibration technique and are reported using the standard delta notation versus Vienna standard mean ocean water (V-SMOW) with a precision of ±0.3‰. Hydrogen isotope ratios were measured using the continuous-flow high-temperature reduction method using a temperature conversion elemental analyzer (TC/EA) coupled to a Delta Plus XL Thermo-Finnigan mass spectrometer. Carbon results are reported in ‰ relative to Peedee belemnite (PDB) with a precision of ±0.3‰. These data were combined with a comprehensive compilation of water and carbonate stable isotope data from the region (summarized in Supplemental Table 1 [see footnote 1]).

Six tephra samples were sent to the U.S. Geological Survey (USGS) laboratory for chemical correlation. Three yielded sufficient glass shards to provide correlation information. Nine major elements (SiO2, Al2O3, FeO, MgO, MnO, CaO, TiO2, N2O, and K2O) were analyzed using the electron microprobe. The raw data were then recalculated to 100% on a fluid-free basis. Coefficient analysis was performed on the chemical data, and the normalized values were compared to a database (∼5800) of chemical “fingerprints” in the USGS Tephrochronology Project reference database. For complete tephrochronology interpretation, independent age control, stratigraphic position, field, petrologic characteristics, and mineralogy were also considered. Supplemental Table 22 summarizes raw and normalized chemical data, comparative chemical data, and lists of chemical correlatives.

U/Pb dating of carbonates may offer a new tool for refining the age and duration of deposition of the Hualapai Limestone and Bouse Formation. Samples were cut out with a rock saw, etched in weak (∼1 N) HNO3 for <30 s, washed with 18 MΩ H2O, dried, and then broken into small pieces with an agate mortar and pestle. Five representative pieces (each with a mass of 30–70 mg) were dissolved separately in 15 N HNO3, spiked with 205Pb, fluxed for >1 h, dried, dissolved in 6 N HCl, dried again, dissolved in 2 N HCl, dried once more, and then dissolved in 1 N HBr. Each sample was then cleaned, and U, Th, and Pb were separated using anion resins. U, Th, and Pb ratios were measured with Neptune MC-ICP-MS. All appropriate ratios were corrected for excess 206Pb from initial 234U by assuming an initial δ234U value of 750‰. Raw data were reduced using PBDATA (Ludwig, 1993) and a Pb blank of ∼10 pg. U-Pb ages were generated using ISOPLOT (Ludwig, 2001). Ages are reported as U-Pb concordia-linear dates, and data reported in Supplemental Table 33.

U-Pb analyses of individual detrital zircons were obtained by laser ablation–inductively coupled plasma–mass spectrometry (ICP-MS) at the University of Arizona LaserChron Center. Laser ablation was conducted with an Excimer laser beam diameter of 30 or 35 µm and a pulse frequency of 7 Hz. Measurements were performed with the Nu high-resolution (HR) MC-ICP-MS system (Gehrels et al., 2008; Gehrels and Pecha, 2014). Analysis sites were randomly targeted. Most zircon yielded U-Pb results that were concordant to within 10%. Complete data tables, evaluation of the 49,127 secondary standard data points, and the basis for making statistical comparisons between samples are provided in Supplemental Table 44. All location data for new geochronology samples and measured sections are in Supplemental Table 55.

Strontium Concentration and Isotopic Ratios

Calcite contains a wealth of semi-independent tracers. The 87Sr/86Sr ratios combined with Sr concentration, [Sr], in waters can faithfully record the nature of source waters and mixtures of different water sources because [Sr] abundance is variable, 87Sr/86Sr ratios are precisely measured, and because different source regions have characteristically different ranges of values. The 87Sr/86Sr ratios in carbonates tend to be a direct reflection of paleowater chemistry because Ca and Sr readily substitute in carbonate rocks. This is demonstrated by data that show similar 87Sr/86Sr values in modern waters and their associated travertine deposits (Crossey et al., 2006).

The 87Sr/86Sr ratio provides a key data set with which to evaluate modern waters and modern to ancient carbonate deposits (Spencer et al., 2013, and references therein), because 87Sr/86Sr is a conservative tracer, and values measured in carbonates mimic the 87Sr/86Sr composition of the water from which they were deposited (Ingram and DePaolo, 1993; Spencer and Patchett, 1997; Spencer et al., 2001; Li et al., 2008a; Crossey et al., 2006). However, we also note two significant caveats and potential problems with using 87Sr/86Sr as a single tracer without other information: (1) Addition of a small volume of high-87Sr/86Sr, high-[Sr] water to low-[Sr] water can have strong leverage on Sr isotope composition; and (2) well-mixed waters should have relatively narrow ranges of 87Sr/86Sr (0.0002 for Salton Sea; Li et al., 2008b), whereas lake and river systems with distinct subbasins and substantial input from different source waters often show greater variation; for example, Great Salt Lake waters vary by 0.006 (Fig. 2; Hart et al., 2004) and the San Francisco Bay estuary system varies by 0.004 (Ingram and Sloan, 1992; Ingram and DePaolo, 1993).

Figure 2 provides a synthesis of new and previously published data for 87Sr/86Sr in modern waters and Pliocene carbonates of the Colorado River system. The 87Sr/86Sr ratios from modern waters in endogenic springs sourced in granitic basement of western Grand Canyon are as high as 0.735. These spring waters mix within the regional groundwater system to produce the extremely wide range of 87Sr/86Sr measured in modern springs and travertines of Grand Canyon (0.708–0.735; Crossey et al., 2006). Speleothem values are 0.713–0.723 (Hill et al., 2008); these are from mammillary coatings on caves interpreted to form near the water table in phreatic groundwater cave systems. Hualapai Limestone 87Sr/86Sr ranges from 0.7114 to 0.7190, and Bouse carbonates are 0.7100–0.07125 (Spencer and Patchett, 1997; Roskowski et al., 2010). Values distinctly higher than marine values of 0.709 provide evidence that most or all of the Bouse and Hualapai Limestones were deposited in lacustrine rather than marine waters (Fig. 2; Spencer and Patchett, 1997). New results are discussed more fully later herein for modern waters and carbonates from the Hualapai Limestone and Bouse Formation.

Oxygen and Carbon Isotopes

The δ18O and δ13C values in travertines exhibit the widest variation in isotopic character of any rock type (Sharp, 2007). Prior work on δ18O and δ13C in waters and calcite of the Colorado River system includes studies by Faulds et al. (2001), Gross et al. (2001), Poulson and John (2003), and Roskowski et al. (2010).

Figure 3 shows covariation trends of δ18O and δ13C from previously published studies of carbonates precipitated from selected lakes and estuaries in the western United States. Marls that precipitate in the water column and settle as pelagic rain have long been used as proxies for the water chemistry; stromatolites, oncolites, and sometimes shells and tufas can also be used to estimate δ18O of waters, although their δ13C can be shifted due to fractionation effects (Talbot, 1990). The different covariation trends in closed basins can be used to characterize the hydrologic setting of internally drained basins (Talbot, 1990). The δ18O of lake waters primarily reflects elevation (latitude) of recharge and degree of evaporation (i.e., whether the lake is open or closed). The δ13C in southwestern lakes reflects surface water and groundwater mixing, and then local degassing and related fractionation of δ13C, primary biogenic activity, and decreasing surface area/depth of lakes (Fig. 3). As examples pertinent to this study, Great Salt Lake and Lake Bonneville data show a covariation trend (r2 = 0.89) that reflects a similar mix of recharge waters over the past 17 Ma and strong evaporative fractionation (Talbot, 1990). Shells from Holocene Lake Cahuilla, which formed by flooding of dominantly Colorado River water into the Salton Trough intermittently in the past 20 ka, show a covariation trend (r2 = 0.79) that records evaporation in a closed lake (Li et al., 2008a). Shells from the modern San Francisco Bay (and dating to 2 ka) show weaker covariation (r2 = 0.5) interpreted to reflect changes in salinity due to variable mixing of riverine water and seawater in an open (flow-through) marine estuary (Ingram et al., 1996).

A realignment of covariation trends can record changes in basin chemistry through time. For example, a marked up-section change is seen in the late Miocene (9–6 Ma) Ridge basin (Talbot, 1994). This basin is pertinent because it restores to a geographic position adjacent to the Eastern Transverse Ranges (Fig. 1), and its upper part was contemporaneous with the lower Imperial Formation (e.g., Dorsey et al., 2007). The lower part of the section has nearly constant δ18O, reflecting source waters with variation in δ13C interpreted to record rapid flow-through in an open lake. The upper covariation segment (r2 = 0.84) is interpreted to record evaporation in a closed lake (Fig. 3; Talbot, 1994). Another example of processes that change the position of collinear arrays is suggested by the difference in published stable isotope values between carbonates of the Hualapai Limestone and Bouse Formation (Roskowski et al., 2010), which reveal an overall shift to the lower right (Fig. 3). This shift may provide information about residence and surface area/depth ratios (e.g., Talbot, 1994) in successively lower lake basins of a developing Colorado River–fed chain of lakes, as discussed in more detail later herein.

Hydrochemistry of Modern Spring Waters

The geochemistry of modern waters in the Grand Canyon region is summarized in Figures 2 and 4, and data are listed in Supplemental Table 1 (see footnote 1; new data are in bold). Figure 2 shows 87Sr/86Sr values of the Colorado River, its tributaries, and springs, including labeled springs and spring groups that characterize the largest volume of groundwater inputs such as Fence Springs and Blue Springs of eastern Grand Canyon, Havasu Spring of central Grand Canyon, and Lava Warm Spring and Travertine Grotto Spring of western Grand Canyon. The highly variable 87Sr/86Sr values in waters are interpreted to reflect complex mixing of deeply sourced (endogenic) waters, karst aquifer groundwater, and meteoric (epigenic) waters (Crossey et al., 2006). Ranges of 87Sr/86Sr values get wider and values are more radiogenic in the western Grand Canyon (Crossey et al., 2009). [Sr] in Colorado River surface water increases from <1 ppm in the eastern Grand Canyon to 1.2 ppm in Lake Mead. Spring waters range from <1 ppm to as much as 15 ppm (Supplemental Table 1 [see footnote 1]). Compiled surface-water (tributary) inputs like the San Juan, Little Colorado, and Virgin Rivers are variable but significantly less radiogenic than groundwater inputs.

Figure 4 shows previously published and new δ18O and δ13C stable isotope data from Grand Canyon area waters (blue symbols). Measured water values (relative to SMOW) were converted to values for the carbonate (in PDB) that would precipitate from the water at 15 °C based on the models from Demeny et al. (2010) for δ18O and from Romanek et al. (1992) for δ13C. Similar to 87Sr/86Sr, δ18O and δ13C are highly variable in Grand Canyon groundwaters due to mixing. For example, a dominantly epigenic (meteoric) spring (Thunder River [TR]) differs from springs with strong endogenic component (Travertine Grotto [TG]). High-volume carbonic springs tend to get more negative in δ13C from east to west (Blue [B], Havasu [H], Lava [L], Travertine Grotto [TG]) at similar δ18O values (Fence Spring does not follow this pattern). Water-chemistry data have been used to model δ13C variation in terms of different sources of carbon for carbonic springs (Crossey et al., 2009). The main sources are: Ccarb = carbon derived from dissolution of carbonate such as from the karst aquifer; Corg = organic-sourced carbon such as from soil; and Cendo = endogenic sources of carbon such as from the mantle and magmatic sources. For the springs plotted in Figure 4, the Cendo component increases from east to west as follows: Thunder River = 0%; Fence = 27%; Havasu = 42%; Lava = 65%; and Travertine Grotto = 84%, in general agreement with 87Sr/86Sr data indicating increasing endogenic fluid inputs in western Grand Canyon springs. Endogenic inputs into the Colorado Plateau groundwater system are geochemically potent, but small volume, such that spring chemistry is quite variable from region to region (Crossey et al., 2006).

Paleogroundwater δ18O and δ13C from Travertine and Speleothems

Figure 4 also plots δ18O-δ13C values from young (<500 ka) Grand Canyon travertines that are associated spatially with the springs. Travertines show groupings that have heavier δ13C than their proximal springs, which we interpret to be due to degassing of the light 12C during travertine formation (rather than evaporation, as samples are from high-flow spring vents). The general similarity of arrays between water and travertine for a given region, and the differences between eastern and western Grand Canyon regions, shows that travertine can preserve paleohydrology information back to more than 500 ka.

Kwagunt travertines were deposited ca. 100–400 ka in eastern Grand Canyon (Pederson et al., 2002), and new isotopic data show a covariation array that is interpreted to reflect seasonal variation of groundwaters that were similar to Blue Springs, but with ∼8‰–10‰ heavier δ13C due to degassing during travertine formation. Havasu modern and Quaternary travertines (data from O’Brien, 2002; O’Brien et al., 2006) show ∼4‰ heavier δ13C relative to modern Havasu Spring for the same reason. Travertine Grotto travertines (new data) are also ∼5‰ heavier than the carbonates expected to be deposited by Travertine Grotto spring. Nevertheless, the full range of modern spring waters overlaps strongly with the range of Quaternary travertines and a new cave speleothem value from near Rampart Cave of westernmost Grand Canyon (Fig. 4), supporting the use of δ18O-δ13C values in carbonates as a proxy for the hydrochemistry of the water from which they were deposited, as has also been well established in other studies (Talbot, 1990).


The Hualapai Limestone (Blackwelder, 1934; Lucchitta, 1966, 1972; Blair and Armstrong, 1979) was deposited during the interval 12–6 Ma (Spencer et al., 2001) and is preserved in outcrops that roughly encompass the region of Lake Mead, from immediately west of the mouth of Grand Canyon to the upper end of Black Canyon (Fig. 5). These freshwater carbonates rest on locally derived fanglomerate and red siltstone that have variably been called Muddy Creek Formation (Lucchitta, 2013; see nomenclature discussion in Young, 2008), rocks of Grand Wash Trough (Bohannon, 1984), and rocks of Overton Arm (Beard et al., 2007). In contrast to Young (2008), we use the local names and recommend that propagation of the name Muddy Creek Formation outside of the Virgin River depression should be abandoned because of the sedimentary and tectonic complexity of internally drained basins of the Grand Wash region (see Dickinson et al., 2014). Here, we use the term Hualapai Limestone for the interbedded carbonates and red siliciclastic rocks above the first appearing limestone in each measured section.

Our new measured sections (for locations, see Fig. 5; Supplemental Table 5 [see footnote 5]) show that carbonates are interbedded with red siltstones in all sections (Fig. 6; Faulds et al., 2001; Wallace et al., 2005; Pearce, 2010) in a series of fault-bounded basins (Lucchitta, 1972; Seixas et al., 2015). The thickness of the Hualapai Limestone varies in different basins and reaches a maximum of 210 m in eastern Grand Wash basin near the mouth of Grand Canyon. This location makes it useful for understanding the history of integration of the Colorado River through western Grand Canyon.

Grand Wash basin contains the thickest section (up to 210 m) of Hualapai Limestone. The age constraints on deposition are shown in measured sections of Figure 6. An age of 13.11 Ma from an ash bed in red siltstones ∼100 m below the first Hualapai carbonate provides a maximum age on the Hualapai carbonate deposition. A new ash bed sampled in Grapevine Wash ∼8 m above the base of the Hualapai Limestone (Fig. 7) is dated at 12.07–11.31 Ma (see following), which now provides the best estimate for the onset of carbonate deposition in Grand Wash Trough. A 10.94 Ma ash age (Wallace et al., 2005) at the red siltstone to carbonate transition in a thinner section at Airport Point (Fig. 6) indicates a time-transgressive onset of carbonate deposition in different parts of the basin. An ash near the top of the Hualapai Limestone in the Grapevine Wash section was dated as 7.43 ± 0.22 Ma (Faulds et al., 2001). If deposition rates were constant, this suggests the 212 m of section could have been deposited between 12 Ma and 7.43 Ma at an average rate of 45 m/Ma Thus, available data suggest that carbonate deposition lasted from 12 to ca. 7 Ma in Grand Wash basin. Measured sections within and between basins are correlated from the exposed tops of the sections (Fig. 6). Temple basin (Mine Road section, Fig. 6) contains the location of the 5.97 Ma ash (Spencer et al., 2001) used for the age of the upper Hualapai Limestone. We recognize that there is a lack of geochronology to prove synchronous deposition, and, even if tops were contemporaneous at ca. 6 Ma, the potential for differential erosion makes this an approximate correlation. If spillover and downward integration models are applied, sedimentation rates may have been faster in Gregg and Temple basins than in Grand Wash basin.

Measured sections in the Hualapai Limestone (Fig. 6; Pearce, 2010) suggest wedge-shaped basin geometries with thickest accumulations on the immediate downthrown side of normal faults. Thickness and facies variations and soft sediment deformation features suggest that the carbonates were deposited syntectonically during slip on the Grand Wash, Wheeler, and Lakeside Mine normal faults (Pearce, 2010; Seixas et al., 2015). We summarize Hualapai Limestone chemostratigraphy and geochronology from east to west using multiple measured sections in each of three adjoining basins: (1) Grand Wash basin, bounded on the east by Grand Wash fault (GW of Fig. 1); (2) Gregg basin, bounded on the east by Wheeler fault (W of Fig. 1); and (3) Temple basin, bounded on the east by the Lakeside Mine fault (LM of Fig. 1). For present purposes, the Hualapai Limestone outcrops in Detrital basin are considered a western part of Temple basin that extends west to a drainage divide at the upper end of Black Canyon (Fig. 1; Roskowski et al., 2010; Spencer et al., 2013). Sedimentation rates were likely faster in deeper depocenters adjacent to active normal faults, and fanning-dip geometries need to be considered in correlation of measured sections.


As summarized in Supplemental Table 2 (see footnote 2), a tephra sample collected from near the base of the Hualapai Limestone in Grand Wash basin (JLP-08–54, Fig. 6) is located several meters above the carbonate–red siltstone contact (Fig. 6). This tephra is chemically similar to tephra derived from the Snake River caldera hotspot of southern Idaho, with closest matches having ages in the range of 12.07–11.31 Ma based on the 40Ar/39Ar dated Trapper Creek section (Supplemental Table 2 [see footnote 2]). The sandstone sample is ∼2 m below the tephra dated by Wallace et al. (2005) as 11.08 ± 0.27 Ma. The age range of 12.07–11.31 Ma overlaps with this ash but perhaps extends the range of carbonate deposition back to ca. 11.3–12 Ma.

Two additional ashes were sampled. One sample is from near the base of the Hualapai Limestone in the Temple basin, Mine Road section (JLP 08–109), which is ∼20 m below a 5.97 Ma ash (Spencer et al., 2001), and one is from the base of the Detrital Wash section (JLP-08–115; Fig. 6). These samples did not yield geochemical correlations with the USGS tephrochronology database, but they had close chemical correlations with each other (similarity coefficient of 0.947; see Supplemental Table 2 [see footnote 2]). This ash thus provides a local time marker that indicates that Detrital basin carbonate deposition began before 6 Ma, when about half (30 m) of the Hualapai Limestone section (carbonates plus red siltstones) had accumulated in the Temple basin. Thus, overall accumulation rates (and thicknesses) in Temple basin were higher than those of Detrital basin. This is compatible with models for increased down-section dips (tilt-fanning) of depositional time lines in each basin and is compatible with integration between basins ca. 6 Ma.

Pilot U-Pb Dating Attempts

Pilot attempts at U/Pb dating of carbonate were performed on several samples from the Hualapai Limestone. The best result came from sample JLP09-129C, which was sampled in Grapevine Wash and consists of a crystalline calcite infilling in micrite (the infilling was analyzed). As shown in Figure 8, this sample yields a three-dimensional isochron age that is consistent, within error, compared to available geochronology. It gives an age of 12.3 ± 3.4 Ma and is stratigraphically between ashes dated at 11 and 7 Ma (Faulds et al., 2001). This age is apparently too old, and its precision too low, to refine this particular section of the Hualapai Limestone. All aliquots yielded positive measured δ234U values of 13‰–27‰, also permissive of an (improbable) age of younger than ca. 2.5 Ma. Our current interpretation is that the positive δ234U values indicate secondary precipitation of young (<2.5 Ma) calcite, which will affect the U/Pb ages in, as yet, unconstrained ways. Until the U and Pb concentrations, isotopic composition, and any secondary calcite are better constrained and additional samples are dated in stratigraphic context, this age should be taken as a promising pilot sample, but not an accurate age.

Detrital-Zircon Data

The histories of present and past drainage systems can be interpreted using U-Pb dating of detrital zircons. Recent advances have emerged from increasingly well-understood detrital-zircon signatures for the region. Pertinent data from various units are shown in Figure 9. Cambrian strata of the Colorado Plateau show the restricted Yavapai-Mazatzal (1.7–1.8 Ga) plus 1.46 Ga “twin age peaks” that reflect erosion from Precambrian basement of the transcontinental arch, whereas Upper Paleozoic strata in the walls of Grand Canyon show a strong peak at 1.1 Ga and numerous Paleozoic grains (numbers 1 and 2 of Fig. 9; Gehrels et al., 2011). The Paleocene–Eocene Music Mountain Formation was derived from Precambrian rocks from the Mogollon highlands to the south of the Colorado Plateau and has peaks of 1.70 and 1.41 Ga (Dickinson et al., 2012).

A test of the age, provenance, and degree of connectedness of Grand Wash Trough basinal sediments was conducted for this study using new detrital-zircon data. These data also help test models for an “old” Grand Canyon of 70–50 Ma (Wernicke, 2011; Flowers and Farley, 2012, 2013) or 17 Ma (Polyak et al. 2008) versus a dominantly “young” western Grand Canyon (Lucchitta, 1966; Karlstrom et al., 2013, 2014). This test relies on the notion that the detrital-zircon makeup of the rocks of Grand Wash Trough that were deposited in the path of the Colorado River from 25 to 6 Ma would be predicted to have appreciable Colorado Plateau–western Grand Canyon–derived detritus in any viable “old” canyon model. New detrital-zircon data are presented in Supplemental Table 4 (see footnote 4).

The 25–17 Ma Rainbow Garden Member of the Horse Spring Formation of the Grand Wash Trough has peaks at 1.70 and 1.1 Ga as well as a syndepositional peak at 20–23 Ma (number 4, Fig. 9). It is interpreted to have been derived from nearby exposures of Lower Mesozoic and Upper Paleozoic strata in adjacent fault blocks with the notable addition of volcanic detritus from the Caliente volcanic field (ca. 19–28 Ma) to the north, but not to have input from the paleo–Colorado River from the east (Lamb et al., 2015). A sample from the 16–14 Ma Thumb Member of the Horse Spring Formation of the Grand Wash Trough (number 5, Fig. 9) has peaks at 1.65, 1.42, 1.05 Ga, and 610 and 415 Ma peaks. An absence of young grains is compatible with models for local derivation in alluvial fans derived from adjacent fault blocks (Beard, 1996).

Four detrital-zircon samples were collected from the rocks of Grand Wash Trough (Bohannon, 1984) that directly underlie the Hualapai Limestone in Grand Wash Trough (number 6, Fig. 9; Kimbrough et al., 2015). The lowest sample (K09-Hual-13) was collected from a sandstone bed that immediately underlies the 13.11 ± 0.08 Ma ash near Pierce Ferry (Fig. 5). This sample is ∼100 m below the basal Hualapai Limestone and is directly in the path of the modern Colorado River. Its detrital-zircon spectrum (90 grains) has a bimodal pattern with peaks of 1.38 Ga and 1.72–1.74 Ga (Supplemental Table 3 [see footnote 3]). Three other samples were collected from the western part of Grand Wash basin, just east of the Wheeler fault and immediately below the basal Hualapai Limestone in the thinner stratigraphic sections west of Airport Point (Fig. 6). K09-Hual 20 was collected in coarse fanglomerate that underlies an ∼10-m-thick Hualapai Limestone section near the divide between Grand Wash and Gregg basins. K09-Hual 21 was collected near the junction of the South Cove and Pearce Ferry roads in red siltstone that underlies a thin Hualapai Limestone section. An additional sample (32606-1) was collected in approximately the same location as K09-Hual 21 (Kimbrough et al., 2015; Supplemental Table 3 [see footnote 3]). The latter three samples are within the upper 20–30 m of the estimated 6–7 Ma top of the Hualapai Limestone in Grand Wash basin.

K09-13 has a well-constrained age (13 Ma), and all of these lower Grand Wash Trough samples (320 grains; number 6 of Fig. 9) represent deposition of red siltstones and fanglomerates in the Grand Wash Trough from 13 to ca. 7 Ma. Collectively, they show a 1.39 Ga peak that is significantly younger than the 1.45–1.4 Ga peak seen in many successions. Pearce (2010) and Pearce et al. (2011) attributed this to derivation from the 1370–1380 Ma Peacock Mountains Granite (Albin et al., 1991) and other ca. 1.37–1.38 Ga granites to the south (black x’s in Fig. 1), but not from western Grand Canyon, as no Paleozoic grains are present. However, a new U-Pb age for the Gold Butte granite near Gold Butte (black star in Fig. 1) gives an age of 1372 ± 12 Ma (Fig. 9C) such that derivation from the local basement to the west seems more likely. The four Grand Wash Trough samples are statistically incongruent (except 326061 and HUAL13, which are virtually identical at P = 0.99, and Hual-13 and Hual-21, which are closely similar at P = 0.37), but collectively, we use them to provide a net composition for the 12–7 Ma pre–Hualapai Limestone rocks of Grand Wash Trough. These spectra are consistent with the Muddy Creek constraint (Lucchitta, 2013) that, between 13 and 7 Ma, limited, or no, far-traveled detritus from the Colorado Plateau or Grand Canyon reached Grand Wash sediments that were deposited in the direct path of the modern Colorado River where it emerges from Grand Canyon. This conclusion is strengthened by the paucity of 1.0–1.2 Ga, Paleozoic, and Mesozoic grains that are common in both the Paleocene and modern Colorado River. Further, these sediments are unlikely to have been derived from the Music Mountain Formation because of the differences in the “twin age peaks.”

We analyzed two samples from a previously unstudied ∼40-m-thick carbonate–red siltstone section in the northern Grand Wash Trough (number 7 of Fig. 9) that has interbedded ash layers of 12–11.7 Ma (Billingsley et al., 2004) and correlates with Hualapai Limestone (Felger and Beard, 2010; Figs. 1 and 5). These samples show detrital-zircon peaks somewhat similar to those of southern Grand Wash Trough, including its young grains, but the northern samples differ in having numerous grains of 1.2–1.0 Ga and Paleozoic age. This spectrum is similar to the Rainbow Garden and Thumb Members (numbers 4 and 5 of Fig. 9) and could reflect derivation from the same sources and perhaps local recycling of these nearby 25–14 Ma units.

Figure 9 (number 8) shows a spectrum from two samples of Hualapai Limestone from the Temple and Gregg basins. The Gregg basin sample is from near the base of the Hualapai Limestone, whereas the Temple basin sample is in the middle of the section and ∼10 m below the 5.97 Ma tephra (Fig. 6). These spectra are similar enough (P = 0.7) to consider them to have been derived from the same population (70% probability). This Temple-Gregg curve generally resembles those in the Rainbow Garden and Thumb Members and the northern Grand Wash Trough, as well as Pliocene–Holocene Colorado River sands (Fig. 9; Kimbrough et al., 2011, 2015), except that the Hualapai Limestone in Gregg and Temple basins lacks age peaks younger than 235 Ma (middle Triassic). Additional analyses are needed to preclude this being a sampling bias, but data so far suggest that no strata younger than Moenkopi Formation were present in the provenance. This is unlike the early Colorado River signature (number 10 of Fig. 9) and suggests that ca. 6 Ma sand in Gregg-Temple basins was derived from Paleozoic and basal Mesozoic strata (Moenkopi Formation contains abundant 245 Ma grains) that overlie Precambrian basement along the flank of the Wheeler Ridge tilt-block between the two basins, thus implying a local provenance west of the Grand Wash Cliffs.

We interpret these data to suggest that, before integration of the different Hualapai Limestone subbasins, sediment was delivered to Grand Wash Trough from the south and west. When Grand Wash spilled over to Gregg basin, progradation from the south was diminished, which allowed local sediment from the closer provenance of Wheeler Ridge to advance southward into Gregg and Temple basins. Alternatively, Gregg and Temple basin were derived from a different provenance than Grand Wash Trough throughout their deposition. The absence of grains younger than 235 Ma in Gregg-Temple basins suggests that none of the Hualapai basins was receiving sand from the Colorado River until post-Hualapai time. The Pliocene Colorado River samples jointly contain 21% grains younger than 215 Ma, and the modern Colorado River sands jointly contain 12% grains younger than 230 Ma. The combined data suggest that Grand Wash basin was internally drained from 25 to ca. 7 Ma. This expands the time frame of the Muddy Creek constraint (Lucchitta, 2013) during which Grand Wash Trough was receiving only local detritus prior to the arrival of the Colorado River 5–6 Ma (Lucchitta, 1966, 2013; Beard, 1996; Pederson, 2008). These detrital-zircon data are less consistent with a proposed proto–western Grand Canyon (Young, 2008; Wernicke, 2011; Flowers and Farley, 2012; Hill and Polyak, 2010, 2014) because such a paleocanyon-paleoriver system should have conveyed Paleozoic detritus from its canyon walls to Grand Wash Trough. Thus, both the geochemical and detrital-zircon data support a model for hydrochemically and depositionally separate, internally drained, fault-bounded basins prior to ca. 6 Ma (Spencer et al., 2013; Faulds et al., 2001) in which the lake/marsh water came dominantly from Colorado Plateau groundwater, but the sediment came from local fault blocks to the south and west.

New Isotopic and Chemostratigraphic Data

Geochemical analyses of Hualapai Limestone were also conducted to evaluate chemostratigraphy within and between individual basins. As noted by other authors (Faulds et al., 2001; Wallace et al., 2005; Spencer et al., 2001), the general absence of evaporite facies between 12 and 6 Ma (except at the base) suggests persistent throughput of water in all three lake/marsh basins during carbonate deposition. Thus, the chemostratigraphy can reveal the nature of waters and temporal trends in water chemistry during this extended 12–6 Ma time interval.

Figure 10 shows up-section variation in 87Sr/86Sr, δ18O, and δ13C. In our best-constrained section (eastern Grand Wash section; black squares), 87Sr/86Sr is fairly consistent (∼0.719) from the base to ∼60 m below the top, and then it decreases up section to values of 0.711. This is similar in the other Grand Wash section (yellow squares). Gregg basin and Temple basin show 87Sr/86Sr values of ∼0.713–0.715 with no marked up-section variation. In δ18O, Grand Wash basins (black, purple, yellow squares) and especially the eastern Grand Wash section (black) show an initial change to more negative (lighter) δ18O starting about −60 m and then a shift toward heavier δ18O in the upper 20 m. Gregg basin (orange circles) is similar and shows an increase in variability and a trend toward heavier δ18O in the upper 50 m, especially in the Hualapai Wash section. Temple basin shows fairly uniform δ18O. There are similar stratigraphic trends in δ13C. Grand Wash basin, and especially in the eastern Grand Wash measured section, shows increasingly negative δ13C values up section starting at ∼60 m from the top. Gregg basin shows this same trend, but perhaps starting ∼80 m from the top. Temple basin also shows lighter δ13C up section in the upper ∼25 m.

The δ18O and δ13C stable isotopic values are fairly constant in the lower 150 m of the Grand Wash basin (perhaps 12 to ca. 8 Ma), whereas isotopic character changes, and higher variability is seen in the upper ∼50 m. We interpret these up-section trends to mean that the Grand Wash basin received progressively more meteoric input (hence lower 87Sr/86Sr) during deposition of the upper ∼50 m of the section, likely corresponding to the time interval from ca. 8 to 6 Ma. The change toward more negative δ18O and more negative δ13C that takes place at about the same stratigraphic level (−60 m) in multiple sections is also consistent with freshening of lake/marsh waters, because lighter δ18O is consistent with addition of high-elevation Colorado Plateau recharge, and lighter δ13C is consistent with more open-system behavior (Talbot, 1990). The change toward heavier δ18O in the uppermost tens of meters of the eastern Grand Wash and Hualapai Wash sections is consistent with increased evaporation due to higher lake levels during integration of the basins.

Figure 11 is a plot of δ13C versus δ18O for the different basins of the Hualapai Limestone. Grand Wash basin δ18O and δ13C values overlap strongly with modern travertines being deposited by Havasu Spring (O’Brien et al., 2006), and they are similar to a western Grand Canyon speleothem sampled near Rampart Cave. The different subbasins of the Hualapai Limestone overlap but have somewhat different stable isotope signatures, as discussed in the following. A new data set from northern Grand Wash Trough (blue symbols in Fig. 11) overlaps with that of southern Grand Wash Trough and suggests they may have been part of the same or similar lake/marsh system formed between the Grand Wash and Wheeler faults.

The 87Sr/86Sr range for Grand Canyon travertines (0.7099–0.7195) is also similar to the Hualapai Limestone (0.7114–0.7195). Stable isotope and 87Sr/86Sr values differ markedly from values of the modern Colorado River near Lake Mead (87Sr/86Sr of 0.7098–0.7108, δ18O ∼−13‰). We conclude that, prior to ca. 8 Ma, Hualapai Limestone precipitated from Colorado Plateau groundwaters similar to those in western Grand Canyon. However, after ca. 8 Ma, Hualapai Limestone deposition reflects increased proportions of higher-elevation Colorado Plateau groundwater. Both source waters were quite distinct from Colorado River water. The marked difference in both 87Sr/86Sr and stable isotope composition between Colorado River water and Hualapai Limestone argues against models for infiltration and piping or syphoning of paleo–Colorado River water through caves and sinkholes as a major step in integration of the Colorado River across the Kaibab uplift and through Grand Canyon (Hunt, 1956; Hill et al., 2008; Hill and Polyak, 2010, 2014; Pederson, 2008).

We conclude that the Hualapai Limestone of Grand Wash subbasin was a spring-fed lacustrine system that, especially in its upper ∼50 m, records evolving regional groundwater flow paths from the Colorado Plateau. Grand Wash, Gregg, and Temple basins have subtle differences, yet they are similar enough that we postulate partial hydrochemical connectedness in the 1–2 Ma leading up to integration of the basins. The variability of Grand Wash and Gregg basin waters suggests poorly mixed waters (spatially and/or temporally), and less variability in Temple basin suggests better mixed waters. By ca. 6 Ma (−50 m in Temple basin section), the similarity of isotopic trends and values in all three basins is compatible with an interpretation that hydrochemical integration had been established. The observed gradual chemostratigraphic changes over the million-year time scale are not well explained by lake spillover models for an abrupt freshening of Lake Hualapai (Scarborough, 2001).

The discovery of meter-scale rimstone dams at the divide between Grand Wash and Gregg basins at western Airport Point (Fig. 12A) provides evidence for meter-scale surface-water flow depth and water turbulence to allow degassing of carbonic waters. The scale of the rimstone dams is similar to structures seen in Havasu Creek (Fig. 12B). Sr data reveal values of ∼0.713–0.715 in the thickest part of the carbonate section, in travertine veins along the Wheeler fault, and for some samples near the top of the section near the pour-over divide. However, samples in adjacent measured sections (Hualapai Wash vs. Spring Wash) show widely disparate values of ∼0.7114 and 0.7195 (respectively). This variation near the top of the Gregg basin sections may record the mixing of spillover waters (from Grand Wash basin; 0.712) with fault- derived spring waters along the Wheeler fault (0.714), and still more radiogenic spring waters from the Wheeler fault (up to 0.7195) that are similar to springs in modern western Grand Canyon. The δ18O and δ13C values in Gregg basin show hydrochemical separation from Grand Wash values in upper parts of the carbonate section, with the heavier δ18O waters interpreted to record increased evaporation as lake surface area increased. More homogenized waters in Temple basin may suggest integration of all three basins by ca. 6 Ma. The combined data are compatible with models proposing partially disconnected hydrologic basins, each characterized by different δ18O-δ13C mixing arrays before 6 Ma, that became integrated across a local spillover path between the two basins near the end of Hualapai Limestone deposition during downward integration (Blackwelder, 1934; House, et al., 2008; Spencer et al., 2013).

Temple basin (Mine Road section of Pearce, 2010) contains the location of the 5.97 Ma ash (Spencer et al., 2001) used for the age of the upper Hualapai Limestone. This section overlies evaporates, and the carbonates may be younger than those of other basins (Roskowski et al., 2010). Chemostratigraphy (Fig. 9) shows little variation in δ18O and an up-section 2‰ increase in δ13C at the top of the section. The combined lower variability of δ18O, δ13C, and 87Sr/86Sr data suggests that this basin was relatively well mixed during Temple basin time, compared to other basins, perhaps suggesting a shorter history of deposition here that was more narrowly linked in time to the integration event.


The Bouse Formation of the Lower Colorado River region (Metzger, 1968) consists of a thin basal carbonate (meter to decimeter scale) that overlies Miocene fanglomerates and in turn is overlain by up to ∼230 m of siliciclastic claystone, mudstone, and sandstone. The basal carbonate is interpreted to record first arrival of Colorado River waters, and the siliciclastic deposits are interpreted to record arrival of Colorado River sediment in a delta that prograded into the lacustrine (Spencer et al., 2008, 2013; Pearthree and House, 2014) and/or perhaps intermittently marine basins (Buising, 1990; McDougall and Miranda Martinez, 2014). The Bouse Formation was deposited in several basins: Las Vegas, Cottonwood, Mojave, Havasu, Bristol, northern Blythe (Parker), southern Blythe, and Yuma basins (Figs. 1 and 2; Spencer et al., 2013). Bouse Limestone is geochemically different than Hualapai Limestone in terms of its 87Sr/86Sr (Fig. 2), and each basin has subtle differences that reflect their own local hydrochemical settings. The age of the Bouse Formation is constrained to be younger than the 5.6–5.84 Ma tuff of Wolverine Creek that underlies it in the Mojave basin, and it is directly dated at 4.83 ± 0.01 Ma by the age of interbedded tuffs in Bristol basin (Harvey, 2014). Dorsey et al. (2007, 2011) concluded that Colorado River sediment first arrived in the Salton Trough at Split Mountain Gorge (location SM, restored in Fig. 1) at 5.3 Ma.

The thin basal Bouse carbonate rests on older fanglomerates across a wide range of elevations within individual basins, and it is interpreted to be a depositional drape that formed as the lakes quickly filled (e.g., House et al., 2008; Spencer et al., 2008, 2013). Elevations of the highest outcrops are interpreted to record lake highstands, which differ from lake to lake across bedrock divides as shown in Figure 1 and the inset to Figure 2: (1) Las Vegas basin (650 m), with divides at upper Black Canyon and modern-day Hoover Dam, (2) Cottonwood basin (560 m) extending to the Pyramid divide (Pearthree and House, 2014), (3) Mojave basin (560 m) extending south to Topoc divide, (4) Havasu basin (360 m) extending to a divide at Parker; (5) Bristol basin (330 m), a NW arm of Blythe basin (Spencer et al., 2013; Miller et al., 2014), and Blythe basin (330 m), which extended to the Chocolate Mountains divide; and (6) Yuma basin, where the top of the inferred Bouse Formation is as much as 1 km below sea level (Olmsted et al., 1973). The Imperial Formation (Fish Creek–Vallecito basin and SM in Fig. 1) was deposited in the northern Gulf of California and is restored by ∼230 km of slip across the San Andreas fault system to its approximate 5–6 Ma position. Similarly, the Ridge basin is restored to its original position near the modern Salton Sea. This offset estimate assumes an average slip rate of 45 km/Ma from 5.3 Ma to the present, as indicated by global positioning system (GPS) velocities and geology-based reconstructions (Oskin and Stock, 2003; Plattner et al., 2007; Darin and Dorsey, 2013).

There are alternatives and cautions to the currently popular interpretation that the elevation of Bouse carbonate outcrops reflects a nearly synchronous bathtub-like carbonate lake coating and that the highest outcrops represent relict lake highstands (e.g., Spencer et al., 2013; Pearthree and House, 2014). First, the age and duration of Bouse deposition in different depocenters are incompletely known. Second, simple calculations based on measurements of possible annual varves in Bouse marl suggest that Bouse deposition could have lasted anywhere from 10 to 300 ka (Homan and Dorsey, 2013). Third, there is well-documented post-Bouse movement on some basin-bounding faults at the divides (Miller et al., 2014; Seixas et al., 2015) and within Bouse basins (Metzger et al., 1973; Miller et al., 2014), suggesting that the original extent and distribution of lakes cannot be precisely reconstructed by contouring the modern elevation of the highest Bouse outcrops (Fig. 1). The amount of post-Bouse offset on faults and related changes in elevation remain undocumented and untested.

In the following sections, we combine results of new analyses with published data for multiple tracers (87Sr/86Sr, [Sr], δ13C, and δ18O) and apply them to Bouse carbonates to test the similarity or difference between basins. The 87Sr/86Sr data are discussed first because of their importance in overturning (Spencer and Patchett, 1997) earlier marine models for parts of the Bouse Formation (Lucchitta, 1979; Buising, 1990).

87Sr/86Sr Isotopic Data

Bouse 87Sr/86Sr ratios range from 0.713 to 0.720 (Fig. 2). Each basin has internal variation of ∼0.001–0.002, suggesting incomplete mixing within each lake. Most Bouse basins are similar in their variability to, for example, the Lake Bonneville system (0.7115–0.7124; Hart et al., 2004). The range and values of 87Sr/86Sr differ subtly but perhaps significantly between basins. For example, Sr isotope ratios in Frenchman Mountain carbonates of the Las Vegas basin (Roskowski et al., 2010) are higher than those of Colorado River water (Gross et al., 2001) but similar to the range seen in Lake Mead spring waters (Supplemental Table 1 [see footnote 1]), and they may have been elevated by nearby radiogenic springs (Fig. 2). Lake-margin travertines near the highstand of Cottonwood basin (new data) are significantly more radiogenic than basin-center carbonates (Roskowski et al., 2010), suggesting spring inputs. Amboy outcrops are more radiogenic than most Bouse carbonates, also suggesting spring inputs in different subbasins.

Where measured sections in individual locations are available (Blythe, Amboy, Frenchman Mountain), there is no observed up-section change in 87Sr/86Sr (thickest section is ∼25 m), suggesting rapid deposition. Where samples are available from elevation transects that may record lake filling of a single basin, Blythe basin shows no relationship between 87Sr/86Sr and elevation (Roskowski et al., 2010), whereas Mojave basin, like Cottonwood basin, shows a slight increase in 87Sr/86Sr at higher elevation. Regionally, the maximum 87Sr/86Sr value in each basin decreases to the south (by ∼0.0003 per basin), suggesting a southward increase of low-87Sr/86Sr waters (e.g., compatible with increased Colorado River input). All values are higher than marine values, although secondary carbonates within the Bouse Formation of southern Blythe basin (Spencer and Patchett, 1997) are intermediate between marine and marl values (Fig. 2). Marine values of 0.709 were measured in shells from a unit correlated as Bouse Formation in the subsurface of the Yuma basin south of the Chocolate Mountains divide (Eberly and Stanley, 1978). With this one exception, all Bouse carbonates are more radiogenic than modern Colorado River water, and most are more radiogenic than the <10 ka Lake Cahuilla tufas that were fed by the ancestral Colorado River (Li et al., 2008a). Thus, we infer that Colorado River water likely mixed with higher-87Sr/86Sr Lake Hualapai water to produce Bouse waters.

Modeling of 87Sr/86Sr and [Sr]

Although 87Sr/86Sr values have figured prominently in studies concluding that the Hualapai Limestone and much of the Bouse Formation are nonmarine (e.g., Spencer and Patchett, 1997), no study has yet considered the importance of strontium concentration [Sr] in mixing models. Figure 13A shows the wide range of 87Sr/86Sr and [Sr] in waters (small symbols) and carbonates (large symbols) of the lower Colorado River region. Colorado River water [Sr] averages ∼1 ppm in eastern Grand Canyon and reaches ∼1.2 ppm in Lake Mead. Grand Canyon springs average ∼0.5 ppm; Lake Mead area springs average ∼3.7 ppm. Evaporative lake waters can have similar [Sr] to the mean of river inflow (e.g., 1.6 and 2.0 ppm for Great Salt Lake and river values, respectively; Hart et al., 2004). However, there is often significant [Sr] enrichment in some saline lakes, for example, ∼24 ppm for modern Salton Sea water, which formed in 1905 when Colorado River water with a mean river input = 2.7 ppm entered the basin (Li et al., 2008b). Geothermal spring inputs also can have very high [Sr]: ∼600 ppm in the Salton Sea (Doe et al., 1966). Carbonates that form from saline lake waters can be significantly enriched in [Sr]; for example, Lake Cahuilla tufas that were deposited between 400 and 20 ka by waters dominated by Colorado River input have a mean of ∼2300 ppm, with shells approaching [Sr] = 6000 ppm (Li et al., 2008b). The modern Gulf of California seawater has [Sr] = 8.4 ppm; San Francisco Bay marine-estuary waters are variable with mean [Sr] = ∼4 ppm (Ingram and DePaolo, 1993). Mean [Sr] for Hualapai Limestone is ∼170 ppm; mean [Sr] for Bouse Formation marls is ∼400 ppm, and mean [Sr] for Bouse shells is ∼900 ppm. Figure 13B shows a north to south decrease in 87Sr/86Sr and increase in [Sr] from Hualapai to Bouse basins.

Figure 14 is the first published effort to model the variability of 87Sr/86Sr values within and between basins and the importance of [Sr] in determining the mixing proportions needed to explain the observed 87Sr/86Sr of different parts of the system. The y axis shows the extreme 87Sr/86Sr variability (0.008) in the Hualapai Limestone and more limited variability in the Bouse Formation (∼0.0015). The x axis is the fraction of the more radiogenic end member for a suite of mixing models. Note that mixing of waters with the same [Sr] (1:1) gives linear mixing, whereas mixing end members with variable [Sr] produces curved models. However, [Sr] is unknown for the waters that deposited past carbonates, so Figure 14 models observed 87Sr/86Sr in terms of different [Sr] assumptions based on end-member values from our data. This provides a speculative, plausible set of mixing models that could be related to downward integration of the early Colorado River, as follows (numbers are keyed to Fig. 14).

(1) Early Hualapai water: We use a spring end-member 87Sr/86Sr value of 0.720 (not our extreme value of 0.735 of western Grand Canyon) to postulate the spring-fed nature of the Grand Wash basin lake/marsh and to explain the fairly consistent values of 0.7193–0.7198 in the early part of the Hualapai Limestone in Grand Wash basin (box 1 of Fig. 14).

(2) Late Hualapai water: Early Hualapai lake/marsh waters are modeled to have mixed with increasing proportions of Colorado Plateau groundwater (CPGW of Fig. 14 with 87Sr/86Sr = 0.712), including an increased component of high-elevation recharge to explain the more negative δ13O in the upper part of the section (Fig. 10). This end member is similar to the mean value of Grand Canyon springs (87Sr/86Sr = 0.7122; [Sr] = 2.6 ppm) and to the high-volume Havasu spring (87Sr/86Sr = 0.7116). This mixing could explain the variability and progressive freshening of 87Sr/86Sr observed in the upper 50 m of the Hualapai Limestone (87Sr/86Sr = 0.715–0.712, represented by 2A to 2B in Fig. 13) in terms of the arrival of increasing amounts of Colorado Plateau groundwater. [Sr] for Hualapai-depositing waters is unknown, but we note that the carbonates deposited in early versus later Hualapai Limestone of Grand Wash basin have similar [Sr] (150 vs. 220 ppm, respectively), so we favor a near 1:1 [Sr] ratio mixing model. In the 1:1 model, 2A values would be explained by ∼40% early Hualapai water mixing with 60% arriving Colorado Plateau groundwater and 2B by ∼10% early Hualapai water mixing with 90% arriving Colorado Plateau groundwater at the very top of the section. Alternatively, if [Sr] values in the early Grand Wash lake/marsh were 4–10 times higher than arriving groundwater, for example, due to evaporative concentration (Li et al., 2008b), a mix of 5%–15% early Hualapai water and 85%–95% Colorado Plateau groundwater would be needed to explain late Hualapai 87Sr/86Sr values of 0.712–0.715 (2C in Fig. 14).

Similar models using arriving Colorado River values of 0.710 (and [Sr] = 1–2 ppm) might also explain the progressive change in the upper 50 m of the section toward decreasing 87Sr/86Sr and increasing [Sr] in the carbonates. However, the gradual change in the top 50 m of the Grand Wash sections suggests a change over several million years, which is more consistent with the interpretation that the arrival of Colorado Plateau groundwater (rather than the Colorado River) led to downward integration from Grand Wash basin to Gregg and Temple basins.

The ranges of 87Sr/86Sr in Gregg and Temple basins (∼0.715–0.0711) are not as large as Grand Wash basin, but they are large enough to suggest incomplete mixing in the lakes/marsh systems; we use a value of 0.712 as a reasonable latest Hualapai Limestone value based on the up-section decrease in 87Sr/86Sr in all basins (Fig. 10). In 2B and 2C, the large proportion of arriving fresher groundwater suggests that groundwater sapping (Pederson, 2001) may have played a significant role in focusing erosion near springs and led to the integration of the Colorado River across upstream drainage divides and through existing paleocanyon segments to instigate downward integration of Hualapai Limestone basins (Karlstrom et al., 2014).

(3) Bouse Water and arrival of the Colorado River (87Sr/86Sr = 0.710): We model the first arrival of far-traveled Colorado River water (as opposed to Colorado Plateau groundwater) to have mixed with late Hualapai water (∼0.712) at 6–5 Ma to produce a few values of 0.712–0.711 in the uppermost Hualapai Limestone and to have initiated downward integration from Temple basin to the Cottonwood basin of the Bouse Formation (Fig. 1; 87Sr/86Sr values of 0.7115–0.7117; 3A in Fig. 14). Strontium concentration in the waters is unknown, but Figure 14 shows that if [Sr] concentrations in Temple basin waters were subequal to [Sr] in the Colorado River (1–2 ppm), addition of a relative proportion of 20%–40% Colorado River water with 60%–80% late Hualapai water could produce the Cottonwood basin values (3A of Fig. 14). More likely, based on modern compositions, [Sr] in the arriving Colorado River water would have been 4–10 times lower than late Hualapai water such that the mixture may have been 20%–45% late Hualapai water mixing with 55%–80% Colorado River (3B of Fig. 14). These models are compatible with the notion that a major water-delivery event of Colorado River water (up to 80% of the total water) initiated downward integration of the Hualapai basin to Bouse subbasins (Fig. 1). Increasing proportions of newly arriving Colorado River water during a sequence of lake filling and overflow events may explain the progressive decrease in maximum 87Sr/86Sr values downstream in the Bouse basins, for example, a mean of ∼0.7108 for Mojave basin (3C of Fig. 14). The models suggest that the range of Bouse 87Sr/86Sr values (0.710–0.7115) can be explained by progressively increased proportions of Colorado River water that mixed with Hualapai basin waters during southward-migrating integration of Bouse basins.

(4) Possible marine mixing in Blythe basin (southern Bouse Formation): In this section, we define “Bouse water” as a mixture of first-arriving Colorado River water and older Hualapai basin waters, in which the relative proportion of the two components was rapidly evolving and highly variable in space and time. Previous workers have stressed that 87Sr/86Sr ratios in carbonates of the southern Blythe basin (0.710–0.711) are significantly more radiogenic than seawater (0.709) and therefore refute a marine influence on Bouse depositional environments (Spencer and Patchett, 1997; Roskowski et al., 2010). To address this question, model 4 in Figure 14 explores the possibility that 87Sr/86Sr ratios in Bouse carbonates could reflect mixing of “Bouse waters” with seawater in a marine-estuary environment (Crossey et al., 2011; McDougall Miranda Martinez, 2014; Miller et al., 2014). If [Sr] in the arriving Bouse water was significantly lower than seawater (∼1–2 ppm for Bouse water, compared to 8 ppm for seawater), a mixture with 85%–95% Bouse water and 5%–15% seawater would produce the observed southern Bouse values of 0.710–0.711 (4A, Fig. 14). If [Sr] of the arriving Bouse water was subequal to that of seawater (∼8 ppm), the observed Sr isotopic ratios could be produced with a mixture of 55%–75% Bouse water and 25%–45% seawater (4B of Fig. 14). If [Sr] in the arriving Bouse water was ∼4 times higher than seawater, a mixture of 25%–40% Bouse water with 60%–75% seawater would produce the Sr isotopic ratios observed in Bouse carbonates (4C of Fig. 14).

Thus, for the southernmost Blythe basin where marine fossils are present, the elevated Bouse 87Sr/86Sr values (0.710–0.711) relative to seawater (0.709) can plausibly be explained with ∼5%–50% seawater mixed with arriving Bouse water (itself a mixture of earliest Colorado River water and older Hualapai basin waters). This modeling indicates that it is neither difficult nor implausible to generate radiogenic values of 0.710 in a Bouse Formation marine-estuary system, contrary to the assertions of previous studies (Spencer and Patchett, 1997; Roskowski et al., 2010; Spencer et al., 2008, 2013). Our mixing models are essentially a “reverse analog” of the modern San Francisco Bay estuary, where waters of the arriving Sacramento and San Joaquin Rivers are much less radiogenic (87Sr/86Sr = 0.7065; [Sr] = 0.13 ppm) and mix with seawater (87Sr/86Sr = 0.7092; [Sr] = 7.9 ppm) in the estuary to produce surface-water values in the estuary that show a wide spread of 87Sr/86Sr values (Fig. 2). This variation was interpreted by Ingram and Sloan (1992) to reflect density stratification as well as seasonal and millennial changes in runoff into the San Francisco Bay. A next step in evaluating this hypothesis for the Bouse Formation is to reconcile the diversity of planktic foraminifers (McDougall and Miranda Martinez, 2014) with possible mixtures of marine and Bouse waters, resulting salinities, and information from other chemical tracers.

In addition to mixing between Bouse waters and seawater, input of endogenic and geothermal groundwaters could also have exerted an important influence on elevated Sr isotopic ratios (Crossey et al., 2011). For example, Salton Sea geothermal sites have 87Sr/86Sr = 0.7110–0.7115, and [Sr] is up to 600–900 ppm (Doe et al., 1966), indicating that geothermal brines could have played a role in mixing. Shallow groundwater in the Imperial Valley is predominantly from the Colorado River, but mixing is seen between these waters and geothermal brines in the Salton Sea region, and seawater mixing with groundwater is documented in coastal springs in Baja, California (Coplen and Kolesar, 1974). In a similar fashion, geothermal groundwaters with high 87Sr/86Sr and [Sr] could have had large leverage in mixing with other water sources in the southern Blythe basin (Crossey et al., 2011).

Overall, our models are compatible with Bouse limestones that formed in a mixture of Lake Hualapai water with voluminous first-arriving Colorado River water (in agreement with Spencer et al., 2013) but not by Colorado River water itself (Patchett and Spencer, 2001). Due to a general lack of data on [Sr] composition of the Pliocene lake waters, the models are permissive of alternate hypotheses (Miller et al., 2014): Bristol basin may have been an arm of Lake Blythe with local spring input that explains higher 87Sr/86Sr, and/or southern Blythe basin could have seen intermittent mixing of marine water (25%–75%) and still deposited carbonates with radiogenic values of 0.710–0.711. Similar mixing models could also be generated that would be compatible with no marine mixing and, instead, an increase in the percentage of Colorado River water through time (Spencer et al., 2013). Thus, our plausibility arguments to try to reconcile Bouse 87Sr/86Sr, marine paleontology (without reliance on avian transport), and stable isotope data (see following) via an estuarine mixing of nonmarine and marine waters are aimed at focusing additional efforts to test alternative models. Additional 87Sr/86Sr analyses, with stratigraphic context, should be focused on outcrops with high diversity of planktic foraminifers below ∼110 m elevation (110 m contour of Fig. 1), as well as subsurface Bouse samples (MacDougal and Miranda Martinez, 2014) to continue to look for any 87Sr/86Sr values between 0.709 and 0.710. Marine values of 0.709 are found in the subsurface of the Yuma basin in rocks correlated with Bouse Formation (Eberly and Stanley [1978],inSpencer and Patchett, 1997) indicating that Gulf of California seawater was at least as far north as latitude 32.7°N in the Colorado River deltaic system at ca. 5–6 Ma (Fig. 1).

Stable Isotopes

Figure 15 shows a compilation of new and published δ18O and δ13C values for Bouse Formation carbonates, in comparison with stable isotopes of the Hualapai Limestone (yellow shading). Bouse samples are highly variable and fall into several groups, with interesting differences between basins. Working downstream, Cottonwood basin is displaced markedly downward (∼4‰ lighter in δ13C and δ18O) relative to Temple basin and most of the Hualapai Limestone but overlaps with samples from the top of Gregg basin. This is compatible with spillover and downward integration of Gregg basin waters to Temple, and then to Cottonwood basin, with considerable evaporation to explain the ∼6‰ heavier δ18O in Cottonwood basin (2–3 of Fig. 15). Mojave basin displays weak covariance of δ18O and δ13C values, which have a similar negative slope to Grand Wash and Gregg basin (interpreted as a spillover trend), and mean values are ∼4‰ lighter in δ18O (3–4 of Fig. 15), compatible with a higher proportion of arriving high-elevation Colorado River recharge during Cottonwood to Mojave spillover (House et al., 2008). Havasu basin values are several per mil heavier in δ13C (4–5 of Fig. 15), compatible with, for example, CO2 degassing during spillover, increased residence time, and/or higher primary organic productivity in Lake Havasu (Fig. 3). Blythe basin values overlap with those of the Cottonwood and Mojave basins, but input source waters (lower left of the covariation trend; Talbot, 1990) are heavier in δ18O (6–7 of Fig. 15), compatible with increased downstream evaporation. Large variation in stable isotopes within basins suggests incomplete mixing. Thus, the overall shift to heavier δ18O and lighter δ13C from Hualapai to the southern Blythe basin (schematically shown in Fig. 3) suggests independent variations between isotopic systems, which is common in open, short-residence-time lakes (Talbot and Kelts, 1990).

Figure 16 shows details of δ18O and δ13C covariation for the southern Blythe basin. Marls form a linear array (R2 = 0.66) with a slope similar to the Lake Cahuilla and Great Salt Lake arrays (Fig. 3) and compatible with evaporation. However, this is anchored near (0,0), which is typical of marine estuaries (e.g., San Francisco Bay) but rare in closed lake systems (e.g., Talbot, 1990; Davis et al., 2012). Fossils in the southern Blythe basin tend to be heavier in δ13C than the marls and up to 10‰ heavier than San Francisco Bay mussels, reflecting organic fractionation and different mixes of waters (Talbot and Kelts, 1990).

Figure 17 is a plot of 87Sr/86Sr versus δ18O. The multitracer approach should be able to produce compatible mixing models for both of these tracers, because both are robust indications of the composition of waters that deposited the carbonates. These data are compatible with the downstream integration sequence described herein for 87Sr/86Sr, and they provide another way to assess the chemostratigraphy trends in Figure 10, thus providing additional insight to possible extent of evaporation in different parts of the system. With numbers keyed to Figure 17, our interpretation is: (1) The most likely source of early Hualapai waters was groundwater from the western Grand Canyon region, where the most radiogenic springs reach values of 87Sr/86Sr = 0.72–0.735, and δ18O is similar to western Grand Canyon travertines (δ18O = −9‰ to −11‰) but not to high-elevation recharge or Colorado River water (δ18O = −12‰). (2) Integration of basins of the Hualapai Limestone system involved progressive freshening of waters due to arrival of Colorado Plateau groundwater with δ18O = −9‰ to −13‰, dissimilar to Colorado River water at Lees Ferry (δ18O = −15‰; 87Sr/86Sr < 0.7105) but similar to mean Grand Canyon and Lake Mead spring water (δ18O = −11.5‰; 87Sr/86Sr = 0.7124). (3) Arrival of the Colorado River abruptly introduced large volumes of low-87Sr/86Sr and negative-δ18O water that drove spillover and downward integration of Temple basin water to Bouse basins. Waters became progressively less radiogenic in 87Sr/86Sr to the south due to higher proportions of river water during downward integration, but they became heavier in δ18O due to increased evaporation. The variability of both parameters suggests incomplete mixing in different basins and/or enough time for local fractionation and mixing processes to affect Bouse basins differently. (4) At the southern end, a near-seawater δ18O = 0‰ could reflect evaporation in a closed lake system (except note that Lake Cahuilla does not get as heavy in δ18O) and/or mixing of seawater in the southern paleolake Blythe (e.g., Ingram et al., 1996; Sampei et al., 2005; Augley et al., 2007). Similarly, 87Sr/86Sr mixing with ∼50% seawater could plausibly explain both the 87Sr/86Sr and near-zero δ18O in the southernmost Blythe basin.

The geochemistry of modern springs in Grand Canyon shows significant variability within a given spring group, with up to 0.02 variation in 87Sr/86Sr, interpreted to reflect differential mixing of deeply derived (endogenic) fluids with meteoric (epigenic) recharge. The ∼2‰–4‰ covarying arrays of δ13C and δ18O are interpreted to reflect seasonal fractionation and groundwater mixing. Western Grand Canyon springs have higher 87Sr/86Sr ratios (0.710–0.735) relative to eastern springs (0.707–0.715). Quaternary (<1 Ma) travertine and speleothems (2–4 Ma) associated with Grand Canyon springs and groundwaters show ranges of 87Sr/86Sr, δ13C, and δ18O in carbonate that overlap with modern water values. These data from young water-carbonate systems provide justification for using carbonate values as a proxy for the water from which they were deposited.

The Hualapai Limestone was deposited from 12 Ma (new ash age) to 6 Ma and is up to 210 m thick. It exhibits variation stratigraphically and between fault-bounded basins. (1) Grand Wash basin has 87Sr/86Sr, δ13C, and δ18O values and ranges similar to western Grand Canyon springs and travertines, and we infer that this basin records a spring-fed lake/marsh system sourced in large part by western Colorado Plateau groundwaters. Up-section chemostratigraphic trends suggest a gradual increase (probably over 1–2 Ma) in groundwaters that delivered recharge from higher elevations. A lack of evaporites (after ca. 13 Ma) indicates continuous flow through Grand Wash and Gregg basins, with a flow path postulated to have been focused southwards along the Grand Wash fault to Red Lake evaporite basin. Meter-scale rimstone dams suggest surface flow and lake spillover were active during downward integration to Gregg basin. (2) Gregg basin shows an up-section increase in δ18O and decrease in 87Sr/86Sr and δ13C, which are interpreted to reflect up-section freshening of waters. (3) Temple basin has less pronounced up-section variation in 87Sr/86Sr, δ13C, and δ18O, suggesting the arrival of combined waters from Hualapai and Gregg basins at ca. 6 Ma. These geochemical data suggest hydrologically distinct basins before ca. 6 Ma. Detrital-zircon spectra from red siltstones underlying the Hualapai Limestone in Grand Wash and Temple basins (Dickinson et al., 2014) also suggest that they were separate basins prior to ca. 6 Ma. We envision a set of fault-segmented internally drained lakes that were integrated ca. 6 Ma to form a single Hualapai lake basin.

We propose that groundwater sapping was an important mechanism for the eventual integration of the Colorado River from the Colorado Plateau to Grand Wash Trough ca. 6 Ma. Groundwater flow and spring discharge can influence erosion and facilitate reorganization of surface-water drainage (Pederson, 2001). Potentially related mechanisms include karst piracy (e.g., transport of Colorado River water through cave passages under the Kaibab uplift; Hill and Polyak, 2014), karst piping (e.g., infiltration of surface water from Hualapai Plateau to Grand Wash through caves; Hunt, 1969; Pederson, 2001), and more general karst connections (where caves connect surface-water drainages beneath a topographic divide; Hill et al., 2008). Our groundwater sapping model does not involve Colorado River surface water moving through caves to explain deposition of the Hualapai Limestone. In contrast to Hunt (1969, 1969), Pederson (2001), and Hill and Polyak (2014), we envision pre–Colorado River movement of groundwater similar to today’s Redwall Muav aquifer. Thermochronologic data suggest that the Kaibab uplift had already been breached 25–15 Ma (Lee et al., 2013; Karlstrom et al., 2014). The ensuing 12–6 Ma deposition of Hualapai Limestone is interpreted here in terms of groundwater flow, first from the western Colorado Plateau (highly radiogenic early Hualapai water) and then from higher-elevation Colorado Plateau sources (late Hualapai water). Deposition of the Hualapai Limestone necessitated high-pCO2, highly radiogenic, and heavy δ18O waters, very different than the modern Colorado River water. The groundwater sapping mechanism, operating over millions of years (12–6 Ma), could have influenced Grand Wash lake throughput and evolving hydrochemistry, focused springs into existing basins and paleocanyons (Karlstrom et al., 2014), and led to drainage reorganization and eventual integration of the Colorado River in its present path 5–6 Ma. Geochemical data seem best explained by a fairly abrupt arrival of voluminous Colorado River that mixed with late Hualapai water and caused downward propagation of lake basins of the Bouse Formation between 5.8 and 4.8 Ma.

Carbonates of the Bouse Formation (5.8–4.8 Ma) were also deposited in several basins separated by bedrock divides (Fig. 1). The basal carbonate is interpreted to record continued downward integration of the Colorado River in water-delivery events resulting from a combination of groundwater sapping and spillover of higher basins. Upper siliciclastic units of the Bouse Formation contain Colorado River sands, and the unit overall is interpreted (in agreement with models of Spencer et al., 2013) to record a downward-integrating lake system that led to final integration of the Colorado River system to the Gulf of California. Hydrochemical tracers suggest several new aspects to this model: (1) Las Vegas has 87Sr/86Sr values compatible with mixing of Hualapai and local spring waters with first-arriving Colorado River water. (2) Cottonwood basin carbonates are more radiogenic (higher 87Sr/86Sr) than carbonates of Las Vegas basin, suggesting addition of radiogenic spring waters similar to modern hot springs below Hoover Dam. The δ18O value in Cottonwood basin is similar to that of Gregg basin but heavier than Temple basin carbonates, suggesting that Cottonwood basin carbonates precipitated from Hualapai-like waters that had seen additional evaporation and/or springs inputs. (3) Mojave basin has less radiogenic 87Sr/86Sr and lighter δ18O than Cottonwood, compatible with increased proportion of arriving Colorado River water. (4) Lake Havasu basin carbonates have δ13C and δ18O values similar to values produced by Hualapai waters mixing with Colorado River waters. (5) Bristol basin is more radiogenic in 87Sr/86Sr than any of the Bouse basins and was likely hydrochemically distinct and/or poorly mixed with the remainder of Blythe basin. (6) Carbonates of the southern Blythe basin have heavier δ18O and lighter δ13C than other Bouse basins and a stable isotope covariation trend similar to that of the modern San Francisco Bay marine estuary. This is compatible with either continued downstream evaporative modification of Bouse waters, and/or intermittent mixing with marine water. The 87Sr/86Sr ratios in Bouse carbonates and shells in all basins are more radiogenic than marine values; hence, Bouse carbonate chemistry is distinctly nonmarine. Bouse carbonates have higher 87Sr/86Sr than the modern or the 10 ka Colorado River water that fed Lake Cahuilla. Each basin shows different δ13C and δ18O mixing trends and a different range of 87Sr/86Sr values relative to other basins. This variability suggests different proportions of mixing among Colorado River water, Hualapai basin water, local spring input, and local meteoric input in hydrologically distinct subbasins, in contrast to patterns expected for well-mixed lakes. The only Bouse paleobasin with possible intermittent marine mixing (based on 87Sr/86Sr, and δ13C and δ18O mixing trends, sedimentology, and paleontology) is the southernmost Blythe basin.

Sr mixing models (87Sr/86Sr vs. [Sr]) can explain the wide range in composition of western Colorado Plateau waters by variable mixing of different end-member waters. Model curves yield mixing proportions needed to explain the observed 87Sr/86Sr and [Sr] values in various parts of the system. Early Hualapai limestone values (e.g., 0.719) are best explained as being sourced from western Colorado Plateau groundwater. Younger Hualapai values (0.712) can be modeled by a groundwater delivery system adding ∼70%–95% higher-elevation Colorado Plateau groundwater similar in composition to Havasu spring groundwater (0.712). This gradual change is recorded in ∼50 m of carbonate deposition (1–2 Ma) and hence is unlikely to be a result of arrival of the Colorado River.

Bouse limestone is modeled by mixing first-arriving Colorado River water with late Hualapai (upper 50 m) lake basin water (0.715). Bouse basins show a downstream decrease in maximum 87Sr/86Sr values and progressively lighter δ18O values, both compatible with increasing proportions of Colorado River water relative to initial Hualapai water as the Colorado River became integrated through progressively lower basins. Significant variability of 87Sr/86Sr and δ18O in Bouse basins suggests incomplete mixing during downward integration. The chemistry of Bristol basin carbonates seems to require mixing of radiogenic spring waters (0.713) with local meteoric water or Colorado River water. For the southernmost Blythe basin, where marine fossils and sedimentary structures are present, the moderately elevated Bouse 87Sr/86Sr values (0.710–0.712) relative to seawater (0.709) can be explained by adding ∼1%–8% endogenic spring/groundwater (e.g., by transport up along faults) to marine water, or by adding ∼50% marine water to an average Bouse value of 0.711. We therefore conclude that elevated 87Sr/86Sr cannot be used to categorically negate the possible influence of marine waters in a marine-estuary environment for the southern Bouse Formation. Nevertheless, our hydrochemical data and modeling generally support a model in which the northern Bouse basins were nonmarine lakes that underwent downward integration due to rapid arrival of Colorado River waters 6–5 Ma. Future work is needed to refine mixing models that include endogenic inputs, surface-water mixing during downward lake basin integration, evolving conditions within each isolated lake system, and the possibility of intermittent mixing with marine waters in the southernmost Blythe basin.

In summary: (1) Hualapai Limestone was deposited by a groundwater system similar to the modern groundwater system such that it records evolving groundwater flow paths on the Colorado Plateau. (2) Increased groundwater input to younger Hualapai basins and groundwater sapping ca. 6–7 Ma may have helped to focus erosion and facilitate integration of the Colorado River through western Grand Canyon at ca. 6 Ma. (3) Grand Wash, Gregg, and Temple basins were hydrochemically and sedimentologically distinct until deposition of uppermost carbonates at ca. 6 Ma; spillover of Grand Wash to Gregg basins is supported by rimstone dams at the spillover point and by the hydrochemistry of δ13C and δ18O data, which is compatible with spillover of the combined Grand Wash and Gregg basins to Temple basin ca. 6 Ma. (4) The hydrochemistry of Bouse carbonates can be explained by spillover of Lake Hualapai basin waters that were mixed with early-arriving Colorado River water plus potential groundwater connections. Each Bouse basin has carbonates with distinct hydrochemistry, suggesting incomplete mixing of waters with multiple sources during downward integration of the early Colorado River. Our synthesis highlights the importance of recognizing multiple water sources and mixing scenarios that could have produced the geochemical and isotopic signatures of carbonates in the Hualapai Limestone and Bouse Formation.

This project benefited from partial financial support from National Science Foundation (NSF) grants EAR-0838575 from the Hydrologic Sciences Program, EAR- 1119629, -1242028, and -1348007 from the Tectonics Program, and DGE-0538396 from the NSF GK-12 program. For tephrochronology, Holly Olson performed the laboratory processing and petrography, Dave Wahl conducted the instrumental and computer analysis, and Elmira Wan interpreted the data. Detrital-zircon samples were run by Mark Pecha and the University of Arizona LaserChron laboratory, and we acknowledge NSF grant EAR-1338583 for support of the Arizona LaserChron Center. W.R. Dickinson helped analyze and interpret the detrital-zircon data. Mark Holland created the detrital-zircon figure, and Jason Ricketts helped with drafting of other figures. We thank Viorel Atudori and the University of New Mexico Stable Isotope Laboratory for stable isotope data. The manuscript benefited from discussions and reviews from W.R Dickinson, David Miller, and an anonymous reviewer, as well as from Guest Associate Editor Andres Aslan.

1Supplemental Table 1. Water and carbonate chemistry from the Grand Canyon–Colorado River region. If you are viewing the PDF of this paper or reading it offline, please visit http://dx.doi.org/10.1130/GES01073.S1 or the full-text article on www.gsapubs.org to view Supplemental Table 1.
2Supplemental Table 2. U.S. Geological Survey tephrochronology report. If you are viewing the PDF of this paper or reading it offline, please visit http://dx.doi.org/10.1130/GES01073.S2 or the full-text article on www.gsapubs.org to view Supplemental Table 2.
3Supplemental Table 3. Data from U-Pb dating of carbonate. If you are viewing the PDF of this paper or reading it offline, please visit http://dx.doi.org/10.1130/GES01073.S3 or the full-text article on www.gsapubs.org to view Supplemental Table 3.
4Supplemental Table 4. U-Pb geochronologic data from detrital zircons and Gold Butte Granite. If you are viewing the PDF of this paper or reading it offline, please visit http://dx.doi.org/10.1130/GES01073.S4 or the full-text article on www.gsapubs.org to view Supplemental Table 4.
5Supplemental Table 5. Location data for geochronology samples and measured sections. If you are viewing the PDF of this paper or reading it offline, please visit http://dx.doi.org/10.1130/GES01073.S5 or the full-text article on www.gsapubs.org to view Supplemental Table 5.