Research ArticleEARTH SCIENCES

Frost for the trees: Did climate increase erosion in unglaciated landscapes during the late Pleistocene?

See allHide authors and affiliations

Science Advances  27 Nov 2015:
Vol. 1, no. 10, e1500715
DOI: 10.1126/sciadv.1500715


Understanding climatic influences on the rates and mechanisms of landscape erosion is an unresolved problem in Earth science that is important for quantifying soil formation rates, sediment and solute fluxes to oceans, and atmospheric CO2 regulation by silicate weathering. Glaciated landscapes record the erosional legacy of glacial intervals through moraine deposits and U-shaped valleys, whereas more widespread unglaciated hillslopes and rivers lack obvious climate signatures, hampering mechanistic theory for how climate sets fluxes and form. Today, periglacial processes in high-elevation settings promote vigorous bedrock-to-regolith conversion and regolith transport, but the extent to which frost processes shaped vast swaths of low- to moderate-elevation terrain during past climate regimes is not well established. By combining a mechanistic frost weathering model with a regional Last Glacial Maximum (LGM) climate reconstruction derived from a paleo-Earth System Model, paleovegetation data, and a paleoerosion archive, we propose that frost-driven sediment production was pervasive during the LGM in our unglaciated Pacific Northwest study site, coincident with a 2.5 times increase in erosion relative to modern rates. Our findings provide a novel framework to quantify how climate modulates sediment production over glacial-interglacial cycles in mid-latitude unglaciated terrain.

  • Climate
  • geomorphology
  • frost processes
  • erosion
  • LGM
  • Paleoclimate
  • soil production
  • frost-cracking


Modern environmental challenges, including soil sustainability, carbon cycling, habitat rehabilitation, and natural hazard assessment, require the establishment and testing of models that describe how past and future climate changes influence bedrock weathering, soil development, and erosion rates (14). The geologic record demonstrates that during the late Cenozoic (characterized by gradual cooling and intensified orbital-scale climate fluctuations), rates of rock exhumation and sedimentation accelerated, particularly in mid-latitude regions (5, 6). Iconic glacial features, such as moraine deposits and U-shaped valleys, record climate-driven changes in erosion rates and processes, but in more widespread and societally relevant unglaciated settings, hillslopes and rivers typically lack obvious climate signatures. This information gap impedes the development of mechanistic theory for how climate modulates landscape fluxes and form.

Abundant paleoenvironmental records demonstrate that unglaciated regions experienced profound climate changes through the late Pleistocene–Holocene transition, but studies quantifying how environmental variables affect erosion and weathering rates in these settings often marginalize or even forego consideration of the role of past climate regimes. Similarly, investigations of orogenic evolution often use Holocene erosion rates for comparison with rates of long-term tectonic forcing; such endeavors may bias the assessment of the relative balance of uplift and denudation if climate modification of erosion is substantial (710). Robust testing of climate-erosion linkages requires data sets with increased temporal continuity to augment terrace sequences and other landforms frequently interpreted as climatic signatures (11). Although changes in precipitation, which regulate river discharge and vegetation coverage, are often invoked to explain terrace formation in unglaciated settings (710), alternative mechanisms have not been well explored (12). In contrast to precipitation-driven erosion linkages, temperature-driven erosional mechanisms that promote rapid denudation in modern alpine settings (13, 14), such as frost cracking and frost creep, may have been more pervasive during the late Pleistocene (15). Unfortunately, few studies have attempted to decipher the previous extent and vigor of these processes [for example, Kirkby (16) and Savi et al. (17)]. Because hillslope and drainage network response time scales are similar to orbital (Milankovitch) cycle time scales (18, 19), the use of a space-for-time substitution framework to infer climate controls on landscape evolution is challenging (11, 2025). Until recently, we have lacked methods and applicable study sites to quantify paleoerosion rates and test models that explicitly connect climate with landscape dynamics. For these reasons, the climatic inheritance of the late Pleistocene in unglaciated landscapes remains poorly constrained.

Organic and inorganic deposits archived in sedimentary basins can record climatic controls on geomorphic processes and rates and enable us to decipher how past climates influence erosion rates, particularly given the advent of isotopic tools such as cosmogenic nuclides. In tectonically active settings, such as the unglaciated Oregon Coast Range (OCR) of the U.S. Pacific Northwest, paleoarchives are rarely preserved. One notable exception is Little Lake in Oregon, a small landslide-dammed remnant of a much larger paleolake (~6 km2 catchment) that has been extensively studied by paleoecologists (Supplementary Materials) (2628). Little Lake is located in the temperate portion of the OCR at latitude 44.2°N, more than 400 km to the south of the maximum extent of the Cordilleran ice sheet (~47.2°N; Fig. 1, inset). Fossil plant communities from Little Lake sediment cores chronicle open canopy forests during Marine Isotope Stage (MIS) 3 [50 to 26 thousand years ago (ka)]; a cold, drier, and nutrient-poor parkland setting during the last glacial (MIS 2; 26 to 13 ka); and a shift toward warm, moister conditions and temperate, closed-canopy Douglas-fir forests at 13 ka with the onset of the modern interglacial interval (MIS 1) (26, 27).

