Annual Findings

Our findings are summarized annually in reports provided to NSF.  Each individual report for years 2008 to 2016 is available for viewing here:  Reports/Proposals.

Below is a select list of our findings provided to the scientific community through journal publication.  Each article is summarized with an abstract and supporting figure.   The selection will be periodically updated.   

Ecology REU student Shelley Pickett discusses a research protocol with Chris Duffy during the summer of 2010.

2017

Understanding watershed hydrogeochemistry: 2. Synchronized hydrological and geochemical processes drive stream chemostatic behavior

A schematic representation of processes in different modules (different colors) of RT-Flux-PIHM. The model allows systematic
understanding of coupled processes at the grid and watershed scales. Note that three types of flow contribute to the stream discharge
QD: overland flow (surface runoff, Qs), subsurface lateral flow in soil (QL), and deep groundwater flow (QG).

Why do solute concentrations in streams remain largely constant while discharge varies by orders of magnitude? We used a new hydrological land surface and reactive transport code, RT-Flux-PIHM, to understand this long-standing puzzle. We focus on the nonreactive chloride (Cl) and reactive magnesium (Mg) in the Susquehanna Shale Hills Critical Zone Observatory (SSHCZO). Simulation results show that stream discharge comes from surface runoff (Qs), soil lateral flow (QL), and deeper groundwater (QG), with QL contributing >70%. In the summer, when high evapotranspiration dries up and disconnects most of the watershed from the stream, Cl is trapped along planar hillslopes. Successive rainfalls connect the watershed and mobilize trapped Cl, which counteracts dilution effects brought about by high water storage (Vw) and maintains chemostasis. Similarly, the synchronous response of clay dissolution rates (Mg source) to hydrological conditions, maintained largely by a relatively constant ratio between ‘‘wetted’’ mineral surface area Aw and Vw, controls Mg chemostatic behavior. Sensitivity analysis indicates that cation exchange plays a secondary role in determining chemostasis compared to clay dissolution, although it does store an order-of-magnitude more Mg on exchange sites than soil water. Model simulations indicate that dilution (concentration decrease with increasing discharge) occurs only when mass influxes from soil lateral flow are negligible (e.g., via having low clay surface area) so that stream discharge is dominated by relatively constant mass fluxes from deep groundwater that are unresponsive to surface hydrological conditions.

Water Resources Research, doi:10.1002/2016WR018935

Understanding watershed hydrogeochemistry: 1. Development of RT-Flux-PIHM

A schematic representation of processes in RT-Flux-PIHM. Flux-PIHM simulates the hydrological and land-surface dynamics (precipitation,
canopy interception, infiltration, recharge, overland flow, subsurface lateral flow, river flow, and surface energy balance) at the watershed scale using the semidiscrete finite volume method. 

Model development in hydrology and geochemistry has been advancing separately with limited integration. We developed a watershed hydrogeochemical code RT-Flux-PIHM to understand complex interactions between hydrological processes (PIHM), land-surface processes (FLUX—Noah Land Surface Model), and multicomponent subsurface reactive transport (RT). The RT module simulates geochemical processes including aqueous complexation, surface complexation, mineral dissolution and precipitation, and cation exchange. The RT module is verified against the widely used reactive transport code CrunchFlow. The code uses semidiscrete finite volume method and irregular gridding and offers data harvesting capabilities from national databases. The application of RT-Flux-PIHM is demonstrated in the Susquehanna Shale Hills Critical Zone Observatory (SSHCZO). We aim to understand key processes that govern hydrogeochemical dynamics of the nonreactive chloride and reactive magnesium. Simulation results indicate that watershed characteristics, in particular topography, dictate the spatial distributions of water content and soil dissolution rates. Ion exchange provides buffering capacities and leads to a hysteresis loop of concentration and discharge relationship of magnesium, which differs from the open hysteresis of chloride.  RT-Flux-PIHM offers physics-based modeling capabilities to integrate the vast amount of water and chemistry data that have now become available, to differentiate the relative importance of competing processes, and to test hypotheses at the interface of hydrology and geochemistry.

Water Resources Research, doi:10.1002/2016WR018934

Combined soil-terrain stratification for characterizing catchment-scale soil moisture variation

Maps of terrain attributes derived from the digital elevation model created from 1-m resolution LiDAR: (a) local slope [m m− 1], (b) depth to bedrock (DtB) [cm], (c) vertical distance to stream (VDS) [m], (d) natural log-transformed upslope contributing area (UCA) [ln (m2)], (e) topographic wetness index (TWI) [−], and (f) surface curvature [−].

Soil properties and terrain characteristics influence spatiotemporal patterns of soil moisture across a watershed. To improve the predictive power of landscape hydrologic models, it is essential to consider both soil and terrain attributes when stratifying a catchment into similar hydrologic functional units. In this study, we developed and validated a new catchment-scale stratification scheme for the Shale Hills watershed by combining soil and terrain attributes in an attempt to delineate soil-landscape units with similar soil moisture dynamics. Terrain was combined with soils information by first using a Random Forest supervised classification algorithm to predict a detailed soil map using 47 field soil samples and terrain variables derived from 1-m LiDAR. A slope class map generated from the LiDAR-derived digital elevation model (DEM) was overlaid on the predicted soil map to delineate areas of similar slope value across the catchment. We compared the performance of this new stratification scheme with two classical stratification schemes, a soil map developed from detailed field survey and a landform unit map based on the DEM, for estimating soil moisture time-series across the forested watershed. The combined soil-terrain method outperformed classical stratification schemes in estimating soil moisture time-series over a 4-year period. Our results demonstrate that combining soil and terrain attributes can help improve the stratification of a catchment into similar soil hydrologic functional units, which is valuable to distributed hydrology modeling and other applications.

Geoderma, doi: 10.1016/j.geoderma.2016.09.031

2016

Rapid tree water transport and residence times in a Pennsylvania catchment

Vessel size distribution for one branch sample from each tree species with counts of vessels in each size category. Between five and nine cross sections were analysed per branch. Total number of vessels analysed:  QUPR, n = 463; QURU, n = 315; CATO, n = 268, ACSA, n = 573.

Trees are responsible for the majority of precipitation recycling over land and can affect soil water storage, stream flow, and ground water recharge. Historically, water has not often been limiting in eastern U.S. forests. As a result, little work has been done to understand the timing of water use by vegetation in these systems. We used deuterium tracer, sap flux, and anatomical techniques to study tree water transport on a forested ridge top in central Pennsylvania. Three trees of each of the dominant ring-porous species, Carya tomentosa (mockernut hickory), Quercus prinus (chestnut oak), and Quercus rubra (red oak), and the diffuse-porous Acer saccharum (sugar maple), were studied. We hypothesized that tracer velocity would be higher in the ring-porous species because of their greater vessel diameters and water transport efficiency, but that tracer residence time would be largely dependent on tree size.  The tracer traveled at velocities of 1 to 18 m d-1 with maximum deuterium concentration in the crowns of the study trees being reached between 1 and 12 days after injection. Tracer residence time ranged from about 5 to 22 days with no evidence of longer residence times in larger trees.  There was also no evidence of a relationship between tracer velocity and calculated xylem specific conductivity, which varied by nearly an order of magnitude between species. However, the soil-to-leaf driving force for water transport may be a strong determinant of tracer velocity across species, and shows promise as a proxy for sap velocity in hydrologic modeling applications.

Ecohydrology, doi: 10.1002/eco.1747

Mineralogical Transformations and Soil Development in Shale across a Latitudinal Climosequence

Bulk soil mineralogy as a function of depth across the climosequence. Vermic/HIV includes vermiculite and hydroxy-interlayered vermiculite minerals.

