Cliff Retreat Rates Associated with a Low-Level Radioactive Waste Disposal Facility in Los Alamos, New Mexico, USA

Abstract

We present an analysis and interpretation of potential cliff stability at a low-level waste disposal facility at Los Alamos National Laboratory, New Mexico, using cliff morphologic and fracture characteristics coupled with carbon-14 surface exposure dating. Our study is important as it directly bears on the licensing criteria for low-level radioactive waste sites. We find that future characteristic cliff failures will likely not breach disposal pits and shafts over the 1000-year minimum regulatory period. Further, we find, using a multivariate regression model, that slope angle and cliff face aspect are sub-equal in importance to predict regions of high risk of failure when combined with surface exposure ages and assuming that old exposure ages are most indicative of stability (instability) and therefore can aid decision making in final design implementation.

1. Introduction

The need for safe and stable waste disposal facilities in well-characterized geologic repositories will increase as the demand for alternative fuel and energy sources increases and as the U.S. Department of Energy (DOE) continues its environmental cleanup mission. Waste disposal facilities must adhere to specific national regulations regarding geological conditions at a site, including distance from cliff edges. In the United States, such regulations are provided in the Resource Conservation and Recovery Act [1], overseen by the Environmental Protection Agency. Separate mandates for DOE facilities are implemented through DOE Order 435.1, Radioactive Waste Management [2]. Safe operation and future closure of the Material Disposal Area G (MDA G) low-level radioactive waste (LLRW) site on the Pajarito Plateau in northern New Mexico (Figure 1) are regulated by DOE Order 435.1, which requires that potential exposure for members of the public over the next 1000 years be less than 10 millirem yr⁻¹ at the fence line [2]. This effectively requires that disposed waste be isolated from the facility surface for at least 1000 years of storage time before any possible exposure. Characterizing the types and rates of current and future geologic changes at this site is critical to ensure adherence to these requirements.

Cliff retreat rates in coastal regions are well characterized and studied in part due to higher population density [3,4,5,6,7], but cliff retreat rates in arid environments are less thoroughly studied [8,9]. Existing and proposed waste disposal facility locations in the southwestern and western U.S. (e.g., Yucca Mountain) are considered viable due to arid conditions and low population density; however, these sites require detailed characterization to fully assess long-term viability for long-term waste storage.

For example, the Pajarito Plateau of northern New Mexico is characterized by long, narrow “finger” mesas and incised canyons in Quaternary volcanic deposits. The mesa on which MDA G is located is referred to as Mesita del Buey (MdB; Figure 1B). MdB is a relatively narrow (120 to 400 m), gently slo** finger mesa that is lined by sub-vertical cliffs on the southside and less steep slopes on the northside (Figure 1C). Fracturing around most of the perimeter of MdB, particularly the south-facing cliffs, introduces the potential for rock falls. Large, detached blocks on the canyon floor and lower cliff slopes provide evidence of cliff failure and estimates of characteristic block size and thus cliff retreat rate.

We aim to determine the impact of future cliff retreat within a 1000-year timeframe on the LLRW facility’s ability to operate within performance objectives. Our case study of the MDA G LLRW facility in northern New Mexico incorporates field observations, image analysis, targeted geologic sampling, cosmogenic nuclide exposure dating of exposed cliff faces, and geospatial analyses to quantitatively assess slope stability in an arid environment. Our slope stability assessments characterize current conditions and failure mechanisms, and forecast future rates of cliff retreat.

2. Background

2.1. Study Site

MDA G is located on the eastern margin of the Pajarito Plateau near Los Alamos, New Mexico (Figure 1A). The Pajarito Plateau is bounded on the west by the Jemez Mountains and on the east by the Rio Grande rift [10,11]. Surface and near-surface geology of the Pajarito Plateau is dominated by the Bandelier Tuff, a volcanic ash-flow tuff erupted in two primary pulses from the nearby Valles Caldera. The first pulse deposited the Otowi Member of the Bandelier Tuff (Qbo) ca. 1.61 Ma [12], while the second pulse deposited the multiple flow facies of the Tshirege Member of the Bandelier Tuff (Qbt) ca. 1.26 Ma [13]. Qbt caps most of the mesas on the Pajarito Plateau and is the predominant cliff-forming unit along the canyons of the plateau. All pits and shafts at the LLRW site are excavated into Qbt [14,15,16,17]. In general, mesa tops are low relief surfaces with thin, poorly developed soil covers. Immediately below the soil or exposed at the surface is one of the members of Qbt. The climate of the Pajarito Plateau is tied to its elevation, with a mean annual temperature in Los Alamos (~2100 m asl) of 9.5 °C and 440 mm precipitation, while in White Rock (~1910 m asl), the mean annual temperature is 10.2 °C with 330 mm precipitation, resulting in a Koppen classification of humid continental, however, with strong seasonality tied to the annual summer monsoon.