Fig. 1 Map of the Little Lake catchment and sample sites.

Core data location for this study is marked with an asterisk. Previous paleoecology data collected in the fens near the Little Lake outlet are marked with a polygon. Modern stream sample locations are delineated with stars. Map shows only a portion of the larger landslide-dammed paleolake deposit that extends to the east of Triangle Lake. (Inset) Map showing the extent of ice sheets 21 ka, location of study area (dark-gray polygon), and Little Lake, OR. Modern analog ecosystem location is identified with a blue square. Continental extent 21 ka is outlined in black, and modern continental extent is outlined in gray.

To optimize paleoenvironmental reconstructions, previous paleoecology studies used core data from fine-grained deposits along the distal fen of the modern lake (26, 27). Not only is Little Lake well suited for paleoenvironmental analysis, but the catchment and sedimentary system also features the prerequisites necessary and favorable for using detrital 10Be to calculate paleodenudation rates (2932). A suitable paleoerosion field site: (i) spans more than one climatic interval; (ii) is found within a single quartz-rich lithology and a quiescent sedimentary environment, such as a deep sea basin or small lake; (iii) features direct hillslope-to-basin sediment dispersal and deposition pathways; (iv) exhibits ample sand-sized sediment; and (v) contains abundant paleoclimate indicators, such as pollen and fossils, to infer millennial-scale climate variations. To quantify how climate paced erosion rates from the Last Glacial Maximum (LGM) to the present day, we measured detrital 10Be concentrations from river-derived fine sand archived in paleolake sediments (collected via drilling) and modern-day river deposits collected upstream of the paleolake (Fig. 1).

We collected our sediment core upstream of the modern lake with the primary goal of deriving erosion rates from the quartz-rich deposits (Fig. 1 and Supplementary Materials). Specifically, we used 10Be to derive erosion rates and fossil-rich lake deposits within our core to derive sediment ages based on a depth-age model (table S1). As a secondary goal, we sought to refine previous Little Lake paleoclimate reconstructions to better characterize the MIS 2 ecosystem and associated erosion mechanisms. Previous Little Lake glacial climate reconstructions were unable to identify spruce pollen to the species level, leaving unresolved whether the climate at Little Lake was closer to that at the modern Pacific Northwest Cascades, the Olympics, coastal British Columbia, or maritime southeastern Alaska. This lack of resolution implied that MIS 2 mean annual temperatures were 7° to 11°C cooler than present-day temperatures (27). Although our core contained abundant sandy sediment, making it ideal for 10Be analysis (Fig. 1 and Materials and Methods), our site was less optimal for preserving continuous late Pleistocene–Holocene deposition, as the paleolake location transitioned from a depositional setting to a stable and intermittently erosional one by ~20 ka. Nonetheless, we collected, described, and analyzed 60 m of lake sediments extending from 50 ka to the late Pleistocene.

Here, we report sedimentary, geochemical, and paleoecological analyses from the LGM portion of our new Little Lake sediment core (Fig. 2) that enable us to explicitly connect late Pleistocene–Holocene climate change with the rates and mechanisms of surface processes. For the purposes of this study, we focus on sedimentary, geochemical, and erosion rate data relevant to LGM climate reconstructions from ~21 ka. Consistent with predictions of a frost weathering model combined with a down-scaled LGM climate reconstruction informed by Little Lake plant macro fossil data, we propose that frost-dominated erosional processes were pervasive and vigorous during the LGM in unglaciated mid-latitude settings. These processes may have been active in regions commonly assumed to be devoid of geomorphic inheritance from the LGM. Our results document a substantial shift in erosion rates and mechanisms attributable to climate-ecosystem change.

Fig. 2 Compilation of Little Lake core observations and data.

P. sitchensis (Sitka spruce) and A. lasiocarpa (subalpine fir) co-occurrence was observed at 29.5 to 22.6 ka based on depth-age model. Numbers to the left of the tree icons are macrofossil counts for each species at each interval. For simplicity, we only do not include macrofossil counts at 29.5 ka (three P. sitchensis and four A. lasiocarpa). Percentages of clay, sand, and silt in the core are based on visual observations. The entire core sequence consists of millimeter- to centimeter-scale laminated lacustrine deposits, with a significant reduction in fine-scale laminations and an increase in grain size ~26 ka, the start of the last glacial (MIS 2) interval. Median calibrated ages are based on the CLAM model best fit. Although the core extends over 50 ky, for the purpose of this study we only present data that are relevant to the time interval of interest.