To investigate factors controlling soil formation, we established a climosequence as part of the Susquehanna-Shale Hills Critical Zone Observatory (SSHCZO) in central Pennsylvania, USA. Sites were located on organic matter-poor, iron-rich Silurian-aged shale in Wales, Pennsylvania, Virginia, Tennessee, Alabama, and Puerto Rico, although this last site is underlain by a younger shale. Across the climosequence, mean annual temperature (MAT) increases from 7 to 24°C and mean annual precipitation (MAP) ranges from 100 to 250 cm. Variations in soil characteristics along the climosequence, including depth, morphology, particle-size distribution, geochemistry, and bulk and clay mineralogy, were characterized to investigate the role of climate in controlling mineral transformations and soil formation. Overall, soil horizonation, depth, clay content, and chemical depletion increase with increasing temperature and precipitation, consistent with enhanced soil development and weathering processes in warmer and wetter locations. Secondary minerals are present at higher concentrations at the warmest sites of the climosequence; kaolinite increases from <5% at northern sites in Wales and Pennsylvania to 30% in Puerto Rico. The deepest observed weathering reaction is plagioclase feldspar dissolution followed by the transformation of chlorite and illite to vermiculite and hydroxy-interlayered vermiculite. Plagioclase, although constituting <12% of the initial shale mineralogy, may be the profile initiating reaction that begins shale bedrock transformation to weathered regolith. Weathering of the more abundant chlorite and illite minerals (~70% of initial mineralogy), however, are more likely controlling regolith thickness. Climate appears to play a central role in driving soil formation and mineral weathering reactions across the climosequence.

Soil Science Society of America Journal, doi:10.2136/sssaj2015.05.0202

Fully-coupled hydrologic processes for modeling landscape evolution

Definition sketch of 3-D control volume on hillslope, where z(x, y, t) ¼ ground surface elevation (m), e(x, y, t) ¼ bedrock-regolith interface elevation (m), h(x, y, t) ¼ regolith thickness in vertical (m) (slope-normal thickness, H, is equal to h cos q, where q is the landscape surface slope), U(x, y, t) ¼ rock uplift rate (m yr1), Bw(x, y, t) ¼ bedrock weathering rate (m yr1) in the vertical direction, qs ¼ the surface flux of regolith sediment by overland flow (m2 yr1), and qc ¼ the lateral volumetric regolith flux by creep and tree throw (m2 yr1) entering through the sides of the control volume. The bottom of the stream may consist of alluvial deposits on bedrock.

Although current landscape evolution models can predict landscapes with specific concave-convex slopes, regolith thicknesses, drainage densities and relief, these models rarely include realistic groundwater and overland flows, and channel-hillslope interactions. To overcome the potential drawbacks, this study couples hydrologic processes with hillslope and channel sediment transport processes to form a new hydrologic-morphodynamic model (LE-PIHM) for regolith formation and landscape evolution. Two scenarios with and without groundwater flow are presented to demonstrate the importance of this coupling. Comparison of the steady state landforms indicates that hillslopes are steeper and relief is higher with groundwater flow. The sensitivity of the solution to mesh geometry is tested and it is shown that model simulations maintain the characteristic features of a landscape over a reasonable range of maximum area and minimum interior angle. To predict long-term landscape change, a morphological acceleration technique is presented and a method for choosing an optimal morphological scale factor is introduced.

Environmental Modelling & Software, doi: 10.1016/j.envsoft.2016.04.014

Designing a suite of measurements to understand the critical zone

Map of Shavers Creek watershed, highlighting (a) topography derived from airborne lidar, (b) geology (Berg et al., 1980), and (c) land use (Homer et al., 2015). In moving from “measure everything everywhere” (our paradigm in the 8 ha Shale Hills catchment; SH) to “measure only what is needed” in the Shavers Creek watershed (164 km2), we chose to investigate two new first-order subcatchments: a
forested sandstone site (along Garner Run, marked GR) and an agricultural calcareous shale site (to be determined). In addition, three sites on Shavers Creek have been chosen as stream discharge and chemistry monitoring sites (marked SCAL – Shavers Creek above lake; SCBL – Shavers Creek below lake; and SCO – Shavers Creek outlet). 

