DRC-2025-000618

~
Div of Waste Management
and Radiation Control
FEB 2 1 2025
--------ENERGYSOLUTIONS ================
February 21, 2025 CD-2025-036
Mr. Doug Hansen, Director
Division of Waste Management and Radiation Control
P.O. Box 144880
Salt Lake City, UT 84114-4880
Subject: Federal Cell Facility Application: Responses to the Director's Requests for Information
DRC-2024-004817 and DRC-2024-006804
Dear Mr. Hansen,
EnergySolutions hereby responds to the Utah Division of Waste Management and Radiation Control's
March 8, 2024 (DRC-2024-004817) and September 16, 2024 (DRC-2024-006804) Requests for
Infom1ation (RFI) on our Federal Cell Facility Application. A response is provided for each request using
the Director's assigned reference number.
Each RFI letter is addressed individually, with the text of the RFI quoted in bold italics followed by its
response. Multiple follow-up RFI letters are addressed in this response, with each letter preceded by a
solid line. The following attachments are also provided:
I. _README _ SIBERIA _FILES, a guide to electronic files provided in support of this work; and
the associated electronic files
Appendix 0:
Erosion Modeling
■ 0-2: (from DRC-2024-006804): Data from Siberia rerun necessary for confirmation of RFJ closure.
An electronic copy of the datafiles necessary to rerun and confirm the SIBERIA analysis is
provided with the compact disc attached to this response document. SIBERIA files are placed
within the following directory tree:
CD: \REFERENCES\Electronic files
Included within this directory is a file labeled _README_SIBERIA_FILES.pdfthe further
description of the location and content of each file associated with the SIBERIA analysis.
■ 0-3: (from DRC-2024-006804): Data from Siberia rerun necessary for confirmation of RFJ closure.
An electronic copy of the datafiles necessary to rerun and confirm the SIBERIA analysis is
provided with the compact disc attached to this response document. SIBERIA files are placed
within the following directory tree:
CD: \REFERENCES\Electronic files
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
[phone] • Fax: [phone] • www .energysolutions.com
Page I of77
~
ENERGY SOLUTIONS
Mr. Doug Hansen
CD-2025-037
February 21, 2025
Included within this directory is a file labeled _README_SIBERIA_FILES.pdfthe further
description of the location and content of each file associated with the SIBERIA analysis.
• 0-4: (from DRC-2024-006804): Data from Siberia rerun necessary for confirmation of RF/ closure.
An electronic copy of the datafiles necessary to rerun and confirm the SIBERIA analysis is
provided with the compact disc attached to this response document. SIBERIA files are placed
within the following directory tree:
CD: \REFERENCES\Electronic files
Included within this directory is a file labeled _README_SIBERIA_FILES.pdfthe further
description of the location and content of each file associated with the SIBERIA analysis.
• 0 -5: (from DRC-2024-006804): Data from Siberia rerun necessary for confirmation of RF/ closure.
An electronic copy of the datafiles necessary to rerun and confirm the SIBERIA analysis is
provided with the compact disc attached to this response document. SIBERIA files are placed
within the following directory tree:
CD: \REFERENCES\Electronic files
Included within this directory is a file labeled _README_SJBERIA_FILES.pdfthe further
description of the location and content of each file associated with the SIBERIA analysis.
• 0 -6: (from DRC-2024-006804): Data from Siberia rerun necessary for confirmation. of RF/ closure.
An electronic copy of the datafiles necessary to rerun and confirm the SIBERIA analysis is
provided with the compact disc attached to this response document. SIBERIA files are placed
within the following directory tree:
CD: \REFERENCES\Electronic files
Included within this directory is a file labeled _README_SIBERIA_FILES.pdfthe further
description of the location and content of each file associated with the SIBERIA analysis.
• 0-7: (from DRC-2024-006804): Data from Siberia rerun necessary for confirmation of RF/ closure.
An electronic copy of the datafiles necessary to rerun and confirm the SIBERIA analysis is
provided with the compact disc attached to this response document. SIBERIA files are placed
within the following directory tree:
CD: \REFERENCES\Electronic files
Included within this directory is a file labeled _README_SIBERIA_FILES.pdfthe further
description of the location and content of each file associated with the SIBERIA analysis.
299 South Mam Street, Suite 1700 • Sal l Lake City, Utah 84111
[phone] • Fax: [phone] • www.energysolutions.com
Page 2 of77
~
ENERGYSOLUTIONS
Mr. Doug Hansen
CD-2025-037
February 21, 2025
■ 0-8: (from DRC-2024-006804): Data from Siberia rerun necessary for confirmation of RF/ closure.
An electronic copy of the datafiles necessary to rerun and confirm the SIBERIA analysis is
provided with the compact disc attached to this response document. SIBERIA files are placed
within the following directory tree:
CD: \REFERENCES\Electronic files
Included within this directory is a file labeled _README_SIBERJA_FILES.pdfthe further
description of the location and content of each file associated with the SIBERIA analysis.
■ 0-31: (from DRC-2024-006804): Data from Siberia rerun necessary for confirmation of RF/ closure.
An electronic copy of the datafiles necessary to rerun and confirm the SIBERIA analysis is
provided with the compact disc attached to this response document. SIBERIA files are placed
within the following directory tree:
CD: \REFERENCES\Electronic files
Included within this directory is a file labeled _README_SIBERJA_FILES.pdfthe further
description of the location and content of each file associated with the SIBERIA analysis.
■ 0-32: (from DRC-2024-006804): Data from Siberia rerun necessary for confirmation of RF/ closure.
An electronic copy of the datafiles necessary to rerun and confirm the SIBERIA analysis is
provided with the compact disc attached to this response document. SIBERIA files are placed
within the following directory tree:
CD: \REFERENCES\Electronic files
Included within this directory is a file labeled _README_SIBERJA_FILES.pdfthe further
description of the location and content of each file associated with the SIBERIA analysis.
■ 0-33: (from DRC-2024-006804): Data from Siberia rerun necessary for confirmation of RF/ closure.
An electronic copy of the datafiles necessary to rerun and confirm the SIBERIA analysis is
provided with the compact disc attached to this response document. SIBERIA files are placed
within the following directory tree:
CD: \REFERENCES\Electronic files
Included within this directory is a file labeled _README_SIBERJA_FILES.pdfthe further
description of the location and content of each file associated with the SIBERIA analysis.
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
[phone] • Fax: [phone] • www.eoergysolutioos.com
Page 3 of77
~
ENERGYSOLUTJONS
Mr. Doug Hansen
CD-2025-037
February 21 , 2025
■ 0-36: (from DRC-2024-006804): Data from Siberia rerun necessary for confirmation of RFI closure.
An electronic copy of the datafiles necessary to rerun and confirm the SIBERIA analysis is
provided with the compact disc attached to this response document. SIBERIA files are placed
within the following directory tree:
CD : \REFERENCES\Electronic files
Included within this directory is a file labeled _README_SJBERIA_FILES.pdfthe further
description of the location and content of each file associated with the SIBERIA analysis.
■ 0-39: (from DRC-2024-006804): Data from Siberia rerun necessary for confirmation of RFI closure.
An electronic copy of the datafiles necessary to rerun and confirm the SIBERIA analysis is
provided with the compact disc attached to this response document. SIBERIA files are placed
within the following directory tree:
CD: \REFERENCES\Electronic files
Included within this directory is a file labeled _README_SJBERIA_FILES.pdfthe further
description of the location and content of each file associated with the SIBERIA analysis.
Appendix 0:
SIBERIA Modeling
■ 0-38.a: (from DRC-2024-004817): The long-term stability of the Federal Cell Facility (FCF) and the
requirement for ongoing maintenance after closure must be based upon analyses of active natural
processes. The analyses must provide reasonable assurance for long-term stability of the FCF, and
that upon closure of the FCF (JO CFR Part 61.13), ongoing active maintenance will not be
necessary.
Questions remain regarding the use of SIBERIA as a tool to evaluate the long-term stability/erosion
of the FCF. The processes leading to rill-gully channelization of landforms are stochastic and driven
by threshold precipitation events at the Clive site.
The SIBERIA model for the FCF uses the same precipitation magnitude in each time step, yet
gullying is a threshold dominated event driven by particularly large storm events. It appears to the
Division that the SIBERIA model is not incorporating the stochastic-threshold processes needed to
begin rill-gully formation and not channelization and may lack the resolution to capture these
events. As for calibration of the SIBERIA model, the RHEM model itself does not seem to result in
significant channeling, gullying or headward erosion, and using a RHEM model to calibrate a
SIBERIA model could lead to a poorly calibrated model and does not give realistic results. The
calibration of the model by matching sediment flux also does not capture the spatial/temporal
distribution of erosion that is important to the long-term stability of the FCF. The approach of
Tucker (Tucker 2004) using a divergence of the sediment flux might provide a more realistic model.
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
[phone] • Fax: (801 ) 880-2879 • www.eaergysolutioas.com
Page 4 of77
~
ENERGY SOLUTIONS
Mr. Doug Hansen
CD-2025-037
February 21, 2025
Given the concerns with gullying, the Division requests that Energy Solutions either modify the
SIBERIA model to account for stochastic precipitation or develop an updated model that can address
the concerns outlined.
The SIBERIA model has been modified to account for stochastic precipitation and to explore
reasonable parameter space for critical parameters. This modification allows SIBERIA to reflect
differing magnitudes of precipitation in each time step by applying a different bl value at each
time step. The differing b 1 values are based on multiple time series of individual storm events.
The modified SIBERIA model was applied both to the engineered cover design for the Federal
Cell Facility (FCF) and to an unarmored landform representing the geometry of the FCF. This
provides an opportunity to compare SIBERIA results with observed conditions for similar
unprotected slopes at the Clive site.
This response document presents erosion simulation results that follow the modeling approach
proposed in June 2024 and documented in "RFl 38a Modeling Proposal," included as Appendix 1
to this response. In addition to the parameters b I, m I, n, I, and QsHold identified in the modeling
proposal, this response considers hillslope diffusivity, dZ, as a parameter to explore across a
reasonable parameter space.
In order to address the concerns raised in RFI O-38.a, an approach containing two key features
was proposed:
1. Reevaluate the distribution of parameters in the current parameter optimization worliflow
based on literature information.
2. Stochastically modify parameters in SJBERJA at each timestep.
Each of these two key features is addressed below. The Executive Summary of this response
document provides a high-level overview of methods, results, and discussion. Sections 2.0, 3.0,
and 4.0 detail methods implemented to incorporate these two key features and corresponding
results. Finally, the following appendices accompany this work: Appendix I provides the "RFI
38a Modeling Proposal," Appendix 2 provides more detail for mathematical derivations, and
Appendix 3 provides additional figures and tables to support the analyses presented in the main
body of this work.
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
[phone] • Fax: [phone] • www.energysolutions.com
Page 5 of77
~
ENERGYSOLUTIONS
1.0 Executive Su mmary
Mr. Doug Hansen
CD-2025-037
February 21, 2025
In order to address the concerns raised in RPI O-38.a, an approach containing two key features
was proposed:
1. Reevaluate the distribution of parameters in the current parameter optimization workflow
based on literature information.
2. Stochastically modify parameters in SIBERJA at each timestep.
The first step in our work was a review of the range of parameter values from the literature.
While the SIBERIA user manual (Willgoose 2005) does provide ranges for values and some
discussion of the dependence among parameters, Neptune also considered parameter distributions
and corresponding ranges in Landlab (Barnhart et al. (2020). Ranges in Landlab were developed
using a comprehensive literature review representing a wide range of environmental and
hydrological conditions and spatial scales. While all landscape evolution models use similar
fundamental equations, parameterizations differ among codebases. Because of this, we
analytically related the parameters considered in Landlab to those used in SIBERIA. In this way
we were able to translate ranges of values from Barnhart et al. (2020) to corresponding values for
the parameters used by SIBERIA. Details of these relationships are provided in Appendix 2.
This work further inform the range of values to be used for the b I parameter in SIBERIA in a
manner that accounts for the effects of stochastic storm events. To accomplish this, the Cligen
stochastic weather generator was used to simulate 15 different 300-year time series of individual
storm events based on the USDA ARS climate file for Dugway, UT. One example of the ability
of the simulated series to depict large stochastic events consistent with the weather in the area is
the simulation of the maximum storm event. One would expect that a 300-year record (Cligen)
would represent larger extreme storm events than a 56-year record (NOAA historical data) and
this is indeed the case here. Across the 15 different 300-year times series, the range of maximum
storm events per series was 1.85 to 2.76 inches. The NOAA Dugway station has a maximum
daily precipitation event of 1.46 inches for the period 1950-2006. As a function of the use of the
300-year time frame, the stochastic storm series used to estimate values for the b 1 parameter in
this work reflect meaningfully greater extreme storm events than what has been observed in the
56-year historical record for the NOAA Dugway site.
As a means to scale daily precipitation statistics to the model space considered in this analysis,
each 300-year simulated series generated using Cligen was then used to drive erosion events in
the Rangeland Hydrology and Erosion Model (RHEM) for 1000 realizations. The equations that
relate the estimated sediment yield via advection and diffusion were optimized to estimate values
for b 1, m 1, and n 1. A different distribution of 1000 values for each of these three parameters was
developed corresponding to each of the 15 storm series. For each parameter-b 1, m 1, and n ]-
the 1000 values corresponding to each of the 15 different storm series were pooled together to
create empirical distributions for b I, m I, and n I. Each of these distributions consists of 15,000
distinct values. These empirical distributions were then subsequently used to provide values to
initialize SIBERIA runs and project potential future erosion at the Federal Cell Facility (FCF).
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
[phone] • Fax: [phone] • www.energysolutions.com
Page 6 of77
~
ENERGY SOLUTIONS
Mr. Doug Hansen
CD-2025-037
February 21 , 2025
In accordance with the modeling plan, SIBERJA was modified so that a different randomly
selected value ofbl (drawn from the empirical distribution with 15,000 values in it) was used
each year throughout the entire time period of a simulation. Prior to running a realization of
SIBERIA, a random sample of b 1 values ( of length corresponding to the number of years of the
simulation) was generated and written to a file read by the modified version of SIBERIA at the
start of a simulation. In this way, a different value ofbl is used each year by this modified
version of SIBERIA. Values form 1, n I, QsHold, and dZ were randomly sampled from the
respective empirical distributions at the beginning of the realization and held constant through
time for each realization. This was consistent with the "RFI 38a Modeling Proposal."
In order to assess the ability of SIBERIA to depict erosion given these parameter distributions
and this novel modeling approach of varying bl at each (annual) timestep, two different cases of
cover condition were subjected to erosion for 500-year time periods. The first case is a cover
condition that represents the armored engineered cover consistent with the design, materials, and
geometry of the FCF. The second case is an unarmored landform representing the geometry of the
FCF-in other words, an unprotected clay pile of identical dimensions to the FCF without gravel
in the surface layer of the top slope and without rip rap as the surface layer on the side slopes.
The evaluation of the armored engineered cover was more comprehensive in terms of the number
of realizations generated and assessed since this is the geometry being considered for impacts of
projected erosion through time. The comparison of these two cases of cover condition serves two
purposes: I) Simulations of projected erosion on the unarmored version of the FCF provide a
reality check on the ability of this approach using SIBERIA to depict erosion as observed on
similar-scale, unarmored landforms at the site, and 2) The comparison of these two cases of cover
condition provides an estimate of the benefit of the armored engineered cover with respect to
reduction in projected erosion impacts.
Multiple transects perpendicular to the primary direction of flow were assessed for erosion depth.
For both cases (i.e., the engineered cover and the unarmored landform), the most eroded transect
among those evaluated was the one 140 m down from the ridge (closest to the edge of the top
slope) on the west side (140-West). After 100 years of simulated erosion on the engineered cover
for the FCF, the median (across 10,000 simulations) 9Y" percentile depth ofa gully along the
140-West transect is projected to be 7.4 cm (2.9 inches). To clarify, this was calculated by
computing the 95th percentile of the erosion depths recorded along the transect at I 00 years and
then repeating this calculation for each of the I 0,000 simulations. The median of those I 0,000
different 95th percentile values was then calculated and reported. At 500 years of simulated
erosion on the FCF engineered cover, the median (across 10,000 simulations) 9511' percentile
depth ofa gully along the 140-West transect was projected to be 13.6 cm (5 .35 inches). For
context, this depth of 13.6 cm represents less than half the thickness of the top layer of the
engineered cover. In other words, the extreme gully depths captured by the 95th percentiles are
projected to penetrate less than 25% of the way through the two feet of soil in the cover profile
between the surface and the frost protection layer. The majority of gullies are much shallower.
For comparison, after I 00 years of simulated erosion on the unam10red landform, the median
(across 10,000 realizations) 95 th percentile depth of a gully along the 140-West transect was
projected to be 28.9 cm ( 11 .4 inches). At 500 years of simulated erosion on the unarmored
landform, the median (across I 0,000 realizations) 95th percentile depth of a gully along the
140-West transect was projected to be 113.4 cm (44.6 inches). The unarmored landform
299 So uth Main Street, Suite 1700 • Salt Lake City, Utah 84111
(801 ) 649-2000 • Fax: (801 ) 880-2879 • www.energysolutions.com
Page 7 of77
~
ENERGY SOLUTIONS
Mr. Doug Hansen
CD-2025-037
February 21, 2025
accordingly displays roughly an order of magnitude more erosion than the engineered cover
under identical modeling conditions.
The comparison of results from the two cases for the transect with the largest projected gully
depths demonstrates the impact of the engineered cover proposed for the FCF on potential erosion
in the future. The SIBERIA results also demonstrate that projected erosion will not reach a depth
that would impact the frost protection barrier within 500 years. No meaningful differences in
projected gully spacing or width due to erosion were observed between the two cover cases.
The results of this analysis are generally consistent with 500-year versions of the simulations used
in DU PA v3.0. One key structural difference between the modeling approach implemented in
this work and that of the approach used in the DU PA v3.0 is the implementation of time-varying
bl values for the results presented here. When higher values ofbl occur earlier in a simulation,
concentrated flow paths with larger drainage area are established early and erosion is amplified
along those flow paths relative to the approach where bl is constant through time.
2.0 Addressing the Two Key Features
2.1 Reevaluate the distribution of parameters in the current parameter
optimization workflow based on literature information
Proposed approach
This work consists of comparing ranges of parameters in SIBERIA to ranges in other models to
determine how the input distributions for SIBERIA should be modified. We primarily focused on
ranges of values provided in Barnhart et al. (2020). Appendix 2 of this document details the
crosswalk between the parameters presented in Barnhart et al. (2020) and those of SIBERIA. For
the sake of clarity, this document uses the parameterization presented in the SIBERIA user
manual.
Ultimately, the literature review was used as an initial starting point to develop the range of
parameters to be considered.
In accordance with the proposed approach, the following parameters I had ranges characterized
using literature review:
nl
bl
ml
QsHold
dZ
= transport-limited slope exponent
= transport-limited stream power erosion coefficient
= transport-limited drainage area exponent
= transport-limited stream power erosion threshold
= hillslope diffusivity
Parameters b 1, m 1, and n 1 are addressed together, followed by QsHold and dZ.
1 Diffusivity, dZ, was also explored as characterized in Section 2.1.3 below.
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
[phone] • Fax: [phone] • www.energysolutions.com
Page 8 of77
~
ENFRGYSOLUTIONS
2.1.1 Develop distributi ons for bl, ml, an d nl
Mr. Doug Hansen
CD-2025-037
February 21, 2025
Our previous work with SIBERIA on other sites showed that there is correlation between
parameter sets that correspond to a given magnitude of projected sediment flux. Hence, there is
some degree of exchangeability with respect to pairs of the values of parameters and the behavior
of the model in terms of projections of flux. This is particularly true for the b 1 and m 1
parameters, such that pairs of these two parameters with meaningfully different individual values
ofbl and ml can yield similar results in terms of projected erosion (i.e., equifinality).
As a consequence, if values for b 1 and m 1 are independently sampled from a wide range of initial
starting values, projections of erosion can be consistent with sites that have substantially more
precipitation and potential erosion. Ultimately, this suggests that pulling parameter values
independently from a large range of values developed solely by the literature review may result in
the application of erosion processes that do not reflect the range of actual potential erosion forces
from the site.
While varying b 1 at each timestep provides a mechanism to more realistically depict the impact
of stochastic storms on gully formation processes, given that the rest of the parameters are fixed
for each specific realization, this situation has the potential to exacerbate concerns regarding
impacts associated with equifinality. That is, realizations are more likely to select combinations
of parameters that do not correspond to the behavior of the area around Clive.
In order to assess the potential impact of pulling parameter values independently using the range
of values solely determined by Ji terature review, we went through a process of parameter
estimation informed by use of the Rangeland Hydrology and Erosion Model (RHEM) for a site
with a much wetter climate. Specifically, we generated storm series using Cligen for the climate
station parameter fi le from Mena, Arkansas. The Mena site has roughly an order of magnitude
more annual precipitation (approximately 150 cm/yr) than the 19.9 cm/year reported for the
Dugway climate station. When the approach detailed below (i.e., using Cligen-generated storm
series with the output from RHEM to optimize in the context of the advection and diffusion
equations) is implemented, the fo llowing parameter values were estimated for bl, m 1, and nl
(9.64E-4, 1.28, 0.53). Figure 1 contrasts sediment yield predicted by RHEM using the Mena, AR
storm series (blue, orange, and green lines representing sediment yield from hillslope profiles
with slope 0.02, 0.024, and 0.03, respectively) with predictions using the Dugway, UT storm
series (red lines). Volumetric sediment yield from the Mena, AR storm series is roughly an order
of magnitude greater than it is using the Dugway, UT storm series. Applying parameters based on
precipitation typical of Mena, AR to the FCF is clearly inappropriate.
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
[phone] • Fax: (801 ) 880-2879 • www.energysolutions.com
Page 9 of77
0.07
0.06
0.05
;;;-0.04 .s
ii; 0.03
0.02
0.01
0.00
0 20
~
ENERGY SOLUTIONS
40 60
area (m2)
80 100
Mr. Doug Hansen
CD-2025-037
February 21, 2025
120
Figure 1. Volumetric sediment yield predicted by RHEM using a storm series based on the Cligen climate
station parameter file for Mena, AR (blue: slope 0.02; orange: slope 0.024; green: slope 0.03) compared to the
full range of sediment yield predicted across 15,000 storm series based on the Dugway, UT climate parameters
(red). Cover conditions correspond to realization 548.
For each of the I 5,000 iterations of the Cligen-RHEM parameter estimation procedure driven by
stochastic storm series, a different vector ( e.g., b I, m 1, n 1) of parameter estimates is obtained.
Figure 2 shows a scatterplot of the 15,000 pairs of b 1 and m 1 that were obtained through this
process, for Dugway, UT. While the value ofbl from the Mena, AR analysis is well within the
range ofbl values obtained from the Cligen-RHEM parameter estimation approach for the
Dugway site, the value of m I from the Mena, AR analysis is well above the range of the largest
values.
1e-03
.0
5e-04
1.15 1.20
m
• Location
• Dugway, UT
• Mena.AR
1.25
Figure 2. Scatterplot of pairs of values of bl and ml obtained from the Cligen-RHEM optimization approach
using stochastically generated storm series from the Dugway, UT climate station (average annual precipitation
of approximately 20 cm). The pair of values obtained for the Mena, AR station (average annual precipitation of
approximately 150 cm) is placed for reference.
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
[phone] • Fax: (801 ) 880-2879 • www.energysolutions.com
Page IO of77
~
ENER.GYSOLUTIONS
Mr. Doug Hansen
CD-2025-037
February 21, 2025
The location of the (b I, m I) point that corresponds to the Mena, AR analysis could be within the
(b 1, m 1) parameter space, if this parameter space was defined solely by the use of a range of
values from the literature. This depends to some extent on where the center of the distribution is
identified to sit within that range. Using the b 1, m 1 pair that corresponds to the Mena, AR site
would grossly overestimate potential erosion at the Clive site, since this corresponds to annual
average precipitation that is an order of magnitude larger than in the region around the Clive site
(as represented by the use of the Dugway, UT station in this analysis). By using the synthetic
storm series that depicts individual events in a manner consistent with the historical record at the
Dugway, UT site to refine the parameter space for b 1 and m 1, this simulation approach depicts
the effects of stochastic storm events that reflect the historical record for the area.
It is challenging to correlate values for b I to precipitation statistics due to the lack of data linking
precipitation directly to observed erosion near the site. Such data has simply not been developed
to date. One way to link daily precipitation records to erosion is through RHEM. Hernandez et al.
(2017) presents the most recent version of RHEM, which is able to ingest a daily precipitation
record for a given area and output simulated estimates of runoff depth and erosion sediment yield.
By default, RHEM uses a 300-year daily precipitation record to predict runoff, sediment yield,
and other related quantities on an annual basis along with the magnitude of storms of selected
return periods.
In the work presented here, the stochastic nature of storm events is represented through the use of
Cligen version 5.3 (USDA ARS 2008) to generate multiple synthetic storm series that are based
on the daily historical data from the Dugway, UT climate station (station# 422257). The
following information (emphasis added) is taken directly from the Cligen website
(https://catalog.data.gov/dataset/cligen-bbee2),
"Cligen is a stochastic weather generator which produces daily estimates of
precipitation, temperature, dewpoint, wind, and solar radiation for a single geographic
point, using monthly parameters (means, SD's, skewness, etc.) derived from the historic
measurements. Unlike other climate generators, it produces individual storm parameter
estimates, including time to peak, peak intensity, and storm duration, which are
required to run the WEPP and the WEPS soil erosion models. Station parameter files to
run Cligen for several thousand U. S. sites are available for download from this website:
(https://www.ars.usda.gov/midwest-area/west-la{ayette-in/national-soil-erosion-
research/docslwepp/cligen/) also data and software to build station files for international
sites."
In the work presented here, Cligen was used to generate 15 different storm series, each of length
300 years based on the Dugway, UT station file. To clarify, each of the 15 distinct synthetic storm
series used in the results presented here depicts daily individual storm events in a manner that is
consistent with the observed data from the Dugway, UT, USDA-ARS climate station
(Ut422257.par). An example of the daily precipitation output that is aggregated from the
individual synthetic storm events depicted in Cligen from storm series G is presented in Figure 3.
In this case, Cligen provides a total of seven storm events that exceed the historical maximum for
Dugway, UT.
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
[phone] • Fax: [phone] • www.energysolutions.com
Page 11 of77
50
~
ENERGY SOLUTIONS
Mr. Doug Hansen
CD-2025-037
February 21, 2025
40 ___________ • __ • ____ -: 0 ___ -•-Hist~Ml!,!t:roum ____ .• • ________ -.• .. _________ •
10
2i)
·oo 20, JOO
Figure 3. Daily precipitation depth (mm) from synthetic storm series G generated by Cligen 5.3. The synthetic
storm series is based on information contained in the observed data from the Dugway, UT climate station
(Ut422257.par). This daily precipitation data is used as an input into the Rangeland Hydrology and Erosion
Model (RHEM). Historical maximum from the Dugway NOAA station 56-year record is shown in the horizontal
red dashed line.
We then used the results of the Cligen-RHEM approach for parameter estimation to further refine
the range and locate the center of the di stributions. This was done to incorporate the use of
precipitation statistics to calculate/scale the values of the b 1 parameter to simulate the effects of
stochastic storm events. The parameter distributions developed as a result of the combination of
the process of literature review and Cligen-RHEM parameter estimation are presented in Table 1.
Table 1. Parameters used in this modeling work.
Parameter Units Barnhart Range SIBERIA Range
et al.
(2020)
Discharge-sediment L3-2m1 /T Kt [1 .0E-04, 1.0E-02] /31 [2.73E-04,
coefficient 1.44E-03]
(transport-limited)
Discharge-sediment Unitless mt [5.2E-02, 2.4E-01] m1 [1 .12E+00,
exponent for discharge 1.23E+00]
(transport-limited)
Discharge-sediment Unitless nt [1 ,2] n1 [1.14E+00,
exponent for slope 1.46E+00]
(transport-limited)
Transport threshold L3-2m1 /T Wt,c Not reported qst, Qt, At [9.29E-08,
(transport-limited) 4.63E-04]
Diffusive coefficient on L2 /T D [1 .0E-07, 1.0E-01] dZ [3.34E-05,
slope 3.46E-04]
L indicates a unit of length; T indicates a unit of time.
In Table l , parameter units are expressed in terms oflength (L) and time (T) dimensions,
followed by symbols representing the parameters in Barnhart et al. (2020), the range of values
reported by Barnhart et al. (2020), and the symbols representing the parameters in SIBERIA. The
final column reports the range of values drawn from the empirical distributions informed by
RHEM (beta_l, m_l , n_l, and dZ) or developed by Monte Carlo sampling of variables
contributing to the sediment transport threshold as described below. The distributions developed
from the ranges in the final column of Table 1 are shown in Figure 4.
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
[phone] • Fax: [phone] • www.energysolutions.com
Page 12 of77
300000
g'200000
"' ::, C' ., it 100000
0.0005
jl
2000 !1
!1 ~ 1500 :11
C .,
~ 1000 ~ u.. 500 I
0
0.50 0 75
~
ENERGY SOLUTIONS
1000
750 ~ C (I) 500 ::,
C' e u.. 250
0.0010 0.0015
b1 Value
600 .li• ,,
!1
g 400 :I, .,
::, er it 200
1 00 1 25 1 50 -14
n1 Value
1 15 1 20
m1 Value
-12 -10
Log(QsHold Value)
jl
f1
!:
1.25
-8
Mr. Doug Hansen
CD-2025-037
February 21, 2025
Figure 4. Histograms of values of bl, ml, n 1, and QsHold obtained from tbe Cligen-RHEM optimization
approach using stochastically generated storm series from the Dugway, UT climate station (average annual
precipitation of approximately 20 cm). The value for each of the parameters that correspond to the Mena, AR
site is shown with the vertical red dashed line.
The RHEM model can ingest a synthetic storm series generated by Cligen and project erosion for
a region with a defined slope, length, and foliar and ground cover conditions. Foliar cover is
specified as an area proportion of plant communities defined as bunch grass, sod grass, forb, and
shrub. Ground cover is specified as an area proportion of rock, litter, plant basal stem, and
biocrust. One thousand realizations of foliar and ground cover conditions were drawn from
statistical distributions as described in Section 4.2 of Neptune (2023). These were combined with
the 15 Cligen storm series to generate a set of 15,000 reference datasets of sediment yield varying
with area and slope by RHEM. Mass-based RHEM sediment yield output is converted to
volumetric sediment flux through bulk density. An example of this is plotted in Figure 5 for one
randomly selected cover conditions realization. Application of the 15 different stochastic storm
series has the effect of broadening the distribution of sediment yield compared with the Cligen
default random seed.
0.00100
0.00115
0.001~0
0.00175
0.00100
0.00075
O.OD050
0.00015
0.0:+000
' 70 T <O ' 60 M ' lC.D :20
Figure 5. Sediment yield predicted by RHEM for 15 different stochastic storm series applied to hillslope profiles
with the same cover conditions. Results from the Cligen default random seed are plotted in black. Blue: slope
0.02; Green: slope 0.024; Orange: slope 0.03. Cover realization 305, units are m3.
299 South Main Street, Suite 1700 • Sall Lake City, Utah 84111
[phone] • Fax: [phone] • www .energysolutions.com
Page 13 of77
~
ENERGYSOLUT/ONS
Mr. Doug Hansen
CD-2025-037
February 21, 2025
A curve-fitting analysis was implemented to evaluate the relationships between drainage area
(D _ area), slope gradient (gslp ), and sediment transport (sym3) using advective and diffusion
equations (Will goose et al. 1991 a, 1991 b ). In this analysis, D _ area (calculated as slope length,
assuming a I-meter-wide profile) and slope gradient serve as the inputs to the advective and
diffusion equations. Slope length (D _area) was analyzed at levels of 5, I 0, 20, 30, 40, 50, 60, 70,
80, 90, 100, 110, and 120 square meters, and slope gradient was evaluated at 2.0, 2.4, and 3.0
percent. All 39 unique combinations of D _ area and gslp were weighted equally with respect to
their impact on the optimization to obtain estimates of the parameters of interest (i.e., b 1, m 1, and
nl).
The curve-fitting process uses nonlinear least-squares optimization with the Levenberg-
Marquardt algorithm to minimize the residual sum of squares (RSS) between two different sets of
sym3 values. The first set of sym3 values represent sediment transport simulated by RHEM based
on combinations of drainage area (D _ area) and slope gradient (gslp ). The second set of sym3
values are computed using the advective and diffusion models using the parameters b 1, m 1, n 1,
and dZ. This optimization estimates the parameters of the advective and diffusion models. For the
advection model, the parameters include b I, m 1, and n 1, which describe how sediment transport
scales with slope length and slope gradient.
The diffusion model extends this by adding a parameter, dZ, to account for sediment diffusion
effects. This analysis ensures consistency in parameter estimates across all combinations of slope
length, slope gradient, and draws. The approach used here estimates the parameters for each
combination of series and draws. The quality of the match between the calculated curve and the
set of sediment yields predicted by RHEM is evaluated by calculating the norm of the residual
vector for both the advective and diffusion models. Finally, we compare the fitted parameters and
residual norms across series and draws to ensure consistency among results .
This analysis was implemented using the curve_fit function from the scipy.optimize library in
Python, although a Nelder-Mead (Press et al. (2007), Section I 0.5) algorithm implemented in the
minimize function in the scipy.optimize module returns the same values.
The advective model is mathematically expressed as:
sym3 = bl • (D_area)m1 • (gslp)n1
The diffusion model extends this by including the diffusion parameter dZ and is expressed as:
sym3 = bl • (D_area)m1 • (glsp)TT1 + glsp • dZ
In these models the sediment yield predicted by RHEM in units Mg ha-1 y-1 are converted to area
by assuming hillslope profiles of 1 meter width and soil bulk density of 1.3 6 Mg m-3.
2.1.2 Developing the Distribution for QsHold
The distribution for sediment flux threshold values was developed from the equation for
transport-limited sediment flux threshold in Barnhart et al. (2020), equation 30.
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
[phone] • Fax: [phone] • www.energysolutions.com
Page 14 of77
~
ENERGY SOLUTIONS
Mr. Doug Hansen
CD-2025-037
February 21, 2025
A Monte Carlo scheme is used to sample values for parameters describing the underlying
sediment properties and hydrological scaling relationships2. On the basis of flow per unit width,
the formula is:
Where
Pw is density of water [M L·3]
R b is nondimensional buoyant density [-]
g is the gravitational constant [L·1 T 2]
Ds is sediment grain size [L]
r; is critical Shields stress [-]
and the nondimensional buoyant density is defined as:
Where
Ps is density of sediment [M L·3]
After multiplying by channel width with a dependence on flow, the formula becomes:
Where
Cw is the discharge-channel width scaling exponent [-]
kw is the discharge-channel width coefficient [L3-2cwrcw]
Cq is the discharge-drainage area scaling exponent[-]
kq is the discharge-drainage area coefficient[ L 1-3cq r-1]
A is drainage area [L2]
R b is nondimensional buoyant density [-]
g is the gravitational constant [L·1 T 2]
Ds is sediment grain size [L]
r; is critical Shields stress [-]
These parameters were sampled across ranges reported in Barnhart et al. (2020) and listed in
Table 2.
2 Monte Carlo simulation was based on unifom1ly sampling parameter value ranges given in Barnhart et
al. (2020).
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
(80 I) 649-2000 • Fax : (80 I) 880-2879 • www.energysolutions.com
Page JS of77
~
ENERGYSOLUTIONS
Mr. Doug Hansen
CD-2025-037
February 21, 2025
Table 2. Value ranges for parameters in the transport-limited sediment threshold formula of
Barnhart et al. (2020), equation 30.
Parameter Units Symbol Range
Discharge-channel width Unitless Cw (0 .37, 0.62]
scaling exponent
Discharge-channel width * kw (5.5, 9.1]
coefficient
Discharge-drainage area Unitless Cq (0.7, 1.0]
scaling exponent
Discharge-drainage area * kq [1 .0E-03, 30]
coefficient
Area m A (0, 120]
Nondimensional buoyant Unitless Rb 1.65
density
Gravitational constant m-1 s·2 g 9.81
Sediment grain size m Ds [6.0E-06, 8.0E-
06]
Critical Shields stress Unitless r * C (0 .3,0.8]
Density of water Mg m·3 Pw 1.0
Density of sediment Mg m·3 Ps 2.65
* units for a given realization depend on values drawn for Cw and Cq as described in Table 1.
QsHold values were initially assessed by Monte Carlo simulation in which the range of values of
constituent parameters given in Table 2 were sampled. Multiple distributions of QsHold values
were obtained by applying selected discrete values of drainage area. Median values from these
distributions were an order of magnitude smaller than the median value of sediment fluxes
predicted by RHEM for all of the 15 ,000 5-meter profiles. This is illustrated in Figure 6. These
small values for QsHold are consistent with a simplifying step described in Barnhart et al. (2020)
which requires the threshold to be small compared to the fluvial sediment flux. Therefore, the
QsHold SIBERIA parameter di stribution was developed using a Monte Carlo sampling approach.
10000~-------------------
60000
50000
~ 40000
] 30000
20000
10000
0.00000
,, ,,
'' I I I I
I I
' ' ' ' ' ' ' ' ' ' I I
' ' I I ,"" ... ,,
I I I \ I I I \ I I I \
I I f \ I I I \ ' \, ,: \..~ __ ,,, .... ---------·. -.............. __ _
0 00005 0.00010 0.00015 0.00020 0.00025
OsHold I SY (m')
RHEM, all storms
SY(m3)
5 m2
••••••• 10m2
••••••• 20m2
Monte Carlo QsHold
--60m2
--80m2
--100 m2
--12om2
Figure 6. Comparison of sediment flux thresholds calculated for selected drainage areas folJowing BARNHART
(solid lines) with sediment yields predicted by RHEM from 15 stochastic storm series and 1,000 cover conditions
realizations (dashed lines).
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
[phone] • Fax: [phone] • www.energysolutions.com
Page 16 of77
~
ENERGYSOLUTIONS
Mr. Doug Hansen
CD-2025-037
February 21, 2025
Larger values of QsHold were assessed as described in Rogers (2023 b) and resulted in slightly
increased incision into the top slope of the Federal Cell design over a very small areal proportion
through I 0,000 years. In that analysis, values of the threshold were based on RHEM-modeled
average sediment yield from I 0-meter profiles with cover conditions predicting maximum
erosion. That value was 0. I 2 Mg ha·1, equivalent to 1.2E-04 Mg off I 0-meter profiles or
equivalently 8.8E-03 m3. The mass value was assigned as QsHold input to SIBERIA. Test cases
were run by multiplying by factors of 10 along with the zero-threshold case. The effects reported
from that analysis are enumerated in Table 3-3 in Appendix 3 and are plotted in Figure 7.
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
[phone] • Fax: [phone] • www.energysolutions.com
Page 17 of77
g
g ;;;
-0,4
-0.5
.~ -0.6
" e
-0.7
-0.8
0.000
-0 025
-0 050
g
~ -0.075 ~
; -0.100
'6 ~
-0.125
-0.150
-0.075
-0.100
-0.125
' >,
~ -0.150
t:
I;; -0.175
-0.200
-0.225
~
E ERGYSOLUTIONS
•---------------------------------------1,
10-1 10-~
OSHold (Tl
I \ I I I I I I
\ I I •
,.
' ' I I ' ,,
,' ' ' ' ' ' I '
-· ,' ,,' ,'
·---------------------------... -----,' I
•----------------------------•
10-1 10-~
OSHold (Tl
' I I I ,'
,'
I / ' ,,, ... , ... •
,.
' I I I I I
I f' I I I I I I I I
I I I I I I
I I A I ,, I
_........ I •---------------------------...--,'
·---------------------------~
10-1 10-5
OSHold (Tl
I I
I
,//.
-· ,,'
10-3
a.
b.
C.
Mr. Doug Hansen
CD-2025-037
February 21 , 2025
Figure 7. Effect on erosion across the top slope as QsHold increases from zero to 0.12 Mg as reported in Rogers
(2023h). Zero threshold is represented by lE-10. Black: 2023 RHEM-based parameters. Brown: 2021 analog site
parameters (Erosion v2.0). a. Maximum incision into the top slope. b. Median incision. c. Average annual
sediment yield.
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
[phone] • Fax: [phone] • www.energysolutions.com
Page 18 of77
~
ENERGYSOLUTIONS
Mr. Doug Hansen
CD-2025-037
February 21, 2025
Maximum incision into the top slope appears insensitive to the sediment threshold until
exceeding 0.001 Mg for parameters developed on the lower gradients of the top slope. For
parameters incorporating higher slopes as the 2021 analog site simulations used (0.02-0.16), the
break point is 0.0 I Mg.
The decrease in median incision and sediment yield accompanying increased maximum incision
depicted in Figure 7 indicates erosion is highly concentrated into channels as higher threshold
values further suppress hillslope erosion relative to channel erosion. While threshold effects may
act to preserve older terrain along high-order channels at the catchment scale (Tucker 2004), for
realistic values of QsHold there is little relative effect across terrain with low-order channels
along a catchment divide such as the FCF.
2.1.3 Diffusivity
Barnhart et al. (2020) report a range of values of a Soil Creep Rate Coefficient3 D between 10-7
and I 0-1 m2/y. This range across six orders of magnitude was derived from a comprehensive
literature review and addresses studies across a wide range of geomorphic settings. There is no
site-specific data available to help constrain values of dZ.
Hancock et al. (2000), working with runoff and sediment data from a 50-year-old existing waste
rock mound, estimated dZ for SIBERIA simulations by trial and error adjustments to match the
observed terrain. Slopes involved were greater than 0.5, which is much steeper than the FCF side
slope of 0.2 and over an order of magnitude greater than the FCF top slope of 0.024. The final
value of dZ was not reported in Hancock et al. (2000). On a different waste rock mound
constructed 5-8 years prior to the study, diffusion was found to contribute little to evolved terrain
(Evans et al. 2000).
Fleming and Johnson (I 975) report soil creep rates between 0.008 and 0.0 I 3 m2 y·1 at the top of a
1.5-meter-thick silty clay soil layer in Menlo Park, CA based on tiltmeter measurements.
However, this was on slopes of 0.12 0.14, again a poor match for the FCF top slope of 0.024 and
its gravel-amended clay surface. The soil in this study had both a higher proportion of sand than
the material proposed for the FCF top slope and its clay was predominantly montmorillonite
which exhibits shrink-swell behavior. Nothing otherwise directly applicable to the top slope
material was found during a literature review.
Previous SIBERIA modeling based on RHEM varied advective parameters while holding dZ
constant at I E-04 m2 y·1 (Neptune 2023). Rogers (2023a) reported the effects on a domain with
0.5-meter grid cells using dZ values of IE-05, IE-04 (default), and lE-03 m2 y·1 while
maintaining advective parameters corresponding to simulations representing the minimum,
median, and maximum predicted erosion based on the 1,000 cover conditions realizations (see
Appendix 3, Table 3-4).
Effects of using these values on the median erosion case with a 1.5-meter resolution DEM were
also reported. In that analysis, average annual sediment yield increased as dZ increased for a
3 The term "D" in Barnhart et al. (2020) corresponds to the "dZ" parameter in SIBERIA, and to "Dz" in
comments on the modeling proposal.
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
[phone] • Fax: [phone] • www.eoergysolutions.com
Page 19 of77
~
ENERGYSOLUTIONS
Mr. Doug Hansen
CD-2025-037
February 21 , 2025
given set of advective parameters while incision into lower cover layers on the top slope
increased by a minor amount. Similar results were reported using these values of dZ on the
original 1.5-meter DEM for simulations derived from cover conditions representing the median
predicted erosion case (see Appendix 3, Table 3-5).
Results from the entire suite of tests augmented with subsequent tests based on the maximum
predicted erosion parameters and a 1.5-meter DEM with random roughness of zero mean and
0.002 m standard deviation are presented in Figure 8. Appendix 3 includes tabulated results in
Table 3-5 . Area proportions were recalculated to reflect an extended region of interest on the top
slope (12.93 ha compared to 12.78 ha, sampling closer to the edges).
M I
-0.25
-0.50
-0.75
M';:,.., -1.00
I 11)
-C
I--1.25
>-Vl
-1.50
-1.75
-2.00
10-7 10-6 10-5 10-4 10-3 10-2 10-1
dZ (m2y-l)
Figure 8. Sediment yield (SY) at different values of dZ, in SIBERIA, realization d725 from Rogers (2023b).
Note that Cover Layer depth intervals were redefined with higher resolution after the response of
Rogers (2023a). Variation in sediment yield and area proportion of cover layer incision in Table
3-5 are quantified on the basis of the revised depth intervals.
Table 3-5 indicates that by 10,000 years the maximum incision depth on the top slope was
insensitive to dZ until exceeding 1.0·3 m2/y with an apparent threshold between 1.0·3 m2/y and I .o-
2 m2/y. Maximum incision depth and average annual sediment yield increased rapidly with
increasing dZ above the threshold.
The measure used to inform the Probabilistic PA Model-area proportion incised into different
cover layers-also rapidly degraded when dZ increased above IE-03 m2 y·1. At IE-03 m2 y·1 the
Frost Protection Layer (Cover Layer 5) at depth 0.61 m was exposed over two of 57,459 model
domain cells; at IE-02 m2 y·1 the Frost Protection Layer was exposed over 1244 of the cells while
Page 20 of77
299 South Maio Street, Suite 1700 • Salt Lake City, Utah 84111
(8 01 ) 649-2000 • Fax: [phone] • www.eoergysolutions.com
~
ENERGY SOLUTIONS
Mr. Doug Hansen
CD-2025-0 37
February 21, 2025
incision advanced into Cover Layer IO as erosion occurred at a rate that tripled sediment yield
over the I 0,000-year period.
RHEM simulations exhibit diffusive effects as positive offsets from the x-axis in plots of
sediment flux at small areas, proportional to slope (Figure 9). Estimates of dZ were obtained by
nonlinear curve fitting as described in Section 2.1.1. A maximum profile length of 30 meters was
used, selected by trial and error to separate the diffusive and advective components. The resulting
range of dZ was between 3.4E-05 and 3.SE-04 m2 y·1.
0.00012
0.00010
0.00008
1 0.00006
1i;
0.00004
0.00002
0.00000
0
0.0006
0.0005
0.0004
1 0.0003
1i;
0.0002
0.0001
0.0000
0 10 15
Area (m'l
20
Area (m'l
6
25 30
10
a.
35 40
b.
Figure 9. Reference sediment yields modeled by RHEM as in Figure 4, viewed across smaller drainage areas to
highlight diffusive effects. Results are from 15 stochastic storm series on hill slopes with the same foliar and
ground cover conditions. Results from the Cligen default random seed are plotted in black. Blue: slope 0.02;
Green: slope 0.024; Orange: slope 0.03. a. 5 -10-meter segment exhibiting diffusive effects in the predictions.
b. 5 -40 meter segment with combined fluvial and diffusive effects apparent. Cover realization 305, units m3.
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
[phone] • Fax: (801 ) 880-2879 • www.energysolutions.com
Page 21 of77
~
E ERGYSOLUTIONS
Mr. Doug Hansen
CD-2025-037
February 21, 2025
The increased erosion predicted by RHEM using the Mena, AR climate station parameter file
displayed a significantly larger diffusion effect compared to that predicted using the Dugway, UT
parameters in Cligen as shown in Figure I 0. The best-fit estimate of dZ for the Mena, AR-based
RHEM sediment yields was 0.01 m2 y-1, compared to a mean value of 2.38E-04 m2 y-1 for erosion
based on Dugway storm series. The empirical distribution obtained using the 15 Dugway UT
storm series was sampled for the I 0,000 simulations in this work.
0.0175
0.0150
0.0125
0.0100
1 0.0075
Iii
0.0050
0.0025
0.0000
-0,0025
0 5 10 15 20
area (m2)
25 30 35 40
Figure 10. Sediment yield predictions by RHEM for cover conditions realization 548 for profiles up to 40 meters
long. Mena, AR Cligen cl imate parameters blue: slope 0.02; orange: slope 0_024; green: slope 0.03. Red: Full
range of sediment yields across slopes 0.02-0.03 using Dugway, UT Cligen climate parameters.
299 South Main Street, Suite 1700 • Salt Lake City , Utah 84 111
(80 I) 649-2000 • Fax: (80 I) 880-2879 • www.energysolutions.com
Page 22 of77
~
E ERGYSOLUTIONS
Mr. Doug Hansen
CD-2025-037
February 21 , 2025
Figure 11 displays the distribution used in this work, with the Mena, AR valu e as a dashed red
line.
10000
7500
>. (.)
C
Q) 5000 ::J CT ~ LL
2500
0
0.0000 0.0025 0.0050 0.0075 0.0 00
Diffusivity (D)
Figure 11. Histograms of values of diffusivity obtained from the Cligen-RHEM optimization approach using
stochastically generated storm series from the Dugway, UT climate station (average annual precipitation of
approximately 20 cm). The value for the parameter that corresponds to the Mena, AR site (average annual
precipitation of approximately 150 cm) is shown with the vertical red dashed line.
2.2 Stochastically modify parameters in SIBERIA at each timestep
The SIBERIA source code was modified to read values for b 1 for each simulation year from a
text file. That is, for a SIBERIA run covering a 500-year time span, there is a vector of b 1 values
that has length 500. The elements of this vector are independent random samples from the
distribution ofbl shown in Figure 4 above. For a given element in a vector, that value is used to
inform SIBERIA for the corresponding year in the simulation. That is, the 24th element of the b 1
vector will be the va lue that SIBERIA uses for the 24th year of a given run.
Since there are likely interactive effects on the simulation of projected erosion among the b I, n I,
m 1, dZ, and QsHold parameters, I 0,000 simulations with randomly selected n 1, m I , dZ, and
QsHold parameters were performed. As outlined in the "RF] 38a Modeling Proposal," these
parameters did not change at each timestep and instead were specified through random selection
from their respective distributions at the beginning of each realization.
Given this, a high-level representation of the steps implemented for a single SIBERIA realization
for 500 years is:
I) Generate a random sample of length 500 from the distribution of bl. Each element of this
vector is the value that inform s SIBERIA in a given timestep.
2) Select a single random value from the distribution ofnl.
3) Select a single random value from the distribution of QsH old.
4) Select a single random value from the distribution of ml.
5) Select a single random value from the distribution of dZ.
299 South Main Street, Suite 1700 • Salt Lake City , Utah 84111
(801 ) 649-2000 • Fax: (801 ) 880-2879 • www.energysolutions.com
Page 23 of77
~
ENERGYSOLUTIONS
Mr. Doug Hansen
CD-2025-037
February 21, 2025
In summary, each of the 10,000 realizations of SIBERIA have distinct, randomly generated
vectors of values for b I. The length of the vector containing values for b I is defined by the
duration of the simulation (e.g., 500 years). The other four parameters remain constant within
each realization, and vary across the 10,000 realizations. An example of a 500-year b I time series
is shown in Figure 12.
00012
0.0010
..
100 200 .IOO 400 500
year
Figure J 2. 500-year time series of bl values applied as realization 1800. Black: ]-year bl value. Orange: J 0-year
moving window median value.
Little difference would be observed between simulated terrain evolved by the default constant b 1
value and landforms evolved using a stochastic b 1 time series generated with the same mean
value, except when m 1 appreciably differs from 1. Since m 1 ranges between 1.12 and 1.24,
stochastic b 1 time series can enhance erosion when higher values of b I occur earlier in the
simulation. Earlier establishment of concentrated flow paths may then erode more quickly than in
the default case during timesteps with lower b I values. The converse may be true to an extent as
well, although higher b 1 values later in the simulation may be sufficient to enhance erosion in the
channel network established by the time they occur.
2.2.1 Engineered cover
In order to assess the ability of SIBERIA to depict erosion given these parameter distributions
and the approach of varying b 1 at each (annual) timestep, two different cases of cover condition
were subjected to erosion for 500-year time periods. The first case is a cover condition that
represents the armored engineered cover consistent with the design, materials, and geometry of
the FCF. The second case is an unarmored landform representing the geometry of the FCF-in
other words, an unprotected clay pile without gravel in the surface layer of the top slope and
without rip rap as the surface layer on the side slopes. The evaluation of the armored engineered
cover was more comprehensive in terms of the number of realizations assessed since this is the
geometry being evaluated for potential erosion impacts through time. The comparison of these
two cover conditions serves two purposes: 1) Simulations of projected erosion on the unarmored
version of the FCF provide a reality check on the ability of the distributions used to depict erosion
at the site, and 2) The comparison of these two cases provides an estimate of the impact of the
armoring with respect to reductions in projected erosion impacts.
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
[phone] • Fax: [phone] • www.energysolutions.com
Page 24 of77
~
ENERGYSOLUTIONS
Mr. Doug Hansen
CD-2025-037
February 21, 2025
The embankment footprint and engineered cover design, dimensions, and grades modeled are
identical to that evaluated in NAC-0017 _R7, Erosion Modeling for the Clive DU PA (Neptune
2023), with design details provided in NAC-0019 _RS, Embankment Modeling for the Clive DU
PA (Neptune 2021).
To test the stochastic b I version of SIBERIA, four different I 0,000-year simulations were
performed using parameter values based on the default Cligen storm series in which days with
less than 0.3 mm precipitation had been removed. Three scenarios were randomly selected from
the set of 1,000 cover conditions realizations-d0725, d 1576, and d4315. Their original 500-
element vectors ofb I values were extended by randomly selected an additional 9500 values with
replacement from the initial set of curve-fit results.
Realization d5326 was added to run as a 10,000-year test since it had the highest mean b I value
over 500 years. For this simulation the 9500 values were sampled with replacement from the
initial 500-element vector to preserve the mean value. This more closely matched a conventional
SIBERIA simulation based on the equivalent b I value of an annual effective geomorphic event.
Each of these simulations ran to completion at 10,000 years.
2.2.2 Unarmored landform
An unarmored landform simulation is configured so that the entire model domain uses the same
b 1 vector used in the corresponding engineered cover simulation. In this context, the four 10,000-
year simulations described for the engineered cover were performed concurrently to provide
contrasting long-term results. This was done concurrently with the 4 different 10,000-year
simulations using the engineered cover configuration. Each of these simulations also ran to
completion at I 0,000 years.
After the final set of I 0,000 500-year simulations were completed on the engineered cover model
domain, a set of l 00 of these scenarios was randomly selected to run with the unarmored
landform configuration. This was done to provide a basis for statistical comparison with the
engineered cover configuration.
3.0Results
As contemplated in the "RFI 38a Modeling Proposal," results are evaluated in terms of gully
depth as well as gully width and spacing. These runs are not directly comparable with those used
in DU PA v3.0, since the DU PA uses l 0,000-year results and the present cases were run for 500
years. Instead, a set of I 00 realizations from those used in the DU PA v3.0 were rerun to obtain
500-year output so equivalent metrics could be compared.
3.1 Gully Depth
The SIBERIA model was used to simulate projected erosion on the FCF over a period of 500
years. Projections of gully depth were generated by assessing the surface heights across multiple
transects of the top slope area for each realization. Confidence intervals were calculated using
quantiles to represent the variability in incision depth metrics across all simulations for a given
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
[phone] • Fax: [phone] • www.energysolutions.com
Page 25 of77
~
E ERGYSOLUTIONS
Mr. Doug Hansen
CD-2025-037
February 21 , 2025
cover case. Intervals were constructed to reflect the plausible range of outcomes based on the
simulations. To achieve this, the 2.5th and 97 .5th percentiles of the simulated depth metrics are
calculated.
The transects were placed along the west face of the top slope at di stances approximately I 0, 30,
50, 75, 100, 120, and 140 meters from the central ridge to capture channel development across
that widest surface (labeled I0TSW, 30TSW, ... , 120TSW, and 140TSW). Figure 13 shows the
transects in planfonn, overlaid onto delta-Z at 500 years for realization 1800. This realization was
randomly selected from the set of I 00 realizations run with the unannored Iandfonn
configuration. These transects were placed below the central ridge with a similar length of 208.5
m to capture effects of flow along the longest hill slope profiles on the top slope. An additional
different transect (TSE) was placed on the east face at a distance from the edge approximately
consistent with the 140W transect; roughly l I 7 meters from the central ridge and just above the
edge with the same length as the transects on the west face. Transects were placed on the north
(TSN) and south (TSS) faces of the top slope at a similar distance from the edge, with equal
length of 232.5 meters extending across the faces .
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
(80 I) 649-2000 • Fax: [phone] • www .~Dergysolutions.com
Page 26 of77
~
E ERGYSOLUTIONS
400
350
300
250
200
LU (f) f-
150 s: s: s: s:
(f) (f) (f) (/) f-f-f-f-l{) 000 ,-... l{) (") .,....
100
TSS
so
0
0 so 100 150 200
..J en en .. w
250 300
Mr. Doug Hansen
CD-2025-037
February 21 , 2025
0.100
0.075
0.050
0.025
0.000
-0.025
-0.050
-0.075
-0.100
Figure 13. Transect locations for gully width and depth analysis, overlaid onto delta-Z at 500 years for
realization 1800. Arrows indicate direction of plotting. The color scale of the elevation change raster is
normalized to interval 1-0.1, 0.1) meters but erosion and deposition may exceed those boundaries.
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
(801 ) 649-2000 • Fax: [phone] • www.energysolutions.com
Page 27 of77
a.
b.
~
ENKR.GYSOLUTIONS
Mr. Doug Hansen
CD-2025-037
February 21 , 2025
Elevation changes across the transects on the west face of the top slope are plotted in Figure 14.
This figure presents a series of paired cross sections of elevation change due to gullying for the
engineered structure, with sections A and C, representing approximately 10 and 140 meters from
the ridgeline on the top slope west (TSW), respectively. For the unarmored landform, sections B
and D represent the same transects. Each color in the figure represents a different year within the
500-year analysis. Positive values for elevation change represent areas of deposition; and
negative values represent erosion. Note that the vertical scale is exaggerated several orders of
magnitude relative to the horizontal scale, in order to depict the complete transect on one figure
while still showing a sense of the magnitude and spacing of gullies. Note also that the scale for
the engineered cover is smaller within each pair of plots than it is for the unarmored landform ,
which shows deeper incision at each transect on the unarmored landform.
!:;~~IL--,--Wl-~-~. _)'T~_fr_"'~'1'-~-~~~----'i
~ :::: I
..:.i -0.005 :,
-0.010
_ 0.4
E ~ 0.2 ,;; UJ
0,0
~ ~ = =
I
o ~ ~ ~o 2~
distance (m)
100 200 300 400 ~o
::f jl'----,--"flll~r,n ~~~"rn~~------r-~-'~~~N~------,---'~u I
100 200 300 400 ~o
100 200 300 400
distance (m)
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
[phone] • Fax: [phone] • www.energysolutions.com
~o
Page 28 of77
C.
d.
~
ENERGY SOLUTIONS
Mr. Doug Hansen
CD-2025-037
February 21 , 2025
! ~jl..______,___"'li_r,ri~n~~r~m ~m,,_,M~, 11~~,m_,1
''"'----.----'''·. I
•001 '_.,,
!:i -0.50
~ -0.75
-1.00
100
100
200 300 400 500
distance (ml
200 JOO 400 500
distance (m)
Figure 14. Elevation change along transects on the west face of the top slope on the engineered cover and
unarmored landform. The top frame in each subplot shows elevation change on the engineered cover. Where
present the middle frame repeats elevation change on the engineered cover with they range scaled to match that
of elevation change on the unarmored landform. The bottom frame shows elevation change on the unarmored
landform. a. Transect I0TSW. b. transect 140TSW. c. Transect WSSU. d. Transect WSSL. Blue: year 5;
Orange: year 50; Green: year 100; Red: year 300; Purple: year 500.
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
[phone] • Fax: [phone] • www.energysolut.ions.com
Page 29 of77
~
ENERGY SOLUTIONS
Table 3 and Table 4 summarize key metrics for each timestep:
Mr. Doug Hansen
CD-2025-037
February 21, 2025
• Median 95th Percentile (95th): Represents the median of the 95th percentile incision
depth across realizations, capturing extreme values in erosion depth.
• Mean of Means (MoM): The average of mean depths across realizations, showing the
central tendency of erosion.
• Median of Medians (Med-Med): Represents the central tendency of median values
across realizations, providing a robust measure of typical erosion depth.
• Upper and Lower Confidence Intervals: Calculated for both mean and median depths,
these intervals quantify the range of plausible erosion outcomes.
For the engineered cover condition proposed for the FCF, 15,000 realizations were analyzed.
Summary statistics of elevation change were calculated for each of the seven transects on the
west face of the top slope, in addition to one transect from each of the north, east, and south faces
of the top slope. The transects from the north, east, and south top slopes are the same distance
from the shoulder as the 140m west slope transect (140TSW). Additionally, 2 west side slope
transects were analyzed (WSSU and WSSL).
Table 3 highlights the impact of the engineered cover on projected erosion. The deepest gullies
form on the west top slope at the 140m transect, though similar to the transect placed in a
comparable location on the east slope. With the proposed engineered cover, the maximum 95th
percentile gully depth on the top slope at 500 years is 13.6 cm (roughly 5.5 inches).
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
[phone] • Fax: [phone] • www.energysolutions.com
Page 30 of77
~
ENERGY SOLUTIONS
Mr. Doug Hansen
CD-2025-037
February 2 I, 2025
Table 3. Summary statistics from projected erosion of gully depths (cm) for the engineered cover
for the FCF.
Time
500
500
500
500
500
500
500
500
500
500
500
500
Transect 95th MoM UCL mean LCL mean Med-Med UCL med LCL med
10m (west) -0.72 -0.35 -0.62 -0.17 -0.31 -0.59 -0 .17
30m (west) -1 .36 -0.41 -0.76 -0.2 -0.32 -0.62 -0 .17
50m (west) -1 .7 -0.53 -0.99 -0 .25 -0.39 -0.79 -0.21
75m (west) -1.84 -0.62 -1 .16 -0 .31 -0.48 -1 .0 -0.29
100m (west) -2.44 -0.78 -1.49 -0.3 -0.43 -0.89 -0.2
120m (west) -4 .59 -1 .32 -2.05 -0.72 -0 .67 -1 .06 -0.39
140m (west) -13.6 -4.36 -5.24 -3.51 -1 .99 -2.8 -1 .29
upper side slope -0 .97 9.53 5.07 13.3 5.83 2.72 9.77
(west)
upper side slope -1.5 0.57 0.3 0.83 0.13 0.07 0.27
(west)
east -13.5 -4.66 -5.58 -3.71 -2 .32 -3.16 -1 .71
north -9.77 -3.62 -4.75 -2 .58 -2.36 -3.16 -1.47
south -12.1 -4.0 -5.12 -2 .79 -2 .1 -3.17 -1 .37
In contrast to the engineered cover case, I 00 realizations were analyzed for the unarmored
landfonn cover condition. The same seven transects on the west top slope, along with one
transect from each of the east, north, and south top slopes, were extracted for both the unarmored
and engineered cover conditions. Similarly, the 2 west side slope transects were extracted for both
cover conditions.
The unarmored landform results provide context for the engineered cover results. Table 4
summarizes the results from the I 00 realizations of unarmored landform simulations. Complete
tables with centennial results up to 500 years are provided in Appendix 3. Without the engineered
cover, the 95th percentile gully depth on the top slope at 500 years is estimated to be 113.4 cm
(roughly 44 inches). This is approximately an order of magnitude larger than the same statistic for
the engineered cover, demonstrating the effectiveness of the engineered cover in reducing
eros10n.
299 South Main Street, Suite 1700 • Salt Lake City, Utah 8411 I
(801 ) 649-2000 • Fax: [phone] • www .energysolutioas.com
Page 31 of77
~
ENERGYSOLUTIONS
Mr. Doug Hansen
CD-2025-037
February 21, 2025
Table 4. Summary statistics from projected erosion of gully depths ( cm) for unarmored landform
cover condition.
Time
500
500
500
500
500
500
500
500
500
500
500
500
Transect 95th MoM UCL mean LCL mean Med-Med UCL med LCL med
10m (west) -0.71 -0.33 -0.48 -0.25 -0.31 -0.46 -0.24
30m (west) -1.38 -0.4 -0.58 -0.29 -0.32 -0.49 -0.27
50m (west) -1.68 -0.51 -0.74 -0 .36 -0.38 -0.57 -0.3
75m (west) -1 .84 -0.6 -0.87 -0.44 -0.48 -0.73 -0.38
100m (west) -3.91 -1.03 -1.72 -0.62 -0.43 -0.73 -0.34
120m (west) -27.1 -5.76 -9.63 -3.6 -1 .26 -2.4 -0.79
140m (west) -113.4 -32.1 -40.1 -27.8 -12.7 -16.5 -10.7
upper side -115.6 -28.4 -33.6 -25.6 -12.7 -15.3 -11.4
slope (west)
lower side -58.5 -14.1 -1 7.0 -12.4 -6 .78 -8.14 -5.84
slope (west)
east -110.1 -27.6 -33.8 -23.5 -11.0 -14.6 -8 .7
north -70.3 -18.7 -23.3 -15.8 -7.10 -9.51 -5.93
south -84.7 -19 .6 -24.0 -16.8 -7.72 -10.4 -6.28
Figure 15 shows the evolution of the 95th percentile gully depths through time across the seven
transects sampled on the west face as well as the east, north, and south transects4. Results are
shown for both cases of the engineered cover proposed for the FCF and the unarmored landform.
In addition to showing the differences between the two cases for the covers, Figure 15 also shows
that, for the engineered cover, the rate of change of erosion through time starts to stabilize around
year 300. Initially high rates of erosion are characteristic of many LEMs as initial roughness in
the DEM is smoothed away and/or fine soil and sediment particles are preferentially removed
leaving less mobile material subject to the remainder of the simulation (Hancock et al. 2011 ;
Hancock and Willgoose 2021).
4 Engineered cover results reflect I 0,000 realizations; unarmored landform results reflect I 00 realizations.
Page 32 of77
299 South Main Street, Suite 1700 • Salt Lake City, Utah 8411 I
[phone] • Fax: [phone] • www.energysolutions.com
~
ENERGY SOLUTIONS
10m (west) 30m (west) 50rn (west)
Mr. Doug Hansen
CD-2025-037
February 21, 2025
75m (west)
00 ~ -0 5
-1 0
·1 5
00
-0 2
0.4
0.6
00
-0 5
, 0
00
1.0
1.5
100m (wes1) 120m (wes1) 140rn(west) WSSU
:~-1:~~~~~ I 2 00 -
-3 -20 90 90
·120
WSSI east norlh 80Ulh
=~~~~~i ~ -60
Time (years)
Figure 15. Time series of the 95th percentiJe gully depths for the different transects sampled along the top slope
as welJ as 2 transects from the west side slope. Engineered cover is blue; unarmored landform is red.
For transects across the majority of the top slope area (i .e., through I 00m west of the crest), the
engineered cover and unarmored landforms show similar gully depths. As the transect approaches
the break between the top slope and side slopes, results begin to diverge dramatically.
3.2 Gully Width and Spacing
Gully boundaries were identified as contiguous points where the depth exceeded the 75th
percentile of erosion for each timestep, transect, and scenario. The depth measurements were
recorded at 1.5-meter intervals along each transect, which informed the calculation of gully width
and spacing. For each gully, the width was calculated as the distance between its starting and
ending points, adjusted by 1.5 meters to account for the spacing intervals. Gully spacing was
determined as the distance between the end of one gully and the start of the next along a transect.
In general, there were negligible differences found between the two cover cases for statistics that
characterize the spacing and width of gullies. These differences may be attributed to the
methodology used to define and measure gully boundaries. Specifically:
• The 75th percentile threshold was calculated for each combination oftimestep, transect, and
scenario. This dynamic thresholding approach adjusts to the distribution of erosion depths for
each scenario, which may reduce the sensitivity of the results to the differences in cover type.
• Measurements were taken at 1.5-meter intervals, representing the grid resolution of the
analysis. While this resolution captures the general pattern of erosion, it may inherently
smooth smaller-scale variations in gully formation, particularly in cases where the spacing or
width of gullies is near the resolution limit.
299 South Main Street, Suite 1700 • Sall Lake City, Utah 84111
[phone] • Fax: [phone] • www .energysolutions.com
Page 33 of77
~
ENERGY SOLUTIONS
Mr. Doug Hansen
CD-2025-037
February 21 , 2025
Figure 16 and Figure 175 show the mean and 95u, percentile through time for projections of gull y
spacing between the two cover types for the transects along the west face. Figure 18 and Figure
19 show the mean and 95th percentile through time for projections of gully width between the two
cover types for th e transects along the west face.
10m (wesl) 30m (west) 50m (west) 75m (weal)
0 § ~ g ~ ~ 0 § ~ g ~ ~ 0 § ~ g ~ ~ 0 § ~ g ~ ~
Time (years)
Figure J 6. Time series of the mean gully spacing for the different transects sampled along the top slope and west
side slopes. Engineered cover is blue; unarmored landform is red. The behavior between the two cover cases
converges through time for top slope transects.
5 Engineered cover results reflect I 0,000 realizations; unarmored land form resu lts reflect I 00 realizations.
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
[phone] • Fax: [phone] • www.energysolutions.com
Page34 of77
200
175
i 150
" E 12.5
rs
50
25
100m (west}
WSSL
t\_
0 ~ ~ ~
30
25
20
15
10
~ §
~
ENERGY SOLUTIONS
1?0M(WH I) 140m fWHI)
Time (years)
wssu
Mr. Doug Hansen
CD-2025-037
February 21 , 2025
Figure 17. Time series of the 95th percentile of gully spacing for the different transects sampled along the top
slope and west side slopes. Engineered cover is blue; unarmored landform is red.
10m(wMt) 30m (west) 50m (west) 75m (west)
35
34
3'
32
31
30
100m (west) 120m(west) 140m (weso wssu
36
35
f? 3' .!! " E 33
3'
WSSL east north 90Ulh
3 20
3 15
3 10
3 05
3 00
Time (years)
Figure 18. Time series of the mean of gully width for the different transects sampled along the top slope and
west side slopes. Engineered cover is blue; unarmored landform is red.
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
[phone] • Fax: [phone] • www .energysolutions.com
Page 35 of77
,om (we,o
r-
100m (weal)
5 75
S 50 L ~ 525 .,
~ 5.00
4 75
• 50
WSSL
~
ENERGY SOLUTIONS
30m (west) 50m (west)
•• I ~ 4 2
40
36
120m (weal) 140m (west)
45
55 ~ .a so
35 0
4 0
30
ea!i-t noM
4 5
4 0
35
30
Mr. Doug Hansen
CD-2025-037
February 21, 2025
75,m (west)
~
wssu
~
south
4 5 F 4 50 I L :P-4 25 56 •o
4 00 52
35
3 75 4.8
30 3 50
o~~~~g o~~g~~
lime (years)
Figure 19. Time series of the 95th percentiJe of gully width for the different transects sampled along the top
slope and west side slopes. Engineered cover is blue; unarmored landform is red.
3.3Comparison with DU PA v3.0
In DU PA v3 .0, I 0,000-year simulati ons were used to construct a lookup table in the probabilistic
PA model. Each of 1,000 simulations corresponded to one row of the lookup table. Each row
related the proportional area of the top slope eroded into cover layers defined by depth intervals.
Since these 500-year simulations cannot be directly compared to the 10,000-year results, a set of
I 00 realizations were rerun to 500 years to obtain equi valent metrics. Summary statistics of
incision depth and sediment yield are compared in Table 5, while the 500-year version of the
lookup table is presented in Table 6.
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
(801 ) 649-2000 • Fax: (80 I) 880-2879 • www.eoergysolutioos.com
Page 36 of77
~
ENERGYSOLUTIONS
Mr. Doug Hansen
CD-2025-037
February 21, 2025
Table 5. Cumulative distributions of maximum incision and average annual sediment yield
comparison for 500-year output from stochastic bl simulations and 500 randomly selected
scenarios with DU PA v3.0 parameters.
max incision (m) Annual SY (Mg ha-1)
quantile DU PA v3.0 Stoch. b1 DU PA v3.0 Stoch. b1
0 0.24 0.20 0.33 0.13
0.05 0.26 0.25 0.34 0.20
0.1 0.26 0.25 0.35 0.22
0.15 0.26 0.26 0.36 0.23
0.2 0.26 0.26 0.36 0.24
0.25 0.27 0.26 0.37 0.25
0.3 0.27 0.26 0.37 0.26
0.35 0.28 0.26 0.37 0.27
0.4 0.28 0.27 0.38 0.28
0.45 0.28 0.27 0.38 0.29
0.5 0.28 0.27 0.38 0.29
0.55 0.28 0.27 0.39 0.30
0.6 0.29 0.27 0.39 0.31
0.65 0.29 0.27 0.39 0.32
0.7 0.29 0.28 0.39 0.33
0.75 0.29 0.28 0.40 0.34
0.8 0.30 0.28 0.40 0.35
0.85 0.30 0.28 0.40 0.37
0.9 0.31 0.28 0.41 0.39
0.95 0.32 0.29 0.42 0.42
1 0.35 0.33 0.44 0.70
At 500 years the current simulations project generally lower average annual sediment yield and
shallower channel incision than values obtained from the 100 randomly selected DU PA v3.0
simulations as listed in Table 5 and plotted in Figure 20.
299 South Main Street, Suite I 700 • Salt Lake City, Utah 84111
[phone] • Fax: [phone] • www.energysolutions.com
Page 37 of77
1.0
o.e
0.6
0.4
0.2
0.0
0.20 0.22
1.0
o.e
0.6
0.4
0.2
0.0
0.2
~
E ERGYSOLUTIONS
0.24 0.26 0.28 0.30
lmax incision! (m)
0.3 0.4 0.5
ISYI (Mg na-1 y-1)
0.32 0.34
0.6 0.7
Mr. Doug Hansen
CD-2025-037
February 21 , 2025
Figure 20. Cumulative distributions of absolute value of maximum incision and average annual sediment yield
listed in Table 5. Blue: DU PA v3.0 parameters (100 randomly selected instances). Orange: 2024 stochastic b]
simulations (complete set of 10,000).
In the lookup table format at 500 years the I 0,000 2024 simulations exhibit a significantly lower
proportion of top slope area incised into cover layer 2 and cover layer 3, the deepest layer
exposed by 500 years. An increase in net deposition across the top slope is apparent as well
(cover layer 0). Table 6 illustrates this by proportional area lookup table values across the 10,000
realizations of the current work (2024) and the I 00 simulations representing DU PA v3 .0 results
at 500 years.
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
(8 01) 649-2000 • Fax: (801 ) 880-2879 • www .energysolutions.com
Page 38 of77
::::::=-
ENERGYSOLUTJONS
Mr. Doug Hansen
CD-2025-037
February 21, 2025
Table 6. Comparison of proportions of top slope area incised into cover layers on the top slope at
500 years between current work ("Stoch. bl") and DU PA v3.0 parameters.
version 10 11 12 13 14 15
min Stoch. B1 0.010825 3.07E-01 0.107451 0 0 0
DU PA v3.0 0.048835 1.69E-01 0.297064 0 0 0
25th pctile Stoch. B1 0.135784 0.543044 2.28E-01 0 0 0
DU PA v3.0 0.125072 5.19E-01 0.326028 0 0 0
median Stoch. B1 0.162377 0.565734 0.271768 0 0 0
DU PA v3.0 0.131 381 0.525209 0.342905 0 0 0
75th pctile Stoch. B1 0.186359 0.584926 3.21 E-01 0 0 0
max
DU PA v3.0 0.139691 0.53335 3.55E-01 0 0 0
Stoch. B1 0.28288 0.618302 0.681843 3.48E-05 0 0
DU PA v3.0 0.155885 0.547051 0.765903 0.016168 0 0
The decrease was determined to derive from different curve fitting methods used to discover
parameter values for b 1, m 1, and n 1. The current work used nonlinear curve fitting implemented
by the Python scipy curve_fit function to obtain values for the advective parameters as described
in Section 2.1.1. The values obtained provided excellent matches to the reference sediment flux
datasets generated by RHEM.
In DU PA v3.0 the curve fitting algorithm was based on a Nelder-Mead method (Press et al.
(2007), Section 10.5) which fit average annual sediment fluxes obtained from
SIBERJA/landscape evolution model simulations to the reference sediment yields modeled by
RHEM. Following Smith and Benson (2016), each iteration of the Nelder-Mead method was
based on a 60-year simulation -the minimum duration suggested by Nearing (2004) for
characterizing erosion characteristics of a site using a continuous model such as WEPP.
For a given set ofreference sediment flux values three different curve-fitting or optimization
algorithms converge to a solution with the same values form I and n I-nonlinear curve-fitting
with the Python curve_fit function, Nelder-Mead optimization implemented in the Python
minimize function, and the Nelder-Mead implementation driving SIBERIA simulations
mentioned above. All three methods consistently characterize the scaling relationship between
sediment yield, drainage area, and slope.
The two Python functions return the same best matching value for the rate coefficient b 1 which
differs from the SIBERIA-Nelder-Mead result. The 60-year average annual sediment flux
calculation incorporates high erosion rates at the start of the simulation as initial roughness in the
DEM is smoothed (Hancock et al. 2011) transitioning to lower erosion rates toward the end of the
period. Since the 60-year average flux rate tends lower than direct RHEM projections the DU PA
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
[phone] • Fax: [phone] • www.energysolutions.com
Page39 of77
~
ENFRGYSOLUTIONS
Mr. Doug Hansen
CD-2025-037
February 21, 2025
v3.0 Nelder-Mead algorithm raises the value ofbl relative to the direct-fit value obtained by the
Python routines match the reference sediment yields. In the current work the direct-fit Dugway-
based values ofbl are approximately 26 percent lower than those reported in DU PA v3.0.
For comparison, Figure 21 shows the parameter distributions used in DU PA v3.0 alongside the
distributions from the current work. The b I distribution, while slightly wider in the current
analysis, reflects scaling with precipitation and the inclusion of 14 additional storm series beyond
the DU PA v3.0 parameterization. Similarly, distributions of m I and n I are widened slightly, as
they too are based on the same expanded precipitation record.
1 00 1 00
i 075 i 075
C C Q) Q)
O 050 O 050 'O 'O Q) Q) ~ 025 ~ 025
Cl) Cl)
0 00 0 00
-3.6 -3 4 -3 2 -3 0 -2 8 115 1 20 1 25
Log10(b1 Value) m1 Value
1 00 !!1 1.00 !i i i , 1 1 ~ 0 75 ~I -~ 0.75 !1 ~ ::El Cl) ::E 1 Q) C:
I C1> O 050 O 050 I
'O ,,
Q) C1> ~ 0.25 'iii
Cl) ~ 0.25
0 00 0.00 I
0 50 0.75 1.00 1 25 1 50 -6 -4 -2
n1 Value Log10(QsHold Value)
Figure 21. Parameter distributions comparing DU PA v3.0 with current work. DU PA v3.0 is orange, current
distributions are blue. While QsHold was set to zero in DU PA v3.0, discrete values tested in Rogers (2023b) are
shown with the vertical green dashed lines. The value for each of the parameters that correspond to the Mena,
AR site are shown with the vertical red dashed line.
4.0Summary
In order to address the concerns raised in RFI 38a, an approach containing two key features was
proposed:
1. Reevaluate the distribution of parameters in the current parameter optimization worliflow
based on literature iriformation.
2. Stochastically modify parameters in SJBERJA at each timestep.
Two different cases of cover condition were subjected to erosion for 500-year time periods using
the SIBERIA landscape evolution model. The first case is a cover condition that represents the
armored engineered cover consistent with the design, materials, and geometry of the FCF. The
second case is an unam10red landform representing the geometry of the FCF -in other words, an
unprotected clay pile without gravel in the surface layer of the top slope and without rip rap as the
surface layer on the side slopes. The evaluation of the armored engineered cover was more
299 South Main Street, Su.ite 1700 • Salt Lake City, Utah 84 111
[phone] • Fax: [phone] • www.energysolutions.com
Page 40 of77
;::::::=-
E ERGYSOLUTIONS
Mr. Doug Hansen
CD-2025-037
February 21, 2025
comprehensive in terms of the number of realizations generated and assessed since this is the
geometry being considered for impacts of projected erosion through time. The comparison of
these two cases of cover condition serves two purposes: I) Simulations of projected erosion on
the unarmored version of the FCF provide a reality check on the ability of this approach using
SIBERIA to depict erosion as observed on similar-scale, unarmored landforms at the site, and
2) The comparison of these two cases of cover condition provides an estimate of the benefit of the
armored engineered cover with respect to reduction in projected erosion impacts.
Multiple transects perpendicular to the primary direction of flow were assessed for erosion depth.
For both cases (i.e., the engineered cover and the unarmored landform), the most eroded transect
among those evaluated was the one 140 m down from the ridge (closest to the edge of the top
slope) on the west side ( I 40-W est). After I 00 years of simulated erosion on the engineered cover
for the FCF, the median (across I 0,000 simulations) 95th percentile depth of a gully along the
140-West transect is projected to be 7.4 cm (2.9 inches). To clarify, this was calculated by
computing the 95th percentile of the erosion depths recorded along the transect at I 00 years and
then repeating this calculation for each of the I 0,000 simulations. The median of those I 0,000
different 95th percentile values was then calculated and reported. At 500 years of simulated
erosion on the FCF engineered cover, the median (across I 0,000 simulations) 95th percentile
depth of a gully along the 140-West transect was projected to be I 3.6 cm (5.35 inches). For
context, this depth of 13.6 cm represents less than half the thickness of the top layer of the
engineered cover. In other words, the extreme gully depths captured by the 95th percentiles are
projected to penetrate less than 25% of the way through the two feet of soil in the cover profile
between the surface and the frost protection layer. The majority of gullies are much shallower.
For comparison, after I 00 years of simulated erosion on the unarmored landform, the median
(across I 0,000 realizations) 95 th percentile depth of a gully along the I 40-West transect was
projected to be 28.9 cm (11.4 inches). At 500 years of simulated erosion on the unarmored
landform, the median (across 10,000 realizations) 95t1, percentile depth of a gully along the
140-West transect was projected to be I 13.4 cm (44.6 inches). The unarmored landform
accordingly displays roughly an order of magnitude more erosion than the engineered cover
under identical modeling conditions.
The comparison of results from the two cases for the transect with the largest projected gully
depths demonstrates the impact of the engineered cover proposed for the FCF on potential erosion
in the future. The SIBERIA results also demonstrate that projected erosion will not reach a depth
that would impact the frost protection barrier within 500 years. No meaningful differences in
projected gully spacing or width due to erosion were observed between the two cover cases.
The results of this analysis are generally consistent with 500-year versions of the simulations used
in DU PA v3.0. One key structural difference between the modeling approach implemented in
this work and that of the approach used in the DU PA v3.0 is the implementation of time-varying
bl values for the results presented here. When higher values ofbl occur earlier in a simulation,
concentrated flow paths with larger drainage area are established early and erosion is amplified
along those flow paths relative to the approach where b 1 is constant through time.
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
(801 ) 649-2000 • Fax : [phone] • www.energysolutions.com
Page 41 of77
~
ENFRGYSOLUTIONS
Appendix 1: RFI 38a Modeling Proposal
Mr. Doug Hansen
CD-2025-037
February 21, 2025
[NOTE: This appendix is a direct quote of a pre-existing document from June 2024.)
The purpose of this document is to outline a proposed modeling approach that is intended to
address concerns raised in RFI 38a. For clarity the text of the RFI 38a is provided below.
• The long-term stability of the Federal Cell Facility (FCF) and the requirement
for ongoing maintenance after closure must be based upon analyses of active
natural processes. The analyses must provide reasonable assurance for long
term stability of the FCF, and that upon closure of the FCF (10 CFR Part 61.13),
ongoing active maintenance will not be necessary.
• Questions remain regarding the use of SIBERJA as a tool to evaluate the long-
term stability/erosion of the FCF. The processes leading to rill-gully
channelization of landforms are stochastic and driven by threshold precipitation
events at the Clive site.
• The SIBERJA model for the FCF uses the same precipitation magnitude in each
time step, yet gullying is a threshold dominated event driven by particularly large
storm events. It appears to the Division that th e SIBERJA model is not
incorporating the stochastic-threshold processes needed to begin rill-gully
formation and not channelization and may lack the resolution to capture these
events. As for calibration of the SJBERJA model, the RHEM model itself does not
seem to result in significant channeling, gullying or headward erosion, and using
a RHEM model to calibrate a SIBERJA model could lead to a poorly calibrated
model and does not give realistic results. The calibration of the model by
matching sediment flux also does not capture the spatial/temporal distribution of
erosion that is important to the long-term stability of the FCF. The approach of
Tucker (Tucker 2004) using a divergence of the sediment flux might provide a
more realistic model.
• Given the concerns with gullying, the Division requests that Energy Solutions
either modify the SIBERJA model to account for stochastic precipitation or
develop an updated model that can address the concerns outlined.
The approach presented here is intended to reflect the content of a discussion with representatives
from UDEQ, SC&A, EnergySolutions, and Neptune on 3 June 2024. In this discussion, options
were considered and explored to ensure that the path forward for the development of a proposed
modeling approach did not contain any fatal flaws that would prohibit the ability for a proposed
modeling approach to meaningfully address the concerns raised in RFI38a. It is understood that
while the details of this proposed modeling approach await review and approval by a subject
matter expert, Dr. Brian Yanites, subcontracted by SC&A, there was agreement among the parties
on the call on 3 June 2024 that, in principle, the two components of the high-level approach
proposed collectively have the potential to meaningfully address the concerns raised in RFI 38a.
In order to address the concerns raised in RFI 38a, an approach is proposed below that has two
key features:
• Reevaluate the distribution of parameters in the current parameter optimization workflow
based on literature information.
• Stochasticall y modify parameters in SIBERIA at each timestep
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
[phone] • Fax: [phone] • www.energysolutions.com
Page 43 of77
~
ENERGY SOLUTIONS
Mr. Doug Hansen
CD-2025-037
February 21 , 2025
1) Reevaluate the distribution of parameters in the current parameter optimization workflow
based on literature information.
Proposed approach
This work will consist of comparing ranges of parameters in SIBERIA to ranges in other models
to determine how the input distributions for SIBERIA should be modified.
Parameters below are represented with their designation from the SIBERIA manual.
nl
bl
ml
QsHold
= transport-limited slope exponent
= transport-limited stream power erosion coefficient
= transport-limited drainage area exponent;
= transport-limited stream power erosion threshold
For each of the parameters nl, bl, ml, QsHold, a range of parameter values will be determined
using a literature review. Probability distributions will be fit to these values such that lower and
upper quantiles will be constrained with the smallest and largest values from review of the
literature. Samples from these probability distributions will then be used to initialize a set of runs
in SIBERIA.
2) Stochastically modify parameters in SIBERIA at each timestep
We will develop a script that allows for the b 1 parameter in SIBERIA to vary at each annual
timestep for the duration of a given run. That is , for a SIBERIA run covering a 500 year time
span, there would be a vector of bl values that has length 500. The elements of this vector will be
independent random samples from the distribution of bl. For a given element in a vector, that
value will be used to inform SIBERIA for the next year in the simulation. That is, the 24th
element of the b 1 vector will be the value that SIBERIA uses for the 24th year of a given run.
Impacts of varying nl, ml, and QsHold will be explored as well. Unlike bl these parameters
will not change at each timestep and will instead be specified through random selection from their
respective distributions at the beginning of each run. Conceptually, there does not seem to be a
justification for a conceptual model where the physical characteristics that these parameters are
meant to represent would vary through time.
Given this, a high level representation of the steps needed for a single SIBERIA run for X years
would be:
1) Generate a random sample of length X from the distribution of b 1. Each element of this
vector will be the value that informs SIBERIA in a given timestep.
2) Select a single random value from the distribution ofnl.
3) Select a single random value from the distribution of QsHold.
4) Select a single random value from the distribution ofml.
299 South Main Street, Suite 1 700 • Sall Lake City, Utah 84111
[phone] • Fax: (80 I) 880-2879 • www.energysolutions.com
Page 44 of77
~
ENERGY SOLUTIONS
Mr. Doug Hansen
CD-2025-037
February 21, 2025
The ability to assess the impact of the uncertainty in bl within a single run will be a function of
the length of the time period being simulated. For longer time periods, much of the impact from
varying bl will be seen within a given simulation. Since there may be interactive effects on
simulation output among the bl, nl, ml, and QsHold parameters, 10,000 simulations with
randomly selected nl, ml, and QsHold parameters will be performed as well. Each of those
10,000 runs of SIBERIA will have distinct, randomly generated vectors of values for bl. The
length of the vector containing values for in bl will be defined by the duration of the simulation
(e.g., 500 years).
Results
Results will be presented at multiple time steps for transects drawn across the cover. This will
include depiction of results for gully depths and width through time as well as the horizontal
spacing. Particular emphasis will be placed on the first 500 years. As such, projected gully
geometries will be presented at centennial timesteps. Reasonable erosion maxima as represented
by the features of the gullies described above will be characterized using an 95th UCL across the
suite of 1,000 simulation runs that are performed.
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
[phone] • Fax : [phone] • www.energysolutions.com
Page 45 of77
~
E ERGYSOLUTJONS
Mr. Doug Hansen
CD-2025-037
February 2 1, 2025
Appendix 2: Relationships between parameters in Barnhart et al. (2020) and SIBERIA
The parameters used in the SIBERIA runs were determined based on mathematical values from
literature, constrained by model output from RHEM version 2.3.
Notation
This section li sts variables used in this writeup. We used the notation from Barnhart et al. (2020)
for most variables. Some variables are only found in the SIBERIA user manual (Willgoose 2005),
in which case we used that notation.
From Barnhart et al. (2020)
K -Erodibili ty coefficient
q* -Einstein bed load number (dimensionless)
r * -Nondimensional boundary shear stress (Shields stress)
T -Dimensional boundary shear stress
T~ -Critical Shields stress
Pw -Density of water
Rb -Nondim ensional buoyant density
g -The acceleration due to gravity
D5 -Grain size
q -Discharge per unit width of flow
Q -Geomorphically effective discharge
W -Channel width (meters)
S -Slope, defined where downwards is positive (gradient).
Sthreshold -Slope threshold for diffusion calculations
Cw -Positive constant relating geomorphically effective discharge to channel width
Cq -Positive constant relating drainage area to geomorphically effective discharge
kq -Positive constant relating drainage area to geomorphically effective discharge
A -Drainage area
kw -Positive constant relating geomorphically effective discharge to channel width
q5 -Volumetric sediment discharge per unit width of flow (S 5.1 .2)
Q5 -Total volumetric sediment flux
k1 -Constant coefficient corresponding to a friction formula; helps decompose shear stress (r)
a -Exponent on area corresponding to a friction formula ; helps decompose shear stress (r)
(unitless)
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
[phone] • Fax : [phone] • www.energysolutions.com
Page 46 of77
::::::=-
ENERGYSOLUTIONS
Mr. Doug Hansen
CD-2025-037
February 21, 2025
/3 -Exponent on slope corresponding to a friction formula; helps decompose shear stress (r)
p -Exponent used for the power law of sediment transport
D -Diffusion rate constant
Dzn -Exponent of nonlinearity for diffusion (unitless)
Dzt -Threshold for diffusion (same units as q5)
From SIBERIA User's Manual
z -elevation (positive upwards)
x, y -Horizontal directions
t -Time
U -Tectonic uplift per time
Pb -Bulk density of soil
q5 -Sediment transport (by mass)
D -diffusivity
Y -Defines whether the channel point is a channel (Y close to 1) or a hi II slope (Y close to 0)
w -Channel width
S -Slope
q -Discharge per unit width
/31 -Coefficient of equation relating sediment discharge to water discharge and slope (variable in
space)
m1 -Exponent of water discharge in equation relating sediment discharge to water discharge and
slope (not variable in space or time)
n1 -Exponent of slope in equation relating sediment discharge to water discharge and slope (not
variable in space or time)
/33 -Coefficient in equation relating discharge to area (variable in space)
m 3 -Exponent of area in equation relating discharge to area (not variable in space or time)
/34 -Coefficient in equation relating width to discharge (variable in space)
m 4 -Exponent of discharge in equation relating width to discharge (not variable in space or time)
/3 -Coefficient in the sediment runoff equation
m -Exponent of area in sediment runoff equation
K -Erodibility
Ot -Reduction factor in sediment transport in hillslopes vs. channels
Qc -Channel discharge
Qt -Threshold below which fluvial sediment transport does not occur (default 0, in SIBERIA as
QsHold)
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
[phone] • Fax: [phone] • www.energysolutions.com
Page 47 of77
~
ENERGYSOLUT/ONS
Ranges of Parameters from Barnhart et al. (2020)
Mr. Doug Hansen
CD-2025-037
February 21 , 2025
To determine the general ranges for parameters, literature values from Barnhart et al. (2020) were
used. The parameters used in SIBERIA (Will goose 2011) were derived in this document from
quantities from which ra nges of values from the literature are known. Simulations can then be
used to give a plausible range of values for the parameters that are used in SIBERIA. These
ranges were used for the SIBERIA parameter QsHold in our simulations.
The following equations are found in Barnhart et al. (2020), based on the Meyer-Peter Muller
formula of Wong and Parker (2006):
From these equations, the following are derived:
q Q
w
Barnhart et al. (2020) provide the dimensional fom1 as:
Using the approximati on r :::: k1qa sP , and multiplying by W to get Qs, it is derived that
Treating the second term in the parentheses as negligible and substituting the expressions for q
and W into the equation results in:
Thus, following Barnhart et al. (2020), an equation of the fom1 Q5 = KtAmt5nt has been derived,
where
299 South Main Street, Suite 1700 • Salt Lake City, Utah 8411 I
(801 ) 649-2000 • Fax: [phone] • www.energysolutioos.com
Page 48 of77
;::::::=-
E ERGYSOLUTIONS
3.97 kiw+l.Sa(l-Cw) k}.S
Kt= P1s Rbgk1sa-1
mt = CqCw + l.Sacq(l -Cw)
nt = l.Sf].
Mr. Doug Hansen
CD-2025-037
February 21, 2025
This expression is used in simulations to understand the range that these coefficients may have in
SIBERIA. The simulated values of Kt, mt, and nt were compared to the predictions for fl, m , and
n1 , respectively, from RHEM (see the section "Moving from RHEM to SIBERIA" below).
The only parameter used explicitly from Barnhart et al. (2020) was the threshold QsHold
(notated as Wt c in Barnhart et al. (2020)), which was found to be
In SIBERIA, this value will be internally scaled for the modeling.
Values per Unit Width
These expressions can also be found for the discharge per unit width. Not multiplying by W and
ignoring the threshold as negligible yields
Thus, getting the equations into the proper form yields
qs = Kt,qAmt,q5nt,q
3.97 k}.S k~.Sa(l-cw)
Kt,q = --p-~-_s_R_b_g_k_~-.s-a-
mt,q = l.Sacq (1 -Cw)
nt,q = l.Sf].
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
(80 1) 649-2000 • Fax: [phone] • www.energysolutions.com
Page 49 of77
~
ENERGY SOLUTIONS
The threshold tenn for the by-width value can be similarly derived as
Wt,q,c
Cligen
Mr. Doug Hansen
CD-2025-037
February 21 , 2025
Cligen v5.3 (USDA ARS 2008) was used to generate additional climate parameters. Fifteen stonn
series for 300 years were created for Dugway, UT using the following seeds:
• the default seed (no command line argument specifying a seed),
• 469 and 666, two testing seeds, and
• 12 additional seeds generated from numpy.random.default_ mg(seed=412).
Moving from RHEM to SIBERIA
In RHEM (USDA ARS SWRC 2015), 15 climate scenarios generated in Cligen were run with 13
different lengths (5, I 0, 20, 30, 40, 50, 60, 70, 80, 90, I 00, 110, and 120 m) and three different
slopes (2%, 2.4%, and 3%). The width was set to I meter in all runs of this experiment. These
scenarios are run I 000 times each for a total of 39,000 runs for each climate scenario. This also
results in 15,000 datasets with 39 rows each that can each be used for one fit of the parameters to
be passed to the SIBERIA modeling. These datasets are used to fit the SIBERIA parameters bl,
ml, nl, and dZ. The SIBERIA parameters b3 and m3 were set to 1 because this investigation is
focused on erosion, not hydrology.
The governing equation for SIBERIA is
oz 1
-= U + -"ii· q + D"i/2 z. ot Pb s
This equation explains landscape elevation change over time (oz/ot) as a function of tectonic
uplift (U), the sediment loss to flow ("i/ • q5 / Pb), and the sediment loss to diffusion (D"i/2 z). The
tectonic uplift is ignored in the sequel because it does not contribute to sediment flow.
From these outputs, the Levenberg-Marquardt algorithm was fit to find the least squares
predictors for the parameters in the equations
299 South Main StTeet, Suite 1700 • Salt Lake City, Utah 84111
[phone] • Fax: (801 ) 880-2879 • www.energysolutions.com
Page SO of77
:;:::::::::=-
ENERGYSOLUTIONS
Mr. Doug Hansen
CD-2025-037
February 21 , 2025
From the SIBERIA user manual (Willgoose (2005), p. 79), the hydraulic flow can be ignored
when focused on the sediment flow (advection) with the equation
/3 /3m1(1-m4)
with /3 = 1 3/3-;'1 and m = m1 m3 (1 -m4)*. This is the equation used for advection.
Diffusive sediment transport is modeled as a diffusivity constant multiplied by the slope, yielding
a full sediment transport equation of qs = f]Amsn1 + DS. The expression DS can be derived
from the SIBERIA user manual (Willgoose (2005), p. 8, eq. 2.3.2) expression
{
DSDzn Sthreshold -------Dzt
q _ Sthreshold -S s -
o
DSDznSthreshold ------>Dzt
Sthreshold -S
DSDznS ' threshold ------:'.SDzc
Sthreshold -S
and assumptions that St » S, Dzn = 1, and Dzt is small. The assumptions of Dzn = 1, and Dzt =
0 are characteristic of simple diffusion, and the assumption that St » S corresponds to the
assumption that slopes are not incredibly steep, which is reasonable given that we are considering
2-3% grades.
To derive distributions for the SIBERIA parameters bl, ml, nl, and dz from the RHEM output,
multiple fits are used from the SIBERIA output. To find the advective term, the equation qs =
f]Amsn1 is fit to all of the data, to derive estimates for the parameters bl (/3), ml (m), and nl
(n1). The equation qs = f]Amsn1 + DS is fit to the RHEM output for each of the 15,000 using
only the runs with A :'.S 30 to find the diffusion term dZ (D). Both fits use the Levenberg-
Marquardt algorithm to estimate the parameters of interest. Only the fits with A :'.S 30 were used
to find dZ (D) because diffusivity is dominant at shorter areas. The values of the parameters bl
(/3), ml (m), and nl (n1) were taken from the fit with all of the data beca use advection is
dominant on larger areas.
This results in I 000 fits for each of the 15 climate scenarios, with the dataset used in each fit
being the 39 combinations of length and slopes combined with the RHEM output in each
situation. The climate scenarios are combined to increase the variance of the sampling
distributions, and the empirical distributions from the 15,000 fits of these parameters are then
used as input di stributions in SIBERIA.
Running SIBERIA
Once the parameters are obtained, SIBERIA is run I 0,000 times. The empirical distributions (i.e.,
raw fitted values) of the 15,000 fits to data from RHEM are then put into SIBERIA as the
distributions for the parameters ml (m), nl (n1) and dZ (D). The parameter bl (/31) is drawn in
each year of the SIBERIA runs from the empirical distribution of the 15,000 RHEM fits for bl.
Output is obtained for I, 10, 26, 50, I 00, 200, 300, 400, and 500 years. The simulation timesteps
were I-year for the first IO years, 2-years until the 50th simulation year, and I 0-years after that.
The numerical timesteps were 0.025 yrs (40 timesteps per year).
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
(80 1) 649-2000 • Fax: [phone] • www.energysolutions.com
Page 51 of77
:;;::::::=-
ENERGYSOLUTIONS
Notes
Mr. Doug Hansen
CD-2025-037
February 21 , 2025
• The SIBERIA User Manual (Willgoose 2005) contains a different derivation. Their equations
6.2. 1 and 6.2.2 are
Using these equations, they obtain that ( eq. 6.1.3)
which is then noted to have the fom1 ( eq. 6.1.4)
/Ji Qm15n1 > Qt
p1 Qm15n1 ::; Qt'
(/Ji/J"(;3)Am1m35n1 > Qt
(/Ji/J"(;3 )Am1m35n1 ::; Q/
f]Amsn1 > Qt
f]Amsn1 ::; Qt'
In these equations, the term Q5 is never defined, and in equations 6.1.3 and 6.1.4, Qt is replaced
by At for unknown reasons. Also, it is unclear why the SIBERIA User Manual conflates Qc
(defined as the discharge in a channel) and q (defined as the discharge per unit width) in these
equations, only explaining that "it should be noted that these are the mean annual equations." Due
to this lack of clarity, an alternative derivation is presented.
Table 2-1 compiles parameters used in Barnhart et al. (2020) and in SIBERIA.
299 South Main Street, Suite 1700 • Sall Lake City, Utah 8411 I
[phone] • Fax: (80 I) 880-2 879 • www.energysolutions.com
Page 52 of77
~
ENERGY SOLUTIONS
Mr. Doug Hansen
CD-2025-037
February 21 , 2025
Table 2-1. Parameters with variable names from Barnhart et al. (2020) and from the SIBERIA
User Manual. Units are expressed in terms of length (L), time (T), and mass (M).
Parameter
Horizontal distances
Elevation
Time
Diffusivity
Channel vs. hillslope
Channel/flow width
Density of water
Bulk density of soil
Gravitational acceleration
Nondimensional buoyant
density
Grain size
Friction factor; Manning's
roughness formula
Friction factor; Darcy-
Weisbach roughness
formula
Nondimensional boundary
Shields stress
Boundary shear stress
Dimensionless critical
shear stress
Critical shear stress
Erosion efficiency
constants
Water velocity
Critical water velocity
Erosion exponent
Dimensionless sediment
discharge per unit width
Volumetric discharge
Units Barnhart Range
L x,y
L z
T t
L2 /T
Unitless
L w
M/L3 (kg/m3) Pw 1000
M/L3 (kg/m3) Ps [1200,3300]
L/T2 (m/s2) g 9.81
M/L3 Rb [0.2,2.3]
L Ds
M/(Lll/5 . T7f5) kr
M/(L7f3 . r4/3) kr
Dimensionless ,·
M/(L • T2) T
Dimensionless ,· C [0.03 ,0.08]
L2 /T 'c 8[0 .6,240)
ke?
L/T u
L/T u· C
Unitless a {1.5,3}
Dimensionless q*
L3 /T Q
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
[phone] • Fax: [phone] • www.energysolutions.com
SIBERIA Range
x,y
z
t
y [0,1]
w
Pb
p [2.5, 3.5]
QC
Page 53 of77
Parameter
Volumetric discharge per
unit width
Sediment transport
Magnitude of sediment
transport
Volumetric sediment flux
Volumetric sediment flux
per unit width
Slope of channel
Drainage Area
Area-discharge coefficient
Area-discharge exponent
Discharge-width coefficient
Discharge-width exponent
Manning's n
~
ENERGY SOLUTIONS
Units Barnhart
L2 /T q
M/T
M/T
L3 /T Qs
L2/T qs
Dimensionless s
L2 A
L3-2qc /T kq
Unitless Cq
L1 -3cwr cw kw
Unitless Cw
T /L1/3 nm
Range
(0, 1]
(0.001 ,30]
(0.70,1.0]
(5.5,9.1]
(0.37,0.62]
(0.01 ,0.1]
Discharge-stress coefficient ML2a -1 /T kt (23,91600]
Discharge-stress exponent Unitless a (0.60,0.67]
for discharge
Discharge-stress exponent Unitless {3 (0.66,0.70]
for slope
Discharge-sediment L 3-2m1 /T Kt b[1 E-04,
coefficient 0.01]
(transport-limited)
Discharge-sediment Unitless mt (0.052 ,0.24]
exponent for discharge
(transport-limited)
Discharge-sediment Unitless nt (1,2]
exponent for slope
(transport-limited)
Transport threshold Wt,c
Tectonic uplift
Channel initiation function
a detachment-limited transport
b range of Kt assumes Darcy-Weisbach friction formula
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
(801 ) 649-2000 • Fax: (801 ) 880-2879 • www.energysolutions.com
Mr. Doug Hansen
CD-2025-037
February 21 , 2025
SIBERIA Range
q
qsx, qsy
qs
s (0,1]
A
{33
m 3
/34
m 4
/31
m1
n 1
qst, Qt, At
u
a
Page 54 of77
~
ENERGYSOLUT/ONS
Mr. Doug Hansen
CD-2025-037
February 21, 2025
An important thing to note in Table 2-1 is that the sediment discharge is expressed in terms of
volume in Barnhart et al. (2020), but sediment discharge is expressed in terms of mass in
SIBERIA.
Table 2-2 lays out equations used to model erosion in the respective codes.
299 South Main Street, Suite 1700 • Salt Lake City, Utah 8411 I
[phone] • Fax: [phone] • www.energysolutions.com
Page 55 of77
~
ENERGY SOLUTIONS
Mr. Doug Hansen
CD-2025-037
February 21 , 2025
Table 2-2. Equations used to model erosion in Barnh art et al. (2020) and the SIBERIA User
Manu al.
Eq #
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
Barnhart SIBERIA
E= ke1(T-Tc)a
E = k ez (TU -TcUDa
qs =
• qs q = D5✓ RdgD5
T r * =
PwRdgDs
q* ~ (T* -T~)3/2
q* ~ (r*)3
T~kr q asfJ
cJz 1 -a =U+-'il·qs +D'il2z t Pb
cJY ( a) at = f dt, Y, at
a= {35qms5ns
Q= k ACq q Q -/3 Am3 C -3
W= k wQCw W = /34Q";4
QC q= -w
qs = {/31 qmi s;1 -qst f31qm15n1 > qst
/3 qm15n1 < q 1 -st
Qs K Amt5nt -W -{{3Am5n1 -A t {3Amsn1 > At -t t,c --qs -0 {3Amsn1 ~ At w w
{
DS0 zn Sthreshold DSDzn Sthreshold
-Dzt
qs = Sthreshold -S S threshold -S
DS0zn Sthreshold 0
Sthreshold -S
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
[phone] • Fax: [phone] • www.energysolutions.com
> Dzt
~ Dzt
Page 56 of77
~
ENERGYSOLUTIONS
Mr. Doug Hansen
CD-2025-037
February 21, 2025
In Table 2-2, equations 1-7 set up the scenario in Barnhart et al. (2020) and equations 8-10 help
set up the situation in the SIBERIA user manual. Equation 11 relates discharge to area, and
equation 12 is the geometric relationship between width and discharge. Equation 13 defines the
discharge per unit width. The equivalent expression does also work in the Barnhart paper, but it is
not explicitly defined.
Equation 14 is the interim step to getting to the final form of equation 15 in SIBERIA. Equation
15 is the key equation that is used for fluvial flow in SIBERIA. Equation 16 represents the
diffusive sediment flow in SIBERIA and can be simplified to q5 = DS under simple diffusion.
In practice, total diffusion is measured in SIBERIA, which involves adding the left-hand sides of
equations 15 and 16, while setting the thresholds equal to Oto get the equation q5 = f]Amsn1 +
DS.
It is also worth noting that the SIBERIA user manual has multiple forms of these equations that
seem to conflict. For equation 11 , the form from SIBERIA that we chose to include was equation
2.1.7 in the SIBERIA user manual. However, equation 4.3.1 is q = f]3Am3 even though q is the
discharge (of water) per uni t width while Qc is the total discharge of water (for the entire width of
the channel); equation 6.2.1 is Q = f]3A1113 , where Q is never explicitly defined in the SIBERIA
user manual. Similarly, equation 14 has the alternate form
Here, Q is never defined, and it is unclear why the threshold term changes from q5t to Qt to At in
the different forms of equation 14 and equation 15.
Table 2-3 provides summary statistics for the parameters evaluated.
Table 2-3. Summary statistics for the parameters. QsHold(Mg) values are from 10,000 random
samples from the distribution in Barnhart et al. (2020). The other four parameters (bl, ml, nl, and
dZ) are informed by the 15,000 iterations of the Cligen-RHEM parameter estimation approach.
count
b1 15000
m1 15000
n1 15000
dZ 15000
QsHold(Mg) 10000
mean std min 25.00% 50.00%
6.59E-04 1.59E-04 2.73E-04 5.45E-04 6.47E-04
1.18E+00 1.52E-02 1.12E+00 1.17E+00 1.18E+00
1.31 E+00 4.67E-02 1.14E+00 1.28E+00 1.31 E+00
2.37E-04 2.46E-05 3.41 E-05 2.19E-04 2.35E-04
1.41E-05 2.64E-05 9.29E-08 1.66E-06 4.70E-06
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
[phone] • Fax : [phone] • www.energysolutions.com
75.00% max
7.61 E-04 1.44E-03
1.1 9E+00 1.23E+00
1.34E+00 1.46E+00
2.53E-04 3.45E-04
1.47E-05 4.63E-04
Page 57 of77
~
ENERGY SOLUTIONS
Appendix 3: Supporting Tables and Figures
Mr. Doug Hansen
CD-2025-037
February 21, 2025
Table 3-1 and Table 3-2 present centennial summary statistics for each transect under unarmored
landform and engineered cover conditions, respectively. The median 95th percentile (95th)
highlights the extremes by capturing the median of the 95th percentile incision depth across
realizations, while the mean of means (means) provides a measure of the central tendency by
averaging mean depths across realizations. The median of medians (Med) offers a robust
indicator of typical erosion depth by focusing on the central tendency of median values, and the
upper and lower confidence limits (UCL and LCL) for both the mean and median depths quantify
the plausible range of erosion outcomes. Together, these metrics enable a direct comparison of
how the engineered cover mitigates erosion relative to the unarmored landform over the
centennial timescale.
Table 3-1. Summary statistics from projected erosion of gully depths for unarmored landform
cover condition.
Time Transect
100 10m (west)
100 30m (west)
100 50m (west)
100 75m (west)
100 1 00m (west)
100 120m (west)
100 140m (west)
100 east
100 north
100 south
200 10m (west)
200 30m (west)
200 50m (west)
200 75m (west)
200 100m (west)
200 120m (west)
200 140m (west)
200 east
200 north
200 south
95th means UCL mean LCL mean Med
-0.15 -0.06 -0.09 -0.05 -0.06
-0.34 -0.06 -0.1 -0.04 -0.05
-0.55 -0.08 -0.13 -0.06 -0.07
-0.83 -0.13 -0.17 -0.1 -0 .11
-0.85 -0.1 -0.15 -0.07 -0.05
-0.85 -0.14 -0.32 -0.07 -0 .07
-28.9 -7.2 -9.87 -5.62 -0 .75
-25.0 -5.78 -8.0 -4.28 -1 .03
-14.2 -2.65 -3.8 -1.94 -0.49
-16.2 -3.23 -4.61 -2.35 -0 .54
-0.28 -0.13 -0.18 -0 .1 -0.12
-0.67 -0.15 -0.22 -0 .1 -0 .12
-0.95 -0.19 -0.27 -0.14 -0.15
-1.11 -0.24 -0.33 -0.18 -0.2
-1.24 -0.22 -0.34 -0.16 -0 .12
-3.72 -0.82 -1.64 -0.44 -0 .26
-59.4 -14.9 -19.5 -12.4 -2 .64
-51 .8 -12.6 -16.2 -10.2 -2 .31
-34.4 -7.09 -9.61 -5.68 -1.6
-41.2 -8.07 -10.5 -6.57 -1.71
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
(80 1) 649-2000 • Fax: [phone] • www .energysolutioos.com
UCL med LCL med
-0.09 -0.05
-0.09 -0.04
-0.1 -0.06
-0.15 -0.09
-0.07 -0.04
-0.11 -0.06
-1 .05 -0.62
-1.19 -0.86
-0.76 -0.37
-0.73 -0.42
-0.19 -0.1
-0.17 -0.09
-0.2 -0.12
-0.27 -0.15
-0.19 -0.09
-0.35 -0.17
-3.58 -2 .25
-3.01 -1 .75
-2 .03 -1 .31
-2.41 -1.45
Page 58 of77
Time Transect
300 1 Om (west)
300 30m (west)
300 50m (west)
300 75m (west)
300 100m (west)
300 120m (west)
300 140m (west)
300 east
300 north
300 south
400 10m (west)
400 30m (west)
400 50m (west)
400 75m (west)
400 100m (west)
400 120m (west)
400 140m (west)
400 east
400 north
400 south
500 10m (west)
500 30m (west)
500 50m (west)
500 75m (west)
500 100m (west)
500 120m (west)
500 140m (west)
500 east
500 north
500 south
~
ENERGY SOLUTIONS
95th means UCL mean LCL mean Med
-0.41 -0.19 -0.28 -0.15 -0.2
-1 .03 -0.23 -0.33 -0.17 -0 .18
-1 .32 -0.3 -0.41 -0.21 -0 .22
-1.32 -0.35 -0.49 -0.27 -0 .29
-1 .64 -0.41 -0 .72 -0.25 -0.21
-10.3 -2 .04 -3 .98 -1.14 -0.43
-82 .8 -21 .3 -27.3 -18.1 -5.29
-74.5 -18.2 -22.9 -15.1 -4.24
-49.0 -11 .3 -14.8 -9.3 -2 .97
-59.4 -12.3 -15.4 -10.4 -3.76
-0.55 -0.26 -0 .38 -0.2 -0.24
-1 .17 -0.31 -0.45 -0.23 -0.27
-1.52 -0.39 -0 .57 -0.29 -0.31
-1 .6 -0.47 -0.66 -0.35 -0.39
-2.63 -0.69 -1.15 -0.4 -0.32
-18.9 -3.75 -6.69 -2.23 -0.74
-100.7 -27.0 -34.0 -23.2 -8.62
-94.6 -23.2 -28.5 -19.5 -7.16
-59.8 -15.2 -19.2 -12.7 -5.04
-73.6 -16.1 -19.8 -13.7 -5.77
-0.71 -0.33 -0.48 -0 .25 -0.31
-1.38 -0.4 -0 .58 -0.29 -0.32
-1 .68 -0.51 -0 .74 -0.36 -0.38
-1 .84 -0.6 -0.87 -0.44 -0.48
-3.91 -1.03 -1.72 -0 .62 -0.43
-27.1 -5.76 -9.63 -3.6 -1 .26
-113.4 -32.1 -40.1 -27.8 -12.7
-110.1 -27.6 -33.8 -23.5 -11 .0
-70.3 -18.7 -23.3 -15.8 -7.1
-84.7 -19.6 -24.0 -16.8 -7.72
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
[phone] • Fax: [phone] • www .energysolutions.com
Mr. Doug Hansen
CD-2025-037
February 21, 2025
UCL med LCL med
-0.28 -0.15
-0.29 -0.15
-0.32 -0.18
-0.4 -0.23
-0.33 -0.15
-0.67 -0.32
-7.06 -4.51
-5.88 -3.54
-4.07 -2.42
-4.77 -3.02
-0.36 -0.2
-0.37 -0.2
-0.45 -0.23
-0.56 -0.31
-0.51 -0.23
-1 .3 -0.52
-12 .0 -7.16
-9.79 -5.81
-6.5 -4.2
-7.45 -4.8
-0.46 -0.24
-0.49 -0.27
-0.57 -0.3
-0.73 -0.38
-0.73 -0.34
-2.4 -0.79
-16.5 -10.7
-14.6 -8.7
-9.51 -5.93
-10.4 -6.28
Page 59 of77
~
ENERGY SOLUTIONS
Mr. Doug Hansen
CD-2025-037
February 21, 2025
Table 3-2. Summary statistics from projected erosion of gully depths for engineered cover condition
at the FCF.
Time Transect
100 10m (west)
100 30m (west)
100 50m (west)
100 75m (west)
100 100m (west)
100 120m (west)
100 140m (west)
100 east
100 north
100 south
200 10m (west)
200 30m (west)
200 50m (west)
200 75m (west)
200 100m (west)
200 120m (west)
200 140m (west)
200 east
200 north
200 south
300 10m (west)
300 30m (west)
300 50m (west)
300 75m (west)
300 100m (west)
300 120m (west)
300 140m (west)
300 east
95th means UCL mean LCL mean Med
-0.15 -0.07 -0.11 -0.03 -0.06
-0.34 -0.07 -0.13 -0.03 -0 .05
-0.56 -0.09 -0.17 -0.04 -0 .07
-0.85 -0.13 -0.22 -0.07 -0 .11
-0.86 -0.11 -0.2 -0.04 -0 .05
-0.83 -0.14 -0.38 -0.05 -0 .07
-7.41 -2.0 -2.86 -1 .08 -0.43
-6.61 -1.82 -2.81 -0.84 -0.88
-3.71 -0.96 -1.72 -0.32 -0.27
-4.42 -1.08 -1.91 -0.4 -0 .34
-0.29 -0.13 -0.23 -0.07 -0.12
-0.7 -0.15 -0.28 -0.07 -0.12
-0.93 -0.2 -0.35 -0.09 -0 .15
-1.11 -0.25 -0.42 -0.14 -0.2
-1 .23 -0.24 -0.48 -0.11 -0 .11
-2.37 -0.48 -0.95 -0.12 -0.21
-11.1 -3.06 -3.87 -2.17 -0.73
-10.3 -3.08 -4.03 -1 .99 -1 .26
-6 .77 -1 .99 -2 .96 -1.09 -0.9
-8 .02 -2.21 -3.23 -1.19 -0.84
-0.42 -0.2 -0 .36 -0.1 -0 .2
-1 .03 -0.24 -0.43 -0.11 -0.18
-1 .32 -0.3 -0 .55 -0 .15 -0.22
-1.33 -0.36 -0 .64 -0.19 -0.29
-1.6 -0.4 -0.86 -0.17 -0.21
-3.6 -0.81 -1 .35 -0.3 -0.39
-12.5 -3.63 -4.45 -2 .8 -1 .14
-12 .0 -3.8 -4.71 -2.76 -1.64
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
[phone] • Fax: [phone] • www.eoergysolutions.com
UCL LCL
med med
-0.11 -0 .04
-0.11 -0.02
-0.12 -0 .05
-0.18 -0 .06
-0.1 -0 .04
-0.16 -0 .05
-0.6 -0 .26
-1 .11 -0 .37
-0.68 -0 .13
-0.68 -0.11
-0.23 -0.06
-0.24 -0.06
-0.27 -0.09
-0.34 -0.11
-0.24 -0.06
-0.4 -0 .09
-1 .14 -0 .54
-1 .64 -0.95
-1 .51 -0.37
-1 .26 -0.49
-0.34 -0.1
-0.34 -0.1
-0.42 -0.12
-0.54 -0.16
-0.44 -0.1
-0 .62 -0.16
-1 .75 -0.75
-2.17 -1.25
Page 60 of77
Time
300
300
400
400
400
400
400
400
400
400
400
400
500
500
500
500
500
500
500
500
500
500
~
ENERGY SOLUTIONS
Transect 95th means UCL mean LCL mean
north -8.54 -2.68 -3.72 -1 .69
south -10.2 -2 .98 -4.05 -1 .87
10m (west) -0 .55 -0 .27 -0.49 -0.14
30m (west) -1 .17 -0 .33 -0.59 -0.16
50m (west) -1 .52 -0.41 -0.76 -0 .2
75m (west) -1 .58 -0.49 -0.89 -0.26
100m (west) -2 .0 -0 .58 -1 .2 -0.23
120m (west) -4.09 -1.08 -1 .7 -0.51
140m (west) -13.2 -4.03 -4.88 -3 .2
east -12.9 -4.29 -5.19 -3 .3
north -9 .3 -3.2 -4.29 -2.17
south -11 .3 -3.54 -4.64 -2.38
10m (west) -0.72 -0.35 -0.62 -0.17
30m (west) -1 .36 -0.41 -0.76 -0.2
50m (west) -1 .7 -0.53 -0.99 -0.25
75m (west) -1 .84 -0.62 -1 .16 -0.31
100m (west) -2.44 -0.78 -1.49 -0.3
120m (west) -4.6 -1 .32 -2 .05 -0.72
140m (west) -13.6 -4.36 -5.24 -3.51
east -13.5 -4.66 -5.58 -3.71
north -9.77 -3.62 -4.75 -2.58
south -12.1 -3.99 -5.12 -2.79
Med
-1.44
-1 .25
-0.26
-0.27
-0.31
-0.39
-0.31
-0.54
-1 .59
-1 .98
-1 .94
-1.66
-0 .31
-0.32
-0 .39
-0.48
-0.43
-0 .67
-1 .99
-2.32
-2 .36
-2 .1
Mr. Doug Hansen
CD-2025-037
February 21 , 2025
UCL LCL
med med
-2 .27 -0.81
-1 .99 -0.82
-0.46 -0.13
-0.48 -0.13
-0.6 -0.17
-0.75 -0.21
-0.66 -0.13
-0.83 -0.28
-2.28 -1 .01
-2 .69 -1 .5
-2 .73 -1 .16
-2 .64 -1 .1
-0.59 -0 .17
-0.62 -0.17
-0.79 -0.21
-1 .0 -0.29
-0.89 -0.2
-1 .06 -0.39
-2 .8 -1 .29
-3.16 -1 .71
-3.16 -1.47
-3.17 -1 .37
Table 3-3, adapted from Rogers (2023b), presents erosion markers for the top slope of the FCF at
10,000 years as affected by varying QsHold on a 1.5-meter resolution model DEM as discussed
in Section 2.1.2.
Table 3-4 and Table 3-5 present effects of varying dZ on erosion markers across the top slope of
the FCF at 10,000 years, adapted from Table 4 and Table 5 in Rogers (2023a). Table 3-3 presents
results on a 0.5-meter resolution model DEM while the DEM resolution reflected in Table 3-5
was 1.5 meters. Quantities in these tables are discussed in Section 2.1.3.
299 South Maio Street, Suite 1700 • Salt Lake City, Utah 841 11
(801 ) 649-2000 • Fax: [phone] • www.energysolutions.com
Page 61 of77
~
ENERGY SOLUTIONS
Mr. Doug Hansen
CD-2025-037
February 21 , 2025
Table 3-3. Effects of large values of QsHold on markers of erosion across the top slope of the FCF at 10,000 years (Table 2 in Rogers
(2023b)).
QsHold Scenario calibrate (Mg)
Draw 725, 2023 RHEM 0
Draw 725, 2023 RHEM 0.00012
Draw 725, 2023 RHEM 0.0012
Draw 725, 2023 RHEM 0.012
Draw 725, 2023 RHEM 0.12
Draw 725, 2021 LANL_analog 0
Draw 725, 2021 LANL_analog 0.00012
Draw 725, 2021 LANL_analog 0.0012
Draw 725, 2021 LANL_analog 0.012
Draw 725, 2021 LANL_analog 0.12
10 11 12 13
0 0 0.995 0.005
0 0 0.994 0.006
0 0 0.987 0.013
0 0.007 0.882 0.111
0.103 0.217 0.606 0.061
0.098 0.008 0.838 0.043
0.095 0.016 0.831 0.043
0.099 0.059 0.778 0.047
0.136 0.13 0.662 0.05
0.295 0.366 0.299 0.017
299 South Mam Street, Suite 1700 • Salt Lake City, Utah 84111
(80 I) 649-2000 • Fax: (80 I) 880-2879 • www.energysolutions.com
15 Maxlncision annual SY 14 through (m) (Mg ha·1) 115
0 0 -0.386 -0.24
0 0 -0.387 -0.239
0 0 -0.374 -0.225
0.001 0 -0.501 -0.192
0.012 0 -0.558 -0.104
0.013 0 -0.608 -0.171
0.015 0 -0.612 -0.169
0.017 0 -0.605 -0.163
0.021 0 -0.613 -0.151
0.023 0 -0.673 -0.067
Page 62 of77
~
ENERGYSOLUTIONS
Mr. Doug Hansen
CD-2025-037
February 21, 2025
Table 3-4. Effect of varying dZ on erosion markers across the top slope of the FCF at 10,000 years for three cover conditions realizations
using a model DEM with 0.5-m grid spacing. From Table 4 of Rogers (2023a).
Cover RSYith GridXy Dz (m2 y·1)
realization (m)
738 0 0.5 1E-03
0.5 0
1.5 0
0.5 1E-05
316 50 0.5 1E-03
0.5 0
1.5 0
0.5 1E-05
725 100 0.5 1E-03
0.5 0
1.5 0
0.5 5
11 12 13 14
0 0.986608 0.012698 0.000693
7.83E-06 0.999419 0.000573 0
0 0.999366 0.000634 0
0 0.999292 0.000708 0
0 0.966505 0.029777 0.003673
0 0.995838 0.003959 0.000204
0 0.998345 0.001655 0
0 0.993711 0.006264 2.54E-05
0 0.925585 0.062737 0.011353
0 0.974912 0.023519 0.001569
0 0.994717 0.005283 0
0 0.987199 0.012801 0
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
[phone] • Fax: [phone] • www.energysolutions.com
15 16 Max SY
incision (Mg ha·1y"1)
(m)
1.96E-06 0 -0.61536 -0.19449
0 0 -0.38611 -0.16662
0 0 -0.34778 -0.17147
0 0 -0.34045 -0.16332
4.50E-05 0 -0.66797 -0.23463
0 0 -0.51855 -0.20299
0 0 -0.38623 -0.21027
0 0 -0.4873 -0.20258
0.000323 1.96E-06 -0.776 -0.27736
0 0 -0.56616 -0.23979
0 0 -0.38635 -0.23981
0 0 -0.41577 -0.2351
Page 63 of77
~
ENERGY SOLUTIONS
Mr. Doug Hansen
CD-2025-037
February 21 , 2025
Table 3-5. Effect of varying dZ on erosion markers across the top slope of the FCF at 10,000 years. Expanded from Rogers (2023a). 57459
DEM cells on the top slope of the FCF at 1.5-m grid spacing.
sen rough_st dZ (m2 y· t (y) max min Sedimen 10 11 12 13 14
dev(m) ') incision incision/ t Yield
(m) max (Mg ha·'
depositio y·')
n(m)
d725 0.002 0.1 10000 -5.00354 -0.43836 -2.14868 0 0 0 0.001566 0.06375
d725 0.002 0.01 10000 -2.262085 -0.12805 -0.42846 0 0 0.631929 0.24226 0.079396
d725 0.002 0.001 10000 -0.666992 0.200439 -0.13715 0.110635 0.012148 0.873231 0.00322 0.000731
d725 2.00E-03 1.00E-04 10000 -0.227173 0.289551 -0.10611 0.169025 0.013679 0.817296 0 0
d725 0.00E+00 1.00E-04 10000 -0.231445 0.285522 -0.10418 0.171253 0.014184 0.814563 0 0
d725 0.00E+00 1.00E-05 10000 -0.227905 0.292114 -0.11042 0.159592 0.013001 0.827407 0 0
d725 0.00E+00 1.00E-06 10000 -0.22937 0.291138 -0.11184 0.161106 0.011887 0.827007 0 0
d725 0.00E+00 1.00E-07 10000 -0.290894 0.296997 -0.10954 0.1653 0.012252 0.822447 0 0
1The top of cover layer 5 is the top of the Frost Protection Layer
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84 111
[phone] • Fax: [phone] • www.energysolutions.com
115 16 17 18 19
0.08587 0.094815 0.110601 0.105397 0.073287
0.02165 0.008023 0.005273 0.003777 0.002576
3.48E-05 0 0 0 0
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
Page 64 of77
110
0.051567
0.001932
0
0
0
0
0
0
:::::=-
ENERGYSOLUTIONS
Mr. Doug Hansen
CD-2025-037
February 2 I, 2025
A complete set of plots of net elevation change at simulation year 500 for scenario 1800 is shown
in Figure 3-1 through Figure 3-19.
0.100
-·-.. T.,-
0.075
400 NSSL
TSN
-'
350 VJ .. ~ VJ 0.050 w ~ ~ • -/ .
300
0.025
250
. 0.000 1 ..
200
w (/)
f--0.025
150 3:: 3:: 3:: 3:: 3:: 3:: (/) (/) (/) (/) (/) (/) f-f-f-f-f-f-00 ll} 000 NO .... ll} C') ~
100 -0.050
TSS
50
u -0.075
0
0 50 100 150 200 250 300
-0.100
Figure 3-1. Locations of transects analyzed and plotted in this report. Transect locations for gully width and depth
analysis are overlaid onto delta-Z at 500 years for realization 1800. Arrows indicate direction of plotting. The color scale
of the elevation change raster is normalized to interval 1-0.1, 0.11 meters but erosion and deposition may exceed those
boundaries.
Elevation change along the transects are plotted in Figure 3-2 through Figure 3-19. The top frame
in each figure is sampled from the engineered cover and scaled vertically to the range of elevation
change through 500 years. The middle frame is the same except the vertical scale is based on
elevation change through 500 years for the unarmored landform. The lower frame is sampled
from the unarmored landform with the appropriate vertical scale.
For all figures colors indicate the following: Blue: year 5; Orange: year 50; Green: year I 00; Red:
year 300; Purple: year 500.
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
[phone] • Fax: [phone] • www.energysolutions.com
Page 65 of77
~
ENERGYSOLUTIONS
Mr. Doug Hansen
CD-2025-036
February 21, 2025
!:~~1 ~1 -o.ooJ 'VV • ~,, lv,.r'1~ f V , .
'-----,-------+-, --~. -~.--~. -
0 50 100 150 200
! ;::1 I
0 50 100 150 200
1 ;~:I I
0 50 100 150 200
distance (m)
Figure 3-2. Transect JOTSW.
! ~;~I I
0 50 100 150 200
!;~1 ~~1
0 50 100 150 200
! ;~1 I
• 0 so 100 150 20C
distance (m)
Figure 3-3. Transect 30TSW.
299 South Main Street, Suite 1700 • Salt Lake City, Utah 841 11
[phone] • Fax: [phone] • www.energysolutions.com
Page 66 of77
~
ENERGY SOLUTIONS
Mr. Doug Hansen
CD-2025-036
February 21, 2025
! ~~~l I
0 50 100 150 200
! :::! I
0 50 -= =
!:::1 ~~ I
0 50 100 150 200
distance (ml
Figure 3-4. Transect 50 TSW.
! :~~1 I
0 50 100 150 200
!~~1 I
0 SO 100 150 200
!~] I
a so 100 iso 200
distance (m)
Figure 3-5. Transect 75 TSW.
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
(80 1) 649-2000 • Fax: [phone] • www.energysolutions.com
Page 67 of77
~
ENERGY SOLUTIONS
Mr. Doug Hansen
CD-2025-036
February 21 , 2025
!J~I ~~~~w~ !
b ~ -~ ~
!::[ ~~v~ •
0 ~ -= =
!:::f ~rn·vvrm I
0 so 100 150 200
distance (ml
Figure 3-6. Transect 100 TSW.
! ::1 !
b so 100 150 200 ! ill '-1'• .. ,.. ... , ,...., pi-,., ... , .. ,,,..,.,,..,,.,..,,,,...,....,,,,..,,.,,., ,.....,,,,,,.. ... ,1.,,,,,.,,... •
LU -0.4
-0.S L---,,------r-, --~. -~.------r-, ~
O so 100 150 200
!l)I ~lwrrr .. rrnn, ·,rrvrn~v I
O ~ 100 150 200
distance (m)
Figure 3-7. Transect 120TSW.
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
[phone] • Fax: [phone] • www.energysolutions.com
Page 68 of77
~
ENERGYSOLUT/ONS
Mr. Doug Hansen
CD-2025-036
February 21, 2025
;~j I
0 ~ = = -
0.01 ''' ...... .,, .. , • w ,. ♦-¥-.. W•t • .. ,.4't ... , ........... , ... " I
!:: .
~ -= -
:~~1 :r1·rnmvm,1r 1,miry I
0 ~ -= -distance (m)
Figure 3-8. Transect 140TSW.
~11 'wwr,rir1~1r1m111~~~'dl'F I
0 100 200 300 400 500
distance (m)
Figure 3-9. Transect WSSU.
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
(80 1) 649-2000 • Fax: [phone] • www.energysolutions.com
Page 69 of77
~
E ERGYSOLUTIONS
Mr. Doug Hansen
CD-2025-036
February 21, 2025
0.15
E 0.10
~ a.as
V,
w 0.00 i
100 200 300 400 500
I
-. ~--11 Ma -.S a: : ~ I ::::: , ~ -0.75
-1.00
'---,----~--~. ---~--~. ---~~
•• Q t
0.00
I -a.2s
,':! -a.so
~ -0.75
-1.00
0 100 200 300 400
100 200 300
distance (ml
Figure 3-10. Transect WSSL.
400
500
500
, ~~~I I
0 50 100 150 200
a.01, .. ,.,.,v •""' ,,.,.., V''"'"',,.,...,..,,,,-.,11..,.,,,,, ..... .,,,, ... ,,,, i
E --0s
L, , , ,
0 50 100 150 200
!::1 ''l 1irrnwrvw11nrr~ I
0 50 100 150 200
distance (m)
Figure 3-1 l. Transect E.
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
(80 1) 649-2000 • Fax: [phone] • www.energysolutions.com
Page 70 of77
0.3
f 0.2
~ O.l
"' w
a.a
a.a
f ~ -05
~ -1.0
-1.5
0.10
f ~ 0.05
~ 0.00
-0.05
a.a
f -0.2
~ -0.4
~ -0.6
-0.8
a.a
K -0.2
~ -0.4
~ -0.6
-0.8
~
ENERGYSOLUTIONS
Mr. Doug Hansen
CD-2025-036
February 21, 2025
,.J~JM~~1,,
0
a
t
J
a
a
I
a
100
100
100
200 300
200 300
200 300
distance (m)
Figure 3-12. Transect ESSU
400
400
400
4•~w~~-4-,,
100 200 300 400
100 200 300 400
100 200 300 400
distance (m)
Figure 3-13. Transect ESSL.
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
[phone] • Fax: [phone] • www.energysolutions.com
500
500
500
\
500
500
500
Page 71 of77
~
ENERGYSOLUTIONS
Mr. Doug Hansen
CD-2025-036
February 21 , 2025
0.000 r-----... ~~~~--~~. ~~~-.~ . ..--~-111 1
_ -0.025
E li' -0.050
vi -0.075
LlJ -0.100
-0.125 ~----------~----~-----~---~
50 100 150 200
o.o 1 •1 ;;.;4.;;,..,~.-ioii+ii-i+iii•iiiii1 ei-i,ii¥iiiiif~•ii♦♦ii<i.,.iiiiJii♦<i;iiiii*iiiii•iit"'Q"•"•t"Y'..,'"t"''ftiiiii1 f''"'l¥Wini,•111-ii9ii,-.,,iii"Mi;iii••t11•••--1
-0.2
:[ -0 4
~ -0.6
Ill W -0.8
-1.0
50 100 150 200
1 ~l ~!t'·rr~,,,,,iiiwiiiiiiii,1rtiiiiiiiiiii~~~iiiiiiiiijii,ij~m~~~~~~, ...
0.10 E li' 0.05
Ill LlJ
0.00 .I I II
0
0.00
E -0.25
li' -0.50
Ill -0.75 LlJ
-1.00
0
0.00
I -0.25
.!:l -0.50
,:.; -0.75 :::,
-1.00
0
50 100
distance (m)
Figure 3-14. Transect N.
so 100 150
.. JI ft ..._
so
so
100 150
100 150
distance (m)
Figure 3-15. Transect NSSU.
150
200
200
200
299 Soulh Main Street, Suite 1700 • Salt Lake City, Utah 84111
[phone] • Fax: (801 ) 880-2879 • www .energysolutions.com
200
250
• I •
250
250
Page 72 of77
~
ENERGY SOLUTIONS
1 ~:1 I ::i .... .
'---,---~ -__ _,
50 100
UO 200 2SO
ISO
d1!>1ance Im)
200 '"
Figure 3-16. Transect NSSL.
)00
~,,~~rw "'F I
50 100 >00
E -0!)
001 • ►++•• ···••t4ff0\if'◄*f' fO\,i :e,;o:ft+W ♦♦--. * I
LIO.____.__~---.------...--
100 ,00 ! :: I ··•vrr~f wmrymrrvr, fYP ffi" I
O 50 100 ,oo
d¥$lance (m}
Figure 3-17. Transect S.
299 South Main Street, Suite 1700 • Salt Lake City, Utah 8411 I
[phone] • Fax: [phone] • www.energysolutions.com
Mr. Doug Hansen
CD-2025-036
February 21 , 2025
Page 73 of77
~
ENERGY SOLUTIONS
Mr. Doug Hansen
CD-2025-036
February 21 , 2025
! _::1
w -1.0.
: ................. ~--...... :. ·:JI. ....... " .... " " • • ........ .
~ = = -=
:::1 '"'~irr'ffmYfTfITTwr~~ ... I
0 ~ = = -= distance (m)
Figure 3-18. Transect SSSU.
1~1 ........... *-~~~~ ..... ~ ...... , I
0 so 100 150 200 250 300
0.0 I ~11111111111111111111 •• , ... 11111•"4 ~ll ... ~--~JJ../111• ·--· ~I,# -Utl it!: tt#ifll11 I I .-,11111 1 Ill JU Ill
-0.2
-0.4
-0.6
~--~. --~-~--~--~--.,.---'.
~ = = -= -
~1 ~~t@"I
0 50 100 150 200 250 300
distance (m)
Figure 3-1 9. Transect SSSL.
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
[phone] • Fax: [phone] • www.energysolutions.com
Page 74 of77
References
~
ENERGY SOLUTIONS
Mr. Doug Hansen
CD-2025-036
February 21 , 2025
Barnhart, K.R., et al., 2020. Inverting Topography for Landscape Evolution Model Process
Representation: 3. Determining Parameter Ranges for Select Mature Geomorphic Transport Laws
and Connecting Changes in Flu vial Erodibility to Changes in Climate, Journal of Geophysical
Research: Earth Surface I 25 ( e20 I 9JF005287) 1-34 doi: I 0. I 029/20 I 9JF005287
Evans, K.G., et al., 2000. Post-Mining Landform Evolution Modelling: I. Derivation of Sediment
Transport Model and Rainfall-Runoff Model Parameters, Earth Surface Processes and Landforms
25 (2000) 743-763
Fleming, R.W., and A.M. Johnson, 1975. Rates of Seasonal Creep of Silty Clay Soil, Q. JI Engng Geo! 8
(1975) 1-29
Hancock, G.R., et al., 20 I 1. An Evaluation of Landscape Evolution Models to Simulate Decadal and
Centennial Scale Soil Erosion in Grassland Catchments, Journal of Hydrology 398 (2011) 171 -
183
Hancock, G.R., et al., 2000. Medium-term Erosion Simulation of an Abandoned Mine Site Using the
SIBERIA Landscape Evolution Model, Aust. J. Soil Res. 38 (2000) 249-263
Hancock, G.R., and G.R. Willgoose, 2021. Predicting Gully Erosion Using Landform Evolution Models:
Insights from Mining Landforms, Earth Surface Processes and Landforms 46 (15) 3271 -3290
doi: 10.1002/esp.5234
Hernandez, M., et al., 2017. The Rangeland Hydrology and Erosion Model: A Dynamic Approach for
Predicting Soil Loss on Rangelands, Water Resources Research (53) 1-24 doi:
I 0.1002/20 l 7WR02065 I
Nearing, M.A., 2004. Capabilities and Limitations of Erosion Models and Data, 13th International Soil
Conservation Organisation Conference, Conserving Soil and Water for Society: Sharing
Solutions, Brisbane, July 2004
Neptune, 2021. Embankment Modeling for the Clive DU PA, Clive DU PA Model v2. 0, NAC-0019 _RS,
Neptune and Company Inc., Los Alamos NM, June 2021
Neptune, 2023. Erosion Modeling for the Clive DU PA, Clive DU PA Model v3.0, NAC-00l 7_R7 ,
Neptune and Company Inc., Lakewood CO, December 2023
Press, W.H., et al., 2007. Numerical Recipes, 3rd Edition: The Art of Scientific Computing, Cambridge
University Press, New York NY
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
[phone] • Fax: [phone] • www.energysolutions.com
Page 75 of77
~
ENERGY SOLUTIONS
Mr. Doug Hansen
CD-2025-036
February 21 , 2025
Rogers, V.C ., 2023a. Responses to Federal Cell Facility Application Request for Information -DRC-
2022-023940, DRC-2022-022282, CD-2023-025, EnergySolutions LLC, Salt Lake City UT, May
2023
Rogers, V.C ., 2023b. Federal Cell Facility Application: Responses to the Director 's Request for
lriformation -DRC-2023-066226, CD-2023-249, EnergySolutions LLC, Salt Lake City UT,
December 2023
Smith, C.L., and C.H. Benson, 2016. Influence of Coupling Erosion and Hydrology on the Long-Term
Performance of Engineered Surface Barriers, NUREG/CR-7200, United States Nuclear
Regulatory Commission (NRC), Washington DC, May 2016
Tucker, G.E., 2004. Drainage Basin Sensitivity to Tectonic and Climatic Forcing: Implications of a
Stochastic Model for the Role of Entrainment and Erosion Thresholds, Earth Surface Processes
and Landforms 29 (2004) 185-205
USDA ARS, 2008. Cligen 5.3, U.S. Department of Agriculture, Agricultural Research Service, available
from https://www.ars.usda.gov/midwest-area/west-lafayette-in/national-soil-erosion-
research/docs/wepp/cligen/
USDA ARS SWRC, 2015. Rangeland Hydrology and Erosion Model (RHEM), version 2.3, U.S.
Department of Agriculture, Agricultural Research Service, Southwest Watershed Research Center,
available from https://dss.tucson.ars.ag.gov/rhern/
Willgoose, G., 2005. User Manual for SIBERIA (Version 8.30), Telluric Research, Scone Australia
Willgoose, G., 2011. SIBERIA, version 8.3.3, available from https://doi.org/10.1594/IEDA/100163
Willgoose, G., et al., 1991 a. A Coupled Channel Network Growth and Hillslope Evolution Model , 2.
Nondimensionalization and Applications, Water Resources Research 27 (7) 1685-1696
Willgoose, G., et al., 1991b. A Coupled Channel Network Growth and Hillslope Evolution Model, 1.
Theory, Wa ter Resources Research 27 (7) 1671 -1684
Wong, M., and G. Parker, 2006. Reanalysis and Correction of Bed-Load Relation of Meyer-Peter and
Muller Using Their Own Database, Journal of Hydraulic Engineering 132 ( 11 )
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
[phone] • Fax: [phone] • www.energysolutions.com
Page 76 of77
~
E ERGYSOLUTJONS
Mr. Doug Hansen
CD-2025-036
February 21 , 2025
If you have further questions regarding these responses to the director's requests of DRC-2024-0048 J 7
and DRC-2024-006804, please contact me at (80 I) 649-2000.
Sincerely,
Vern C.
Rogers
Vern C. Rogers
D,g1taUy signed by Vern C. Rogers
DN: cn=Vern C. Rogers,
o=EnergySolutions, ou=Waste
Management Division,
email•[e-mail],
c=US
Date: 2025.02.21 12:09:49-07'00'
Director, Regulatory Affairs
enclosure
I certify under penalty of law that this document and all attachments were prepared under my direction or supervision in accordance with a system
designed to assure that qualified personnel properly gather and evaluate the information submitted Based on my inquiry of the person or persons
who manage the system, or those persons directly responsible for gathering the information, the information submitted is, to the best of my
knowledge and belief. true, accurate, and complete. I am aware that there are significant penalties for submittingfalse information, including the
possibility of fine and imprisonment for knowing violations.
299 South Main Street, Suite 1700 • Salt Lake City, Utah 84111
[phone] • Fax: [phone] • www.energysolutions.com
Page 77 of77