Little Lake sediment core

In the sediment core, we observe a transition from finely laminated red, brown, and gray lacustrine clay, silt, and sands to coarse lacustrine blue-gray sand deposits at ~26 ka, coincident with the LGM transition from open canopy forests to a cold and dry parkland setting (Fig. 2) (26, 27). These distinct subangular sand deposits persist throughout the glacial interval and transition to layered, poorly sorted lacustrine deposits during the forested Holocene. The open canopy forest of pollen zone LL-1 (42.5 to 27 ka) and the subalpine parkland of LL-2 (27 to 13 ka) (thousands of calibrated years median before present) correspond to the latter part of MIS 2 and MIS 3, respectively (26, 27).

From an analysis of needle fragments in our sediment core, we refine previous Little Lake pollen-based vegetation reconstructions with macrofossil assemblages that indicate the co-occurrence of Picea sitchensis (Sitka spruce) and Abies lasiocarpa (subalpine fir) from ~29.5 to 22.6 ka (Fig. 2). These two species are rarely found together today at an elevation similar to Little Lake, with the exception of cold parkland settings in southeastern Alaska (latitude 55.6°N) (33). With our identification of the co-occurrence of Sitka spruce and subalpine fir from the end of MIS 3 through MIS 2 [over an interval of 7 thousand years (ky)], we can now better constrain the late Pleistocene climate with a modern analog, a geographically restricted region near Hyder, AK (33). Pollen data from this time period show other species indicative of cool or cold maritime settings, including Tsuga mertensiana (mountain hemlock), Tsuga heterophylla (western hemlock), and herbaceous pollen types (Supplementary Materials) (27). Unlike the forested ecosystems of MIS 3 and MIS 1, subalpine parkland settings are composed of scattered trees or small forested patches amidst open mountain meadows. The presence of these species starting at ~29.5 ka suggests that our 10Be-derived erosion rates from the LGM reflect climate and ecosystem controls on soil production and erosion processes that had persisted for millennia.

LGM frost weathering and increased sediment production and erosion

LGM 10Be-derived erosion rates remained constant at 0.2 ± 0.01 mm/year (mean ± SE; n = 3), which is 2.5 times faster than present-day catchment erosion rates that average 0.08 ± 0.01 mm/year (n = 2) (Fig. 2 and table S2). The modern 10Be-derived erosion rates at Little Lake are similar to modern erosion rates measured throughout the OCR, which range from ~0.05 to 0.14 mm/year (3436). (We have adjusted reported values to account for overestimation of muon production in fast-eroding settings; Supplementary Materials.) In addition, the depth-age model linear regression for the glacial interval (extending from −3.22 to −29.69 m) suggests that sediment accumulation was rapid, at a rate of 6 mm/year (r2 = 0.88) (table S1).

Given the following observations relevant to the LGM at our study site: (i) the presence of tree species characteristic of periglacial settings, (ii) plant species and sediment color indicative of nitrogen-poor regolith (26), (iii) a significant increase in grain size relative to modern size, and (iv) a 2.5 times increase in erosion rates relative to modern rates, we propose that LGM frost weathering and transport processes in our study area contributed to increased sediment production and erosion. Our findings challenge the oft-cited assertion that the OCR approximates steady-state conditions, such that despite stochastic soil production and transport processes, both erosion and landscape form are time-invariant over hillslope response time scales of 50 to 100 ky (3, 18, 34, 37, 38). Rather, our observed 2.5 times late Pleistocene to Holocene erosion rate contrast implies that climate-driven erosion fluctuations can be significant in unglaciated settings despite the lack of an obvious or diagnostic topographic signature. Climate simulations, a maritime setting, and plant macrofossil evidence (27) indicate moderately drier conditions (LGM mean precipitation ~1200 mm/year) with patchy snow cover during the LGM in the OCR, which likely precludes increased discharge and runoff as a means to increase sediment flux (fig. S1). Although we cannot rule out LGM changes in discharge frequency or magnitude because effective runoff predictions from water balance calculations are poorly known, evidence from field observations and paleoclimate simulations suggests that significant swaths of land south of late Pleistocene ice sheets were subject to vigorous frost weathering, frost heave, and solifluction during the LGM (16, 3942).

Paleoclimate simulations and a frost weathering model