Many scientists have begun to refer to the earth surface environment from the upper canopy to the depths of bedrock as the critical zone (CZ). Identification of the CZ as an integral object worthy of study implicitly posits that the study of the whole earth surface will provide benefits that do not arise when studying the individual parts. To study the CZ, however, requires prioritizing among the measurements that can be made – and we do not generally agree on the priorities. Currently, the Susquehanna Shale Hills Critical Zone Observatory (SSHCZO) is expanding from a small original focus area (0.08 km2, Shale Hills catchment), to a larger watershed (164 km2, Shaver's Creek watershed) and is grappling with the prioritization. This effort is an expansion from a monolithologic first-order forested catchment to a watershed that encompasses several lithologies (shale, sandstone, limestone) and land use types (forest, agriculture). The goal of the project remains the same:  to understand water, energy, gas, solute, and sediment (WEGSS) fluxes that are occurring today in the context of the record of those fluxes over geologic time as recorded in soil profiles, the sedimentary record, and landscape morphology.

Given the small size of the Shale Hills catchment, the original design incorporated measurement of as many parameters as possible at high temporal and spatial density. In the larger Shaver's Creek watershed, however, we must focus the measurements. We describe a strategy of data collection and modeling based on a geomorphological and land use framework that builds on the hillslope as the basic unit. Interpolation and extrapolation beyond specific sites relies on geophysical surveying, remote sensing, geomorphic analysis, the study of natural integrators such as streams, groundwaters or air, and application of a suite of CZ models. We hypothesize that measurements of a few important variables at strategic locations within a geomorphological framework will allow development of predictive models of CZ behavior. In turn, the measurements and models will reveal how the larger watershed will respond to perturbations both now and into the future.

Earth Surface Dynamics, doi:10.5194/esurf-4-211-2016

2015

Reliance on shallow soil water in a mixed-hardwood forest in central Pennsylvania

Depiction of Shale Hills catchment and related plant and soil traits with relative slope position and relative elevation.

We investigated depth of water uptake of trees on shale-derived soils in order to assess the importance of roots over a meter deep as a driver of water use in a central Pennsylvania catchment. This information is not only needed to improve basic understanding of water use in these forests but also to improve descriptions of root function at depth in hydrologic process models. The study took place at the Susquehanna Shale Hills Critical Zone Observatory in central Pennsylvania. We asked two main questions: (i) Do trees in a mixed-hardwood, humid temperate forest in a central Pennsylvania catchment rely on deep roots for water during dry portions of the growing season? (ii) What is the role of tree genus, size, soil depth and hillslope position on the depth of water extraction by trees? Based on multiple lines of evidence, including stable isotope natural abundance, sap flux and soil moisture depletion patterns with depth, the majority of water uptake during the dry part of the growing season occurred, on average, at less than ∼60 cm soil depth throughout the catchment. While there were some trends in depth of water uptake related to genus, tree size and soil depth, water uptake was more uniformly shallow than we expected. Our results suggest that these types of forests may rely considerably on water sources that are quite shallow, even in the drier parts of the growing season.

Tree Physiology, doi: 10.1093/treephys/tpv113

Frequency and Control of Subsurface Preferential Flow: From Pedon to Catchment Scales

Measured (M) and interpolated (I) preferential flow (PF) occurrence frequency based on the data from the period May 2011 to June 2012 for all 35 sites, with all the available data in each soil profile used in estimation (with the bar height representing the size of PF frequency, which ranged from <1% at Site 31 to 70.1% at Site 67). The interpolation map was generated using regression kriging.

Quantitative assessment of frequency and control of preferential flow (PF) across the landscape has been largely lacking. Previous work evaluated PF occurrence at 10 sites along a hillslope in the Shale Hills Catchment using soil moisture response to 175 precipitation events. We expanded the analysis to include (i) 237 additional events to test the temporal consistency and predictability of PF occurrence and (ii) 25 additional sites to upscale to the entire catchment. The results showed considerable temporal consistence in both frequency and main controls of PF at the hillslope scale, attributed largely to statistical stability of precipitation patterns during the 6.5-yr monitoring and relatively stable subsurface PF paths. Generally, PF tended to occur more often in response to intense rainfalls and favored conditions at dry hilltop or wet valley sites. When upscaling to the catchment, topographic controls became more evident, leading to the identification of a hidden subsurface PF network. Higher frequency of PF occurred at the hilltop (average 46%) and the valley floor (average 41%), while the overall average frequency for swales was 26% and that for planar and convex hillslopes was 18%. Soil-terrain attributes provided a limited estimation (R2 = 0.43–0.48) of PF occurrence, suggesting complexities involved in PF dynamics. This study confirmed that the initiation and persistence of PF were controlled by interactions among landforms, soils, initial moisture conditions, precipitation, and seasons. Further investigations of these key controls can lead to improved understanding and modeling of PF from pedon to catchment scales.

Soil Science Society of America Journal, doi: 10.2136/sssaj2014.08.0330

Cyber-Innovated Watershed Research at the Shale Hills Critical Zone Observatory

PIHM structure and processes. The upper subfigure shows the hydrological processes of PIHM at a cross section of a watershed. The lower subfigure shows the spatial structure of PIHM. Blue lines represent the stream channels, and triangles represent the catchment domain.

Cyberinfrastructure is enabling ever more integrative and transformative science. Technological advances in cyberinfrastructure have allowed deeper understanding of watershed hydrology by improved integration of data, information, and models. The synthesis of all sources of hydrologic variables (historical, real time, future scenarios, observed, and modeled) requires advanced data acquisition, data storage, data management, data integration, data mining, and data visualization. In this context, cyber-innovated hydrologic research was implemented to carry out watershed-based historical climate simulations at the Shale Hills Critical Zone Observatory. The simulations were based on the assimilation of data from a hydrologic monitoring network into a multiphysics hydrologic model (the Penn State Integrated
Hydrology Model). We documented workflows for the model application and applied the model to short-time hyporheic exchange flow study and long-term climate scenario analysis. The effort reported herein demonstrates that advances in cyberscience allows innovative research that improves our ability to access and share data; to allow collective development of science hypotheses; and to support building models via team participation. We simplified communications between model developers and community scientists, software professionals, students, and decision makers, which in the long term will improve the utilization of hydrologic models for science and societal applications.

IEEE SYSTEMS JOURNAL, doi: 10.1109/JSYST.2015.2484219

Natural and anthropogenic processes contributing to metal enrichment in surface soils of central Pennsylvania

Conceptual diagram of metal cycling in soils at SSHO

Metals in soils may positively or negatively affect plants as well as soil micro-organisms and mesofauna, depending on their abundance and bioavailability. Atmospheric deposition and biological uplift commonly result in metal enrichment in surface soils, but the relative importance of these processes is not always resolved. Here, we used an integrated approach to study the cycling of phosphorus and a suite of metals from the soil to the canopy (and back) in a temperate watershed. The behavior of elements in these surface soils fell into three categories. First, Al, Fe, V, Co, and Cr showed little to no enrichment in the top soil layers, and their concentrations were determined primarily by soil production fluxes with little influence of either atmospheric inputs or biological activity. Second, P, Cu, Zn and Cd were moderately enriched in surface soils due to a combination of atmospheric deposition and biological uplift. Among the metals we studied, Cu, Zn and Cd concentrations in surface soils were the most sensitive to changes in atmospheric deposition fluxes. Finally, Mo and Mn showed strong enrichment in the top soil layer that could not be explained strictly by either current atmospheric deposition or biological recycling processes, but may reflect both their unique chemistry and remnants of past anthropogenic fluxes. Mn has a long residence time in the soil partly due to intense biological uplift that retains Mn in the top soil layer. Mo, in spite of the high solubility of molybdate, remains in the soil because of strong binding to natural organic matter. This study demonstrates the need to consider simultaneously the vegetation and the soils to understand elemental distribution within soil profiles as well as cycling within watersheds.

Biogeochemistry, doi: 10.1007/s10533-015-0068-5

2014

Designing a Suite of Models to Explore Critical Zone Function

The Critical Zone (CZ) incorporates all aspects of the earth's environment from the vegetation canopy to the bottom of groundwater. CZ researchers target processes that cross timescales from that of water fluxes (milliseconds to decades) to that of the evolution of landforms (thousands to tens of millions of years). Conceptual and numerical models are used to investigate the important fluxes: water, energy, solutes, carbon, nitrogen, and sediments. Depending upon the questions addressed, these models must calculate the distribution of landforms, regolith structure and chemistry, biota, and the chemistry of water, solutes, sediments, and soil atmospheres. No single model can accomplish all these objectives. We are designing a group of models or model capabilities to explore the CZ and testing them at the Susquehanna Shale Hills CZ Observatory. To examine processes over different timescales, we establish the core hydrologic fluxes using the Penn State Integrated Hydrologic Model (PIHM) – and then augment PIHM with simulation modules. For example, most land-atmosphere models currently do not incorporate an accurate representation of the geologic subsurface. We are exploring what aspects of subsurface structure must be accurately modelled to simulate water, carbon, energy, and sediment fluxes accurately. Only with a suite of modeling tools will we learn to forecast – earthcast -- the future CZ.

Procedia Earth and Planetary Science, doi: 10.1016/j.proeps.2014.08.003

Aspect-dependent variations in regolith creep revealed by meteoric 10Be

Although variations in insolation and emergent feedbacks among soil moisture, vegetation, and soil cohesion are commonly invoked to explain topographic asymmetry that depends on aspect, few studies have directly quantified the efficiency of regolith transport along hillslopes of opposing aspect. We utilize meteoric 10Be concentrations in regolith (n = 74) to determine mass flux along equatorial-facing and polar-facing hillslopes in three forested upland watersheds in and adjacent to the Susquehanna Shale Hills Critical Zone Observatory in central Pennsylvania (USA). In combination with regolith depth measurements and high-resolution topography, these fluxes allow us to evaluate transport rate laws and the efficiency of regolith creep. Concentrations of meteoric 10Be in regolith along six separate transects imply that regolith flux is similar along all hillslopes, despite differences in topographic gradient and regolith thickness. Comparison of flux with regolith depth and topographic gradient reveals that transport depends on regolith depth, and that regolith creep is twice as efficient along low-gradient, south-facing slopes with thin regolith as compared to steep, north-facing slopes mantled with thicker regolith. We suggest that the observed topographic asymmetry in these watersheds has evolved over geologic time as a result of differences in the frequency of freeze-thaw events between hillslopes of opposing aspect.

Geology doi:DOI: 10.1130/G35357.1

Parameter estimation of a physically-based land surface hydrologic model using the ensemble Kalman Filter: A synthetic experiment

Flowchart of Flux-PIHM data assimilation framework for parameter and state estimation.

This paper presents multiple parameter estimation using multivariate observations via the ensemble Kalman filter (EnKF) for a physically based land surface hydrologic model. A data assimilation system is developed for a coupled physically based land surface hydrologic model (Flux-PIHM) by incorporating EnKF for model parameter and state estimation. Synthetic data experiments are performed at a first-order watershed, the Shale Hills watershed (0.08 km2). Six model parameters are estimated. Observations of discharge, water table depth, soil moisture, land surface temperature, sensible and latent heat fluxes, and transpiration are assimilated into the system. The results show that, given a limited number of site-specific observations, the EnKF can be used to estimate Flux-PIHM model parameters. All the estimated parameter values are very close to their true values, with the true values inside the estimated uncertainty range (1 standard deviation spread). The estimated parameter values are not affected by the initial guesses. It is found that discharge, soil moisture, and land surface temperature (or sensible and latent heat fluxes) are the most critical observations for the estimation of those six model parameters. The assimilation of multivariate observations applies strong constraints to parameter estimation, and provides unique parameter solutions. Model results reveal strong interaction between the van Genuchten parameters α and β, and between land surface and subsurface parameters. The EnKF data assimilation system provides a new approach for physically based hydrologic model calibration using multivariate observations. It can be used to provide guidance for observational system designs, and is promising for real-time probabilistic flood and drought forecasting.

Water Resources Research doi:10.1002/2013WR014070

The role of macropores and multi-resolution soil survey datasets for distributed surface–subsurface flow modeling

Five timestamps are selected: A (pre-storm condition), B (peakflow reached), C (at the beginning of recession), D (in the middle of recession), and E (after recession). 

Simulated spatial and response to a precipitation event for Cases I to III showing shallow groundwater table contours of pre-storm condition (column A), successive changes in saturated storage after precipitation event (B–A), (C–A), (D–A), (E–A). 

The average percentage changes in saturated storage relative to the total pore volume during runoff event.

Distributed watershed-scale modeling is often used as a framework for exploring the heterogeneity of runoff response and hydrologic performance of the catchment. The objective of this study is to apply this framework to characterizing the impacts of soil hydraulic properties at multiple scales on moisture storage and distributed runoff generation in a forested catchment. The physics-based and fully-coupled Penn State Integrated Hydrologic Model (PIHM) is employed to test a priori and field-measured properties in the modeling of watershed hydrology. PIHM includes an approximate representation of macropore flow that preserves the water holding capacity of the soil matrix while still allowing rapid flow through the macroporous soil under wet conditions. Both phenomena are critical to the overall hydrologic performance of the catchment. Soils data at different scales were identified: Case I STATSGO soils data (uniform or single soil type), Case II STATSGO soils data with macropore effect, and Case III field-based hydropedologic experiment revised distributed soil hydraulic properties and macropore property estimation. Our results showed that the Case I had difficulties in simulating the timing and peakflow of the runoff responses. Case II performed satisfactorily for peakflow at the outlet and internal weir locations. The distributed soils data in Case III demonstrated the model ability of predicting groundwater levels. The analysis suggests the important role of macropore flow to setting the threshold for recharge and runoff response, while still preserving the water holding capability of the soil and plant water availability. The spatial variability in soil hydraulic properties represented by Case III introduces an additional improvement in distributed catchment flow modeling, especially as it relates to subsurface lateral flow. Comparison of the three cases suggests the value of high-resolution soil survey mapping combined with a macropore parameterization can improve distributed watershed models.

Journal of Hydrology doi:10.1016/j.jhydrol.2014.02.055

2013

Spatiotemporal patterns of water stable isotope compositions at the Shale Hills Critical Zone: Linkages to subsurface hydrologic processes

Monthly amount-weighted isotope composition in precipitation, the kriged map of soil water isotope composition through depth and time, and monthly precipitation amount.  Red values correspond to less-depleted isotope composition while blue values correspond to more-depleted values.  Black circles correspond to catchment-averaged lysimeters (see text).  Ten-day moving averages of the groundwater level in the valley floor are plotted as a black line. Kriging uncertainty contours are plotted over the same space-time domain (see supplemental materials).  Ordinary kriging and kriging variance were computed in Matlab.

To better understand flow pathways and patterns in the subsurface, a stable isotope monitoring network was established at the Susquehanna-Shale Hills Critical Zone Observatory (SSHCZO). Soil water samples were collected approximately biweekly using suction-cup lysimeters installed at multiple depths along four different transects in the catchment. Groundwater and stream water were collected daily in the valley using automatic samplers, while precipitation samples were collected automatically on an event basis. The 3+ years (2008–2012) of monitoring data showed strong seasonal precipitation isotope compositions, which were imprinted in seasonal patterns of soil water at different spatial locations and depths. The groundwater isotope composition remained relatively constant throughout the year and closely matched the yearly amount-weighted precipitation average, suggesting groundwater received recharge water in each season, although recharge mechanisms differed between growing and nongrowing seasons. Soil water samples showed clear attenuation with depth, with the largest variability in the shallow soil water (≤30 cm) mirroring precipitation inputs, moderate variability in the intermediate depths (40–100 cm), and the least variability in the deep soil water (≥120 cm) where the average remained near the groundwater average. Soil water isotope composition profiles also provided clear evidence for preferential flow occurring both laterally and vertically in different seasons and at various soil depths in the catchment. Putting all together, the extensive dataset of soil water isotopic compositions obtained in this study have provided a number of insights into complex subsurface hydrologic processes that are transferable to other similar landscapes.

Vadose Zone Journal doi:10.2136/vzj2013.01.0029

Climate dependence of feldspar weathering in shale soils along a latitudinal gradient

The average integrated rate of Na loss (QNa) plotted versus MAT and MAP for transect sites. The 3D curve was calculated using fit parameters obtained from fitting Eq. (11) to data. Transect sites are projected onto the surface using MAT and MAP values listed in Table 1 of the article. 

Although regolith, the mantle of physically, chemically, and biologically altered material overlying bedrock, covers much of Earth’s continents, the rates and mechanisms of regolith formation are not well quantified. Without this knowledge, predictions of the availability of soil to sustain Earth’s growing population are problematic. To quantify the influence of climate on regolith formation, a transect of study sites has been established on the same lithology – Silurian shale – along a climatic gradient in the northern hemisphere as part of the Susquehanna Shale Hills Critical Zone Observatory, Pennsylvania, USA. The climate gradient is bounded by a cold/wet end member in Wales and a warm/wet end member in Puerto Rico; in between, mean annual temperature (MAT) and mean annual precipitation (MAP) increase to the south through New York, Pennsylvania, Virginia, Tennessee and Alabama. The site in Puerto Rico does not lie on the same shale formation as the Appalachian sites but is similar in composition. Soils and rocks were sampled at geomorphologically similar ridgetop sites to compare and model shale weathering along the transect. Focusing on the low-concentration, non-nutrient element Na, we observe that the extent and depth of Na depletion is greater where mean annual temperature (MAT) and precipitation (MAP) are higher. Na depletion, a proxy for feldspar weathering, is the deepest reaction documented in the augerable soil profiles. This may therefore be the reaction that initiates the transformation of high bulk-density bedrock to regolith of low bulk density. Based on the shale chemistry along the transect, the time-integrated Na release rate (QNa) increases exponentially as a function of MAT and linearly with MAP. NY, the only site with shale-till parent material, is characterized by a QNa that is 18 times faster than PA, an observation which is attributed to the increased surface area of minerals due to grinding of the glacier and kinetically limited weathering in the north. A calculated apparent Arrhenius-type temperature dependence across the transect (excluding NY) for the dissolution of feldspar (Na depletion) is 99 ± 15 kJ mol-1, a value similar to field-measured values of the activation energy (14–109 kJ mol-1) or laboratory-measured values of the enthalpy of the albite reaction (79.8 kJ mol-1). Observations from this transect document that weathering losses of Na from Silurian shale can be understood with models of regolith formation based on chemical and physical factors such that weathering progresses from kinetically limited sites (Wales to AL) to the transport-limited site in Puerto Rico. Significant advances in our ability to predict regolith formation will be made as we apply more quantitative models to such transect studies on shales and other rocks types.

Geochimica et Cosmochimica Acta 122 (2013) 101–126  doi:10.1016/j.gca.2013.08.001

Spatial Distribution of Tree Species Governs the Spatio-Temporal Interaction of Leaf Area Index and Soil Moisture across a Forested Landscape

Temporal dynamics of leaf area index (L: m2 m-2) across Susquehanna Shale Hills CZO during 2010.

Temporal dynamics of surface (10 cm) volumetric soil water content (θ: m3 m-3) across Susquehanna Shale Hills CZO during 2010.

Quantifying coupled spatio-temporal dynamics of phenology and hydrology and understanding underlying processes is a fundamental challenge in ecohydrology. While variation in phenology and factors influencing it have attracted the attention of ecologists for a long time, the influence of biodiversity on coupled dynamics of phenology and hydrology across a landscape is largely untested. We measured leaf area index (L) and volumetric soil water content (h) on a co-located spatial grid to characterize forest phenology and hydrology across a forested catchment in central Pennsylvania during 2010. We used hierarchical Bayesian modeling to quantify spatio-temporal patterns of L and h. Our results suggest that the spatial distribution of tree species across the landscape created unique spatio-temporal patterns of L, which created patterns of  water demand reflected in variable soil moisture across space and time. We found a lag of about 11 days between increase in L and decline in h. Vegetation and soil moisture become increasingly homogenized and coupled from leaf-onset to maturity but heterogeneous and uncoupled from leaf maturity to senescence. Our results provide insight into spatiotemporal coupling between biodiversity and soil hydrology that is useful to enhance ecohydrological modeling in humid temperate forests.

PLoS ONE 8(3): e58704. doi:10.1371/journal.pone.0058704

Regolith production and transport in the Susquehanna Shale Hills Critical Zone Observatory, Part 1: Insights from U-series isotopes

Measured (234U/238U) and (230Th/238U) activity ratios for regolith samples from NPRT, NPMS, and NPVF plotted against regolith depth.

To investigate the timescales of regolith formation on hillslopes with contrasting topographic aspect, we measured U-series isotopes in regolith profiles from two hillslopes (north- vs. south-facing) within the east-west trending Shale Hills catchment in Pennsylvania. This catchment is developed entirely on the Fe-rich, Silurian Rose Hill gray shale. Hillslopes exhibit a topographic asymmetry: the north-facing hillslope has an average slope gradient of ~20°, slightly steeper than the south-facing hillslope (~15°). The regolith samples display significant U-series disequilibrium resulting from shale weathering. Based on the U-series data, the rates of regolith production on the two ridge tops are indistinguishable (40 ± 22 versus 45 ± 12 m/Ma). However, when downslope positions are compared, the regolith profiles on the south-facing hillslope are characterized by faster regolith production rates (50±15 to 52±15 m/Ma) and shorter durations of chemical weathering (12±3 to 16±5 ka) than those on the north-facing hillslope (17±14 to 18±13 m/Ma and 39±20 to 43±20 ka). The south-facing hillslope is also characterized by faster chemical weathering rates inferred from major element chemistry, despite lower extents of chemical depletion. These results are consistent with the influence of aspect on regolith formation at Shale Hills; we hypothesize that aspect affects such variables as temperature, moisture content and evapotranspiration in the regolith zone, causing faster chemical weathering and regolith formation rates on the south-facing side of the catchment. The difference in microclimate between these two hillslopes is inferred to have been especially significant during the peri-glacial period that occurred at Shale Hills at least ~15 ka before present. At that time, the erosion rates may also have been different from those observed today, perhaps denuding the south-facing hillslope of regolith but not quite stripping the north-facing hillslope. An analysis of hillslope evolution and response timescales with a linear mass transport model shows that the current landscape at Shale Hills is not in geomorphologic steady state (i.e., so-called dynamic equilibrium), but rather is likely still responding to the climate shift from the Holocene periglacial to the modern, temperate conditions.

Journal of Geophysical Research - Earth Surface doi: 10.1002/jgrf.20037

Above- and belowground controls on water use by trees of different wood types in an eastern US deciduous forest

Map of the study site showing locations of trees in which sap flow was measured, locations from which soil water content and water potential data were obtained (SMS), and the location of the instrument tower from which relative humidity, air temperature and solar radiation data were obtained (Met Tower). Species abbreviations: ACSA, Acer saccharum; LITU, Liriodendron tulipifera; PIVI, Pinus virginiana; QUAL, Quercus alba; QUPR, Quercus prinus; QURU, Quercus rubra; TSCA, Tsuga canadensis.

Stomata control tree transpiration by sensing and integrating environmental signals originating in the atmosphere and soil, and co-occurring species may differ in inherent stomatal sensitivity to these above- and belowground signals and in the types of signals to which they respond. Stomatal responsiveness to environmental signals is likely to differ across species having different types of wood (e.g., ring-porous, diffuse-porous and coniferous) because each wood type differs in the structure, size and spatial distribution of its xylem conduits as well as in the scaling of hydraulic properties with stem diameter.  The objective of this study was to evaluate the impact of variation in soil water availability and atmospheric evaporative demand on stomatal regulation of transpiration in seven co-occurring temperate deciduous forest species representing three wood types. We measured whole-tree sap flux and soil and atmospheric variables in a mixed deciduous forest in central Pennsylvania over the course of a growing season characterized by severe drought and large fluctuations in atmospheric vapor pressure deficit (D). The relative sensitivity of sap flux to soil drying was ~2.2–2.3 times greater in the diffuse-porous and coniferous species than in the ring-porous species. Stomata of the ring-porous oaks were only about half as responsive to increased D as those of trees of the other two wood types. These differences in responsiveness to changes in the below and aboveground environment implied that regulation of leaf water potential in the ring-porous oaks was less stringent than that in the diffuse-porous angiosperms or the conifers. The results suggest that increases in the frequency or intensity of summer droughts in the study region could have multiple consequences for forest function, including altered successional time courses or climax species composition and cumulative effects on whole-tree architecture, resulting in a structural and physiological legacy that would restrict the ability of trees to respond rapidly to more favorable growth conditions.

Tree Physiology 00, 1–12 doi:10.1093/treephys/tpt012

2012

Fe cycling in the Shale Hills Critical Zone Observatory, Pennsylvania: An analysis of biogeochemical weathering and Fe isotope fractionation

Contour maps showing the distribution of concentrations of HCl-extractable Fe (a) and δ56Fe values in bulk Fe (b) across the south planar transect. Small circles indicate locations from which samples were obtained. Soil depth in cm on y-axis, while distance along ridge on x-axis.

During weathering, Fe in primary minerals is solubilized by ligands and/or reduced by bacteria and released into soil porewaters. Such Fe is then removed or reprecipitated in soils. To understand these processes, we analyzed Fe chemistry and isotopic composition in regolith of the Shale Hills watershed, a Critical Zone Observatory in central Pennsylvania overlying iron-rich shale of the Rose Hill Formation. Elemental concentrations were measured in soil from a well-drained catena on a planar hillslope on the south side of the catchment. Based upon X-ray diffraction and bulk elemental data, loss of Fe commences as clay begins to weather ∼15 cm below the depth of auger-refusal. More Fe(III) was present than Fe(II) in all soil samples from the ridge top to the valley floor. Both total and ferrous iron are depleted from the land surface of catena soils relative to the bedrock. Loss of ferrous Fe is attributed mostly to abiotic or biotic oxidation. Loss of Fe is most likely due to transport of micron-sized particles that are not sampled by porous-cup lysimeters, but which are sampled in stream and ground waters. The isotopic compositions (δ56Fe, relative to IRMM-014) of bulk Fe and 0.5 N HCl-extracted Fe (operationally designed to remove amorphous Fe (oxyhydr)oxides) range between −0.3‰ and +0.3‰, with Δ56Febulk-extractable values between ∼0.2‰ and 0.4‰. Throughout the soils along the catena, δ56Fe signatures of both bulk Fe and HCl-extracted Fe become isotopically lighter as the extent of weathering proceeds. The isotopic trends are attributed to one of two proposed mechanisms. One mechanism involves Fe fractionation during mobilization of Fe from the parent material due to either Fe reduction or ligand-promoted dissolution. The other mechanism involves fractionation during immobilization of Fe (oxyhydr)oxides. If the latter mechanism is true, then shale – which comprises one quarter of continental rocks – could be an important source of isotopically heavy Fe for rivers.

Geochimica et Cosmochimica Acta, Volume 99, 15 December 2012, Pages 18–38, doi:10.1016/j.gca.2012.09.029

Changing controls of soil moisture spatial organization in the Shale Hills Catchment

Plot of the temporal evolution of averaged soil moisture content within the three soil–landform units over 91 measurement days from 2004 to 2008 at each measurement depth. Shaded areas indicate dry summer periods (06/16-10/02/2005, 06/30-08/31/2006, 06/06-10/06/2007, and 05/27-10/19/2008). Topmost panel is depth to water table (W.T.) at site #15 located in the valley (B.G.S. = below ground surface). Second row panels are potential evapotranspiration calculated using the climate station data. Rainfall is represented as the cumulative precipitation in the preceding 3 days (3 day antecedent precipitation index, API3).

Modeling hydrological processes often requires the identification of dominant controls on soil moisture spatial organization under different climatic conditions at various soil depths. In this study, we utilized a four-year database consisting of soil moisture measurements at 106 locations from near-surface down to 1.1 m depth across a forested catchment in central Pennsylvania, USA. Our objectives were to 1) compare the spatial organization of soil moisture within different soil–landform units and its temporal persistence at different depths under varying catchment wetness conditions and 2) investigate correlation strength between soil moisture content and 11 soil–terrain attributes and the temporal change of such correlation. Our results showed that the catchment's near-surface (< 0.3 m) soil moisture organization exhibited clear seasonal trends: during winter through early summer, areas of high soil moisture were concentrated within convergent landforms; while during summer through early fall, soil moisture was more uniformly distributed throughout this complex terrain catchment. This suggests that under dry conditions soil moisture removal (primarily through evapotranspiration) had a significant influence on the organization of near-surface soil moisture, while topography was an important control on soil moisture spatial organization under wet conditions. Subsurface (> 0.3 to 1.1 m) soil moisture organization, however, exhibited increasing temporal persistence with depth, and subsoil moisture above the catchment-wide average was concentrated within convergent landforms under both wet and dry conditions. Topographic wetness index, slope, depth to bedrock, and percent (by weight) clay and rock fragment were significant (p < 0.05) factors influencing soil moisture content on at least 80% of 91 measurement days analyzed for all soil depths. Intermediate depths (> 0.3 to < 0.7 m) exhibited the highest coefficient of determination (R2) in linear regression for topographic wetness index, suggesting that lateral subsurface flow may be an important driver of soil moisture dynamics at these depths in this catchment. Mean R2 values for slope, depth to bedrock, and percent clay and rock fragment increased with increasing depth, confirming the importance of deep soil moisture storage on subsoil moisture organization. We conclude that the controls on this catchment's soil moisture spatial organization at the near-surface (< 0.3 m) fluctuates seasonally between evapotranspiration and topography; that at intermediate depths (0.3 to 0.7 m) the soil moisture organization is controlled significantly by lateral subsurface flow; and that the organization at deeper depths (> 0.7 m) becomes more temporally persistent and is primarily a function of both topography and soil depth.