2.2. Geologic Setting

Cooling unit 1 (Qbt1) of the Tshirege Member is a 12–30 m thick succession of ash-flow tuffs deposited over a widespread area of the Pajarito Plateau and is divided into units 1g and 1v [18]. Unit 1v is further divided into 1v-c and 1v-u, only the latter of which is exposed at MdB. Unit 1v-u is a non-welded to moderately welded, powdery, white, vapor-phase altered unit. Pumice clasts typically make up 30% to 50% of the unit and are commonly up to 6 cm in diameter. Pumice accumulation zones are common near the top of this unit at MdB; the high concentration of pumice allows the top of this unit to weather easily in the presence of water, leading to instability in the overlying cliffs [18].

Cooling unit 2 (Qbt2) is also a succession of ash-flow tuffs, 2–24 m thick [18]. Unit 2 is a more densely welded cliff-forming unit compared to unit 1. At MdB, unit 2 is partially welded, with pumice fragments that are smaller (<2 cm) and less abundant (2% to 15%) than those observed in unit 1v-u. The base of Qbt2 is commonly marked by an abundance of pumice-containing lapilli up to 15 cm in diameter. This unit has greater phenocryst content, mainly quartz and sanidine [18], and is highly fractured, with near-vertical fractures resulting from contraction and brittle failure during cooling of the tuff, which introduced the pre-existing instability of the cliffs [19]. Failure along these cooling fractures results in large coherent simple block and cliff failures rather than wedge-type failures.

2.3. Previous Work

Previous Pajarito Plateau erosion studies focused on (1) fracture quantification and characterization [19], (2) mass wasting events [20], and (3) cosmogenic-nuclide-derived erosion rate estimates [21,22]. Wohletz [19] documented a background linear density of approximately 20 fractures per 30 m intervals, increasing to values in excess of 50 fractures per 30 m intervals in locations coincident with mapped traces of local faults. Fractures showed average strikes of either N35W or N47E, average dips between 75N and 82N, and average apertures of 0.7 cm.

Distinct differences in mass wasting events exist between the north and south rims of Pajarito Mesa, located 3 km west of MdB [20]. The north rim displays large-scale mass movement features in a zone typically 30–60 m wide, while the south rim is dominated by the failure of narrow fracture-bounded tuff blocks, with an average block thickness of 1.0–1.3 m and a maximum block thickness of 6.1 m (Figure 1). The minimum average rate of cliff retreat is 200 m Myr⁻¹, assuming that the initial incision of the canyons occurred soon after the emplacement of the Tshirege Member of the Bandelier Tuff [20]. However, Reneau and McDonald [20] hypothesized that cliff retreat typically occurs by discrete failures of varying size and time intervals; therefore, the use of average erosion rates could be misleading for short time periods.

Additionally, Poths and Goff [21] used cosmogenic Ne-21 to date three samples of Bandelier Tuff on the Pajarito Plateau: two surface samples from a mesa top and one shielded sample from a road cut. The mesa top samples yielded erosion rates of 180 and 280 m Myr⁻¹. Similarly, Albrecht et al. [22] found similar erosion rates of 50 to 110 m Myr⁻¹ using Be-10 and Al-26.

3. Materials and Methods

We focused data collection and analyses on the southern edge of MdB and the south-facing cliffs for two reasons: (1) some shafts and portions of disposal pits are located within 20 to 25 m of the southern edge of the mesa compared to >50 m from the northern edge; and (2) the southern edge is characterized by steep, near-vertical cliffs more prone to failure compared to the gently sloped northern edge. We note that the siting criteria for the disposal units required that they be a minimum of 15 m from cliff edges, although that minimum was never met.