Because our paleoarchive analyses indicate colder LGM climes relative to modern climes (Fig. 3), we combined a process-based frost weathering model with paleoclimate reconstructions to predict the spatial extent of periglacial processes in the late Pleistocene. To quantify the intensity of frost weathering at geomorphically relevant scales in our study area, we downscaled (to 90 m) modern PRISM data and CMIP5 (Coupled Model Intercomparison Project Phase 5)/PMIP3 (Palaeoclimate Modelling Intercomparison Project Phase 3) simulations to generate mean monthly near-surface air temperature data (4345). Generally, LGM simulations have a “warm bias” in MAT, primarily generated by insufficient winter cooling (44). We selected the most appropriate CMIP5/PMIP3 models by comparing paleotemperature simulations downscaled to Little Lake with mean monthly temperature data from Hyder, AK, a representative modern environmental analog (33), based on the LGM macrofossil flora documented in our sediment core (Fig. 3). Only two climate models produce annual temperature cycles similar to the modern analog. Although the Community Earth System Models (COSMOS) simulations provide a better fit in the warmer months, the Model for Interdisciplinary Research on Climate (MIROC) simulation more closely conforms to the paleoecological constraints imposed by our Little Lake study site (Fig. 3).

Fig. 3 Annual temperature curves.

Comparison of annual temperature curves based on mean monthly temperature data for modern Little Lake (a; dashed green line); the range of data from downscaled paleosimulations with temperatures above the frost-cracking window (b; gray band with CCSM4, GISS-E2-R, ISPL-CM5A-LR, MPI-ESM-P, and MRI-CGM3 delineated with light-gray lines); COSMOS-ASO (c; solid gray line); Hyder, AK, climate station data (d; dashed blue line); and MIROC-ESM data (e; blue line). All paleomodel data are from CMIP5/PMIP3 simulations (Materials and Methods).

Theoretical and experimental studies highlight the complex physics of premelted films and frost-cracking processes (4650). The interplay of variable permeability, fluctuations in water availability, rock-air temperature coupling, porosity, vegetation, aspect, and cloud cover conspires to control the intensity of segregation ice growth and rock damage (1). Here, we adopt a modified version of a field-verified, frost-cracking model (51), applicable to our regional study domain, as a proxy for frost-driven weathering during the LGM. With this approach, we intend to account for hillslope- and watershed-scale effects. In our model, the amplitude of seasonal air temperature variations and MAT are the dominant factors controlling the vigor of frost weathering at a given location (Materials and Methods). A contour plot of frost-cracking intensity as a function of amplitude and MAT thus provides a simple framework for predicting frost weathering intensity across a range of conditions (Fig. 3 and fig. S2). Zones of high frost-cracking intensity are predicted with MAT values above and below 0°C, where annual variability is high (which generates steep temperature gradients when rock passes through the frost-cracking window). As MAT increases above 0°C (or decreases below 0°C), enhanced frost weathering is predicted, given concurrent increases in the annual temperature amplitude (Fig. 4). Differences in summer insolation and in the degree of ocean and land surface interactions contributed to greater interannual temperature variation during the Pleistocene (52), and thus, frost-cracking potential is amplified during glacial intervals even when MAT does not decrease significantly.

Fig. 4 Contour map of frost-cracking intensity.

Model output of frost-cracking intensity for peak amplitude ranging between 0 and 15 and for MAT ranging between −15°C and 15°C, representing Earth’s range of annual temperature variability. MAT and amplitude values for Little Lake at 21 ka and for the present are delineated by “a” and “b,” respectively.

Late Pleistocene climate reconstruction and frost weathering

In western Oregon, bordered by the Pacific Ocean, elevations range from 0 to 2143 m, and LGM MATs decrease with elevation (Fig. 5, A and B) and latitude (Fig. 5 and fig. S3). The amplitude of annual temperature cycles increases to the east, with distance from the buffering influence of the Pacific Ocean; topography also has an effect, with low amplitudes found in high-elevation areas and with high amplitudes observed in valleys (Fig. 5C). Hence, the highest LGM frost-cracking intensity values are not at higher elevations, as is predicted for the modern climate, but rather in low-elevation settings as a result of greater winter-to-summer temperature fluctuations. Only the southwestern corner of the OCR (42°N to 43°N) does not exhibit frost weathering during the LGM, with the exception of mountain crests, where MATs (~2°C) promote frost cracking despite low temperature amplitude (Figs. 4 and 5). Together, the influences of ocean, topography, and ice-sheet proximity during the LGM combine to generate predicted frost weathering across 90% of our study domain (Fig. 5D). At Little Lake, our model predicts relatively high frost-cracking intensity coincident with our observed LGM acceleration in erosion rate. According to our model, Little Lake frost-cracking intensity values averaged 352 ± 44°C·day (mean ± SD) during the LGM. However, Little Lake is not exceptional, as ~25% of the LGM OCR is predicted to have had frost-cracking intensity values equal to or greater than those at Little Lake and 40% of the OCR had frost-cracking intensity values of 200°C·day or higher (fig. S4). We present these values to highlight the likelihood of frost cracking as a significant LGM geomorphic agent in this unglaciated region. Nonetheless, the physics of frost cracking and rock damage remain under investigation (1, 53), and our model results thus serve as a framework for considering the likelihood of frost cracking rather than a precise map of erosion rate variations. For example, inland of buffering ocean influence, vigorous frost weathering is also predicted across swaths of mid-latitude North America and other continental settings during the LGM, consistent with mapped periglacial features (3941).