Geoderma 173-174:289-302, doi:10.1016/j.geoderma.2011.11.003

Insights into the weathering of black shale: Cu isotopes and concentrations in the Marcellus Formation shale, Huntingdon County, Pennsylvania (U.S.A.).

Cu isotope compositions of samples of soil (from ridgetop sites RT1, RT2, and valley floor site VFS), parent rocks (grey bar, calculated as the average Cu isotope measurements from parent materials n=7), and pore waters (lysimeter samples).

The complex weathering processes which govern the production of soil from bedrock have proven difficult to understand for many lithologies. Weathering of black shale is of particular interest because it releases organic carbon and heavy metals as solutes and therefore impacts the health of terrestrial and aquatic ecosystems. To understand black shale weathering, a geochemical survey was initiated for soils developed on shales of the Marcellus Formation at a zero-order catchment at a satellite site of the Susquehanna/Shale Hills Critical Zone Observatory located in Jackson Corner, Pennsylvania. This formation is an organic- and metal-rich, carbonaceous shale that underlies much of New York, Pennsylvania, Ohio and West Virginia. In this paper, we focus on the effects of weathering on variations of Cu isotopes in the shale. Cu concentration data for soil were normalized using Ti concentrations to document the mobility of Cu relative to bedrock. At both the ridgetop and valley floor, depletion profiles for Cu are documented in the soils. The Cu in the soils is depleted in 65Cu (average δ65Cu = − 0.5‰ ± 0.2) compared to the parent material (average δ65Cu = 0.03‰ ± 0.15). Consistent with loss of Cu from soils, the pore waters contain 10 ppb Cu on average and are enriched in the heavy isotope (average value δ65Cu = 1.14‰ ± 0.44). Rayleigh fractionation models using the concentration and isotope data of the soils are consistent with pyrite weathering and loss of Cu from the ridgetop, but downslope transport and Cu re-precipitation at the valley floor.

