IDENTIFYING TOPOGRAPHIC SIGNATURES OF RIVER NETWORK
DEVELOPMENT IN VARIED TERRAINS OF THE PACIFIC NORTHWEST,
USA
by
NATHANIEL THOMAS KLEMA
A DISSERTATION
Presented to the Department of Earth Sciences
and the Division of Graduate Studies of the University of Oregon
in partial fulfillment of the requirements
for the degree of
Doctor of Philosophy
September 2023
DISSERTATION APPROVAL PAGE
Student: Nathaniel Thomas Klema
Title: Identifying Topographic Signatures of River Network Development in Varied
Terrains of the Pacific Northwest, USA
This dissertation has been accepted and approved in partial fulfillment of the
requirements for the Doctor of Philosophy degree in the Department of Earth
Sciences by:
Leif Karlstrom Chair
Joshua Roering Core Member
Kathy Cashman Core Member
Sarah Cooley Institutional Representative
and
Krista Chronister Vice Provost for Graduate Studies
Original approval signatures are on file with the University of Oregon Division of
Graduate Studies.
Degree awarded September 2023
2
© 2023 Nathaniel Klema
This work is licensed under a Creative 
Commons Attribution License 
3
DISSERTATION ABSTRACT
Nathaniel Thomas Klema
Doctor of Philosophy
Department of Earth Sciences
September 2023
Title: Identifying Topographic Signatures of River Network Development in Varied
Terrains of the Pacific Northwest, USA
Geomorphology leverages the competition between processes of uplift and
erosion to infer geologic time-scale forcing on the earth’s upper crust. While this
framework has revolutionized our understanding of global mountain building, it
has mostly been applied to relatively simple tectonic landscapes and few tools exist
for applying such analysis to more complex volcanic terrains. Such settings, and
particularly continental volcanic arcs, make up some of the earths most geologically
active landscapes and are thus of broad scientific interest (Hilley et al. (2022)).
In this work we focus on varied terrains of the pacific northwestern region of the
USA to develop tools and workflows that can applied to the study of the Cascade
Volcanic arc.
In Chapter II we derive a land form classification scheme based on classical
surface theory that can be used as an objective means of studying topography. We
focus this analysis on the Oregon Coast range, a region studied for it’s assumed
steady-state characteristics. We show the efficacy of our curvature-based approach
in identifying process regimes that define the evolution of fluvial landscapes, and
show how land surface geometry varies as a function of drainage area.
4
In Chapter III we apply the methods used in Chapter II to probe timescales
of fluvial development in the Central Oregon Cascades. Here a well known state-
shift in hydrology is driven by weathering processes in aging volcanic bedrock and
so the degree of fluvial development can inform both bedrock age and hydrologic
regime (Jefferson, Grant, Lewis, and Lancaster (2010)). We show that this
evolution is marked by an increase in the magnitude of land surface curvatures
providing an objective means of studying this state-shift that does not depend
on modern catchment delineations. We show that a full transition from volcanic
constructional to fluvially dissected topography occurs over ∼ 1 My, which is
consistent with past estimates (Jefferson et al. (2010)).
In Chapter IV we use an interdisciplinary approach to show how magmatism
is expressed in mountain range scale land forms of the Columbia River Gorge where
a cross section of the Cascade Arc has been exposed by persistent incision of the
Columbia River. We combine modelling of fluvial knickpoints with a flexural uplift
model applied to a geologic structural datum to show the relative influence of
intrusive and extrusive magmatism to land form development. We also show how
our model for the magmatic origin of this section of the Cascades is supported by
heat flow, bouguer gravity, and surface wave tomography data.
Together the chapters of this dissertation show how through
multidisciplinary approaches geomorphology can aid in disentangling the landscape
evolution history of complex volcanic terrains.
5
ACKNOWLEDGEMENTS
I would like to thank my entire committee, and specifically Leif Karlstrom,
for the incredible inspiration and mentorship through this process. I would also
like to thank my parents, Barbara and Thomas Klema, for instilling in me passion
and curiosity for the natural world, and to my siblings Matthew Klema, Julia
Klema, and Mariah Richstone for teaching me how to explore that passion and
curiosity and supporting me along the way. I am ever grateful for Laurie Williams,
Randy Palmer, Gerald Crawford, and Gary Gianinni for helping me build the
necessary skills, and Karl Karlstrom, Gene Humphreys, Jim O’Connor, Charles
Cannon, Gordon Grant, Ray Wells, Chengxin Jiang, and Brandon Schmandt for
the conversations and guidance along the way. There are too many people to
name who have been integral to this journey but to name a few in no particular
order Trevor Doty, Ben Luck, Somer Morris, Louis Geltman, Brooke Hunter,
Eric Levinson, Jamie Wright, Eric Tomczak, Darby McAdams, Cooper Lambla,
Charles King, Kristin Alligood, the Thunderdome, Allison Kubo, Will Struble, Dan
O’Hara, John Ditmars, The Tenney family, Bozo Cardozo, Andrew McEwen, Middy
Tilghman, Chris Harper, Zebulon, SZ4Grads, Jim Snyder, Ed Skrzypkowski, Marla
Trox, Sandy Thoms, Dave Stemple, and the Karlstrom and Roering lab groups
were all essential to getting me to the finish line. Last but definitely not least thank
you to Kira Tenney for her patience, love, and support. I could not have done it
without you.
6
For my parents Barbara and Thomas Klema. Thank you for giving me the curiosity
and creativity I needed to reach this point.
7
TABLE OF CONTENTS
Chapter Page
I. INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . 25
II. THE INTRINSIC AND EXTRINSIC GEOMETRY OF
LANDSCAPES . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.1. Connecting Landscape Form to Process . . . . . . . . . . . . . 29
2.2. Connecting topographic curvature to classical surface theory . . . . 35
2.3. Oregon Coast Range . . . . . . . . . . . . . . . . . . . . . . 42
2.4. Calculation of curvature . . . . . . . . . . . . . . . . . . . . 42
2.4.1. Deriving the first and second fundamental forms . . . . . . 42
2.4.2. Finding the principal curvatures . . . . . . . . . . . . . 44
2.4.3. Defining the Mean and Gaussian Curvatures and
Surface Shape Classes . . . . . . . . . . . . . . . . . . 48
2.5. Calculation of Curvature on Gridded Data Sets . . . . . . . . . . 53
2.5.1. Spectral Filtering . . . . . . . . . . . . . . . . . . . . 54
2.5.1.1. Minimizing Long-Wavelength Contamination . . . . 56
2.5.1.2. Calculation of amplitude spectrum . . . . . . . . 57
2.5.1.3. Removal of High-Frequency Noise . . . . . . . . . 59
2.5.2. Scale Selection . . . . . . . . . . . . . . . . . . . . . 61
2.6. Alignment of Slope and Principal Curvature Vectors . . . . . . . . 62
2.7. Connections to the Laplacian . . . . . . . . . . . . . . . . . . 67
2.8. Land Surface Shape-Class Distribution . . . . . . . . . . . . . . 72
2.9. Geometric View of the Oregon Coast Range . . . . . . . . . . . 74
2.9.1. Gaussian curvature-based landscape partitioning . . . . . . 76
2.9.1.1. Region 1: Drainage areas less that 3.3 ×
102 m2 . . . . . . . . . . . . . . . . . . . . 76
8
Chapter Page
2.9.1.2. Region 2: Drainage areas between 3.3 ×
102 m2 and 1.3× 103 m2 . . . . . . . . . . . . . 78
2.9.1.3. Region 3: Drainage areas between 1.3 ×
103 m2 and 4.5× 103 m2 . . . . . . . . . . . . . 80
2.9.1.4. Region 4: Drainage areas between 4.5 ×
103 m2 and 9.8× 105 m2 . . . . . . . . . . . . . 82
2.9.1.5. Region 5: Drainage areas greater than
9.8× 105 m2 . . . . . . . . . . . . . . . . . . 84
2.9.2. Mean curvature based landscape partitioning . . . . . . . . 86
2.9.2.1. Mean curvature inflection point . . . . . . . . . . 86
2.9.2.2. Mean curvature based channel
threshold: Drainage areas greater than
1× 105 m2 . . . . . . . . . . . . . . . . . . . 90
2.10. Other Potential Applications . . . . . . . . . . . . . . . . . . 92
2.10.1. Identifying debris-flow source regions . . . . . . . . . . . 92
2.11. Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . 94
2.11.1. Bridge . . . . . . . . . . . . . . . . . . . . . . . . . 94
III. CURVATURE SIGNATURES OF VOLCANIC BEDROCK
AGING IN THE CENTRAL OREGON CASCADES . . . . . . . . . 97
3.1. The Cascades Chronosequence . . . . . . . . . . . . . . . . . 98
3.1.0.1. Groundwater Dominated High Cascades . . . . . . 98
3.1.0.2. Surface Water Dominated Western Cascades . . . . 101
3.2. Signatures of Topographic Development . . . . . . . . . . . . . 101
3.2.1. Curvature based surface classification . . . . . . . . . . . 104
3.2.2. Timescales of Fluvial Development . . . . . . . . . . . . 105
3.3. Further Directions . . . . . . . . . . . . . . . . . . . . . . . 107
3.4. Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . 110
3.5. Bridge . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
9
Chapter Page
IV. THE MAGMATIC ORIGIN OF THE COLUMBIA RIVER
GORGE, USA . . . . . . . . . . . . . . . . . . . . . . . . . . 112
4.1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . 112
4.2. The Columbia River Gorge . . . . . . . . . . . . . . . . . . . 112
4.3. Transcrustal magmatic structure . . . . . . . . . . . . . . . . 115
4.4. Long-term balance of intrusions versus eruptions . . . . . . . . . 116
4.5. Topographic response to magmatism . . . . . . . . . . . . . . . 120
4.6. Entwined upper crustal deformation, arc mountain
building, and magma ascent . . . . . . . . . . . . . . . . . . 123
APPENDIX: APPENDIX A - CHAPTER IV SUPPLEMENARY
MATERIAL . . . . . . . . . . . . . . . . . . . . . . . . . . . 126
A.0.1. Structure and timing of Gorge deformation . . . . . . . . . 126
A.0.2. Topographic analysis and fluvial erosion . . . . . . . . . . 129
A.0.3. Interpretation of Bouguer gravity anomaly . . . . . . . . . 134
A.0.4. Elastic deformation/knickpoint response model . . . . . . . 137
A.0.5. Isostatic response to erosion by the Columbia River . . . . . 141
A.0.6. Heat flow . . . . . . . . . . . . . . . . . . . . . . . 147
REFERENCES CITED . . . . . . . . . . . . . . . . . . . . . . . . 154
10
LIST OF FIGURES
Figure Page
1. Figure from Hack, Seaton, and Nolan (1957)
showing the power-law relationship between stream
length and upstream drainage area, now known to
as Hack’s Law, in several catchments in Maryland
and Virginia. . . . . . . . . . . . . . . . . . . . . . . . . . . 34
2. Error in curvature defined as second partial
derivative as a function of slope along a 1-d curve. . . . . . . . 38
3. Difference between distances and angles measured
on a map projection versus on the surface. A. Map
projection of DEM including map grid defined by E-W and
N-S lines with grid spacing dx. The red line corresponds
to the rectangular outline of panel B. B. DEM viewed as a
2-d manifold embedded in a 3-d space. Dashed lines how a
locally defined u-v coordinate system that follows x and y
curves on the map projection, but which are not orthogonal
or of equal length due to surface distortion. E and G are
coefficients of the first fundamental form, and ds is the
displacement vector that results moving one grid space
along each of these coordinate vectors . . . . . . . . . . . . . . . . 40
4. Curvature as a function of the angle between the
osculating plane and a reference direction. Here the
angle θ = 0 is associated with the first principal curvature
(k1). k2 is the second principal curvature and KM is the
Mean curvature about which the value oscillates. . . . . . . . . . . . 49
5. Shape classes into which points on the surface can
be sorted based on the signs on the Mean (KM) and
Gaussian (KG) curvatures. In this analysis we focus on
those classes that can be assigned based on raw curvature
values, which are synformal saddles, antiformal saddles,
basins, and domes, and do not include developable surfaces
or perfect saddles. . . . . . . . . . . . . . . . . . . . . . . . . 52
11
Figure Page
6. Interpolation artifacts in 1-m LiDAR data. A.
The location of the LiDAR swath within the study area.
B. LiDAR hillshape showing regions dominated by
interpolation artifacts over length scales of hundreds of
meters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
7. Comparison of windowing techniques used in
preprocessing DEMs. A. Raw DEM. B. Mirrored
data with Tukey window. The red box is the orignal DEM
extent. C. Mirrored with hanning window. D. Hanning
window applied to DEM without mirroring. E. Profiles of
topography for each of the preprocessing methods compared
to the original DEM. . . . . . . . . . . . . . . . . . . . . . . . 58
8. Identification of gridding artifacts in the 8.1 m DEM
data used in this study. A. Mean squared amplitude
vs. radial frequency of topography. Colored boxes highlight
frequency ranges shown in panels C, D, and E. B. Raw
DEM. C.DEM bandpass filtered to admit wavelengths
between 30 and 50 m showing the gibbs phenomena
thought to produce the spike in the power spectrum at
these wavelengths. E. DEM high-pass filtered to admit
wavelengths less than 30 m. Straight line are artifacts from
the gridding process. . . . . . . . . . . . . . . . . . . . . . . . 60
9. Surface geometry metrics binned by upstream
drainage area for a range of low-pass filter cut-
offs between 50 m and 500 m calculated on 50 m
intervals. A. Gradient. B. Gaussian curvature (KG).
C. Mean Curvature (KM). D. First principal curvature
(k1). E. Second principal curvature (k2). F. Surface area
proportionality factor (α) . . . . . . . . . . . . . . . . . . . . . 63
10. Alignment between principal curvatures and slope
of the tangent plane measured via the surface
normal N. A. DEM low-pass filtered to 200 m. B. Slope
of tangent plane. C. Scalar product of the slope vector and
first principal curvature vector (k1).D.Scalar product of the
slope vector second principal curvature vector (k2). . . . . . . . . . . 66
12
Figure Page
11. Principal curvatures binned by upstream drainage
area for DEM filtered to 200 m. Left axis measures
mean curvature values shown by solid lines. Right
axis measures the scalar product between the principal
curvatures and slope shown by the dotted lines. . . . . . . . . . . . 67
12. Comparison of different slope metrics binned by
drainage area. . . . . . . . . . . . . . . . . . . . . . . . . . 68
13. Map-view of surface curvature metrics. A.
Topography. B. Gaussian Curvature (KG). C. Difference
between mean curvature calculated using the principal
curvatures and Mean curvature calculated using the
components of the laplacian operator. D. Mean curvature (KM). . . . . 70
14. Surface curvatures binned by upstream drainage area. . . . . . 71
15. Comparison between mean curvature calculated
using the principal curvatures (KM), and mean
curvature calculated using the components of the
laplacian (KMLP). A. Calculated mean curvatures binned
by upstream drainage area. B. Difference between KM
and KMLP binned by KM . C. Difference between KM and
KMLP binned by slope . . . . . . . . . . . . . . . . . . . . . . 73
16. Shape class distributions as a function of drainage
area. A. Probability density function of upstream drainage
areas within the study area. B. Probability density
functions of shape classes as a function of drainage area.
C. Conditional PDFs of shape classes at a given drainage area. . . . . . 75
17. Compilation of surface geometry data used in this
study. A. PDF of upstream drainage area. B. Slope
metrics binned by upstream drainage area. C. Surface
curvatures binned by upstream drainage area. D. Principal
curvature values and slope alignment binned by upstream
drainage area. E. Conditional PDFs of shape classes as a
function of upstream drainage area. F. Map-view of shape
class distribution. G. PDFs of slope metrics. H. Pie chart
distribution of shape classes. . . . . . . . . . . . . . . . . . . . . 77
13
Figure Page
18. Surface geometry data for points in the landscape
with upstream drainage areas less than 3.3 × 102
m2. The red box in each of the area plots highlights the
region of interest. A. PDF of upstream drainage area. B.
Slope metrics binned by upstream drainage area. C. Surface
curvatures binned by upstream drainage area. D. Principal
curvature values and slope alignment binned by upstream
drainage area. E. Conditional PDFs of shape classes as a
function of upstream drainage area. F. Map-view of shape
class distribution. G. PDFs of slope metrics. H. Pie chart
distribution of shape classes. . . . . . . . . . . . . . . . . . . . . 79
19. Surface geometry data for points in the landscape
with upstream drainage areas between 3.3 × 102
m2 and 1.3 × 103 m2. The red box in each of the
area plots highlights the region of interest. A. PDF of
upstream drainage area. B. Slope metrics binned by
upstream drainage area. C. Surface curvatures binned
by upstream drainage area. D. Principal curvature values
and slope alignment binned by upstream drainage area. E.
Conditional PDFs of shape classes as a function of upstream
drainage area. F. Map-view of shape class distribution. G.
PDFs of slope metrics. H. Pie chart distribution of shape classes. . . . 81
20. Surface geometry data for points in the landscape
with upstream drainage areas between 1.3 × 103
m2 and 4.5 × 103 m2. The red box in each of the
area plots highlights the region of interest. A. PDF of
upstream drainage area. B. Slope metrics binned by
upstream drainage area. C. Surface curvatures binned
by upstream drainage area. D. Principal curvature values
and slope alignment binned by upstream drainage area. E.
Conditional PDFs of shape classes as a function of upstream
drainage area. F. Map-view of shape class distribution. G.
PDFs of slope metrics. H. Pie chart distribution of shape classes. . . . 83
14
Figure Page
21. Surface geometry data for points in the landscape
with upstream drainage areas between 4.5 × 103
m2 and 9.8 × 105 m2. The red box in each of the
area plots highlights the region of interest. A. PDF of
upstream drainage area. B. Slope metrics binned by
upstream drainage area. C. Surface curvatures binned
by upstream drainage area. D. Principal curvature values
and slope alignment binned by upstream drainage area. E.
Conditional PDFs of shape classes as a function of upstream
drainage area. F. Map-view of shape class distribution. G.
PDFs of slope metrics. H. Pie chart distribution of shape classes. . . . 85
22. Surface geometry data for points in the landscape
with upstream drainage areas greater than 9.8 × 105
m2. The red box in each of the area plots highlights the
region of interest. A. PDF of upstream drainage area. B.
Slope metrics binned by upstream drainage area. C. Surface
curvatures binned by upstream drainage area. D. Principal
curvature values and slope alignment binned by upstream
drainage area. E. Conditional PDFs of shape classes as a
function of upstream drainage area. F. Map-view of shape
class distribution. G. PDFs of slope metrics. H. Pie chart
distribution of shape classes. . . . . . . . . . . . . . . . . . . . . 87
23. Surface geometry data for points in the landscape
where the Mean curvature (KM) is majority
positive. The red box in each of the area plots highlights
the region of interest. A. PDF of upstream drainage area.
B. Slope metrics binned by upstream drainage area. C.
Surface curvatures binned by upstream drainage area. D.
Principal curvature values and slope alignment binned by
upstream drainage area. E. Conditional PDFs of shape
classes as a function of upstream drainage area. F. Map-
view of shape class distribution. G. PDFs of slope metrics.
H. Pie chart distribution of shape classes. . . . . . . . . . . . . . . 88
15
Figure Page
24. Surface geometry data for points in the landscape
where the Mean curvature (KM) is majority
negative. The red box in each of the area plots highlights
the region of interest. A. PDF of upstream drainage area.
B. Slope metrics binned by upstream drainage area. C.
Surface curvatures binned by upstream drainage area. D.
Principal curvature values and slope alignment binned by
upstream drainage area. E. Conditional PDFs of shape
classes as a function of upstream drainage area. F. Map-
view of shape class distribution. G. PDFs of slope metrics.
H. Pie chart distribution of shape classes. . . . . . . . . . . . . . . 89
25. Distribution of surface geometry metrics for regions
defined by the Mean curvature inflection point.
A. Map-view of landscape partitioned about the point of
inflection in KM . B. Distribution of tangent slope in regions
of both negative and positive average KM . C. Distribution
of Mean curvature in regions of both negative and positive
average KM . D. Distribution of gaussian curvatures in
regions of both negative and positive average KM . . . . . . . . . . 91
26. Surface geometry data for points in the landscape
with upstream drainage area greater than 1× 105 m2
as determined by the global minimum in the mean
curvature. The red box in each of the area plots highlights
the region of interest. A. PDF of upstream drainage area.
B. Slope metrics binned by upstream drainage area. C.
Surface curvatures binned by upstream drainage area. D.
Principal curvature values and slope alignment binned by
upstream drainage area. E. Conditional PDFs of shape
classes as a function of upstream drainage area. F. Map-
view of shape class distribution. G. PDFs of slope metrics.
H. Pie chart distribution of shape classes. . . . . . . . . . . . . . . 93
27. Map of colluvial hollows inferred from our analysis. Purple
outlines mark basin structures that only containing points
with areas less than an assigned threshold of 1 × 105 m2.
These are overlain on a shape-class map of region 4, which
we interpret to represent the uppermost reaches of the
channel network . . . . . . . . . . . . . . . . . . . . . . . . . 95
16
Figure Page
28. Map of the central Oregon Cascades. A. Topography.
Red boxes indicate areas of focused analysis in Figs. 30 and
33. B. Age of volcanic bedrock units (from Sherrod and
Smith (2000)) . . . . . . . . . . . . . . . . . . . . . . . . . . 99
29. Characteristic hydrographs of the spring dominated
High Cascades, and the surface run-off dominates
Western Cascades (from Deligne et al. (2017) which is a
modified figure from Tague and Grant (2004) . . . . . . . . . . . . 100
30. Compilation of surface geometry data used in
Chapter II applied to the Western Cascades. A.
PDF of upstream drainage area. B. Slope metrics binned
by upstream drainage area. C. Surface curvatures binned
by upstream drainage area. D. Principal curvature values
and slope alignment binned by upstream drainage area. E.
Conditional PDFs of shape classes as a function of upstream
drainage area. F. Map-view of shape class distribution. G.
PDFs of slope metrics. H. Pie chart distribution of shape classes. . . . 102
31. Plot of drainage density versus bedrock age in the
upper McKenzie River Basin (from Jefferson et al. (2010)) . . . . . 103
32. Maps of surface bedrock age in the Oregon
Cascades averaged over different spatial scales. . . . . . . . . . 104
33. Surface metric distributions for sections of the
Oregon Coast Range, Western Cascades, and
High Cascades. A. Oregon Coast Range topography.
B. Western Cascades topography. C. High Cascades
topography. D. Gaussian curvature PDFs. E. Mean
curvature PDFs. F. Slope PDFs. G. α PDFs. . . . . . . . . . . . . 106
34. Gaussian curvature of the Cascades chronosequence.
A. Map-view of calculated Gaussian curvature. B. Gaussian
curvature binned by bedrock age. C. Gaussian curvature
binned by UTM easting. . . . . . . . . . . . . . . . . . . . . . 108
35. Mean curvature of the Cascades chronosequence. A.
Map-view of calculated Mean curvature. B. Mean curvature
binned by bedrock age. C. Mean curvature binned by UTM easting. . . 108
36. Tangent slope of the Cascades chronosequence. A.
Map-view of calculated tangent slope. B. Tangent slope
binned by bedrock age. C. Tangent slope binned by UTM easting. . . . 109
17
Figure Page
37. α of the Cascades chronosequence. A. Map-view
of calculated area proportionality factor (α). B. Area
proportionality factor (α) binned by bedrock age. C. Area
proportionality factor (α) binned by UTM easting. . . . . . . . . . . 109
38. Overview of the Columbia River Gorge. A. Structural
setting, showing GNSS velocities with respect to stable
North American reference frame (McCaffrey, King, Payne,
and Lancaster (2013)) (black arrows), Yakima Folds (orange
lines), the modern and paleo (‘Bridal Veil’) Columbia River
(blue and pink lines), and extent of the Cascades arc (red
dashes). B. Symbols are Bridal Veil Channel hyaloclastite
deposits and best-fitting knickpoint locations in Columbia
tributaries. Maximum topographic elevations within 10 km
south of the Columbia (black line) with prominent faults
(dashed lines) (Mcclaughry, Wiley, Conrey, Jones, and Lite
(2012)). C. Southward orthoview of fluvial catchments,
hyaloclastite deposits, and major faults. D. Orthoview of
Oneonta Creek showing the location of the trunk stream
knickpoint used in this study (large blue square), analogous
knickpoints in smaller tributary channels, and waterfalls
associated with layered basalt stratigraphy (black stars).
E. Regional map of the Columbia river drainage network,
showing the Euler pole for Cascadia forearc block relative to
North America. Orange lines are compressional structures
while yellow lines mark arc-marginal normal faults. F.
Longitudinal profile using area-normalized upstream
distance χ for Oneonta creek showing the location of the
best fitting slope-break knickpoint (Supplement). Lithologic
units are labeled. . . . . . . . . . . . . . . . . . . . . . . . . 114
18
Figure Page
39. Geophysical anomalies across the Gorge and
their relationships to BV channel deformation.
A. Best-fitting uplift model from maximum likelihood
estimation with 95% confidence model envelope. Light
blue squares show knickpoint elevations inferred from a
transient landscape evolution model forced by flexural uplift
constrained by hyaloclastite deformation (pink squares). B.
Interpolated Bouguer gravity and surface heat flow spanning
the Gorge. C. Average surface wave velocity anomaly.
The white polygon marks inferred 2.5 − 4% partial melt
(Supplement). Black dots are earthquake hypocenters for
all events underlying the Gorge in the Pacific Northwest
Seismic Network catalog since 1970. D. Hyaloclastite
elevation (squares) and knickpoint elevation (diamonds)
vs knickpoint location χkn. Red dashed line shows uplift
predicted by fluvial model using slope-area regression, while
shaded swath shows model predictions from a coupled
flexural uplift/landscape evolution model (Supplement).
Mismatch between χkn and knickpoint elevation likely
implies time-evolving hydraulic conductivity/erodibility of
young lava flows. E. Regional map of shear wave seismic
velocity at 30 km depth showing swath in panel C. . . . . . . . . . . 117
40. Inference of transcrustal magmatic structure. A.
Volume weighted Gaussian kernel density of Quaternary
edifices within 10 km of the study area. B. Intrusive and
extrusive magmatic topography, using best fitting flexural
uplift model. Dashed lines indicate major fault systems that
bound the Hood River Valley (Mcclaughry et al. (2012)).
C. Best fitting plate model defining elastic thickness of
the upper crust. Alternating tensile and compressional
deviatoric stresses in the plate are inferred to direct shallow
magma ascent and faulting. D. Ratios of intrusive:extrusive
(I:E) magmatism calculated using topography (assuming
north-south continuity in flexural uplift), knickpoint
positions, and total calculated volumes of erupted and
intruded magma. . . . . . . . . . . . . . . . . . . . . . . . . . 119
19
Figure Page
A.41.Photos of hyaloclastite deposits associated with
the BV channel. A. Hyaloclastite (orange/yellow) with
interbedded pillow complexes seen from I-84 about 4 km
west of Hood River, OR. near Ruthton Point B. Lens of
hyaloclastite material with channel capping LKT flows in
Ainsworth State Park near Dodson, OR. C. Exposures of
hyaloclastite material uplifted nearly 1000 m above the
modern Columbia River (visible on right side of panel) seen
from the Pacific Crest Trail near Cascades Locks, OR. . . . . . . . . 128
A.42.Raw Bridal-Veil Channel elevation data. A. Mean
elevations of hyaloclastite deposits and minimum elevations
of LKT flows that overlie the BV Channel. B. Map view of
hyaloclastite and LKT polygons. . . . . . . . . . . . . . . . . . . 129
A.43.16 Columbia River tributaries used in this study.
A. Blue lines are longitudinal profiles extracted from 3m
LIDAR data plotted against an area-normalized χ-scale
calculated using best fitting concavity value of θ = 0.455
and A 2o = 1 m . Red dashed lines are piece-wise linear fits
that best capture slope-break knickpoints with best fitting
knickpoint locations marked by blue squares on each profile.
B. PDF of θ values calculated using best-fitting knickpoint
decompositions for each value. . . . . . . . . . . . . . . . . . . 132
A.44.Slope-area plot for fluvial channels downstream
of the knickpoint. Points represent raw slope-area
data for all points within fluvial channels downstream
of paleosurface elevations implied by the flexure model.
The blue line and shaded region represents a model fit to
equation A.1 using the range of ks values from the coupled
grid-search inversion and concavity values derived in the
knickpoint fitting process. The red dashed line is a fit to the
raw slope-area data using the robustfit function in Matlab. . . . . . . 134
A.45.Bouguer gravity and heat flow data. A. Bouguer
gravity anomaly (Bonvalot et al. (2012)) with locations of
quaternary volcanic vents, the modern Columbia River,
and the extent used to generate the cross-arc gravity profile
shown in Figure 39.B. B. Heat flow data from Williams and
DeAngelo (2008) (C. F. Williams and DeAngelo (2008))
interpolated onto a 2-minute grid showing measurement stations. . . . . 135
20
Figure Page
A.46.Interpretation of admittance between gravity and
hyaloclastite deformation. A. Detrended Bouguer
gravity plotted versus hyaloclastite deformation where
the slope of the linear fit between the datasets is the
gravitational admittance over a wavelength of ∼ 80 km.
B. Regime diagram showing expected ranges of admittance
from different loading scenarios. The black line represents
isostatic compensation at the base of the crust. Points
above the line represent a flexural response to buried
loading while points below the line represent a flexural
response to surface loading. . . . . . . . . . . . . . . . . . . . . 136
A.47.Posterior PDFs and covariance plots from maximum
likelihood estimation. D is flexural rigidity, Fnet is
the total force acting on the plate, and ks is the fluvial
proportionality constant that maps hyaloclastite uplift onto
knickpoint positions. Orange bars represent 95% confidence
intervals for parameter values. . . . . . . . . . . . . . . . . . . . 140
A.48.Possible values for effective elastic thickness (Te) as
a function of flexural rigidity (D), Young’s Modulus (E)
and Poisson’s Ratio (ν). A. Variation in Te as a function of
E. The solid black line shows Te for our best-fitting value
of flexural rigidity (D) and the dashed black lines represent
the 95% confidence interval. The blue and red dashed lines
show variation in Te that result from varying Poisson’s
Ratio where blue and red represent ν = 0.3 and ν = 0.1
respectively. B. Variation in Te for all values of D explored
in the maximum likelihood estimation, where the solid black
line represents the best-fitting model and the dashed lines
are 95% confidence intervals. . . . . . . . . . . . . . . . . . . . . 141
A.49.Isostatic uplift from post 3.5 Ma incision implied
by our uplift model. A. Elevation difference between
topography and best fitting flexure model projected north
and south along the arc. B. Erosion implied by north-
south model projection. C. Maximum flexural uplift
from erosional unloading shown in panel B. D. Range of
model fits to hyaloclastite data and range of values for
uplift caused by erosion within the north-south range of
hyaloclastite deposits (shown by dashed black lines). . . . . . . . . . 143
21
Figure Page
A.50.Predicted Vp, Vs, and ρ with respect to the
variations of pressure and temperature condition
for the solid phase using PerpleX, both for basalts
(A-C) and dacites (D-E). The red and yellow stars
define the maximum and minimum value of each parameter
to define a range used in the calculation. The black lines in
B and E separate the 2D space into several domains where
slightly different mineral assemblages are stable. See Table
S2 for the detailed mineral assemblages for each domain. . . . . . . . 148
A.51.Modelled Vs and melt fraction relationship for
basalts (A) and dacites (B). The color of the curves
represents the different aspect ratios assumed during the
calculation, with the textural equilibrium regime of aspects
ratio between 0.1-0.15 highlighted in yellow. The solid lines
are from calculations using the maximum values of elastic
parameters (red stars in Figure A.50) and the dashed lines
are from the minimum values. The grey-shaded region
indicates Vs range observed in this study. . . . . . . . . . . . . . . 149
A.52.Regional maps of Vs at 20 km (A), 25 km (B), 30
km (C), and 35 km (D) depths, respectively. The
red-filled stars denote the volcanoes of Mount St. Helens
(MSH), Mount Adams (MA), and Mount Hood (MH).
The grey-shaded region highlights the corridor where
Vs values are averaged to produce the transect shown in
Figure 39C. The color of the curves represents the different
aspect ratios assumed during the calculation, with the
textural equilibrium regime of aspects ratio between 0.1-0.15
highlighted in yellow. The solid lines are from calculations
using the maximum values of elastic parameters (red stars
in Figure A.50) and the dashed lines are from the minimum
values. The grey-shaded region indicates Vs range observed
in this study. . . . . . . . . . . . . . . . . . . . . . . . . . . 150
22
Figure Page
A.53.Illustrative conductive geotherms with and without
radiogenic heat production (H = 9 × 10−100 W/kg). Curve
colors correspond to assumed surface heat flow, spanning
the range observed in the Gorge. The black curve is Dacite
solidus with 5 weight percent water (Blatter, Sisson, and
Hankins (2017)), computed in PerpleX. The curves are
lower bounds for geothermal gradient as they neglect heat
associated with intrusions, and hydrothermal circulation
in the shallow crust (Ingebritsen, Sherrod, and Mariner
(1989)), supporting our seismic inference of lower-crustal
melt at present. . . . . . . . . . . . . . . . . . . . . . . . . . 153
23
LIST OF TABLES
Table Page
A.1.Composition of typical basalts and dacites at Mount
St. Helens used to calculate seismic velocity. The
basaltic composition is from Blatter et al., 2017, and dacitic
composition is from Leeman et al., (2005). . . . . . . . . . . . . . . 147
A.2.The mineral assemblages stabilized at different P-T
conditions. The number and corresponding domains are
specified in Figure A.50 panels B and E. . . . . . . . . . . . . . . . 151
24
CHAPTER I
INTRODUCTION
Geomorphology is the study of how the geometry of topographic surfaces
relates to the geologic processes that shape them. Much of the modern science is
rooted in the observation that land is sculpted through competition between forces
that uplift the landscape - such as tectonics, volcanism, and mantle dynamics -
and erosional processes that modulate topographic growth by transporting uplifted
material to sediment sinks where it is more gravitationally stable. Most erosion
processes are gradient driven, and so all else equal higher landscape uplift rates,
which generate more relief, cause faster erosion and overall steeper topography
(Whipple, DiBiase, and Crosby (2013)). This basic conceptual framework has
inspired generations of geomorphologists to find quantitative methods for relating
topographic form to erosion rates with the goal of mapping uplift variation and
thus better understanding spatial and temporal patterns of deeper forcing on the
landscape (Wobus et al. (2006)). This in turn requires disentangling effects of
climatic, lithologic, and biologic variation on topographic development and so most
workers have focused on landscapes where some variables can be assumed constant.
As a result there are few applications of geomorphology to terrains built through
combinations of tectonics and magmatism.
Volcanic arcs comprise some of the earths most geologically active terrains,
and have been identified as priority settings for focused research efforts in coming
decades (Hilley et al. (2022)). Here repaving of topography by eruptive products
makes it difficult to disentangle geologic histories using traditional mapping
techniques; thus the development of geomorphology tools that could be applied
to understand geologic histories of arcs could revolutionize our understanding of
25
such settings. One reason it is challenging to do geomorphology studies in volcanic
settings is that most studies surface process models leverage geometric relationships
that are directly associated with a certain process, and there are relatively few
metrics of topographic form that are decoupled from process and thus can be used
to quantify local variation in lithology, hydrology, and uplift rate characteristic of
volcanic settings. The primary goal if this work is to develop quantitative tools and
approaches of topographic analysis that can be used in settings where assumptions
common to most geomorphology studies break down towards use of topography as
a dataset for studying the evolution of the Cascades volcanic arc.
In Chapter II we look at a section of the well studied Oregon Coast Range
and perform a rigorous calculation of topographic surface characteristics that
builds directly on the canonical work by Carl Frederick Gauss (Gauss (1902)).
We demonstrate how this approach provides an abundance of information not
previously used in landscape evolution contexts, and show how the curvature
of the landscape evolves with drainage area thus putting these novel metrics in
the context of standard approaches of geomorphology. In doing so we derive a
quantitative metric of landscape form that is independent of process, but can be
easily linked to process transitions throughout this characteristic fluvial landscape.
While this chapter has no explicit connection to the volcanic landscape of the
Cascade Volcanic Arc we use the Coast Range as a test location for defining the
tool and demonstrate its connections to more standard slope, area, and curvature
metrics common in the literature.
In Chapter III we apply some of the tools derived in Chapter II to the
Oregon Cascades, where a hydrologic state-shift has been connected to a decrease
in permeability over time in basalts sourced from the arc (Tague and Grant
26
(2004)). While the timescale of this state-shift has already been explored using
morphological approaches (Jefferson et al. (2010); Luo, Grudzinski, and Pederson
(2010)), a dependence on drainage density in these approaches limits their
applicability to specific regions where fluvial processes are present and hydraulic
watersheds are plausibly defined by topography. Since young volcanic landscapes,
such as the High Cascades, are groundwater dominated with no consistent
connection between topography and hydrology (Gannett, Lite, Morgan, and Collins
(2001)) we argue that our curvature based approach, which makes no assumptions
about process, provides a more general approach for exploring the complex coupling
of hydrology and topography in volcanic landscapes. We confirm past estimates for
a timescale of ∼ 1 My for this transition to take place.
Chapter IV is an example of how geomorphology can be used to probe
dynamics of arc construction, and focuses on the topographic development of the
Columbia River Gorge. Here we use an approach that combines the analysis of
fluvial knickpoints with flexural modelling of the upper crust as constrained by
deposits of a Columbia River paleo-channel. We correlate the observed deformation
with heat flow, gravity, and seismic tomography data, all of which are consistent
with ongoing mid-crustal magmatism as a driver of uplift in the Columbia River
Gorge. The cross section provided by the Columbia allows us to partition uplift
into relative contributions from intrusive and extrusive magmatism and show
that the longer-wavelength intrusive component exerts a stronger control on the
topographic response, likely due to the higher permeability of erupted deposits.
This dissertation serves as a template for how concepts of geomorphology
can be built upon to provide a means of better understanding complex volcanic
settings, however more broadly is an example of how interdisciplinary research
27
can be used to tackle challenging problems in the earth sciences. The technical
methods used in Chapters II and III combine geomorphology with formal surface
theory, while Chapter IV combines technical analysis of several different datasets,
which together paint a holistic picture of volcanic mountain building that cannot be
discerned from any one of the data sets alone.
28
CHAPTER II
THE INTRINSIC AND EXTRINSIC GEOMETRY OF LANDSCAPES
2.1 Connecting Landscape Form to Process
The majority of research on erosion dynamics has been done in fluvially
dissected landscapes where sediment transport is driven by river networks.
Quantitative models of erosion in these settings are inherently scale dependant with
different models granting different insight. For example, at the scale of a mountain
range average elevations tracks deformation of the crust and individual erosional
features appear as noise about this long wavelength trend. For some land forms at
this scale erosion has been modeled simply as a diffusive process governed by the
heat equation
zt = κ ∇2s z (2.1)
where z is elevation, ∇2z = zxx + zyy is the Laplacian of the surface, and κs
is an empirical diffusivity governing the rate of topographic relaxation (Ruh
(2020); Watts (2001)) (note that in this work we use subscript notation for partial
derivatives, which is consistent with the curvature literature). This approach allows
for the calculation of the landscape relaxation time following orogenic events and
can easily admit basic uplift forcings such as the isostatic response of the crust to
erosional unloading (Watts (2001)). This is useful for understanding long timescale
cycles of uplift and erosion and the connections between ancient tectonic events and
the modern rock record, however does little to put modern topography in context
of these orogenic cycles. Modeling erosion as purely diffusive does little to elucidate
the diverse physical mechanisms at play, and therefore is limited in its predictive
power.
29
On the other end of the spatial scale there is a rich body of literature
focused on building mechanistic models of individual erosion processes at the
outcrop/field scale down to the scale of individual sediment particles. For example
within fluvial channels the popular stream power family of models (Bagnold (1966);
Howard and Kerby (1983); Seidl, Dietrich, and Kirchner (1994)) show the potential
of average shear stresses within rivers to mobilize sediment, while an adjacent
literature focuses on connections between fluid turbulence and the development
of observed bed structures (A. A. Hurst, Anderson, and Crimaldi (2021); Li,
Venditti, Rennie, and Nelson (2022); MacVicar and Roy (2007); Stoesser, Kara,
MacVicar, and Best (2010); Venditti (2007); Venditti et al. (2014); Wrick and
Pasternack (2008)). The last decades have also seen the development of hillslope
sediment transport models with both continuum mechanics (Gabet (2003); Lamb,
Scheingross, Amidon, Swanson, and Limaye (2011); Roering, Kirchner, and
Dietrich (1999); Roering, Perron, and Kirchner (2007)) and statistical mechanics
(Furbish, Haff, Dietrich, and Heimsath (2009); S. G. Williams and Furbish (2021))
approaches, as well as models describing landsliding (Booth, Roering, and Rempel
(2013)), debris flow processes (Struble, Mcguire, Mccoy, Barnhart, and Marc
(2023)), gully retreat (Flores-Cervantes, Istanbulluoglu, and Bras (2006)), plunge-
pool formation (Scheingross and Lamb (2017)) and a host of other location specific
modes of sediment mobilization and transport.
Outcrop scale models and observations are useful because they can
calibrated with laboratory scale physical models (Baynes et al. (2018); Deshpande,
Furbish, Arratia, and Jerolmack (2021); Gabet (2003); Lamb and Dietrich (2009);
Pasternack, Ellis, Leier, Vallé, and Marr (2006)), and also be connected to
locally derived erosion rates using tools like thermochronology (Nie et al. (2018);
30
Simoes, Braun, and Bonnet (2010)) and cosmogenic nuclide analysis (Gayer,
Mukhopadhyay, and Meade (2008); Hippe et al. (n.d.)). In addition they have
provided the foundation for a growing ability to quantify the dynamic stability
of landscapes towards the forecasting of geologic hazards. The weakness of this
approach is that individual mechanistic models are difficult to connect across
process regimes. Studies tend to compartmentalize erosion into isolated regions that
don’t reflect the continuum nature of landscapes, or capture multiscale feedbacks
that may play a large roll in topographic development.
Between the observational scales defined above is a range of scales at which
fluvial landscapes are defined by dendritic networks of fluvial channels divided by
ridge networks to which they are connected by hillslopes and first order channels.
These are the scales at which it is easiest to connect topographic form to the
temporal evolution of the landscape and deep uplift forcings, which is often done
through implementation of landscape evolution equations of the basic form
zt = U − E (2.2)
where z = z(x, y, t) is elevation, U = U(x, y, t) is uplift rate, and E = E(x, y, t)
is the erosion rate Braun and Willett (2013); Whipple and Tucker (1999). In many
studies erosion is broken into two terms defining fluvial and hillslope processes.
The fluvial component is generally defined along a river longitudinal profile using
a version of the stream power model written as E(s) = KA(s)m|S(s)|n where
s = s(x, y) is an along stream coordinate, K,m, and n are empirical constants
that capture effects of lithology, hydrology and climate, A(s) is drainage area
upstream of a point in the channel (itself a proxy for discharge), and S(s) is
along-channel slope (Whipple and Tucker (1999)). The hillslope flux is usually
derived by defining a relationship between slope and sediment flux. While steep
31
hillslopes are best modeled using a non-linear relationship (Roering et al. (1999,
2007)), the simplest models assume sediment flux is a linear function of slope in
which case hillslope erosion can be written as E = ∇ · q = D∇2z where q is
slope dependent sediment flux and D is a diffusivity that quantifies the efficiency
of hillslope transport processes (Barnhart et al. (2020)). Note that in this case
curvature ∇2z is defined as positive downwards, which differs from much of the
geomorphology literature (Roering et al. (2007)) but will be consistent with the
derivation of invariant curvature metrics in section 2.3.
With the above definitions of erosion equation 2.2 can be written in terms of
competition between uplift and erosion as
zt = U −KAm|∇z|n −D∇2z. (2.3)
Modelling at the landscape scale has proven powerful for understanding timescales
of topographic development, and factors that determine how uplift (or climate)
perturbations propagate through the landscape. It has the benefit of being able
to encompass transitions in process dominance between channels and hillslopes,
but has the disadvantage of blurring over transitional processes that are likely
important drivers of landscape adjustment. Landscape evolution governed by
equation 2.3 depends strongly on parameter values that are difficult to constrain
and likely not constant throughout the landscape as well as in time (Snyder,
Whipple, Tucker, and Merritts (2000)). Thus landscape evolution model results are
nearly always non-unique, and less mechanistic than may be desired in many cases:
it is challenging to calibrate semi-empirical fluvial and hillslope erosion laws over
the full range of relevant scales (Venditti, Li, Deal, Dingle, and Church (2019)).
A commonality between landscape evolution models at all scales, is that
sediment transport rates are usually assumed to be recorded in metrics of land
32
surface geometry, such as slope and curvature. This builds on ideas established
in the late 19th century when work by G.K. Gilbert and W.M. Davis suggested
connections between hillslope convexity and rates and styles of denudation in
mountain terrains (Davis (1892); Gilbert (1877)). These early observations, and
subsequent papers (Fenneman (1908); Gilbert (1908)) showed that both slope and
curvature near drainage divides is related to erosion rates. These early works noted
spatial variation in the dominance of different erosion processes and began the
effort, ongoing today, to partition the landscape into areas defined by particular
mechanisms of erosion.
It is generally accepted that fluvial networks are organized by the
accumulation of drainage area (and thus discharge), which gives rise to a rich and
complex tiling of drainage basins that can appear to have sometimes fractal (Claps
and Oliveto (1996)) or even periodic structure (Perron, Kirchner, and Dietrich
(2009)). This quality was first defined by John Tilton Hack (Hack et al. (1957))
who showed that the lengths of stream channels in select regions of Virginia and
Maryland are at any point in the channel related to upstream drainage area (A) via
the power-law expression L = cAb where c and b are empirical constants (Figure 1).
This relationship has proven remarkable in its ability to describe the
geometry of networks ranging from the laboratory to continental scale with an
average value for the scaling exponent of b ≈ 0.6, though the value does exhibit
spatial variation (Cheraghi, Rinaldo, Sander, Perona, and Barry (2018)). Hack
also noted a power-law relationship between drainage area and channel slope,
which later work by J.J. Flint (Flint (1974)) would suggest is indicative of an
equilibrium tendency. In modern studies this relationship - termed Flints Law -
is often leveraged for the sake of landscape process partitioning, and the evaluation
33
Figure 1. Figure from Hack et al. (1957) showing the power-law
relationship between stream length and upstream drainage area, now
known to as Hack’s Law, in several catchments in Maryland and
Virginia.
34
of the degree of equilibrium of fluvial networks. In particular, considerable effort
has gone into identifying the locations and causes of scaling breaks in binned plots
of slope-area data, which have been correlated both to process transitions, and the
transmission of signals of dynamic adjustment in landscapes (Montgomery and
Foufoula-Georgiou (1993); Stock and Dietrich (2003)). In this work we quantify
the evolution of topographic curvature in such drainage area space, seeking better
resolution on the location of process transitions and to gain insight in to the
organizational structure of fluvial landscapes.
2.2 Connecting topographic curvature to classical surface theory
The term “curvature” is not associated with a single geometric entity, but
is rather an entire class of mathematical operations that describe the degree to
which a surface (or more generally, a manifold) deviates from being planar. It is
an invariant quality of a surface, independent of chosen coordinates and unchanged
under linear transformations. Measures of curvature are often classified as either
extrinsic or intrinsic, where extrinsic measures depend on an externally defined
reference frame (for example, a higher dimensional space in which the surface is
embedded) while intrinsic measures do not. A common conceptualization of this
idea states that intrinsic curvature must be built on information that would be
readily available to tiny inhabitants of a surface who have no external perspective,
or knowledge of any additional dimensions within which the surface is embedded
(Needham (2021); Struik (1950)).
Drawing connections between extrinsic and intrinsic curvature has occupied
some of the greatest scientific minds of the last several centuries (Needham
(2021); Penrose (2004)) and has led to significant scientific discoveries in that
time. Perhaps the most striking example is Albert Einsteins 1915 hypothesis that
35
gravity reflects the intrinsic curvature of 4-dimensional space-time bent by local
concentrations of mass/energy (Einstein (1982)). This realization revolutionized
cosmology, heralded the modern era of satellite technology, and showed that
while our 3-d human perspective struggles to interpret the 4-d space we inhabit
much is learned from considering the intrinsic properties of systems embedded
in higher dimensional spaces. For the purpose of geomorphology, a combined
extrinsic/intrinsic perspective of surfaces viewed as 2-d manifolds embedded in a
3-d space allows us to derive a rich description of topographic form. We will show
that standard surface characterization from 2-d map projections lose significant
information compared to a more rigorous curvature analysis.
There are a variety of methods for curvature calculation presented in the
geomorphology literature. For example Shary (1995) used a similar conceptual
framework to that of this chapter, but with a more explicit connection to the
gravity field, to derive 12 curvature classes that they argued could be used to
classify the landscape and possibly gain insight into the degree of landscape
equilibrium. Passalacqua, Do Trung, Foufoula-Georgiou, Sapiro, and Dietrich
(2010) used curvature of constant elevation lines in combination with drainage
area thresholding to construct a method for channel extraction. Minár, Evans, and
Jenčo (2020) presented an extensive list of possible land surface curvature metrics
and proposed possible links to topographic equilibrium. Finally, Schmidt, Evans,
and Brinkmann (2003) derived the same invariant metrics presented in this chapter
for 2-d polynomial fits of topography. The goal of our work is not to present new
definitions of curvature, but rather present a system of curvature metrics that are
useful in landscape evolution contexts. In landscape evolution studies “curvature”
usually refers to an approximation given by the Laplacian operator, either applied
36
directly to DEMs or 2-d polynomial fits of elevation data (M. D. Hurst, Mudd,
Walcott, Attal, and Yoo (2012); Perron et al. (2009); Roering et al. (1999); Schmidt
et al. (2003)). Such measures, which come directly from mechanistic derivations
assuming diffusive behavior, are extrinsically defined and thus depend on the choice
of coordinate system, and orientation within that reference frame. This means that
such methods admit systematic error that scales with slope and deformation of
the surface. In 1-dimension slope dependent error of a second partial derivative
curvature (d2z/dx2), shown in Figure 2 can be easily quantified through comparison
to the standard calculus definition of curvature (k)
d2z/dx2
k = (2.4)
[(1 + dz/dx)2]3/2
with the error given as (Bergbauer and Pollard (2003) and references therein)
k − d2z/dx2
error(%) = × 100 = 1− [1 + (dz/dx)2]3/2 × 100. (2.5)
k
In 2-dimensions error is much more complicated and also depends on the curvature
and orientation of the surface relative to the defined coordinates Bergbauer and
Pollard (2003). Thus curvature methods that depend on externally defined global
coordinate systems are limited in their utility on complex surfaces motivating our
more formal definition of land surface curvatures.
Real Earth surface topography is nearly always approximated by a gridded
discrete set of elevation points called a Digital Elevation Model (DEM). To
understand topography in the terms of classical surface theory it is necessary
to make a conceptual shift from viewing a DEM as a flat map of topography
to viewing it as a set of irregularly spaced data points sampling a 2-d surface
embedded in a 3-dimensional space. This is an important distinction because while
map projections of DEMs assume uniform distances and angles between grid points
37
Figure 2. Error in curvature defined as second partial derivative as a
function of slope along a 1-d curve.
38
it is mathematically impossible to map a surface onto a lower dimensional manifold
without introducing distortion, a fact proven by mathematician Leonard Euler
in 1775 (Needham (2021)). It follows that surface metrics calculated on gridded
elevation data admit this distortion and introduce an additional source of error that
is difficult to quantify.
The effects of the projection process can be seen in Figure 3, which
compares the distances and angles of a map projection (Figure 3.A) to those of the
same grid lines overlain on a 3-d model of the surface (Figure 3.B) . It can be seen
that the E-W and N-S grid lines of the map are not in fact perpendicular on the
surface nor are grid cells composed of uniform dimensions. In other words, if one
were to walk from point p along vectors defined by the grid lines they would not
cover the distance dŝ to arrive at point q̂ as implied by the map, but would rather
travel the distance ds arriving at point q.
Defining relationships between surfaces and their map representation was a
goal of Carl Frederick Gauss, who roughly fifty years after Euler’s treatment of the
problem would publish a revolutionary approach for extracting a combination of
intrinsic and extrinsic surface information that could be applied within Cartesian
coordinate systems (Gauss (1902)). One reason this is challenging is that while
Euclidean spaces exhibit global parallelism (the quality that parallel lines remain
parallel for ever; Figure 3.A), curved surfaces do not, and so cannot be defined by
any single set of global coordinates. Gauss got around this by everywhere defining a
local coordinate system comprised of arbitrary intersecting curves (the dashed lines
in Figure 3.B) and using them to define variation in the surface about a point.
Because coordinate curves can themselves be defined in terms of a Cartesian
reference frame, Gauss was able to relate intrinsic distortion of the surface to
39
ALorem ipsum dolor sit B
q
q q
s qd ds
p dx
= E1
/2dx
p du 
Figure 3. Difference between distances and angles measured on a map
projection versus on the surface. A. Map projection of DEM including map
grid defined by E-W and N-S lines with grid spacing dx. The red line corresponds
to the rectangular outline of panel B. B. DEM viewed as a 2-d manifold embedded
in a 3-d space. Dashed lines how a locally defined u-v coordinate system that
follows x and y curves on the map projection, but which are not orthogonal or
of equal length due to surface distortion. E and G are coefficients of the first
fundamental form, and ds is the displacement vector that results moving one grid
space along each of these coordinate vectors
40
dy
dv = G1/2dy
extrinsically defined curvatures through two quadratic partial differential equations
known as the first and second fundamental forms. Through this approach Gauss
was able to derive a measure of curvature that is independent of the shape of the
surface in space, an observation that he highlighted in his Theorema Egregium
(Remarkable Theorem), and which showed the ability of his method to extract
truly intrinsic surface qualities. This discovery inspired an effort to derive a fully
general description of surfaces, which led to the proof of the fundamental theorem
of surface theory in 1867 by Pierre Ossian Bonnet (Cogliati and R (2022)), which
states that all qualities of a surface, with the exception of its orientation in space,
can be derived using the 6 coefficients of the first and second fundamental forms.
During the same time period Bernhard Riemann, a student of Gauss, set
the stage for modern differential geometry by generalizing Gauss’ ideas to higher
dimensional spaces and showing that characterizing curvature of an n-dimensional
manifold requires 1 n2(n2 − 1) distinct components. Thus while a single measure
12
of curvature can be sufficient to characterize a 2-d map projection of topography
(or curvature in a temperature field as in the heat equation), fully characterizing
topographic curvature in 3-d space requires six distinct components, which are
the coefficients of the first and second fundamental forms. The field of differential
geometry was again revolutionized by the development of tensor calculus in 1901 by
Gregorio Ricci and Tullio Levi-Civita (Needham (2021)). Their methods combined
with the theory of Riemannian manifolds has allowed the field to evolve well
beyond the work of Gauss, however we find the original theory is sufficient to shed
significant light on the structure of topography. We therefore root our analysis in
the fundamental forms and a relatively simple classical definition of topographic
curvature.
41
2.3 Oregon Coast Range
We test our method of geometric classification in the Oregon Coast Range,
USA, in the Umpqua River Basin near Reedsport, Oregon. The Oregon Coast
Range is one of the most heavily studied terrains in the field of geomorphology
and is generally considered to be an archetypal example of a landscape where
uplift and erosion are balanced over long timescales (a common definition of
‘equilibrium’ in landscape evolution studies)(Reneau and Dietrich (1991)). The
Coast Range is mostly composed of a suite of accreted Eocene turbidites known
as the Tyee Formation that is fairly uniform and gently sloped, mostly at angles
less than 20 degrees, and so there is thought to be little variation in the lithologic
control on erosion (Roering, Kirchner, and Dietrich (2005)). It has a temperate
climate with precipitation rates that are fairly uniform throughout. This spatially
uniform landscape forcing results in pervasive ridge-valley topography, thought to
be indicative of steady fluvial erosion processes .
2.4 Calculation of curvature
In this section we derive the coefficients of Gauss’ two fundamental forms,
and use them to derive curvature metrics used in our topographic analysis. While
there are many different methods for deriving such measures curvature we base our
derivation on that of Struik (1950), which is the same derivation used in similar
work presented in Bergbauer and Pollard (2003) and Mynatt, Bergbauer, and
Pollard (2007). Everything in the following sections on deriving curvature can be
found in chapter 2 of Struik (1950) unless otherwise noted.
2.4.1 Deriving the first and second fundamental forms. Surface
curvature at a given point can be defined by considering the 1-d curvature of paths
over the surface. Infinitesimal displacements along such a path trace out small
42
sections of curve (ds) that lie completely within the so-called osculating plane,
which also contains the unit tangent (t) and normal (N) vectors of the curve. We
first define a position on the surface in 3-d Euclidean space as
r = r(u, v) = r1e1 + r2e2 + r3e3, (2.6)
where the parameters u and v represent two intersecting curves on the surface that
define a local coordinate system, and e1, e2, and e3 are basis vectors representing a
Cartesian reference frame. The choice of u and v is completely arbitrary except for
the requirement
ru × rv ̸= 0 (2.7)
which ensures that a normal vector can be defined.
A tangent vector (t) and unit normal vector (N) can then be defined in
terms of the u, v curves as
t = dr/ds (2.8)
and
ru × rv
N = . (2.9)
|ru × rv|
respectively, where dr = rudu + rvdv is the change in position resultant of a small
displacement on the surface and the magnitude of the displacement is given by the
Pythagorean Theorem as
I = ds2 = dr · dr = Edu2 + 2Fdudv +Gdv2. (2.10)
Equation 2.10, known as the first fundamental form (I) or surface metric formula,
is a measure of distances on the surface, where E = ru · ru, F = ru · rv, and
G = rv · rv (known as the metric coefficients) quantify the proportionality of real
distances to distances on a 2-d map representation (Needham (2021)).
43
The curvature, defined as the change in the tangent t along ds, can be
decomposed into two vector components as
dt
= κN+ κgg (2.11)
ds
where κ is the normal component of curvature, found by projecting the full
curvature onto the unit normal vector (N), and κg is the geodesic curvature which
measures meandering of the path along the surface and points in a direction
orthogonal to the path within the tangent plane. Since it is always possible to
define ds along a geodesic curve such that κg = 0 Equation 2.11 reduces to
dt dt
= κN = ( ·N)N. (2.12)
ds ds
Leveraging the fact that N · t = 0 we can write
dt
κ = ·N = − dNt · (2.13)
ds ds
where dN = Nudu+Nvdv. Combining Equations 2.8 and 2.13 yields the expression
for normal curvature
dr · dN edu
2 + 2fdudv + gdv2 II
κ = = = . (2.14)
ds ds Edu2 + 2Fdudv +Gdv2 I
In the above equation the expression in the numerator
II = dr · dN = edu2 + 2fdudv + gdv2 (2.15)
is the second fundamental form (II) with coefficients e = ruu ·N, f = ruv ·N, and
g = rvv · N that are measures of directional curvatures that are in turn scaled by
the coefficients of the first fundamental form (I) to define the directional change in
orientation of the surface.
2.4.2 Finding the principal curvatures. In the previous section a
general expression for curvature was defined relative to a single arbitrary path over
the surface. The existence of infinite paths through any given point means there are
44
infinite potential values for curvature, however there exist orthogonal paths along
which the maximum and minimum surface curvatures, known as are found. This
observation can be attributed to Leonard Euler, who showed that as the osculating
plane rotates about an axis defined by the surface normal vector N the curvature
varies systematically as
κ(θ) = k 21 cos θ + k2 sin
2 θ (2.16)
where k1 and k2 are the curvature extrema known as the principal curvatures. The
above equation, known as Euler’s Theorem, shows that if the principal curvatures
are calculated, they can then be used to understand curvature along any path.
To find the directions of the two principal curvatures we first define λ =
dv/du allowing us to rewrite Equation 2.14 in terms of a single variable as
e+ 2fλ+ gλ2
κ = . (2.17)
E + 2Fλ+Gλ2
Since the principal curvatures represent extreme values they must correspond to
directions where dκ/dλ = 0, and so we differentiate Equation 2.17 with respect to λ
which gives
dκ
= (E + 2Fλ+Gλ2)(f + gλ)− (e+ 2fλ+ gλ2)(F +Gλ) = 0. (2.18)
dλ
This can be put in the quadratic form
(Fg −Gf)λ2 + (Eg −Ge)λ+ (Ef − Fe) = 0 (2.19)
with roots that correspond to the principal curvature directions. The magnitudes
of the principal curvatures can be found through a similar approach. Given the
condition dκ/dλ = 0 equations 2.17 and 2.18 can be combined to give a more
simple expression for the curvature
f + gλ
κ = . (2.20)
F +Gλ
45
Recognizing that
E + 2Fλ+Gλ2 = (E + Fλ) + λ(F +Gλ) (2.21)
and
e+ 2fλ+ gλ2 = (e+ fλ) + λ(f + gλ) (2.22)
equation 2.18 can be rearranged to show
f + gλ e+ fλ
= (2.23)
F +Gλ E + Fλ
and so it follows that
II f + gλ e+ fλ
κ = = = . (2.24)
I F +Gλ E + Fλ
If we rewrite the two expressions for curvature give above as
f − κF + λ(g − κG) = 0 (2.25)
and
e− κE + f − κF = 0 (2.26)
and multiply both equations through by du (remembering that λ = dv/du) then we
arrive at a system of equations in terms of our original u-v coordinate system
    e− κE f − κFdu 0=  . (2.27)
f − κF g − κG dv 0
This has a non-trivial solution if and only if the determinant of the coefficient
matrix is zero. Taking the determinate gives
(EG− F 2)κ2 − (gE − 2fF + eG)κ+ (eg − f 2) = 0, (2.28)
46
a quadratic equation in κ of which the roots are the principal curvatures of the
surface. We sort these principal curvatures by magnitude, classifying the higher
magnitude curvature as k1, and the lower magnitude curvature as k2.
It should be noted that in more modern literature the information contained
in the fundamental forms is often presented in terms of a geometric entity known as
the Shape Operator, which is the second order tensor given by taking the negative
gradient of a normal vector taken along an arbitrary tangent. This is written
mathematically as
S(t) = −∇tN (2.29)
where t is the tangent vector (O’Neill (2006)). In terms of our u-v coordinate
curves the first and second fundamental forms can be defined via this approach
as I(u,v) = u · v and II(u,v) = S(u) · v (Needham (2021)).
In tensor form, also more common in the modern literature, the first
fundamental form is described by the so-called metric tensor
E F

I =  , (2.30)
F G
the second fundamental form the surface curvature tensor
 e fII =  , (2.31)
f g
and the shape operator is given by
S = I−1II. (2.32)
The eigenvalues of the shape operator tensor are the principle curvatures and the
eigenvectors define the directions of the paths associated with these curvatures.
47
The Gaussian and Mean curvatures, which we will define in the following section,
are the determinate and trace of this matrix respectively.
An approach to surface curvature calculation for gridded DEM data that
utilizes the shape operator can be found in Mynatt et al. (2007) and Pearce, Jones,
Smith, McCaffrey, and Clegg (2006) which apply a similar approach to ours in the
context of structural geology. While the mathematical simplicity of this approach
is appealing, it is considerably more computationally expensive to calculate the
shape operator. This in practice limits the scale over which the methods can be
applied. For example when curvatures were calculated on a 1000× 1000 DEM using
a MATLAB code similar to that employed by Mynatt et al. (2007) the calculation
took 13 seconds, while our code that instead finds the roots of the quadratic
equations defined above gave identical results in 0.11 seconds, which makes the
tool much more practical for landscape scale analysis.
2.4.3 Defining the Mean and Gaussian Curvatures and Surface
Shape Classes. Once the principal curvatures are found using one of the
approaches outlined above they can be used to calculate two useful measures of
surface curvature, namely the Mean and Gaussian curvatures, which themselves
measure different characteristics of the surface. The mean curvature follows directly
from Euler’s Theorem, and is the value about which the curvature oscillates in the
periodic function defined by equation 2.16, and shown in figure 4.
We calculate the Mean curvature (KM) from the principal curvatures as
k1 + k2
KM = . (2.33)
2
It is an extrinsic measure of curvature, and so is an measure of the shape and
orientation of the surface in reference to our defined coordinate space, but would
not be directly observable to inhabitants of the surface. One reason the Mean
48
Figure 4. Curvature as a function of the angle between the osculating
plane and a reference direction. Here the angle θ = 0 is associated with the
first principal curvature (k1). k2 is the second principal curvature and KM is the
Mean curvature about which the value oscillates.
49
curvature is useful is that it actually does not require any knowledge of the
principal curvatures, but can be calculated as the average curvature of any two
orthogonal paths through a given point. We will leverage this in section 2.7 to
assess error in curvature calculated using the Laplacian operator, which should
theoretically give a value twice that of the Mean curvature at any given point, as
long as the coordinate paths are orthogonal and slope error (figure 2; equation 2.5)
is accounted for.
The Gaussian curvature is the remarkable measure of intrinsic curvature
that Gauss’ Theorema Egregium noted as a measure of surface deformation that is
invariant under isometries, which are spatial transformations that conserve distance
(Needham (2021)). This means its value does not actually depend on the shape of
the surface, but captures a more subtle quality, which is the degree of stretching or
bending required to deform a flat plane to fit a given surface. To a tiny inhabitant
of the surface the Gaussian curvature defines the distances between points, and the
rate of divergence or convergence of locally parallel geodesic paths, and is thus an
extremely fundamental property of surfaces. In terms of the principal curvatures it
is calculated simply as
KG = k1k2. (2.34)
We note that while we will consistently refer to this curvature as the Gaussian
curvature, in the literature you will see it called the intrinsic curvature, total
curvature, or sometimes just the curvature (Needham (2021)).
The Mean and Gaussian curvature together can be used to simply classify
the geometry about a point into the 8 distinct shape classes shown in Figure 5.
Since the Gaussian curvature is the product of the two principal curvatures it will
only be positive in instances where k1 and k2 have the same sign. This tells us
50
that when KG is positive a surface is either a dome or basin, though because the
Gaussian curvature is isometrically invariant it does not alone determine which.
If KG is negative then k1 and k2 necessarily have opposing signs and so the
surface is locally a saddle, though again with an orientation in space that cannot
be uniquely determined by the Gaussian curvature. If either k1 or k2 is equal
to zero then KG is also zero. This case itself comprises its own class of surfaces,
known as developable surfaces, which can be formed from a plane without altering
surface area and thus are intrinsically flat. Relatively developable surfaces can
be extracted from the landscape by assigning a zero value to curvatures under
a defined threshold (Mynatt et al. (2007)). Curvature thresholding to extract
developable forms could be a promising way of quantifying the contributions of
different landforms to overall curvature distributions, however in the interest of
minimizing variables we do not explore this quality of surfaces here.
The sign of the mean curvature tells us the orientation of a shape in space,
and so allows us to develop KG based classifications in a landscape reference frame.
KM is positive in two cases; when either both k1 and k2 are positive, or the higher
magnitude curvature is positive. This means that points in the landscape with
KM > 0 are locally either domes or antiformal saddles. Similarly if KM is negative
then the orientation of the surface must be dominantly concave up, and so the
surface is either a basin or synformal saddle.
In rare cases where k1 and k2 are equal and opposite the surface is a perfect
saddles. This too comprises its own special class of surfaces, known as minimal
surfaces, which have been largely tied to the optimization of surface area in
self-organized systems such as soap foams (Osserman (1969)). The existence of
approximate perfect saddles could be shown in the landscape by rounding the
51
Figure 5. Shape classes into which points on the surface can be sorted
based on the signs on the Mean (KM) and Gaussian (KG) curvatures. In
this analysis we focus on those classes that can be assigned based on raw curvature
values, which are synformal saddles, antiformal saddles, basins, and domes, and do
not include developable surfaces or perfect saddles.
52
principal curvatures to lower precision, which has the possibility of revealing
process scales for which minimal surface theory could be applied, however we do
not endeavour to explore this aspect of the data. Instead we focus on the 4 basic
shape classes that can be assigned based on raw curvature calculations.
2.5 Calculation of Curvature on Gridded Data Sets
In order to calculate the curvatures of topography it is necessary to first
select a DEM dataset, and then smooth the data in order to remove errors
introduced in the gridding process. Since DEMs are themselves an interpolation
of other forms of elevation data, artifacts are known to arise and are sensitive
to data density, the type of raw data used, and the interpolation method used in
gridding (Reuter, Hengl, Gessler, and Soille (2009)). Removing such contamination
is essential for our curvature calculations since they are computed on a grid of
neighboring cells, and are thus sensitive to the highest frequency content on a given
surface.
While it is well established by the Shannon-Nyquist sampling theorem that
discretely sampled data cannot resolve any structure with an effective wavelength
less than twice the sample spacing (Stewart and Podolski (1998)), gridding artifacts
stemming from data interpolation represent a less obvious constraint on the scale
of structures that can be resolved in a given DEM. In the past decade, DEMs
generated with LiDAR (Light Detection and Ranging) have become the preferred
data source for connecting morphometric analysis to landscape evolution models
due to their superior resolution and lack of obvious gridding artifacts relative to
other elevation datasets.
While these factors can rightfully be taken as evidence of the superior
quality of LiDAR over older data products, it is also important to recognize that
53
LiDAR DEMs are built through interpolation techniques that are themselves a
form of low-pass filtering, and which are known to introduce errors dependent on
slope, terrain roughness, and sample density, in the original point cloud (Bater and
Coops (2009); Bui and Glennie (2023); Reuter et al. (2009); S. Smith, Holland, and
Longley (2004)).
Within our study area there is 1-m resolution LiDAR data available
through the Oregon Department of Geology and Mineral Industries
(https://gis.dogami.oregon.gov/maps/lidarviewer/), however close examination
of this data reveals interpolation artifacts that we interpret as being indicative of
poor resolution in the source data. These artifacts can be seen as the triangular
structures in Figure 6.B, which shows the hillshade for an unfiltered swath of
LiDAR within the study area. While these artifacts, which are common in high
curvature areas and persist for hundreds of meters, can easily be smoothed over it
is not clear over what scale this data is actually representative of a topographic
surface. While it is beyond the scope of this work to quantify how this type of
error propagates into the calculation of landscape form, it does show that LiDAR
grid spacing is not always reflective of the useful resolution of the data (T. Smith,
Rheinwalt, and Bookhagen (2019)). Given the scale of these interpolation artifacts
we opt to perform our analysis on 8.1 m resolution DEM data freely available
through the National Map (https://apps.nationalmap.gov/downloader/).
2.5.1 Spectral Filtering. To filter the DEM we use the Discrete
Fourier Transform (DFT; also a contribution of Gauss), which decomposes
discretely sampled signals into sums of trigonometric functions. This allows for
evaluation of spatial signals in terms of their wavenumber content, and in the
context of smoothing, the ability to eliminate wavenumbers higher than a defined
54
Figure 6. Interpolation artifacts in 1-m LiDAR data. A. The location of
the LiDAR swath within the study area. B. LiDAR hillshape showing regions
dominated by interpolation artifacts over length scales of hundreds of meters.
cutoff leaving behind longer wavelength structures. In this way a low-pass filter
can be constructed that preserves the amplitude of signals at wavelengths longer
than that of gridding noise, however this approach requires careful processing
and interpretation, as calculation of the Fourier spectrum can easily generate
spurious artifacts at a range of scales. Fourier methods have already been applied
in geomorphology towards the identification of characteristic process scales
(Perron, Kirchner, and Dietrich (2008)), landform identification (Booth, Roering,
and Perron (2009)), and to assess topographic controls on the development
of volcanic topography (Richardson and Karlstrom (2019)). We build on this
literature to develop a low-pass filtering technique that is slightly more robust
against long-wavelength contamination, and show the utility of this method for
removing erroneous high wavenumber signals introduced in the process of gridding
topographic data (Reuter et al. (2009)).
55
2.5.1.1 Minimizing Long-Wavelength Contamination. One
challenge of applying Fourier analysis to DEMs, is that the method fits the
surface with a sum of periodic functions which collectively must go to zero at
the edges. Since topography at the scales of interest is not bounded by regions
of zero elevation it is required to taper the data on the margins such that this
boundary condition is met. Most geomorphology studies to date accomplish this by
convolving the DEM grid with a 2-d raised cosine (aka Hanning window), such that
the resulting topography is equal to its actual value only in center of the grid, and
is elsewhere subdued by an amount that scales with distance towards the margins.
This approach not only diminishes the spectral power of relevant landscape
features, but introduces artificial long-wavelength signals that are difficult to
differentiate from real structures in the resulting power spectrum. These effects
of windowing are generally not acknowledged in the geomorphology literature,
however this issue is well described in the signal processing literature (ex. Harris
(1978)), as well as the geophysical literature (where spectral filtering is common for
computing gravity/topography admittance). For example Mcnutt (2005) showed
that contamination can be minimized by first reflecting the topographic grid in
each direction, then tapering the data only in the reflected portions, which are
outside the limits of the original DEM. This avoids the introduction of distortion
within the study area through the windowing process, and only introduces artificial
signals with significant power that are at or above the scale of the original DEM
extent. After removing high frequency content and performing a reverse Fourier
transform the original data can then be extracted yielding a more accurate low-pass
filtered representation of the surface.
56
We adopt this mirroring approach and apply a Tukey (Harris (1978))
window, which consists of a boxcar function convolved with a cosine taper along
the margins. We apply the default Tukey window in the window2 function in
Matlab. We find this low-pass filtering approach more compelling than the cosine
taper method as distortion is not introduced in the pre-processing step. Figure
7 compares our pre-processing approach with the default window used in the
2dSpecTools toolbox for Matlab Perron et al. (2008), and also shows the effects
of applying the window used by (Perron et al. (2008)) to a mirrored version of
the data. It can be seen that both methods employing the hanning window cause
significant distortion of the data that can be expected to modify the spectral
content over relevant length scales.
2.5.1.2 Calculation of amplitude spectrum. Once the data has
been mirrored the Discrete Fo∑urier∑Transform (DFT) is calculated asNx−1Ny−1 k m
Z(k , k ) = z(m∆x, n∆y)e−2πi(
kxm+ y )
Nx Nyx y (2.35)
m=0 n=0
where Nx and Ny are the number of grid cells in each direction, m and n
are array indices, ∆x and ∆y are the grid spacings in each direction, and kx and
ky are the wavenumbers in the respective x and y directions (Perron et al. (2008)).
Each value in the output array given by the above equation is associated with a
frequency in each the x and y directions with magnitudes
kx ky
fx = , fy = . (2.36)
Nx∆x Ny∆y
These frequencies can then be used to√define a radial frequency as
fr = (f
2
x + f
2
y ). (2.37)
57
Elevation (m) Elevation (m) Elevation (m)
-200 -100 0 100 200 -200 -100 0 100 200 -200 -100 0 100 200
A B C
Elevation (m) 200
-200 -100 0 100 200 E
150
Detrended DEM
D Hanning window100 Mirrored with tukey window
Mirrored with hanning window
50
0
-50
-100
-150
4.2 4.21 4.22 4.23 4.24 4.25 4.26 4.27
x 105
Easting (m)
Figure 7. Comparison of windowing techniques used in preprocessing
DEMs. A. Raw DEM. B. Mirrored data with Tukey window. The red box is the
orignal DEM extent. C. Mirrored with hanning window. D. Hanning window
applied to DEM without mirroring. E. Profiles of topography for each of the
preprocessing methods compared to the original DEM.
58
Elevation (m)
We calculated the DFT periodogram as
1
Pf (kx, ky) = |Z(kx, ky)|2 (2.38)
N2 2xNy
which can be plotted against radial frequency to provide an estimate of the
landscape power spectrum as seen in figure 8.A.
2.5.1.3 Removal of High-Frequency Noise. The primary goal of
the filtering process is to remove high wavenumber artifacts of the gridding process.
Since DEMs are themselves interpolations of other forms of elevation data they are
prone to gridding artifacts that vary between datasets (Reuter et al. (2009)). In
addition to the gridding artifacts themselves it is expected that such discontinuities
in the data give rise to Gibbs phenomena, which is high wavenumber “ringing”
around discontinuities in the data that comes from trying to fit a discontinuity with
layered high-wavenumber signals. This is especially of concern for the methods used
in this study, as curvature is calculated on a grid of neighboring DEM cells, and so
is itself a measure of the highest frequency content in a given DEM.
In the particular dataset used in this study the grid artifacts, and resultant
Gibbs phenomena can easily be seen in the amplitude spectrum. In figure 8.A the
farthest right bump in the power spectrum is the result of gridding artifacts as can
be seen in figure 8.E, in which the raw topography (panel B) is high pass filtered
to only admit wavelengths under 30 m. The stripes that comprise most of the
signal at these wavelengths is characteristic of gridding artifacts (Bater and Coops
(2009)). The bump at slightly lower frequencies we attribute to Gibbs phenomena
that stem from the gridding artifacts. Features associated with this peak can be
seen in 8.D which shows the result of bandpass filtering the raw DEM to isolate
wavelengths between 30 m and 50 m. The interpretation of these features as Gibbs
phenomena comes from the observation that structures at these scales follow the
59
Figure 8. Identification of gridding artifacts in the 8.1 m DEM data used
in this study. A. Mean squared amplitude vs. radial frequency of topography.
Colored boxes highlight frequency ranges shown in panels C, D, and E. B. Raw
DEM. C.DEM bandpass filtered to admit60wavelengths between 30 and 50 m
showing the gibbs phenomena thought to produce the spike in the power spectrum
at these wavelengths. E. DEM high-pass filtered to admit wavelengths less than 30
m. Straight line are artifacts from the gridding process.
orientation of grid lines, but truly defining the source of these signals requires
further analysis. This could be accomplished through testing of different windowing
methods, which should conserve real structures, but alter the expression of such
ringing signals (Harris (1978)). We conclude that this data is not suitable for the
calculation of curvatures at wavelengths of less than 50 m and so take this as the
minimum low-pass filter cut-off for accurate curvature calculation, however we will
do most of our analysis on DEMs filtered to 200 m, and so do not further explore
the sources of these high frequency signals. The DEM filtered to a lowpass filter
cut-off of 50 m to eliminate these signals can be seen in figure 8.C.
2.5.2 Scale Selection. Once the landscape has been filtered all of the
curvature metrics outlined in sections 2.4.2 and 2.4.3 are calculated for every pixel.
The values for each of the curvature metrics are binned by drainage area in order
to couch our curvature analysis in a context consistent with that of the literature
for identifying process transitions. We first look for dominant trend in the evolution
of each metric, and iteratively smooth the landscape to greater wavelengths to get
a sense of how our measures of curvature vary with scale. To eliminate the high
wavenumber gridding artifacts identified in the previous section we first smooth the
landscape to 50 m, then increase the low-pass filter cutoff by increments of 50 m up
to 600 m.
The results can be seen in Figure 9, which shows the land surface gradient,
Gaussian and Mean curvatures, as well as the two principal curvatures and area
proportionality factor (α) binned by drainage area. The area proportionality factor
(α) measures the proportionality of surface area on the landscape to its map-
view representation, and can be calculated in terms of the coefficients of the first
61
√
fundamental form as α = EG− F 2. It is thus a measure of distortion inherent to
the map projection.
We see that general trends in all of these curvature metrics are similar across
this range of filter cutoffs, making the selection of a particular scale for further
analysis somewhat arbitrary. However we note that at cutoffs greater than ∼ 200
m peaks in the curves start to get pulled to longer wavelengths suggesting a shift
in the structures being resolved at those scales. This is most easily seen in the
main peaks of the Gaussian curvature in Figure 9.B. We take this as an upper
end of the range of length scales for primary landscape features, and so do our
focused analysis on a DEM filtered to 200 m. This filter scale allows us to identify
certain landscape features, however inhibits our ability to see the effects of finer
scale processes. For example signatures of processes, such as tree throw, cannot be
resolved and would require analysis of finer datasets and/or less aggressive filters.
2.6 Alignment of Slope and Principal Curvature Vectors
Upon calculating the principal curvatures, they are sorted such that the first
principal curvature (k1) is defined as the curvature with the greatest magnitude
regardless of sign as explained in section 2.4.2. This is not a unique choice (Figure
4 shows the first principal curvature as the more positive of the curvatures) but we
find this way of defining the principal curvatures is slightly easier to interpret in
terms of landform development.
Since we also are able to define the orientation of the principal curvature
vectors via Equation 2.19 we can then look at how each principal curvature aligns
with slope throughout the landscape. It is common in geomorphology studies
not to define the method of slope calculation, however many studies use either
the gradient operator (∇z(x, y) = ∂z + ∂z ) (Roering et al. (1999); Scherler
∂x ∂y
62
Figure 9. Surface geometry metrics binned by upstream drainage area for
a range of low-pass filter cut-offs between 50 m and 500 m calculated on
50 m intervals. A. Gradient. B. Gaussian curvature (KG). C. Mean Curvature
(KM). D. First principal curvature (k1). E. Second principal curvature (k2). F.
Surface area proportionality factor (α)
63
and Schwanghart (2019)), or an 8-point neighborhood gradient which calculates
the slope between a pixel and its lowest immediate neighbour (Montgomery
and Foufoula-Georgiou (1993); Schwanghart and Scherler (2014)). Both these
methods can only resolve slope contributions in directions defined by the DEM
grid structure and so do not have sufficient aspect resolution for correlation with
principal curvature vectors. For this reason we instead use the slope of the tangent
plane at each pixel defined via the unit normal vector given by (Equation 2.9). If N
is written in vector form as N = N1e1 + N2e2 + N3e3 then the magnitude of the
slope is given by (π )
ST = tan − sin−1(N3) (2.39)
2
and the orientation of the normal vector on the projected map can be written in
unit vector form as
N√1e1 +N2e2ST = . (2.40)
N21 +N
2
2
The degree of alignment between a principal curvature vector and this measure of
slope can be found using the scalar product where ST · ki = 1 if the slope and
curvature are perfectly aligned and ST · ki = 0 when they are orthogonal.
Viewing this alignment metric in map view, as seen in Figure 10, reveals
intuitive trends that can help inform the location of process transitions. The
second principal curvature (k2) is the least curved path over the surface and so
will always track the long axis of a given structure, which will only align with
the tangent slope in areas of symmetric convergence and divergence in directions
orthogonal to k2. Indeed Figure 10.D shows maximum alignment of k2 and ST
along ridgelines and channels, and in the hillslopes and colluvial hollows that mark
the termination of those structures. The first principal curvature (k1) is aligned
64
with slope in transition regions that are neither especially convergent or divergent.
Sharp transitions in the alignment of the principal curvatures thus reflect transition
areas.
We can see the same alignment trend in area space, though with less detail.
Figure 11 shows both the average magnitude of each principal curvature, and
our alignment metric binned as a function of drainage area. The sharp increase
in alignment between ST and k2 at drainage areas slightly higher than 10
3 m2
follows the transition of both principal curvatures to being dominantly negative,
itself a signature of flow convergence that we will show is associated with a major
shift in landscape concavity. We will suggest that this increase in alignment is
associated with colluvial hollows at the head of the drainage network, and thus is a
signature of convergent processes that mark the onset of hydrologic channelization
in landscapes.
Since we rely on a non-standard slope metric for this analysis it is important
to compare this metric to others commonly used in the field. Figure 36 shows
comparison between ST , the gradient operator (∇z), the 8-point neighborhood
gradient calculated using the gradient8 function in TopoToolbox, and the gradient
calculated in 1-dimension along all channels above an area threshold of (104) m2,
which are extracted using the STREAMobj object class in TopoToolbox. We
note that while ST and the 8-point neighborhood gradient are well aligned in the
steepest parts on the landscape, the 8-point neighborhood gradient predicts much
lower slopes at high drainage areas than the other metrics. This is relevant as it is
the metric often used to define process transitions in slope area space (Montgomery
and Foufoula-Georgiou (1993)). The gradient operator predicts systematically
higher slopes than the other metrics, which is important to note as it is often used
65
Figure 10. Alignment between principal curvatures and slope of the
tangent plane measured via the surface normal N. A. DEM low-pass
filtered to 200 m. B. Slope of tangent plane. C. Scalar product of the slope
vector and first principal curvature vector (k1).D.Scalar product of the slope vector
second principal curvature vector (k2).
66
Figure 11. Principal curvatures binned by upstream drainage area for
DEM filtered to 200 m. Left axis measures mean curvature values shown by
solid lines. Right axis measures the scalar product between the principal curvatures
and slope shown by the dotted lines.
in hillslope models (Roering et al. (2007)). The tangent slope is well aligned with
the slope of channels in the upper reaches of the network, however these metrics
also diverge at drainage areas that seem to correspond with a transition zone
already identified in the literature as connected to debris flow processes (Stock
and Dietrich (2003)). We have no reason to think any one of these slope metrics
superior to the others, however do note that each way of calculating slope extracts
different information from the landscape and will lead to different conclusions if
used to define the area thresholds of process transitions.
2.7 Connections to the Laplacian
The distribution of both the Mean (KM) and Gaussian (KG) curvatures
are shown in map-view (Figure 13) and in area space (Figure 14) giving different,
but complementary, views of how these curvature metrics relate to structures in
the landscape. As outlined in section 2.4.3, the general shape of the surface is
characterized by the Gaussian curvature, which takes on positive values when the
67
Figure 12. Comparison of different slope metrics binned by drainage area.
68
sign of the two principal curvatures are the same and negative values when they
are not. Thus basins and domes are marked by positive Gaussian curvature and
negative values indicate saddle structures, where the orientation of these shape
classes is recorded in the sign of the mean curvature. In Figure 13.B we see that
the landscape is composed of a somewhat tessellated pattern of oscillating Gaussian
curvatures while Figure 13.D shows that the implied alternating domal/saddle
structures have orientation that is dictated by their location either in the divergent
(ridge/hilltop) or convergent (valley network) portion of the landscape. The area
space perspective shown in Figure 14 gives a more blurred representation, and
suggests that the landscape dominantly consists of sources and sinks represented
by domes and basins with the exception of a transitional area around where the
mean curvature switches from positive to negative. We will examine this area space
trend in greater detail in section 2.9, and show that the curvature evolution shown
in Figure 14 is useful in defining the evolution of process regimes in the Oregon
Coast Range.
The Mean curvature also provides a way of quantifying complicated
error inherent to common map view Laplacian methods. While in this study
KM is calculated at point using the two principal curvatures, Euler’s Theorem
(Equation 2.16) actually shows that it can be found using the curvatures of any two
orthogonal paths. This means that if one is able to apply the Laplacian operator
along coordinate vectors that are orthogonal on the surface and account for slope
error, then the resulting measure of curvature relates to KM via the relation
z + z ∇2xx yy z
KM = = . (2.41)
2 2
69
Figure 13. Map-view of surface curvature metrics. A. Topography. B.
Gaussian Curvature (KG). C. Difference between mean curvature calculated using
the principal curvatures and Mean curvature calculated using the components of
the laplacian operator. D. Mean curvature (KM).
70
Figure 14. Surface curvatures binned by upstream drainage area.
The magnitude of the difference between the mean curvature and that given by the
calculation of directional derivatives along the map grid can be then written as
∇2z
error = |KM − |. (2.42)
2
Figure 13.C shows the distribution of error in map-view. We see that
maximum error in the Laplacian calculation occurs along ridge and channel
networks where the magnitudes of both the Gaussian and Mean curvatures are
highest. This is not surprising since these are the regions in the landscape with the
most complexity, and thus the most distortion introduced through projection onto
the map grid. This is also reflected in the area-space representation given in Figure
15.A, which shows that agreement between the Mean curvature calculated using
the principal curvatures and that calculated using the Laplacian are highest in
areas where the magnitude of both curvatures are lowest. Perhaps non-intuitively
slope is shown to have a lesser effect on the error in Laplacian curvature (Figure
15.C) and in fact error is minimized at higher slope. This is contrary to the general
71
expectation (Bergbauer and Pollard (2003)), however makes sense given the strong
control exerted by surface curvature itself on error in the Laplacian (Figure 15.B).
We note that this analysis only pertains to the simplest possible
implementation of the Laplacian operator calculated on our filtered DEM surface,
and thus should not be viewed as a measure of error in curvatures calculated in
geomorphology studies that use the Laplacian as most studies implement more
sophisticated pre-processing steps (M. D. Hurst et al. (2012)), and often calculate
curvature along 1-d profiles (Roering et al. (2007); Struble and Roering (2021)).
Rather this is a simplistic measure of the relative contributions of slope and
map distortion to general errors in 2-dimensional partial derivative curvature
calculations. It should also be noted that, while the magnitude of the curvature
is very different, the general trend in area-space is the same, and so it is likely
that general inferences of connections between curvature and erosion rate in the
literature (Roering et al. (2007); Struble and Roering (2021)) are not changed by
this observation.
2.8 Land Surface Shape-Class Distribution
As outlined in section 2.4.3, the Mean and Gaussian curvatures can together
be used to sort the landscape into four distinct shape classes. These amount to
a way of partitioning the landscape that provides a more intuitive framework for
discussing geomorphic structures than curvature metrics themselves. We look at
the distribution of these shape classes in the landscape by generating probability
density functions of each shape class as a function of contributing drainage area.
The initial PDFs are shown in Figure 16.B, which shows that probabilities are
strongly weighted by the distribution of drainage area (shown in Figure 16.A) and
72
Figure 15. Comparison between mean curvature calculated using the
principal curvatures (KM), and mean curvature calculated using the
components of the laplacian (KMLP). A. Calculated mean curvatures binned
by upstream drainage area. B. Difference between KM and KMLP binned by KM .
C. Difference between KM and KMLP binned by slope
.
73
are thus not a good reflection of shape class probabilities existing in a particular
area-dependant process regime.
For this reason we calculate the conditional probability of each shape class
at a given drainage area. We first define the probabilities of each shape class as a
function of drainage area (Figure 16.B), which can be thought of as the probability
of both the shape class and drainage area occurring simultaneously. More formally
this is the intersection of the sets containing shape classes and drainage area values
which is written P (C ∩ A). We can also define the probability of each drainage
area value (P (A); Figure 16.A), which allows us to calculate the conditional
probability of a shape class at a given drainage area (P (C|A); (Figure 16.C)) using
the probability axiom (Dekking (2010))
P (C ∩ A)
P (C ∩ A) = P (C|A)P (A) → P (C|A) = . (2.43)
P (A)
We find this version of the shape class metric the most useful for interrogating the
connection between surface geometry and surface processes, and explicitly connect
it to the other quantitative metrics of topographic form in section 2.9.
2.9 Geometric View of the Oregon Coast Range
Having defined in the previous sections several quantitative measures
of land surface geometry we now look holistically at these metrics as a suite
of complementary data that inform how the geometry of topographic surfaces
varies as a function of catchment area in fluvially dissected landscapes. We first
examine the study area as a whole, then decompose the landscape into regions
defined by transition points in the invariant measures of curvature, which we argue
are indicative of process transitions that define how signals of fluvial dissection
propagate through connected catchments.
74
Figure 16. Shape class distributions as a function of drainage area. A.
Probability density function of upstream drainage areas within the study area. B.
Probability density functions of shape classes as a function of drainage area. C.
Conditional PDFs of shape classes at a given drainage area.
75
Figure 17 shows a compilation of the geometric data and will serve as a
template for the decomposition of the landscape in the following sections. Panels
A -E are the area-space representations of the geometry metrics presented in the
previous section and do not contain any new information. Panel F shows the map-
view distribution of shape classes, while Panel H present this same information as a
pie chart. We see that at this spatial scale the landscape is evenly divided between
structures that are concave up and those that are concave down, which is reflective
of a symmetry in mean curvature that we will show in more detail in section X.
Panel G shows a probability density function of slopes and is helpful in drawing
connections between our analysis and the broader geomorphology literature, in
which process regions are often defined in terms of average slope (Roering et al.
(2007); Stock and Dietrich (2003)). We will first look at a landscape decomposition
based on the Gaussian curvature.
2.9.1 Gaussian curvature-based landscape partitioning. We
take a detailed look at how the Gaussian curvature changes with increasing
drainage area, and how it relates to process regimes already identified in the
literature.
2.9.1.1 Region 1: Drainage areas less that 3.3 × 102 m2 .
The first threshold we apply in our landscape analysis corresponds with the first
inflection point in the Gaussian curvature, which occurs at drainage areas of ∼ 330
m2. This region, highlighted in Figure 18, is comprised of the highest elevations
in the landscape and is associated with the network of peaks and ridges that
divide adjacent catchments and, in terms of plan-form area, make up 17% of
the landscape (Panel A). In map view this network is seen to consist of almost
uniformly spaced domes (peaks) connected via antiformal saddles that make up
76
Figure 17. Compilation of surface geometry data used in this study. A.
PDF of upstream drainage area. B. Slope metrics binned by upstream drainage
area. C. Surface curvatures binned by upstream drainage area. D. Principal
curvature values and slope alignment binned by upstream drainage area. E.
Conditional PDFs of shape classes as a function of upstream drainage area. F.
Map-view of shape class distribution. G. PDFs of slope metrics. H. Pie chart
distribution of shape classes.
77
the ridge network (Panel F). Together these two shape classes comprise over 99% of
this landscape region (Panels E and H).
Binned by area the two principal curvatures have fall on the same trend
in this region (Panel D), which we find noteworthy since in surface theory areas
where the two principal curvatures are the same represent special geometric
entities, known as an umbilic points, which are source or sink regions that are
deeply connected to the overall structure and topology of vector and/or scalar
fields (Needham (2021)). While is is difficult to see a direct connection between
these binned curvature values and the structure of individual points we find
this connection to geometric theory interesting, especially in a portion of the
landscape where the gradient vector field is characteristically divergent and serves
as a primary boundary in the organization of fluvial catchments. In this region
slope increases with drainage area, likely corresponding to convex rollovers onto
hillslopes.
This is a region of the landscape where curvature is often calculated on 1-d
profiles selected such that they are orthogonal to ridge-lines. The slope distribution
given in Panel G can be roughly correlated with Figure 2 to get a broad sense of
the potential uncertainty associated with such calculations.
2.9.1.2 Region 2: Drainage areas between 3.3 × 102 m2 and
1.3 × 103 m2. As drainage are increases beyond the 330 m2 threshold that
defined region 1 the Gaussian curvature transitions to being negative and remains
so up to drainage areas of ∼ 1300 m2. This region, which we correlate to side
slopes, is highlighted in Figure 19, and makes up 58% of our study area. Here
the landscape is strongly transitional, encompassing the only inflection point
in the mean curvature (Panel C), and a zone of divergence in the two principal
78
Figure 18. Surface geometry data for points in the landscape with
upstream drainage areas less than 3.3 × 102 m2. The red box in each of the
area plots highlights the region of interest. A. PDF of upstream drainage area. B.
Slope metrics binned by upstream drainage area. C. Surface curvatures binned by
upstream drainage area. D. Principal curvature values and slope alignment binned
by upstream drainage area. E. Conditional PDFs of shape classes as a function of
upstream drainage area. F. Map-view of shape class distribution. G. PDFs of slope
metrics. H. Pie chart distribution of shape classes.
79
curvatures (Panel D). Accordingly this is the most geometrically complex part of
the landscape. Broadly this region encompasses a transition from the dominance of
domes to that of basins, however in more detail the landscape is primarily saddles,
with an equal partitioning of synforms and antiforms (Panels F and H). This is
also the steepest section of the landscape and thus accounts for the majority of
landscape relief and course-scale roughness (Panels B and G). The distribution of
gradient values found in this region (Panel G) agrees with values found for Coast
Range hillslopes in other studies (Roering et al. (2007)).
2.9.1.3 Region 3: Drainage areas between 1.3 × 103 m2 and
4.5 × 103 m2. Following a second inflection point in area space the Gaussian
curvature sharply increases to a maximum at ∼ 4500 m2. The location of this local
maximum, which is seen in Figure 9 to be fairly consistent over a range of low-pass
filter cutoffs, defines the upper bound for the next region we examine, which is
shown in Figure 20. In this region the combination of dominantly positive Gaussian
curvature, negative mean curvature, and a decrease in slope with added drainage
area are all consistent with colluvial deposits at the base of hillslopes, where debris
from steeper upper reaches collects on the margins of channel networks as slope
decreases thus limiting the ability for material to carry momentum downhill. The
landscape is now strongly convergent, as shown by the strong dominance of concave
geometries in Panels E, F, and H. This region also aligns with a sharp increase
in alignment between the second principal curvature (k2) and the tangent slope,
with a corresponding decrease in alignment between slope and the first principal
curvature (k1).
This trend, seen in Panel C, is what we expect for the onset of
channelization and so we interpret this as a signal of first order channels at the
80
Figure 19. Surface geometry data for points in the landscape with
upstream drainage areas between 3.3 × 102 m2 and 1.3 × 103 m2. The red
box in each of the area plots highlights the region of interest. A. PDF of upstream
drainage area. B. Slope metrics binned by upstream drainage area. C. Surface
curvatures binned by upstream drainage area. D. Principal curvature values
and slope alignment binned by upstream drainage area. E. Conditional PDFs of
shape classes as a function of upstream drainage area. F. Map-view of shape class
distribution. G. PDFs of slope metrics. H. Pie chart distribution of shape classes.
81
head of the fluvial network. This region accounts for 18% of the study area, and so
has a similar contribution to the total landscape as the ridge-valley network defined
by our analysis as region 1.
2.9.1.4 Region 4: Drainage areas between 4.5×103 m2 and 9.8×
105 m2. The next defined region has a range of drainage areas spanning 3 orders of
magnitude over which geometric change happens gradually. Figure 21 shows that
the magnitude of positive Gaussian curvature decreases towards another inflection
point at A = 9.8 × 105 m2 which is indicative of a gradual decrease in the relative
proportion of basins to synformal saddles (Panel E). This trend intuitively implies
a growing influence of channels in defining landscape curvature and so we associate
this region with channels at the head of the fluvial network. This is consistent with
the mapview extent of this region (Panel F) which shows that these areas exist at
the bottoms of valleys in an area in, or directly adjacent to, main channels.
As drainage area increases the contribution to the total composition of the
landscape approaches zero with this entire region only making up 7% of the study
area (Panel A). This region contains on average the highest magnitude negative
mean curvatures, implying that it is extremely convergent, which taken with the
fact that it underlies a region defined by decreasing slope and thus the collection
of colluvium (region 3) can be taken to mean that it is a prime region for the rapid
accumulation of both hillslope debris and overland flow and thus is expected to be
a region where debris flow processes dominate. This is in broad agreement with
other literature focused on debris flow processes (Penserini, Roering, and Streig
(2017); Stock and Dietrich (2003); Struble et al. (2023)). This region also contains
the major divergence in the values of the different slope metrics. We include on
Panel G the average gradient calculated in 1-d along the longitudinal profiles of
82
Figure 20. Surface geometry data for points in the landscape with
upstream drainage areas between 1.3 × 103 m2 and 4.5 × 103 m2. The red
box in each of the area plots highlights the region of interest. A. PDF of upstream
drainage area. B. Slope metrics binned by upstream drainage area. C. Surface
curvatures binned by upstream drainage area. D. Principal curvature values
and slope alignment binned by upstream drainage area. E. Conditional PDFs of
shape classes as a function of upstream drainage area. F. Map-view of shape class
distribution. G. PDFs of slope metrics. H. Pie chart distribution of shape classes.
83
streams within the study region extracted using the STREAMobj object class from
TopoToolbox using a threshold drainage area for channelization of 4500 m 2.
2.9.1.5 Region 5: Drainage areas greater than 9.8 × 105 m2.
Our final point of differentiation based on the evolution of the Gaussian curvature
is located at its final transition from mostly positive to mostly negative, which
occurs at a drainage areas of 9.8× 105 m2, or roughly 1 km2. In the geomorphology
literature this is a commonly applied threshold area for the onset of fluvial
processes, and –perhaps reassuringly– this is reflected in the continued evolution
on the curvature metrics. The increasingly negative Gaussian curvature (Panel
C) is an expression of a greater tendency towards synformal saddles, which is the
morphology of a stereotypical channel structure.
In spite of this we see that actually the fluvial channels are almost evenly
partitioned into basins and saddles (Panel H), which reflects fine-scale oscillations
in the Gaussian curvature that are obscured by the binning process, but are
observable in the map-view representation of the data (Panel F). This points out
a major weakness in area space representations of landscape data, as binning over
orders of magnitude does not allow for the observation of fine-scale structure that is
likely important for understanding processes, or the interactions between different
sections of the landscape.
We also note that in this region the different metrics for slope give very
different distributions. While it seems clear that the 1-d channel gradient is the
best for understanding the slope of channels it is less clear how to best connect this
measure to the measure of tangent slope on which we depend for our correlation
between slope and principal curvature vectors. This points out a major challenge
for slope-area analysis, which is that it may not be possible to capture the subtlety
84
Figure 21. Surface geometry data for points in the landscape with
upstream drainage areas between 4.5 × 103 m2 and 9.8 × 105 m2. The red
box in each of the area plots highlights the region of interest. A. PDF of upstream
drainage area. B. Slope metrics binned by upstream drainage area. C. Surface
curvatures binned by upstream drainage area. D. Principal curvature values
and slope alignment binned by upstream drainage area. E. Conditional PDFs of
shape classes as a function of upstream drainage area. F. Map-view of shape class
distribution. G. PDFs of slope metrics. H. Pie chart distribution of shape classes.
85
of the full range of processes using a single metric. It is thus a bit unclear exactly
what is being shown by slope-area plots across process regimes with variable surface
geometry. Finally we note that despite fluvial channels being the connection
between the landscape and base level that likely sets the pace of erosion, fluvial
channels make up less that 1% of the landscape by area.
2.9.2 Mean curvature based landscape partitioning.
2.9.2.1 Mean curvature inflection point. In the previous section
we laid out one possible approach for classifying landscape process transitions in
terms of inflection points and extrema of the Gaussian curvature. The distribution
of Mean curvatures in area space is much less complex, but can also provide
interesting insight into the organizational structure of the landscape, though only
with resolution allowed by averaging into area bins.
The Mean curvature has a single inflection point, which occurs at drainage
areas of ∼ 750 m2 within the transitional hillslope regime classified above as
region 2. We decompose the landscape into two regions separated at this point
with results shown in Figures 24 and 23. The existence of linear hillslope regions is
well documented in the literature, and their relationship to hilltop curvatures have
convincingly been tied to rates and style of erosion (Roering et al. (1999, 2007)).
In the context of our work we see that in fact an inflection point, and thus region
of linearity, is required for maintaining continuity between dominantly convex and
concave regions of the landscape outlined above.
Not surprisingly this inflection point is aligned with the highest slopes in
the landscape, which fall at the same drainage area in all three of the slope metrics
evaluated in this study. The slope at this point in the landscape likely represents
an average angle of repose, above which slope are gravitationally unstable, and
86
Figure 22. Surface geometry data for points in the landscape with
upstream drainage areas greater than 9.8 × 105 m2. The red box in each of
the area plots highlights the region of interest. A. PDF of upstream drainage area.
B. Slope metrics binned by upstream drainage area. C. Surface curvatures binned
by upstream drainage area. D. Principal curvature values and slope alignment
binned by upstream drainage area. E. Conditional PDFs of shape classes as a
function of upstream drainage area. F. Map-view of shape class distribution. G.
PDFs of slope metrics. H. Pie chart distribution of shape classes.
87
Figure 23. Surface geometry data for points in the landscape where the
Mean curvature (KM) is majority positive. The red box in each of the area
plots highlights the region of interest. A. PDF of upstream drainage area. B.
Slope metrics binned by upstream drainage area. C. Surface curvatures binned by
upstream drainage area. D. Principal curvature values and slope alignment binned
by upstream drainage area. E. Conditional PDFs of shape classes as a function of
upstream drainage area. F. Map-view of shape class distribution. G. PDFs of slope
metrics. H. Pie chart distribution of shape classes.
88
Figure 24. Surface geometry data for points in the landscape where the
Mean curvature (KM) is majority negative. The red box in each of the area
plots highlights the region of interest. A. PDF of upstream drainage area. B.
Slope metrics binned by upstream drainage area. C. Surface curvatures binned by
upstream drainage area. D. Principal curvature values and slope alignment binned
by upstream drainage area. E. Conditional PDFs of shape classes as a function of
upstream drainage area. F. Map-view of shape class distribution. G. PDFs of slope
metrics. H. Pie chart distribution of shape classes.
89
below which material will tend to collect as colluvium. This is consistent with the
curvature data which shows that at length scales represented by our analysis the
average geometries upslope of this demarcation line are convex (Figure 23) and
thus steepen downhill, while below this threshold geometry tends to be concave (
Figure 24) with decreasing slope downhill. This can be seen in Panel E, F, and H
of Figures 23 and 24 where domes and antiformal saddles give way to basins and
synformal saddles roughly at this point.
This inflection point also aligns with the most probable drainage area in the
study area and, though area is not normally distributed, the landscape is equally
partitioned about this inflection point such that 50% of points in the landscape
sit both above and below this threshold. This remarkable landscape symmetry
is further explored in Figure 25 which shows a map-view representation of the
landscape divided into concave and convex regions as defined by this area threshold
and shows the distributions of tangential slope (ST ), Mean curvature (KM), and
Gaussian curvature (KG) in these respective portions of the landscape. We see that
there are nearly identical distributions of slope above and below this area threshold,
and that the Gaussian curvature distributions are also nearly aligned while the
distributions of mean curvatures roughly mirror each other.
This work is still too preliminary to make any sort of statement as to the
meaning of this inflection point, but this strong signal of symmetry and geometric
balance in the landscape seems profound enough to in itself warrant further study
that explores how this inflection point varies between landscapes.
2.9.2.2 Mean curvature based channel threshold: Drainage
areas greater than 1 × 105 m2. Other than the inflection point discussed
in the previous section the Mean curvature (KM) only has one major transition
90
Figure 25. Distribution of surface geometry metrics for regions defined by
the Mean curvature inflection point. A. Map-view of landscape partitioned
about the point of inflection in KM . B. Distribution of tangent slope in regions
of both negative and positive average KM . C. Distribution of Mean curvature in
regions of both negative and positive average KM . D. Distribution of gaussian
curvatures in regions of both negative and positive average KM
91
point, which is a minimum at ∼ 1 × 105 m2. While this turning point is not very
pronounced in topography filtered to 200 m, it can be seen in figure 9 that the
existence of an extreme value at approximately these drainage areas is a robust
feature of the data at a range of scales. Figure 26 shows the results of imposing
a threshold at this drainage area. As would be expected this region is composed
of synformal saddles and basins consistent with the channel network, and with a
relative distribution that is between that of regions 4 and 5 defined by the Gaussian
curvature.
While partitioning the landscape this way does not stand out from the other
channel network thresholds in terms of curvatures, it is notable that it aligns with
a scaling break in the slope-area data, which has often been associated with process
transitions (Montgomery and Foufoula-Georgiou (1993)), and so this point may
represent the onset of fluvial processes. Conceptually this makes sense as we would
expect narrow channels at this point in the landscape, which would correspond
with high Mean curvature values.
2.10 Other Potential Applications
2.10.1 Identifying debris-flow source regions. One way to define
the heads of channel networks is through the identification of colluvial hollows
that make up the transition regions between hillslopes and first order debris flow
channels. Colluvial hollows are of broad interest since they often are areas where
hillslope sediments collect that are ultimately mobilized in debris flows (Stock
and Dietrich (2003); Struble et al. (2023)), and so they have direct connection
to debris flow hazards. Morphologically colluvial hollows tend to be regions of
omnidirectional upward concavity, and so we expect them to present in our shape
classification as basins, and more specifically the basins containing the smallest
92
Figure 26. Surface geometry data for points in the landscape with
upstream drainage area greater than 1 × 105 m2 as determined by the
global minimum in the mean curvature. The red box in each of the area
plots highlights the region of interest. A. PDF of upstream drainage area. B.
Slope metrics binned by upstream drainage area. C. Surface curvatures binned by
upstream drainage area. D. Principal curvature values and slope alignment binned
by upstream drainage area. E. Conditional PDFs of shape classes as a function of
upstream drainage area. F. Map-view of shape class distribution. G. PDFs of slope
metrics. H. Pie chart distribution of shape classes.
93
accumulated area. We can identify potential colluvial hollows by selecting basin
structures on the basis of accumulated area. Since channelization is found in
section 2.9 to become established at drainage areas between 4.5 × 103 m2 and
9.8 × 105 m2 we assign a threshold drainage area of 1 × 105 m2 and extract all
basin structures that do not include points with drainage areas greater than the
threshold. The results of this exercise are shown in Figure 27.
2.11 Conclusions
In this study we have applied our method of surface classification to the
simplest possible landscape in the hopes of identifying structures that define
steady-state fluvial topography. It is certain that the magnitudes and distributions
of each of the metrics defined here vary between landscapes, and depend on a
host of variables including erosion rate, climate, lithology, ecology, etc. A natural
next step in developing these tools is to explore how variation in these factors is
expressed in curvature, and how signals of topographic disequilibrium are recorded
in the system of land surface curvatures.
In addition we have compared both curvature and slope calculated using
formal surface theory with approximations common in the geomorphology
literature. We see that slope and curvature dependent errors in these metrics can
be avoided with our more rigurous approach to classi
2.11.1 Bridge. In Chapter II we derived a system of land surface
curvatures and showed how they are expressed in a landscape defined by fluvial
equilibrium. In Chapter III we show the applicability of this methodology to more
complicated landscapes by applying it to a bedrock chronosequence in the Oregon
Cascades, where time dependent hydrology in volcanic bedrock is recorded in the
94
Figure 27. Map of colluvial hollows inferred from our analysis. Purple outlines
mark basin structures that only containing points with areas less than an assigned
threshold of 1 × 105 m2. These are overlain on a shape-class map of region 4, which
we interpret to represent the uppermost reaches of the channel network
.
95
degree of surface dissection by fluvial channels. We show that this hydrologic and
topographic transition is reflected in curvature metrics derived in Chapter II.
96
CHAPTER III
CURVATURE SIGNATURES OF VOLCANIC BEDROCK AGING IN THE
CENTRAL OREGON CASCADES
Along most continental subduction margins the highest and most
geologically active topography is located along the volcanic arc front, however there
exist few tools for studying the long-term dynamics of such terrains. As a result
the potential of geomorphology as a framework for studying volcanic landscape
evolution has not been fully realized. However a growing body of work suggests
that with the proper tools and approaches geomorphology has the potential to
identify signatures of volcanic forcing on the landscape and thus could revolutionize
our understanding of the dynamics of volcanic landscape evolution and thus
magmatic processes (O’Hara, Karlstrom, and Roering (2019); O’Hara, Klema,
and Karlstrom (2021); Singer et al. (2018)). One reason this is challenging is that
basaltic bedrock has a non-trivial time-dependent permeability structure and so
general approaches in geomorphology that correlate erosional efficiency to upstream
drainage area are not easily applied. In this chapter we show how curvature metrics
derived in Chapter II can be applied to a basaltic bedrock chronosequence in the
central Oregon Cascades to infer timescales of a well-documented state-shift from
groundwater to surface water dominated hydrology, which is correlated with the
evolution of volcanic constructional topography to a fluvially dissected landscape
with signatures of equilibrium. Since this approach does not require assumptions
regarding the partitioning of meteoric water between the surface and subsurface
it provides an objective way of quantifying volcanic topography towards better
understanding the mechanisms and timescales of volcanic landscape evolution.
97
3.1 The Cascades Chronosequence
In the central Oregon Cascades bedrock age increases westward from the arc
axis (Figure 28) reflecting an eastward migration of the main locus of volcanism
over timescales of tens of Myr (du Bray and John (2011)). Arc front migration
relative to the plate boundary is common in subduction zones (Dickenson and
Snyder (1979); Karlstrom, Lee, and Manga (2014)), and in the context of the
Cascades is likely due to the clockwise rotation of the Pacific Northwest fore-arc
block relative to North America (R. E. Wells and McCaffrey (2013)) combined with
potential climatic and magmatic drivers (Hildreth and Wilson (2007); Karlstrom et
al. (2014); V. A. Muller, Sternai, Sue, Simon-Labric, and Valla (2022)). As a result
the Oregon Cascades are often referenced in terms of two geologic subprovinces; the
older Western Cascades, and the younger High Cascades where volcanic bedrock is
associated with the modern arc (du Bray and John (2011)). While this bedrock
chronosequence is defined on the basis of detailed field mapping (Sherrod and
Smith (2000)) fundamental differences between the Western and High Cascades
are also reflected in both hydrology and topography, which themselves are tightly
coupled. Thus the Cascades chronosequence is an ideal natural laboratory for the
study of time dependent hydrologic processes in volcanic landscapes, and how the
are recorded in the evolution of topography.
3.1.0.1 Groundwater Dominated High Cascades. Young basaltic
bedrock has high bulk permeability meaning that meteoric water in the High
Cascades is mostly absorbed into shallow aquifers and there is little surface runoff
(Saar and Manga (2004)). As a result the High Cascades contain few fluvial
channels and topography is mostly defined by primary volcanic structures and
glacial valleys carved during the LGM (Deligne et al. (2017); Sweeney and Roering
98
Figure 28. Map of the central Oregon Cascades. A. Topography. Red boxes
indicate areas of focused analysis in Figs. 30 and 33. B. Age of volcanic bedrock
units (from Sherrod and Smith (2000))
(2017)). Where stream channels exist in the High Cascades they are mostly spring
fed with damped hydraulic responses to precipitation indicating the dominance of
stable groundwater systems in modulating their discharge (Manga (1999); Tague
and Grant (2004)). This can be seen in Figure 29, which is taken from Deligne
et al. (2017) and compares an annual hydrograph of the McKenzie River to that
of the Little North Santiam River. The more muted response of the McKenzie
to precipitation and run-off events is characteristic of spring fed channels (Tague
and Grant (2004)). In addition there also exist in the High Cascades ephemeral
channels that have rapidly incised into lava flows and show the relevance of
standard geomorphic processes in the development of volcanic terrains (Sweeney
and Roering (2017)), however these channels are only locally similar to those of
non-volcanic landscapes, possibly representing an intermediate stage of fluvial
landscape development.
99
Figure 29. Characteristic hydrographs of the spring dominated High
Cascades, and the surface run-off dominates Western Cascades (from
Deligne et al. (2017) which is a modified figure from Tague and Grant (2004)
100
3.1.0.2 Surface Water Dominated Western Cascades. Over
time there is a decrease in the bulk permeability of basaltic bedrock as micro-pores
are clogged with fine sediment such as glacial till (Sweeney and Roering (2017)),
volcanic ash (Jefferson et al. (2010)), and soil particulate formed through en situ
chemical weathering (Lohse and Dietrich (2005)). The diminishing infiltration
capacity corresponds with an increase in overland flow and a correlated increase
in the ability of fluvial channels to mobilize sediment (Jefferson et al. (2010)). As a
result the aging landscape is increasingly defined by channel networks with slope-
area scaling relationships and curvature distributions similar to those defined in
Chapter II for the Oregon Coast Range. Figure 30, shows the same presentation
of surface metrics given in Chapter II, but calculated on a section of the Western
Cascades in which bedrock has an age of 20 ± 2.9 My (Sherrod and Smith (2000)).
It can be seen that while the magnitudes of these metrics - and the exact locations
of inflection points and extrema - are slightly different than those of the Coast
Range, the landscape has similar area-space trends and geometric distributions.
This implies the Western Cascades can be considered an equilibrium landscape
despite its complex initial conditions and evolutionary trajectory.
3.2 Signatures of Topographic Development
The idea of leveraging the delayed development of fluvial channels in the
Cascades as a means of understanding the hydraulic state-shift is not new and
builds directly on the work of Jefferson et al. (2010), who used a combination of
drainage density, bedrock age, and stream base flows to show that processes driving
the hydraulic state-shift occur over 100 ky timescales, with the full transition to
fluvial topography taking ∼ 1 Ma. Figure 31 shows this result in the form of a
basin averaged drainage density plotted against basin averaged bedrock age.
101
Figure 30. Compilation of surface geometry data used in Chapter II
applied to the Western Cascades. A. PDF of upstream drainage area. B.
Slope metrics binned by upstream drainage area. C. Surface curvatures binned by
upstream drainage area. D. Principal curvature values and slope alignment binned
by upstream drainage area. E. Conditional PDFs of shape classes as a function of
upstream drainage area. F. Map-view of shape class distribution. G. PDFs of slope
metrics. H. Pie chart distribution of shape classes.
102
Figure 31. Plot of drainage density versus bedrock age in the upper
McKenzie River Basin (from Jefferson et al. (2010))
A spatial correlation between bedrock age and drainage density can be
observed throughout the Cascades, and is seen in Figure 32 which shows map-
view representations of bedrock age and drainage density (calculated off of USGS
Bluelines) averaged over several spatial scales using a simple moving window
filter. The problem with this method for inferring bedrock age and hydrologic
characteristics is that it requires making a choice as to what defines a channel,
which is challenging in the complex topography of volcanic landscapes. This was
demonstrated by Luo and Stepinski (2008), who compared different methods of
drainage network extraction in the Cascades and showed that drainage density
is highly dependent on classification method. They also found that a curvature
based algorithm that does not depend on any drainage area assumptions (Molloy
and Stepinski (2007)) most accurately recreates the channel network as defined
by field mapping. Still it does not provide a tool for more detailed morphometric
analysis in un-channelized portions of the landscape along the arc axis. This
103
Figure 32. Maps of surface bedrock age in the Oregon Cascades averaged
over different spatial scales.
relates to another challenge in applying the methods of Jefferson et al. (2010) more
generally, which is that both drainage density and bedrock age are averaged within
topographic watersheds, which may not always align with hydraulic divides in
volcanic landscapes where relict topography is overprinted with permeable surface
lithologies (Gannett et al. (2001)).
3.2.1 Curvature based surface classification. The surface
classification method derived in Chapter II carries no inherent assumptions about
process and so can be used as an objective measure of topographic form across the
entire Cascades chronosequence. Area-space representations are not appropriate
in the High Cascades due to disconnects between topographic and hydrologic
watersheds (Gannett et al. (2001)), and so we focus on the surface evolution
relative to bedrock age. We first explore the fundamental differences in topographic
104
form by comparing the distributions of surface geometry metrics in the High and
Western Cascades to those of the Oregon Coast Range. In this section we focus
on the Mean curvature (KM) , Gaussian curvature (KG), the tangent slope of the
surface (ST ), and the surface area proportionality constant (α). α is a measure of
the change in area between the surface and its map projection that follows directly
from the first fundamental form, and which we interpret to be a local measure of
topographic relief.
Figure 33 shows PDFs for each of these metrics calculated on 10 m DEMs
filtered to 200 m for direct comparison with the analysis in Chapter II. Even this
very general presentation of the geometric data gives insight into the fundamental
differences between fluvial and volcanic-constructional landscapes. While curvatures
are normally distributed in all of the landscapes the narrow distributions in the
High Cascades indicate much lower magnitude curvatures than the Western
Cascades, which have wider distributions that are more similar to those of the
Coast Range. This possibly suggests that the widening of curvature distributions is
a temporal signature of landscapes trending towards fluvial equilibrium. A similar
trend is seen in the distributions of both slope and α, with overall lower magnitude
values in the High Cascades, and higher magnitude values in the Western Cascades
and Coast Range.
3.2.2 Timescales of Fluvial Development. We now apply our
analysis to the full Cascades chronosequence to give a more detailed view of how
the hydrologic state-shift is recorded in land surface geometry. Surface metrics are
calculate throughout the full study area and are binned both by UTM Easting and
by age with the results shown in Figures 34-35. We first note that this analysis
confirms the result from Jefferson et al. (2010) that the topographic complexity
105
Figure 33. Surface metric distributions for sections of the Oregon Coast
Range, Western Cascades, and High Cascades. A. Oregon Coast Range
topography. B. Western Cascades topography. C. High Cascades topography. D.
Gaussian curvature PDFs. E. Mean curvature PDFs. F. Slope PDFs. G. α PDFs.
106
associated with fluvial dissection emerges in the Cascades when volcanic bedrock
ages reach about 1 Ma. This same timescale is reflected in the increased magnitude
of all geometric parameters around that time. In terms of easting the transition
is more gradual, but we note a sharp increase coincident with major arc bounding
normal faults (R. M. Conrey, Taylor, Donnelly-Nolan, and Sherrod (2002)). This
possibly indicates some structural control on the state-shift, which could be due
to the build up of young basaltic lavas along the margins of the arc (Sherrod and
Smith (2000)). The bell shape of the surface metrics binned by easting is not easily
explained, and could possibly be connected to subsidence of the Willamette Valley
(Mitchell, Vincent, Weldon, and Richards (1994)) or orographic precipitation.
However there is also similarity in the shape and wavelength of these trends to
hyaloclastite deformation in the Columbia River Gorge Tolan and Beeson (1984),
which we explore in depth in Chapter IV. We note that flexural deformation of
bedrock has been identified in limited sections of the Western Cascades (Lopez
and Meigs (2016); Pesek, Perez, Meigs, Rowden, and Giles (2020)), and so there
is a possibility that these surface geometry trends record signatures of structural
deformation that are not widely recognized.
3.3 Further Directions
In this chapter we have demonstrated a simple way in which surface
geometry can be used to study the evolutionary trajectory of arc topography in
the Oregon Cascades. This analysis was conducted over large spatial scales in an
attempt to extract statistical characteristics of topography, however geometry is
inherently scale dependent and so there is much to be learned by such analysis at a
range scales. For example it is likely possible to identify signatures of topographic
forms not easily described by standard geomorphology tools such as stratocones,
107
Figure 34. Gaussian curvature of the Cascades chronosequence. A. Map-
view of calculated Gaussian curvature. B. Gaussian curvature binned by bedrock
age. C. Gaussian curvature binned by UTM easting.
Figure 35. Mean curvature of the Cascades chronosequence. A. Map-view
of calculated Mean curvature. B. Mean curvature binned by bedrock age. C. Mean
curvature binned by UTM easting.
108
Figure 36. Tangent slope of the Cascades chronosequence. A. Map-view
of calculated tangent slope. B. Tangent slope binned by bedrock age. C. Tangent
slope binned by UTM easting.
Figure 37. α of the Cascades chronosequence. A. Map-view of calculated area
proportionality factor (α). B. Area proportionality factor (α) binned by bedrock
age. C. Area proportionality factor (α) binned by UTM easting.
109
lava flows, and glacial landforms. We have also stopped short of describing the
evolution of the Cascades in terms of the shape classes derived in Chapter II,
or looked in detail at how the distributions and orientations of the principal
curvatures evolve in aging volcanic bedrock. The different shape in the curvature
distributions between the Western and High Cascades also make this an ideal
setting for the application of curvature thresholds in defining surface shape classes,
which could likely be used to more rigorously define the geometric evolution of the
Cascade Arc.
3.4 Conclusion
In this chapter we have shown how the transition from a volcanic
constructional to a fluvially dissected landscape in the Oregon Cascades is recorded
in land surface geometry. Our approach confirms past estimates of a state-shift
from ground to surface water dominated hydrology occurring over ∼ 1 My with a
correlated development of ridge-valley topography. Next steps in this work include
more sophisticated classification including an analysis of shape class distributions as
was done for the Oregon Coast Range in Chapter II.
3.5 Bridge
In Chapters II and III. we explored a new method for land surface curvature
classification based on classical surface theory. We now use a more traditional
geomorphology approach of knickpoint analysis to disentangle the volcanic history
of the Cascade Arc seen in cross-section in the Columbia River Gorge. While this
setting is complicated in the context of landscapes in which knickpoint analysis is
generally applied, we show that an interdisciplinary approach leveraging geologic
and geophysical data gives considerable insight into the evolution of mountain
range scale topography in this setting allowing for unique insights into long-
110
term trends in magmatic processes and how magmatism contributes to landscape
evolution in arc settings. This work is taken from a manuscript in review, and so
is comprised of main body text in addition to supporting information contained in
Appendix A.
111
CHAPTER IV
THE MAGMATIC ORIGIN OF THE COLUMBIA RIVER GORGE, USA
Reproduced with permission from Klema, N.T, Karlstrom, L., Cannon,
C., Jiang, C., O’Connor, J., Wells, R., Schmandt, B., (in review). The Magmatic
Origin of the Columbia River Gorge, USA. Science Advances
4.1 Introduction
Volcanic arcs represent the highest, steepest, and the most rapidly
changing topography at subduction zones. Localized and unsteady uplift driven
by magmatism, superimposed on background tectonic deformation, makes these
terrains an outstanding laboratory for the study of volcanic mountain building
(Gregory-Wodzicki (2000); Lee, Thurner, Paterson, and Cao (2015)). For over
100 years the Columbia River Gorge (hereafter the Gorge), an enigmatic high
relief landscape that has shaped human society in the region (Lyman (1963)), has
been recognized as an ideal setting for the study of arc processes (Bretz (1917);
I. A. Williams (1923)). As the largest river on Earth to cross an active volcanic arc,
the Columbia River has dissected the Cascades for at least 17 Myr (Evarts, Conrey,
Fleck, and Hagstrum (2009)), exposing a cross section of the Cascade Range and
serving as a near-sea-level base level driving erosion of the surrounding landscape.
4.2 The Columbia River Gorge
As the regional topographic low, the Columbia River corridor has received,
conveyed, and been diverted by lava flows from Columbia River Flood Basalt
volcanism to the east and locally from the Cascade volcanic arc that it bisects
(Reidel, Camp, Tolan, Kauffman, and Garwood (2013)). Remnants of these
volcanic deposits record subsequent deformation of the arc by their position relative
to the modern Columbia River (O’Connor et al. (2021)). An especially prominent
112
structural datum is marked by the most recent Columbia River diversion in the
Gorge at 3.5−3.0 Ma, when arc-derived lava flows and associated hyaloclastite filled
the existing valley and displaced the river northward to its modern course (Beeson
and Tolan (1990)). Hyaloclastite deposits filling this “Bridal Veil” (BV) paleo-
channel (Figure 38) record as much as ∼ 900 m of subsequent uplift in a broad,
asymmetric anticline spanning ∼ 80 km of the arc (Beeson and Tolan (1990)). A
superposition of this uplift and overlying erupted deposits forms the iconic high-
relief topography of the Gorge, which locally rises to 1512 m at the shield volcano
Mount Defiance.
Clockwise rotation of the forearc relative to North America (Figure 38)
since at least the early Miocene (R. E. Wells and McCaffrey (2013)) provides a
framework to analyze the underlying drivers of this uplift. North of the Columbia
River, tectonic rotation produces compression, mainly manifest by east-west
oriented folds and thrust faults (Reidel et al. (2013)). South of the Columbia
River the regime is mostly extensional as evidenced by north-south oriented normal
faults. Rotational shear is accommodated by right lateral displacements on north to
north-northwest strike-slip faults (R. M. Conrey et al. (2002); R. E. Wells (1998);
R. E. Wells and McCaffrey (2013)). The Columbia River thus sits at a transition
point where tectonic strain is dominantly north-south, making crustal shortening an
unlikely source of the broad uplift centered on the Gorge. Instead, we hypothesize
that the uplift is due primarily to magmatic crustal inflation. We will show that
late Pliocene and Quaternary deformation of BV channel deposits is a flexural
response to magmatic intrusions below an elastic upper crust. Magmatic forcing
is active at present, as inferred from gravity, heat flow, and seismic tomography
coincident with uplift, and patterns of Quaternary volcanism. Current Columbia
113
Figure 38. Overview of the Columbia River Gorge. A. Structural setting,
showing GNSS velocities with respect to stable North American reference frame
(McCaffrey et al. (2013)) (black arrows), Yakima Folds (orange lines), the modern
and paleo (‘Bridal Veil’) Columbia River (blue and pink lines), and extent of
the Cascades arc (red dashes). B. Symbols are Bridal Veil Channel hyaloclastite
deposits and best-fitting knickpoint locations in Columbia tributaries. Maximum
topographic elevations within 10 km south of the Columbia (black line) with
prominent faults (dashed lines) (Mcclaughry et al. (2012)). C. Southward
orthoview of fluvial catchments, hyaloclastite deposits, and major faults. D.
Orthoview of Oneonta Creek showing the location of the trunk stream knickpoint
used in this study (large blue square), analogous knickpoints in smaller tributary
channels, and waterfalls associated with layered basalt stratigraphy (black stars).
E. Regional map of the Columbia river drainage network, showing the Euler
pole for Cascadia forearc block relative to North America. Orange lines are
compressional structures while yellow lines mark arc-marginal normal faults.
F. Longitudinal profile using area-normalized upstream distance χ for Oneonta
creek showing the location of the best fitting slope-break knickpoint (Supplement).
Lithologic units are labeled.
114
tributary drainages record transient topographic adjustment to this spatially
variable magmatic uplift field in the form of fluvial knickpoints that mirror uplift
of the BV paleo-channel.
4.3 Transcrustal magmatic structure
The mechanical response of Earth’s crust to an applied load depends on
the magnitude and spatial wavelengths of forcing as well as rheological properties
that are strongly time and temperature dependent (Audet, Jellinek, and Uno
(2007); Burov and Diament (1995)). In the case of an elastic response to vertical
loading the lateral wavelength of deformation is set by a flexural rigidity that
scales with the effective elastic thickness of the crust (Te) cubed. In a maturing arc
setting, heat released from generations of magmatic intrusions ascending through
a transcrustal transport network (Sparks and Cashman (2017)) likely induces a
shallowing of Te as crustal heating promotes ductile mechanical response at depth
(Karlstrom, Paterson, and Jellinek (2017)). Intrusions dilate overlying elastic crust,
producing a mechanical response similar to a thin elastic plate with thickness Te set
by this rheologic transition (O’Hara et al. (2021); Wilson and Russell (2020)).
The magnitude and pattern of uplift recorded in BV channel hyaloclastites
is well-modeled with Te = 4 − 8 km subject to a net upward force of ∼ 34 × 106
N distributed over a wavelength of 110 − 130 km. Removal of crust by erosion
by the Columbia River and its tributaries contributes to this force imbalance, but
not enough for isostatic rebound to be a primary driver of uplift (section A.0.5).
We estimate the removal of 270 km3 of sediment by erosion in tributaries of the
Columbia that enter within the Gorge. For values of Te that fit the deformation
wavelength of hyaloclastite deposits the implied removal of mass would contribute
∼ 120m of flexural uplift at most. This unloading, however, is offset by a
115
comparable volume of erupted material deposited on top of the paleo-surface by
Cascade Range volcanoes. We therefore neglect erosional unloading relative to
magmatic inflation in forcing hyaloclastite uplift.
Inference of intrusive magmatic forcing as a mechanism for bedrock uplift
is then augmented by a broad suite of indirect observations. Heat flow and
gravity data implicate elevated temperatures and silicate melt under the area
of maximum hyaloclastite uplift both at present and likely quasi-continuously
since the onset of hyaloclastite deformation (Figure 39). Bouguer gravity data
(Bonvalot et al. (2012)), although not corrected for locally steep topography or
laterally heterogeneous subsurface structures, shows a negative anomaly over the
arc consistent with flexural uplift due to a buried load (section A.0.3). Similarly,
regionally interpolated heat flow measurements (Ingebritsen and Mariner (2010))
exhibit elevated near-surface temperatures consistent with a long-lived magmatic
heat source centered on the arc.
High resolution seismic tomography using Rayleigh and Love waves (Jiang,
Schmandt, Abers, Kiser, and Miller (2023)) shows a region of low shear-wave
velocities (Vs) at depths > 20 km beneath the arc with magnitudes that likely
require both partial melt and elevated temperatures (Supplement). The seismic
anomaly is aligned with hyaloclastite deformation in the Gorge and extends
northwest and southeast connecting Mt. St. Helens and Mt. Hood at 20-30 km
depth (Figure 39; Jiang et al. (2023)).
4.4 Long-term balance of intrusions versus eruptions
Broad covariance along the Cascade Range between crustal geophysical
anomalies and high topography associated with Quaternary volcanic vents (O’Hara,
Karlstrom, and Ramsey (2020)) suggests that magmatism contributes to uplift
116
Figure 39. Geophysical anomalies across the Gorge and their relationships
to BV channel deformation. A. Best-fitting uplift model from maximum
likelihood estimation with 95% confidence model envelope. Light blue squares
show knickpoint elevations inferred from a transient landscape evolution model
forced by flexural uplift constrained by hyaloclastite deformation (pink squares).
B. Interpolated Bouguer gravity and surface heat flow spanning the Gorge. C.
Average surface wave velocity anomaly. The white polygon marks inferred 2.5 − 4%
partial melt (Supplement). Black dots are earthquake hypocenters for all events
underlying the Gorge in the Pacific Northwest Seismic Network catalog since
1970. D. Hyaloclastite elevation (squares) and knickpoint elevation (diamonds)
vs knickpoint location χkn. Red dashed line shows uplift predicted by fluvial model
using slope-area regression, while shaded swath shows model predictions from
a coupled flexural uplift/landscape evolution model (Supplement). Mismatch
between χkn and knickpoint elevation likely implies time-evolving hydraulic
conductivity/erodibility of young lava flows. E. Regional map of shear wave seismic
velocity at 30 km depth showing swath in panel C.
117
throughout the arc. However understanding the partitioning between intrusive and
extrusive (i.e., eruptions) modes of magmatism is challenging (White, Crisp, and
Spera (2006)), so rates of topographic growth from intrusive magmatic processes
are difficult to quantify. The cross-section of the Cascade Range created by the
Columbia River provides a unique opportunity to calculate the ratio of intrusive
to extrusive magmatism (I:E) on million-year timescales, and assess the relative
contribution of each to arc mountain building.
A lower bound estimate for the volume of intruded magma can be calculated
as the amount of accommodation space created by hyaloclastite deflection between
the western margin of uplift and the fault-bounded eastern margin of uplift (Figure
40). Total extruded volcanic output is then the difference between elevation of the
paleosurface and the modern topography. If we assume an onset time between
3.0 − 3.5 Ma (section A.0.1) this gives average intrusive and extrusive magma
flux estimates of 10 − 13 km3/My/km and 3 − 4 km3/My/km, respectively.
Accounting for spatially variable topography within a 10 km swath south of the
Gorge (assuming intrusive uplift is unchanged) gives a distribution with 95%
confidence interval of I:E = 2.9+3.6−2.4 (Figure 40).
We can compare this to previous estimates of magmatic flux required to
produce heat flow measured in the Cascades (Ingebritsen, Mariner, and Sherrod
(1994)) of 9 − 33 km3/My/km. Taken with eruptive flux estimates of 3 − 6
km3/My/km (Sherrod and Smith (2000)) this gives I:E= 1.5 − 11. Topographically
derived I:E values for the Gorge are a lower bound, because they do not account
for deep intrusions which generate longer wavelength deformation near the surface
(O’Hara et al. (2021)), and only account for magma underlying the deformed paleo-
surface.
118
Figure 40. Inference of transcrustal magmatic structure. A. Volume
weighted Gaussian kernel density of Quaternary edifices within 10 km of the study
area. B. Intrusive and extrusive magmatic topography, using best fitting flexural
uplift model. Dashed lines indicate major fault systems that bound the Hood River
Valley (Mcclaughry et al. (2012)). C. Best fitting plate model defining elastic
thickness of the upper crust. Alternating tensile and compressional deviatoric
stresses in the plate are inferred to direct shallow magma ascent and faulting.
D. Ratios of intrusive:extrusive (I:E) magmatism calculated using topography
(assuming north-south continuity in flexural uplift), knickpoint positions, and total
calculated volumes of erupted and intruded magma.
119
Seismic tomography of present day crustal structure suggests ongoing
intrusion and provides modern context for these flux estimates. The volume
of stored partial melt implied by seismic anomalies depends on lower crustal
compositions, temperature, and grain boundary melt geometry (Takei (2002)).
Petrologic constraints on storage temperatures for Mt. St. Helens magmas suggest
lower crustal temperatures in the range of 750 − 950 C (Wanke (2019)). A magma
storage zone with this temperature range is consistent with observed heat flow
at the arc front, which also suggests sustained heat from partial melt at depth
(section A.0.6). Assuming textural equilibrium of melt in the deep crust, seismic
tomography suggests ∼ 4 km3 magma per kilometer of arc currently underneath the
Gorge.
4.5 Topographic response to magmatism
We can gain additional insight into the time-dependent rates of magmatism
by leveraging the stable, near-sea-level base level imposed by the Columbia River
as it transects the Cascade Range. Actively eroding landscapes tend towards
dynamic equilibrium conditions wherein surface erosion and uplift are balanced and
fluvial channels maintain power-law scaling between channel slope and upstream
drainage area (Flint (1974)). An increase in uplift rate generates an upstream-
propagating kinematic wave of incision that manifests as a step-change in this
slope-area relationship often referred to as a ‘slope-break knickpoint’ (Whipple
et al. (2013)). Thus we expect that a past increase in the rate of magmatic uplift
should be recorded in the slope-area relations of tributaries entering the Columbia
River.
We focus on channels that enter from the south side of the Gorge where
southward tilting bedrock limits deep-seated landsliding that is prevalent on the
120
north side of the river (Pierson, Evarts, and Bard (2016)) (section A.0.2). Here
16 prominent tributary streams entering the Columbia River provide a natural
landscape evolution experiment in which adjacent streams are forced by a known
uplift perturbation beginning at the same time (deposition of BV deposits), but
with a magnitude that varies significantly between tributaries. In all catchments,
channels cut through similar layered basalt stratigraphy consisting of Cascades arc
volcanic rock overlying Columbia River Flood Basalts and lenses of hyaloclastite
(Beeson and Tolan (1990)). Variability between basalt layers causes sharp localized
knickpoints in stream profiles (the iconic waterfalls of the Columbia Gorge),
however at the catchment scale most tributaries contain a single slope-break
knickpoint that separates each longitudinal profile into two segments with constant
concavity but different slope-area scaling (Figure 38).
The upstream position of these knickpoints tracks hyaloclastite deformation
(Figure 1), implicating a geomorphic adjustment to magmatic uplift. Knickpoint
propagation speed depends on upstream drainage area (A) (Berlin and Anderson
(2007)), and so to compare different sized catchments directly we normalize
upstream distance along each channel with respect to drainage area and examine
knickpoint positions relative to the non-dimensional upstream coordinate χ
(Royden and Taylor Perron (2013)). The ‘stream power’ model for fluvial bedrock
erosion in its simplest form predicts that a knickpoint associated with an increase
in uplift rate occurring at time τo is given by χkn = Uτo/ks (section A.0.2) where U
is the average uplift rate and ks is an empirical constant that scales the power law
relationship between upstream drainage area and channel slope, dz/dx = ksA(x)
−θ
(Flint (1974)). Maintaining continuity in the elastically deformed upper crust
requires that τo (assumed to coincide with diversion of the Columbia at 3.0-3.5
121
Ma) is the same in all catchments. The stream power prediction for the position
of knickpoints associated with the hyaloclastite uplift (why) can then be written as
χkn = why/ks, which matches the data with ks = 140±17 for all catchments (Figure
39).
The lateral position of knickpoints is less well correlated to knickpoint
elevations in channels (Figure 39.D), which reflect variable Pliocene-Quaternary
arc lavas overlying the paleosurface. This observation is consistent with episodic
deposition of volcanic bedrock, which leads to strong hydromechanical layering and
hydraulic anisotropy that impacts bedrock erosion long after primary emplacement
(Jefferson et al. (2010); Sweeney and Roering (2017)). We infer that young lava
flows in the Gorge produce topographic relief but have negligible influence on time-
averaged rates of channel incision in underlying strata.
Inference of constant ks across Gorge tributaries is consistent with fluvial
erosion patterns in other high-relief landscapes. Although this parameter is
sometimes assumed to covary with uplift rate (Whipple et al. (2013)), a global
compilation of fluvial catchments in active tectonic settings (Hilley et al. (2019))
suggests a threshold in channel steepness that creates an upper limit on ks in
rapidly uplifting landscapes. Increase of coarse sediment supply with uplift,
resulting in transport limited conditions that modulate the efficiency of fluvial
bedrock erosion (Lai, Roering, Finnegan, Dorsey, and Yen (2021)), is one way limits
on ks might be achieved.
The inferred relationship between knickpoint locations and hyaloclastite
uplift shows that fluvial profiles record the relative contribution of intrusion
versus eruption to topographic development. With thickness of erupted deposits
as the difference between observed channel elevations and uplift associated with
122
the knickpoints, the ratio of intrusive to extrusive magmatic uplift can then be
calculated directly from fluvial channels as I:E= ksχkn/zkn = 1.4 ± 0.5. The
alignment of this value with I:E calculated through other means (Figure 3D)
suggests that long-term magma dynamics can be encoded in river networks, as
has been inferred for tectonic and/or climatic forcing in non-magmatic landscapes
(Schoenbohm, Whipple, Burchfiel, and Chen (2004); Wobus et al. (2006)).
4.6 Entwined upper crustal deformation, arc mountain building, and
magma ascent
Direct constraints on the mechanical response of the upper crust to intrusive
magmatism can explain several other primary features of the central Cascades arc.
In particular, plate-like flexure of the elastic upper crust generates alternating
tensile and compressional stresses (Hieronymus and Bercovici (2001)) that may
explain prominent normal faults along the Hood River Fault Zone (Figure 38).
These extensional features are consistent with north-south maximum horizontal
compressive stress, and align with peaks in deviatoric stresses from our best-fitting
flexural model (Figure 40). To the west of the arc front, subsidence of the Portland
basin occurs coincident with the implied flexural bulge of our model (R. Wells et al.
(2020)). While flexure cannot itself explain the magnitude of subsidence there (or
in the Hood River Valley to the east), the resulting topography would route fluvial
sediments and erupted material to arc margins, causing further loading of the plate
and increased subsidence (Garcia-Castellanos (2002)).
Unusual spatial patterns of magmatism in the vicinity of the Columbia
river, the widest point of the Quaternary Cascades arc (Hildreth (2007)), may
also reflect the influence of shallow crustal stresses implied by our model. Within
the Gorge, local maxima in both the volume-weighted Gaussian kernel density
123
of volcanic vents (O’Hara et al. (2020)) and the thickness of associated erupted
deposits are located east and west of maximum hyaloclastite layer deformation
(Figure 40). This is consistent with compressional flexural stresses in the lower half
of an elastic upper crust that route ascending magma towards off-axis minima in
differential stresses (Hieronymus and Bercovici (2001)). The most active Holocene
Cascades volcano, Mount St. Helens, sits notably off-axis in the forearc northwest
of the Gorge. Arc-centered flexure can contribute to observed lateral migration of
ascending magma around inherited upper crustal structures (Bedrosian, Peacock,
Bowles-Martinez, Schultz, and Hill (2018)).
Outside the gorge, local widening of the Quaternary arc is associated with
the Boring and Simcoe distributed vent fields (Figure 38) (Hildreth (2007)). The
Columbia River itself could play a role in this local widening of arc magmatic
surface expression, as surface unloading from persistent fluvial incision generates
topographic stresses known to influence magma ascent pathways in other settings
(J. R. Muller, Ito, and Martel (2001)). Current topography represents a ∼ 25
MPa topographic load difference between the Gorge and adjacent arc (Garcia,
Sandwell, and Luttrell (2015)). Quaternary average erosion/sediment production
rates implied by projecting the hyaloclastite layer over the Gorge are 90 km3/Myr.
This missing volume represents total magmatic topographic construction, at least
270 km3 or ∼ 6 − 9 km3/My/km based on ∼ 10 − 15 km north-south extent of the
Gorge. This roughly matches inferred deep magmatic influx rates.
Finally, transcrustal evidence ranging from geomorphology of river channels
to seismic tomography of the crust suggest persistent, long-lived magmatism that
is not directly associated with any major stratovolcano. Seismic and geodetic
data show that arc magma reservoirs are commonly offset from their associated
124
volcanoes (Lerner et al. (2020)), suggesting that sustained magma storage between
adjacent volcanoes (as inferred here) is more common than previously recognized
Bedrosian et al. (2018). In general, this implies that arc magma transport networks
are not always strongly aligned with volcanic edifices. Our work suggests that
intrusive magmatism has a broad organizing influence (Karlstrom, Dufek, and
Manga (2009)) spanning surface topographic development, crustal stress state, and
magma ascent pathways.
125
APPENDIX
APPENDIX A - CHAPTER IV SUPPLEMENARY MATERIAL
A.0.1 Constraints on structure and timing of Gorge
deformation. The Columbia River has crossed the Cascades arc for at least 17
Myr though its exact location has shifted several times as the result of volcanic
infilling of the channel. In the Gorge laharic breccias of the Eagle Creek Formation
uncomformably underlie early Columbia River Basalt (CRB) flows, suggesting
a dominance of deposition by fluvial tributary systems draining local volcanic
terrain in the early Miocene (Korosec (1987)). The first known avulsion of the
main Columbia channel occurred through emplacement of the Grande Ronde
Basalts between 16.0 and 15.5 Ma (Reidel and Tolan (2013a)). The ∼150,400 km3
of material associated with this period of CRB volcanism erupted from fissures
in eastern Oregon and inundated the established Columbia River network. This
displaced the river and forced it on a more southern course across the Cascades
along a developing fold structure known as the Mosier-Bull Run syncline (Reidel
and Tolan (2013b)). Subsequent emplacement of the Wanapum basalts forced the
river back north into the Bridal Veil Channel by 12 Ma, where it persisted until it
was pushed to its modern course during a widespread pulse of mafic volcanism in
the northern Oregon Cascades between ∼4.4 and 2.1 Ma (R. Conrey, Sherrod, Uto,
and Uchiumi (1996); Tolan and Beeson (1984)).
In the mid Pliocene a series of voluminous lavas, mostly with low potassium
tholeiitic basalt (LKT) compositions, were erupted from the area south of the
modern Columbia River. Some of these lavas flowed into the ancestral channel of
the Columbia River and were rapidly chilled and fragmented during interaction
with water generating hyaloclastite (Mcclaughry et al. (2012); Swanson (1986)).
126
The Bridal Veil channel was eventually filled by hyaloclastite and associated
LKT lava flows (Figure A.41) and the river was diverted north to its present
course. 40Ar/39Ar age determinations for LKT lava flows inter-bedded with, and
overlying, the hyaloclastite (Fleck, Hagstrum, Calvert, Evarts, and Conrey (2014);
Mcclaughry et al. (2012)) suggest diversion of the channel between 3.5 and 3.0
Ma (Table S1) thus determining our age brackets for the beginning of the current
phase of uplift. During this time arc-marginal normal faulting associated with the
northward propagation of an intra-arc rift began in the Hood River Fault Zone
(R. Conrey, Sherrod, Hooper, and Swanson (1997); Mcclaughry et al. (2012)).
127
Figure A.41. Photos of hyaloclastite deposits associated with the BV
channel. A. Hyaloclastite (orange/yellow) with interbedded pillow complexes
seen from I-84 about 4 km west of Hood River, OR. near Ruthton Point B. Lens
of hyaloclastite material with channel capping LKT flows in Ainsworth State Park
near Dodson, OR. C. Exposures of hyaloclastite material uplifted nearly 1000 m
above the modern Columbia River (visible on right side of panel) seen from the
Pacific Crest Trail near Cascades Locks, OR.
128
Figure A.42. Raw Bridal-Veil Channel elevation data. A. Mean elevations
of hyaloclastite deposits and minimum elevations of LKT flows that overlie the BV
Channel. B. Map view of hyaloclastite and LKT polygons.
A.0.2 Topographic analysis and fluvial erosion. Today the
anticlinal uplift of the hyaloclastites (Figure A.42) is superimposed on a slight
(2 − 3◦) north-south downward tilt observable in CRB deposits and underlying
volcaniclastic sediments. As a result, the north side of the Columbia River is prone
to deep-seated landslides that tend to fail along slip surfaces at unconformable
contacts at the base of the CRB and move into the accommodation space created
by Columbia incision (Pierson et al. (2016)). On the south side of the Columbia,
these slip surfaces are buttressed by southward tilt so deep-seated mass wasting
events are not common. For this reason we focus our topographic analysis on the
south side of the Columbia River. Of the seventeen largest catchments on the south
side of the Columbia one (Tanner Creek) is eliminated as its longitudinal profile
129
has a large knickpoint near its confluence with the Columbia that is inconsistent
with the dominant pattern of transient adjustment. This knickpoint is likely related
to partial filling of the Tanner Creek channel by basalt flows observed near the
channel mouth, but compositional affinity and age are poorly resolved in current
maps limiting our ability to confirm this explanation.
It has been documented that mature fluvial networks exhibit inverse power-
law scaling between channel slope and upstream drainage area at a given point in
the channel that is often interpreted as indicative of a steady-state balance between
uplift and erosion (Flint (1974); Whipple et al. (2013)). This is expressed as
dz
= k A(x)−θs , (A.1)
dx
where x is the along-channel spatial coordinate (increases upstream), θ is a
constant that sets longitudinal curvature of the channel, A(x) is drainage area
upstream of position x, and ks is an empirical constant. Changes in uplift rate
generate knickpoints that travel up through the channel network until the entire
system re-equilibrates to the new uplift conditions. These “slope-break” knickpoints
(Whipple et al. (2013)) manifest as discontinuities in the relationship between
channel slope and drainage area given by Eq. A.1.
We identify knickpoints in the trunk streams of sixteen Columbia tributaries
using 3 m resolution LiDAR data compiled by the Oregon Department of Geology
and Mineral Industries (https://www.oregongeology.org/lidar/) using the default
flow routing algorithm in MATLAB-based Topotoolbox (Schwanghart and Scherler
(2014)). Each profile is smoothed using a 100m moving window to eliminate
channel features associated with lithologic variation (waterfalls), as well as DEM
error inherent to gridded elevation data sets in steep terrain. We use a coordinate
transformation to remove the effects of drainage area (and concavity) from
130
longitudinal profiles in Eq. A.1 which allows us to write
z(x) = ksχ(x), ∫ (A.2)s θ
where z(x) is elevation relative to base level and χ(x) = Ao θ dx is an area-so A(x)
normalized upstream coordinate (Perron and Royden (2013)). Ao is an arbitrary
constant introduced to give χ units of length. In these coordinates longitudinal
profiles (z(x) versus χ(x)) of equilibrium channels plot as straight lines with slope
dz/dχ = ks, and slope-break knickpoints manifest as changes to linear channel
slopes (Royden and Taylor Perron (2013)). While it is common in the literature to
assign Ao a value of 1×106 m2 (Perron and Royden (2013)) we instead choose Ao =
1 m2 in which case the slope of the χ-space longitudinal profile is equivalent to ks
used in slope-area analysis (Whipple et al. (2013)), allowing us to directly compare
two common approaches for classifying fluvial profiles (Figure A.44). To constrain
θ we perform a grid search over concavity values ranging from 0.3 to 0.6 where for
each θ value we find a piece-wise least squares linear fit of all profiles assuming
a single knickpoint in each. This gives knickpoint locations for all catchments as
well as the value of θ that best linearizes the stream profiles (FigureA.43). We find
θ = 0.46+0.01−0.02 which is consistent with expectation for fluvial erosion (Whipple
(2004)).
We assume a simple ‘stream power’ type erosion law in which fluvial
erosion is assumed proportional to channel slope raised an exponent (n) (Whipple
and Tucker (1999)). Although a range of values have been proposed,it has been
suggested from hydraulic geometry relationships that n ≈ 1 in bedrock channels
(Venditti et al. (2019)). Assuming this linear form of the stream power model,
knickpoint upstream positions in χ-space scale with total cumulative uplift since
the onset of increased uplift rate that initiated the knickpoint (Goren, Fox, and
131
Figure A.43. 16 Columbia River tributaries used in this study. A. Blue
lines are longitudinal profiles extracted from 3m LIDAR data plotted against an
area-normalized χ-scale calculated using best fitting concavity value of θ = 0.455
and Ao = 1 m
2. Red dashed lines are piece-wise linear fits that best capture
slope-break knickpoints with best fitting knickpoint locations marked by blue
squares on each profile. B. PDF of θ values calculated using best-fitting knickpoint
decompositions for each value.
Willett (2014)). Thus for a given uplift rate (U) and onset time (τ) Eq. A.2 can be
written as Uτ = ksχ.
For a given catchment we then expect the χ-position of a transient signal
associated with the onset of hyaloclastite uplift to be χ = w/ks where w is the
magnitude of hyaloclastite uplift. By computing hyaloclastite displacement for each
catchment from our flexural model we can find ks values that best map this uplift
onto knickpoint positions, thus building a coupled model between hyaloclastite
uplift and the transient fluvial response. We find that knickpoint positions can
be fit via this model with ks = 140 ± 17 for all catchments and that error is
minimized when the component of uplift associated with erupted material is not
included in uplift term of the knickpoint model. In other active tectonic settings,
132
numerous studies (Adams, Whipple, Forte, Heimsath, and Hodges (2020); Gallen
and Fernández-Blanco (2021)) vary the slope exponent n to help explain the fluvial
response to uplift gradients. Here, strong physical constraints on spatially varying
uplift across the Gorge reveal that it is possible to fit 16 tributary channel profiles
with fewer parameters using the simpler linear stream power model.
Slope-Area Analysis
The assumption that ks is constant across strongly (factor of 4 − 5) varying
uplift requires justification, as uplift rate is commonly linked to ks in studies
of transient landscape evolution (Kirby and Whipple (2012)). We can test the
robustness of constant ks and θ between tributaries by independently constraining
their values through slope-area analysis (Wobus et al. (2006)). We first calculate
slope and upstream drainage area for every pixel in the LiDAR dataset using
algorithms in Topotoolbox (Schwanghart and Scherler (2014)). We posit the
likelihood of spatially (but not temporally) constant uplift in each channel, which
means that only sections of channel below the knickpoint reflect the fluvial response
to modern uplift rates. We clip off channel reaches above the elevations implied
by the hyaloclastite uplift model for a given channel. We further filter this data
to only include points in the landscape with upstream drainage areas greater than
2× 106 km2.
In this case we do not limit our analysis to knickpoint positions in the larger
trunk streams, but include all channel points within the study area which yields a
much larger data set. The resulting slope-area distribution of fluvial channels below
the lithologic horizon of the hyaloclastite band is fit using the robustfit function in
MATLAB. This gives best fitting values of ks = 132 and θ = 0.45, consistent within
133
uncertainty with our values derived from the regression of longitudinal profiles
outlined above (Figure A.44).
Figure A.44. Slope-area plot for fluvial channels downstream of the
knickpoint. Points represent raw slope-area data for all points within fluvial
channels downstream of paleosurface elevations implied by the flexure model. The
blue line and shaded region represents a model fit to equation A.1 using the range
of ks values from the coupled grid-search inversion and concavity values derived
in the knickpoint fitting process. The red dashed line is a fit to the raw slope-area
data using the robustfit function in Matlab.
A.0.3 Interpretation of Bouguer gravity anomaly. The low in
the Bouguer gravity anomaly collocated with hyaloclastite uplift (FigureA.45A
and Figure 2B of main text) implies a flexural response to a buried load. To show
134
this we interpolate detrended Bouguer gravity data (Bonvalot et al. (2012)) onto
hyaloclastite outcrop positions and calculate the admittance between gravity and
hyaloclastite deformation as the slope of a linear fit between the two datasets
(FigureA.46A). Similar Bouguer gravity lows in other volcanic settings have been
interpreted as indicative of isostatic compensation of a magmatically thickened
crust (Perkins et al. (2016)), and so we compare this value to expectation for
isostatic equilibrium. FigureA.46B shows that the calculated admittance in the
Columbia Gorge is indicative of a flexural response to buried loading (Forsyth
(1985)). While this data informs the development of our structural model inherent
parameter trade-offs in gravity models limit the extent to which we can constrain
the depth and composition of the buried load (Blakely (1994)).
Figure A.45. Bouguer gravity and heat flow data. A. Bouguer gravity
anomaly (Bonvalot et al. (2012)) with locations of quaternary volcanic vents, the
modern Columbia River, and the extent used to generate the cross-arc gravity
profile shown in Figure 39.B. B. Heat flow data from Williams and DeAngelo
(2008) (C. F. Williams and DeAngelo (2008)) interpolated onto a 2-minute grid
showing measurement stations.
135
Figure A.46. Interpretation of admittance between gravity and
hyaloclastite deformation. A. Detrended Bouguer gravity plotted versus
hyaloclastite deformation where the slope of the linear fit between the datasets
is the gravitational admittance over a wavelength of ∼ 80 km. B. Regime diagram
showing expected ranges of admittance from different loading scenarios. The black
line represents isostatic compensation at the base of the crust. Points above the
line represent a flexural response to buried loading while points below the line
represent a flexural response to surface loading.
136
A.0.4 Elastic deformation/knickpoint response model. We build
a forward model that assumes east-west uplift of hyaloclastite deposits through the
Columbia River Gorge can be approximated as deformation of a thin elastic plate
subject to basal pressure from magmatic intrusions. Assuming no east-west lateral
forcing on the plate (as implied by the tectonic field) the thin-plate flexure equation
in 1-d is
d4w
D + ρlwg = q(x). (A.3)
dx4
Here w is deflection of the plate, the second term on the left is a restoring force
that scales with uplift, q(x) is a pressure distribution from magmatic intrusion at
the base of the plate, and D is flexural rigidity defined as
ET 3
D = e . (A.4)
12(1− ν2)
where E is the Young’s Modulus, ν is Poisson’s Ratio, and Te is the elastic
thickness of the plate. The pressure distribution acting on the plate q(x) is modeled
as a Gaussian function
Fnet − (x−x
2
c)
q(x) = √ e 2σ2 , (A.5)
σ 2π
where σ2 is the variance and xc is the center of the pressure distribution. Since
plates deforming within their flexural waveband are insensitive to the shape of the
load (Watts (2001)) we assign a constant variance of σ2 = 25 km and only vary
the total force Fnet acting on the plate. We solve Eq. A.3 in the Fourier domain
(Sandwell (1984)) defining the Fourier∫transform pair∞
f(x) = F (s)e2πisxds (A.6)
∫ −∞∞
F (s) = f(x)e−2πisxdx. (A.7)
−∞
137
Taking the Fourier transform gives an expression for displacement of the plate as a
function of wavenumber for the pressure distribution Q(s)
Q(s)
W (s) = , (A.8)
(2πs)4D + ρlg
which can be inverted to find w(x).
The mean hyaloclastite uplift value within each catchment is then mapped
onto knickpoint locations using equation A.2, rearranged such that χi = wi/ks for
each basin (where i = 1 : 16).
With this forward model we perform a grid search over plausible values of
flexural rigidity (5 × 1020Nm≤ D ≤ 12 × 1020 Nm), total force on the plate (31 ×
106N≤ F 6net ≤ 38×10 N), and the fluvial proportionality constant (110m ≤ ks ≤ 170
m). This we pose as a Bayesian inverse problem with a solution of the form
σ(m) = kα(m)L(m) (A.9)
which states that the posterior probability density function σ(m) for a parameter
distribution m is the product of the a priori probability distribution (α(m)), a
likelihood function L(m) and a normalization constant k. We assume uniform
priors and so the posterior probability distribution is proportional to the likelihood
function defined as ∏n [ ]1
L(m) = exp − (gi(m)− d )Ti Σ−1i (g2 i
(m)− di) k (A.10)
i=1
where d is our data vector, g(m) is the output of a forward model for a given
input parameter set m, Σ−1 is the determinate of the covariance matrix storing
data uncertainty values (Mosegaard and Tarantola (1995); Tarantola (2004)). In
our model d = [d Thy,dkn] where dhy consists of mean hyaloclastite elevations
in 20 km bins, and dkn contains the best-fitting knickpoint locations for the 16
138
fluvial tributaries. The covariance matrix is a diagonal matrix in which uncertainty
in each data point in the data vector d is stored along the diagonal. For dhy the
covariance is approximated as the full range about the mean in each of the bins,
while for dkn we calculate a 95% confidence interval for knickpoint locations given
our best fitting value of the area constant (θ = 0.46). The PDFs generated
through this inversion approach (FigureA.47) define a family of models capable
of explaining hyaloclastite uplift in combination with the transient topographic
response observed in fluvial tributaries.
With flexural rigidity (D) constrained by the above maximum likelihood
estimation we calculate effective elastic thickness Te of the deformed upper crust by
rearranging equation A.4 to ( ) 1
12D(1− ν2) 3
Te = . (A.11)
E
FigureA.48 shows the results of varying E from (1 − 8) × 1020Pa for possible
lithologies underlying Gorge topography Turcotte and Schubert (2014) and vary
Poisson’s Ratio ν between 0.1− 0.3.
139
Figure A.47. Posterior PDFs and covariance plots from maximum
likelihood estimation. D is flexural rigidity, Fnet is the total force acting on
the plate, and ks is the fluvial proportionality constant that maps hyaloclastite
uplift onto knickpoint positions. Orange bars represent 95% confidence intervals for
parameter values.
140
Figure A.48. Possible values for effective elastic thickness (Te) as a function
of flexural rigidity (D), Young’s Modulus (E) and Poisson’s Ratio (ν). A. Variation
in Te as a function of E. The solid black line shows Te for our best-fitting value
of flexural rigidity (D) and the dashed black lines represent the 95% confidence
interval. The blue and red dashed lines show variation in Te that result from
varying Poisson’s Ratio where blue and red represent ν = 0.3 and ν = 0.1
respectively. B. Variation in Te for all values of D explored in the maximum
likelihood estimation, where the solid black line represents the best-fitting model
and the dashed lines are 95% confidence intervals.
A.0.5 Isostatic response to erosion by the Columbia River. We
estimate the volume of material removed by erosion since the start of hyaloclastite
uplift by projecting our best fitting flexure model north and south across the
141
gorge and differencing with modern topography (FigureA.49A). This gives a rough
estimate for the full isostatic load where negative values represent the amount of
material removed by erosion (FigureA.49B). We calculate the amount of uplift
caused by removal of this material (FigureA.49C and FigureA.49D) using a 2D
thin-plate flexure equation (Watts (2001)). There are additional loads pushing back
down on the plate (erupted material) that will offset some of this uplift, however
since topography overlying the projected flexure model includes structures not
related to magmatism we do not attempt to account for this type of loading in our
force balance and instead take this to be an upper bound on isostatic uplift since
deposition of the hyaloclastites.
142
Figure A.49. Isostatic uplift from post 3.5 Ma incision implied by our
uplift model. A. Elevation difference between topography and best fitting flexure
model projected north and south along the arc. B. Erosion implied by north-south
model projection. C. Maximum flexural uplift from erosional unloading shown in
panel B. D. Range of model fits to hyaloclastite data and range of values for uplift
caused by erosion within the north-south range of hyaloclastite deposits (shown by
dashed black lines).
Seismic tomography
To derive the relationship between Vs and melt fraction, we approximate the
mid-lower crustal magmatic system as a complex of crystals and melt, and apply
143
the effective-medium theory (Berryman (1980)) to estimate its elastic properties,
including compressional wave speed (V p), shear wave speed (V s) and density
(ρ). This follows a two-step approach adopted by a previous study (Maguire et
al. (2022)).
The first step is to estimate the elastic properties of the solid and melt phase
separately. For the solid phase, we use the PerpleX package (Connolly (2009)) to
compute V p, V s, and ρ considering both basaltic and dacitic bulk compositions,
which represent the two major types of magmas erupted at the Mount St. Helens
(Wanke, Clynne, von Quadt, Vennemann, and Bachmann (2019); Wanke, Karakas,
and Bachmann (2019)). We use the oxidates reported by Blatter et al., (2017)
(Blatter et al. (2017)) for the dacitic composition and those from Leeman et
al., (2005) (Leeman, Lewis, Evarts, Conrey, and Streck (2005)) for the basaltic
composition. Table A.1 summarises the compositions of major oxidates used in
this study. For the PerpleX calculations, we account for a relatively large range
of temperature conditions between 750°C – 1000°C; while we vary the pressure
condition in a normal range of 700 MPa – 900 MPa for such depths.
Though some higher melting temperatures are reported from petrological
analysis in the mid-lower crust of the region, such as Mount St. Helens (Wanke,
Clynne, et al. (2019)), we are considering the average temperature condition for a
15 km thick column. Figure A.50 shows the variations of V p, V s, and ρ relative
to the change of pressure (P) and temperature (T) condition, respectively, which
outlines several regimes where slightly different mineral assemblages are stable (see
Table A.2 for details). The boundaries between the domains are associated with
mineral change and they together define a complex relationship between the elastic
parameters and the P-T condition.
144
To approximate the range of each elastic property of the solid phase, we
use the corresponding values of each elastic parameter at the conditions where the
minimum and maximum values of the sum of V p, V s, and ρ are achieved. The
preference for using the sum of the three parameters over the individual parameter
space is to consider their inter-connected nature. The red and yellow filled stars
in Figure A.50 mark the maximum and minimum values of each parameter,
respectively. For the melt phase, we follow the approach of Ueki and Iwamori
(2016) (Ueki and Iwamori (2016)) to estimate the properties of V p and ρ (with its
Vs being zero) at the same P-T conditions where maximum and minimum values
of elastic parameters are achieved (marked by the red and yellow stars in Figure
A.50).
The method of Ueki and Iwamori (2016) is effective in deriving the
relationships between pressures, temperatures, and H2O concentrations of various
hydrous melts from ultramafic to felsic compositions at pressures of 0–4.29 GPa. To
account for the generally hydrous conditions of the erupted magmas (Blatter et al.
(2017)), 5wt% H2O is added to the reported oxidates before they are normalized
for the melt phase calculation. Increasing the H2O will generally make the melt
phase more hydrous, thus reducing the Vp and density of the melt phase, the Vs
of the two-phase aggregate as well as the volume of melt. However, the resulted
change is minor. The melt-phase calculation is made for both basaltic and dacitic
compositions.
In the second step, we use the estimated elastic properties of the solid and
melt phase to compute the velocity of a two-phase aggregate by initially assuming
air-filled pore spaces. We then correct the velocity by substituting it in the liquid
phase following Gassmann’s equation (Gassmann (1951)). The calculations are
145
made via the open-source package of ElasticC (Paulatto et al. (2022)), which is
available at http://github.com/michpaulatto/ElasticC. Note that our approach in
this step yields the low frequency or relaxed seismic velocity of the partially molten
aggregate.
In the above procedures, one key parameter that will significantly
influence the melt volume estimation later is the aspect ratio, which describes the
geometrical organization of melt relative to a usually assumed elliptical pore space
along grain boundaries (Holness (2006)). Smaller aspect ratios mean less melt is
required to wet grain boundaries compared to more equant pockets isolated near
grain boundary junctions. Maguire et al. (2022) report that the predicted melt
volume for the shallow magmatic reservoir beneath Yellowstone can be up to 4
times more if using an aspect ratio of 1 compared to that from a value of 0.05.
In our calculation, we also vary the aspect ratios between 0.05 and 1 and found
that a similar conclusion holds for our case (Figure A.51). However, theoretical
work (Takei (2002)) indicates that texturally equilibrated partially molten rocks
have aspect ratios usually in the range between 0.1 – 0.15, which is equivalent to
a dihedron angle range of 20◦ − 40◦. We think this could be the case for mid-
lower crustal reservoirs beneath the Gorge, where stored magma is likely not
directly linked to active edifices. We note however that uncertainties exist on
this assumption, and it is also difficult to ascribe one aspect ratio for the system
considering the relatively large geographic distribution of the mid-lower crustal
reservoir (Figure A.52) and potentially complex melt distribution in general.
To estimate the volume of active melts present in the mid-lower crust
currently beneath the Gorge, we use the 3.55 km/s velocity contour as a diagnostic
of the presence of partial melt. This contour corresponds to a melt fraction of 3 –
146
6% for basaltic composition, and slightly lower, 1.5 – 4% for dacitic compositions.
To construct a velocity image for the melt estimation, we average the Vs along an
E-W corridor centered around the Gorge area but with an extension of 0.2° to the
north and south, a width that is close to the resolution of our seismic tomography
model at these depths. The Vs distributions in the mid-lower crust (e.g., 20-35
km) show generally symmetric features within this corridor around the Gorge
area (Figure A.52), justifying the process of averaging. Assuming an average melt
fraction of 6% implies a total volume of active melts of ∼ 3.7 km3 per km arc
distance while an average melt fraction of 1.5% melt indicates a lower volume of
∼ 1 km3 per km arc distance.
Oxide Weight % for basalts Weight % for dacites
SiO2 50.40 65.37
TiO2 1.98 0.54
Al2O3 16.28 17.9
FeO 10.55 3.84
MnO 0.15 0.08
MgO 6.84 1.52
CaO 8.55 4.63
Na2O 3.91 4.78
K2O 1.34 1.29
Table A.1. Composition of typical basalts and dacites at Mount St.
Helens used to calculate seismic velocity. The basaltic composition is from
Blatter et al., 2017, and dacitic composition is from Leeman et al., (2005).
A.0.6 Heat flow. Quantifying magmatic origins of heat flow
anomalies in the central Cascades is complicated by strong near-surface meteoric
circulation (Saar and Manga (2004)). This is reflected by non-conductive
geothermal gradients at the arc axis reflecting lateral advection of heat (Ingebritsen
et al. (1989)). The regional compilation of (Ingebritsen and Mariner (2010)) shows
localized high heat flow coincident with Quaternary volcanic vents throughout the
147
Figure A.50. Predicted Vp, Vs, and ρ with respect to the variations of
pressure and temperature condition for the solid phase using PerpleX,
both for basalts (A-C) and dacites (D-E). The red and yellow stars define
the maximum and minimum value of each parameter to define a range used in
the calculation. The black lines in B and E separate the 2D space into several
domains where slightly different mineral assemblages are stable. See Table S2 for
the detailed mineral assemblages for each domain.
148
Figure A.51. Modelled Vs and melt fraction relationship for basalts (A)
and dacites (B). The color of the curves represents the different aspect ratios
assumed during the calculation, with the textural equilibrium regime of aspects
ratio between 0.1-0.15 highlighted in yellow. The solid lines are from calculations
using the maximum values of elastic parameters (red stars in Figure A.50) and the
dashed lines are from the minimum values. The grey-shaded region indicates Vs
range observed in this study.
149
Figure A.52. Regional maps of Vs at 20 km (A), 25 km (B), 30 km (C),
and 35 km (D) depths, respectively. The red-filled stars denote the volcanoes
of Mount St. Helens (MSH), Mount Adams (MA), and Mount Hood (MH). The
grey-shaded region highlights the corridor where Vs values are averaged to produce
the transect shown in Figure 39C. The color of the curves represents the different
aspect ratios assumed during the calculation, with the textural equilibrium regime
of aspects ratio between 0.1-0.15 highlighted in yellow. The solid lines are from
calculations using the maximum values of elastic parameters (red stars in Figure
A.50) and the dashed lines are from the minimum values. The grey-shaded region
indicates Vs range observed in this study.
150
Domain Basalts Dacites
almandine, spessartine, enstatite, almandine, spessartine, kyanite,
1 diopside, jadeite, albite, sphene, diopside, albite,
sanidine, ilmenite, hercynite sanidine, anorthite, quartz
forsterite, almandine, spessartine, almandine, spessartine, kyanite,
2 diopside, jadeite, albite, diopside, albite, sanidine,
sanidine, ilmenite, hercynite anorthite, quartz, ilmenite
forsterite, fayalite, almandine, almandine,spessartine, enstatite,
3 spessartine, diopside, albite, diopside, albite, sanidine,
sanidine, ilmenite, hercynite anorthite, quartz, ilmenite
fayalite, almandine, spessartine,
4 enstatite, diopside, albite,
sanidine, ilmenite, hercynite
Table A.2. The mineral assemblages stabilized at different P-T conditions.
The number and corresponding domains are specified in Figure A.50 panels B and
E.
arc (with increasing correlation for Holocene age vents, (O’Hara et al. (2020))).
Interpolated heat flow variations around the Gorge (C. F. Williams and DeAngelo
(2008)) are displayed in Figure A.45B and Figure 39B.
Cascades arc-centered heat flow anomalies likely arise from magmatic
intrusions at depth (Blackwell, Steele, Kelley, and Korosec (1990)). Here we show
that the magnitude of peak heat flow is consistent with long-lived magma flux into
the lower/mid crust at depths spanning observed low seismic Vs velocities. While
a joint inversion of geophysical data that accounts for fluid advection of heat is
needed to capture the details of regional heat flow, a simple 1D model similar
to that used in previous work (Till et al. (2019)) is sufficient to argue that the
presence of anomalous heat at the arc axis is likely due to intrusions at depth, and
consistent with the seismically inferred melt at present day.
We estimate Cascades geothermal gradient via solution to
d2T (z)
kT = P (z), (A.12)
dz2
151
with T temperature as a function of depth z, subject to T (0) = 13 C (yearly
average for Portland, OR) and ksdT/dz|z=0 = qsurf , kT is the effective conductivity
in the near-surface. We vary as through the range of heatflow observed in the
gorge, qsurf ∈ [45, 85] mW/m2. P (z) is a radiogenic heat production term, of the
generic form
P (z) = ρH0 exp (−z/λ), (A.13)
with H0 the radiogenic heat production and λ a scale length for concentration in
the shallow crust (Turcotte and Schubert (2014)), with ρ an average crustal density.
We experiment with and without this term to demonstrate its effects qualitatively.
Figure SA.53 shows several representative conductive geotherms that
demonstrate the likelihood of melt as inferred from seismic tomography.
152
0
with radiogenic heat
without radiogenic 
5 heat
q s u  r f= 45 mW/m 2
q s u  r f= 55 mW/m 2
q s u  r f= 65 mW/m 2
10 q s u  r f= 75 mW/m 2
q s u  r f= 85 mW/m 2
Mt. St. Helens dacite
solidus, ρ=2800 kg/m 3
(Blatter et al., 2017)
15
20
25
30
0 200 400 600 800 1000 1200
Temperature (C)
Figure A.53. Illustrative conductive geotherms with and without radiogenic
heat production (H0 = 9 × 10−10 W/kg). Curve colors correspond to assumed
surface heat flow, spanning the range observed in the Gorge. The black curve is
Dacite solidus with 5 weight percent water (Blatter et al. (2017)), computed in
PerpleX. The curves are lower bounds for geothermal gradient as they neglect
heat associated with intrusions, and hydrothermal circulation in the shallow crust
(Ingebritsen et al. (1989)), supporting our seismic inference of lower-crustal melt at
present.
153
Depth (km)
REFERENCES CITED
Adams, B. A., Whipple, K. X., Forte, A. M., Heimsath, A. M., & Hodges, K. V.
(2020). Climate controls on erosion in tectonically active landscapes. Science
Advances , 6 (42), eaaz3166. doi: 10.1126/sciadv.aaz3166
Audet, P., Jellinek, A. M., & Uno, H. (2007). Mechanical controls on the
deformation of continents at convergent margins. Earth and Planetary
Science Letters , 264 (1-2), 151–166. doi: 10.1016/j.epsl.2007.09.024
Bagnold, R. A. (1966). An Approach to the Sediment Transport Problem from
General Physics. USGS Professional Paper , 42. doi:
10.1017/S0016756800049074
Barnhart, K. R., Tucker, G. E., Doty, S. G., Shobe, C. M., Glade, R. C., Rossi,
M. W., & Hill, M. C. (2020). Inverting Topography for Landscape Evolution
Model Process Representation: 3. Determining Parameter Ranges for Select
Mature Geomorphic Transport Laws and Connecting Changes in Fluvial
Erodibility to Changes in Climate. Journal of Geophysical Research: Earth
Surface, 125 (7), 1–34. doi: 10.1029/2019JF005287
Bater, C. W., & Coops, N. C. (2009). Evaluating error associated with
lidar-derived DEM interpolation. Computers and Geosciences , 35 (2),
289–300. doi: 10.1016/j.cageo.2008.09.001
Baynes, E. R., Lague, D., Attal, M., Gangloff, A., Kirstein, L. A., & Dugmore,
A. J. (2018). River self-organisation inhibits discharge control on waterfall
migration. Scientific Reports , 8 (1), 1–8. doi: 10.1038/s41598-018-20767-6
Bedrosian, P. A., Peacock, J. R., Bowles-Martinez, E., Schultz, A., & Hill, G. J.
(2018). Crustal inheritance and a top-down control on arc magmatism at
Mount St Helens. Nature Geoscience, 11 (11), 865–870. Retrieved from
http://dx.doi.org/10.1038/s41561-018-0217-2 doi:
10.1038/s41561-018-0217-2
Beeson, M. H., & Tolan, T. L. (1990). The Columbia River Basalt Group in the
Cascade Range: A Middle Miocene Reference Datum for Structural
Analysis. Journal of Geophysical Research, 95 , 19547–19559.
Bergbauer, S., & Pollard, D. D. (2003). How to calculate normal curvatures of
sampled geological surfaces. Journal of Structural Geology , 25 (2), 277–289.
doi: 10.1016/S0191-8141(02)00019-6
154
Berlin, M. M., & Anderson, R. S. (2007). Modeling of knickpoint retreat on the
Roan Plateau, western Colorado. Journal of Geophysical Research: Earth
Surface, 112 (3), 1–16. doi: 10.1029/2006JF000553
Berryman, J. G. (1980). Long-wavelength propagation in composite elastic media
ii. ellipsoidal inclusions. The Journal of the Acoustical Society of America,
68 (6), 1820–1831.
Blackwell, D. D., Steele, J. L., Kelley, S., & Korosec, M. A. (1990). Heat flow in
the State of Washington and thermal conditions in the Cascade Range.
Journal of Geophysical Research, 95 (B12). doi: 10.1029/jb095ib12p19495
Blakely, R. J. (1994). Extent of partial melting beneath the Cascade Range ,
Oregon : Constraints from gravity anomalies and ideal-body theory. Journal
of Geophysical Research, 99 , 2757–2773.
Blatter, D. L., Sisson, T. W., & Hankins, W. B. (2017). Voluminous arc dacites as
amphibole reaction-boundary liquids. Contributions to Mineralogy and
Petrology , 172 , 1–37.
Bonvalot, S., Balmino, G., Brias, A., Kuhn, M., Peyrefitte, A., Vales, N., . . .
Sarrailh, M. (2012). World Gravity Map (Vol. 1; Tech. Rep.).
Booth, A. M., Roering, J. J., & Perron, J. T. (2009). Automated landslide
mapping using spectral analysis and high-resolution topographic data: Puget
Sound lowlands, Washington, and Portland Hills, Oregon. Geomorphology ,
109 (3-4), 132–147. Retrieved from
http://dx.doi.org/10.1016/j.geomorph.2009.02.027 doi:
10.1016/j.geomorph.2009.02.027
Booth, A. M., Roering, J. J., & Rempel, A. W. (2013). Topographic signatures and
a general transport law for deep-seated landslides in a landscape evolution
model. Journal of Geophysical Research: Earth Surface, 118 (2), 603–624.
doi: 10.1002/jgrf.20051
Braun, J., & Willett, S. D. (2013). A very efficient O(n), implicit and parallel
method to solve the stream power equation governing fluvial incision and
landscape evolution. Geomorphology , 180-181 , 170–179. Retrieved from
http://dx.doi.org/10.1016/j.geomorph.2012.10.008 doi:
10.1016/j.geomorph.2012.10.008
Bretz, J. H. (1917). The Satsop Formation of Oregon and Washington. The
Journal of Geology , 25 (5), 446–458. doi: 10.1086/622509
155
Bui, L. K., & Glennie, C. L. (2023). Estimation of lidar-based gridded DEM
uncertainty with varying terrain roughness and point density. ISPRS Open
Journal of Photogrammetry and Remote Sensing , 7 (December 2022),
100028. Retrieved from https://doi.org/10.1016/j.ophoto.2022.100028
doi: 10.1016/j.ophoto.2022.100028
Burov, E. B., & Diament, M. (1995). The effective elastic thickness (Te) of
continental lithosphere: what does it really mean? Journal of Geophysical
Research, 100 (B3), 3905–3927. doi: 10.1029/94JB02770
Cheraghi, M., Rinaldo, A., Sander, G. C., Perona, P., & Barry, D. A. (2018).
Catchment Drainage Network Scaling Laws Found Experimentally in
Overland Flow Morphologies. Geophysical Research Letters , 45 (18),
9614–9622. doi: 10.1029/2018GL078351
Claps, P., & Oliveto, G. (1996). I q = l / 3. , 32 (10), 3123–3135.
Cogliati, A., & R, R. (2022). The Origins of the Fundamental Theorem of Surface
Theory. HIstoria Mathematica(61), 45–79.
Connolly, J. (2009). The geodynamic equation of state: what and how.
Geochemistry, geophysics, geosystems , 10 (10).
Conrey, R., Sherrod, D., Uto, K., & Uchiumi, S. (1996). Potassium-Argon Ages
from Mount Hood Area of Cascade Range, Northern Oregon.
Isochron/West , 63 (63), 10–20.
Conrey, R., Sherrod, D. R., Hooper, P. R., & Swanson, D. A. (1997). Diverse
primitive magmas in the Cascade arc, Northern Oregon and Southern
Washington. Canadian Mineralogist , 35 (2), 367–396.
Conrey, R. M., Taylor, E. M., Donnelly-Nolan, J. M., & Sherrod, D. R. (2002).
North-central Oregon Cascades: Exploring petrologic and tectonic intimacy
in a propagating intra-arc rift. Field guide to geologic processes in Cascadia,
36 (36), 47–90. Retrieved from
http://www.academia.edu/download/41265842/North-central Oregon
Cascades Exploring 20160116-15874-1mrpyao.pdf%0Apapers3://
publication/uuid/F680CAC7-FB0A-4F9B-BC6A-8A4212B606AA
Davis, A. W. M. (1892). The Convex Profile of Bad-Land Divides. Science,
20 (508), 27–28.
Dekking, M. (2010). A modern introduction to probability and statistics:
understanding why and how. London: Springer.
156
Deligne, N. I., Mckay, D., Conrey, R. M., Grant, G. E., Johnson, E. R., O’Connor,
J., & Sweeney, K. (2017). Field-trip guide to mafic volcanism of the Cascade
Range in Central Oregon—A volcanic, tectonic, hydrologic, and geomorphic
journey. U.S. Geological Survey Scientific Investigations Report , 110.
Retrieved from http://pubs.er.usgs.gov/publication/sir20175022H
doi: 10.3133/sir20175022H
Deshpande, N. S., Furbish, D. J., Arratia, P. E., & Jerolmack, D. J. (2021). The
perpetual fragility of creeping hillslopes. Nature Communications , 12 (1),
1–7. Retrieved from http://dx.doi.org/10.1038/s41467-021-23979-z
doi: 10.1038/s41467-021-23979-z
Dickenson, W., & Snyder, W. (1979). Geometry of Subducted Slabs Related to San
Andreas Transform Author ( s ): William R . Dickinson and Walter S .
Snyder Source : The Journal of Geology , Vol . 87 , No . 6 ( Nov ., 1979 ),
pp . 609-627 Published by : The University of Chicago Press Stable.
Geology , 87 (6), 609–627.
du Bray, E. A., & John, D. A. (2011). Petrologic, tectonic, and metallogenic
evolution of the Ancestral Cascades magmatic arc, Washington, Oregon, and
northern California. Geosphere, 7 (5), 1102–1133. doi: 10.1130/GES00669.1
Einstein, A. (1982). The Special and General Theory of Relativity. Man and the
Universe, 1 (series , ??–??
Evarts, R., Conrey, R. M., Fleck, R. J., & Hagstrum, J. (2009). Volcanoes to
Vineyards: Geologic Field Trips through the Dynamic Landscape of the
Pacifi c Northwest. Volcanoes to Vineyards: Geologic Field Trips through the
Dynamic Landscape of the Pacific Northwest: Geological Society of America
Field Guide 15 , 015 (32), 253–270. Retrieved from
http://geoprisms.org/~geo/images/stories/documents/Alaska/
Boring FT guide%28GSA2009%29.pdf doi: 10.1130/2009.fl
Fenneman, N. (1908). Some Features of Erosion by Unconcentrated Wash. The
Journal of Geology , 16 (8), 746–754.
Fleck, R. J., Hagstrum, J. T., Calvert, A. T., Evarts, R. C., & Conrey, R. M.
(2014). 40Ar/39Ar geochronology, paleomagnetism, and evolution of the
boring volcanic field, Oregon and Washington, USA. Geosphere, 10 (6),
1283–1314. doi: 10.1130/GES00985.1
Flint, J. J. (1974). Stream gradient as a function of order, magnitude, and
discharge. Water Resources Research, 10 (5), 969–973. doi:
10.1029/WR010i005p00969
157
Flores-Cervantes, J. H., Istanbulluoglu, E., & Bras, R. L. (2006). Development of
gullies on the landscape: A model of headcut retreat resulting from plunge
pool erosion. Journal of Geophysical Research: Earth Surface, 111 (1), 1–14.
doi: 10.1029/2004JF000226
Forsyth, D. W. (1985). Subsurface loading and estimates of the flexural rigidity of
continental lithosphere. Journal of Geophysical Research, 90 (B14),
12623–12632. doi: 10.1029/jb090ib14p12623
Furbish, D. J., Haff, P. K., Dietrich, W. E., & Heimsath, A. M. (2009). Statistical
description of slope-dependent soil transport and the diffusion-like
coefficient. Journal of Geophysical Research: Earth Surface, 114 (4), 1–19.
doi: 10.1029/2009JF001267
Gabet, E. J. (2003). Sediment transport by dry ravel. Journal of Geophysical
Research: Solid Earth, 108 (B1), 1–8. doi: 10.1029/2001jb001686
Gallen, S. F., & Fernández-Blanco, D. (2021). A New Data-Driven Bayesian
Inversion of Fluvial Topography Clarifies the Tectonic History of the Corinth
Rift and Reveals a Channel Steepness Threshold. Journal of Geophysical
Research: Earth Surface, 126 (3), 1–29. doi: 10.1029/2020JF005651
Gannett, M., Lite, K., Morgan, D., & Collins, C. (2001). Ground-Water Hydrology
of the Upper Deschutes Basin, Oregon; Water-Resources Investigations
Report 00–4162. , 72.
Garcia, E. S., Sandwell, D. T., & Luttrell, K. M. (2015). An iterative spectral
solution method for thin elastic plate flexure with variable rigidity.
Geophysical Journal International , 200 (2), 1012–1028. doi:
10.1093/gji/ggu449
Garcia-Castellanos, D. (2002). Interplay between lithospheric flexure and river
transport in foreland basins. Basin Research, 14 (2), 89–104. doi:
10.1046/j.1365-2117.2002.00174.x
Gassmann, F. (1951). Elastic waves through a packing of spheres. Geophysics ,
16 (4), 673–685.
Gauss, C. F. (1902). General Investigations of Curved Surfaces of 1827 and 1825.
Nature, 66 (1709), 316–317. doi: 10.1038/066316b0
Gayer, E., Mukhopadhyay, S., & Meade, B. J. (2008). Spatial variability of erosion
rates inferred from the frequency distribution of cosmogenic 3He in olivines
from Hawaiian river sediments. Earth and Planetary Science Letters ,
266 (3-4), 303–315. doi: 10.1016/j.epsl.2007.11.019
158
Gilbert, G. K. (1877). Geology of the Henry Mountains. U.S. Geographical and
Geological Survey of the Rocky Mountain Region, 196.
Gilbert, G. K. (1908). The Convexity of Hilltops. The Journal of Geology , 17 (4),
344–350.
Goren, L., Fox, M., & Willett, S. D. (2014). Tectonics from fluvial topography
using formal linear inversion: Theory and applications to the Inyo
Mountains, California. Journal of Geophysical Research F: Earth Surface,
119 (8), 1651–1681. doi: 10.1002/2014JF003079
Gregory-Wodzicki, K. M. (2000). Uplift history of the Central and Northern Andes:
A review. Bulletin of the Geological Society of America, 112 (7), 1091–1105.
doi: 10.1130/0016-7606(2000)112⟨1091:UHOTCA⟩2.0.CO;2
Hack, J. T., Seaton, F. A., & Nolan, T. B. (1957). Studies of Longitudinal Stream
Profiles in Virginia and Maryland. , 45–97.
Harris, F. J. (1978). On the Use of Windows for Harmonic Analysis with the
Discrete Fourier Transform. Proceedings of the IEEE , 66 (1), 51–83. doi:
10.1109/PROC.1978.10837
Hieronymus, C. F., & Bercovici, D. (2001). A theoretical model of hotspot
volcanism: Control on volcanic spacing and patterns via magma dynamics
and lithospheric stresses. Journal of Geophysical Research, 106 (B1),
683–702.
Hildreth, W. (2007). Quaternary magmatism in the Cascades - Geologic
perspectives. US Geological Survey Professional Paper(1744). doi:
10.3133/pp1744
Hildreth, W., & Wilson, C. J. (2007). Compositional zoning of the bishop tuff.
Journal of Petrology , 48 (5), 951–999. doi: 10.1093/petrology/egm007
Hilley, G. E., Brodsky, E. E., Roman, D. C., Shillington, D., Brudzinski, M., Behn,
M. D., & Tobin, H. (2022). SZ4D Implementation Plan. Stanford Digital
Repository .
Hilley, G. E., Porder, S., Aron, F., Baden, C. W., Johnstone, S. A., Liu, F., . . .
Young, H. H. (2019). Earth’s topographic relief potentially limited by an
upper bound on channel steepness. Nature Geoscience, 12 (10), 828–832.
Retrieved from http://dx.doi.org/10.1038/s41561-019-0442-3 doi:
10.1038/s41561-019-0442-3
159
Hippe, K., Jansen, J. D., Skov, D. S., Lupker, M., Ivy-ochs, S., Kober, F., . . .
Egholm, D. L. (n.d.). Cosmogenic in situ 14 C- 10 Be reveals abrupt Late
Holocene soil loss in the Andean Altiplano. Nature Communications(2021).
Retrieved from http://dx.doi.org/10.1038/s41467-021-22825-6 doi:
10.1038/s41467-021-22825-6
Holness, M. (2006). Melt–solid dihedral angles of common minerals in natural
rocks. Journal of Petrology , 47 (4), 791–800.
Howard, A. D., & Kerby, G. (1983). Channel changes in badlands. Geological
Society of America Bulletin, 94 (6), 739–752. doi:
10.1130/0016-7606(1983)94⟨739:CCIB⟩2.0.CO;2
Hurst, A. A., Anderson, R. S., & Crimaldi, J. P. (2021). Toward Entrainment
Thresholds in Fluvial Plucking. Journal of Geophysical Research: Earth
Surface, 126 (5), 1–19. doi: 10.1029/2020jf005944
Hurst, M. D., Mudd, S. M., Walcott, R., Attal, M., & Yoo, K. (2012). Using hilltop
curvature to derive the spatial distribution of erosion rates. Journal of
Geophysical Research: Earth Surface, 117 (2), 1–19. doi:
10.1029/2011JF002057
Ingebritsen, S. E., & Mariner, R. H. (2010). Hydrothermal heat discharge in the
Cascade Range, northwestern United States. Journal of Volcanology and
Geothermal Research, 196 (3-4), 208–218. Retrieved from
http://dx.doi.org/10.1016/j.jvolgeores.2010.07.023 doi:
10.1016/j.jvolgeores.2010.07.023
Ingebritsen, S. E., Mariner, R. H., & Sherrod, D. R. (1994). Hydrothermal systems
of the Cascade Range, north-central Oregon (Vol. 1044-L; Tech. Rep.). doi:
10.3133/pp1044l
Ingebritsen, S. E., Sherrod, D. R., & Mariner, R. H. (1989). Heat flow and
hydrothermal circulation in the cascade range, north-central Oregon.
Science, 243 (4897), 1458–1462. doi: 10.1126/science.243.4897.1458
Jefferson, A., Grant, G. E., Lewis, S. L., & Lancaster, S. T. (2010). Coevolution of
hydrology and topography on a basalt landscape in the Oregon Cascade
Range, USA. Earth Surface Processes and Landforms , 35 (7), 803–816. doi:
10.1002/esp.1976
Jiang, C., Schmandt, B., Abers, G. A., Kiser, E., & Miller, S. (2023). Segmentation
and radial anisotropy of the deep crustal magmatic system beneath the
Cascades arc. Geochemistry, Geophysics, Geosystems , 24 . doi:
10.1029/2022GC010738
160
Karlstrom, L., Dufek, J., & Manga, M. (2009). Organization of volcanic plumbing
through magmatic lensing by magma chambers and volcanic loads. Journal
of Geophysical Research: Solid Earth, 114 (10), 1–16. doi:
10.1029/2009JB006339
Karlstrom, L., Lee, C. T. A., & Manga, M. (2014). The role of magmatically driven
lithospheric thickening on arc front migration. Geochemistry, Geophysics,
Geosystems Research, 2655–2675. doi: 10.1002/2014GC005355.Received
Karlstrom, L., Paterson, S. R., & Jellinek, A. M. (2017). A reverse energy cascade
for crustal magma transport. Nature Geoscience, 10 (8), 604–608. doi:
10.1038/NGEO2982
Kirby, E., & Whipple, K. X. (2012). Expression of active tectonics in erosional
landscapes. Journal of Structural Geology , 44 , 54–75. Retrieved from
http://dx.doi.org/10.1016/j.jsg.2012.07.009 doi:
10.1016/j.jsg.2012.07.009
Korosec, M. A. (1987). Geologic map of the Hood River quadrangle, Washington
and Oregon. Olympia, Wash: Washington Dept. of Natural Resources,
Division of Geology and Earth Resources.
Lai, L. S. H., Roering, J. J., Finnegan, N. J., Dorsey, R. J., & Yen, J. Y. (2021).
Coarse sediment supply sets the slope of bedrock channels in rapidly
uplifting terrain: Field and topographic evidence from eastern Taiwan.
Earth Surface Processes and Landforms , 46 (13), 2671–2689. doi:
10.1002/esp.5200
Lamb, M. P., & Dietrich, W. E. (2009). The persistence of waterfalls in fractured
rock. Bulletin of the Geological Society of America, 121 (7-8), 1123–1134.
doi: 10.1130/B26482.1
Lamb, M. P., Scheingross, J. S., Amidon, W. H., Swanson, E., & Limaye, A.
(2011). A model for fire-induced sediment yield by dry ravel in steep
landscapes. Journal of Geophysical Research: Earth Surface, 116 (3), 1–13.
doi: 10.1029/2010JF001878
Lee, C. T. A., Thurner, S., Paterson, S., & Cao, W. (2015). The rise and fall of
continental arcs: Interplays between magmatism, uplift, weathering, and
climate. Earth and Planetary Science Letters , 425 , 105–119. Retrieved from
http://dx.doi.org/10.1016/j.epsl.2015.05.045 doi:
10.1016/j.epsl.2015.05.045
Leeman, W. P., Lewis, J. F., Evarts, R. C., Conrey, R. M., & Streck, M. J. (2005).
Petrologic constraints on the thermal structure of the cascades arc. Journal
of Volcanology and Geothermal Research, 140 (1-3), 67–105.
161
Lerner, A. H., O’Hara, D., Karlstrom, L., Ebmeier, S. K., Anderson, K. R., &
Hurwitz, S. (2020). The Prevalence and Significance of Offset Magma
Reservoirs at Arc Volcanoes. Geophysical Research Letters , 47 (14). doi:
10.1029/2020GL087856
Li, T., Venditti, J. G., Rennie, C. D., & Nelson, P. A. (2022). Bed and Bank Stress
Partitioning in Bedrock Rivers. Journal of Geophysical Research: Earth
Surface, 1–17. doi: 10.1029/2021jf006360
Lohse, K. A., & Dietrich, W. E. (2005). Contrasting effects of soil development on
hydrological properties and flow paths. Water Resources Research, 41 (12),
1–17. doi: 10.1029/2004WR003403
Lopez, W., & Meigs, A. (2016). Post 6.3 Ma erosionally driven rock uplift and
landscape evolution ofthe Cascade Range in the Pacific northwest
(Unpublished doctoral dissertation). Oregon State University.
Luo, W., Grudzinski, B. P., & Pederson, D. (2010). Estimating hydraulic
conductivity from drainage patterns-a case study in the Oregon Cascades.
Geology , 38 (4), 335–338. doi: 10.1130/G30816.1
Luo, W., & Stepinski, T. (2008). Identification of geologic contrasts from landscape
dissection pattern: An application to the Cascade Range, Oregon, USA.
Geomorphology , 99 (1-4), 90–98. doi: 10.1016/j.geomorph.2007.10.014
Lyman, W. D. (1963). The Columbia River : its history, its myths, its scenery, its
commerce ([4th ed.]. ed.). Portland, Or: Binfords & Mort.
MacVicar, B. J., & Roy, A. G. (2007). Hydrodynamics of a forced riffle pool in a
gravel bed river: 1. Mean velocity and turbulence intensity. Water Resources
Research, 43 (12). doi: 10.1029/2006WR005272
Maguire, R., Schmandt, B., Li, J., Jiang, C., Li, G., Wilgus, J., & Chen, M. (2022).
Magma accumulation at depths of prior rhyolite storage beneath yellowstone
caldera. Science, 378 (6623), 1001–1004.
Manga, M. (1999). On the timescales characterizing groundwater discharge at
springs. Journal of Hydrology , 219 (1-2), 56–69. doi:
10.1016/S0022-1694(99)00044-X
McCaffrey, R., King, R. W., Payne, S. J., & Lancaster, M. (2013). Active tectonics
of northwestern U.S. inferred from GPS-derived surface velocities. Journal of
Geophysical Research: Solid Earth, 118 (2), 709–723. doi:
10.1029/2012JB009473
162
Mcclaughry, J. D., Wiley, T. J., Conrey, R. M., Jones, C. B., & Lite, K. E. (2012).
Digital geologic map of the Hood River Valley, Hood River and Wasco
counties, Oregon (Tech. Rep.). Portland, Or.
Mcnutt, M. (2005). Influence Of Plate Subduction On Isostatic Compensation in
Northern California. , 11 (4), 253.
Minár, J., Evans, I. S., & Jenčo, M. (2020). A comprehensive system of definitions
of land surface (topographic) curvatures, with implications for their
application in geoscience modelling and prediction. Earth-Science Reviews ,
211 (September 2019). doi: 10.1016/j.earscirev.2020.103414
Mitchell, C. E., Vincent, P., Weldon, R. J., & Richards, M. A. (1994). Present-day
vertical deformation of the Cascadia Margin, Pacific northwest, United
States. Journal of Geophysical Research, 99 (B6). doi: 10.1029/94jb00279
Molloy, I., & Stepinski, T. F. (2007). Automatic mapping of valley networks on
Mars. Computers & Geosciences , 33 (6), 728–738. Retrieved from https://
www.sciencedirect.com/science/article/pii/S0098300406002159 doi:
https://doi.org/10.1016/j.cageo.2006.09.009
Montgomery, D. R., & Foufoula-Georgiou, E. (1993). Channel network source
representation using digital elevation models. Water Resources Research,
29 (12), 3925–3934. doi: 10.1029/93WR02463
Mosegaard, K., & Tarantola, A. (1995). Monte Carlo sampling of solutions to
inverse problems. Journal of Geophysical Research, 100 (B7). doi:
10.1029/94jb03097
Muller, J. R., Ito, G., & Martel, S. J. (2001). Effects of volcano loading on dike
propagation in an elastic half-space. Journal of Geophysical Research,
106 (B6), 11101–11113.
Muller, V. A., Sternai, P., Sue, C., Simon-Labric, T., & Valla, P. G. (2022).
Climatic control on the location of continental volcanic arcs. Scientific
Reports , 12 (1), 1–13. Retrieved from
https://doi.org/10.1038/s41598-022-26158-2 doi:
10.1038/s41598-022-26158-2
Mynatt, I., Bergbauer, S., & Pollard, D. D. (2007). Using differential geometry to
describe 3-D folds. , 29 , 1256–1266. doi: 10.1016/j.jsg.2007.02.006
Needham, T. (2021). Visual Differential Geometry and Forms: A Mathematical
Drama in Five Acts. Princeton: Princeton University Press.
163
Nie, J., Ruetenik, G., Gallagher, K., Hoke, G., Garzione, C. N., Wang, W., . . . Liu,
S. (2018). Rapid incision of the Mekong River in the middle Miocene linked
to monsoonal precipitation. Nature Geoscience, 11 (12), 944–948. Retrieved
from http://dx.doi.org/10.1038/s41561-018-0244-z doi:
10.1038/s41561-018-0244-z
O’Connor, J. E., Wells, R. E., Bennett, S. E., Cannon, C. M., Staisch, L. M.,
Anderson, J. L., . . . Evarts, R. C. (2021). Arc versus river-The geology of
the Columbia River Gorge. GSA Field Guides , 62 (05), 131 -186. doi:
10.1130/2021.0062(05)
O’Hara, D., Karlstrom, L., & Roering, J. J. (2019). Distributed landscape response
to localized uplift and the fragility of steady states. Earth and Planetary
Science Letters , 506 , 243–254. Retrieved from
https://doi.org/10.1016/j.epsl.2018.11.006 doi:
10.1016/j.epsl.2018.11.006
O’Hara, D., Klema, N., & Karlstrom, L. (2021). Development of magmatic
topography through repeated stochastic intrusions. Journal of Volcanology
and Geothermal Research, 419 , 107371. Retrieved from
https://doi.org/10.1016/j.jvolgeores.2021.107371 doi:
10.1016/j.jvolgeores.2021.107371
O’Neill, B. (2006). Elementary differential geometry (Rev. 2nd e ed.). Burlington,
MA ;: Academic Press, an imprint of Elsevier.
Osserman, R. (1969). A survey of minimal surfaces. New York: Van Nostrand
Reinhold Co.
O’Hara, D., Karlstrom, L., & Ramsey, D. W. (2020). Time-evolving surface and
subsurface signatures of Quaternary volcanism in the Cascades arc. Geology ,
48 (Xx), 1–6. doi: 10.1130/g47706.1
Passalacqua, P., Do Trung, T., Foufoula-Georgiou, E., Sapiro, G., & Dietrich,
W. E. (2010). A geometric framework for channel network extraction from
lidar: Nonlinear diffusion and geodesic paths. Journal of Geophysical
Research, 115 (F1), 1–18. doi: 10.1029/2009jf001254
Pasternack, G. B., Ellis, C. R., Leier, K. A., Vallé, B. L., & Marr, J. D. (2006).
Convergent hydraulics at horseshoe steps in bedrock rivers. Geomorphology ,
82 (1-2), 126–145. doi: 10.1016/j.geomorph.2005.08.022
Paulatto, M., Hooft, E. E., Chrapkiewicz, K., Heath, B., Toomey, D. R., & Morgan,
J. V. (2022). Advances in seismic imaging of magma and crystal mush.
Frontiers in Earth Science, 10 .
164
Pearce, M. A., Jones, R. R., Smith, S. A., McCaffrey, K. J., & Clegg, P. (2006).
Numerical analysis of fold curvature using data acquired by high-precision
GPS. Journal of Structural Geology , 28 (9), 1640–1646. doi:
10.1016/j.jsg.2006.05.010
Penrose, R. (2004). The road to reality : a complete guide to the laws of the
universe (1st Americ ed.). New York: A.A. Knopf.
Penserini, B. D., Roering, J. J., & Streig, A. (2017). A morphologic proxy for
debris flow erosion with application to the earthquake deformation cycle,
Cascadia Subduction Zone, USA. Geomorphology , 282 , 150–161. Retrieved
from http://dx.doi.org/10.1016/j.geomorph.2017.01.018 doi:
10.1016/j.geomorph.2017.01.018
Perkins, J. P., Ward, K. M., De Silva, S. L., Zandt, G., Beck, S. L., & Finnegan,
N. J. (2016). Surface uplift in the Central Andes driven by growth of the
Altiplano Puna Magma Body. Nature Communications , 7 , 13185–13195.
Retrieved from http://dx.doi.org/10.1038/ncomms13185 doi:
10.1038/ncomms13185
Perron, J. T., Kirchner, J. W., & Dietrich, W. E. (2008). Spectral signatures of
characteristic spatial scales and nonfractal structure in landscapes. Journal
of Geophysical Research: Earth Surface, 113 (4), 1–14. doi:
10.1029/2007JF000866
Perron, J. T., Kirchner, J. W., & Dietrich, W. E. (2009). Formation of evenly
spaced ridges and valleys. Nature, 460 (7254), 502–505. doi:
10.1038/nature08174
Perron, J. T., & Royden, L. (2013). An integral approach to bedrock river profile
analysis. Earth Surface Processes and Landforms , 38 (6), 570–576. doi:
10.1002/esp.3302
Pesek, M. E., Perez, N. D., Meigs, A., Rowden, C. C., & Giles, S. M. (2020).
Exhumation Timing in the Oregon Cascade Range Decoupled From
Deformation, Magmatic, and Climate Patterns. Tectonics , 39 (9), 1–21. doi:
10.1029/2020TC006078
Pierson, T. C., Evarts, R. C., & Bard, J. A. (2016). Landslides in the western
Columbia Gorge, Skamania County, Washington. USGS Scientific
Investigations Map, 3358 , 25. Retrieved from
http://pubs.er.usgs.gov/publication/sim3358 doi: 10.3133/sim3358
165
Reidel, S. P., Camp, V. E., Tolan, T. L., Kauffman, J. D., & Garwood, D. L.
(2013). Tectonic evolution of the Columbia River flood basalt province.
Special Paper of the Geological Society of America, 497 (June 2014),
293–324. doi: 10.1130/2013.2497(12)
Reidel, S. P., & Tolan, T. L. (2013a). The Grande Ronde Basalt , Columbia River
Basalt Group. Geological Society of America Special Papers , 497 (05),
117–153. doi: 10.1130/2013.2497(05).For
Reidel, S. P., & Tolan, T. L. (2013b). The late Cenozoic evolution of the Columbia
River system in the Columbia River flood basalt province. Special Paper of
the Geological Society of America, 497 (November 1987), 201–230. doi:
10.1130/2013.2497(08)
Reneau, S. L., & Dietrich, W. E. (1991). Erosion rates in the southern oregon coast
range: Evidence for an equilibrium between hillslope erosion and sediment
yield. Earth Surface Processes and Landforms , 16 (4), 307–322. doi:
10.1002/esp.3290160405
Reuter, H. I., Hengl, T., Gessler, P., & Soille, P. (2009). Preparation of DEMs for
geomorphometric analysis. Developments in Soil Science, 33 (C), 87–120.
doi: 10.1016/S0166-2481(08)00004-4
Richardson, P., & Karlstrom, L. (2019). The multi-scale influence of topography on
lava flow morphology. Volcanology Bulletin, 1–42.
Roering, J. J., Kirchner, J. W., & Dietrich, W. E. (1999). Evidence for nonlinear,
diffusive sediment transport on hillslopes and implications for landscape
morphology. Water Resources Research, 35 (3), 853–870. doi:
10.1029/1998WR900090
Roering, J. J., Kirchner, J. W., & Dietrich, W. E. (2005). Characterizing structural
and lithologic controls on deep-seated landsliding: Implications for
topographic relief and landscape evolution in the Oregon Coast Range, USA.
Bulletin of the Geological Society of America, 117 (5-6), 654–668. doi:
10.1130/B25567.1
Roering, J. J., Perron, J. T., & Kirchner, J. W. (2007). Functional relationships
between denudation andhillslope form and relief. Earth and Planetary
Science Letters , 264 (1-2), 245–258. doi: 10.1016/j.epsl.2007.09.035
Royden, L., & Taylor Perron, J. (2013). Solutions of the stream power equation
and application to the evolution of river longitudinal profiles. Journal of
Geophysical Research: Earth Surface, 118 (2), 497–518. doi:
10.1002/jgrf.20031
166
Ruh, J. B. (2020). Numerical modeling of tectonic underplating in accretionary
wedge systems. Geosphere, 16 (6), 1385–1407. doi: 10.1130/GES02273.1
Saar, M. O., & Manga, M. (2004). Depth dependence of permeability in the
Oregon Cascades inferred from hydrogeologic, thermal, seismic, and
magmatic modeling constraints. Journal of Geophysical Research: Solid
Earth, 109 (4), 1–19. doi: 10.1029/2003JB002855
Sandwell, D. T. (1984). Thermomechanical evolution of oceanic fracture zones.
Journal of Geophysical Research, 89 (B13), 11401–11413. doi:
10.1029/JB089iB13p11401
Scheingross, J. S., & Lamb, M. P. (2017). A Mechanistic Model of Waterfall Plunge
Pool Erosion into Bedrock. Journal of Geophysical Research: Earth Surface,
122 (11), 2079–2104. doi: 10.1002/2017JF004195
Scherler, D., & Schwanghart, W. (2019). Identification and ordering of drainage
divides in digital elevation models. Earth Surf. Dynam. Discuss.,
2019 (October), 1–35. Retrieved from https://
www.earth-surf-dynam-discuss.net/esurf-2019-51/%0Ahttps://
www.earth-surf-dynam-discuss.net/esurf-2019-51/esurf-2019-51.pdf
doi: 10.5194/esurf-2019-51
Schmidt, J., Evans, I. S., & Brinkmann, J. (2003). Comparison of polynomial
models for land surface curvature calculation. International Journal of
Geographical Information Science, 17 (8), 797–814. doi:
10.1080/13658810310001596058
Schoenbohm, L. M., Whipple, K. X., Burchfiel, B. C., & Chen, L. (2004).
Geomorphic constraints on surface uplift, exhumation, and plateau growth
in the Red River region, Yunnan Province, China. Bulletin of the Geological
Society of America, 116 (7-8), 895–909. doi: 10.1130/B25364.1
Schwanghart, W., & Scherler, D. (2014). TopoToolbox 2 – MATLAB-based
software for topographic analysis and modeling in Earth surface sciences.
Earth Surface Dynamics , 2 , 1–7. doi: 10.5194/esurf-2-1-2014
Seidl, M. A., Dietrich, W. E., & Kirchner, J. W. (1994). Longitudinal Profile
Development into Bedrock: An Analysis of Hawaiian Channels. The Journal
of Geology , 102 (4), 457–474.
Shary, P. A. (1995). Land surface in gravity points classification by a complete
system of curvatures. Mathematical Geology , 27 (3), 373–390. doi:
10.1007/BF02084608
167
Sherrod, D. R., & Smith, J. G. (2000). Geologic Map of Upper Eocene to Holocene
Volcanic and Related Rocks of the Cascade Range. Ground Water , 80225.
doi: 10.3133/ds842
Simoes, M., Braun, J., & Bonnet, S. (2010). Continental-scale erosion and
transport laws: A new approach to quantitatively investigate macroscale
landscapes and associated sediment fluxes over the geological past.
Geochemistry, Geophysics, Geosystems , 11 (9). doi: 10.1029/2010GC003121
Singer, B. S., Le Mével, H., Licciardi, J. M., Córdova, L., Tikoff, B., Garibaldi, N.,
. . . Feigl, K. L. (2018). Geomorphic expression of rapid Holocene silicic
magma reservoir growth beneath Laguna del Maule, Chile. Science
Advances , 4 (6), 1–11. doi: 10.1126/sciadv.aat1513
Smith, S., Holland, D., & Longley, P. (2004). The importance of understanding
error in Lidar digital elevation. International Archives of the
Photogrammetry, Remote Sensing and Spatial Information Sciences ,
35 (January 2004), 996–1001. Retrieved from https://www.isprs.org/
proceedings/xxxv/congress/comm4/papers/488.pdf
Smith, T., Rheinwalt, A., & Bookhagen, B. (2019). Determining the Optimal Grid
Resolution for Topographic Analysis on an Airborne Lidar Dataset. Earth
Surface Dynamics Discussions , 1–25. doi: 10.5194/esurf-2018-96
Snyder, N. P., Whipple, K. X., Tucker, G. E., & Merritts, D. J. (2000). Stream
profiles in the Mendocino triple junction region, northern California. GSA
Bulletin, 112 (8), 1250–1263. doi: 10.1130/0016-7606(2000)112⟨1250
Sparks, R. S. J., & Cashman, K. V. (2017). Dynamic magma systems: Implications
for forecasting volcanic activity. Elements , 13 (1), 35–40. doi:
10.2113/gselements.13.1.35
Stewart, S. A., & Podolski, R. (1998). Curvature analysis of gridded geological
surfaces. Geological Society Special Publication, 127 , 133–147. doi:
10.1144/GSL.SP.1998.127.01.11
Stock, J., & Dietrich, W. E. (2003). Valley incision by debris flows: Evidence of a
topographic signature. Water Resources Research, 39 (4). doi:
10.1029/2001WR001057
Stoesser, T., Kara, S., MacVicar, B., & Best, J. L. (2010). Turbulent Flow over a
mildly sloped pool-riffle sequence. Proceedings of the IAHR River Flow 2010
Conference, Bundesanstalt fuer Wasserbau, 409–417.
168
Struble, W. T., Mcguire, L. A., Mccoy, S. W., Barnhart, K. R., & Marc, O. (2023).
Debris-Flow Process Controls on Steepland Morphology in the San Gabriel
Mountains , California Journal of Geophysical Research : Earth Surface. ,
1–29. doi: 10.1029/2022JF007017
Struble, W. T., & Roering, J. J. (2021). Hilltop curvature as a proxy for erosion
rate: Wavelets enable rapid computation and reveal systematic
underestimation. Earth Surface Dynamics , 9 (5), 1279–1300. doi:
10.5194/esurf-9-1279-2021
Struik, D. J. D. J. (1950). Lectures on classical differential geometry. Cambridge,
Mass: Addison-Wesley Press.
Swanson, R. D. (1986). A stratigraphic-geochemical study of the Troutdale
Formation and Sandy River Mudstone in the Portland Basin and lower
Columbia River Gorge (Doctoral dissertation). Retrieved from
http://pdxscholar.library.pdx.edu/open access etds/3720/
Sweeney, K. E., & Roering, J. J. (2017). Rapid fluvial incision of a late Holocene
lava flow: Insights from LiDAR, alluvial stratigraphy, and numerical
modeling. Bulletin of the Geological Society of America, 129 (3-4), 500–512.
doi: 10.1130/B31537.1
Tague, C., & Grant, G. E. (2004). A geological framework for interpreting the
low-flow regimes of Cascade streams, Willamette River Basin, Oregon.
Water Resources Research, 40 (4), 1–9. doi: 10.1029/2003WR002629
Takei, Y. (2002). Effect of pore geometry on vp/vs: From equilibrium geometry to
crack. Journal of Geophysical Research: Solid Earth, 107 (B2), ECV–6.
Tarantola, A. (2004). Inverse problem theory. Philadelphia, PA: Society for
Industrial and Applied Mathematics. Retrieved from
http://www.ipgp.fr/~tarantola/Files/Professional/SIAM/
InverseProblemTheory.pdf%0Apapers2://publication/uuid/
62A7A164-F0E3-4298-AA2B-222FF673BF4E
Till, C. B., Kent, A. J., Abers, G. A., Janiszewski, H. A., Gaherty, J. B., & Pitcher,
B. W. (2019). The causes of spatiotemporal variations in erupted fluxes and
compositions along a volcanic arc. Nature Communications , 10 (1).
Retrieved from http://dx.doi.org/10.1038/s41467-019-09113-0 doi:
10.1038/s41467-019-09113-0
Tolan, T. L., & Beeson, M. H. (1984). Intracanyon Flows of the Columbia River
Basalt Group in the Lower Columbia River Gorge and Their Relationship
To the Troutdale Formation. Bulletin of the Geological Society of America,
95 (4), 463–477. doi: 10.1130/0016-7606(1984)95⟨463:IFOTCR⟩2.0.CO;2
169
Turcotte, D. L., & Schubert, G. (2014). Geodynamics (Third edit ed.). Cambridge,
United Kingdom: Cambridge University Press.
Ueki, K., & Iwamori, H. (2016). Density and seismic velocity of hydrous melts
under crustal and upper mantle conditions. Geochemistry, Geophysics,
Geosystems , 17 (5), 1799–1814.
Venditti, J. G. (2007). Turbulent flow and drag over fixed two- and
three-dimensional dunes. Journal of Geophysical Research: Earth Surface,
112 (4), 1–21. doi: 10.1029/2006JF000650
Venditti, J. G., Li, T., Deal, E., Dingle, E., & Church, M. (2019). Struggles with
stream power: Connecting theory across scales. Geomorphology(xxxx),
106817. Retrieved from
https://linkinghub.elsevier.com/retrieve/pii/S0169555X19302892
doi: 10.1016/j.geomorph.2019.07.004
Venditti, J. G., Rennie, C. D., Bomhof, J., Bradley, R. W., Little, M., & Church,
M. (2014). Flow in bedrock canyons. Nature, 513 (7519), 534–537. doi:
10.1038/nature13779
Wanke, M. (2019). Geochemical and petrological diversity of mafic magmas from
Mount St . Helens. Contributions to Mineralogy and Petrology , 174 (1), 1–25.
Retrieved from http://dx.doi.org/10.1007/s00410-018-1544-4 doi:
10.1007/s00410-018-1544-4
Wanke, M., Clynne, M. A., von Quadt, A., Vennemann, T. W., & Bachmann, O.
(2019). Geochemical and petrological diversity of mafic magmas from mount
st. helens. Contributions to Mineralogy and Petrology , 174 , 1–25.
Wanke, M., Karakas, O., & Bachmann, O. (2019). The genesis of arc dacites: the
case of mount st. helens, wa. Contributions to Mineralogy and Petrology ,
174 , 1–14.
Watts, A. B. (2001). Isostasy and flexure of the lithosphere. Cambridge, United
Kingdom: Cambridge University Press.
Wells, R., Haugerud, R., Niem, A., Niem, W., Ma, L., Evarts, R., . . . Sawlan, M.
(2020). Geologic map of the Greater Portland Metropolitan Area and
surrounding region, Oregon and Washington. U.S. Geological Survey,
Scientific Investigations Map SIM 3443 .
Wells, R. E. (1998). Fore-arc migration in Cascadia and its neotectonic significance.
Geology , 26 (8), 759–762. doi:
10.1130/0091-7613(1998)026⟨0759:FAMICA⟩2.3.CO;2
170
Wells, R. E., & McCaffrey, R. (2013). Steady rotation of the Cascade arc. Geology ,
41 (9), 1027–1030. doi: 10.1130/G34514.1
Whipple, K. X. (2004). Bedrock rivers and the geomorphology of active orogens.
Annual Review of Earth and Planetary Sciences , 32 , 151–185. doi:
10.1146/annurev.earth.32.101802.120356
Whipple, K. X., DiBiase, R. A., & Crosby, B. T. (2013). Treatise on
Geomorphology. Academic Press , 9 , 550–573. doi:
10.1016/B978-0-12-374739-6.00254-2
Whipple, K. X., & Tucker, G. E. (1999). Dynamics of the stream-power river
incision model: Implications for height limits of mountain ranges, landscape
response timescales, and research needs. Journal of Geophysical Research:
Solid Earth, 104 (B8), 17661–17674. Retrieved from
http://doi.wiley.com/10.1029/1999JB900120 doi:
10.1029/1999JB900120
White, S. M., Crisp, J. A., & Spera, F. J. (2006). Long-term volumetric eruption
rates and magma budgets. Geochemistry, Geophysics, Geosystems , 7 (3).
doi: 10.1029/2005GC001002
Williams, C. F., & DeAngelo, J. (2008). Mapping geothermal potential in the
Western United States. Transactions - Geothermal Resources Council , 32 ,
155–161.
Williams, I. A. (1923). The Columbia River Gorge: Its Geologic History Interpreted
from the Columbia River Highway. Portland, Or.: James, Kerns & Abbot.
Williams, S. G., & Furbish, D. J. (2021). Particle energy partitioning and
transverse diffusion during rarefied travel on an experimental hillslope.
Earth Surface Dynamics , 9 (4), 701–721. doi: 10.5194/esurf-9-701-2021
Wilson, A. M., & Russell, J. K. (2020). Glacial pumping of a magma-charged
lithosphere: A model for glaciovolcanic causality in magmatic arcs. Earth
and Planetary Science Letters , 548 , 116500. Retrieved from
https://doi.org/10.1016/j.epsl.2020.116500 doi:
10.1016/j.epsl.2020.116500
Wobus, C., Whipple, K. X., Kirby, E., Snyder, N., Johnson, J., Spyropolou, K., . . .
Sheehan, D. (2006). Tectonics from topography: Procedures, promise, and
pitfalls. Special Paper of the Geological Society of America, 398 (04), 55–74.
doi: 10.1130/2006.2398(04)
171
Wrick, J. R., & Pasternack, G. B. (2008). Modeling energy dissipation and
hydraulic jump regime responses to channel nonuniformity at river steps.
Journal of Geophysical Research: Earth Surface, 113 (3), 1–25. doi:
10.1029/2007JF000873
172