Fig. 5 Elevation, MAT, amplitude, and frost-cracking intensity.

Maps showing elevation, MAT, annual amplitude, and frost-cracking intensity 21 ka. All data are overlain on a present-day hillshade map delineating the continental extent of the OCR study area. Downscaled paleodata are adjusted for continental extent 21 ka. White circle delineates the Little Lake study area. (A) Present-day study area elevation. (B) Mean annual temperature 21 ka based on downscaled MIROC model output. (C) Amplitude values (half the temperature range) 21 ka based on mean monthly temperature data. (D) Frost-cracking intensity 21 ka in the OCR.


Our theoretical framework emphasizes the role of frost processes in bedrock weathering because frost-driven sediment transport in cold climes tends to be vigorous and thus unlikely to limit LGM sediment production at our study site (13, 14, 54). Furthermore, our study area features steep hillslope gradients (0.70 ± 0.2, mean ± SD) and narrow, confined valley bottoms, implying rapid sediment transport and minimal opportunity for storage. In modern settings with steep slopes, seasonal frost, and MAT similar to Little Lake’s LGM MAT of ~0°C, solifluction rates average 6 cm/year and can approach ~20 cm/year (14) compared to typical soil creep velocities of ~1 cm/year in forested, temperate settings (54, 55). Therefore, we expect that frost-driven LGM soil transport at Little Lake likely exceeded biogenically driven rates characteristic of the modern closed-canopy forest ecosystem. In addition, our core data (albeit limited in extent) imply more than 25 m of sediment accumulation over the span of the LGM (table S2). When this volume is distributed over the 5 × 106 m2 contributing area of the Little Lake source catchment, it implies exhumation and transport of nearly 5 m of sediment during the LGM. Given the shallow attenuation depth for cosmic rays at our site (~80 cm; Supplementary Materials), we thus conclude that LGM regolith production would have been rapid to accommodate observed paleoerosion rates.

Assuming that hillslope transport processes during the LGM can be represented with a slope-dependent transport model, our measured paleoerosion rates enable us to infer climate-driven changes in model parameters. Soil transport parameters based on a nonlinear model were previously calibrated for OCR hillslopes using Holocene (closed-canopy forest) erosion rates (Materials and Methods) (55). Because the hillslope transport coefficient (K) (L2/T) varies proportionally with erosion rate, our observed 2.5 times increase in erosion suggests that LGM K values increased by 2.5 times relative to modern values. Alternatively, different soil column properties during the LGM may have influenced the critical slope parameter (Sc), which may have contributed to the increased fluxes required to generate our observed late Pleistocene–Holocene erosion rate contrast. Furthermore, for a given soil production relationship, increased erosion typically implies decreasing soil thickness. Although the OCR soil production function estimated by Heimsath et al. (34) for the Holocene closed-canopy forest ecosystem can account for our observed LGM erosion rates, it would require thin soils that may not have been able to accommodate fluxes necessary to produce the observed sedimentary record observed in Little Lake. Instead, we propose that rapid LGM soil production rates may be associated with an alternative parameterization or formulation, perhaps similar to that proposed by Anderson et al. (1).

Our results suggest that increased late Cenozoic erosion in mid-latitude settings (5, 56) may be facilitated by frost-derived sediment fluxes emanating from extensive unglaciated terrain. In addition, our findings encourage a reassessment of conceptual models for how glacial-interglacial climate change manifests in unglaciated settings (11, 15, 57). Rather than precipitation and precipitation-driven controls on vegetation, which are commonly invoked to explain terrace formation or sediment pulses, temperature may be the dominant control on LGM increases in sediment production and transport in many settings (11, 56). By combining mechanistic theory with temperature reconstructions informed by the geologic record and physical geography, we provide a means to resolve the extent to which unglaciated landscapes deviate from modern process mechanisms and rates as a result of climate fluctuations.

Our analyses suggest that frost weathering could be a key control on increased LGM sediment production. Estimates of sediment flux and our paleoclimate reconstruction suggest favorable conditions for frost heave, solifluction, and rapid erosion at Little Lake and similar montane settings. At Little Lake, where measured LGM 10Be-derived erosion rates are 2.5 times greater than Holocene rates, our model predicts moderate to high frost-cracking intensity relative to other regions of our study domain (Fig. 5D). Furthermore, our observed erosion rate contrast is likely a minimum value because of isotopic inheritance that occurs due to transient erosion. Given lower erosion rates during the modern forested climate interval, we speculate that erosion rates were also lower during MIS 3, particularly from ~50 to ~39.5 ka based on paleoclimate reconstructions (Supplementary Materials) (26, 27, 32, 58). If erosion rates are highly sensitive to simulated variations in frost-cracking intensity and sediment transport, our model predicts nonuniform erosion across the OCR during the LGM. This would contrast with the relatively uniform 10Be-derived erosion rates observed for the modern OCR (this study) (3436). The consistency of these modern rates is likely modulated by tree-driven bedrock-to-soil production, colluvial processes, and shallow landsliding associated with closed-canopy coniferous forest stands that characterize the region today (34, 59). Field studies that identify the signature of late Pleistocene frost processes and erosion rates will refine process models and enable us to “map” the vigor of erosional mechanisms in past climes. In contrast to the piedmont region of the eastern United States, where airborne light detection and ranging (LiDAR) data and field observations reveal pervasive periglacial landforms across unglaciated hillslopes (42, 60), such diagnostic features are elusive in the heavily bioturbated and more rapidly eroding steeplands of western Oregon.