Chemical Geology 304–305:175–184, doi: 10.1016/j.chemgeo.2012.02.015.

2011

Hydraulic patterns and safety margins, from stem to stomata, in three eastern US tree species

Measured water potential of leaves predawn, midday, and stem water potential measured midday.  ACRU is A. rubrum, LITU is L. tulipifera, and PIVI is P. virginiana.  Vertical error bars represent standard error. 

Adequate water transport is necessary to prevent stomatal closure and allow for photosynthesis. Dysfunction in the water transport pathway can result in stomatal closure, and can be deleterious to overall plant health and survival. Although much is known about small branch hydraulics, little is known about the coordination of leaf and stem hydraulic function. Additionally, the daily variations in leaf hydraulic conductance (Kleaf), stomatal conductance and water potential (ΨL) have only been measured for a few species. The objective of the current study was to characterize stem and leaf vulnerability to hydraulic dysfunction for three eastern US tree species (Acer rubrum, Liriodendron tulipifera and Pinus virginiana) and to measure in situ daily patterns of Kleaf, leaf and stem Ψ, and stomatal conductance in the field. Sap flow measurements were made on two of the three species to compare patterns of whole-plant water use with changes in Kleaf and stomatal conductance. Overall, stems were more resistant to hydraulic dysfunction than leaves. Stem P50 (Ψ resulting in 50% loss in conductivity) ranged from −3.0 to −4.2 MPa, whereas leaf P50 ranged from −0.8 to −1.7 MPa. Field ΨL declined over the course of the day, but only P. virginiana experienced reductions in Kleaf (nearly 100% loss). Stomatal conductance was greatest overall in P. virginiana, but peaked midmorning and then declined in all three species. Midday stem Ψ in all three species remained well above the threshold for embolism formation. The daily course of sap flux in P. virginiana was bell-shaped, whereas in A. rubrum sap flux peaked early in the morning and then declined over the remainder of the day. An analysis of our data and data for 39 other species suggest that there may be at least three distinct trajectories of relationships between maximum Kleaf and the % Kleaf at Ψmin. In one group of species, a trade-off between maximum Kleaf and % Kleaf at Ψmin appeared to exist, but no trade-off was evident in the other two trajectories.