3.1. Cliff and Fracture Characterizations

In September 2014, Lidar for the Pajarito Plateau was collected. Lidar resolution was five points per square meter, and all points were referenced to a single control point located near the Los Alamos airport. The data were post-processed to a digital elevation model (DEM) with 0.3 m resolution.

We paired the Lidar DEM with aerial photographs to map the trace of mesa top fractures and measure the distance of each fracture from the mesa edge, and then verified each fracture with field measurements. The strike, dip, aperture, and presence of fill for each fracture was also collected and entered into a geographic information system (GIS) database. We also used aerial imagery to document additional fractures that were not accessible in the field (e.g., some fractures immediately adjacent to cliff edges were unsafe to approach). In these cases, the position and orientation of fractures were captured directly within the GIS database.

We used mesa top fractures adjacent to the cliff edge to calculate fracture density and orientation relative to the cliff edge. Fracture density was calculated within ESRI ArcGIS using the Kernel Density tool with a search radius of 20 m to produce a raster with 1 m cell size. The intersection angles between each fracture and the cliff edge were also calculated in GIS; fractures that did not intersect the cliff were extrapolated to the cliff edge such that the strike of the fracture remained unchanged to calculate the intersection angle.

3.2. Surface Exposure Dating

Nine samples were analyzed for carbon-14 surface exposure dating. These samples were collected exclusively from Qbt2 because it is the uppermost unit at MdB and poses the greatest cliff failure risk. Sample locations were selected based on accessibility and lack of archaeological sites, and to sample a representative range of cliff face aspects. Samples measuring approximately 30 cm × 30 cm × 3 cm (L × W × thickness) were collected using a battery-powered angle grinder with a dry diamond blade. Sample locations were determined using hand-held GPS units, elevations derived from the Lidar-derived DEM, and topographic shielding collected using a hand-held sighting inclinometer approximating the topography.

Quartz phenocrysts were isolated by crushing and sieving to isolate the 250–500-micron-size fraction. Iron within the sanidine phenocrysts and groundmass meant that quartz was efficiently isolated using induced-roll magnetic mineral separation. Samples were then etched in a diluted mixture of hydrofluoric and nitric acids until visual purity was achieved, generally requiring four 24 h periods of etching in which the acid solution was changed, and the sample was rinsed with 18.2 MΩ water after each 24 h period.

Carbon was extracted from the samples using the Tulane University Carbon Extraction and Graphitization System following procedures outlined in [23]. Carbon isotope ratios were measured at the Lawrence Livermore National Laboratory Center for Accelerator Mass Spectrometry (Livermore, CA, USA) (¹⁴C/¹³C) and the UC-Davis Stable Isotope Facility (Davis, CA, USA) (¹³C/¹²C). The former ratios were normalized to the oxalic acid II secondary standard (NIST 4990C), while stable carbon isotope ratios were measured relative to NBS-19 (NIST RM 8544) and back calculated to the Peedee belemnite (PDB) carbon isotope standard material. Data reduction followed standard practices for in situ ¹⁴C measurements [24].

¹⁴C exposure ages were calculated using the scaling scheme of Lifton et al. [25] and the geomagnetic model of Lifton [26], the combination of which is commonly referred to as LSDn scaling. Production rates are derived from measurements of CRONUS-A intercomparison material [24,27] measured at Tulane. All reported uncertainties are at the 1-sigma level unless otherwise noted. We further assume that the measured ¹⁴C represents that accumulated in a free-air cliff face, rather than that accumulated during progressive fracture aperture increases prior to complete failure. Additionally, because the samples are from near the cliff base, there is sufficient shielding to essentially completely shield sample surfaces from significant production of ¹⁴C by muons. Thus, the ages most accurately represent the timing of block failure and exposure of the sampled location.

4. Results

The slope angle of the south-facing Qbt2 cliffs along the entire length of MdB is generally greater than 45° (Figure 2). The median length of each fracture is 9.4 m, the median distance of each fracture center point to the cliff edge is 3.4 m, and average fracture aperture is typically greater than 5 cm (Figure 3). While zones of high fracture density exist throughout the entire southern edge of MdB, fracture density tends to be highest—on the order of 10 or more fractures per 100 m²—on the southwest-facing slopes (Figure 3). Fractures exhibit variable intersection angles with the cliff edge, although clustering of fractures parallel to sub-parallel (<45°) to the cliff face are present on the southwest-facing slopes (Figure 4).