Informed by regional paleoecology data, modeling the efficacy of frost processes (tempered by elevation and proximity to large water bodies and ice sheets among other factors) provides a framework to assess the legacy of past climates on modern processes, such as soil development and ecosystem dynamics. For example, our frost weathering model predicts that rock damage can be significant several meters below the surface. This implies that modern critical zone characteristics, particularly shallow bedrock properties such as fracture density, may contain the signature of late Pleistocene processes and thus influence modern processes. Our results suggest that many mountainous regions and broad swaths of continental landscapes proximal to major ice sheets likely experienced accelerated sediment production via frost processes during glacial intervals, inviting a re-evaluation of how past climate regimes imprint modern landscapes.


Experimental design

The combined sedimentology, paleoerosion, and paleovegetation findings in the Little Lake sediment core motivated us to consider a sediment production and transport system unlike the one active today. Our objective was to link geomorphic processes and erosion rates for past climate regimes. More specifically, we sought to incorporate a sophisticated climate perspective that honors how climate plays out across relevant geomorphic spatial scales. Although most studies predict paleotemperatures by using lapse rates and a uniform temperature anomaly (17, 61, 62), we coupled a mechanistic process model that accounted for geographical variability with a regional paleoclimate reconstruction constrained by our analog landscape to test our hypothesis that frost weathering conditions prevailed in the unglaciated OCR 21 ka.

Sediment core