Tree Physiol (2011) 31 (6): 659-668, doi: 10.1093/treephys/tpr050

Quantifying Solute Transport at the Shale Hills Critical Zone Observatory

Measured breakthrough curves with advection-dispersion equation (ADE), mobile-immobile model (MIM), and continuous time random walk (CTRW) solutions for soil cores from (A) 0- to 0.2-, (B) 0.6- to 0.8-, (C) 1.6- to 1.8-, (D) 2.3- to 2.5-m depths.  EC = effluent electrical conductivity, Br = Br- concentration.

We collected and analyzed Br breakthrough curve (BTC) data to identify the parameters controlling transport from a series of soil cores and a field-scale tracer test at the Shale Hills Critical Zone Observatory (SH-CZO) in central Pennsylvania. The soil cores were retrieved from a continuous hole that extended through the soil profile to quantify also how solute transport behavior changes with depth and weathering. Additionally, we performed a field-scale doublet tracer test to determine transport behavior in the weathered shale bedrock. Hydraulic conductivity and porosity were as low as 10−15 m s−1 and 0.035, respectively, in the shale bedrock and upward of 10−5 m s−1 and 0.45, respectively, in the shallow soils. Bromide BTCs demonstrated significant tailing in soil cores and field tracer experiments, which does not fit classical advection–dispersion processes. To quantify the behavior, numerical simulation of solute transport was performed with both a mobile–immobile (MIM) model and a continuous-time random walk (CTRW) approach. One-dimensional MIM modeling results yielded low mass transfer rates (<1 d−1) coupled with large immobile domains (immobile/mobile porosity ratio of 1.5–2). The MIM modeling results also suggested that immobile porosity was a combination of soil texture, fractures, and porosity development on shale fragments. One-dimensional CTRW results yielded a parameter set indicative of a transport regime that is consistently non-Fickian within the soil profile and bedrock. These modeling results confirm the important role of preferential flow paths, fractures, and mass transfer between more- and less-mobile fluid domains and advance the need to incorporate a continuum of mass transfer rates to more accurately quantify transport behavior within the SH-CZO.

Vadose Zone Journal (2011)  Vol. 10 No. 3, p. 843-857, doi: 10.2136/vzj2010.0130

Hot Spots and Hot Moments of Dissolved Organic Carbon Export and Soil Organic Carbon Storage in the Shale Hills Catchment

Interpolated maps of soil organic carbon (SOC) content (% w/w) in the (a) A and (b) B soil horizons, respectively, and SOC storage (g cm-2) within ≤ 1.1-m solum based on 56 sampling sites shown on each map.