The results of the ¹⁴C surface exposure analyses are shown in Figure 5,

Table 1. Surface exposure dating results. Additional geochemical details can be found in Supplementary Material Table S1.

Sample	Latitude (dd)	Longitude (dd)	Elevation (m)	Shielding	Concentration (104 atoms g−1)	Age	Erosion Rate
						(ka)	(m Myr−1)
15-MDB-01-SS	35.82813	−106.23586	1998	0.5332	7.28 ± 0.12	3.47 ± 0.07	319 ± 6.9
15-MDB-02-SS	35.82895	−106.23636	2000	0.4463	19.20 ± 0.27	21.28 ± 1.38	10 ± 1.5
15-MDB-03-SS	35.82897	−106.23637	2000	0.4463	19.20 ± 0.26	21.29 ± 1.38	10 ± 1.5
15-MDB-04-SS	35.82862	−106.23799	1997	0.6083	16.20 ± 0.23	7.99 ± 0.19	73 ± 2.4
15-MDB-05-SS	35.82836	−106.23859	2001	0.6843	5.76 ± 0.10	2.07 ± 0.04	576 ± 12.2
15-MDB-08-SS	35.82912	−106.24004	2004	0.6371	8.69 ± 0.14	3.63 ± 0.07	281 ± 6.0
15-MDB-10-SS	35.82996	−106.24197	2010	0.8007	18.10 ± 0.25	6.75 ± 0.15	93 ± 2.6
15-MDB-12-SS	35.83016	−106.24307	2012	0.799	9.17 ± 0.14	3.15 ± 0.06	325 ± 6.6
15-MDB-14-SS	35.83113	−106.24581	2014	0.4833	4.70 ± 0.09	2.17 ± 0.05	615 ± 14.2

All samples have an assumed density of 2 g cm⁻³ and thickness of 2.5 cm.

, and Table S1. Exposure ages range from ca. 2 to 22 ka. The youngest surface exposure dates are generally associated with prominent, south- to southwest-facing cliff features (i.e., “noses”), while older surface exposure data are generally correlated with side drainages and less-prominent cliff position.

5. Discussion

5.1. Fracture Density, Orientation, and Slope Aspect Thresholds

While fracture characteristics, exposure age results, and slope analyses provide individual metrics of slope stability at MdB, we employed a multi-modal data fusion approach to synthesize the potential risk for cliff failure. We integrated slope and fracture characteristics with the exposure age results using a multivariate regression analysis to (1) identify which parameters exert the greatest control over cliff retreat using age as a surrogate and (2) predict the high-risk locations for cliff failure.

First, we used ArcGIS to convert slope and fracture characteristic datasets from vector-based datasets (points and lines) to raster-based datasets (grid cells); then, we reclassified rasters into integer datasets such that a value of 1 represents the least desirable condition and a value of 9 represents the most desirable condition (

Table 2. Raster reclassification parameters.

Raster	Condition When Integer = 1 (Least Desirable)	Condition When Integer = 9 (Most Desirable)
Fracture orientation	Fractures sub-parallel to cliff face	Fractures orthogonal to cliff face
Fracture density	High fracture density	Low fracture density
Slope aspect	S45° E to S45° W	N45° E and N45° W
Slope angle	Steep slopes (>60°)	Less steep to shallow slopes (<60°)

).