For our study, we sited our core location upstream of the modern lake to satisfy the following criteria and considerations. Given our primary goal of deriving erosion rates from the quartz-rich deposits, we sought to maximize the occurrence of hillslope-derived deposits and thus set our core location in the valley axis proximal to the sediment source area (Fig. 1). The samples required sufficient quartz mass to obtain erosion rates over short (<1000 years) time intervals from sediments with size ranging from 0.25 to 2 mm within the 63.5-mm-diameter core. We used mud rotary drill rigs to collect near-continuous core samples from the paleolake surface through 63 m of sediments to the original valley surface (Fig. 2). After removing the intact 63.5-mm-diameter cores from their 15.2-cm-long metal casings, we split the cores in lengthwise sections and visually described all sections, noting lamination spacing, grain sizes and depositional features (for example, silt with minor sand and vivianite or coarse sand with angular gravels and wood fragments), sediment color, and obvious macrofossils. As part of the drilling process, we preserved the bottom of the core material caught in the “core catcher” and labeled samples from the bottom of the core with “tip.” Although the disturbed samples from the core catcher were unsuitable for stratigraphic description, they contained useable fossils and quartz grains for age and erosion rate analysis. We also noted obvious landslide deposits characterized by poorly sorted coarse angular gravel deposits frequently containing large wood fragments. We avoided any potential landslide deposit for 14C dating and 10Be samples. Using 14C from 17 plant macrofossils, we constructed a depth-age model. To construct the depth-age model, we used a monotonic spline fit to the measured sample depths and the best-modeled median calibrated ages (Little Lake Paleoclimate Archives, table S1), generated with the CLAM model (version 2.2; (63). We also fit a linear regression to the measured depths and median calibrated ages to derive a sediment accumulation rate.

10Be erosion rates

The interpretation of paleoerosion rates requires consideration of cosmogenic nuclide accumulation during hillslope erosion, sediment transport, and sediment deposition. The Little Lake watershed is steep, highly dissected, and subject to colluvial mass wasting processes. In addition, the area features a uniform sandstone lithology and thus is well suited for inferring erosion rates from the cosmogenic nuclide 10Be. The topography suggests minimal potential for sediment storage upstream of the lake deposits (Fig. 1). Topographic analysis of our LiDAR digital elevation model data and sedimentological evidence in the core support a deepwater setting, which would rapidly attenuate postdepositional production by secondary cosmic rays; thus, we have not corrected for nuclide production during or after sediment deposition.

We estimated erosion rates with spatially averaged production rates determined using LiDAR-derived basin hypsometry (table S2). We divided the basin into nine equally spaced elevation bins and calculated production rates by nucleon spallation for each using the CRONUS online calculator (64).The spatially averaged production rate was obtained by weighting each bin according to its fractional area. We calculated production by muon reactions using local production rates following the depth curves of Granger and Smith (65) adjusted for muon cross sections reported by Balco et al. (66) and assuming steady-state erosion. Ignoring muons would lower the calculated erosion rates by 8%.

We interpolated our depth-age model to estimate dates for the paleoerosion rate data. Modern erosion rates were determined from in-stream sediments collected from channels above the influence of the paleolake. All 10Be ratios calibrated to 07KNSTD (67), measured at the PRIME Laboratory (Purdue University). Cosmogenic nuclide analysis requires subtraction of a process blank to account for the small amount of 10Be introduced in the laboratory. At the time of our sample processing, there was an unusually high process blank attributable to contaminated reagents in the chemistry laboratory. Although the blank concentration was high, it was reproducible for blanks before and after each batch of samples, providing confidence in our measurements. This blank subtraction, together with the low 10Be concentrations in the samples, has led to unusually high uncertainties for our erosion rates (10 to 17%). The change in erosion rates that we observe from LGM to modern far exceeds these uncertainties. Table S2 contains more details on the erosion rate calculations.

Erosion rate and hillslope sediment flux calculations

We used a 1-m LiDAR-derived digital elevation model to calculate the contributing watershed area. We used the depths and 14C ages from our sediment core over the LGM interval to calculate sediment thickness and deposit time. For contributing area and hillslope gradient quantifications, we excluded the paleolake deposit area from the topographic analysis. For the nonlinear sediment transport calculation, we usedEmbedded Image(1)where qs is sediment flux, K is diffusivity, ∇z is hillslope gradient, and Sc is the critical hillslope gradient (55).

Frost-cracking model

Annual air temperature variability, which results from local insolation patterns and the strength of ocean-atmosphere and land-surface interactions, produces annual temperature curves significantly better represented by two harmonics rather than a simple sine wave when considering how temperature varies across topography (68, 69). Thus, we integrated harmonic functions for the one-dimensional (1D) surface heat-flux boundary conditions to capture geographic controls on temperature variability. We use an analytical solution for 1D heat conduction (68, 69)Embedded Image(2)where T is daily temperature, z is depth, t is time, MAT is the mean annual temperature calculated from monthly mean near-surface air temperatures, α is the thermal diffusion coefficient, Py is the time period for the curve (we use an annual cycle), and A1, B1, A2, and B2 are coefficients of a Fourier series fit to monthly temperature data extracted from paleoclimate simulations. In the present context, A1 and B1 jointly control the amplitude and phase of the annual temperature cycle (that is, the first harmonic of the monthly temperatures), whereas A2 and B2 control the asymmetry of the cycle. We set α to 0.01 cm2/s. Previous 1D heat conduction models disallowed frost cracking when surface temperatures froze water at ≤0°C (1, 51, 62). However, field and experimental evidence suggests that frost cracking persists even when the surface temperature is below 0°C (13, 70). Thus, we follow previous methods in calculating frost-cracking intensity (1, 51), with a modification that relaxes the restriction against frost cracking when the surface temperature is ≤0°C. We now set the criterion to T > 0°C at the surface or at depth to allow for bidirectional freezing and frost cracking in warm permafrost settings.

Although snow cover can modify surface temperature dramatically, we have chosen to ignore it in part because the maritime influence limits snow accumulation in our study domain. Moreover, the coarse-grained nature of the current model, applied over grid cells of 90 m2 (and thus over slopes with both N and S aspect and across ridges and valleys), implies large variations in the offset between air and surface temperatures due to a variety of effects (for example, aspect, vegetation, wind, and topography). Over landscape scales relevant to geomorphic change, variations in air temperatures are expected to be diagnostic of the surface temperature variations that control temperatures at depth.

In addition to snow, sediment cover will modify predicted frost-cracking intensity values by reducing the thickness of bedrock subject to frost cracking. In calculating a depth-integrated value, we have chosen to ignore sediment cover for the purpose of this study. As OCR hillslopes are steep, frost-generated regolith (unconsolidated material above bedrock or saprolite) transport rates would have been rapid (14), and regolith cover is not likely to limit the pace of sediment production over geomorphic time scales. Modern sediment cover in the OCR is variable as a result of the stochastic nature of biotic disruption (34, 60); however, local stochastic changes in soil depth (for example, induced by bioturbation) equilibrate rapidly (71, 72). Modern soils are generally thin on noses and side slopes, with average depths of less than 0.5 m. On the basis of estimated LGM sediment fluxes and solifluction rates in similar settings (14), we expect that sediment depths would have been, at most, equal or, more likely, thinner during the LGM. As we are not using the frost-cracking index to calculate absolute values but rather to demonstrate geographic variation and the overall expanse of frost-cracking intensity, we have taken the most direct approach to assessing frost weathering across our study domain.

As our paleoclimate simulations depend on downscaled modern PRISM climate data (PRISM Climate Group, Oregon State University, (see Paleoclimate Reconstruction), which record near-surface air temperatures, we use air surface rather than surface (skin) temperatures, as the alternative would require a land-surface model that can predict deviations between air and surface temperatures. We chose to avoid generating secondary effects that are sensitive to poorly constrained and heterogeneous factors such as snow cover, fracture density, vegetation, cloud cover, and permeability.

Paleoclimate reconstruction

We downloaded climate model output for the LGM (lgm) and modern control (piControl) simulations from the CMIP5 Web sites (terms in italics are the database field names). The CMIP5 models selected included CCSM4, COSMOS-ASO, GISS-E2-R, ISPL-CM5A-LR, MIROC-ESM, MPI-ESM-P, and MRI-CGM3. We used the last 100 years of monthly data for near-surface air temperature (tas). Because these models have a very coarse grid resolution (0.5 to 3°), we relied on interpolated downscaled modern data to generate geomorphically relevant downscaled paleoclimate reconstructions. We obtained 30-year monthly mean maximum and minimum temperature data (tmin and tmax) from the 30-s (~800 m) PRISM climate data set. From these, we generated monthly average temperatures. The monthly average temperature data were then used to calculate local topographic lapse rates for each month of the year (73). The lapse rates and CMIP5 data were both (bilinearly) interpolated onto a 90-m digital elevation model, applying the lapse rates to generate elevationally adjusted monthly temperatures on the grid (73). Using the seven CMIP5/PMIP3 simulations and the ensemble average, we generated monthly lgm minus piControl long-term mean differences (or “anomalies”). We then added these long-term mean differences to the present-day 90-m grid to produce a downscaled 90-m map of simulated LGM monthly temperatures.

We compared model output at location 44.18°N, 123.56°W, a representative location in the Little Lake watershed, midway between the valley floor and ridgetop elevations, at 400 m, to present-day mean monthly temperature data for a similar representative location (based on vegetation) near Hyder, AK, our analog landscape, at 400 m (55.75°N, 130.49°W), using ClimateWNA (ClimateWNA, University of British Columbia, (74).

Model implementation

We calculated LGM frost weathering intensity across western Oregon by inputting the mean monthly temperature data derived from the downscaled MIROC paleoclimate simulations into our frost-cracking model (Figs. 3 and 5B) and by calculating frost-cracking intensity for each 90-m grid cell. Using the downscaled paleoclimate simulations, we calculated MAT and the harmonic shape coefficients by fitting a two-term Fourier series to the monthly temperature data and the median Julian day within each month for each grid node within our study domain. We then used Eq. 2 to calculate frost-cracking intensity. For computational efficiency, we subsampled the 90-m grid at a 270-m interval.


Supplementary material for this article is available at

Materials and Methods

Fig. S1. Precipitation anomaly map.

Fig. S2. Additional information on frost-cracking intensity phase map (Fig. 3).

Fig. S3. Temperature anomaly map.

Fig. S4. Frost-cracking intensity distributions for Little Lake and the OCR.

Table S1. 14C data and depths used in depth-age model.

Table S2. Cosmogenic nuclide data and calculated erosion rates.

References (7580)

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank K. Izumi for assisting with the preparation of the paleoclimate simulations, H. Grist for assistance with core and 10Be processing, and A. Noble and J. Noble for drilling and sampling access on their property. Helpful comments from the assistant editor K. Whipple and reviews from A. Heimsath, K. Ferrier, and one anonymous reviewer significantly improved the manuscript. Funding: J.A.M., J.J.R., D.G.G., and D.E.G. were supported by NSF EAR-0952186. Author contributions: J.A.M. directed the drilling operations, oversaw all core extraction core logging and preliminary 10Be processing, and developed the reformulated frost-cracking equation after consultation with P.J.B., whereas D.G.G. identified all macrofossils from the core. J.A.M. and D.G.G. cleaned all fossils used for 14C analysis. D.G.G. developed the depth-age model with assistance from J.A.M., whereas S.J.P. generated downscaled elevation-adjusted modern monthly temperature curves for the study area. S.P. and P.J.B. produced maps of LGM monthly temperatures for the study area. T.C.H. provided the original model code for generating depth-integrated frost cracking, modified by J.A.M. Theoretical and model implications of frost weathering effectiveness were analyzed and described by A.W.R., whereas D.E.G. supervised all 10Be processing and analysis. P.J.B. and T.C.H. were involved in the study design. J.A.M. and J.J.R. designed the study and analyzed the data. J.A.M. and J.J.R. wrote the paper. All authors discussed the results and commented on the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.
View Abstract

Navigate This Article