Dissolved organic carbon (DOC) export from watersheds and soil organic carbon (SOC) storage are intimately linked in the terrestrial carbon cycle. However, predictions of hot spots and hot moments of DOC and SOC in watersheds remain uncertain because of high spatiotemporal variability and changing controls. In this study, we investigated the linkage between SOC storage and landform units across the 7.9-ha Shale Hills Critical Zone Observatory (CZO) and its implications for potential hot spots of DOC.  We also examined the trends of DOC in soil pore water along two hillslopes of contrasting soils and topography and the impacts of rainfall, stream discharge, and stream temperature on DOC export to identify possible hot moments. Based on the SOC distribution throughout the entire catchment, swales (particularly south-facing swales) were identified as hot spots because they exhibited significantly higher SOC storage and more active hydrology as compared tothe rest of the catchment. Along the two hillslopes reported here, average soil pore water DOC concentrations were noticeably higher (35 ± 12%) along the swale as compared to the planar hillslope. Soil pore water DOC concentrations were elevated at the soil–bedrock interface at the ridgetop and at the Bw–Bt horizon interface in the valley floor, suggesting transport-driven hot spots along restrictive layer interfaces. Stream water DOC concentration at the catchment outlet averaged 6.2 ± 5.3 mg L−1 from May 2008 to October 2010, which was significantly correlated with stream discharge and stream water temperature. Transport-driven hot moments of stream water DOC were observed during the periods of snowmelt and late summer to early fall wet-up, which together contributed ~55% of total stream water DOC exported in 2009. This reflected the control of DOC export by flushing (linked to discharge) and biological activity (related to temperature) and its variation during different seasons of a year. This study showcased the impacts of complex soil and topography interactions—coupled with changing weather and seasonal biological activity—on the spatiotemporal dynamics of DOC export in a temperate forested catchment and its link to SOC distribution.

Vadose Zone Journal (2011) Vol. 10 No. 3, p. 943-954, doi: 10.2136/vzj2010.0149

Soils Reveal Widespread Manganese Enrichment from Industrial Inputs