Initially, we assigned all four raster datasets listed in Table 2 an equal weight percent of importance. The correlation between the surface exposure age and the weighted overlay analysis is roughly linear, with younger dates corresponding to lower values (i.e., less desirable cliff characteristics and greater risk of failure) and older dates corresponding to higher values (i.e., more desirable cliff characteristics and less risk of failure). While there is a general correlation between cosmogenic nuclide age and the weighted overlay analysis when all variables are weighted equally, we also adjusted the weight percent of each raster to better match the cosmogenic nuclide age results using a least squares minimization (Figure 6). In essence, we aimed to maximize the correlation such that young dates (i.e., recent failures) correspond to the least desirable locations (high fracture density, etc.) and old dates (i.e., older failures) correspond to more desirable locations (low fracture density, etc.). Since no single slope or fracture characteristic correlates significantly with surface exposure ages, integrating these data and their respective weights into a single dataset is critical for identifying locations that have the greatest likelihood of failure. We used the resulting least-squares-derived weightings (

Table 3. Weights for overlay analysis. Resulting weights are not subject to requirement of summing to 1.

Fracture Orientation	Slope Angle	Aspect	Fracture Density
0.3545	0.303	0.2794	−0.0811

) to calculate at-risk locations for future cliff failures (Figure 7). Interestingly—within the current dataset—fracture density exerts the least control in terms of cliff retreat events; rather, slope angle and fracture orientation dominate our cliff stability metric.

Our analysis assumes that younger exposure ages represent the locations most at-risk for future cliff failures because their recent failure may imply greater failure recurrence. However, the opposite is also plausible—younger ages might represent locations least at-risk since they experienced a cliff failure event in geologically recent times, thus resetting their failure cycle. In the absence of additional exposure ages and/or modern-day cliff failure results, we assume younger exposure ages are indicative of higher-risk cliff locations.

5.2. Cliff Failure Recurrence Interval and Time to Facility Exposure

We are ultimately interested in the estimated time it will take for any given fracture to fail and when the cliff edge may eventually reach waste storage shafts and pits. The mesa top pits and shafts were sited no closer than 15 m (and generally greater than 30 m) from any cliff edge. Thus, we need to know the quasi-recurrence interval for cliff failure based on surface exposure ages. Here, we refer to the quasi-recurrence interval as we do not have information on the magnitude of cliff failure for a given period. Instead, we model the quasi-recurrence interval by calculating the difference in time between any two of our observed exposure ages using a bootstrap Monte Carlo approach. Doing so yields a histogram of recurrence intervals (Figure 8), showing a bimodal distribution dominated by the two oldest exposure ages from the reentrant valley on the east side of MdB. The resulting median recurrence interval is 3.7 kyr. We argue this approach has the further benefit of simulating an ability to sample a large number of sites for surface exposure ages, including additional sites within the reentrant drainages of MdB.

Next, if we assume that the distance to the nearest fracture from any given cliff edge location defines the next blockfall at that location, we can infer that the median distance from a fracture to a cliff edge will take 3.7 kyr to collapse. In the MdB case, this means that on average, a 3.38 m block will fall every 3.7 kyr. Assuming a worst-case scenario of 15 m to the nearest disposal pit or shaft from a cliff edge yields a time to disposal unit exposure of 16.4 kyr. It is important to note that the above relies on the median recurrence interval, while the use of the mean recurrence interval increases the time scale of exposure to 27 kyr. In either case, there is evidence from the exposure ages and fracture characteristics that implies generally stable cliff faces on the prominent southern edge, and even more stability along reentrant drainage edges. Further exposure dating from the cliff faces and fallen blocks would aid in the development of a probabilistic cliff failure model based on time–magnitude relationships.

6. Conclusions

At the MdB low-level waste repository, large tabular boulders on the canyon floor adjacent to the cliff face indicate that cliff retreat is dominated by episodic rockfall caused by the dislodgement of one or more fracture-bounded blocks. Such behavior thus forms the basis for future cliff failure and informs management decisions associated with the LLRW.

Direct exposure of waste through cliff failure and subsequent breaching of the disposal pits and shafts is extremely unlikely over the next 1000 years and even the next 10,000 years. Surface exposure ages provide constraints on the probability of future failures at this location. At this time, we cannot make any conclusions regarding correlation between surface exposure ages (as a surrogate for cliff failure events) and seismic events. But the prevalence of multiple cliff failure events spanning the Holocene and the state of knowledge regarding Holocene seismic event timing imply that fracture and slope characteristics [28], rather than seismic events and ground motion, are the primary driver of cliff retreat in arid environments such as the Pajarito Plateau. The weighted overlay analysis results for future cliff failure, shown in Figure 7, will be used for the practical application of informing the design of the final earthen site closure cover. For example, the information may be used to bolster cover stability or increase set back in those areas most susceptible to cliff failure.