Normalized Mn concentrations (τZr,Mn, eq 1) are plotted as different symbols versus depth for 21 ridgetop soil cores at the Susquehanna Shale Hills Observatory (SSHO).  Here, depth in the soil is normalized so that 1 ) depth of the bedrock-soil interface. Soils document Mn enrichment (τZr,Mn > 0) near the surface that approaches underlying parent rock composition (τZr,Mn ) 0) at depth. Surface soils are up to 13 times more enriched in Mn than parent. Error bars represent the propagated uncertainty in chemical analyses ((3%). Where no errors are shown, the bars are the size of symbols or smaller.

It is well-known that metals are emitted to the air by human activities and subsequently deposited to the land surface; however, we have not adequately evaluated the geographic extent and ecosystem impacts of industrial metal loading to soils. Here, we demonstrate that atmospheric inputs have widely contaminated soils with Mn in industrialized regions. Soils record elemental fluxes impacting the Earth’s surface and can be analyzed to quantify inputs and outputs during pedogenesis. We use a mass balance model to interpret details of Mn enrichment by examining soil, bedrock, precipitation, and porefluid chemistry in a first-order watershed in central Pennsylvania, USA. This reveals that 53% of Mn in ridge soils can be attributed to atmospheric deposition from anthropogenic sources. An analysis of published data sets indicates that over half of the soils surveyed in Pennsylvania (70%), North America (60%), and Europe (51%) are similarly enriched in Mn. We conclude that soil Mn enrichment due to industrial inputs is extensive, yet patchy in distribution due to source location, heterogeneity of lithology, vegetation, and other attributes of the land surface. These results indicate that atmospheric transport must be considered a potentially critical component of the global Mn cycle during the Anthropocene.

Environmental Science & Technology 45 (1):241-247, doi:10.1021/es102001w

2010

An object-oriented shared data model for GIS and distributed hydrologic models

Distributed physical models for the space–time distribution of water, energy, vegetation, and mass flow require new strategies for data representation, model domain decomposition, a priori parameterization, and visualization. The geographic information system (GIS) has been traditionally used to accomplish these data management functionalities in hydrologic applications. However, the interaction between the data management tools and the physical model are often loosely integrated and nondynamic. This is because (a) the data types, semantics, resolutions, and formats for the physical model system and the distributed data or parameters may be different, with significant data preprocessing required before they can be shared; (b) the management tools may not be accessible or shared by the GIS and physical model; and (c) the individual systems may be operating system dependent or are driven by proprietary data structures. The impediment to seamless data flow between the two software components has the effect of increasing the model setup time and analysis time of model output results, and also makes it restrictive to perform sophisticated numerical modeling procedures (real-time forecasting, sensitivity analysis, etc.) that utilize extensive GIS data. These limitations can be offset to a large degree by developing an integrated software component that shares data between the (hydrologic) model and the GIS modules. We contend that the prerequisite for the development of such an integrated software component is a ‘shared data model’, which is designed using an object-oriented strategy. Here we present the design of such a shared data model taking into consideration the data type descriptions, identification of data classes, relationships, and constraints. The developed data model has been used as a method base for developing a coupled GIS interface to Penn State Integrated Hydrologic Model (PIHM), called PIHMgis.

International Journal of Geographical Information Science 24 (7):1061-1079, doi:10.1080/13658810903289460.

Mineral weathering and elemental transport during hillslope evolution at the Susquehanna/Shale Hills Critical Zone Observatory

Located in the uplands of the Valley and Ridge physiographic province of Pennsylvania, the Susquehanna/Shale Hills Critical Zone Observatory (SSHO) is a tectonically quiescent, first-order catchment developed on shales of the Silurian Rose Hill Formation.  We used soil cores augered at the highest point of the watershed and along a subsurface water flowline on a planar hillslope to investigate mineral transformations and physical/chemical weathering fluxes. About 25 m of bedrock was also drilled to estimate parent composition. Depletion of carbonate at tens of meters of depth in bedrock may delineate a deep carbonate-weathering front. Overlying this, extending from ˜6 m below the bedrock–soil interface up into the soil, is the feldspar dissolution front. In the soils, depletion profiles for K, Mg, Si, Fe, and Al relative to the bedrock define the illite and chlorite reaction fronts. When combined with a cosmogenic nuclide-derived erosion rate on watershed sediments, these depletion profiles are consistent with dissolution rates that are several orders of magnitudes slower for chlorite (1–5 x 10-17 mol m-2 s-1) and illite (2–9 x 10-17 mol m-2 s-1) than observed in the laboratory. Mineral reactions result in formation of vermiculite, hydroxy-interlayered vermiculite, and minor kaolinite. During weathering, exchangeable divalent cations are replaced by Al as soil pH decreases.  The losses of Mg and K in the soils occur largely as solute fluxes; in contrast, losses of Al and Fe are mostly as downslope transport of fine particles. Physical erosion of bulk soils also occurs: results from a steady-state model demonstrate that physical erosion accounts for about half of the total denudation at the ridgetop and midslope positions. Chemical weathering losses of Mg, Na, and K are higher in the upslope positions likely because of the higher degree of chemical undersaturation in porewaters. Chemical weathering slows down in the valley floor and Al and Si even show net accumulation. The simplest model for the hillslope that is consistent with all observations is a steady-state, clay weathering-limited system where soil production rates decrease with increasing soil thickness.

Geochimica et Cosmochimica Acta 74 (13):3669-3691, doi:10.1016/j.gca.2010.03.036

Dynamical modelling of concentration–age–discharge in watersheds

There is now a wide literature on the use of tracer age and transit time distributions to diagnose transport in environmental systems. Theories have been proposed using idealized tracer age modelling for ocean ventilation, atmospheric circulation, soil, stream and groundwater flow. Most approaches assume a steady flow regime and stationarity in the concentration (tracer) distribution function for age, although recent work shows that this is not a necessary assumption. In this paper, dynamic model for flow, concentration, and age in volume-averaged and a spatially distributed watershed system are derived in terms of the moments of the underlying distribution function for tracer age, time, and position. Several theoretical and practical issues are presented: (1) The low-order moments of the age distribution function are sufficient to construct a dynamical system for the mean age and concentration under steady or transient flow conditions. (2) Solutions to the coupled system of equations for flow, concentration and age show that ‘age’ of solutes stored within the watershed or leaving the watershed is a dynamic process which depends on flow variations as well as the solute or tracer dynamics. (3) Intermittency of wetting and drying cycles leads to an apparent increase in the tracer age in proportional to the duration of the ‘dry’ phase. (4) The question of how mobile/immobile flow may affect the age of solutes is examined by including a low permeable, passive store that relaxes the well-mixed assumption. (5). A spatially distributed advective and dispersive transport solution for age evolution over a simple 1-D hillslope is developed to demonstrate the age theory for a distributed source of water and tracer, and the solution is shown to have very similar input–output behaviour when compared to the volume-average model for comparable parameters.

Hydrological Processes Journal 24:1711-1718, DOI: 10.1002/hyp.7691

Cu isotope compositions of samples of soil (from ridgetop sites RT1, RT2, and valley floor site VFS), parent rocks (grey bar, calculated as the average Cu isotope measurements from parent materials n=7), and pore waters (lysimeter samples).

Plot of the temporal evolution of averaged soil moisture content within the three soil–landform units over 91 measurement days from 2004 to 2008 at each measurement depth. Shaded areas indicate dry summer periods (06/16-10/02/2005, 06/30-08/31/2006, 06/06-10/06/2007, and 05/27-10/19/2008). Topmost panel is depth to water table (W.T.) at site #15 located in the valley (B.G.S. = below ground surface). Second row panels are potential evapotranspiration calculated using the climate station data. Rainfall is represented as the cumulative precipitation in the preceding 3 days (3 day antecedent precipitation index, API3).

Measured water potential of leaves predawn, midday, and stem water potential measured midday.  ACRU is A. rubrum, LITU is L. tulipifera, and PIVI is P. virginiana.  Vertical error bars represent standard error. 

Measured breakthrough curves with advection-dispersion equation (ADE), mobile-immobile model (MIM), and continuous time random walk (CTRW) solutions for soil cores from (A) 0- to 0.2-, (B) 0.6- to 0.8-, (C) 1.6- to 1.8-, (D) 2.3- to 2.5-m depths.  EC = effluent electrical conductivity, Br = Br- concentration.

Map of the study site showing locations of trees in which sap flow was measured, locations from which soil water content and water potential data were obtained (SMS), and the location of the instrument tower from which relative humidity, air temperature and solar radiation data were obtained (Met Tower). Species abbreviations: ACSA, Acer saccharum; LITU, Liriodendron tulipifera; PIVI, Pinus virginiana; QUAL, Quercus alba; QUPR, Quercus prinus; QURU, Quercus rubra; TSCA, Tsuga canadensis.

Measured (234U/238U) and (230Th/238U) activity ratios for regolith samples from NPRT, NPMS, and NPVF plotted against regolith depth.

Temporal dynamics of leaf area index (L: m2 m-2) across Susquehanna Shale Hills CZO during 2010.

Interpolated maps of soil organic carbon (SOC) content (% w/w) in the (a) A and (b) B soil horizons, respectively, and SOC storage (g cm-2) within ≤ 1.1-m solum based on 56 sampling sites shown on each map.

Temporal dynamics of surface (10 cm) volumetric soil water content (θ: m3 m-3) across Susquehanna Shale Hills CZO during 2010.

Contour maps showing the distribution of concentrations of HCl-extractable Fe (a) and δ56Fe values in bulk Fe (b) across the south planar transect. Small circles indicate locations from which samples were obtained. Soil depth in cm on y-axis, while distance along ridge on x-axis.

Normalized Mn concentrations (τZr,Mn, eq 1) are plotted as different symbols versus depth for 21 ridgetop soil cores at the Susquehanna Shale Hills Observatory (SSHO).  Here, depth in the soil is normalized so that 1 ) depth of the bedrock-soil interface. Soils document Mn enrichment (τZr,Mn > 0) near the surface that approaches underlying parent rock composition (τZr,Mn ) 0) at depth. Surface soils are up to 13 times more enriched in Mn than parent. Error bars represent the propagated uncertainty in chemical analyses ((3%). Where no errors are shown, the bars are the size of symbols or smaller.

The average integrated rate of Na loss (QNa) plotted versus MAT and MAP for transect sites. The 3D curve was calculated using fit parameters obtained from fitting Eq. (11) to data. Transect sites are projected onto the surface using MAT and MAP values listed in Table 1 of the article. 

Monthly amount-weighted isotope composition in precipitation, the kriged map of soil water isotope composition through depth and time, and monthly precipitation amount.  Red values correspond to less-depleted isotope composition while blue values correspond to more-depleted values.  Black circles correspond to catchment-averaged lysimeters (see text).  Ten-day moving averages of the groundwater level in the valley floor are plotted as a black line. Kriging uncertainty contours are plotted over the same space-time domain (see supplemental materials).  Ordinary kriging and kriging variance were computed in Matlab.

Flowchart of Flux-PIHM data assimilation framework for parameter and state estimation.

Five timestamps are selected: A (pre-storm condition), B (peakflow reached), C (at the beginning of recession), D (in the middle of recession), and E (after recession). 

Simulated spatial and response to a precipitation event for Cases I to III showing shallow groundwater table contours of pre-storm condition (column A), successive changes in saturated storage after precipitation event (B–A), (C–A), (D–A), (E–A). 

The average percentage changes in saturated storage relative to the total pore volume during runoff event.

Conceptual diagram of metal cycling in soils at SSHO

PIHM structure and processes. The upper subfigure shows the hydrological processes of PIHM at a cross section of a watershed. The lower subfigure shows the spatial structure of PIHM. Blue lines represent the stream channels, and triangles represent the catchment domain.

Measured (M) and interpolated (I) preferential flow (PF) occurrence frequency based on the data from the period May 2011 to June 2012 for all 35 sites, with all the available data in each soil profile used in estimation (with the bar height representing the size of PF frequency, which ranged from <1% at Site 31 to 70.1% at Site 67). The interpolation map was generated using regression kriging.

Depiction of Shale Hills catchment and related plant and soil traits with relative slope position and relative elevation.

Map of Shavers Creek watershed, highlighting (a) topography derived from airborne lidar, (b) geology (Berg et al., 1980), and (c) land use (Homer et al., 2015). In moving from “measure everything everywhere” (our paradigm in the 8 ha Shale Hills catchment; SH) to “measure only what is needed” in the Shavers Creek watershed (164 km2), we chose to investigate two new first-order subcatchments: a
forested sandstone site (along Garner Run, marked GR) and an agricultural calcareous shale site (to be determined). In addition, three sites on Shavers Creek have been chosen as stream discharge and chemistry monitoring sites (marked SCAL – Shavers Creek above lake; SCBL – Shavers Creek below lake; and SCO – Shavers Creek outlet). 

Definition sketch of 3-D control volume on hillslope, where z(x, y, t) ¼ ground surface elevation (m), e(x, y, t) ¼ bedrock-regolith interface elevation (m), h(x, y, t) ¼ regolith thickness in vertical (m) (slope-normal thickness, H, is equal to h cos q, where q is the landscape surface slope), U(x, y, t) ¼ rock uplift rate (m yr1), Bw(x, y, t) ¼ bedrock weathering rate (m yr1) in the vertical direction, qs ¼ the surface flux of regolith sediment by overland flow (m2 yr1), and qc ¼ the lateral volumetric regolith flux by creep and tree throw (m2 yr1) entering through the sides of the control volume. The bottom of the stream may consist of alluvial deposits on bedrock.

Bulk soil mineralogy as a function of depth across the climosequence. Vermic/HIV includes vermiculite and hydroxy-interlayered vermiculite minerals.

Vessel size distribution for one branch sample from each tree species with counts of vessels in each size category. Between five and nine cross sections were analysed per branch. Total number of vessels analysed:  QUPR, n = 463; QURU, n = 315; CATO, n = 268, ACSA, n = 573.

Maps of terrain attributes derived from the digital elevation model created from 1-m resolution LiDAR: (a) local slope [m m− 1], (b) depth to bedrock (DtB) [cm], (c) vertical distance to stream (VDS) [m], (d) natural log-transformed upslope contributing area (UCA) [ln (m2)], (e) topographic wetness index (TWI) [−], and (f) surface curvature [−].

A schematic representation of processes in RT-Flux-PIHM. Flux-PIHM simulates the hydrological and land-surface dynamics (precipitation,
canopy interception, infiltration, recharge, overland flow, subsurface lateral flow, river flow, and surface energy balance) at the watershed scale using the semidiscrete finite volume method. 

A schematic representation of processes in different modules (different colors) of RT-Flux-PIHM. The model allows systematic
understanding of coupled processes at the grid and watershed scales. Note that three types of flow contribute to the stream discharge
QD: overland flow (surface runoff, Qs), subsurface lateral flow in soil (QL), and deep groundwater flow (QG).