PLACE IN RETURN BOX to remove this checkout from your record. TO AVOID FINES return on or before date due. MAY BE RECALLED with earlier due date if requested. DATE DUE DATE DUE DATE DUE JAM 0,2011 'CT 7U I h— l A Q T 01 v‘ V- U! I0 9 1.. ‘ In 37.3 c.1107; 6/07 p:/ClRC/Date0ue.indd-p.1 MAPPING AND MODELING GLACIAL LAKE ALGONQUIN IN NORTHERN MICHIGAN AND WESTERN ONTARIO WITH MODELS OF UNCERTAINTY by Scott A. Drzyzga A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree Of DOCTOR OF PHILOSOPHY Department of Geography 2007 ABSTRACT MAPPING AND MODELING GLACIAL LAKE ALGONQUIN IN NORTHERN MICHIGAN AND WESTERN ONTARIO WITH MODELS OF UNCERTAINTY By Scott A. Drzyzga This dissertation is the first work to map the ”Upper Group" sequence of glacial Lake Algonquin water planes in northern Michigan and western Ontario. Position data were collected from relict shoreline features, which were uplifted during isostatic readjustment after deglaciation. Geostatistical models were built to describe regional uplift trends and reconstruct the now warped water planes associated with the Main, Ardtrea, Upper Orillia and Lower Orillia phases. Each water plane model was used to adjust a digital contemporary elevation model and create a digital pro-glacial elevation model, which reflected both local topography and the regional pattern of isostatic depression that existed during that particular lake phase. Data were subjected to accuracy assessments that permitted models Of uncertainty about input data errors to be included in each landscape reconstruction process. Accordingly, each map depicts the extent of a former lake phase and expresses two components of uncertainty about predicted paleoshoreline locations: one associated with the precision of the predicted water plane elevation; and another associated with co-located uncertainty about error in the digital contemporary elevation model. Important findings of this work include, but are not limited to: a) Main phase and lower shorelines found on Cockburn Island, Ontario, suggest the island was deglaciated before the outlet at North Bay, Ontario was opened; b) relict features at Beaver Island, Michigan, indicate the Main and Ardtrea phases extended into the Michigan basin, but data collected near Little Traverse Bay, Michigan, suggest many ”highest” glacial Lake Algonquin Shorelines in the Michigan basin do not correlate with the Main phase; and c) similarities between isobases mapped during this work and recently published isobases that Show contemporary uplift rates suggest that the regional pattern of isostatic readjustment have not changed much during the last 11,000 years. This study also found evidence of convergence between the Main and Ardtrea water planes, which suggests differential uplift relative to the Fenelon Falls outlet was the cause of shoreline deformation and tilt during these phases. This study does not resolve the question, however, Of where glacial Lake Algonquin stood in the southern parts of the Huron and Michigan basins. Copyright by SCOTT A. DRZYZGA 2006 ACKNOWLEDGEMENTS I owe a special thanks to Ashton Shortridge for accepting the responsibility of advising me during his first and busy years at Michigan State University, and for his time and commitment ever Since. Not only has he encouraged me to think in new ways and supported my professional endeavors, he has encouraged me to express my opinions in my writing and to not conceal them behind the curtain Of absolute Objectivity. I am grateful for each member of my doctoral committee. I thank Randy Schaetzl for recruiting me into his Seminar in Physical Geography; the experience truly inspired this work and provided me with valuable opportunities to publish early in my career. I thank Judy Olson for helping me to think through many mapping problems; not the least was learning to understand how map projections and coordinate systems really work. From Iudy, I also learned the values of knowing your audience and good map design. And, I am grateful for the many conversations I shared with Alan Arbogast; particularly one on a hot day at Ludington State Park. We climbed atop a large sand dune and looked over many others situated between Lake Michigan and Hamlin Lake. Alan talked about how he could see the landscape evolve in his mind’s eye. He could see winds, changes in water level, stabilization and destabilization, and sands burying soils. At the time, I just stood there and wondered because I saw only sand and junipers. That conversation stuck with me, however, because it let me realize I needed to ”see” Lake Algonquin as Alan had ”seen” his dunes if I were to be successful in this work, and for it to be accepted by others. That moment fully arrived while I, alone on an uninhabited island, was standing on a relict offshore bar and looking at a relict wave-cut bluff; I could see those ancient waves break and those waters swash where today no water is present. I see those waves still and have owned my research ever since. I thank Dr. Grahame Larson, Of the Department of Geology, for serving as the Dean of the College of Social Science’s representative at my dissertation defense. His questions helped me to identify and recognize valuable findings that I had overlooked. I enjoyed learning from him. I extend gratitude and thanks to Mr. L. Avra, Of the Huron Timber Company of Thessalon, Ontario; my work at Cockburn Island would not have been possible without his permission or the assistance provided to me by his foreman and crew. I also thank Mr. H. McQuarrie, a retired logger and seasonal resident of Cockburn Island, who gave me his hand-drafted map Of properties, logging roads and foot trails. Funding for fieldwork was provided by the Association Of American Geographers via a Doctoral Dissertation Improvement Grant, and a Graduate Research Fellowship from the Department of Geography at Michigan State University. vi I extend gratitude and thanks tO Drs. P. F. Karrow, W. R. Cowan and C. F. M. Lewis, and to Kevin Kincare and Beth Weisenborn, for Sharing data and insights regarding glacial Lake Algonquin. I thank Rick and janis Herman, my wonderful cousins, for graciously providing me a home while I was in Michigan and conducting field work. I will always remember our morning coffees, and those afternoons when Rick would help me clear my head with a round of golf and a cold beer. Last but not least, I thank Cathryn Dowd for supporting me while we were in Illinois and I was getting my research off the ground, for helping me return to Michigan to conduct fieldwork, and for putting up with me while we were in Pennsylvania and I was trying simultaneously to finish my dissertation and start a new career. I owe a tremendous debt to her. vii TABLE OF CONTENTS LIST OF TABLES ................................................................................ xi LIST OF FIGURES .............................................................................. xiv CHAPTER 1: INTRODUCTION ............................................................ 1 Background about glacial Lake Algonquin ........................................ 2 Statement of problem ................................................................... 4 Statement Of purpose and research Objectives .................................... 5 Definitions of key terms ................................................................ 6 Error .................................................................................. 7 Accuracy ............................................................................. 7 Precision ............................................................................. 7 Uncertainty ......................................................................... 8 Elevation ............................................................................ 8 Location ............................................................................. 9 Position .............................................................................. 9 Justification for research ............................................................... 9 Organization Of the dissertation ...................................................... 11 CHAPTER 2: LAKE LEVELS AND SHORELINES .................................... 14 Introduction ............................................................................... 14 A word about shorelines ............................................................... 16 Glacial Lake Algonquin ................................................................ 17 Relict shoreline indicators and water level proxies .............................. 27 Variation in Great Lakes water levels ............................................... 36 Short-term events .................................................................. 38 Cyclic fluctuations ................................................................ 40 Long-term trends .................................................................. 43 Summary ............................................................................ 44 Summary ................................................................................... 46 CHAPTER 3: METHODS AND FIELD CAMPAIGNS ................................ 48 Introduction ............................................................................... 48 The initial glacial Lake Algonquin dataset ........................................ 48 Opportunities to improve the glacial Lake Algonquin dataset ............... 51 Shoreline identification, bluff selection, and sampling methods ............. 53 Data processing ........................................................................... 58 Summary of prior mapping Campaigns ............................................ 61 Mapping campaigns conducted for this work .................................... 64 viii Cockburn Island, Ontario ....................................................... 67 Sugar Island, Michigan .......................................................... 74 Rexford, Michigan ................................................................ 76 The Munuscong Islands ......................................................... 79 Kincheloe, Michigan .............................................................. 82 Mackinac Island, Michigan ..................................................... 83 The Douglas Lake area, Michigan ............................................. 87 Beaver Island, Michigan ......................................................... 92 From Petoskey to Traverse City, Michigan .................................. 99 The Leelanau Peninsula, Michigan ........................................... 104 From Mullet Lake to Alpena, Michigan ..................................... 108 From Alpena to the Au Sable River, Michigan ............................ 110 From the Au Sable River to Saganing, Michigan .......................... 112 Summary ................................................................................... 116 CHAPTER 4: ACCURACY ASSESSMENTS ............................................. 126 Introduction ............................................................................... 126 Two datasets .............................................................................. 127 Assessment of the glacial Lake Algonquin dataset .............................. 129 Analysis Of potential selection errors ......................................... 130 Ground control ..................................................................... 131 Analysis of measurement errors ............................................... 134 Cumulative elevation errors in paleoshoreline data ..................... 140 Assessment Of the National Elevation Dataset .................................... 142 Geoprocessing the NED ......................................................... 145 Hypsometric analysis ............................................................ 146 Ground control ..................................................................... 150 Select results of the first exploratory analysis .............................. 152 Complete results Of the second exploratory analysis ..................... 157 Covariation with location, elevation, and slope ........................... 164 Spatial autocorrelation ........................................................... 170 Summary ................................................................................... 170 CHAPTER 5: THE SPATIAL MODEL ..................................................... 175 Introduction ............................................................................... 175 Modeling outcomes of a spatial random process ................................ 176 Building the spatial model ............................................................ 187 Models Of elevation errors in the paleoshoreline elevation data ............. 192 Assumed model of selection errors ........................................... 192 Empirical model of measurement errors .................................... 194 Model of cumulative errors ..................................................... 195 Models of differentially upwarped water planes .................................. 197 A note about coordinate transformation ....................................... 201 Main phase ............................................................................ 202 ix Ardtrea phase ........................................................................ 207 Upper Orillia phase ................................................................. 212 Lower Orillia phase ................................................................. 217 Cedar Point phase ................................................................... 221 Models of errors in the digital elevation model .................................. 223 Variography ........................................................................... 225 Making predictions .................................................................. 230 Running the Spatial model ........................................................ 232 Summary ................................................................................... 236 CHAPTER 6: MAPPING THE RESULTS ................................................. 241 Introduction ............................................................................... 241 Summary of results obtained from the spatial model ........................... 241 Maps of glacial Lake Algonquin shorelines ........................................ 246 The Main phase ...................................................................... 253 The Ardtrea phase .................................................................. 256 The Upper Orillia phase ........................................................... 258 The Lower Orillia phase ........................................................... 260 Maps Of glacial Lake Algonquin Shorelines that reveal uncertainty ........ 261 Composites map Of Main Lake Algonquin and the Upper Group... . .. 264 Summary ................................................................................... 279 CHAPTER 7: CONCLUSIONS .............................................................. 283 Introduction ............................................................................... 283 Re-evaluation of research Objectives ................................................ 284 To find and survey relict shoreline features that best serve as 284 To reconstruct water planes associated with glacial Lake Algonquin 290 To reconstruct the landscape associated with each lake phases and 295 To express the degree of uncertainty with which each glacial Lake 296 Future research ........................................................................... 297 Glacial Lake Algonquin ............................................................ 298 Topographic map analysis ........................................................ 300 Visualizing uncertainty ............................................................ 300 Global climate change and sea-level rise ....................................... 301 Final thoughts ............................................................................. 302 APPENDICES .................................................................................... 297 Appendix A: The glacial Lake Algonquin dataset ............................... 298 Appendix B: Scripts written for use with R ....................................... 317 REFERENCES .................................................................................... 326 LIST OF TABLES Table 2.1: Glacial Lake Algonquin lake phases with time frames. Lake phases are listed in descending order in terms of both contemporary elevation and age. Radiocarbon dates are from Karrow et a1. (1975) ............. 25 Table 2.2: Seasonal variation in contemporary Great Lakes water levels (Quinn, 2002; NOAA, 2006) .................................................................. 37 Table 2.3: Potential range of Great Lakes water-level fluctuations (see text for references) .................................................................................... 45 Table 3.1: Approximate elevations for glacial Lake Algonquin water planes and the Nipissing water plane at Sault Ste. Marie and Manitoulin Island, Ontario (Cowan, 1985) ......................................................................... 62 Table 3.2: Lake level elevations for glacial Lake Algonquin water planes and the Nipissing water plane at the southern end of protO-St. Joseph Island, Ontario (Karrow, 1987) ........................................................................ 63 Table 3.3: Approximate elevations Of glacial Lake Algonquin water planes at Cockburn Island, Ontario .................................................................... 72 Table 3.4: Approximate elevations of glacial Lake Algonquin water planes and the Nipissing water plane measured in western Chippewa County ........ 78 Table 3.5: Approximate elevations of glacial Lake Algonquin water planes measured in the Munuscong Islands region ............................................. 82 Table 3.6: Approximate elevations Of glacial Lake Algonquin water planes measured at Mackinac Island, Michigan .................................................. 87 Table 3.7: Approximate elevations of glacial Lake Algonquin water planes and the Nipissing water plane measured in the Douglas Lake area ............... 91 Table 3.8: Approximate elevations of glacial Lake Algonquin water planes measured at Beaver Island, Michigan ...................................................... 99 Table 4.1: Vertical differences between measured and published heights. N GS bench marks are indicated by their Permanent Identifiers (PIDS) ........... 133 xi Table 4.2: Summary statistics that describe the set of differences between measured and published bench mark elevations ....................................... 135 Table 4.3: Results Of Kendall’s non-parametric tests for associations between position and elevation error .................................................................. 139 Table 4.4: Horizontal differences between measured and published NGS benchmark locations ........................................................................... 140 Table 4.5: Summary statistics that describe differences between actual and co-located elevations in the original and corrected DEM surfaces ................. 154 Table 4.6: Summary statistics that describe differences between actual and co-located DEM elevations ................................................................... 159 Table 4.7: Results Of two-sample difference-Of-means tests among DEM subsets ............................................................................................. 163 Table 4.8: Results of Kendall’s tests for association between easy-tO-measure landscape variables and: a) signed error and b) unsigned error .................... 166 Table 5.1: Main phase trend parameters estimated via ordinary least squares. 203 Table 5.2: Summary Of Main surface residuals Obtained via ordinary least squares (OLS), cross-validated universal kriging with a nugget model only (CV.NUG), and cross-validated universal kriging with a Matern model Of spatial structure (CV.MAT) .................................................................. 203 Table 5.3: Ardtrea phase trend parameters estimated via ordinary least squares (OLS), and generalized least squares (GLS) ................................... 209 Table 5.4: Summary statistics for residuals (n=41) Obtained from Ardtrea water plane models predicted via ordinary least squares (OLS), generalized least squares (GLS), and cross-validated universal kriging (CV.NUG) ........... 210 Table 5.5: Upper Orillia phase trend parameters estimated via ordinary least squares ............................................................................................. 214 Table 5.6: Summary statistics for residuals (n=33) Obtained from Upper Orillia elevation trend models obtained via ordinary least squares (OLS), cross-validated universal kriging with pure nugget model (CV.NUG) and cross-validated universal kriging with anisotropic covariance model (CV .MAT) ........................................................................................ 214 xii Table 5.7: Lower Orillia phase trend model parameters estimated via ordinary least squares (OLS), and generalized least squares (GLS) ............... 218 Table 5.8: Summary statistics for residuals (n=16) Obtained from Lower Orillia water plane models predicted via ordinary least squares (OLS), generalized least squares (GLS) and cross-validated universal kriging (CV .UK) ........................................................................................... 219 Table 5.9: Cedar Point phase trend parameters estimated via ordinary least squares (OLS) .................................................................................... 222 Table 6.1: Estimates of total land area under water during each lake phase. Percentages of contemporary land area under water are indicated in parentheses ....................................................................................... 246 Table A1: The glacial Lake Algonquin datatset ......................................... 298 xiii LIST OF FIGURES Figure 1.1: The extents of glacial Lake Algonquin hypothesized by Hough (1958) and Larsen (1987) ...................................................................... 3 Figure 2.1: Terms used to describe morphologic features along the beach profile (adapted from Komar, 1998) ....................................................... 15 Figure 2.2: Baedke et al. (2004) used buried foreshore deposits as shoreline indicators. Figure adapted from descriptions in their text and from Figure 2.6 in Johnston (2004) .......................................................................... 32 Figure 2.3: Schaetzl et al. (2002) used bouldery/ gravelly lags at the bases of wave-cut bluffs as Shoreline indicators. Figure was adapted from their Figure 2b (p.405) ................................................................................ 34 Figure 3.1: The study area explored by Schaetzl et a1. (2002), points in the original glacial Lake Algonquin dataset, and the extent of the Main phase of glacial Lake Algonquin as they had determined it. Figure adapted from their Figure 1 ............................................................................................ 50 Figure 3.2: Map of the current study area with names Of places and features mentioned in the text. County boundaries and names added for reference ......................................................................................... 52 Figure 3.3: Schaetzl et a1. (2002) used bouldery/ gravelly lags at the bases of wave-cut bluffs as shoreline indicators. Figure was adapted from their Figure 2b (p.405) ................................................................................ 54 Figure 3.4: Index map for smaller-scaled maps that show field campaigns and field sites mentioned in this chapter ................................................ 66 Figure 3.5: Relict bluff sites at Cockburn Island, Ontario. Sites are labeled with a unique site identification number and an elevation value (e. g., 642: 200 m) ............................................................................................. 68 Figure 3.6: Image of a relict bluff on the south side of McCaigs Hill, Cockburn Island, Ontario .................................................................... 70 Figure 3.7: Diagram of the relict bluff shown in Figure 3.6. Shoreline positions were measured at gravel lag situated at the base Of the bluff .......... 70 xiv Figure 3.8: Apparent contours across the Main Lake Algonquin water plane. 73 Figure 3.9: Relict bluff sites at Sugar Island, Michigan. Sites are labeled with a unique site identification number and an elevation value (e. g., 642: 200 m). Also shown are the locations of features surveyed by Cowan (1985) ............. 75 Figure 3.10: Relict bluff sites near Rexford, Michigan. Sites are labeled with a unique site identification number and an elevation value (e. g., 642: 200 m).. 77 Figure 3.11: Relict bluff sites in the Munuscong Islands region. Sites are labeled with a unique site identification number and an elevation value (e. g., 642: 200 m) ....................................................................................... 81 Figure 3.12: Relict bluff sites at Mackinac Island. Sites are labeled with a unique Site identification number and an elevation value (e.g., 642: 200 m). . .. 85 Figure 3.13: Relict bluff sites in the Douglas Lake area. Sites are labeled with a unique site identification number and an elevation value (e. g., 642: 200 m)... 90 Figure 3.14: Relict bluff sites on Beaver Island, Michigan. Sites are labeled with a unique site identification number and an elevation value (e. g., 642: 200 m) ............................................................................................. 94 Figure 3.15: Exposed sediment profile in a gravel pit located near the southern end of protO-Beaver Island. A 2-foot carpenter's level is shown for scale ................................................................................................ 98 Figure 3.16: Relict bluff sites near Petoskey, Michigan. Sites are labeled with a unique site identification number and an elevation value (e.g., 642: 200 m)... 100 Figure 3.17: Relict bluff sites on the Leelanau Peninsula. Sites are labeled with a unique site identification number and an elevation value (e. g., 642: 200 m) ............................................ 103 Figure 3.18: Relict bluff sites in the Black Lake area. Sites are labeled with a unique site identification number and an elevation value (e. g., 642: 200 m). . .. 108 Figure 3.19: Relict bluff sites in the Long Lake area. Sites are labeled with a unique site identification number and an elevation value (e. g., 642: 200 m). . .. 109 Figure 3.20: Relict bluff sites near Alpena, Michigan. Sites are labeled with a unique site identification number and an elevation value (e.g., 642: 200 m). . . .. 111 XV Figure 3.21: Relict bluff sites in near Harrisville, Michigan. Sites are labeled with a unique site identification number and an elevation value (e. g., 642: 200 m) ............................................................................................. 112 Figure 3.22: Relict bluff sites in near Oscoda, Michigan. Sites are labeled with a unique site identification number and an elevation value (e.g., 642: 200 m) ............................................................................................. 114 Figure 4.1: Distribution of error in sampled bench mark elevations: a) Stern and leaf plot - solid line denotes the decimal place; and b) q-q plot of standardized errors with a one-to-one line and 90% confidence envelope ...... 135 Figure 4.2: Scatterplots of Signed error in sampled bench mark elevations and position. Errors and elevations are reported in meters. Latitude and Longitude values are reported in decimal degrees .................................... 137 Figure 4.3: Scatterplots of unsigned error in sampled bench mark elevations and position. Errors and elevations are reported in meters. Latitude and Longitude values are reported in decimal degrees .................................... 138 Figure 4.4: Frequency histograms of elevation. Bin width is one meter .......... 148 Figure 4.5: Histograms of: a) 529 check point elevations; and b) all DEM elevations in the study area. The difference between frequency scales is the results of unequal sample sizes ............................................................. 152 Figure 4.6: Scatterplots Of DEM-reported and actual elevations for the (a) original DEM and for (b) the corrected DEM ............................................ 153 Figure 4.7: Scatterplots of DEM error and the linear components Of slope aspect .............................................................................................. 155 Figure 4.8: Map of DEM errors at control point locations. The map reveals no global pattern but suggests error values exhibit some spatial continuity at short distances. The image is presented in color ........................................ 158 Figure 4.9: Histograms of elevation errors in the DEM .............................. 160 Figure 4.10: Q-Q plots of standardized elevation errors in the DEM ............. 161 Figure 4.11: Scatterplots of: a) signed error and slope gradient; and b) unsigned error and slope gradient. Black lines drawn to emphasize zero error ................................................................................................ 169 xvi Figure 4.12: Map and scatterplot of ground control locations used to assess error in the Ontario subset ................................................................... 169 Figure 5.1: The spatial model Schaetzl et a1. (2002) used to reconstruct landscapes in northern Michigan associated with glacial Lake Algonquin. . . 188 Figure 5.2: The spatial model used in this work to digitally reconstruct water planes and landscapes in northern Michigan and western Ontario that were associated with glacial Lake Algonquin. The model was run 100 times for each water plane, which produced a set of 100 results that were summarized to reveal the influences of input data errors ............................................. 191 Figure 5.3: Histograms of errors generated randomly. Error values are binned in 0.25 m intervals .................................................................... 195 Figure 5.4: Histogram of cumulative paleoshoreline elevation errors. Error values are binned in 0.25 m intervals. Frequency scale is different than that used in Figure 5.2 ............................................................................... 197 Figure 5.5: Empirical variograms of Main trend surface residuals generated via ordinary least squares. Solid lines represent fitted variogram models ...... 205 Figure 5.6: The Main Algonquin water plane: a) predicted elevations are contoured with a 10-meter interval; and b) greater uncertainties are expressed with darker gray tones. Sample sites are shown as black dots ....... 207 Figure 5.7: Empirical variogram model fitted to the Spatial structure of Ardtrea water plane residuals generated during ordinary least squares. . . . . . 210 Figure 5.8: The Ardtrea phase water plane: a) predicted elevations are contoured with a 10-meter interval; and b) greater uncertainties are expressed with darker gray tones. Sample Sites are shown as black dots ....... 211 Figure 5.9: Spatial structure of residuals Obtained during attempts to model the Upper Orillia water plane via ordinary least squares ............................ 215 Figure 5.10: The Upper Orillia phase water plane: a) predicted elevations are contoured with a 10-meter interval; and b) greater prediction uncertainties are expressed with darker gray tones. Sample sites are shown as black dots... 217 Figure 5.11: Spatial structure of 01.5 residuals generated during attempts to model the Lower Orillia water plane ...................................................... 219 xvii Figure 5.12: The Lower Orillia phase water plane: a) predicted elevations are contoured with a 10-meter interval; and b) greater prediction uncertainties are expressed with darker gray tones. Sample sites are shown as black dots... 221 Figure 5.13: Empirical semi-variograms and variogram models for DEM subsets ............................................................................................. 227 Figure 5.14: Composite surface of predicted DEM errors. The image in this dissertation is presented in color ............................................................ 231 Figure 5.15: Uncertainty in the surface of predicted DEM errors. Darker tones are associated with larger confidence intervals and hence greater uncertainties. Areas Shown in white were masked from analysis to facilitate computation ..................................................................................... 232 Figure 5.16: Composite surfaces Of DEM errors. Cells shown in gray were masked from analysis to facilitate computation, but later set to zero. The image in this dissertation is presented in color .......................................... 234 Figure 5.17: Contemporary uplift rates (cm / century) derived from the ICE- 3G geophysical model (Tushingham and Peltier, 1991). Figure adapted from Tushingham (1992: 445) ....................................................................... 238 Figure 5.18: Contemporary uplift rates (cm/ century) derived from historical lake gauge data. Figure adapted from Mainville and Craymer (2005: 1077). 239 Figure 6.1: Predicted shorelines of the Main phase Of glacial Lake Algonquin in the study area. Shorelines are represented as blackened areas; greater widths imply greater uncertainty about shoreline location .......................... 248 Figure 6.2: Predicted shorelines of the Ardtrea phase Of glacial Lake Algonquin in the study area. Shorelines are represented as blackened areas; greater widths imply greater uncertainty about shoreline location ................ 249 Figure 6.3: Predicted shorelines of the Upper Orillia phase of glacial Lake Algonquin in the study area. Shorelines are represented as blackened areas; greater widths imply greater uncertainty about shoreline location ................ 250 Figure 6.4: Predicted shorelines of the Lower Orillia phase of glacial Lake Algonquin in the study area. Shorelines are represented as blackened areas; greater widths imply greater uncertainty about shoreline location ................ 251 xviii Figure 6.5: The Main phase Of glacial Lake Algonquin. Darker tones indicated greater uncertainties about the predicted water plane. The image in this dissertation is presented in color ................................................ 265 Figure 6.6: The Ardtrea phase of glacial Lake Algonquin. Darker tones indicated greater uncertainties about the predicted water plane. The image in this dissertation is presented in color ..................................................... 266 Figure 6.7: The Upper Orillia phase of glacial Lake Algonquin. Darker tones indicated greater uncertainties about the predicted water plane. The image in this dissertation is presented in color ..................................................... 267 Figure 6.8: The Lower Orillia phase Of glacial Lake Algonquin. Darker tones indicated greater uncertainties about the predicted water plane. The image in this dissertation is presented in color ..................................................... 267 Figure 6.9: The Main (dark gray), Ardtrea (pink), Upper Orillia (purple), and Lower Orillia (green) phases Of Lake Algonquin. The image in this dissertation is presented in color ............................................................ 269 Figure 6.10: Areas on Cockburn Island that were predicted to be above glacial Lake Algonquin water levels. The extent of this figure corresponds with the extent of Figure 3.5. The image in this dissertation is in color .......... 270 Figure 6.11: Areas in the Munuscong Island region that were predicted to be above glacial Lake Algonquin water levels. The extent Of this figure corresponds with the extent of Figure 3.11. The image in this dissertation is in color ............................................................................................. 271 Figure 6.12: Areas in the Douglas Lake area that were predicted to be above glacial Lake Algonquin water levels. This extent of this figure corresponds with the extent of Figure 3.13. The image in this dissertation is presented in color ................................................................................................ 272 Figure 6.13: Areas near Charlevoix that were predicted to be above glacial Lake Algonquin water levels. This extent Of this figure corresponds with the extent Of Figure 3.16. The image in this dissertation is presented in color ...... 274 Figure 6.14: Areas south of Alpena, Michigan, predicted to be above glacial Lake Algonquin water levels. This extent of this figure corresponds with the extent Of Figure 3.20. The image in this dissertation is presented in color ...... 275 xix Figure 6.15: Areas near Oscoda, Michigan, predicted to be above glacial Lake Algonquin water levels. This extent Of this figure corresponds with the extent of Figure 3.22. The apparent line west Of Oscoda marks the southern limit of the study area. The image in this dissertation is presented in color. . .. 276 Figure 6.16. The Rapid River outvvash fan/ delta during the: a) Main; b) Ardtrea; c) Upper Orillia; and d) Lower Orillia phases. Pictured is a portion of the USGS 1:24,000 Torch River 7.5 Minute Series Quadrangle, but rescaled to 1:63,500 for display. Darkened areas are predicted to have been under water at with 90% confidence. The image in this dissertation is presented in color .............................................................................................. 277 Figure 6.17. Pullup Lake and the surrounding area during the : a) Main; b) Ardtrea; c) Upper Orillia; and d) Lower Orillia phases. Pictured is a composite of the USGS 1:24,000 Gilchrist and Roberts Corner 7.5 Minute Series Quadrangles, but rescaled to 1:63,500 for display. Darkened areas are predicted to have been under water with 90% confidence. The image in this dissertation is presented in color ............................................................ 278 XX CHAPTER 1: INTRODUCTION This dissertation reports research about reconstructing and visualizing a sequence of pro-glacial lakes that existed in the Great Lakes basin more than 10,000 years ago. The sequence, known as glacial Lake Algonquin, cannot be Observed directly, but its existence is inferred from relict Shoreline features perched above contemporary lake levels. The tasks of mapping relict Shoreline features and drawing lines on maps to reconnect them are made difficult by several factors: some features exhibit evidence Of many water levels while others exhibit evidence of only one; long distances separate many relict features; and the earth’s surface has been differentially lifted and tilted in the Great Lakes region, so different features created along the perimeter Of the same water surface might no longer be situated along the same plane of elevation. Accordingly, there is some confusion in the literature regarding the numbers of water levels and lake phases that existed, and which shoreline features belong to which lake phase (Schaetzl et al., 2002). This work enlists the power Of computers and geostatistical methods to explore and analyze paleoshoreline data, to build and test spatial models of differential uplift and tilt, and to reconstruct some of these ancient landscapes. This research is much more than reconnecting Old bits of Shorelines; it is about working with imperfect data and dealing with uncertainties that propagate through spatial models and into results. Moreover, Shorelines are not static like the lines drawn on maps to represent them. They are dynamic over multiple spatial and temporal scales. This research attempts to describe completely what is known and knowable about shoreline locations. Ultimately, this work develops a new way Of looking at the land-water interface and Offers a stochastic representation --- a fundamental breakthrough in the way shorelines can be visualized. Brief background about glacial Lake Algonquin Glacial Lake Algonquin was a sequence of pro-glacial lakes in the Great Lakes basin (Figure 1.1). It maintained relative high water levels from 11,200 to 10,400 yr BP 1 (see Figure 12 in Karrow, 1975) and developed conspicuous shoreline features (Larson and Schaetzl, 2001). As the Laurentide glacier receded from the basin, lower outlets were uncovered and lake levels fell in stages (Larson and Schaetzl, 2001). The land surface rebounded upward as the weight of the ice was removed; the adjustment process is called isostatic rebound. Locations in northern Michigan and Ontario, which are closer to the former centers of ice loading, have rebounded at rates faster than locations in southern Michigan (The Coordinating Committee on Great Lakes Basic Hydraulic and Hydrologic Data, 1977, 2001). Consequently, landforms associated with the sequence have become differentially lifted and tilted. 1 Conventional radiocarbon ages are reported in years before present (yr BP). The datum, 0 yr BP, corresponds tO calendar year 1950 AD. «4‘ _ Approximate ice- 6%),» Extent of current lakes .n—I margin position State or province Extent of Lake Algonquin Au Train- White Fish Sault Ste. Marie - QUEBEC North 0' N L _‘ _____ * " mats. W ...'-V\ \d I. a k 9 ,4 {T "#59,, M X??? ONTARIO I) é Kirkfield Fenelon * Falls Glacial “c. j ‘ 93939324,“ Grand River _5381”3W ‘ -------- lowla l nd . / ,r/ I ,. ONTARIO ( U!“ T "7 Q} 5 e Ontario :4 “ . " MICHIGAN PortHumn M u.s ,..-.~»~~J\ A ,//fi\fi_ f . - ,3; .e /NEw YORK j V'Nr/ E r ‘ /’/ . .. . ‘ (/ K e ,,,,, ~ , “\x o /v,,,./’ ., . INDIANA OHIO >3:\ -.-_./ PENNSYLVANIA Figure 1.1: The extents of glacial Lake Algonquin hypothesized by Hough (1958) and Larsen (1987). Much has been written about how glacial Lake Algonquin was formed (Spencer, 1891; Gilbert, 1898; Goldthwait, 1908; Leverett and Taylor, 1915; and Stanley, 1936), but little agreement exists about the elevations and spatial extents of its water surfaces (Hough, 1958; Futyma, 1981; Larsen, 1987; and see Figure 1.1). Little agreement exists because the earth’s surface has been lifted and tilted over the past 10,000 years and many relict landforms have been altered by natural and human processes; some were buried or reworked by later lake transgressions. Uncertainty about the timing and extents of this lake sequence prompted an extensive survey by Schaetzl et a1. (2002). The authors employed global positioning technology to measure positions along relict Lake Algonquin features. They used a GIS and a digital elevation model to analyze their paleoshoreline dataset and interpolate unknown shoreline locations. In brief, they used a GIS to model the cumulative effects of the isostatic rebound process over space and estimate the spatial extents of several water surfaces. Statement of problem Spatial data quality is a theme that has ebbed and flowed through the course of literature devoted to Geographic Information Systems (GIS) and Geographic Information Science (GISc). Yet, it seems many people are indifferent to this line of work because the same warnings have been issued time and time again. For example, Goodchild and Gopal (1989) claimed virtually all spatial datasets stored in a GIS are, to some extent, contaminated by error. Thirteen years later, Heuvelink and Burrough repeated the same message: ”data stored in a spatial database are rarely, if ever, truly free of error” (Heuvelink and Burrough, 2002: 111). Burrough (1992) once baited the spatial modeling community to include assessments of spatial data quality and error propagation as ingredients of intelligent spatial model design. The community at large did not take the bait and Heuvelink (1998) was compelled to reiterate the dangers in work that utilizes GIS yet fails to assess errors inherent in the data sets used or propagated through the spatial models employed. Together, Heuvelink and Burrough (2002) have argued it is crucial ” to know how accurate the data contained in spatial databases really are" (p.111), because without that knowledge neither the true value of the information derived, nor the correctness of the decisions they support can be evaluated. Assessing spatial data quality and decontaminating spatial datasets are difficult and often tedious and frustrating tasks for in most cases only imperfect knowledge exists about the amount of error contained in any one observation, let alone among any set of observations. Data errors, inaccuracies, or imprecision, even when acknowledged, are often absorbed into analyses; their impacts on results left unexamined and potential consequences ignored. A need exists to demonstrate the value of assessing spatial data quality and uncertainty propagation in work that utilizes GIS. Perhaps Motro summarized this need best when he wrote: Uncertainty permeates our understanding of the real world. The purpose of information systems is to model the real world. Hence information systems must be able to deal with uncertainty (Motro, 1993, cited in Parsons, 1996: 2). Statement of purpose and research objectives The purpose of this research is to demonstrate the value of assessing spatial data quality, incorporating information and uncertainty about input data errors into spatial model design, and visualizing uncertainties propagated through work accomplished within a GIS environment. To achieve this end, this work reconstructs the paleolandscapes of glacial Lake Algonquin, and shows by example how accuracy assessments, uncertainty models, and visualizations of uncertainty can be used to aide and improve such an undertaking. Accordingly, this research has the following objectives: 1. To find and survey relict shoreline features that best serve as indicators of former lake levels; 2. To conflate new paleoshoreline data with data collected during previous mapping campaigns; 3. To reconstruct water planes associated with the various phases of glacial Lake Algonquin; 4. To reconstruct the landscape associated with each lake phase and map the distribution of shorelines; 5. And, to express the degree of uncertainty with which each shoreline location has been mapped. Definitions of key terms Many data quality concepts are commonly, but erroneously, used inter- changeably. This work relies heavily on several of these concepts, so I offer definitions for them below. Error Error is the difference between an attribute value and its measured representation. All measured values contain error because the precisions of human ability and equipment are finite. The term error includes not only mistakes, bias, and blunders, but the statistical concept of variation. Accuracy Accuracy refers to the quality of a measurement and indicates whether or not it contains an acceptable amount of error. A measurement is accurate if it contains an acceptable amount of error (or less); it is inaccurate if it contains more. Although accuracy can be considered a binary condition (i.e., being one of only two possible states), it is used more commonly to refer to the degree of conformity between measured and true values (Wolf and Dewitt, 2000). Such use permits the accuracy of one measurement to be compared to the accuracy of another in relative terms (e. g., measurement x is more accurate than measurement y). Precision Precision refers to the quality of a measurement process and is commonly used to indicate the exactness of values obtained by that process. According to Wolf and Dewitt (2000), precision can be assessed by taking repeated measurements of a single quantity and examining the variance. The process, and the exactness of values obtained from it, is considered imprecise if the variance exceeds a predetermined threshold. The process, and the exactness of values 7 obtained from it, is considered precise if repeated measurements vary less than the threshold. Uncertainty Western science idealizes prefect data and information because they support rational thought and quick decisions, but uncertainty lurks as an intrinsic property of many concepts, ideas, and data (Couclelis, 2003). According to MacEachren (2005: 140), ” when inaccuracy is known objectively, it can be expressed as error; when it is not known, the term uncertainty applies.” Yet, there is no generally agreed-upon definition of uncertainty, so it is most often treated as a multi-faceted characterization that embodies related concepts including error, accuracy, precision, validity, variability, noise, fuzziness, vagueness, completeness, and reliability (Buttenfield, 1993; Pang, 2001; MacEachren, 2005). I use the term to connote the degree of confidence with which something can be known given available data or information. According to Christakos (2000), determinations and reports of uncertainty are the hallmarks of geostatistical analysis and important to processes related to risk assessment and decision-making. This research also relies heavily on the terms elevation, location and position. I define them below. Elevation Elevation is a 1-dimensional measure of vertical distance, or height. All elevations in this work are referenced to the WGS84 ellipsoid and GEOID03 8 geoid models, which are used to specify the familiar concept of mean sea level in and around the United States. Location Location refers to a 2-dimensional measure of horizontal distances between a point and the origin of a coordinate system. The system can be geographic, in which case locations are specified with pairs of longitude and latitude values. All geographic coordinates in this work are referenced to the WG584 and NAD83 geographic coordinate systems, which are used interchangeably. The system can be planar, in which case locations are expressed as pairs of easting and northing values. Unless noted otherwise, all planar coordinates are referenced to the Michigan GEOREF planar coordinate system. Position Position refers to a 3-dimensional specification of a point in space and includes specifications of both elevation and location (see above). Justification for research This research continues the exploration and mapping traditions of Geographic research and offers new and original contributions. For example, this work offers the first known maps of the Upper Orillia and Lower Orillia shorelines (Deane, 1950) of glacial Lake Algonquin. And, this work offers the first reported survey of relict shoreline features at Cockburn Island, Ontario. Other contributions are incremental because they build upon the geomorphologic research and data collection activities initiated more than a century ago by Spencer (1891), Lawson (1891), and others. For example, this research expands upon the work of Schaetzl et al. (2002) and Drzyzga et al. (2002) by offering updated reinterpretations of extant paleoshoreline data, updated water plane reconstructions, and maps of Lake Algonquin shorelines that cover a geographic area broader than they had studied. New data points also fill some of the large data gaps that Cowan (1985) had mentioned; there be fewer dragons on the map now (Blake, 1999). This research also provides the field of Geography with a fundamental contribution - a new way of looking at the land-water interface. A stochastic representation of a shoreline embodies variation in water levels and, hence, the degree of confidence with which the location of a shoreline, as it is traditionally conceptualized, can be known. While the application of this new method is limited in this research to the reconstruction of ancient lakes in the Great Lake basin, it can also be applied to other paleolimnologic reconstructions and adapted for use in contemporary applications that focus on riparian and coastal land use issues. This research is intrinsically geographic because it synthesizes concepts, knowledge, and tools from multiple disciplines to study a portion of the Earth. It employs methods developed in the field of Geographic Information Science to address a long-standing physiographic and paleolimnologic problem. Some geomorphologists might find issue with the lack of stratigraphy or rheology in 10 this work, for example, and geostatisticians will not find any new theories here, just applications of existing methods. But, it is the synthesis and novel application of ideas borrowed from multiple disciplines to address a long- standing problem that forms the basis of my contribution. As an aside, the idea of large pro-glacial lakes and icy waters washing over vast areas of land we now call Ontario or Michigan is both fascinating and unnerving because it forces the mind to visualize a landscape unlike any that can be seen on earth. It prompts consideration of a landscape with large glaciers and without anthropogenic modifications; a cold, wet, windy, and dangerous place. Yet, glacial Lake Algonquin has inherent value beyond prompting the imagination in such wild ways -- it places the contemporary Great Lakes system into a context of landform evolution. Knowing glacial Lake Algonquin existed reminds us that the Great Lakes have evolved through a history of greater great lakes and are evolving still. Organization of the dissertation The work reported here follows two interleaved threads of research: the first thread is associated with Objectives 1 and 2 and focuses on my efforts to find, survey, and correlate relict shoreline features with known glacial Lake Algonquin lake phases; the second thread is associated with satisfying Objectives 3, 4 and 5 and focuses on my efforts to digitally reconstruct glacial Lake Algonquin water planes and landscapes, and to estimate the degree of 11 uncertainty with which interpolated shorelines have been mapped. Because each thread provides a context for examining the other, some reviews of literature, descriptions of methods, and reports of results are arranged in a noncontiguous way. This dissertation is organized into seven chapters. As is traditional, the first chapter, this one, identifies the scope of the research (i.e., background materials, the research problem, research questions, and the context and justification for the work). Chapter 2 reviews the known history of glacial Lake Algonquin and provides a discussion that recognizes uncertainty as an intrinsic property of lake levels and shorelines. I also make a case for incorporating the dynamic nature of shorelines into a study of ancient lakes. Chapter 3 describes methods used and field campaigns conducted to collect data about Lake Algonquin lake phases. I document the methods used survey, process, and classify paleoshoreline data. Because Objectives 1 and 2 needed to be satisfied before work on Objects 3, 4 and 5 could begin, correlations among collected paleoshoreline data and Lake Algonquin lake phases are also reported in Chapter 3. It just made sense to report relict shoreline features and their correlations with lake phases in one place rather than force the reader to mentally dovetail multiple landform descriptions and lake phase correlations that, if reported in the traditional way, would have spanned multiple chapters. Readers familiar with the seminal work of Leverett and Taylor (1915) will recognize my attempt to describe fieldwork, features, and correlations as they 12 had done. The narrative style of writing is less formal than I use in other chapters, but it lends itself to others who might use the text as a field guide during future research. Chapter 4 offers a comprehensive assessment of the principal datasets I employ in my analysis: the glacial Lake Algonquin dataset and a digital elevation model. Relevant literature and methods for assessing the accuracy of these data are also presented. Chapter 5 presents the crux of the work. I report work to build uncertainty models that accommodate errors in my datasets, build geostatistical models of four differentially upwarped water planes, and incorporate all of this information into a spatial model that predicts whether or not locations in the area of interest were above or below each water plane. Relevant literature is also reviewed. Chapter 6 delivers the results; maps that illustrate paleoshorelines associated with the Main and ”Upper Group” phases of glacial Lake Algonquin, as well as areas where the locations of these shorelines remain uncertain. Figure 6.10 is a large-format map that shows the locations of data points in the glacial Lake Algonquin dataset as well as the predicted shoreline locations for the entire Upper Group sequence of water levels. Chapter 7 summarizes the work and identifies goals worth pursuing in future research. Finally, Appendix A lists every record in the glacial Lake Algonquin dataset as well as every record in the supplemental database that holds surveyed features that were retained for future analysis. 13 CHAPTER 2: LAKE LEVELS AND SHORELINES Introduction This chapter provides a context for surveying, mapping, and modeling a sequence of differentially-uplifted water planes associated with glacial Lake Algonquin. I begin the chapter by introducing glacial Lake Algonquin‘and providing an abridged history of processes and events that lead to its formation, lake-level changes through time, and demise. Next, I review work that employs coastal landforms as shoreline indicators and uses measurements taken from these landforms as proxies for former water levels. Special attention is given to the pioneering work of Lawson (1891), who established a precedent for mapping glacial Lake Algonquin features and correlating sequences stranded at different locations. Last, I offer a discussion of variation in contemporary Great Lakes water levels, which is used to infer variation in glacial Lake Algonquin water levels. The discussion serves two purposes: (1) to estimate the precision with which an average water level can be measured by taking a spot measurement from a relict shoreline feature, and (2) to estimate the precision with which two successive lake phases can be discriminated on the basis of such measurements. I adopted Komar’s (1976, 1998) nomenclature for describing shorezone morphology (Figure 2.1). Larsen (1994) employed Komar’s system of names in his work about glacial Lake Algonquin, so it seems reasonable to continue using it. Readers are referred to Komar’s latest text (1998) for complete definitions. 14 V littoral zone 17 -— offshore+ a inshore foreshore *s— backshore —+‘ If berms -\ beach I ' b'Ufi I beach-// scarp longshore bar Figure 2.1: Terms used to describe morphologic features along the beach profile (adapted from Komar, 1998). The landforms illustrated in Figure 2.1 are based on a discrete model of geographic space. A continuous zone of the earth’s surface is partitioned into discrete conceptual entities (e. g., longshore bar, beach face, berm, or bluff) according to sets of relative geomorphometric (e. g., slope aspect, gradient, curvature), spatial (e. g., landward of a, or seaward of b) or sedimentologic characteristics (e. g., parent material, grain size). Discrete conceptual entities can be represented by cartographic objects (i.e., points, lines, and areas) and located with respect to a known spatial reference system. Digital forms of mapped objects can be attributed, stored, edited, manipulated, queried, analyzed and displayed within a GIS environment. 15 A word about shorelines A shoreline is a conceptual entity that generalizes a dynamic intersection among bodies of air, land, and water into a simple geometric construct. The word shoreline connotes an infinitely thin line along the shore, and thus a high- degree of precision, which lends itself to efficient representation in a visually manageable form. This connotation is evident in Komar’s (1998) definition: ”Shoreline: The line of demarcation between the water and the exposed beach. The water’s edge.” (p.47). In one sense, Komar’s definition is useful because it reduces a complex phenomenon to an object that can be discussed and analyzed; it facilitates scientific inquiry and efficient representation on maps or within a GIS environment. The problem with the definition, however, is that practitioners do not always agree on how to identify a shoreline in the field (Graham et al., 2003), and when they do, it is nonetheless difficult to identify and measure consistently (Parker, 2003). This research recognizes that a line object is an imprecise representation of a shoreline because it fails to capture intrinsic variations inherent to the dynamic air, land, and water interface over the long term. Shorelines move landward as Water levels rise, and they move lakeward as water levels fall. The more gentle the beach slope, the greater the horizontal movement. Shorelines do not move horizontally only when beach slope is infinite (i.e., a vertical cliff or wall). This 16 research attempts to mitigate this imprecision problem by employing a representation that captures the range of locations (i.e., the area) over which interface movements occur given variability in water levels. Extant research about contemporary and historic Great Lakes water-level fluctuations can be used to condition an expectation about variability in glacial Lake Algonquin water levels. Glacial Lake Algonquin Much of what is known today about glacial Lake Algonquin, its shorelines, and its water levels was established at the turn of the twentieth century on the basis of pioneering land surveys and geomorphic reconnaissance missions. Spencer (1891a) is credited as the first person to attach the Algonquin moniker to conspicuous bluff relicts and shoreline strands located ”above the eastern rim of Lake Huron” (Schaetzl eta1., 2002:399). The work of Leverett and Taylor (1915), however, is generally considered the first authoritative history of glacial Lake Algonquin. They integrated a set of previous and disparately located works conducted by themselves and others (e. g., Lawson, 1891 ; Spencer, 1891a; Leverett, 1897; Taylor, 1894 and 1897) in order to make generalizations about the glacial and post-glacial history of the Great Lakes region. The authors report included hypotheses and claims about outlet locations, outlet Chronologies, correlations among features at disparate field sites, enduring lake Stages, and general patterns of differential crustal deformation. Extant research 17 (e.g., Farrand and Eschman, 1974; Schaetzl et al., 2002), this work included, continues to recognize the value of Leverett and Taylor’s (1915) synthesis. Two major contributions have synthesized and thus improved the collective understanding of Great Lakes evolution: Hough’s (1958) Geology of the Great Lakes; and Karrow and Calkin’s (1985) Quaternary Evolution of the Great Lakes. The later represents a collection of individual works that update the lake- Ievel and outlet histories in each of the lake basins. More recent works by Pengelly et al. (1997) and Larson and Schaetzl (2001) offer new summaries of Great Lakes evolution. Another review is warranted here. Water basins in the Great Lakes region are largely products of repeated scouring and erosion as continental glaciers advanced, oscillated, and retreated during the last 2.4 million years; ice flows tended to follow bedrock valley systems that existed before glaciation (Larson and Schaetzl, 2001). Many sites in Michigan and Ontario, as well as others in adjacent states, show geomorphic and stratigraphic evidence of multiple glaciations. Large morainic systems (e.g., the Lake Border, Port Huron or Munising morainic systems) mark the limits of particular glacial advances. The last major glacial advance of the Laurentide ice sheet began 26 ka, covered the entire watershed, and reached a maximum extent 18 ka (Larson and Schaetzl, 2001). After oscillating at the maximum extent for a few thousand years, the ice sheet margin retreated northward, ”but that retreat Was interrupted by several readvances at about 15.5, 13.0, 11.8, and 10.0 ka” (Larson and Schaetzl, 2001: 537). As the southern margin of the Laurentide ice 18 sheet retreated, large pro-glacial lakes formed in topographic sinks where waters were impounded by glacial ice to the north and high topography elsewhere. Relict shoreline features (e. g., wave—cut bluffs, beach ridges, spits, and deltas) mark the limits of pro- and post-glacial lakes that formed during glacial retreats. Lake waters were displaced during readvances and shoreline features were overridden or destroyed by glacial ice. The Laurentide ice sheet withdrew from the watershed roughly 9000 years ago (Larson and Schaetzl, 2001). The highest set of relict shoreline features north of Georgian Bay was identified by Spencer (1891a) as remnants of the earliest knowable lake to occupy the Lake Huron, Michigan and Superior basins simultaneously (Eschman and Karrow, 1985; Larsen, 1987). Other highest shoreline features were discovered in the region and reported to occur at prOgressively lower altitudes from north to south (Larsen, 1987). Lawson (1891), Taylor (1894, 1895, 1897), Leverett (1913) and Goldthwait (1908, 1910) surveyed relict shoreline features, made correlations among features stranded at disparate locations, and mapped former lake extents. They also reported lower shorelines that marked lower and younger lakes. According to Leverett and Taylor (1915), the four highest shorelines occurred at progressively lower elevations until they converged at a point south of Traverse City, Michigan. Larsen (1987) described how converging shorelines and spatial patterns of Shoreline deformation and tilt provided grist for Opposing theories of crustal movement during the late 18003. Gilbert (1890), for example, believed tilted 19 shorelines in the Great Lakes region bore similarities to tilted shorelines found elsewhere in Utah and, thus, the earth’s crust in the Great Lakes region had been depressed by the massive weight of a continental ice sheet during glacial advance and rebounded as the weight of ice diminished during retreat. Gilbert (1890) adopted an elastic model of crustal movement; patterns of shoreline deformation and tilt reflected a continuous process of differential isostatic adjustment since glaciation. Goldthwait (1908, 1910) adopted an alternative model -- a rigid model of crustal movement -- that required a stable southern region and a tectonically active northern region; a hinge line model represented the boundary between these two regions. Goldthwait (1908) assumed crustal stability in the southern region allowed the glacial Lake Algonquin shoreline to remain horizontal while episodic tectonic movements somewhere in the northern region forced crustal uplift. The hinge metaphor provided an explanation for how a single shoreline could split into multiple shorelines and diverge with increasing distance, and it suggested rates of episodic movement increased with distance from the hinge line. Leverett and Taylor (1915) adopted the hinge line model for use in their landmark manuscript 1, which, in effect, installed the rigid model of crustal movement as the proper geologic context for conducting 1 Leverett and Taylor (1915) claimed Lake Algonquin attained its greatest extent before a ”great uplift” (p410) had warped the land. 20 subsequent geomorphic research and interpreting relict shorelines during the next fifty years (e.g., Stanley, 1936; Deane, 1950; and Hough, 1958). Not until the 19708 did works like Clarke and Persoage (1970) and the Coordinating Committee on Great Lakes Basic Hydraulic and Hydrologic Data (1977) report the entire Great Lakes watershed was undergoing differential vertical movements associated with postglacial rebound - and not only those locations situated north of a supposed hinge. The spatial isostatic adjustment pattern was smooth, continuous, and tended to increase from southwest to northeast and towards the former centers of ice loading. Locations along the southern shore of Lake Michigan were adjusting at slower rates than locations near the Port Huron outlet - but adjusting nonetheless. The Coordinating Committee published a second report (Coordinating Committee, 2001) confirming the process was ongoing and noting water levels along the southern shores of the Great Lakes were rising (i.e., transgressing) slowly over time. Others (e.g., Tushingham, 1992; Clark et al., 1994; Mainville and Craymer, 2005) have either repeated the experiment or conducted similar work, and all obtained similar results. These reports debunk Goldthwait’s (1908) rigid earth model. Kaszycki (1985) and Larsen (1987) have since offered evidence that conflicts with the hinge line model and suggested glacial Lake Algonquin shorelines plunge below contemporary lake levels in the supposed southern region. The glacial readvance of 11.8 ka covered the northern half of the Michigan basin and the northwestern end of the Huron basin, and it blocked the Straits of 21 Mackinac that, at times, had allowed water to flow from one basin to the other (Hansel and Mickelson, 1988; Larson and Schaetzl, 2001). Broecker and Farrand (1963) reported a radiocarbon date (11,850 yr BP) obtained from a forest near Two Rivers, Wisconsin, that was overridden and buried by this advance. It is generally believed that pro-glacial Lake Chicago formed in the southern Michigan basin and drained via an outlet near Chicago, Illinois (Leverett, 1897; Hansel and Mickelson, 1988), and pro-glacial Lake Algonquin formed in the southern Huron basin and drained via the Fenelon Falls outlet near Kirkfield, Ontario (Finamore, 1985; Hansel and Mickelson, 1988). Differential isostatic rebound changed the relationship between land surface elevations and each field of equal gravitational potential along which each lake surface rested. For example, Hansel and Mickelson (1988) reported lake transgression occurred during the Calumet phase of Lake Chicago as land was uplifted, and Larson and Schaetzl (2001) reported a contemporaneous lake transgression occurred in the southern end of the Huron basin as the Fenelon Falls outlet was uplifted at a faster rate. As glacial ice retreated from the maximum extent of 11.8. ka, lakes in the Michigan and Huron basins expanded into sinks formerly occupied by ice. The lakes were joined when retreating ice uncovered the Indian River lowlands and, slightly later, the Straits of Mackinac (11,200 yr BP according to Larsen, 1987). Upon joining, water levels in glacial Lake Chicago dropped to the glacial Lake Algonquin level (Hansel et al., 1985) and left the outlet at Chicago abandoned 22 (Larson and Schaetzl, 2001). The magnitude of the drop is unknown (Kaszycki, 1985). Continued ice margin retreat allowed glacial Lake Algonquin to enter the Superior basin (Cowan, 1985; Larson and Schaetzl, 2001). The northern extent of glacial Lake Algonquin is not known, however, because many shoreline features in the Superior basin were overridden or destroyed during the next glacial advance (Cowan, 1985; Farrand and Drexler, 1985). Differential isostatic rebound continued to change the relationship between land surface elevations and the field of equal gravitational potential along which the large surface of glacial Lake Algonquin rested. Eschman and Karrow (1985: 89) remarked the ”Main Algonquin beach south of the [Fenelon Falls] outlet represents the maximum of a transgressing water level” whereas north of the outlet it represents just one of several regressing water levels. Larson and Schaetzl (2001) noted debate exists regarding whether or not transgressing waters reached the Port Huron outlet. For example, Pengelly et al. (1997) insisted that the excavation of the lowest portion of the Niagara Gorge (between the Erie and Ontario basins) cannot be explained if the Erie basin did not receive inflow from the Huron basin, but Mothersill and Schurer (2004) claimed Lake Huron waters bypassed Lake Erie, which allowed extensive plant- choked swamps to form in Lake Erie’s western sediment basins. In either case, a combination of continued ice margin retreat and continued uplift ultimately shifted outflow away from the Fenelon Falls outlet and to the Fossmill outlets that spilled water into the Ottawa River valley Eschman and Farrand, 1985). 23 The sequence of descending shorelines in the northern Great Lakes region marks a sequence of progressively younger lake phases. In descending order (Table 2.1), the sequence is comprised of the highest or Main phase of Lake Algonquin (Spencer, 1891a; Leverett and Taylor, 1915); an Upper Group containing the Ardtrea, Upper Orillia and Lower Orillia shorelines / phases that Deane (1950) identified in Ontario’s Lake Simcoe district; a Lower Group containing the Wyebridge, Penetang, Cedar Point and Payette shorelines / phases that Stanley (1936) identified on Ontario’s Penetanguishene Peninsula; and the very low Sheguiandah and Korah shorelines / phases that Hough (1958) identified near Sault Ste. Marie, Ontario. Both Karrow (1987) and Larsen (1987) have noted the controversy that has surrounded these shorelines because it is not yet known whether glacial Lake Algonquin water levels dropped abruptly as new outlets were uncovered, or the land surface was differentially lifted above the field of equal gravitational potential along which the lake surface rested. Karrow (pers. comm.) recently correlated a sequence of descending outlets with the sequence of descending lake phases, but his results remain unpublished. Nevertheless, it seems reasonable to consider that a sequence of discrete short- term outlet change events occurred during the same period as the continuous long-term process of isostatic readjustment; the cumulative impacts of both sequences could produce a very complicated lake-level history. Glacial Lake Algonquin ended when the North Bay outlet was deglaciated (Hansel et al., 1985) and the large lake drained leaving three pools: Lake 24 Chippewa in the Michigan basin (Hough, 1955; Hansel et al., 1985; Larson and Schaetzl, 2001), Lake Stanley in the Huron basin (Hough, 1955; Hansel et al., 1985; Larson and Schaetzl, 2001), and Lake Hough in the Georgian Bay basin (Lewis, 1969; Larson and Schaetzl, 2001). Karrow et al. (1975) reported radiocarbon dates that indicate glacial Lake Algonquin, from the Main phase to the lowest, existed between 11,200 yr BP and 10,400 yr BP. Individual lake phase initiations and durations, however, have not yet been dated (Table 2.1). Table 2.1: Glacial Lake Algonquin lake phases and time frames. Lake phases are listed in descending order in terms of both contemporary elevation and age. Radiocarbon dates are from Karrow et al. (1975). Name of lake phase Time frame The Main phase (Spencer, 1891) 11,200 yr BP - 7 The Upper Group (Deane, 1950) Ardtrea 7 - 7 Upper Orillia 7 - 7 Lower Orillia 7 - 7 The Lower Group (Stanley, 1936) Wyebridge Penetang Cedar Point Payette Sheguiandah (Hough, 1958) . - Korah (Hough, 1958) ? - 10, 0 v w) w: w: °\J I I O w.) -\J 'x) °\) °\J 4:. yrBP In summary, glacial Lake Algonquin was a sequence of pro-glacial lakes in the Great Lakes basin (Figure 1.1). It maintained relative high water levels from 11,200 to 10,400 BP (Karrow et al., 1975; Larson and Schaetzl, 2001). Each 25 water level is assumed to have rested along a field of equal gravitational potential that can be approximated by a horizontal plane (Clarke et al., 1990). Assuming the same kinds of erosional and depositional processes that occur today also occurred during glacial Lake Algonquin times, wave-action at the perimeter of each lake phase formed shoreline features at or near the horizontal plane. Observations of stranded sets of descending shoreline features support the theory that lake levels fell in discrete stages (Lawson, 1891 ; Spencer, 1891b). The sequence of descending shorelines was attributed to a series of lake phases (Leverett and Taylor, 1915; Stanley, 1936; Deane, 1950; Hough, 1958). Observations of tilt among relict shorelines support the theory that water planes, and hence the entire land surface, were differentially uplifted as the glacier retreated and the weight of ice diminished (Gilbert, 1898; Schaetzl et al., 2002). Observations of greater tilt at northern locations and lesser tilt at southern locations support the theory that the former ice sheet originated in the north and, after reaching maxima extents, retreated to the north. Observations of greater tilt rates among higher shorelines and lesser tilt rates among lower shorelines support the theory that rates of differential uplift have decreased over time (Larsen, 1987). Together, these theories and the observations that support them provide an explanation for the origin of the contemporary Great Lakes: they evolved through a history of greater great lakes. They also provide the basis for reconstructing glacial Lake Algonquin water planes. 26 Relict shoreline indicators and water level proxies Lawson’s (1891) study of abandoned shorelines along the north side of Lake Superior seems to have established a precedent for mapping relict shorelines and correlating stranded sequences across space. Although Taylor (18971111) did not agree with many of his interpretations, he did praise ”the high value of Prof. Lawson’s observations.” Lawson’s data collection methods, as well as the statements he made about data uncertainty, have particular relevance to this work about mapping and modeling paleolake levels and shorelines, so several excerpts are examined below. Lawson (1891) was interested in the cumulative effects of isostatic rebound on glacial Lake Algonquin shorelines. His purpose was to ascertain ”what evidences the ancient strand lines of Lake Superior afford of differential crustal movements” (p.185) and to determine whether or not the ancient strand lines ”maintained their original horizontality or had been thrown into inclined altitudes" (p.186). Lawson’s (18912186) objective was to ”secure accurate data as to the altitude of the various strand lines at as many points as possible and to gather evidence which would help to trace particular strands around the coast.” With his purpose and objective in mind Lawson wrote: [It] was clearly recognized at the outset that it would be hopeless, owing to the merely local development of shore features, to settle the question by attempting to trace out continuously any single 27 strand line or any set of such strands. It was, however, believed that the strands of different vertical series could be correlated, unless there had been rather complex movements, by carefully comparing the intervals between them (Lawson, 1891:230-231). He continued: In order to effect such a correlation the most accurate possible data were necessary, and steps were accordingly taken to ascertain the elevations of the various features indicative of shore lines by precise engineering methods. (Lawson, 1891:231) Lawson (1891) was concerned about two types of error in his work: measurement and selection. The quote above implies concern about the taking of accurate and precise measurements, yet throughout his work the author voiced confidence about his team’s ability to conduct an accurate and precise surveyz. Lawson dwelled considerably, however, on selection error, for he was very concerned about his ability to interpret relict shorelines consistently and to select points that indicate former water-levels accurately. He wrote: The figures [i.e., measured elevations -autlzor’s inclusion] do not represent always the actual level of the water. The crest of a beach 2 Unfortunately, he could not reference his survey to an absolute datum, so no way exists to assess the absolute accuracy of his measurement technique. 28 is always above the level of the water line at which it is formed; the rear of a terrace may be sometimes above the actual water line and occasionally it may be below it, so that in nearly all cases a correction has to be made for the actual water level. The amount of this correction is never definitely known so it has been thought best not to apply it, but to give the figures for the actual observations. In this connection it is well to note that many of the apparent discrepancies in strand lines are doubtless due to the varying value of this correction. (Lawson, 1891:231) Lawson was concerned about uncertainty. He tried to mitigate uncertainty about estimating former water plane positions by making observations of modern beach features encountered at the base of each topographic profile he measured. For example, Lawson (1891:231) noted, ” the base of a sea-cliff is commonly one or two feet above the calm water, though occasionally it coincides with the lake- level closely, and the rear of terraces forming at the mouths of streams is usually below lake-level." Ultimately, Lawson’s observations lead him to rely primarily on the bases of stranded sea-cliffs and lake bluffs to be the best shoreline indicators; elevations taken at the bases to be the best proxies for water levels. Prompted by debate and uncertainty in the contemporary literature, Morton and Speed (1998) sought to determine which shore features best serve as monitors of average water levels and shoreline positions. The authors studied a 29 wave-dominated and micro-tidal environment near Galveston, Texas, and measured water levels and positions along morphological features on several sandy beaches. The Great Lakes are considered to be a wave-dominated and atidal environment, so their findings have some application here. Morton and Speed (1998) argued random fluctuations in water level behavior can inhibit the taking of accurate measurements of long-term shoreline positions, but such fluctuations can be mitigated by using ”a stable beach feature” (p.1381) as an indicator. They found that, for shores composed of erodable materials, ” the stability of shoreline features increases landward and the frequency of movement of a shoreline feature increases seaward” (p.1381). For example, the toe of a wave-cut bluff is a more stable indicator of shoreline position than a berm crest or beach ridge along the same shore. The authors also report discrepancies between measured and actual water levels tend to be ” greatest on low-gradient sandy beaches along a microtidal and wave- dominated coast”, and, ” in such settings, the mean highest high water line consistently plots far seaward of the berm crest” (p.1383). Their work suggests that the average water-level, which by definition is lower in elevation than the mean highest high water line, also tends to occur far seaward of the berm crest - likely farther. Given this finding, high-gradient shorezones COmposed of materials that are less erodable than beach sand (e. g., glacial till, gravel, or bedrock) may offer better environments for taking accurate measurements of water levels and shoreline positions than do low-gradient and 30 sandy shorezones. The authors suggest using stable features like erosional bluffs or scarps, or a vegetation line as indicators. In related work, Moore (2000: p.116 after Morton, 1991) argued that ”where beach morphology permits the bluff/cliff toe or a vegetation line may be a better choice for shoreline analysis” than other beach features. The recent findings of Morton and Speed (1998) and Moore (2000) support the observations and choices that Lawson (1891) made a century before. In work on reconstructing late Holocene lake-levels, Thompson (1992) and others (Baedke et al., 2004; Johnston, 2004; and Thompson and Baedke, 1995; 1997; 1999) focused on sandy depositional features in protected post-glacial embayments. Relatively coarse foreshore deposits between successive beach ridges were used as shoreline indicators; elevations taken from them were used as proxies for water levels (Figure 2.2). Baedke and Thompson (2004) provided the most recent description of Thompson’s sampling technique, so I offer a summary of their description here. A search is made along a horizontal strandplain of beach ridges. Soil and buried SEdiments are excavated via coring at ” the change in slope that occurs along the lakeward side of the beach ridge” (Baedke et al., 2004: 440). Core site elevations are surveyed with sub-centimeter precision. Organic deposits are collected from the landward swale of the same beach ridge. In a laboratory, cores are analyzed and summarized in terms of grain size distributions, biological components, and physical structures. Abrupt and gradational grain-size changes within each 31 vertical core are used to infer the contact between relatively coarse foreshore deposits and other deposits. The vertical point of contact is assumed to be the coarsest portion of the horizontal foreshore and related to plunge-point deposits that occur ”at or very near the still-water elevation of the lake” (Baedke et al., 2004: 439). 3 Radiocarbon age determinations are made of the organic deposits. <—Lakeward Landward—’ Beach ridge Beach ridge crest crest Swale Swale watertabze Soil co're contact between coarse foreshore deposits and other sediments Figure 2.2: Baedke et al. (2004) used buried foreshore deposits as shoreline indicators. Figure adapted from descriptions in their text and from Figure 2.6 in Johnston (2004). 3 Neither Thompson (1992) nor Baedke et al. (2004) have offered evidence to suggest the location of the ”change in slope” on the land surface is an accurate indicator of the location of the coarsest portion (a.k.a. plunge-point) of a buried foreshore. Building on Thompson’s early work, Larsen (1994) documented a record of relative lake-level changes among three locations (Whitefish Point, MI; Sturgeon Bay, WI; and the Gary Airport, IN) by correlating dated beach ridge sequences. Like Thompson, Larsen studied beach ridges because they ”preserve continuous geomorphic and stratigraphic records of lake-level change” (p.108). But unlike Thompson (vis a vis Baedke et al., 2004), Larsen recognized elevation data taken from buried geomorphic features should be interpreted as possible outcomes of a fluctuating lake system and not as highly-accurate sedimentary indicators of statistical mean lake levels (+ / - 0.003 m). Moreover, Larsen (1994) claimed buried foreshore deposits usually indicate levels within the upper range of lake levels and, so, not average lake levels. When considered within the context of the range of contemporary Great Lakes water level fluctuations (Table 2.2), elevations taken from individual beach ridges could overestimate average lake-levels by 2.7 meters or more. Larsen’s (1994) argument was later supported indirectly by Argyilan et al. (2005), who claimed beach ridges preferentially reflect high water levels. In work related to this research, Schaetzl et al. (2002) reported concern about their abilities to identify relict shorelines consistently in the field and to select points that indicate former water-levels accurately. Like Lawson (1891), they recognized that different kinds of features, which are formed by different process, have differing abilities to indicate average paleowater planes consistently. For example, they reported offshore bars are typically formed 33 below a lake surface whereas a spit crest might rise several meters above it; elevations taken from these features would underestimate or overestimate average water levels, respectively. Schaetzl et al. (2002) observed that postlacustrine debris is often located on the surfaces of benches just below the bases of wave-cut bluffs and claimed, where present, bouldery/ gravelly lags provide ”the most credible estimation of mean water level" (p.405). Figure 2.3 illustrates their observation. wave—cut bluff 3/// debris (colluvium & slopewash) GPS former mean water level ./ 3" bench bouldery lag ,_... .49., wave-built terrace _' - Figure 2.3: Schaetzl et al. (2002) used bouldery/ gravelly lags at the bases of wave-cut bluffs as shoreline indicators. Figure was adapted from their Figure 2b (p.405). Where present, the bouldery/ gravelly lags at the bases of wave-cut bluffs and the foreshore deposits buried beneath beach ridges have been used effectively as geomorphic indicators of relict shoreline positions. Both are 34 effective indicators because each consists of relatively coarse—sized sediments that physically mark the sorting action of swashing and back swashing waters. Wave-cut bluffs form in erosional environments. Beach ridges form in depositional environments. In terms of surveying and mapping relict shoreline features of glacial Lake Algonquin, however, the bouldery/ gravelly lags at the bases of wave-cut bluffs seem to offer two distinct advantages. First, glacial Lake Algonquin refers to a sequence of descending lake phases that occurred during a period of ”strong (hurricane-force) winds” (Krist and Schaetzl, 2001: 11) and ”exceedingly strong waves and currents” (Krist and Schaetzl, 2001: 13), which suggests wind and wave conditions that favored erosional processes likely prevailed over gentler conditions that favored depositional processes. Consequently, I expected erosional landforms associated with Lake Algonquin shorelines, like wave-cut bluffs, to be more numerous (and easier to find) than depositional landforms associated with Lake Algonquin shorelines, like strandplains of beach ridges. And second, for practical purposes, bases of wave- cut bluffs are more accessible than buried foreshore deposits because wave-cut bluffs typically do not require trenching or deep coring to reach them in the field. Therefore, based on the long tradition of glacial Lake Algonquin research, Lawson's (1891) observations, the recent findings of Morton and Speed (1998) and Moore (2000), and the related work of Schaetzl et al. (2002), this research relies on the bouldery/ gravelly lags that form at the bases of wave-cut bluffs to be the best indicators of glacial Lake Algonquin shorelines. And, using words 35 modified from Lawson (1891), it is understood that measured elevations do not represent actual water levels, but simply the elevations of shore features ”from which the actual water levels may be inferred, within certain considerable limits of error” (Lawson, 1891: 231-232). Variation in Great Lakes water levels Contemporary Great Lakes water levels are generally resistant to change because the bodies of water are deep, have large surface areas, have small drainage basins relative to lake areas, and are controlled by relatively small outlets. Under ceteris paribus conditions, the surface of each lake is expected to rest along a field of equal gravitational potential (Clarke et al., 1990). Yet, variations in precipitation, air temperature, in-flow, and out-flow force the water surface to fluctuate and compensate, so lake gauge stations are used to capture, and statistical techniques are used to summarize, variability in water levels over time. The historical water level record provides grist for scientific inquiry. Quinn (2002), for example, analyzed seasonal variation in the record; Table 2.2 lists average lake levels and portions of Quinn’s summary (Quinn, 2002: Table 1 on page 455). Bishop also (1990) analyzed inter-annual variation and found, ” the ranges in measured mean monthly levels for Lakes Michigan-Huron and Erie are approximately 2.4 m and 1.8 m, respectively" (Bishop, 1990: 406). And Changnon’s (2004) analysis of the lake-level record for Lakes Superior, Michigan- 36 Huron, and Erie found ”record highs or record lows were being set on one or more lakes during nearly half the time over the 141—year period” (Changnon, 2004: 197). So, although generally resistant, Great Lakes water levels are not immune to change and change differently over time. Table 2.2: Seasonal variation in contemporary Great Lakes water levels (Quinn, 2002; NOAA, 2006). 4 Average lake Seasonal Standard La_k_e level (m a.s.l.) range (m) deviation m1) Time period Superior 183.5 0.30 0.11 1980-2000 Michigan 176.5 0.27 0.12 1987—2000 Huron 176.5 0.27 0.12 1987-2000 Erie 174.1 0.46 0.15 1861-2000 Ontario 74.6 0.55 0.20 1960-2000 Fluctuations around, or deviations from, long-term average water levels are generalized into three categories or scales: short-term events (i.e., short time frames that range from less than one hour to several days and, thus, considered as events); cycles (i.e., various periods), and long-term trends (i.e., multi-year). The following sections describe each scale in terms of the contemporary Great Lakes, and attempt to place them within the context of the relatively large spatial and temporal domains once occupied by glacial Lake Algonquin. Estimating a 4 Quinn (2002) did not report average lake levels, so I estimated average lake levels for a twenty-year period between January 1, 1980 and January 1, 2000 using NOAA’s Center for Operational Oceanographic Products and Services data (NOAA, 2006). Data from the lake gauge situated nearest to each outlet was summarized. likely range of glacial Lake Algonquin water-level fluctuations is important because the range can serve as a guiding threshold for discriminating evidences of different lake surfaces in the field. For example, given a field site with two instances of a relict coastal landform and a vertical distance between them that is less than the threshold value, it is reasonable to conclude, based on the observed difference between elevations alone, that both features were shaped by the same lake phase operating around the same mean water level. Conversely, given a field site with two instances of a relict coastal landform but a vertical distance between them that is greater than the threshold value, it is reasonable to conclude, based on the difference between elevations alone, that each feature was shaped by a different lake phase operating around a different mean water level. Short-term events Short—term events are the results of differential wind conditions (i.e., speed, direction, and fetch) or barometric pressures, which force localized and temporary imbalances along a water surface. Energy obtained from the atmosphere is dispersed across the water surface and delivered to the shore in waves. Wood et al. (1994) studied waves in Lake Michigan and reported most are one-meter high or less; the maximum observed was approximately six meters. Gabriel et al. (1997) and Brown et al. (2005) claim storm-generated waves and seiches can upset water levels by as much as several meters for several hours. Such short-term events effect local and temporary deviations from, but no permanent changes to, long-term average water levels. 38 When considered in the context of an entire lake surface and over a long time frame (e.g., decades), chaotic or noisy variations in lake level associated with all short-term events can be summarized over both space and time. Accordingly, the large-space and long-term distribution of water-level fluctuations associated with short-term disturbances on the contemporary Great Lakes is characterized in this work as the result of an unbiased and unstructured random process with a zero mean. Although the true variance is unknown, the quantitative information provided by Wood et al. (1994) suggests the average deviation is one meter or less. Storms and storm-generated waves occurred on glacial Lake Algonquin. Krist and Schaetzl’s (2001) treatment of Lake Algonquin type spits, for example, provides evidence of coastal erosion and longshore transport processes that were associated with ”exceedingly strong waves and currents” (p. 13) driven by ”frequent and intense cyclonic storms and their associated wind fields” (p.14). The modern Great Lakes are not commonly subjected to” strong (hurricane- force) winds” (p.11) as was glacial Lake Algonquin, so this work infers from Krist and Schaetzl’s (2001) findings that average and maximum wave heights were likely greater during glacial Lake Algonquin times than they are today. The large spatio-temporal distribution of water-level fluctuations associated with short-term disturbances on glacial Lake Algonquin can be characterized as the result of an unbiased and unstructured random process with a zero mean and an average deviation that was at least as great or greater than one meter. 39 Cyclic fluctuations Lenters (2001) analyzed the historical lake gauge record (from 1860 to 1998) for each Great Lake. His analysis of intra-annual variations revealed the influences of the annual hydrologic cycle that has been reported by many others (e.g., Fairchild, 1926; Bishop, 1990; and Angel, 1995). Peak lake levels are associated with peaks in spring run-off and precipitation, and typically occur in June (Lakes Huron, Michigan and Superior) or July (Lakes Erie and Ontario); winter lows, which are associated with precipitation lows, precipitation stored on land as snow, and evaporation, typically occur in late March or early April. And, according to Angel (1995), the annual hydrologic cycle generates fluctuations that range from 20 to 40 cm in any given year. Lenters’ (2001) analysis of inter-annual variation also revealed changes in the annual water-level cycle. Each Great Lake has exhibited a trend of increasing high spring levels and decreasing low autumn levels, although only trends associated with Lakes Erie and Ontario could be characterized as statistically significant. Quinn (2002) analyzed water level data for the period between 1860 and 2000, found similar shifts, and concluded that all of the lakes have ”experienced changes in their seasonal cycle over the past 140 years” (Quinn, 2002: 464). Todd Thompson, by himself (Thompson, 1992) and through work with others (Thompson et al., 1991; Thompson and Baedke, 1995, 1997 and1999; Baedke and Thompson, 2000; and Baedke et al., 2004), has developed a research 40 program to investigate periodicity in lake-levels in the upper Great Lakes. The program is designed to focus not only on historical lake gauge data, but on sandy depositional features in protected, post-glacial embayments. Thompson searches the local geologic record for preserved evidence of lake level history. For now, the longest history observed extends to approximately 4,700 BP. Baedke and Thompson (2000) identified quasi-periodic lake-level fluctuations that occur at 29 to 38 yr intervals with an average amplitude of 0.5 meters, and a longer-term fluctuation that occurs over durations between 120 and 180 years and with average amplitudes between 0.5 and 1.5 meters. More recently, Quinn and Sellinger (2006) found supporting evidence for such quasi-periodic extreme lake levels in their analysis of dendrochronologic data, which they linked to coincident changes in air temperature and precipitation. The record of glacial Lake Algonquin water level fluctuations is not known, but it is reasonable to expect that, then, water levels were subject to cyclic deviations forced by cyclic changes in air temperature or precipitation. Moreover, glacial Lake Algonquin was bounded to the north by a retreating ice sheet (Karrow, 1987; Farrand, 1967/ 1998). Fluctuations associated with glacial discharge and advances or retreats in ice-margin position likely occurred and, perhaps, were substantial. Tinkler and Pengally (1995) simulated outbursts from Lake Agassiz into the Upper Great Lakes and found water levels would have risen between 2 and 20 m; the largest would have required nearly three years for the Upper Great Lakes to re-attain equilibrium. Disagreement exists, however, 41 about whether or not outflow from Lake Agassiz actually entered the Upper Great Lakes during Algonquin times (see e. g., the correspondence between Karrow, 2002 and Teller, 2002). Yet, given the large size of the ice sheet and the amount of water it probably relinquished during the summer seasons, it is reasonable to speculate that such variations were greater than, perhaps much greater than, those contained within the historical lake gauge record, Quinn and Sellinger’s (2006) tree-ring record, or Baedke and Thompson’s (2000) beach-ridge record. Cyclic processes are important catalysts for local coastal erosion processes (e.g., Loope and Arbogast, 2000; Zurek et al., 2003). When considered in the global context of the large surface area and the long time frame associated with glacial Lake Algonquin, variation associated with cyclic processes can be summarized over space and time. The large spatiotemporal distribution of cyclic lake—level deviations is characterized in this work as the result of an unbiased and unstructured random process with a zero mean. Given the quantitative evidence provided by Lenters (2001) and Baedke and Thompson (2000), and given the assumption that cyclic water-level fluctuations during glacial Lake Algonquin times were greater than modern fluctuations, then the range of cyclic water-level deviations associated with glacial Lake Algonquin lake phases may easily have exceeded two meters or more. 42 Long-term trends Many studies have examined a long-term trend that impacts lake levels, which is generally thought to be caused by differential vertical movements of the earth’s crust across the Great Lakes basin (e. g., Coordinating Committee of Great Lakes Hydraulic and Hydrologic Data, 2001 and 1977; Tushingham, 1992; Clarke et al., 1994; and Mainville and Craymer, 2005). Quinn and Sellinger (1990) provide a summary of the process: ”A tremendous weight such as a regional glacier can theoretically force localized land masses downward into the crust until a compensatory equilibrium is achieved. Once that weight is removed, as in a glacial retreat, the landmass will rebound upward to once again achieve equilibrium. The Laurentide glacier of the Pleistocene Epoch depressed the crust in the Great Lakes area differentially because the thickest ice was to the northeast Due to the varying thicknesses of ice cover over the Great Lakes region, differential isostatic rebound rates are location specific.” (Quinn and Sellinger, 1990: 134) Differential isostatic rebound continues today with land at northern locations rebounding faster than land at southern locations. Lee and Southam (1994: 409) report, for example, that land at Michipicoten, Ontario is currently rising 0.525 m/ century faster than land at Duluth, MN, and land at Pt. Iroquois, MI is rising 43 0.122 m/ century with respect to Marquette, MI (from Table 1 on p.410). Ongoing crustal movement results in differential changes in the relationship between the land surface that holds water bodies and the fields of gravity that act on the surfaces of those water bodies. Consequently, isostatic rebound results in differential changes in the relationship between the elevations of contemporary shoreline features and relict shoreline features. Although differential isostatic rebound in the Great Lakes region has been well studied, little agreement exists about the behavior of the phenomenon, rates of change in the behavior of the phenomenon since deglaciation, or the amount of cumulative effect it has had on individual water planes and relict shorelines. For the purposes of discussion, however, it is important to note that the isostatic rebound process cannot be characterized as another unbiased and unstructured random process with a mean of zero; rather, as a spatially structured process with a non-stationary mean amount of uplift that increases toward the northeast. Models of isostatic rebound are presented in Chapter 5. Summary Apparent random behavior in water-level fluctuations reflects the superposition of short-term events onto longer-term cyclic processes, which interact in ways that can produce record highs and lows. In terms of the contemporary Great Lakes, the cumulative distribution of water-level fluctuations likely includes deviations from the long-term average that exceed 2.7 meters (see Table 2.3). Quinn (1999) cautions, however, that our understanding 44 of modern lake-level behavior is limited because the current regime may produce ”high and low lake levels with extremes significantly higher than our present measurements would indicate” (p4). Table 2.3: Potential range of Great Lakes water-level fluctuations (see text for references). Process and duration Typical deviation (m) Storms or waves ~from seconds to several hours :10 Seasonal hydrologic cycle ~1 year 0.2 Baedke and Thompson (2000) ~33-year cycle 0.5 ~150-year cycle 1.0 Differential isostatic rebound (varies with location) ? Total i 2.7 m The range of historic water-level fluctuations is important because it suggests the precision with which average water levels (statistical quantities subject to variation) can be estimated from spot elevations (measured quantities) taken from static landforms. As explained above, this research relies on the gravel lags that form at the bases of wave-cut bluffs to be geomorphic indicators of shoreline events because they tend to form at or near the land-water interface and where wave energies are relatively high. Elevations taken from the tops of these lags serve as proxies for average water levels (Schaetzl et al., 2002). Because beach sediments and shoreline features can move in response to water level fluctuations, I must recognize that estimates of glacial Lake Algonquin water 45 levels based on geomorphic proxies are subject to uncertainty. Fortunately, the range of contemporary Great Lakes water-level fluctuations (Table 2.3) provides a convenient (albeit rough) estimate of this uncertainty. This research proceeds by assuming the minimum vertical difference between two gravel lags that can be used to discriminate two average glacial Lake Algonquin water levels is 3 meters. Summary Lake levels and shorelines are results of complex and dynamic interactions among bodies of air, land and water, and gravitational forces. Such interactions can occur over multiple spatial and temporal scales. Changes in lake level, as well as rates of change in lake level, are generally recognized as important controls on shoreline position, coastal processes, and coastal landform changes. In general, shorelines migrate upward and landward as water levels increase, and migrate downward and lakeward as water levels decrease. Statistical techniques are usually required to summarize the distributions of historic or contemporary lake-level fluctuations in terms of average lake levels. Researchers interested in mapping and modeling paleolakes and paleoshorelines must deal with difficulties associated with estimating statistical lake levels by taking elevation measurements from relict shoreline features. The precision of this method can be influenced by the choice of shoreline indicator as 46 well as the legacies of earth processes that have modified the landscape or moved relict shoreline features since they were formed. This research relies on the bouldery/ gravelly lags that formed at the bases of relict wave-cut bluffs to be indicators of former shoreline positions (Schaetzl et al., 2002). It is understood, however, given natural variability among local site conditions, different lags that formed along the edge of the same lake surface might have elevations that vary. The tradition of glacial Lake Algonquin research and the later work of Morton and Speed (1998) and Moore (2000) suggest the precision with which a lake level can be established by the elevation of a wave-cut bluff alone is likely sub-meter, but knowledge about contemporary water-level behavior indicates the precision may be, at the extreme, no better than three meters. Accordingly, this research recognizes imprecision as a form of uncertainty and proceeds by assuming the precision with which a former lake level can be established by the elevation of a gravel lag situated at the base of a wave—cut bluff is one meter, and the minimum vertical distance that can be used to confidently discriminate different lake levels is three meters. 47 CHAPTER 3: METHODS AND FIELD CAMPAIGNS Introduction This chapter describes methods used to collect relict shoreline data and build a glacial Lake Algonquin dataset. I begin by summarizing the work of Schaetzl et al. (2002), which established a relevant methodology and produced an initial dataset, which I inherited. I continue by recognizing opportunities to improve upon the work they had done. Next, I review the prior mapping campaigns led by Cowan (1985) and Karrow (1987) and describe how data they collected were conflated with data collected during this work. I then report new data I collected during a series of field campaigns designed to expand and improve the geographic distribution of sampled shorelines. My report of relict features at Cockburn Island, Ontario, for example, is a significant contribution to the literature for no other report describes these features in detail. And, in light of new data and recently published work, I offer reinterpretations of old data. The initial glacial Lake Algonquin dataset Uncertainty about the elevations, timing and extents of glacial Lake Algonquin lake phases prompted Schaetzl et al. (2002) to conduct an extensive survey of relict coastal landforms in northern Michigan. They targeted more than 200 sites and were guided by descriptions of relict landforms in the literature (e.g., Leverett and Taylor, 1915; Futyma, 1981; Farrand and Drexler, 1985; Schaetzl et al., 2002; and Krist and Schaetzl, 2001) and their interpretations 48 of features on 7 .5-minute USGS topographic maps. The authors investigated sites to determine whether or not each one contains a relict wave-cut bluff suitable for sampling. At suitable sites, Schaetzl et al. (2002) measured bluff positions, rather than accept any previously published elevation values, because the literature provides: a) few detailed descriptions about where previously reported position measurements were actually made; b) few references to the spatial reference systems employed (i.e., either the horizontal or vertical datum was unspecified); and c) no assessments of data quality. Without such metadata, the authors felt they could not assess confidently whether or not any previously reported elevation value was valid, commensurate with other published data, or acceptable for inclusion in further analysis. Their dataset contains 146 records (Figure 3.1). Individually, each record represents the position of a relict shoreline feature. In aggregate, their data reveal the spatial extents of four, possibly five, glacial Lake Algonquin lake phases in northern Michigan (Schaetzl et al., 2002). Schaetzl et al. (2002) concluded that pro- and postglacial lake levels can be effectively identified, reconstructed, and correlated via examinations of elevation data collected from wave-cut bluffs. They claimed their rebound curves and isobase maps are more accurate than those published previously, but suggested more geomorphologic and stratigraphic research is needed beyond their area of study. 49 if I l . > The Main Lake Algonquin shoreline was -- Current shorelines i ,/ not generated for Canada's mainland. Show," of mm <5; Lake Algonquin >§ ONTARIO ' GPSpointnaardod D- n Town 2 g :- Univ «Michigan u Biologic Sutton in Sugar Island 5’ -‘ ’ Neebish Island . r: Approximate ice. margin position State or province bou - PENNSYLVANIA Figure 3.1: The study area explored by Schaetzl et al. (2002), points in the original glacial Lake Algonquin dataset, and the extent of the Main phase of glacial Lake Algonquin as they had determined it. Figure adapted from their Figure 1. 50 Opportunities to improve the glacial Lake Algonquin dataset Several opportunities existed to substantially yet efficiently improve the initial dataset built by Schaetzl et al. (2002); I capitalized on three. First, I expanded the study area (see Figure 3.2) to include relict features at Sault Ste. Marie (Cowan, 1985), St. Joseph Island (Leverett, 1913 and Karrow et al., 1987) and Cockburn Island, Ontario (Leverett, 1913; Chapman and Putnam, 1966), as well as others on Beaver Island, Michigan (Dietrich, 1978). Relict bluff features at Sault Ste. Marie, Ontario and St. Joseph Island are important because they are well-known and uniquely situated between Lakes Huron and Superior; they link the evolutionary histories of these lakes. The correlations Cowan (1985) and Karrow (1987) made between glacial Lake Algonquin phase names and stranded shoreline elevations can be used to revise correlations made by Farrand and Drexler (1985) and Schaetzl et a1. (2002) in northern Michigan. Second, the initial geographic distribution of sampled positions (Figure 3.1) is not oriented along the generally accepted direction of differential isostatic uplift (approx. N 15°E; after Futyma, 1981), but obliquely to it (approx. N 35°W). The consequences that such an orientation has on the process of estimating rebound trend parameters are not known, but it seems reasonable to treat any such parameters as provisional until additional data can be collected along, rather than normal to, the uplift trend. Collecting and integrating data from the islands mentioned above adds important features to the glacial Lake Algonquin dataset and much needed support to unsampled sectors along the former water planes. 51 Sault Ste. Mme. ONTARIO CHIPP-EWA Kincheloa 1 'O ' " n scon MACKINAC "LEM; 5 1M Hansel PRESQUE ISLE \ l 1 ”his" -... .,.,._...._.,.. l “N 1 CHEBOYGAN 1 CHARLEVOIX L._...._-.~..._—.~.—w—a-—..._. .M..,.. E MONTMORENCY ; Ii 073560 y ALPENA ANTRlM 1 KALKASKA . CRAWFORD O-SCOD‘A Ham ALCONA 1 nmvsnss i ngpono l MISSAUKEE ROSCOMMON Jr a OGEMAW l loscoJ fig\ i- ..‘_.....4...‘._,....._... Figure 3.2: Map of the current study area with names of places and features mentioned in the text. County boundaries and names added for reference. 52 Last, Schaetzl et al. (2002) and Drzyzga et al. (2002) used an undersized set of ground control points (n=8) to assess the quality of sampled shoreline elevations (n=146). Drzyzga et a1. (2002) made simplistic assumptions about the distribution of error in the dataset, which may or may not represent the actual distribution of error. So, I collected a larger set of ground control points (n=24). The larger set of ground control points is useful for better assessing the quality of the glacial Lake Algonquin data set and mitigating the effects of an apparent error process (e.g., a systematic bias or spatial pattern). The accuracy assessment is reported in Chapter 4. Shoreline identification, bluff selection, and sampling methods Frank Taylor (1895) stressed the value and importance of studying the ancient islands of glacial Lake Algonquin. Notable examples include the ancient Munuscong Islands in Mackinac County, Michigan (Taylor, 1895) and proto-St. Joseph Island (Karrow, 1987), proto-Cockburn Island (Chapman and Putnam, 1966), and proto-Beaver Island (Dietrich, 1977). Accordingly, I focused my attention and fieldwork on contemporary landforms that may have existed as islands in glacial Lake Algonquin and tried to collect data from as broad an area as possible. I inherited the dataset and employed the same set of shoreline identification and bluff selection methods used by Schaetzl et al. (2002: 405-406). And, like they did, I relied on the bouldery/ gravelly lags that form at the bases 53 of wave-cut bluffs to be geomorphic indicators of ancient shorelines (Figure 3.3). Using the same indicator ensures similar reliability and consistency characteristics between data collected during their campaign and new data collected during mine. Similarities notwithstanding, I provide a summary of methods below because several differences do exist between the GPS leveling technique and data processing methods I used and those that Schaetzl et al. (2002) used. For example, I take advantage of the latest geoid model (GEOID03, Roman et al., 2004) that is a better approximation of mean sea level throughout the United States than the GEOID96 geoid model (Smith and Milbert, 1999) they used. wave—cut bluff f/ debris (colluvium & slopewash) GPS former mean water level . “fig-”WM b h wave-bunt terry , _ // enc ./’/ / (”11-7" _ .6532’”~ bouldery lag Figure 3.3: Schaetzl et al. (2002) used bouldery/ gravelly lags at the bases of wave-cut bluffs as shoreline indicators. Figure was adapted from their Figure 2b (p.405). 54 First, I collected 7.5-minute USGS topographic map sheets for all locations within the study area. I also collected 1:50,000 Canadian Topographic Map sheets to supplement coverage of Cockburn Island, St. Joseph Island, and Sault Ste. Marie, Ontario. Field notes taken during previous campaigns as well as extant literature (e.g., Leverett and Taylor, 1915; Hough, 1958; Dietrich, 1978; and Karrow, 1987) were used as guides for targeting likely sites with relict wave-cut bluffs. Initial selections were limited to only those features that are landward of, and higher than, the generally accepted N ipissing water surface (Goldthwait, 1910). I plotted extant data and new potential bluff targets on topographic map sheets to identify areas with little data. Next, I visited targeted sites. A target was selected for sampling if it exhibited a conspicuous bluff-to-beach morphology and presented a well- defined bluff toe. If a target exhibited a conspicuous bluff-to-beach morphology but not a well-defined toe, then I traveled along the feature (for up to several hundred meters in either direction of the initial site) in search of a ”best” site that contained a well-defined toe and a relative paucity of forest cover (to facilitate radio signal reception from orbiting GPS satellites). Locations that evinced anthropogenic change or bluff failure are not suitable locations for sampling and, so, were not sampled. Third, an initial position was selected along the bluff toe. Soil was excavated with a bucket auger to a maximum depth of 1.5 meters. If a bouldery/ gravelly lag (Figure 3.3) could not be detected, then sampling 55 proceeded down slope at successively lower locations until a buried lag was found. (Fortunately, at many sites a conspicuous boulder lag is present below the bluff.) Field notes were taken at the ultimate location to describe the site and record depth to the buried gravel lag. Elevation measurements, which were taken at the land surface, were adjusted downward to account for the vertical distance between the land surface and the buried lag. Also, when weather permitted, an image of the site was recorded with a digital camera. And fourth, real-time differential Global Positioning System (GPS) technology was employed to measure relict shoreline positions. When used correctly, GPS leveling offers at least five distinct advantages over other field methods: a) on average, position measurements {x,y,z triplets} can be made quickly and with sub-meter errors; b) position values collected among scattered field sites can be georeferenced to the same spatial reference system; c) field equipment is light-weight and portable; d) fieldwork can be conducted effectively and efficiently by one person, which precludes the need for a second person to work the other end of a level line; and e) all position data are recorded and stored in digital form, which facilitates translation into data formats that can be analyzed within a (313 or spreadsheet environments. I planned and organized data collection activities in three important ways to take best advantage of GPS technology. First, the field season was limited to leaf-off conditions (i.e., late fall) because summer-time leaf canopies can obfuscate, perturb, or block GPS radio signals. The winter leaf-off season was 56 also avoided, however, because regional winter snow blankets ground features and inhibits exploration. In the field, the GPS receiver was mounted onto a tripod and raised two meters above the ground. Second, all data were collected during temporal windows that framed optimal GPS satellite constellation configurations (i.e., VDOP S 3.0) 1. GPS receiver software can manipulate and tumble numbers under a plethora of constellation configurations, but relatively few configurations offer geometries sufficient for taking accurate and precise vertical measurements. Trimble’s Quick Plan software (Trimble Navigation Limited, 2002) and current satellite ephemeris data were used to identify and select suitable VDOP windows. Data collection activities were scheduled and limited to these windows only. And third, a point-averaging technique was used to estimate position at each site. The point-averaging technique requires 1 A Dilution of Precision (DOP) value is an estimate of the quality of a GPS position. For a given point in time, the DOP value takes into account the location of each GPS satellite relative to the locations of others in the constellation, and their geometry in relation to the GPS receiver (Reilly, 2002). A low DOP value suggests a relatively high probability of taking an accurate measurement. A low VDOP (Vertical Dilution of Precision) suggests a relatively high probability of taking an accurate vertical measurement. Because this research is concerned with taking accurate vertical measurements, primary emphasis is placed on obtaining low VDOP values during data collection activities. Most satellite-viewer geometries that are optimal for taking accurate vertical measurements also allow taking accurate horizontal measurements. 57 multiple measurements of the same feature to be made. Collecting a sample of measurements and summarizing it in terms of average coordinate values and precisions, rather than collecting a single measurement, reduces the risk of accepting data with large random errors. Schaetzl et al. (2002) used the point averaging technique and collected a sample of no less than 100 position measurements at each bluff site. They achieved reliable results with the technique, but had to discard some valuable data because, for some sites, many of the 100+ measurements were deleted during data processing (pers. comm.) Collecting a larger sample at each site would have helped preserve their data. So, for this work, I collected a larger sample of no less than 300 position measurements at each bluff site. All horizontal locations are referenced to the North American Datum of 1983 (N AD83). All elevations were referenced first to the World Geodetic System 1984 (WGS84/ C3873) ellipsoid model and later converted to orthometric heights (see below). Data processing Data files (one sample of shoreline measurements per file) were processed using Trimble Pathfinder Oflice (GPS) and Microsoft Excel (spreadsheet) software. Data processing followed this model: First, each file was assigned to one of two groups: those with position data that were differentially corrected in the field, and those that were not. The radio link between a GPS receiver and the satellite that delivers OMNIStar’s real-time 58 correction service is sometimes lost behind local obstructions (e. g., a massive bluff or a dense forest condition). Each file with uncorrected data required post- processing using differential-correction data made available via the network of National Geodetic Survey, Continuously Operating Reference Stations (CORS). Next, the spatial distribution of individual position measurements collected at each bluff was examined using a scatterplot matrix. The 3-D distribution was examined and filtered of positions associated with conspicuous multi-path errors (i.e., a linear pattern of sequential measurements that drifted over distances greater than one meter and suggested the incoming radio signal was bouncing off a nearby object) or indicators of poor satellite geometries (i.e., VDOP values > 3.0). The remaining set of acceptable positions was summarized geometrically to support average x, y and z coordinate values and measures of precision (i.e., standard deviation). Elevation values were processed one step further to convert heights referenced to the WGS84 ellipsoid to orthometric heights, which are referenced to the North American Vertical Datum of 1988 (N AVD88) and tied to the familiar notion of mean sea-level. Orthometric height (H) is calculated as: H=h~N [Eq. 1] where: h is the average of sampled heights above the WGS84 ellipsoid as obtained from GPS leveling; and N is the geoidal undulation at that location and 59 obtained from the GEOID03 model (Roman et al., 2004). In Michigan, users of GEOID03 are able to convert GPS-measured ellipsoid heights to orthometric heights to within 2.7 cm of traditional land-based geodetic measurements (Roman etal., 2004: 11, Table 1). Last, field notes, digital images, and feature positions were reviewed a final time. Each elevation value was adjusted downward to account for the tripod height (two meters) any distance between the land surface and the buried gravel lag. Feature positions associated with conspicuous wave—cut bluffs were added as records to the glacial Lake Algonquin dataset and exported to an ASCII text file, which is readable by common GIS or spreadsheet applications. Feature positions not associated with conspicuous wave-cut bluffs are retained in a supplemental database. This research employs the same sets of selection criteria, GPS equipment (T rimble XRS receiver, antenna and software) and real-time differential correction procedures used by Schaetzl et al. (2002). So, by following the same data collection methods and using the same field equipment, it is reasonable to assume that the data collected during this exercise can be conflated with data collected during the previous work. Yet, to ensure data are conflated properly in every manner, I utilized the original and unprocessed set of GPS data files collected by Schaetzl et al. (2002), which I inherited, and processed each according to the model describe above. The only notable difference in methods was that Schaetzl et al. (2002) employed the GEOID96 geoid model (Milbert and 60 Smith, 1996) and this work uses the current GEOID03 geoid model (Roman et al., 2004) to convert heights measured above the ellipsoid to elevations tied to the N AVD88. The current geoid model is considered the most accurate available. Summary of prior mapping campaigns Sly and Lewis (1972), Cowan (1985), and Karrow (1987) collected paleoshoreline data at Manitoulin Island, Sault Ste. Marie and St. Joseph Island, Ontario, respectively, and correlated their features with Algonquin lake phases named by Deane (1950: Ardtrea, Upper Orillia and Lower Orillia), Stanley (1936: Wyebridge, Penetang, Cedar Point and Payette) and Hough (1958: Sheguiandah and Korah). Their campaigns relied on the bases of wave-cut bluffs to be indicators of former lake levels and they employed survey-grade equipment to measure paleoshoreline elevations precisely and accurately. Moreover, features at Sault Ste. Marie and St. Joseph Island are situated near others surveyed in Northern Michigan, so I used them to inform my classification of relict features in Northern Michigan and, in effect, revise some previous classifications proposed by Farrand and Drexler (1985) and Schaetzl et al. (2002). Cowan (1985) mapped Quaternary sediments at Sault Ste. Marie, Ontario, and concluded the sandy terrain situated between the exposed Canadian Shield (to the north) and St. Marys River (to the south) exhibits a descending sequence of stranded shoreline features. I-Ie associated the uppermost (at 309 m) with the Main phase of glacial Lake Algonquin. Although Cowan reported elevations, he 61 did not report the locations of features he surveyed, so readers of his report cannot use it as a field guide. I contacted Dr. Cowan to inquire about the locations of bluff features he surveyed and, graciously, he sent me copies of his field notes and data. He also sent copies of topographic maps, on which he marked the locations of benchmarks, level lines and relict bluffs. Given this wealth of information, I was able to combine his surveyed elevation data (2) with corresponding location data (x,y) that I interpolated from each topographic map. According to Cowan (1985), his elevation data are accurate to within 30 centimeters. The location data I interpolated from topographic maps are accurate to within 25 meters. Table 3.1 summarizes Cowan’s (1985) results (Table 1, p34). Table 3.1: Approximate elevations for glacial Lake Algonquin water planes and the Nipissing water plane at Sault Ste. Marie and Manitoulin Island, Ontario (Cowan, 1985). Elevation at Elevation at Sault Ste. Marie (m) Manitoulin Island (m) Lake phase 309 309 Main 302? - Ardtrea 295 287-296 Upper Orillia ? - Lower Orillia 275? 271-276 Wyebridge 265 267 Penetang 257 259 Cedar Point 247 244 Payette 233i: 233 Sheguiandah 210 212 Korah 198 198 N ipissing 62 Karrow (1987) mapped Quaternary sediments at St. Joseph Island, Ontario and surveyed stranded shoreline features. Specifically, he measured the elevations of wave-cut bluffs and the tops of beach bars. He projected his surveyed elevation data onto a line of profile oriented along the presumed direction of maximum uplift (N15°E; after Futyma, 1981). Karrow compared his projected elevation data to the set of previously named water planes (Hough, 1958; Sly and Lewis, 1972; and Cowan, 1985) and concluded, despite some inconsistencies, that all but the central hill of St. Joseph Island was submerged by the Main phase of glacial Lake Algonquin. I summarize Karrow’s Figure 5 (p116-117) as Table 3.2 to facilitate comparison with data I collected in this work. Table 3.2: Lake level elevations for glacial Lake Algonquin water planes and the N ipissing water plane at the southern end of proto-St. Joseph Island, Ontario (Karrow, 1987). Elevation (m) Association 280 Main Algonquin 27 2? Ardtrea 266? Upper Orillia 258? Lower Orillia 240? Wyebridge 232 Penetang 225 Cedar Point 212 Payette 200 Sheguiandah - Korah 196 Nipissing Readers of Karrow’s (1987) report cannot use it as a field guide, so I contacted him to inquire about the locations of the bluff features he surveyed. 63 He graciously sent me copies of his survey data and marked topographic maps with the locations of relict wave-cut bluffs he surveyed. Given this information, I was able to combine his surveyed elevations (z) with corresponding location data (x,y) that I interpolated from topographic maps (Figure . In related work, Karrow (1988) indicates his surveying method yield elevations with decimeter- level accuracy. The location data I interpolated from topographic maps are accurate to Within 25 meters, which is the precision afforded by the width of a printed contour line at the 150,000 scale. Mapping campaigns conducted for this work Frank Taylor stressed the importance of ancient islands when he wrote: ”. . . [they] can hardly be overestimated, especially if they are situated well out from the mainland. Each one furnishes a new point of support for the restoration of the water plane and must prove valuable in the ultimate study of the earth’s history as disclosed in the deformation of former water levels" (Taylor, 1895: 28) Accordingly, I focused much of my attention and fieldwork on contemporary landforms that may have existed as islands in glacial Lake Algonquin. I visited sites with potential bluff features between June 2000 and December 2003. The subsections that follow document features and observations I made during these visits. In the text I provide maps and refer to field sites using a system of unique 64 site identification numbers (e. g., SID642) and report features elevations in meters. Figure 3.4 can be used to reference each of the smaller—scaled area maps. Last but not least, I report correlations I made between elevations taken from relict wave- cut bluffs and named glacial Lake Algonquin lake phases. 65 ONTARIO l l lPRESQUEISLE l ......_........i. . .._... . m... . .. .1 ’T T M» j CHEBOYGAN l . .I ~ 1 . 1 j ‘- CHARLEVOIX l r: I 1 MONTMORENCY ’ 1 078560 . 1 V ANTRIM l l ‘ I 1 1 ‘ I ______ ,1t~_~_,. l l l L, l i 1 ¥ KALKASKA 1 CRAWFORD 1 OSCODA l J ; j ' * ‘ l l l 1 ‘GRAND , j I ‘ TRAVERSE 1 l l 1-1.--_m..-- ”M V; I - F“ o m. .I-I - NA“: V 'N‘” > "m h‘“ ‘V -- i ‘OSCO cum 3‘.“ r32 WEXFORD l MISSAUKEE ‘ ROSCOMMON f OGEMAW .1 Figureenfifgg Figure 3.4: Index map for smaller-scaled maps that show field campaigns and field sites mentioned in this chapter. Cockburn Island, Ontario Leverett (1913) and Chapman and Putnam (1966) mentioned that a possible islet of glacial Lake Algonquin exists at Cockburn Island, Ontario, but neither they nor subsequent writers have offered a detailed description of the supposed islet or data to support their observations. P. F. Karrow (personal communication, 2003) surveyed a line in 1991 across the northern part of Cockburn Island to the supposed islet, but his data remain unpublished as part of a larger study of nearby Manitoulin Island. I conducted fieldwork on Cockburn Island during October, 2003. Transportation to and from the Island, lodging, and access to an off-road vehicle were provided by Huron Timber Company of Thessalon, Ontario, which owns more than 80 percent of the ~100 km2 island. Cockburn Island has no permanent residents, paved roads, or services of any kind; the sole town of Tolsmaville evinces former life as a port community but is now defunct. Concession roads and logging trails offer limited access to interior portions of the island, but accessing most relict features must be achieved on foot. I encountered Mr. Harold McQuarrie, a seasonal island resident, by chance. He maintains a hand- drafted map of island features that depicts logging trails, property ownership boundaries, concession- and patent-holder names, and basic hydrographic features. Mr. McQuarrie's map (McQuarrie, 2000) greatly facilitated my efforts to traverse the Island and approach bluff targets. Figure 3.5 shows Cockburn Island and identifies field sites visited and discussed in the text. 67 5153:27. . ‘ " 152: 25.27%?” 316: 278.4m ' .<.-'137'272.9m ~ - " 36274.1": Figure 3.5: Relict bluff sites at Cockburn Island, Ontario. Sites are labeled with a unique site identification number and an elevation value (e.g., 642: 200 m). I focused attention on a small upland called McCaigs Hill. McCaigs Hill is the highest feature on Cockburn Island, approximately 2 km2 in size, composed of unconsolidated sediments, and divided roughly into quadrants by two intersecting concession roads. McCaigs Hill exhibits a relatively smooth surface that dips visibly southward; only a shallow saddle, which trends NNE-SSW, and scattered well-rounded boulders interrupt the surface. Soil samples taken along the surface (to a maximum depth of 1.5 meters) revealed red-brown sands and 68 few stones; clay content tended to increase with depth. Much of the upland surface and some of the eastern and southern slopes were once cleared of forest cover and are now covered by grasses; relict landforms in these areas are exposed plainly and clearly. A small area in the northeastern section of McCaigs Hill is known locally as Scotch Block, for a settlement of Scottish families that once inhabited the site. The steep sides of McCaigs Hill exhibit strong bluffs at four distinct elevations (280, 272, 266, and 259 m); two can be traced around the Hill, but all are interrupted in places by localized bluff failures or road cuts. The two highest (~280 m: SID128, 136, 141, 144, 152 and 154; and ~272 m: SID127, 137, 138, 142, 153 and 155) are conspicuous and well-defined (Figures 3.6 and 3.7), particularly in the eastern side of the Hill. Two lower sets of bluffs (at ~259 m: SID126 and 157; and at ~266 m: SID145 and 156) are also conspicuous and well-defined, but on the northwestern side of the Hill only. The bluff at 259 m (height >16 m) is well preserved, but shrouded by dense tree cover along most of its length. All other bluffs are absent where this bluff is strongest; it visibly truncates the next higher bluff (266 m) and, perhaps, an intermediate bench. The bluff changes character rapidly, for it shortens into a low bluff and then degrades into three wave-cut scarps or bars that run in a south-southeasterly direction away from northwestern side of the upland. The bars are composed of rounded gravels and are relatively tall, narrow and compactly spaced at locations near the bluff site, but tend to shorten, widen and separate as they curve away from the Hill. 69 Figure 3.6: Image of a relict bluff on the south side of McCaigs Hill, Cockburn Island, Ontario. lake bottom Figure 3.7: Diagram of the relict bluff shown in Figure 3.6. Shoreline positions were measured at gravel lag situated at the base of the bluff. 7O A pair of lower shorelines is also present (at ~230 m: SID132, 146, 149 and 150; and at ~222 m: SID131, 133, 147, 148, 151 and 158). They are situated several kilometers from McCaigs Hill and outline what must have been younger and larger proto—Cockburn Islands. I characterize these two shorelines as a pair because each occurs in close proximity to the other. For example, at some locations these features present as a set of gravelly bars, with one situated several meters higher than the other. At other locations they appear as a combination of erosional bluff and bar or as a set of erosional bluffs. Incidentally, this pair is the highest set of relict shoreline features connected to bedrock outcrops. I found small dolomite outcrops at the northern and southern tips of this ancient island which might have supplied the parent material for the gravelly bars. Cedar swamps sandwich narrow hardwood stands situated along the sloping sands and gravels between the bluffs at the 222 m and 230 m elevations and might hold organic material suitable for radiocarbon analysis and dating. Lumbermen were harvesting hardwoods in this area and during the time of fieldwork. Their safety concerns prohibited me from accessing features at these worksites, so I consider the survey of relict features as incomplete. Relict bluffs at Cockburn Island have elevations that correspond with the elevations of relict bluffs stranded at Sault Ste. Marie, St. Joseph Island, and Manitoulin Island as well as others in the eastern half of Michigan’s Upper Peninsula. Table 3.3 lists the elevations of relict bluffs at Cockburn Island and the names of coplanar water levels at nearby St. Joseph Island. 71 Table 3.3: Approximate elevations of glacial Lake Algonquin water planes measured at Cockburn Island, Ontario. Approximate elevation (mL Water plane / lake phase 280 Main 272 Ardtrea 266 Upper Orillia 259 Lower Orillia 230 Penetang 222 Cedar Point Cowan (1985) describes relict shoreline features at Sault Ste. Marie, Ontario and argues they are directly comparable with others at Little Current, Ontario (Sly and Lewis, 1972), because respective interpretations of both "Main Algonquin and Nipissing levels are within one meter of each other" (p.36). If Cowan's argument is correct, then both locations have experienced similar rates of, and cumulative amounts of, isostatic recovery over time. Cartographically, an isoline of paleolake-plane elevation can be drawn to represent the apparent connection between these two locations (Figure 3.8, 309 m). Such a line is oriented similarly to other isolines drawn by Goldthwait (1910:38), Leverett and Taylor (1915:429), Hough (1958256) and Clark and Persoage (1970). Cockburn Island is situated south of the isoline drawn between Sault Ste. Marie and Little Current, Ontario (Figure 3.8). Given the generally accepted trend of increasing isostatic rebound across the region (N15°E; Futyma, 1981) and given the relative location of Cockburn Island with respect to this isoline, it is reasonable to expect that intersections between glacial Lake Algonquin water 72 planes and the Cockburn Island land surface occur today at elevations lower than those measured at Sault Ste. Marie and Little Current, ON (309 m, after Cowan, 1985: p. 34). The highest shoreline measured at Cockburn Island is at 280.2 m, which is lower than the highest shorelines found at Sault Ste. Marie and Little Current, and proximate to the highest shorelines found at St. Joseph Island, ON (280 - 285 m, after Karrow, 1987) and in northern Michigan (282 - 285 m, after Schaetzl et al., 2002). The highest shorelines at St. Joseph Island and in northern Michigan are associated with the Main phase of glacial Lake Algonquin, so it seems reasonable to infer from the near-equivalent elevation of the highest shoreline at Cockburn Island that McCaigs Hill, a small upland, was an ancient islet of the Main phase of glacial Lake Algonquin. 2 I Y 1 1" am 32w r i l U) ,_ . . Al Sault '1 T i ' Lake " , ' ' ' Supenim‘ ( K/ " \ fl \ ”fr ;\ . » 'x_ if 1 .x// \i_/"/ ‘ l M 1 7. \ b \_ r‘ ~.\\'/A“’;/ ,’ \.( _.., ‘35“ . l"\,_ 44‘ A, NJ; ‘1 1' "' -... ‘ 4“ 3m\ 3 DDN \s N \/ “u x ’5" Jim. / bx... -/‘ m!) M ,’ i [bx ( I. ‘ Mackinac “ "~'Ex.-..-»\,\L.\ ‘4 ;'/.I ~ 1: g l W, \N Qflsland Cockburn mi. 1. w}! i . l") “a, x'ri‘m- ’ / '1 Michigan ,3 “\\“‘3 ISIand . _L".-., /,.-j / w—w} \lw' L04. Manitoulin liuck j}? V./\.\ (3 filll‘o’] Island \ K ’ 2 saw 32w ‘~/ 1 Figure 3.8: Apparent contours across the Main Lake Algonquin water plane. 73 Interpretation of lower shorelines is problematic, however, because they cannot be correlated consistently with water planes named at nearby places. For example, Schaetzl et a1. (2002) argued that only the Main, Ardtrea, Wyebridge and Payette phases can be interpreted consistently among the data they collected in Michigan. This set of four lake phases is not sufficient for interpreting the six water planes I observed at Cockburn Island. Based on their isoline maps of bluff elevations (Schaetzl et al., 20022411) bluffs at Cockburn Island that occur at 280, 272, and 266 m correlate closely with their Main, Ardtrea and Wyebridge phases, respectively, the bluffs stranded at 259, 230, and 222 m have no correlates. Karrow (1987 :118) reports many relict shoreline features on St. Joseph Island that appear to correspond to Upper Group and Lower Group beaches, but he claims his data are insufficient ”to clearly distinguish particularly consistent and well- developed water planes" below the highest one. Yet, assuming St. Joseph Island is to Cockburn Island as Sault Ste. Marie is to Manitoulin Island, I use Karrow’s lake plane profiles (Figure 5, p116-117) to extrapolate named water levels from St. Joseph Island to Cockburn Island. In descending order, paleoshorelines at Cockburn Island correspond with the Main, Ardtrea, Upper Orillia, Lower Orillia, Penetang and Cedar Point water planes. Sugar Island, Michigan Schaetzl et a1. (2002) surveyed two relict features at Sugar Island, which is situated in the St. Marys River and between Sault Ste. Marie and St. Joseph 74 Island, Ontario (Figure 3.9). The authors measured a high feature twice (240 m: SID395 and 407) and correlated it with the Payette phase. They measured a low feature twice (228 m: SID396 and 453) but correlated one measurement with Lake Minong phase of Lake Superior and attributed the other to an unnamed ”lower” phase of Lake Huron. In light of data available in nearby Canada, I believe the lower feature is best correlated with the Sheguiandah phase. L: 1’0 r3 v ‘1 Figure 3.9: Relict bluff sites at Sugar Island, Michigan. Sites are labeled with a unique site identification number and an elevation value (e.g., 642: 200 m). Also shown are the locations of features surveyed by Cowan (1985). 75 Rexford, Michigan Main Algonquin occupied a ”sizeable part of Whitefish Bay in eastern Lake Superior” (Cowan, 1985: 34). Farrand and Drexler (1985) recognized several post-Algonquin lakes, ”at least through Sheguiandah” (p30), that transgressed into the Superior basin; they correlated features found near Pendills Bay with the Wyebridge and Payette phases. A careful read suggests, however, that Farrand and Drexler (1985) were not fully confident in their assignments of Algonquin phase names to the features they studied (i.e., ”Wyebridge(?)”, p29). In other work, Futyma (1981) identified Main and Battlefield (now deprecated) type features in this same region, but did not correlate features having intermediate elevations with named lake phases. Schaetzl et a1. (2002) surveyed 16 features in western Chippewa County, MI (Figure 3.10) and based their correlations on those named previously by Farrand and Drexler (1985). The highest stranded bluff in this region, located 3 km due east of McNearney Lake, has a toe elevation of 281.9 m (SID393) and is correlated with Main Lake Algonquin; no higher bluff exists in this region. Two nearby features were measured at 262.1 m (SID392) and 259.5 m (SID404) and the authors assigned them to the lower Wyebridge phase. Surprisingly, Schaetzl et a1. (2002) associated features near Dollar Settlement and stranded at ~232 m with Lake Minong. Farrand and Drexler clearly describe a bench near Dollar Settlement ”at 213 to 219 m" to be ”a good candidate for the Minong shoreline” because it occurs in the Pendills Lake embayment immediately west of Nadoway 76 Point and does not occur ”to the east on the Lake Huron slope” (Farrand and Drexler, 1985: 29). Confusion about these features prompted me to reinterpret W111 tcfz sh B a y 1 37127062111 3912mm. , ' . 3741231 .3m .. 392261.111]. . 403:230.2m;w,.- 372:259.7m ‘* ‘- .. : $251311»- .11'0432595m P 1...: _ I. ’4 _F f':_;:,:f 7’ . . .. .= . ., i: ' '1 :5 48153551111? 3 5 A‘ , . 406.272.2111; ‘ . a”: . '. Rexford f . =—‘“‘“ ' - - . i. .7... 5" 4552597111" ,. .3 "7 " . , 1" ._r' . ’9; ;_, ’ Eastern, ‘ Figure 3.10: Relict bluff sites near Rexford, Michigan. Sites are labeled with a unique site identification number and an elevation value (e. g., 642: 200 m). The highest stranded shoreline in this region has an elevation of 281.9 m. Nearby features associated with the Nipissing water plane rest at 198 m. The elevations of Main Algonquin and Nipissing features correspond very closely with those measured at St. Joseph Island and do not correspond with those located, for example, north of Sault Ste. Marie (see Figure 3.2). Assuming St. Joseph Island is to the Pendills Bay area as Manitoulin Island is to Sault Ste. Marie (after Cowan, 1985), I used Karrow’s lake plane profiles (Figure 5, p116- 117) to correlate features at Pendills Bay with those at St. Joseph Island. In my opinion, correlating two sequences of nearby features with near equivalent elevations makes better sense than extrapolating lake levels from either Sault Ste. Marie or the Penetang Peninsula as Farrand and Drexler (1985) seem to have done. I report updated associations in Table 3.4. In effect, I have recognized connections among relict features stranded in the Pendills Bay area and others stranded at St. Joseph Island and Cockburn Island (Figure 3.5). Having surveyed features in the Pendills Bay area, as well as those at Cockburn Island, I believe the relative strengths of bluffs along each sequence are better matched. Table 3.4: Approximate elevations of glacial Lake Algonquin water planes and the N ipissing water plane measured in western Chippewa County. Approximate elevation (m) Water plane / lake phase 280 Main 273 Ardtrea - Upper Orillia 259 Lower Orillia 240 Wyebridge 232 Penetang 198 N ipissing 78 The Munuscong Islands Taylor (1895) reported observations and provisional height measurements taken along a set of Algonquin strands situated ~1.5 km north of Hessel, MI. According to Taylor, the Munuscong Islands were the only landforms that interrupted the highest Algonquin water surface ”between Mackinac and the Canadian Highlands back of Sault Ste. Marie" (p28). Although the author reports that at least three islands remain, he visited only two: the one closest to Hessel and ”irregular in outline” (p.26); and a smaller island situated northeast of Hessel that exhibited ”precipitous cliffs of limestone 75 to 100 feet in height” (p26). Taylor did not visit the largest of the three islands he mentioned, which is situated northwest of Hessel. Eighty-five years later, Futyma (1980) revisited the Munuscong Islands. Unfortunately, his report offers no new information about these landforms. A careful review of Appendix A (Futyma, 1980) indicates the author focused on just the one island closest to Hessel and relied only on Taylor's (1895) and Leverett and Taylor’s (1915) early findings. Twenty years after Futyma, Schaetzl et al. (2002) targeted these islands, were finally able to access and study a portion of the unvisited northwestern island, and provided new information. In all, Schaetzl et al. (2002) measured 10 positions in the Munuscong Island region. They identified shorelines associated with the Main, Ardtrea and Wyebridge phases along the central island closest to 79 Hessel; the Wyebridge and Payette phases along a minor islet; and the Wyebridge, Payette and lower phases along the large northwestern island. I conducted fieldwork in the Munuscong Islands region, visited old sites, and collected five new position measurements (Figure 3.11). I focused most of my attention on the large northwestern island. The upland covers the better part of 9 km2 and was accessible to me only via a poorly maintained two-track named Fish Road. All other points of entry, particularly the one along Criderman Road, are gated and posted. In the southern half of Sec. 14 T43N R2W, I observed two bluffs, one massive (height > 35 m) and one lesser, but only by comparison. Huge limestone boulders litter the bases of these giant bluffs. The sheer sizes of these well-rounded boulders, perched atop beach sand, clearly imply powerful water action. Unable to receive GPS radio signals at positions along the bluff toe, I traced the features northward into the northern half of Section 14. There I was able to measure positions along the upper and lower shorelines. The upper shoreline is stranded at 266.9 m (SID172) and I associate it with the Main phase. The lower shoreline is stranded at the 259.7 m (SID171) plane and, I believe, is associated with the Ardtrea phase. I measured a third shoreline position in Sec. 26 T43N R2W at the 264.3 m plane (SID169), which I believe is also associated with the Main phase. I attempted to maneuver toward the northwestern side of the large island, but nearly every pathway and trail is gated or posted. I believe the northwestern side holds additional massive bluffs, which are visually apparent from nearby Interstate 75, and I hope to study them someday. 80 (is); ‘ ' 9 ,. 166:223.8m 3592221 ~2m‘.167:224.5m 9166:222.3m Chippewa County 3.2323355" I . m 319.241.4111 . . The 172:266-9m_i'171:259.7m Munuscong Islands _ . o449:227.0m 169:262.0m+ 1741250-0'“ _ . 448:249.5m “(J-2650'” 450:258.4m » . ,_ 451:265.7m°. . ,2 . ‘1 ' Mackinac County i “2241.4“..- 3' 3.: 17:256-2m I l 3042442111 l , o 4 6 8Knometer§ J i I . . 3 4 5Miles Hessel Figure 3.11: Relict bluff sites in the Munuscong Islands region. Sites are labeled with a unique site identification number and an elevation value (e. g., 642: 200 111). An apparent fourth and easternmost island of the Munuscong Island chain remains unstudied. Much of the land area is currently held by mining interests; a large quarry occupies large portions of Sections 3 and 10 of T42N R1E and, perhaps, is the site where Taylor (1895226) observed his ”precipitous cliffs of limestone.” Table 3.5 lists the glacial Lake Algonquin water planes I associated with paleoshoreline elevations in this region. 81 Table 3.5: Approximate elevations of glacial Lake Algonquin water planes measured in the Munuscong Islands region. Approximate elevation (m) Water plane / lake phase 265 Main Algonquin 258 Ardtrea 250 Upper Orillia 242 Lower Orillia 220 Penetang Kincheloe, Michigan The Kinross State Correctional Facility (formerly the Kincheloe Air Force Base) is situated atop a raised terrace; a conspicuous wave-cut bluff marks nearly every side. Schaetzl et a1. (2002) surveyed the toe of a substantial bluff on its eastern flank and correlated its 221 m elevation with Lake Minong (SID389). I surveyed the same feature at three locations along 2 km of its length (SID166, 167 and 168) and obtained an average elevation of 223 m (Figure 3.11). Farrand and Drexler (1985), working northwest of Kincheloe, MI in the Pendills Bay area, identified Lake Minong shorelines at elevations between 213 and 219 m and contend a ”Nadoway Point-Gross Cap sill” (p29) retained Minong waters in the Superior Basin until it was downcut (or breached by inflow from Lake Agassiz), thereby lowering water level even further. In light of Kincheloe’s relative location with respect to the proposed Nadoway Point-Gross Cap sill and situation with respect to St. Joseph Island (~30 km), I believe the stranded bluff at Kincheloe is too high to be associated with the Minong phase of Lake Superior. I 82 believe this feature is better correlated with the Cedar Point water level of glacial Lake Algonquin. Mackinac Island, Michigan Leverett and Taylor (1915) identified an ancient island at Mackinac Island, MI, which is uniquely situated near the confluence of Lakes Michigan and Huron, and between the Munuscong Islands in Michigan’s Upper Peninsula and the Pellston/ Riggsville Islands region of the Lower Peninsula. They describe ”splendid displays” of ”all the beaches from the upper Algonquin down” (Leverett and Taylor, 1915: p424) and assigned the now deprecated Upper Group, Battlefield and Fort Brady labels, in descending order, to three groups of successively lower relict features. Stanley (1936) studied relict shoreline features along the Penetanguishene Peninsula (a.k.a. Penetang Peninsula) of Ontario and compared water levels there with those he and Leverett and Taylor (1915) observed at Mackinac Island. I include an excerpt from his report here (Stanley, 1936:1955): On theoretical grounds, all the Algonquin planes here discussed should be separated by similar intervals wherever the total amount of Algonquin deformation is equal, provided that tilting between Algonquin and Nipissing times proceeded always in the same direction. For instance, on Mackinac Island there should be a shoreline about 90 feet lower than the highest Algonquin (809-812) 83 in order to correlate with the Wyebridge shoreline in the Penetang area (which is 90 feet below the Algonquin in that region). The type Battlefield beach (719) on Mackinac Island fulfills this requirement exactly. Stanley’s provisional correlation of Ontario’s Wyebridge type beach and Mackinac Island’s Battlefield beaches required a leap of faith considering that he could not identify the Wyebridge shoreline near Wasaga Beach (see Stanley, 1936: 1938 and Table 3 on page 1951) and his work at Mackinac Island yielded, in his own words, ”poor results” (p1955). He expressed uncertainty regarding ”just how the lower Algonquin beaches of the Mackinac region fit in with those near Penetang” (p1955). Deane (1950: 74), working later and near Lake Simcoe, Ontario, argued Wyebridge features are ”not strong” and indicate ”only a short pause of Lake Algonquin at this level.” He also reported elevations and relative strengths of features formed along four upper lake levels (Main, Ardtrea, Upper Orillia and Lower Orillia) and characterized the Main and Upper Orillia paleoshorelines as strongest. Uncertainty about the elevations of features at Mackinac Island and their correlations with others prompted Schaetzl et al. (2002) to visit the island and survey relict features there; they obtained five paleoshoreline positions (Figure 3.12). Schaetzl et al. (2002) correlated the highest bluff (242.0 m: SID307) with Main Lake Algonquin, a lower observation (236.3 m: SID302) with the Ardtrea 84 phase, and three lower observations (~228.0 m: SID301, 305 and 306) with the Wyebridge phase. They claim paleoshorelines in the vicinity of Mackinac Island correlate well with those near the center of Deane’s (1950) study area. Lake Huron % LII/((3 Slrails of M a (‘/\' i II a (‘ Nfichigan Figure 3.12: Relict bluff sites at Mackinac Island. Sites are labeled with a unique site identification number and an elevation value (e. g., 642: 200 m). Although two paleoshorelines at Mackinac Island apparently share elevations with others located at Wyebridge, Ontrario, Schaetzl et al. (2002) did not apply Stanley or Deane’s typology consistently to a third. For example, the 85 Nipissing plane at Mackinac Island occurs at ~193 m (634 ft) and seems to correspond with Goldthwait (1910) and Stanley’s (1936) Nipissing shoreline at Wyebridge, Ontario (194 m, 637 ft). Mackinac bluffs stranded at ~228.0 m (750 ft) have elevations that correspond with Stanley’s (1936) Wyebridge bluff at Wyebridge. But Schaetzl’s highest feature (242.0 m: SID307), which they associated with the Main phase, does not correlate with Stanley and Goldthwait’s Main elevation at Wyebridge (256 m, 840 ft) -- a difference of nearly 14 meters. Therefore, either the highest bluff at Mackinac Island was either last cut by a lower water level, or lower features at Mackinac Island are not correlated with those at Wyebridge, Ontario. The highest bluff at Mackinac Island has long been identified as belonging to the Main phase of glacial Lake Algonquin (Stanley, 1945; Schaetzl et al., 2002), so I continue to recognize this association. The next lower bluff rests 5-6 m lower than the Main bluff, which is typical of Ardtrea bluffs located near Orillia, Ontario, so I continue to recognize this association. I am compelled, however, to reject the Wyebridge correlation Schaetzl et al. (2002) made with features stranded at 228 m. Conspicuous examples of Upper Orillia shorelines and few instances of Wyebridge shorelines are present at Sault Ste. Marie, Cockburn Island, the Munuscong Islands, the Penetang Peninsula, or in Lake Simcoe District, so it seems unlikely that Wyebridge shorelines would be expressed better at Mackinac Island than those carved elsewhere, or expressed stronger than those carved by the Upper Orillia water plane. So, in light of the elevations 86 and relative strengths of Upper Orillia shorelines elsewhere, and the relative location of Mackinac Island with respect to these Upper Orillia shoreline features, I correlate features at ~228 m on Mackinac Island with the Upper Orillia phase of glacial Lake Algonquin and not with the Wyebridge phase. Table 3.6: Approximate elevations of glacial Lake Algonquin water planes measured at Mackinac Island, Michigan. Approximate elevation (r11 Water plane / lake phase 242 Main 236 Ardtrea 228 Upper Orillia ? Lower (observed but not surveyed) The Douglas Lake area, Michigan Schaetzl et al. (2002:406) began classifying their shoreline data near Douglas Lake and the University of Michigan Biological Station in Cheboygan County, MI ”because only two highly distinct shorelines exist in this region.” According to their work, and their interpretation of Leverett and Taylor’s (1915) previous work, relict shorelines of the Main phase of glacial Lake Algonquin are expected at the 225 m elevation plane and remnants of a lower and later shoreline are expected at ~210 m. Interestingly, Spurr and Zumberge (1956) also argue that two distinct shorelines exist in this region; they accept the Main phase at 738 ft (225 m) but claim residual segments of a pre-Algonquin shoreline are stranded at 724 ft (220 m). Spurr and Zumberge (1956) make no mention of a post-Main shoreline at 210 m. Obviously, these reports differ in their 87 interpretations. I am disinclined to accept Spurr and Zumberge’s (1956) pre- Algonquin lower shoreline hypothesis, however, because I expect that any feature composed of unconsolidated sediments (as are the ones in question) would have been buried, reshaped or reworked under the high waters of the relatively long-standing Main phase. Also, the idea that drowned shorelines simply re-emerged unscathed after several hundred years of subaqueous exposure, especially in an environment that suffered ”the long fetch afforded by Lake Algonquin” and ”frequent and intense cyclonic storms” (Krist and Schaetzl, 2001:14), seems unreasonable. Before conducting fieldwork, I reviewed Schaetzl's data and plotted them on the Pellston, Indianville, and Mullet Lake MI 7.5 minute USGS topographic map sheets. Of the thirteen points they surveyed in this area, four have elevation values near 230 m, six have values near 225 m, one at 221 m, and two others at ~211 m. Given the 3-meter limit of precision I have adopted during this work, the differences among these elevation values suggest that four, not two, distinctive water levels existed in this region. I revisited field sites at locations around Douglas Lake and surveyed other relict shorelines nearby (Figure 3.13). I started at the boat ramp located due south of Pells Island and traveled counterclockwise around the upland known as Pellston Island. In SEi Sec. 30 T37N R3W, I measured the gravel lag at the base of a protruding relict bluff at 222.2 m (SID180). I then followed the base of the upland as it crosses over the Pellston Island spit (Krist and Schaetzl, 2001) and 88 then south and down into SE34- Sec. 25 T37N R4W; where I measured the base of a relict bluff at 223.0 m (SID181). I continued southward and surveyed several positions along the west side of the island (unfortunately, most of these points were later discarded because each sample contained large anomalies, thus indicating the nearby tree stands or the intermittent snowfalls of the week perturbed GPS signal reception). In NE% SW% Sec. 7 T36N R3W, the base of a stranded bluff rests at 223.0 m (SID182). Schaetzl et al. (2002) measured a 221.6 m elevation at a location 1 km east of this point and along the same feature (SID328). The apparent difference in elevation between these two samples (~1.4 m) is less than the critical threshold value (3.0 m) used to distinguish distinctive lake levels, so they are considered to represent that same water plane. I was able to observe several remnants of what I assume to be lower and later shorelines as I traveled around the north end of Burt Lake. For example, as Indian Road approaches Burt Lake, several houses are situated along a bluff cut at 187.6 m (SID184), which Spurr and Zumberge (1956:102) might associate with the Nipissing Great Lakes. On the east side of Burt Lake, a pipeline runs north- south through Sec. 12 T36N R3W and crosses two apparent shorelines. Admittedly, the bluffs here are weak but they are discernible. The higher bluff was measured at 211.8 m (SID185) and the lower at 197.0 m (SID186). Returning to Douglas Lake, starting at South Fish Tail Bay (NE % Sec. 33 T37N R3W) and working westward, a single relict shoreline and bluff can be traced along and above the southern rim of Douglas Lake. I found an instance of 89 a double bluff near Grapevine Point, which is located due north of the University of Michigan Biological Station campus (Sec. 28 T37N R3W). Figure 2 in Spurr and Zumberge (1956) illustrates the feature. The base of the lower bluff is cut at 220.3 m (SID179), which is 3 m higher than the average water level of Douglas Lake. The base of the upper bluff rests at 227.4 m (SID178). ,333:226.1m 1%.; 469227.4in.468;231.2m . 470:230;6mé471:233.4mr ; ' ' , . , 5388:229.7m V , , - 7: * _ .322:224.6m «- A ‘ ' ”3467:2262m , ”V.- Aiunoo u'ésfioqauo 3261210317? . - ,185321130’1 . yi‘1863197flm _ 0323" 3252253 _“ 472:225.8m "3972172111 , 3411225'9'3398221.s'm -2o7.4m 1.0 ~.2. 21' 6 anaemia l . ' . 37912182111 '0 1 . 2 3 4.45m»; e , . . ____ ,7 _._.__ 7,. ...___._1 Figure 3.13: Relict bluff sites in the Douglas Lake area. Sites are labeled with a unique site identification number and an elevation value (e.g., 642: 200 m). 90 I surveyed other relict bluffs north of Douglas Lake (e.g., NW of North Fishtail Bay in Sec. 22 T37N R3W and SE of Lancaster Lake in Sec. 8 of the same township), but had to discard these samples because they contain highly scattered position data. My field notes indicate a heavy snow / rain mix was falling during data collection, so it is likely that precipitation interfered with GPS signal transmissions as water collected atop the antenna. Given the complex glacial history of this former ice-marginal region (Blewett, 1995; Spurr and Zumberge, 1956), it comes as little surprise that so many different relict features over a relatively wide range of elevations exhibit the erosive influences of water. I believe, however, that there is evidence to support six distinctive water levels (see Table 3.7) in this area. Table 3.7: Approximate elevations of glacial Lake Algonquin water planes and the N ipissing water plane measured in the Douglas Lake area. Approximate elevation (m) Water plane / lake phase 230 Ice-contact margin? 225 Main 221 Ardtrea 212 Upper Orillia - Lower Orillia 197 Uncertain 188 Nipissing As mentioned above, I reviewed Schaetzl’s data for this region and found 4 points with elevation values near 230 m (e. g., SID327), which is 5 meters higher than the generally accepted elevation (~225 m) of the Main phase (Spurr and 91 Zumberge, 1956; Schaetzl et al., 2002). Features at this elevation might mark an ice-contact margin described by Spurr and Zumberge (1956). I removed these points from the glacial Lake Algonquin dataset and placed them into the supplemental dataset. Among the areas covered by the Pellston, Indianville, and Mullet Lake, MI topographic map sheets, 8 samples now mark relict bluffs along the 225 m elevation plane. Given the long history of associating features at this elevation to the Main phase of glacial Lake Algonquin, I am compelled to do the same. Smaller sets of lower shoreline relicts are associated with later and lower glacial Lake Algonquin lake phases as indicated in Table 3.7; the lowest position is associated with the Nipissing Great Lakes. Beaver Island, Michigan Beaver Island is a part of Charlevoix County, Michigan. It is situated among a small group of islands in Lake Michigan, approximately 50 km north- northwest of the City of Charlevoix and 55 km west of Mackinac City. Leverett and Taylor (1915:425) described an upland area on Beaver Island as a ”morainic tract about 4 miles long from northeast to southwest and 1% miles wide” with altitudes ”considerably above the highest Algonquin beach.” The authors explored ”unusually well developed” (p.425) beaches on the north, south and west sides of this upland; they make no mention, however, of any relict coastal landforms on the eastern side of Beaver Island. 92 Dietrich (1978) visited Beaver Island fifty years later, examined soils and landforms around Leverett and Taylor’s (1915:425) ”morainic tract” and found features that ”appear to be beach gravels that are well below what are generally recognized as Lake Algonquin levels” (p.287). He reasoned these beach gravels might belong to later and lower glacial Lake Algonquin lake phases because they occur at elevations higher than the expected Nipissing water plane. Furthermore, he claimed the most substantial bluff features, particularly those situated on the southern and western sides of the relict upland, were eroded by Lake Nipissing and not by Lake Algonquin - as Leverett and Taylor (1915) claimed. Dietrich (1978) continued to argue that ” the island’s only prominent wave-cut scarp attributable to Lake Algonquin” exists on the eastern side of the island ”between Iron Ore Bay and Miller Marsh" (p.287). He also speculated that two small inland lakes (Fox Lake and Green’s Lake) existed as bays in glacial Lake Algonquin. According to Dietrich’s (1978) interpretation of Leverett and Taylor’s (1915zp439) isobase map, the highest glacial Lake Algonquin water plane intersects the contemporary land surface between the 215 and 230 m planes of elevation. The highest elevation on the island surface is approximately 245 m, so it is possible that a proto-Beaver Island existed as an island in glacial Lake Algonquin. The apparent conflicts between claims made about the relationship between the highest Lake Algonquin water plane and the locations of relict lacustrine features on Beaver Island prompted me to visit the Island and examine these features. 93 I conducted fieldwork on Beaver Island, Michigan during late October, 2003 (Figure 3.14). I used the Beaver Island North and Beaver Island South 7.5 minute USGS topographic maps and The Wojan-Cashman Map of Beaver Island (Cashman and Wojan, 1977) as field guides. The provisional USGS maps provide good hypsometric information, but the topology of the local road network is better represented by the Wojan-Cashman Map. The morainic tract is situated in the southern half of the modern island. Michigan BeaVe‘r' Island ‘ ;, 62:229.,3m f" , ' , 0915922232111 _. 1 53160325111qu , '_ .i’163:229.0m. Figure 3.14: Relict bluff sites on Beaver Island, Michigan. Sites are labeled with a unique site identification number and an elevation value (e.g., 642: 200 m). 94 I selected an initial set of twelve potential bluff targets, which includes only those locations with elevations above the expected Nipissing water plane (189 m after Dietrich, 1978). Unfortunately, none of the maps I used illustrates the actual state of land ownership fragmentation or the maze of private in- holdings that exists within the Jordan River State Forest. Many in-holdings are fenced, gated and padlocked, so access to relict features on these properties is prohibited without permission. I was able to contact several landowners and obtain permissions to access their lands, but was unable to secure enough permissions to visit all potential bluff targets. Mr. Larry Dawson owns property (Sec. 29 T38N R10W) that abuts Leverett and Taylor’s (1915:425) ”morainic tract.” A moderate bluff (~6 m high) is cut into the northern side of the tract. No other feature is higher at this location. I sampled soil at the top of the bluff and found red-brown till with coarse gravel. Sampling at the base of the bluff, however, revealed a thin layer of organic material and colluvium (~10 cm) over a thicker layer of beach sand (~1 m). The juxtaposition of soil samples taken at such proximate locations suggests this bluff marks the highest remnant shoreline on Beaver Island. A buried lag of small rounded gravel rests within the layer of beach sand at 229.3 m (SID162), which is approximately 6 In higher than the highest shoreline that Dietrich (1978) measured. I followed the bluff eastwardly and down into a wetland where it loses its distinctive toe. I continued to follow the upland tract around the southern rim of 95 Fox Lake where I found a second wave-cut bluff at 223.2 m (SID159). The elevation of this feature is similar to the highest elevation that Dietrich (1978) reported and is approximately 2 meters higher than the average water level of Fox Lake (as indicated on USGS tOpographic maps). The near coincidence between bluff and lake elevations is, perhaps, the reason why Dietrich (1978) speculated Fox Lake once existed as a bay in glacial Lake Algonquin. Sampling at this location again revealed red-brown till at the bluff top and beach sand at the bluff base. The lower relict shoreline lost its distinctive character and graded into a wetland as I traveled southeasterly along the upland tract. I could not distinguish the higher shoreline feature (229.3 m) anywhere along the southern rim of Fox Lake. About 1.2 km farther east, the lower shoreline regains its character in the form of a very large bluff (~15 m). Sampling at the top of the bluff revealed red- brown till composed mainly of sand with some gravel. Large boulders are scattered across the surface of the bluff top. Soil sampling at the base of the bluff (SID160: 225.1 m) revealed a layer of plant detritus, which rests over a layer of peat, which rests over a layer of beach sand. An apparent relict bar sits 50 m ”lakeward” of the relict bluff, which may have acted as a barrier between open water and a lagoon that allowed the peat to form. Vegetation is dense and the ground is very wet along the base of the large bluff. So, I retreated uphill to a trail that parallels the bluff top and continued eastward. After another kilometer or so I left the trail and followed the morainic 96 tract as it turns south. A few hundred meters farther, I found a conspicuous wave-cut bluff cut into the eastern side of the morainic tract at 229.0 m (SID163). Beyond this site the bluff loses its distinctive toe, vegetation becomes denser and the ground gets wetter. I assume this feature is correlated with the high feature I found on Mr. Dawson’s property. I could not distinguish lower shoreline features at this site. The eastern flank of the morainic tract is very difficult to approach from the northern end. I attempted to follow the Dotty Trail as marked on the Wojan- Cashman Map of Beaver Island (Cashman and Wojan, 1977), but found the path flooded and impassible. The western flank is also very difficult to approach because almost every access road is gated and locked, or posted with a hostile message for potential trespassers. Lacking permission, I did not search for relict shoreline features along these two flanks. The eastern and western flanks of the morainic tract converge toward the southern end of proto—Beaver Island, but the tip of the tract has apparently been truncated by Nipissing waters (Dietrich, 1978). I visited the prominent Lake Algonquin scarp Dietrich (1978) described along the eastern side and ”between Iron Ore Bay and Miller Marsh” (p.287) and measured three elevations along one-kilometer of its length (214.0 [SID165], 212.5 [SID161] and 211.4 m [SID164]). Each value is within 3 m of each other, so all are considered marks along one relict shoreline. In my opinion, this feature serves to distinguish a third and lower Algonquin water plane at Beaver Island. 97 As a final note, the supposed morainic tract changes character halfway along its 8 km length; rough hills typical of the northern 5 kilometers give way to a relatively smooth surface typical of the southern 3 kilometers. As mentioned above, soil samples taken above the highest shoreline in the northern zone yielded red-brown tills with gravels. Sediments sampled along the smooth surface near the southern edge, however, did not reveal a red-brown till, but a sequence of well-sorted sands and gravels. A gravel pit at the southern tip of the upland reveals a profile of sorted sediments (Figure 3.15). Figure 3.15: Exposed sediment profile in a gravel pit located near the southern end of proto-Beaver Island. A 2-foot carpenter’s level is shown for scale. In my opinion, the apparent differences (i.e., topographic expression and general sediment structure) between the northern and southern zones of the central upland on Beaver Island indicate the two zones were shaped by different 98 processes. Also, the elevation of the land surface above the gravel pit (at 220 m) is lower than the two highest shoreline elevations I measured at the north end of the upland (at ~229.0 m), so it seems reasonable to suggest the southern zone of the upland was formed, or at least shaped, under lake waters. I measured relict shoreline features at three distinct levels of elevation on Beaver Island, MI (see Table 3.8). The highest shoreline, which is displayed best at Mr. Dawson’s property, is correlated with the Main phase of glacial Lake Algonquin. Paleoshoreline elevations at Beaver Island correlate well with others stranded nearby on Leverett and Taylor’s (1915) Brutus Island (Schaetzl etal., 2002). Evidence of other water planes exists at Beaver Island, particularly along the eastern and western flanks of the central upland, but the beach gravels Dietrich (1978) described are beyond the scope of this research. I believe much more can be learned from further study at Beaver Island. Table 3.8: Approximate elevations of glacial Lake Algonquin water planes measured at Beaver Island, Michigan. Approximate elevation (m) Water plane / Lake phase 229 Main 223 Ardtrea 212 Wyebridge From Petoskey to Traverse City, Michigan The area between Petoskey and Traverse City, MI contains many upland areas that contain long and narrow lakes (e.g., Lake Charlevoix and Walloon, 99 Intermediate and Torch Lakes). The terrain is complex and relatively difficult to interpret. Leverett and Taylor (19152426) speculated glacial Lake Algonquin ”probably” extended up these valleys and left ”distinct shorelines around Pine Lake” (now called Lake Charlevoix - see Figure 3.16) that ”mark the highest Algonquin beach and some lower members of the upper group.” The authors also mentioned relict features at ”Undine and Hortons [sic] Bay” (Leverett and Taylor, 1915: 426). h-k ”1...... 1 Petés a; ’ ,7- --. E ‘ v.5? 7 ;’11o:g13.om ‘ j; c _ z r # 1695213: . 12.. 3 , 1,; i' ‘ ' ' 3 «.3 ,4 3451 lg; J. Figure 3.16: Relict bluff sites near Petoskey, Michigan. Sites are labeled with a unique site identification number and an elevation value (e. g., 642: 200 m). 100 A high bluff parallels the southern limit of Little Traverse Bay west of Petoskey. I studied three sites near Bay Shore and along Pincherry Road (SID108, 109 and 110). All three have ~213 m elevations. This shoreline is the highest in the region and, according to Spurr and Zumberge (1956: Figure 2), it represents the Main phase. I do not agree with this correlation, however, and attribute it to the Ardtrea phase instead. I considered the well known instance of a Main phase bluff that is situated along the north shore of Little Traverse Bay. Leverett and Taylor (1915:424) describe ”the best-preserved parts” of ”the upper group of Algonquin beaches” between Middle Village and Canby, MI. And, I considered that Schaetzl et a1. (2002) measured bluff elevations at two sites near Good Hart, M1: the highest was at 224.4 m (SID313), which they attributed to the Main phase; and the lower one at 221.8 m (SID314), which they attributed to the Ardtrea phase. These features are well preserved (pers. obs.), where not interrupted by the N ipissing transgression and can be traced confidently for fifteen kilometers to Harbor Springs. A site situated approximately 4 km west of Harbor Springs marks the Main plane at 222.4 m (SID315); another situated 4 km west marks the Ardtrea plane at 215.9 m (SID316). Supposed ”Main” shorelines at Bay Shore and Harbor Springs have elevations that differ by 9.4 m but are separated by no more than 10.5 km, thus implying a rate of uplift near 0.9 m / km. An uplift rate as large as 0.9 m/ km is not expected at this location or latitude and, if proposed, would imply some kind of geologic unconformity exists. Lacking evidence of any such unconformity, I do not carry the Main phase across 101 Little Traverse Bay. Instead, I carry only the Ardtrea phase, which implies a 0.2 m/ km rate uplift. Although such a correlation might seem provocative, I am not the first to propose it. Larsen (1987) summarized ”a lack of consistency between the [Michigan and Huron] basins, indicating that the [Main] Algonquin level of Michigan is considered younger than that of Lake Huron” (p20). Between the Undine Cemetery and Horton Bay (SEi Sec. 1 T33N R7W), I measured a well-defined instance of a Nipissing shoreline; a small bluff is cut at 186.3 m (SID113), situated a few hundred meters north of Lake Charlevoix (not shown in Figure 3.16). Goldthwait (1910: Figure 4) reported a Nipissing elevation of 613 ft (186.8 m) near Burgess, MI, which is located only a few miles away. Many houses along Boyne City Road in Horton Bay are built into a high bluff. For example, at Mr. Jim Chellis’ property (SE% NEi Sec. 6 T33N R6W), I measured the old shoreline at 208.4 m (SID112). This shoreline is the highest apparent shoreline in the region and I correlate it with the Upper Orillia. The Horton Bay area exhibits several intermediate shorelines situated between the high bluff at Mr. Chellis’ property and the current water level of Lake Charlevoix. These features are represented by a series of beach ridges that are stranded along what once was a low-gradient shorezone. These intermediate shoreline features were not surveyed because beach ridges are not always stable indicators of former water level elevations. 102 Leverett and Taylor (1915: 426) claimed glacial Lake Algonquin once covered ”Grass, Clam, Torchlight and Elk lakes to the basin of Grand Traverse Bay” (T orchlight Lake is now known as Torch Lake). They offered, however, only one description of a relict landform: A mile and a half south of Elk Rapids the upper beach is marked by a strong gravel ridge and a hooked spit, and one or two of the members below it are faintly developed. It runs thence south close to the shore of East Arm, much of the way as a wave-cut bench and bluff, in which form it passes south into swampy flats (Leverett and Taylor, 1915:426). These features are clearly illustrated on the Elk Rapids, MI and Williamsburg, MI 7.5 minute USGS topographic map sheets. Plotted contour lines indicate, however, that these features occur at elevations near 185 m, which is lower than expectations of glacial Lake Algonquin water levels but consistent with expectations of Nipissing water levels. Field observations confirm the elevations. Further review of the topographic maps prompted me to inquire about other features in the area - particularly those situated along Torch Lake. For example, just north of French Point (Sec. 18 T29N R8W), I observed what I believe to be a set of two relict shorelines as indicated by a double bluff. The lower shoreline is situated at 185.5 m (SID119), which I correlate with the Nipissing Great Lakes. The higher bluff is situated at 193.6 m (SID117 and 118). 103 I traced the higher bluff both northwardly and southwardly and looked for other conspicuous instances to sample. Most of the land area along Torch Lake has been disturbed by residential and seasonal home development, which rendered several observable instances unsuitable for sampling. I did, however, find one instance situated on a vacant lot that had only been staked for development. Here, a relict shoreline is marked by a substantial bluff at 192.2 m (SID116). I found no evidence of higher or lower Algonquin water levels around Torch Lake and remain uncertain about which phase of glacial Lake Algonquin, if any, these points are associated with. So, I placed them into the supplemental database. The Leelanau Peninsula, Michigan Leverett and Taylor (1915) indicated the Algonquin beach is well—defined along the Leelanau Peninsula of Michigan. They identified a relict islet near Omena, a high lake cliff near Northport and several sets of lower Lake Algonquin beach ridges among several other locations. More recently, Kincare and Larson (2002) mapped the surficial geology of Northern Leelanau Peninsula and provisionally marked the locations of the uppermost Algonquin shoreline. The authors claim ”isostatic uplift has elevated Algonquin shorelines from about 197 meters (646 ft) in Leland to just over 200 meters (656 ft) near Lighthouse Point.” They did not, however, find any evidence of ”lower Algonquin strands.” Instead, Kincare and Larson (2002) suggest that the low strands that do exist likely belong to the Algoma and Nipissing Great Lakes. 104 I investigated sites described by Leverett and Taylor (1915) and mapped by Kincare and Larson (2002). Figure 3.17 show the area. The upper Algonquin water plane and the later Nipissing Water plane are expected to be separated by less than 20 meters of elevation. In some locations, the higher and older shoreline was obviously cut away by transgression of the lower. In other locations, both relict shorelines have been cut away by recent wave action and substantial bluffs now mark the contemporary shoreline. Grand Tram: rse B a y Figure 3.17: Relict bluff sites on the Leelanau Peninsula. Sites are labeled with a unique site identification number and an elevation value (e.g., 642: 200 m). 105 In Sec. 15 T30N R11W, I measured the base of a strong wave-cut bluff at 199.5 m (SID101); it is the highest wave-cut feature in the area. As the ground slopes away from the bluff and toward Suttons Bay it exhibits a sequence of minor ridges and terraces. The ridges run parallel to each other and the modern shoreline. While these ridges and terraces may represent agricultural modifications, they may, rather, represent relict beach ridges associated with falling water levels. Additional study is needed. A small orchard is situated at the site and trees have been planted along these terraces. In SW% Sec. 35 T31N R11W, Freeland Road switches direction to ascend a steep bluff. Kincare and Larson (2002) indicate the Algonquin shoreline occurs here. A gravel lag is buried at the base of the bluff and rests at 199.3 m (SID102). Incidentally, this site is one of the last sites along the Leelanau Peninsula that remains forested, wet, and undeveloped for residential or seasonal home uses. I visited the knoll Leverett and Taylor (1915) described east of Omena. While the authors suggested several lower Algonquin members are exhibited here, my observations indicate only one high and conspicuous shoreline is present; it has been truncated by the N ipissing shoreline in two places. I found no evidence of any higher, lower or intermediate shorelines situated around the knoll. I surveyed the northeast side of the knoll and the uppermost shoreline was cut at 196.9 m (SID103); the Nipissing shoreline is at ~186 m (SID119). Michigan Route 22 approaches a glacial Lake Algonquin bluff about two km north of Omena. The relict shoreline rests at 196.4 m (SID104). 106 Leverett and Taylor (1915) described the highest Algonquin beach as it ”turns directly west across the peninsula about a mile north of Northport" (p.426). I visited this area and believe their description is not correct. Although the Algonquin bluff does begin to turn west north of N orthport, it does not proceed directly across the peninsula. Rather, the shoreline traces the back of a small bay (Kincare and Larson, 2002) and returns east just north of Garthes Pond. From Garthes Pond, I observed the bluff to head in a northeasterly direction for another 2.2 km. It is from here that the shoreline turns west and across the peninsula. A Lake Algonquin shoreline is marked at 198.0 m (SID105) by a massive bluff near Christmas Cove (Sec. 21 T32N R11W). The bluff has been truncated by the contemporary shoreline of Lake Michigan. I did not find lower shorelines below the massive bluff and above the Nipissing water plane. The high truncated bluff reappears south of Christmas Cove and near the Onominese Cemetery (Sec. 5 T31N R11W). I measured the shoreline and it is cut at 201.0 m (SID106). This site is cluttered by gravelly colluvium in many places below the steep bluff face, so this elevation is considered high. As mentioned above, Kincare and Larson (2002) indicate the Algonquin shoreline is situated at ”about 197 meters (646 ft) in Leland.” I observed the base of a wave-cut bluff at 194.2 m (SID107: 637 ft) approximately two km east of Leland and near the eastern shore of Lake Leelanau (SW% SW% Sec. 11 T30N R12W). This feature is well-defined for about one km (near Wardens Point) and sheltered from westerly winds. 107 Prom Mullet Lake to Alpena, Michigan Schaetzl et al. (2002) followed Leverett and Taylor’s (1915) descriptions of relict landforms and surveyed many in the area around Black Lake (Figure 3.18). Most have elevations that correspond with others located in the Douglas Lake area. Each water plane tends to decline as one travels south toward Alpena. Such correspondence is expected given that paleoshorelines are oriented parallel to the generally presumed isobase of isostatic rebound (after Hough, 1958). 'u .. co, . . 9 1." .o C m '92 0 , 9 Lake MN .. g Huron . .q _‘_ ‘2 _ -- I, , 368:229.5m ,5 A . Q‘ 5357:2213m '*: q 1;“ ’; ',465:221.7m j 9:. g - 464:217.2m f‘ rfini35622202m 364:227.8m , . 6362:218.5m - - “sz ,..-... 367229.411: a: _ é - x 4 ”2“: :r‘ ’35952095m j _344:219.0m - 7 _ ”" fl J , a 343:223.0m ; ' '- .- 3. . .360:222."7Ijn - . .j <83. '342:225.5m 36122591111 .. ' _ j < ' f. .- _‘ " ,g .4 * ~' 23% ..... :3 ‘7 Q . r. " I a O 'l' 9' V "H ‘3 {a 2: ..-: r-”.,0 2_ 4.7.6? BKHometers —r 3 ”I" — {2 .l' 4:) _ - s ‘ . o 1 2 V‘s/.74.. 5,Miless Figure 3.18: Relict bluff sites in the Black Lake area. Sites are labeled with a unique site identification number and an elevation value (e.g., 642: 200 m). 108 The highest paleoshoreline (227 - 222 m) is conspicuous and easily associated with Main Algonquin. A few instances of a slightly lower water plane are associated with the Ardtrea phase (222 — 217 m). In keeping with my recognition of Orillia water planes in Ontario and elsewhere in Michigan, I am compelled to update several sites that had been attributed to the lower Wyebridge phase (5113399, 410). .3702206.5m .334:210.4m 335:519’6m' .418:193.9m 7 , 420:201;6m ' ' .,,._ - 421:221.6m‘. Figure 3.19: Relict bluff sites in the Long Lake area. Sites are labeled with a unique site identification number and an elevation value (e.g., 642: 200 m). 109 Schaetzl et al. (2002) also surveyed features in the area around Long Lake (Figure 3.19). The authors attributed a single instance of a very low paleoshoreline (201.6 m: SID420) with the Payette phase (Figure 3.19), but I do not agree with this correlation. Interpretation of Karrow (1987, Figure 5) and Stanley (1936: Figure 2) indicates the Payette water plane is expected to dip below the N ipissing water plane well before it reaches the contemporary Michigan shore; consequently, Payette shorelines are not expected to exist landward of Alpena, MI. So, the instance of a very low paleoshoreline must represent a water plane higher than the Payette water plane. From Alpena to the Au Sable River, Michigan The Thunder Bay River area represents a conspicuous break between sets of paleoshorelines stranded north and south of Alpena (Figure 3.20). Leverett and Taylor (1915) expressed great difficulty regarding their attempts to identify distinctive water planes along the flat and sandy terrain, and to connect these two sets. Also, Schaetzl et al. (2002: Figure 4a, p411) show their statistical model of the Main Algonquin water plane did not perform well in this region (as expressed by relatively large positive residuals). The Thunder Bay River area seems to present the same problem I encountered at Little Traverse Bay: if the nearest highest feature to the south is associated with the nearest feature attributed to the Main phase, then the local linear rate of uplift between the two sites is surprisingly high. Perhaps this is no coincidence. 110 ' 410:2142 m Lake Huron Th 21 11 d e r B a .1/ p 25*; {4 :l‘21‘lzgmys5a1asm ‘ 43020936111 ' Figure 3.20: Relict bluff sites near Alpena, Michigan. Sites are labeled with a unique site identification number and an elevation value (e. g., 642: 200 m). I think it is reasonable to speculate that features situated south of Alpena, having been worked by several converging water planes, might not express the elevation of the first water level to act upon it. Therefore, Main Lake Algonquin features may not be present here. As best as I can determine, only two sites south of Alpena express Main Lake Algonquin (214 m: SID431, 435). The balance, I believe, are associated with the Ardtrea and Upper Orillia water planes. Paleoshorelines situated near Harrisville, between Alpena and the Au 111 Sable River, have elevations that correspond very closely with others located south of Little Traverse Bay (Figure 3.21). 'lz: him/'11 ’ lit/tr I. a It 13 H 1.1 r a II Figure 3.21: Relict bluff sites near Harrisville, Michigan. Sites are labeled with a unique site identification number and an elevation value (e. g., 642: 200 m). From Au Sable River to Saganing, Michigan The area between the Au Sable River and Saganing, MI is generally composed of fine lucustrine materials (Farrand and Bell, 1982) and exhibits a topography that is generally flat or slightly rolling. Two instances of end 112 moraines composed of fine-textured till can be found near Turner and Tawas City (Farrand and Bell, 1982), however, fewer instances of conspicuous wave-cut bluffs can be found within the region. Leverett and Taylor (1915: p.419) described a discontinuous beach ridge, which they characterize as the main Lake Algonquin ridge, that begins near a cemetery in Bay City, ”passes about 1% miles southeast of Omer” and ”runs about a mile back from the lake [i.e., Lake Huron] northward to Alabaster, in Iosco County.” They also described a ”lake cliff and upper ridge” that crosses the Au Sable River ”4 miles west of Oscoda” (p.420). Recent editions of 7.5 minute USGS topographic maps indicate the features mentioned by Leverett and Taylor (1915) and situated between Bay City and Alabaster have an approximate present-day elevation of 184 m, which coincides with the expected N ipissing water plane at these locations (Goldthwait, 1910). Field observations confirm this elevation. Therefore, these relict features were likely last shaped by post-glacial Nipissing waters and not by Lake Algonquin waters. Leverett and Taylor’s (1915zp420) lake cliff, however, is well-defined on the Foote Site Village 7.5 minute USGS topographic map and rests at present-day elevations greater than 190 m, which is higher than the expected N ipissing plane at this location (Goldthwait, 1910). Burgis (1977) describes the same shore-cliff in terms of a graded deltaic sequence that formed in relatively high lake waters (she suggests Lakes Grassmere and Lundy) and was later truncated by glacial Lake 113 Algonquin waters. She also illustrates a lower delta (Burgis, 1977; Figure 44) that formed at the former mouth of the Au Sable River and in glacial Lake Algonquin waters. Burgis (1977) did not provide elevation data for these features. Old logging roads provide access to the truncated delta, a portion of which rests in the Au Sable State Forest, north of the Au Sable River, and west of Oscoda, Michigan (Figure 3.22). I made observations along the base of a very strong bluff in Sections 26, 23, 14 and 11 of T24N R8E. For example, in NW} NE% Sec. 26, the bluff is approximately 30 m high, well-preserved and steep. I... a It 1: H u. r 0 u ..... 1.3: :11- , F 'f a..." _.A_ r 4121196109?" . . i *5 11221399311139 , > . g; , it...” ~ v..,.,.z_,..,r . “gas-”lags! g , 71221936111 v. . V. Figure 3.22: Relict bluff sites near Oscoda, Michigan. Sites are labeled with a unique site identification number and an elevation value (e.g., 642: 200 m). Schaetzl et al. (2002) reported an elevation of 196.0 m (SID412) along the base of this bluff at a location just north of Bissonette Road (SW1? NEl: Sec. 26). I report elevation values of 196.3 m (SID124: SE1: SWi Sec. 26) and 196.8 m (SID125: NWi N132} Sec. 26). As one proceeds northward into Sections 23 and 14, however, the bluff degrades into a series of benches. Two upper benches are 114 backed by faint scarps with basal elevations at 211.4 m (SID996) and 204.3 m (SID995). The lowest portion descends into a swamp and no conspicuous wave- cut features can be identified along its slope. The upper bench positions, which are not associated with conspicuous wave-cut bluffs, are retained in the supplemental database. South of the Au Sable River in Sec. 10 of T23N R8E, I made observations along the bases of a weak scarp (SID999: 203.0 m - not shown), a strong scarp (SID122: 196.5 m) and a conspicuous wave-cut bluff (SID123: 193.4 m) roughly 6 m high. The two later observations were made along the same feature but at locations separated by 530 m. At this time, I cannot explain the apparent large difference in elevation (~3.0 m) between two instances of the same feature. Additional sampling and, perhaps, a rigorous examination of soil stratigraphy at these locations may yield some explanation. Yet, because the later two observations are associated with a conspicuous bluff feature, their positions are included in the shoreline dataset - as is the uncertainty that exists between them. The highest of the three positions, which is not associated with a conspicuous bluff feature, is retained in the supplemental database. Given the substantial relief of this bluff (~30 m) and the coincidental higher elevations of scarps found at locations both north (204.3 m) and south (203.2 m) of the Au Sable River, the base of this bluff (~196 m) is likely not correlated with the highest Algonquin water level, but a subsequent and lower one. I associate the base of this bluff with the Upper Orillia water plane. 115 Summary I described the methods I used to collect relict shoreline data and substantially improve the glacial Lake Algonquin dataset created by Schaetzl et al. (2002). I surveyed 88 new relict features in Northern Michigan and Western Ontario and was able to attribute 69 features to known glacial Lake Algonquin lake phases; the balance is reserved in the supplemental database. I retained 128 of the 146 positions Schaetzl et al. (2002) collected. Four elevations taken near Douglas Lake, Michigan, mark anomalous features that rest higher than the Main water plane (~225 m). The other fourteen points not retained were taken from features that resemble offshore bars, barrier bars, or spits and, so, were not collected according to my sampling protocol; I reserve these eighteen points in the supplemental database and they await future review. Also, many of the elevation values I obtained at sites in the Munuscong Islands region, near Douglas Lake and near the Au Sable River compliment the elevation values reported previously by Schaetzl et al. (2002) - so it seems our efforts identified the same features and obtained similar elevation values from them. I conflated well-known paleoshoreline data taken by others at Sault Ste. Marie (10) and St. Joseph Island (19). Together, the glacial Lake Algonquin dataset and the supplemental database hold 243 bluff records; Appendix A presents a complete list of them; Figure 6.10 shows bluff locations with respect to predicted shorelines. Individually, each sample marks a point along a former shoreline of glacial Lake Algonquin (although a few mark N ipissing shorelines). 116 In aggregate, they outline several distinct water planes within the region of interest: Main (57), Ardtrea (41), Upper Orillia (33), Lower Orillia (16), Wyebridge (12), Penetang (13), Cedar Point (17), Payette (7), Sheguiandah (3), Korah (1) and others (43). Apparent differences among water plane elevations are useful for modeling the cumulative effects of differential isostatic rebound (Chapter 5). A major difference exists between the method Schaetzl et al. (2002) used to correlate paleoshorelines with named water planes, and the method I used. Schaetzl et al. (2002) began correlating shoreline data and Lake Algonquin lake phases in the Douglas Lake area, worked their way east and west, and then proceeded north and south. I attempted to follow that procedure, but encountered difficulties when attempting to recognize their set of four lake levels at locations north of the Straits of Mackinac - where additional water planes are clearly evident. Ultimately, I started again but began correlating paleoshoreline data with Lake Algonquin lake phases at Sault Ste. Marie instead. Sly and Lewis (1972), Cowan (1985), and Karrow (1987) correlated features in the region with respect to the Lake Algonquin lake phases named by Deane (1950: The Upper Group containing the Ardtrea, Upper Orillia and Lower Orillia phases), Stanley (1936: The Lower Group containing the Wyebridge, Penetang, Cedar Point and Payette phases) and Hough (1958: Sheguiandah and Korah), so it was relatively easy to interpret my results at Cockburn Island in terms of their findings. Having recognized the connections Cowan (1985) made between Manitoulin 117 Island and Sault Ste. Marie, and the apparent connections that I recognize between Cockburn Island and St. Joseph Island, it became clear that relict features in Michigan’s Upper Peninsula needed to be reclassified. From these northern locations I then progressed southwestwardly and continued to recognize water planes as I went. Relict shoreline features at Cockburn Island are particularly important features because they add much needed support near the northern rims of glacial Lake Algonquin water planes. They also provide clues about former ice margin position (Karrow, 1987). We know that the Superior Lobe of the Laurentide ice sheet retreated northeastwardly after reaching a maximum extent that covered the northern portion of Michigan’s Lower Peninsula (Karrow, 1987; Larson and Schaetzl, 2001). Lakes in the Michigan and Huron basins were joined when retreating ice uncovered the Indian River lowlands and, slightly later, the Straits of Mackinac, which Larsen places at approximately 11,200 yr BP (after Hansel et al., 1985). Water levels in glacial Lake Chicago dropped to the glacial Lake Algonquin level upon joining (Hansel et al., 1985), which left the outlet at Chicago abandoned (Larson and Schaetzl, 2001) and directed all outflow to the Fenelon Falls outlets near Kirkfield, Ontario (Finamore, 1985). Thus began the Main phase of glacial Lake Algonquin. The presence of Main phase shorelines on Cockburn Island suggests that the ice sheet retreated through Cockburn Island sometime during the early portion of this lake phase; early enough to allow bluff cutting and other shoreline development processes to occur. Karrow 118 (1987) concluded that St. Joseph Island was deglaciated ”about 11,000 yr BP or shortly thereafter” (p119). Cockburn Island is apparently situated along the same isobase that runs through St. Joseph Island (Figure 3.8), so both islands were likely deglaciated at the same time. Therefore, locations between the Straits of Mackinac and the Cockburn Island-St Joseph Island region must have been deglaciated within, roughly, a ZOO-year period.2 Main phase and lower shorelines on Cockburn Island tell us that Cockburn Island was deglaciated before the Fossmill outlets at North Bay, Ontario, because lower lake phases are associated with falling water levels brought about by deglaciation of the outlets at North Bay (Kaszycki, 1985; Karrow, 1987). Karrow et a1. (1975) reported radiocarbon dates that indicate glacial Lake Algonquin, from the Main phase to the lowest, existed between 11,200 yr BP and 10,400 yr BP. Relict shoreline features at Cockburn Island are important features because they help to place the ice sheet margin as well as the northern extent Main Lake Algonquin into a spatial context within this very short time frame. 2 Karrow (1987) mapped a separation between areas of the Laurentide ice sheet that are known as Superior ice and Algoma ice (Figure 2, page 114). Superior ice covered the Straits of Mackinac and Algoma ice covered St. Joseph Island and Cockburn Island. This separation complicates the deglaciation history of the northern Lake Huron basin and precludes, for now, speculation about how glacial ice retreated in locations between the Straits of Mackinac and the St. Joseph-Cockburn Islands region. 119 Understanding the relict features at Mackinac Island remains a key to establishing the proper link between the Upper Peninsula and the Lower Peninsula- interpretation of these features, however, remains a crux. Proto- Mackinac Island was small and distant from other Algonquin islands, and it holds relatively few conspicuous paleoshorelines to study (despite their notoriety). Schaetzl et al. (2002) correlated features at 242 m with the Main phase, those at 236 m with the Ardtrea phase, and those at 228 m with the Wyebridge phase; they claimed features at Mackinac Island have elevations similar to others in the Lake Simcoe District of Ontario (Deane, 1950). I was compelled, however, to reinterpret these data as I attempted to integrate new data with old findings. I reviewed Deane’s (1950) work and consider the work of Cowan (1985) and Karrow (1987), as well as my own, and found reasons to believe that features Schaetzl et al. (2002) correlated with the Wyebridge phase rest at elevations better correlated with the earlier Upper Orillia water plane. For example, I considered Deane’s (1950) descriptions of features; he characterized Upper Orillia bluffs as strong and characterized Wyebridge features as ”not strong” and indicative of ”only a short pause of Lake Algonquin at this level” (Deane, 1950: 174). I considered my work at Cockburn Island — the strongest bluffs there occur at elevations that correlate very well with the Upper Orillia water plane and rest at elevations higher than the expected Wyebridge water plane. Also, I did not find features on Cockburn Island that were associated with the Wyebridge water plane. Finally, I considered the relative location of 120 Mackinac Island with respect to other places that mark the Upper Orillia water plane. Ultimately, I correlated features at ~228 m on Mackinac Island with the Upper Orillia phase of glacial Lake Algonquin and not with the Wyebridge phase. More study is needed here. Making correlations between relict shoreline features and Lake Algonquin lake phases is difficult in areas near Little Traverse Bay (Lake Michigan) and Thunder Bay (Lake Huron) because paleoshorelines in these areas are oriented roughly north-south, and interrupted by contemporary water features that are oriented roughly east-west. Ever since Leverett and Taylor (1915), the highest paleoshorelines on either side of these water features have been considered as representatives of the same water level - Main Lake Algonquin. Such connections, however, implicitly require anomalous crustal deformations to have occurred in both bays. In Emmet County, for example, the highest paleoshoreline cut into the north side of Levering Island occurs at 232 m and the highest paleoshoreline cut into the south side of Brutus Island (on the north side of Little Traverse Bay) is 222 m; both shorelines are associated with the Main phase (Leverett and Taylor, 1915; Schaetzl et al., 2002). These two sites are separated by approximately 30 km, thus yielding a reasonable 0.33 m/ km linear rate of elevation change along the Main water plane. Continuing southwardly, if we assume the highest paleoshoreline on the north side of Little Traverse Bay is the same as the highest paleoshoreline on the south side of Little Traverse Bay, then that would require an effective 0.90 m / km linear rate of change along the 121 water plane over a mere 10.5 km span. In other words, assuming the highest paleoshoreline on both sides of Little Traverse Bay are the same shorelines would require an anomaly in the isostatic rebound process at this location. A similar anomaly is required if the highest paleoshorelines north and south of the Thunder Bay River are also considered as the same. Lacking any knowledge of anomalous deformations in these two areas, I must consider that the highest paleoshorelines situated north of Little Traverse Bay / Alpena are not correlated with the highest paleoshorelines south of Little Traverse Bay / Alpena. In related work, Goldthwait (1910) described features along the highest paleoshoreline in the area between Kincardine and Port Elgin, Ontario. He also reported several features occur at anomalously low elevations. Karrow (1988) surveyed this region and found that features previously identified as the highest Lake Algonquin shoreline had actually been eroded and lowered by later water levels. Consequently, the features in question indicate the action of younger lake phases and not the Main phase of glacial Lake Algonquin. I speculate that a similar pattern of erosion occurred in the areas south of Little Traverse Bay and Thunder Bay, Michigan. Both Deane (1950) and Schaetzl et al. (2002) reported Main and Ardtrea water planes seem to converge at roughly this latitude, so it is quite reasonable to speculate that later, lower, and converging water planes continued to erode bluffs that might have once held marks of the Main phase, but no longer do. The Main Lake Algonquin water plane elevation at Port Elgin, Ontario (220 m after Karrow, 1988: p159) correlates well with highest 122 paleoshoreline found north of Alpena and west of Long Lake (221.6 m). The elevations of lower relict features situated south of Port Elgin, Ontario (208 m after Karrow, 1988: p159) correlate well with features situated south of Alpena, Michigan (209 m). In both cases, however, the elevations of lower relict features correlate best with the younger Upper Orillia lake phase and not the higher and older Main or Ardtrea phases. Karrow (1987) and Larsen (1987) have noted the long-standing controversy regarding whether certain water planes in the Lake Algonquin series are parallel, thus implicating changes in outlets were responsible for changes in water levels, or whether they converge toward a single outlet, which would implicate isostatic rebound as the sole cause of shoreline deformation. Stanley (1936) and Hough (1958) suggested the shorelines were parallel; each represented a new outlet. Deane (1950) and Schaetzl et al. (2002), however, found evidence of convergence between the Main and Ardtrea water planes; each shoreline represented differential uplift relative to the same (also uplifting) outlet. Larsen (1987) argued for convergence after considering historic uplift data (Clark and Persoage, 1970). This work finds some evidence of both. Like Dean (1950) and Schaetzl et al. (2002), this study finds evidence of convergence between the Main and Ardtrea water planes. Proceeding roughly southwardly: Main and Ardtrea bluffs at Cockburn Island are vertically separated by approximately 8 m (280 and 272 m, respectively); in the Munuscong Islands region they are separated by ~7 m (265 and 258 m, respectively); on 123 Mackinac Island they are separated by ~6 m (242 and 236 m, respectively); and, in the Douglas Lake area, they are separated by ~4 m (225 and 221 m, respectively). Moreover, the Main and Ardtrea shorelines are difficult to differentiate at locations south of Little Traverse Bay (Lake Michigan) and Thunder Bay (Lake Huron) because they either occur within 3 m of each other or they were reworked or destroyed by younger and lower lake phases. Convergence between the Main and Ardtrea shorelines suggests differential uplift relative to the Fenelon Falls outlet was the cause of deformation and tilt during the Main and Ardtrea phases. Upper Orillia shorelines consistently occur below Ardtrea shorelines at Cockburn Island, the Munuscong Islands, Mackinac Island, and in the Douglas Lake area, but elevation differences between them provide conflicting evidence regarding convergence. Between some locations they appear to converge toward the south, while among other locations there is evidence of parallelism. Few sites exist that mark multiple instances of later and lower water planes. Cockburn Island and the Munuscong Islands region, for example, are the only two locations that hold instances of both the Upper Orillia and Penetang shorelines. The vertical difference between these shorelines at Cockburn Island is 36 m, whereas the difference in the Munuscong Islands region is only 30m. A decreasing separating distance with decreasing latitude suggests convergence again, hence the influence of differential rebound. A comparison between two points, however, is not sufficient for characterizing the spatial patterns of two 124 regional water planes, so no valid conclusion can yet be drawn about the relationship between the Upper Orillia and Penetang water planes. Evidence of convergence or parallelism between particular water planes in northern parts of the Huron basin notwithstanding, these data cannot be used to resolve the question of where Lake Algonquin stood in the southern parts of the Huron and Michigan basins. Like Schaetzl et al. (2002), my study area did not extend into those areas. The new elements of the glacial Lake Algonquin dataset, including those surveyed during this work (n=69) or conflated from previous work (n=29), fill large data gaps (Cowan, 1985) and furnish new points of support for the restoration of glacial Lake Algonquin water planes (Taylor, 1895). 125 CHAPTER 4: ACCURACY ASSESSMENTS Introduction All spatial data are, to some extent, contaminated by error (Foody, 2003; Burrough and McDonnell, 1998; Heuvelink and Burrough, 1993; Openshaw, 1989 and Goodchild and Gopal, 1989). Error implies a discrepancy between an attribute of an entity or phenomenon and its measured representation. A potential danger exists in work that uses spatial data because errors can propagate through analyses, render results uncertain and prompt misleading interpretations (Shortridge, 2001; Heuvelink, 1988). Heuvelink and Burrough (2002: 111) have argued the importance of knowing ”how accurate the data contained within spatial databases really are” so the ”value of the derived information” might be known. Accordingly, the accuracies of the datasets employed in this work need to be assessed so that knowledge about data errors can be used to improve data analysis and qualify the results. I report exploratory data analyses and accuracy assessments of the paleoshoreline dataset and a selected portion of the National Elevation Dataset (NED: Gesch et al., 2002). I use the public-domain R environment for statistical computing and graphics (R Development Core Team, 2005), particularly the geoR (Ribeiro Jr. and Diggle, 2001), gstat (Pebesma 2004), maptools (Bivand, 2005), and RandomFields (Schlather, 2001) packages, to explore and analyze data and conduct statistical tests. 126 Two datasets I employ two datasets: a) a primary set of paleoshoreline positions; and b) a secondary digital elevation model that covers the region of interest. Every elevation value in each dataset, however, is the result of a random measurement process and likely contains error. From the outset I expected errors in elevation to have much greater potentials for contributing uncertainty to spatial model results than errors in horizontal locations because a significant difference in scale exists between the horizontal and vertical dimensions that bound glacial Lake Algonquin shorelines. The region of interest (see Figure 1) includes the confluence of Lakes Michigan and Huron and it extends northward past Sault Ste. Marie, Ontario, southward to the Au Sable River, westward to the Leelanau Peninsula and eastward to Harrisville, MI. The north-south dimension of the study area covers approximately 240 km of terrain whereas the east-west dimension covers 195 km. If a paleoshoreline position was to contain a horizontal location error of 10 m, then the error magnitude would represent no more than 0.005% of the east-west extent or 0.004% of the north-south extent. In other words, the relative size of a 10-meter horizontal error would be very small, and if input to a spatial model, then its influence on model results would also be small. I chose the GPS leveling methods and equipment described in Chapter 3 because they typically yield horizontal coordinates with sub-meter errors. A sub-meter error is an order of magnitude less than a 10-meter error, so its relative influence is expected to be negligible. 127 The contemporary range of elevations associated with stranded glacial Lake Algonquin features (~133 m) is much smaller than either range of horizontal coordinate values, which is the important difference in scale mentioned above.1 Consider a paleoshoreline position that contains a 10 m elevation error, rather than a horizontal error. A 10-meter error in elevation would represent 7.5% of the elevation range. Hypothetically, if a 10-meter elevation error were entered into a spatial model, then it would likely perturb results. By comparison, the relative impact of a 10-meter error in elevation is much greater than the impact of a 10-meter horizontal error given the substantial difference between vertical and horizontal extents. Moreover, the contemporary ranges of elevation associated with individual glacial Lake Algonquin water planes are less than 133 m, so individual water-plane models supported by such erroneous positions are subject to even greater impacts and uncertainties. 2 That is why I used differential GPS technology and designed a strict sampling protocol that limited data collection activities to occur only during temporal 1 The average level of water in the contemporary Michigan-Huron basin is 176 m, below which glacial Lake Algonquin shorelines cannot be observed. The highest elevation associated with a glacial Lake Algonquin feature is approximately 309 m (Cowan, 1985). 2 Schaetzl et al. (2002: Table 4, p412) reported three post-Main Lake Algonquin water plane elevations at Alpena, MI and Sault Ste. Marie, Ontario. The differentially warped Ardtrea, Wyebridge and Payette water planes include elevations that range 75, 77 and 51 meters, respectively, between the two sites. 128 windows that framed optimal satellite constellation geometries suitable for the taking of accurate and precise elevation measurements. Differential GPS technology can yield elevations with sub-meter accuracy when used properly. Yet, proper planning and use notwithstanding, it is common knowledge among the GPS community that GPS technology resolves elevation values less precisely than horizontal coordinate values. So, the set of GPS-measured paleoshoreline data, particularly the set of elevation values, might contain errors that need mitigation. I focus my attention on these. And, it is common knowledge that ”elevations recorded within digital models are fraught with errors” (Fisher, 1998:215), so I also focus on them. Assessment of the glacial Lake Algonquin dataset The task of assessing accuracy is concerned with describing errors in data (Shortridge, 2001). Typical assessments begin with the assumption that each error is the result of a random measurement process, proceed by comparing a set of measured values to a set of known or higher-order values assumed to be accurate, and end with a description of the distribution of error. A higher- quality paleoshoreline dataset does not exist to assess the set of data collected during this work. Yet, rather than assume paleoshoreline data are accurate by virtue of being the highest-order available, I attempt an assessment by recognizing potential sources of error and conceptually or empirically estimating their potential values. Each shoreline sample is subject to two sources of error: a) 129 the limits of precision associated with the method I used to select relict shorelines in the field; and b) the limits of precision associated with the methods I employed to measure and process shoreline positions at selected sites. I address both. Analysis of potential selection errors As explained in Chapter 2, I rely on the gravel lags that form at the bases of wave-cut bluffs to be the set of best indicators of former water-level positions. Morton and Speed (1998) suggest the precision with which a lake level can be established by the elevation of a wave-cut bluff alone is sub-meter, but my knowledge of contemporary water-level behavior suggests the precision might be as great as three meters (see Chapter 2). Such imprecision is supported by Goldthwait (1910: Figure 2) who, when drafting profiles of uplifted water planes, used bands ”6 feet broad according to the vertical scale” to represent ”small scale discordances as must be expected because of the original variation” (p37) in water plane heights. Every field site, save the few noted otherwise, exhibited a conspicuous gravel lag from which the limits of long-term wave-action can be inferred fairly confidently. A possibility exists, however, for some gravel lags to have been deposited at elevations higher than long-term water levels (e.g., during storm events). If I selected such a feature, then a probability exists for the selected elevation to overestimate the former water level. I expect such errors are rare. 130 Based on my experience in the field and following the research of others (e. g., Goldthwait, 1910; Morton and Speed, 1998), I expect most selection errors to be 0 m, the majority of selection error values to fall within the interval of 0 and 1 m, and, on rare occasions, extreme values to fall within the interval of 1 and 3 m. Ground control Measurement errors in the paleoshoreline dataset cannot be assessed directly because no paleoshoreline surveyed in this work has been surveyed independently or by higher-order methods; an indirect assessment method is needed. So, I assess the precision of the data collection and processing methods, from which I infer the accuracy of data they yield. Small sets of National Geodetic Survey (NGS) bench marks were measured during each field campaign. The protocol used to measure and process paleoshoreline position data was used to measure and process bench mark data. After each campaign, each set of bench marks was added to the pool of ground control data. Calculated discrepancies between measured and actual bench mark positions are used to assess the accuracy of paleoshoreline position data and correct any bias. The National Geodetic Survey is the authority for determining latitude, longitude, elevation, crustal velocity and gravity throughout the United States and it serves as the depository for geodetic data of the highest-order (NGS, 1998). Geodetic data represent a network of precisely surveyed control stations with accurate position attributes (Bossler, 1984). Horizontal control stations, also 131 known as monuments, are established to provide users with accurate longitude and latitude values. Vertical control stations, also known as bench marks, are established to provide accurate height (elevation) information. Interestingly, not all bench marks have horizontal data of geodetic-quality and not all monuments have height data of geodetic-quality, but some do have both. The NGS maintains a geodetic datasheet (NGS, 2004) for each monument and bench mark. Datasheets are archived annually on a regional basis and those that represent control stations in Michigan were last archived in November, 2004. Datasheets contained within this archive are used to assess the precision of the data collection methods described in Chapter 3 and to assess, by inference, the accuracy of paleoshoreline position data. Twenty-four different NGS control stations (see Table 4.1) were occupied and measured during my field campaigns, of which ten stations have geodetic longitude and latitude information. During data processing, however, I discovered six station samples to be unfit for use. Three samples (PIDs OJ0706, PK0154 and PK0165) were collected on November 13, 2003 and each contains conspicuous large errors - particularly station PK0154. Such large errors prompted me to review all field notes, field equipment and data processing methods, but I found nothing from these sources to help me explain their occurrences. I discovered, ultimately, the US Coast Guard archive of GPS Status Message Reports (US Coast Guard Navigation Center, 2004), which contains historic status information about the GPS space and control segments. 132 Table 4.1: Vertical differences between measured and published heights. NGS bench marks are indicated by their Permanent Identifiers (PIDs). Observations used to assess shoreline elevation data. GPS height NGS height Difference PID (m) (m) (m) N F0260 259.86 260.26 -0.40 PK0227 193.85 194.28 -0.43 PL0361 192.20 192.21 -0.01 PL0369 190.63 190.61 0.02 PL0370 195.62 195.83 -0.21 PL0389 182.93 182.36 0.57 QJ0086 180.31 179.70 0.61 QJ0151 198.21 196.51 1.70 QK0087 185.55 184.84 0.71 QK0325 255.22 255.55 -0.33 QK0410 177.29 177.56 -0.27 QK0469 194.63 194.52 0.11 QK0470 205.24 206.60 -1.36 QK0484 232.83 232.04 0.79 RJ0161 193.20 192.37 0.83 RJ023O 236.24 234.77 1.47 RJ0644 208.01 208.49 -0.48 RJ0691 234.17 234.68 -0.51 Observations perverted by satellite malfunction (rejected). OJ0706 183.32 184.72 -1.40 PK0165 203.88 206.02 -2.14 PK0154 185.28 190.66 -5.38 Observations taken at bench marks with non-geodetic quality height information (rejected). QK0852 243.16 242.00 1.16 * 267.97 266.70 1.27 * 249.40 248.41 0.99 Several GPS status reports (USCG, 2004) revealed that a GPS satellite malfunctioned on November 13, 2003 and was broadcasting bad data for at least twenty-four hours. The GPS signal was eventually turned off, but not before I 133 collected data. I discarded these samples because I could not remove the influences of bad data obtained from the now decommissioned satellite. This event prompted me to review every archived status report for each day data collection activities occurred, which included the previous work accomplished by Schaetzl et al. (2002). Fortunately, I found no evidence to suggest other bench marks or any paleoshoreline data should be discarded for similar reasons. I rejected three more control station samples because each was obtained at a site having non-geodetic quality information. For example, the datasheet that represents bench mark QK0852 contains elevation data not obtained by rigorous survey, but interpolated from an old topographic map built using the obsolete North American Vertical Datum of 1929 (N AVD29). Unfortunately, this mark is the only documented vertical control station on Beaver Island, MI. I rejected two marks collected by Schaetzl et al., 2002 (marked in Table 4.1 with an asterisk - *) because they were taken at road intersections annotated with spot heights on topographic maps (and rounded to the nearest foot), and not at bench marks with high-order geodetic information. The lesser-qualities of these three control points are insufficient for assessing the greater precision of GPS-based measurement methods or the accuracy of the data they yield, so they are rejected. Analysis of measurement errors Differences between measured and published elevations are listed in Table 4.1 for eighteen bench marks having geodetic-quality height information. 134 Summary statistics (Table 4.2) and a stem-and-leaf plot (Figure 4.1a) show measurement errors are small in magnitude; the distribution appears symmetric and centered near zero. Table 4.2: Summary statistics that describe the set of differences between measured and published bench mark elevations. Sample statistics Value Test for mean at 0 Value n 18 minimum -1.36 m median 0.00 m mean (bias) 0.16 m matched-pairs t p = 0.4010 standard error of mean 0.19 m standard deviation 0.77 m maximum 1.70 m 2 -3 i, 1 -2 E 3 -1 4 c. 0 -0 55443320 ‘2; 0 0166788 E 1 57 in" ’1 2 3 -2 -2 -1 0 1 2 Quantiles of a Normal Distribution (6) (b) Figure 4.1: Distribution of error in sampled bench mark elevations: a) Stern and leaf plot - solid line denotes the decimal place; and b) q-q plot of standardized errors with a one-to—one line and 90% confidence envelope. 135 Results of a matched-pairs t test indicate the mean difference between measured and published elevations (0.16 m) is not statistically different than zero at the 90% confidence level. A q-q plot of standardized error values (Figure 4.1b) and results of a Shapiro-Wilk test for normality (a = 0.10: W=0.9540 and p = 0.4915) provide no significant evidence to suggest the sample distribution of elevation error does not fit a Normal distribution. No elevation difference exceeds 1.7 m, which is less than the three-meter threshold used to distinguish different long-term water levels. As expected, the majority of measured elevations (15 of 18, or 83%) contain sub-meter errors. I used scatterplots to examine associations between error and position (Figures 4.2 and 4.3). I used Kendall’s non-parametric tau statistic to test the strengths of associations because latitude and longitude values are non-linear measures of distance. Kendall's tau does not rely on assumptions about the underlying distribution of data. It is applied to the ranks of data (which mitigates the non-linear distance problem) and takes values ranging between -1 and +1 (which can be interpreted like Pearson’s or Spearman’s statistics). No plot or statistical test result provides strong evidence to suggest errors, either signed or unsigned, are associated with location or elevation. 136 2 A 1 .. g 515' B 0 4‘. , t: .. «r a LU -1 . -2 -3 180 200 220 240 260 Actual elevation (m) (a) 3 3 2 s 2 , a a " as g 33 I' '3' ‘P \E’ {5 {E are E 0 . £1 “#4' “Iii-I. 3 E 0 $5 ilk—hi IIIII I? LL! _1 LU -1 T g: {3: -2 -2 -3 -3 -86.0 -85.0 -84.0 -83.0 44.0 45.0 46.0 47.0 Longitude (dd) Latitude (dd) 0)) (C) Figure 4.2: Scatterplots of signed error in sampled bench mark elevations and position. Errors and elevations are reported in meters. Latitude and Longitude values are reported in decimal degrees. 137 3.0 :5: 2.5 {‘3’ 2.0 g, 1.5 .4 as (D E 1.0 -. e 1:: «ft? m 0.5 I” '13 q” 4, 0.0 i 180 200 220 240 260 Actual elevation (m) (a) 3.0. 3.0 E 2.5-: E 2.5- % 2.0 8 2.0 3 . :5 . (U i (U ' ‘ E 1.0 ,. E 1.0 9 - z... 33 1": .3; e .. a: g :1: a; 0.5: ; t...” a. Lil-J 0.5 .. f s 0.0: 1+ 0.0 is i“ -86.0 -85.0 -s4.0 -83.0 44.0 45.0 46.0 47.0 Longitude (dd) Latitude (dd) 0?) (C) Figure 4.3: Scatterplots of unsigned error in sampled bench mark elevations and position. Errors and elevations are reported in meters. Latitude and Longitude values are reported in decimal degrees. I used Moran’s I statistic (Bailey and Gatrell, 1995: 270) to test for spatial autocorrelation. The calculated statistics for signed and unsigned errors (I = 0.047 and -0.090, respectively) are near zero and indicate no spatial 138 autocorrelation is present, so I infer no spatial structure exists among the errors in these bench mark elevations. Table 4.3: Results of Kendall’s non-parametric tests for associations between position and elevation error. Signed error Unsigned error tau p value tau p value Longitude 0.150 0.410 0.294 0.096 Latitude 0.072 0.709 0.268 0.131 Elevation -0.163 0.369 0.137 0.454 The actual distribution of error is not known, but sampled errors appear small, random, normally-distributed and not a function of location, elevation or autocorrelation. Because the methods used to measure paleoshoreline elevations are the same as those used to measure bench mark elevations, errors in paleoshoreline elevations are inferred to be similarly small, random, normally distributed and spatially unstructured. Calculated differences between measured and published latitude and longitude values (converted to meters) are listed in Table 4.4 for seven monuments having geodetic-quality latitude and longitude information (recall that ten were measured but three were rejected because of a satellite malfunction). All location coordinates contain sub-meter errors. On average, a measured location is typically within 0.5 m of the actual location. 139 Table 4.4: Horizontal differences between measured and published N GS benchmark locations. Error in Error in PID Latitude (m) Longitude (m) Distance (m) NF0260 0.16 0.09 0.19 PL0361 -0.47 0.22 0.52 QK0325 -0.96 0.19 0.98 QK0469 -0.48 -0.41 0.63 QK0470 -0.65 -0.02 0.65 QK0484 -0.30 0.30 0.42 R]0161 -0.03 0.21 0.21 The methods used to measure paleoshoreline locations are the same as those used to measure bench mark locations, so horizontal errors in paleoshoreline locations are inferred to be similarly small. Moreover, the maximum error observed in horizontal coordinate data (0.98 m at mark QK0325) is smaller than the widths of most gravel lags that form along contemporary Great Lakes shorelines (pers. obs.). From this observation, and given the large extent of the study area, random and small horizontal errors likely associated with measured shorelines, while likely not zero, are hereafter considered to be not meaningfully different than zero. In short, they are ignored. Cumulative elevation errors in paleoshoreline data Recall that each shoreline measurement is subject to two sources of error: a) the limits of precision associated with the method I used to select relict shorelines in the field; and b) the limits of precision associated with the methods I employed to measure and process shoreline positions. I have no reason to believe that either one of these sources of error can be influenced by the presence, location, magnitude or direction of the other, so I assume that each error process operates independently of the other. Therefore, the total error in any paleoshoreline elevation is expected to be the simple sum of selection and random measurement errors. I describe in Chapter 5 the process I used to estimate total elevation error in my paleoshoreline elevation data. 141 Assessment of the National Elevation Dataset The National Elevation Dataset (Gesch et al., 2002) is a digital elevation model (DEM) that provides seamless data coverage of the United States. It is produced by the US Geological Survey (USGS) and distributed in raster format. The NED represents bare earth elevations and is the best such dataset available for the area under investigation. Unlike older DEMs (e. g., the USGS 7.5-minute and 3 arc-second models) that contain elevation data tied to different spatial reference systems, the NED contains elevations that are referenced only to the North American Datum of 1983 (NAD1983) and North American Vertical Datum of 1988 (NAVD88) -- and they share a common unit of measure (meter). NED elevations are arranged on a grid and spaced at one arc-second intervals, which, at 46° N latitude, spans approximately 31 meters in the north/ south direction and 22 meters in the east/ west direction. When first released, most NED elevations were derived from reprocessed versions of older DEMs, but increasingly they are being derived from new and more accurate source materials. The NED is now assembled completely from 10-meter and 30-meter resolution sources and new tiles at 10-meter resolution are incorporated as they are produced (Gesch et al., 2001). The NED is updated every two months; this work employs the December 2004 release that contains new elevation data for a portion of northern Lower Michigan. 142 DEMs contain errors in elevation that are the cumulative results of imprecise sampling, measurement and interpolation processes (see e.g., Shortridge, 2001; Fisher, 1998; Bolstad and Stowe, 1994; Brown and Bara, 1994 and Shearer, 1990). 3 The NED has not yet been assessed for accuracy, but according to the USGS (2005), ”a plan for assessing the overall accuracy of the National Elevation Dataset (NED) based on independent high-accuracy geodetic control is currently under development.” The USGS (2005) suggests NED users should consult ” published information on the accuracy of the source digital elevation models (DEM's) from which NED was assembled.” Spatial metadata distributed with NED data contain polygons that outline footprints of source DBMS and attribute tables with statistical summaries of error. A root mean square error (RMSE) quantifies the average mismatch between a small set of ground-based and DEM—reported elevations within each footprint. Such values are not specific to particular locations and hence cannot be used to pinpoint spots where elevations can be considered as reliable or in need of improvement (Shortridge, 2001; Kyriakidis et al., 1999). Furthermore, RMSE values represent mismatches calculated prior to NED processing, which 3 Horizontal locations stored in a digital elevation model are controlled during the production process (i.e., they are selected and not measured), so digital elevation models by design do not contain horizontal errors (Maune, 1996). Admittedly, it is possible for an elevation to be assigned incorrectly to a location, but such violations still are considered elevation errors and not horizontal errors. 143 may be different than actual values. The NED data processing protocol includes: a) the application of large filters to remove production artifacts or smooth discontinuities along input tile boundaries; b) bilinear interpolation methods that coerce data from extant resolutions into available spots along the arc second grid; and c) geodetic transformations that change elevations formerly rounded and stored as integers (either in feet or meters) into decimal meters stored with floating point precision (USGS, 2005). Consequently, a potential mismatch exists between metadata-reported and actual RMSE values. So, while historic RMSE values associated with source files may be the best available indicators of accuracy, actual accuracy is uncertain. The NED also contains elevation data for Canada that are situated near the US/ Canada border and within the one-degree tile structure employed by the USGS (Gesch et al., 2001; Gesch et al., 2002). Canadian elevation data were not collected according to USGS protocols, but likely derived from the Digital Terrain Elevation Dataset (DTED: NIMA, USGS and NGDC, 1997). The DTED was created for large areas of the world by the former Defense Mapping Agency during the 19605 and 19705, and used primarily for military purposes. NED metadata provide scant information about how these data were collected or processed, but metadata-reported RMSE values (10 - 29 m) indicate elevation data for Canada contain large errors relative to those associated with other sources (RMSE 1 - 7 m). Consequently, elevation data for Canada are considered a likely source of uncertainty in model results. 144 Geoprocessing the NED NED data were projected onto the Michigan GEOREF planar coordinate grid and resampled to 45-meter resolution using GIS software. I was initially reluctant to reduce the horizontal precision of NED data, but was compelled by computer memory and storage limits (NED data for the study area consume nearly 0.5 GB of disk space at 30-m resolution). Nevertheless, the new resolution is an improvement over the coarser 90—m data employed in related work by Schaetzl et al. (2002) and Drzyzga et a1. (2002). Also, many grid cells associated with Lakes Michigan, Huron and Superior were masked from analysis because they do not represent land surface information. The geoprocessed version of NED data is hereafter referred to as, simply, ” the DEM.” I divided the DEM into three subsets: Lower Peninsula, Upper Peninsula, and Ontario. I divided Michigan data into two subsets because elevation data for the selected portion of the Lower Peninsula were derived primarily from contour lines on topographic maps, and most elevation data for the selected portion of the Upper Peninsula were derived photogrammetrically. Elevation data for Ontario were derived differently than elevation data for Michigan, so I suspected internal inconsistencies might differ as well. The literature has long recognized that differences among DEM production methods yield DEMs with different error characteristics (e. g., Wood, 1996 or Brown and Bara, 1994), so I sought to examine each subset to ascertain whether or not such differences existed wihin my data. I originally considered subsetting the DEM using source tile polygons 145 that were attributed with different production method data, but doing so would have fragmented the study area into many non-contiguous and hard-to-manage units. More importantly, however, the USGS can update any or all data within a source tile with elevations obtained via methods not consistent with those reported in the spatial metadata file (USGS, 2005). Consequently, footprints of NED source tiles, which indicate the extents of pre-NED production methods, may not coincide with the extents of actual production methods. This complicating factor prompted me to use simple majorities to inform my choice of three subsets. Hypsometric analysis Initial assessment of the DEM and each subset can be achieved by visual analysis of frequency histograms, which can reveal internal inconsistencies if any are present (Hutchinson and Gallant, 2000; Wood, 1996). The frequency distribution of my DEM elevations (Figure 4.4a) is skewed toward lower elevations and it contains several conspicuous spikes. The largest spike is situated between 176 and 181 m and represents a large number of ”water” cells associated with the contemporary Lake Michigan-Huron water plane that remain in the model, as well as others associated with inland lakes Burt, Charlevoix and Mullet. A secondary peak is situated between 213 and 218 m and associated with large homogenous groups of cells that represent the elevations of Hubbard Lake (Lower Peninsula) and North Manistique Lake, Manistique Lake and South 146 Manistique Lake (Upper Peninsula). Several lesser spikes are also evident at 213, 244, 274 and 305 m, but these cannot be explained by known topographic regularities. The Lower Peninsula subset is the largest of the three and its distribution (Figure 4.4b) dominates the distribution of the whole DEM as the two histograms are quite similar. High frequencies of 213, 244, 274 and 305 m elevations in the DEM are not evident in this subset. Other and lesser spikes, however, appear at three-meter intervals. Most 1:24,000 USGS topographic maps have elevation contours mapped at ten-foot intervals (roughly three meters), so it seems that artifacts associated with DEM production methods that utilized digitized contour lines have survived geoprocessing and are observable. The Upper Peninsula histogram (Figure 4.4c) contains many low elevations that represent the characteristic lake plain topography of the area, and a small peak centered near 255 m, which represents several bedrock uplands. Two conspicuous spikes occur at 177 and 217 m: the lower is associated with ”water” grid cells representing contemporary Great Lakes water level elevations; the higher is associated with a large number of small inland lakes that occur at various locations across the lake plain of the eastern Upper Peninsula. The Seney National Wildlife Area, for example, is predominantly swampland and most open water bodies situated within it were mapped with the same elevation value. 147 120 73 70‘ 100 O O O O 2 e 80 6 6 60 C C 9 3 40 a a t: if: 20 O 200 300 400 500 200 300 400 500 Elevation (m) Elevation (m) (a) the DEM (b) Lower Peninsula subset 120' 120 13 100 13 100 O O 8 8 ' 5 803 5 80 l 3 60 3 601 i 5 . g g l a . j J l E} 40 E} 40 . i f “.3 93 l u. 20 u. 20 ' 0‘ 0 200 300 400 500 200 300 400 500 Elevation (m) Elevation (m) (c) Upper Peninsula subset ((1) Ontario subset Figure 4.4: Frequency histograms of elevation. Bin width is one meter. Large spikes at 213, 244, 274 and 305 m in the DEM histogram are not evident in the Upper Peninsula histogram. Yet, other and lesser spikes do appear at three- meter intervals. The majority of elevation data in this subset were derived 148 photogrammetrically yet some were derived from contours, which explains the consistent pattern of overrepresented elevations at three-meter intervals. The Ontario subset (Figure 4.4d) is the smallest of the three but contains the full range of elevation values. Topography in the Ontario region is interesting and unique because it contains both stranded lake plains and a high and exposed portion of the Canadian Shield. Yet, irregular Shield topography does not account for the regular spikes that occur at 213, 244, 274 and 305 m or the intermediate spikes at 229, 259, and 289 m. When meters are converted to imperial feet, these spikes reveal a consistent pattern of overrepresented elevations at 100 and 50 ft intervals, respectively, and suggest a systematic and non-topographic cause (Hutchinson and Gallant, 2000; Wood, 1996). Given that 1:250,000 USGS topographic maps of the region were drafted using a 50-foot contour interval, it seems reasonable to speculate that such mapped contours were used as source materials for DTED elevations, which are the likely sources of NED elevations in Ontario. Again, NED metadata contain no information about the sources used to assemble data for locations in Ontario, so this is only speculation. Regardless of source, however, spikes in Ontario elevation data betray a likely artifact of a DEM production process and they indicate the elevation model over-represents contour elevations within this region. 149 Ground control The task of assessing accuracy is concerned with describing errors in data (Shortridge, 2001). Typical assessments: begin with the assumption that each error is the result of a random measurement process; proceed by comparing a set of measured values to a set of known or higher-order values assumed to be accurate; and end with a description of the error distribution. I assembled a pool of high-order check points consisting of National Geodetic Survey (NGS) and Geodetic Division of Canada (GDC) control stations. Stations were selected if they are located within the area of interest and have geodetic-quality location and elevation information. As mentioned above, not all stations have geodetic- quality location information or geodetic-quality height information, so not all can be used for this purpose. The set of control stations was reduced further to exclude those mounted atop manmade objects (e.g., an antenna, lighthouse, range light, dock, bridge etc.) or those not mounted flush with the ground surface. Geodetic control stations are preferentially situated near engineered works (e.g., a railroad, highway, or commercial waterway) or atop accessible high points (pers. obs. ; Kumler, 1994). Consequently, the set of station locations is not representative of the full range of topographic conditions in the region. I added the set of GPS-measured paleoshoreline data to the pool of check points so that locations distant from engineered works can also be assessed. Paleoshoreline data are suitable for this purpose because they contain sub-meter location errors, typically contain sub-meter errors in elevation, and are, therefore, 150 of a higher-order than DEM data. Like control stations, paleoshoreline positions were also preferentially selected and so, for example, do not include hill tops or river valleys, but their inclusion adds important landscape variation that is missing from the set of geodetic control station locations. Moreover, adding paleoshoreline positions to the pool enables me to ensure that agreement exists between the surface elevations I observed in the field and the surface elevations stored in the DEM (see Chapter 5). A pool of 529 check points (256 NGS control stations, 4 GDC control stations and 269 paleoshoreline positions) was used to quantitatively assess DEM elevations. A comparison of frequency histograms reveals that the distribution of check point elevations (Figure 4.5a) is shaped and skewed like the distribution of DEM elevations (Figure 4.5b) and suggests that the sample of check point elevations is representative of the larger sample held by the DEM surface. The check point distribution lacks, however, large spikes in the lowest elevation bins because the elevations of open water bodies, which are represented in the surface, are not represented in check point data. 151 80 200 29? ”E 60 o 150 V O >, V— 0 v 5 40 5‘ 100 3 C 0' (D e 3 u— 20 E 50 LL 0 0 200 300 400 500 200 300 400 500 Elevation (m) Elevation (m) (a) (b) Figure 4.5: Histograms of: a) 529 check point elevations; and b) all DEM elevations in the study area. The difference between frequency scales is the results of unequal sample sizes. Select results of the first exploratory analysis Ultimately, I conducted two exploratory data analyses of error. The first led me to find a latent horizontal shift in the original NED- a systematic shift that propagated through geoprocessing unnoticed and perturbed my attempts to assess the distribution of error. Rather than report complete sets of results for both analyses, however, I use this section to highlight interesting characteristics of the original DEM that betrayed evidence of the shift. I use the next section to report a complete accuracy assessment of the corrected and properly-situated DEM. Original DEM elevations were assigned to control point locations using nearest-neighbor interpolation and GIS software. A scatterplot (Figure 4.6a) 152 suggested DEM elevations generally agreed with actual elevations. Results of a Kendall’s test for statistical association (a = 0.10: coefficient = 0.897, p < 0.0000) indicated a strong and significant correlation existed between the two sets. 450 5 450 g 400- g 400 5 C C .9 .Q 5 ‘ 5 . 7;, 3001 E 300; E 250; E 250; o 5 o . 200* 200? 200 300 400 200 300 400 Actual elevation (m) Actual elevation (m) (a) Original DEM (b) Corrected DEM Figure 4.6: Scatterplots of DEM-reported and actual elevations of the (a) original DEM and for (b) the corrected DEM. I estimated errors in DEM elevations by subtracting check point elevations from interpolated DEM elevations and then calculated summary statistics to describe the distribution. Table 4.5 contains the first and second sets of summary statistics. For example, the original sample mean difference (3? = -0.55 m) was significantly different than zero at the 90% confidence level. 153 Table 4.5: Summary statistics that describe differences between actual and co- located elevations in the original and corrected DEM surfaces. The original DEM The corrected DEM n 529 529 minimum -26.2 m -18.6 m median -0.30 m ~0.02 m mean (bias) -0.55 m -0.36 m standard error of mean 0.27 m 0.18 m standard deviation 6.79 m 4.07 m maximum 27.8 m 19.8 m Ultimately, I investigated possible covariates with error and was surprised to find significant associations with the linear components of slope aspect.4 Circular measures of aspect {0 S aspect < 360 degrees} were decomposed into the linear U (east-west) and V (north-south) components of a unit vector {-1, ,1} to facilitate linear statistical analysis. A pure south-facing slope, for example, has a V component value of -1, and a pure north-facing slope has a V component of +1. Scatterplots (Figure 4.7) and results of Kendall’s tests for association both suggested error was positively and significantly correlated with the east-west component (tau = 0.199 and p < 0.0000), and negatively and significantly correlated with the north-south component (tau = -0.287 and p < 0.0000). 4 I derived a surface of slope aspects from the grid of DEM elevations and assigned co- located values to control points using GIS software. 154 Error (m) -1.0 -0.5 0.0 0.5 1.0 -1.0 -0.5 0.0 0.5 1.0 Aspect (u component) Aspect (v component) (a) original DEM (b) original DEM 30 30 20 . 20 E ° ‘2 Lu -10 UJ -20 —20 —30 -30 —1.0 —0.5 0.0 0.5 1.0 —1.0 —0.5 0.0 0.5 1.0 Aspect (u component) Aspect (v component) (c) corrected DEM (d) corrected DEM Figure 4.7: Scatterplots of DEM error and the linear components of slope aspect. Although the statistical relationships were weak, the graphical evidence was striking (Figure 4.7a and 4.7b) and indicated elevations in the original DEM tended to overestimate elevations at locations with east- or south-facing slopes, 155 and underestimate elevations at locations with west- or north-facing slopes. Such results indicated clearly that the DEM was not georeferenced prOperly. I visually inspected control point and DEM data simultaneously using GIS software and observed the apparent shift for the first time. Multiple attempts to estimate the correct locations of control points within the footprints of nearby grid cells yielded consistent estimates. On average, elevations in the original DEM were shifted 75 m south and 15 m east of control point locations. So, I corrected the DEM by translating the origin of the original NED surface 75 m north and 15 m west, and reprocessed the model according to the methods described above. Hereafter, any mention of ”the DEM” employed in this work refer to the correctly georeferenced elevation surface. Corrected DEM elevations were assigned to control point locations using nearest-neighbor interpolation and GIS software and I conducted a second exploratory analysis that yielded markedly improved results. For example, the association between error and slope-aspect is now considered random, uncorrelated and not indicative of any systematic bias (see Figure 4.7c and 4.7d). The association between DEM-reported and actual elevation is stronger (see Figure 4.6b) and suggests a better match (a = 0.10: Pearson’s correlation coefficient = 0.9978, p < 0.0000), the sample mean and median values moved toward zero, the two extreme error values have much lesser magnitudes, and the sample variance diminished (see Table 4.5). These new results indicate that the DEM is georeferenced more accurately than before. 156 Complete results of the second exploratory analysis Again, corrected DEM elevations were assigned to control point locations using nearest-neighbor interpolation and GIS software. Errors in elevation were estimated by subtracting check point elevations from interpolated elevations. A map of data postings (Figure 4.8) reveals no conspicuous trend or drift. Yet, many small errors occur near other small errors, many positive errors are situated near other positive errors and large negative errors tend to occur near other large negative errors (although, in some instances, large positive errors are situated near large negative errors). In short, differences from the sample mean (-0.36 m) seem to exhibit spatial continuity or structure. Table 4.6 contains statistical summaries of error, Figure 4.9 shows frequency histograms and Figure 4.10 shows q-q plots for the entire DEM and each subset. The calculated mean of all DEM errors ()7 = -0.36 m) is significantly different than zero at the 90% confidence level. Summary statistics and the frequency histogram (Figure 4.9a) suggest the distribution of error is bell-shaped, symmetric and centered near zero. The q-q plot (Figure 4.10a) reveals, however, larger-than-expected errors in both tails. Results of a Shapiro-Wilk test for normality (a = 0.10: W = 0.917 and p < 0.0000) indicate the standardized distribution of error is generally correlated with, but does not significantly fit a Normal distribution. 157 I:- :? _ I" . O 650 1 ’1'”? a: ”(9 r, . ,1 ‘ . o ”"1" u) 7. g H ‘ 3‘ . . ‘.‘ . ‘ '1 . 6 . a n 600 9 A _ 11 . D > 8 m E 0 4to 8 i5 3 -. '. g 0 to 4 O) O 0.9;? a .4. " ‘4 t0 0 E 550 . «a; o' 0'31 - -8to-4 E ’4: ‘—'~ :0 <-8m Z . .,1 I F - l. O . "'w :1 O f; ,. 1 500 .‘ Q J . . j .1? .1 . 2 H. .a . a _‘ O 0' . e .211 + ‘ 4 . 1:. a 450 1 1. . .. ‘ O Q . l . 3 . 500 550 600 650 700 750 Easting (km) Figure 4.8: Map of DEM errors at control point locations. The map reveals no global pattern but suggests error values exhibit some spatial continuity at short distances. The image is presented in color. Elevations in the Lower Peninsula subset were collected originally from contour lines either digitized or scanned from topographic maps. The Lower Peninsula sample mean of calculated errors (f = -0.06 m) is not significantly different than zero at the 90% confidence level. 158 Table 4.6: Summary statistics that describe differences between actual and co- located DEM elevations. The entire DEM Value Test for mean 95 0 Value n 529 minimum -18.6 m median -0.02 m mean (bias) -0.36 m matched-pairs t p = 0.043 standard error of mean 0.18 m standard deviation 4.07 m maximum 19.8 m Lower Peninsula subset n 348 minimum -13.3 m median 0.05 m mean (bias) -0.06 m matched-pairs t p = 0.718 standard error of mean 0.16 m standard deviation 2.94 m maximum 13.1 m Upper Peninsula subset n 108 minimum -13.6 m median -0.30 m mean (bias) -0.50 m matched-pairs t p = 0.147 standard error of mean 0.34 m standard deviation 3.57 m maximum 7.0 m Ontario subset n 73 minimum -18.6 m median -2.07 m mean (bias) -1.58 m matched-pairs t p = 0.081 standard error of mean 0.90 m standard deviation 7.65 m maximum 19.8 m 159 150 150 E E >. > g 100 1 g 100 a: a) D D a a u“. 50 E 50 0 O -20 -10 O 10 20 -20 -10 0 10 20 Error (m) Error (m) (a) entire DEM (b) Lower Peninsula subset 150 150 E E 5 6 C 100 . C 100 0: c0 3 D 8 8 u: 50 I: 50 0 :41: ".‘l 0 ..... -20 -10 0 10 20 -20 -10 0 10 20 Error (m) Error (m) (c) Upper Peninsula subset ((1) Ontario subset Figure 4.9: Histograms of elevation errors in the DEM. Summary statistics and a frequency histogram (Figure 4.9b) suggest the distribution of error is bell-shaped, symmetric and centered near zero, but may contain larger-than-expected positive errors. 160 Sample quantiles O Sample quantiles O -2 -2 a" _4 .. ~‘,_ -4 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 Theoretical quantiles Theoretical quantiles (a) the entire DEM (b) Lower Peninsula subset 3 2 , g 0 ' g 1 3 r 3 o- , o- _ 2 ’1 :2 0 a ' a. E —2 -- E “3 - m —1 : (I) a) . -3 ' -2 -2 -1 0 1 2 -2 -1 0 1 2 Theoretical quantiles Theoretical quantiles (c) Upper Peninsula subset (d) Ontario subset Figure 4.10: Q-Q plots of standardized elevation errors in the DEM. The Lower Peninsula q-q plot (Figure 4.10b) reveals, however, larger-than- expected errors in both tails. Results of a Shapiro-Wilk test for normality (a = 0.10: W = 0.930 and p < 0.0000) indicate the distribution or error in the Lower Peninsula does not fit a Normal distribution. 161 Elevations in the Upper Peninsula were collected primarily using photogrammetric techniques. Bolstad and Stowe (1994) compared photogrammetrically-derived elevations at 42 NGS bench marks near Blacksburg, VA and reported such values tend to underestimate actual elevations. The sample mean of calculated errors ()7 = -0.50 m) for the Upper Peninsula subset is less than zero and so suggests this model also tends to underestimate actual elevations, but the sample mean is not significantly different than zero at the 90% confidence level. A histogram (Figure 4.9c) suggests the distribution of error is bell-shaped, but has large errors in the lower tail. A q-q plot (Figure 9.10c) indicates the lower half of the distribution is not well-behaved for it contains several values that are larger than expected should errors be distributed as Normal. Results of a Shapiro-Wilk test for normality (a = 0.10: W = 0.930 and p < 0.0000) indicate the distribution of error in the Upper Peninsula does not fit a Normal distribution. Elevations in Ontario, Canada were likely derived from the DTED and coerced into the NED. Summary statistics and a frequency histogram (Figure 4.9d) suggest the distribution of DEM error is skewed left and centered below zero. The Ontario subset contains extreme values and they are the same as those associated with the entire DEM. The sample mean of calculated errors (J—C = -1.58 m) is significantly different than zero at the 90% confidence level and the sample standard deviation (7.65 m) is relatively large. As expected, elevation data for Ontario contain large errors and so should be considered a likely source of 162 uncertainty in model results. A q-q plot (Figure 4.10d) suggests that the distribution of error might not be distributed as normal, but results of a Shapiro- Wilk test for normality (a = 0.10: W = 0.982 and p = 0.364) indicate the distribution or error does fit a Normal distribution. I partitioned the DEM into three subsets based on ancillary information about the original DEM production methods. I expected to find marked differences in the resultant behaviors of their elevations and errors. It seems, however, that subsequent geoprocessing mitigated (or obfuscated!) many of the apparent differences, for only the Ontario subset appears stands out with its relatively large mean and variance. 1 transformed this observation into a hypothesis and conducted three, two-sample difference-of—means tests to assess statistical differences (a = 0.1, m = 0.0 m and assuming unequal variances). All test results (Table 4.7) indicate, however, that no statistical difference exists between the means. Nevertheless, I continued to work with the three subsets because each contains data derived from different source materials. Table 4.7 : Results of two-sample difference-of-means tests among DEM subsets. Lower Peninsula Upper Peninsula Upper t =1.164 Peninsula df = 107 p value = 0.233 Ontario t = 1.612 t =1.126 df = 72 df = 72 p value = 0.201 p value = 0.274 163 Covariations with location, elevation, and slope Bolstad and Stowe (1994) reported large discrepancies between ground- based and DEM-reported elevations tend to occur along steep slopes and ”in high-relief forested areas” (p1329), thus implying DEM errors exhibit spatial continuity and covary with other landscape elements. Such ideas about spatial continuity are not unique for others also have reported spatial patterns of error (e.g., Holmes et al., 2002, Kyriakidis et al., 1999; Hunter and Goodchild, 1997; and Veregin, 1997), but they are important because they provide grist for examining error in terms of location, elevation and other landscape attributes. Specifically, Bolstad and Stowe (1994) claim dense tree canopies negatively impact the precision of some DEM production methods, and photogrammetric image correlations can fail in areas with steep slopes. Such imprecision and failures yield inaccurate data for those locations. From this information, I plotted maps (not shown), reviewed field notes and examined datasheets for check points associated with large errors (> | 8 ml). Indeed, I found many that occur in high- relief or forested areas (NGS datasheets contain descriptions of the approach and area near control stations). Station R]1215 (error = -13.6 m) is situated atop a high and forested knoll, station R]1118 (error = -12.0 m) is mounted along a large outcrop of limestone, and station R]0383 (error = -13.6 m) is situated near the peak of a high and steep hill of sand. Station PL0631 (error = -8.8 m), however, is situated on flat and open ground. 164 Of the thirty-four control points that have large absolute errors, only four are NGS control stations; the others are paleoshoreline sites. Had I not already assessed the quality of my paleoshoreline data (mean = 0.16 m, standard deviation = 0.77 m), then I might have considered this result to be an indictment of the lesser quality of my paleoshoreline data. Instead, I believe that many of the large errors (> I 8 ml) are associated with ground slope and data model imprecision. Elevations in the DEM grid are spaced at 45-meter intervals that span terrain discontinuities and cannot capture finer-scale variations in topography (Maune, 1996; Carter, 1992; Shearer, 1990). Many relict lake bluffs, for example, have steep slopes that rise from toe to top over runs less than 45 meters - all within the footprint of a single grid cell. Consequently, the DEM will, by design, contain errors that covary with such breaks or changes in slope. Check-points have location, elevation and error attributes so relationships among these variables can be assessed easily. To assess the relationship between error and slope, however, I derived surfaces of slope gradients and slope aspects from the grid of DEM elevations and assign co-located values to control points using nearest-neighbor interpolation and GIS software. Any comparison between an elevation error and an elevation-based derivative, like slope gradient or aspect, presupposes extant errors in elevation do not influence the slope values derived for the same locations. Such a presupposition is unreasonable, especially since Carter (1992) demonstrated how gradient and aspect calculations 165 are sensitive to even small elevation errors, but I continued for the sake of exploration. Table 4.8: Results of Kendall’s tests for association between easy-to-measure landscape variables and: a) signed error and b) unsigned error. a) signed error b) unsigned error The entire DEM tau p value tau p value Easting 0.021 0.476 0.126 < 0.001 Northing -0.037 0.204 0.178 < 0.001 Elevation -0.053 0.067 0.074 0.011 Slope gradient 0.033 0.256 0.273 < 0.001 Slope aspect (U) 0.004 0.902 -0.007 0.804 Slope aspect (V) 0.051 0.078 0.035 0.227 Lower Peninsula subset Easting 0.108 0.003 0.018 0.621 Northing -0.003 0.928 0.085 0.018 Elevation -0.062 0.083 0.021 0.563 Slope gradient 0.039 0.275 0.308 < 0.001 Slope aspect (U) 0.024 0.498 0.019 0.598 Slope aspect (V) 0.140 < 0.001 0.050 0.160 Upper Peninsula subset Easting 0.050 0.444 -0.032 0.619 Northing 0.021 0.750 -0.034 0.604 Elevation -0.062 0.344 0.176 0.007 Slope gradient 0.076 0.244 0.331 < 0.001 Slope aspect (U) 0.043 0.512 -0.017 0.799 Slope aspect (V) -0.145 0.026 0.099 0.128 Ontario subset Easting -0.424 < 0.001 0.174 0.029 N orthing 0.287 < 0.001 -0.207 0.010 Elevation -0.006 0.943 -0.046 0.567 Slope gradient -0.168 0.036 0.051 0.520 Slope aspect (U) 0144 0.071 0.034 0.675 Slope aspect (V) -0.003 0.970 0.057 0.478 166 I used scatterplots (not shown) and Kendall’s non-parametric test for association to examine the relationships between error and location, elevation and slope. I examined unsigned errors as well. Table 4.8 lists the results. A strong association between variables is indicated when there is a large probability that the ranks of data are in the same order (tau > 0.667). Of all the relationships tested, none exhibited a strong association. Jumbled pairings (tau < 0.333) or statistically insignificant test results (p > 0.100) eliminate the need to consider most of these relationships further. A few, however, exhibit moderate (0.333 < tau < 0.667) and significant (p < 0.100) associations, so I focused on them. Unsigned error is positively correlated with slope gradient in every set, although results for the Ontario subset are weak and not statistically significant. These results are reasonable because they suggest large errors occur in locations where elevation is hard to measure -- along steep slopes -- and they agree with results reported by Bolstad and Stowe (1994) and Hunter and Goodchild (1997). A scatterplot of this relationship offers conflicting evidence (Figure 4.11) however, for it shows the full range of error occurs at locations with relatively shallow slopes (< 5% gradient). Also, both point clouds show no linear patterns of association. The relationship between error and slope gradient provides grist for general conclusions only. Carter’s (1992) work cautions that gradient calculations are uncertain given errors in the elevation model. From these results, steep slopes can be associated with large DEM errors, but both large and 167 small DEM errors can occur along any slope regardless of gradient. Therefore, an error estimate cannot be reduced if the terrain is known to have a low gradient. The only other moderate correlation among error and location is in the Ontario subset and it indicates error decreases in an easterly direction. Review of ground control data, however, provides a reason for why this test result might not be a valid result. Ground control data used to assess the Ontario portion of the DEM are highly clustered (Figure 4.12a). They contain only four GDC control stations - the balance are paleoshoreline positions collected at field sites. Paleoshoreline data near Sault Ste Marie, Ontario, for example, were collected by Cowan (1985) and his data comprise the northwestern cluster. Figure 4.12b reveals most DEM elevations tend to overestimate elevations in this cluster. Paleoshoreline data at St. Joseph Island were collected by Karrow (1987) and his data comprise the central cluster. DEM elevations tend to overestimate Karrow’s elevations, but not always. And, I collected paleoshoreline data at Cockburn Island, which comprise most of the southeastern cluster. DEM elevations tend to underestimate elevations in this cluster, but not always. In short, the test for association between error and Easting is seemingly influenced by the interrelations among three distant clusters and not among seventy-three well- distributed observations (Figure 4.12b). Consequently, the test statistic has likely been overestimated and the associated p value is too small. I did not to pursue this association further. 168 20 “1'- 20 E a) 15 g A 'D l E g 9'7 .5 g, 10 t m m E ‘9' 5 LE -20 0 0 5 10 15 20 0 5 10 15 20 Gradient (percent) Gradient (percent) (a) (b) Figure 4.11: Scatterplots of: a) signed error and slope gradient; and b) unsigned error and slope gradient. Black lines drawn to emphasize zero error. 20 660 - :E: . . 7.1.. ‘_ , A 10 . C ‘, h " . . 11.. 1 1 , t 620 -. ~_." - - a Li , l, 1112*» g o 1 1 . 4:. '1‘ . .7 at. Z . _. " - -. . -1o ‘ 1 -' @4155. - "‘1‘ ‘ ,1 -=‘.‘x,,, : . > 1 j, . ~1~ ' {r 1, fi —20 620 640 660 680 700 620 640 660 680 700 Easting (km) Easting (km) (a) (b) Figure 4.12: Map and scatterplot of ground control locations used to assess error in the Ontario subset. 169 Spatial autocorrelation The distribution of DEM error (Figure 4.8) seems to exhibit some form of local spatial continuity or structure - meaning two errors located close to each other tend to have more similar values than two errors that are far apart (after Isaaks and Srivastava (1989:50—51). This is structure that cannot be explained by the weak associations with location, elevation or slope. The semi-variogram is a useful tool for exploring and characterizing the way in which the deviations from a mean value at different locations covary or are correlated (Bailey and Gatrell, 1995: 162). Yet successful use and interpretation of a semi-variogram of errors presupposes an understanding of regionalized variables and the assumptions of stationarity. These concepts are much more closely aligned with the spatial models discussed in Chapter 5 than with the accuracy assessments presented here, so I report these results in Chapter 5. Summary I designed my sampling and data processing protocols (Chapter 3, this volume) to ensure that paleoshorelines were measured as accurately as possible given the lesser precision of GPS technology relative to traditional yet more labor-intensive technologies. Careful planning and diligent use has yielded a set of measured positions that closely match actual positions. Measurement errors in elevation are random, typically sub-meter and always less than two meters; the distribution is centered near zero and Normal. Horizontal errors are 170 negligible. There is no evidence of any remaining gross or systematic error that can or should be mitigated. Selection errors in paleoshoreline data present a problem because they cannot be easily identified or redressed. Gravel lags are the best indicators of former shorelines, but using them as indicators of shoreline elevation is imprecise. I suspect such errors are often ignored during spatial modeling and analysis. In this case, however, the range of possible selection errors is large enough to exceed the three—meter critical threshold I adopted for distinguishing different water levels, and so they warrant further consideration. A position with a three-meter selection error is not acceptable for use, but I do not know if any of my data contain such an error. Consequently, my only recourse is to add the information I do have about possible selection errors to my spatial model. In short, I need to simulate their existence in order to ascertain how they might impact my results. The US Coast Guard archive of GPS Status Message Reports (US Coast Guard Navigation Center, 2004) is a valuable resource; I did not know it existed until after I collected data and began analysis. The US Coast Guard posts status reports that list planned outages and system changes as well as recent and unexpected satellite behaviors (such as the one I encountered). This kind of information will impose additional constraints on future planning missions, but it will add value to any data I collect with GPS in the future. 171 The NED is the best available source of ground surface elevation data for the region of interest. Unfortunately, these data contain both large and small errors that are the results of a long and complex geoprocessing history. Analysis of the NED has provided the biggest surprise for I did not expect to find a systematic bias in location. I translated the NED surface 75 m north and 15 m west and, subsequently, removed a systematic pattern of elevation errors and reduced the sample mean, median and variance of errors. Simultaneous visualization of the corrected surface and check point data provides solid evidence of a better fit. I have since revisited all NED-related documents on the USGS website and scoured the literature for any report that might address the shift I observed, but I have found no report that does so. I forwarded my results to the USGS, but have not yet received a reply. Although corrected, the DEM still holds a pool of variations associated with topography, artifacts of DEM production methods, and adulterations imposed by various geoprocessing algorithms. Partitioning the model into three subsets has proven useful for it has allowed me to isolate some of these variations and consider their impacts. The distribution of elevations in the Ontario subset, for example, is clearly dominated by artifacts of a data production process. Had I not subset the model, the non-topographic pattern of spikes in the Ontario histogram would have remained largely hidden in the histogram of DEM elevations. Elevations in the Lower and Upper Peninsula subsets also exhibit artifacts of data production processes, but these appear mild 172 by comparison. By and large, the Lower and Upper Peninsula subsets contain elevations that represent surface topography. Mapped errors in DEM elevations exhibit evidence of local spatial continuity as positive errors tend to occur near other positive errors and negative errors tend to occur near other negative errors. Yet, correlation test results demonstrate DEM error, either signed or unsigned, does not covary with location coordinates, elevation, or slope. Moreover, trend surface analyses (not reported) provided no statistical power to suggest otherwise. The apparent association between unsigned error and slope gradient is acknowledged but the parameters that might define its behavior are uncertain because gradient calculations are so sensitive to errors in elevation (Carter, 1992) and the DEM employed in this work certainly contains errors in elevation. In my opinion, no correlation is sufficiently strong to advocate using knowledge about it to improve estimates of error at unsampled locations. While it would be convenient to have found a strong and significant association between error and location, elevation, or slope - or any other cheap-to—measure landscape variable for that matter - other means of qualifying spatial continuity are necessary. I take up this issue again in Chapter 5. Absolute errors in DEM elevations often exceed three meters and, in Ontario, can reach almost twenty meters. Such errors far exceed the critical threshold I’ve adopted for distinguishing different water levels. The DEM, as is, is not fit for my use. Unfortunately, the current DEM is the best one available, 173 and so I must use it. My only recourse is to add the information I do have about DEM errors to my spatial model. In short, I need to fix the elevations errors where I know they occur, and I need to simulate the presence of error at locations between check point locations. Doing so will enable me to ascertain how they might impact my results. 174 CHAPTER 5: THE SPATIAL MODEL Introduction This chapter describes a spatial model for digitally reconstructing ancient landscapes associated with glacial Lake Algonquin. The model has two parts. The first part reconstructs water planes that have been upwarped by differential isostatic rebound over the past 11,000 years. Initial reconstructions are supported by paleoshoreline data (described in Chapter 3) and parameters that describe the upwarped surface are estimated via least squares. Water planes are reconstructed ultimately via universal kriging, which yields both a prediction surface and an estimate of uncertainty at each predicted location. The second part of the spatial model addresses error and uncertainty about error in a contemporary digital elevation model (DEM), which is used to represent ground surface form for the region of interest. A set of control points and information about elevation errors at control point locations is used to condition the DEM and, via ordinary kriging and Monte Carlo simulation techniques, to generate alternate possible realizations of the ground surface. Finally, results of the two parts (a predicted water plane and a set of alternate possible DEMs) are combined to allow locations within the area of interest to be characterized as: above water plane; below water plane; or uncertain. Four landscapes are reconstructed and they are associated with the Main, Ardtrea, Upper Orillia, and Lower Orillia phases of glacial Lake Algonquin. 175 I offer a review of geostatistical concepts that underlie general spatial prediction and kriging before presenting the spatial model and its two parts. My review summarizes essential points related to the theory of spatial prediction, but does so in a user-friendly format. Formal treatments of these theories and methods are presented by Journel and Huijbregts (1978), Isaaks and Srivastava (1989), Bailey and Gatrell (1995) and Gringarten and Deutsch (2001). Modeling outcomes of a spatial random process A field is a conceptual model that is useful for representing an attribute of a phenomenon that varies over space. Elevation, for example, is an attribute of terrain and it varies over space. Similarly, relative humidity and air temperature are attributes of the atmosphere that vary over space. A field model typically represents space by employing a tessellation of discrete spatial elements (e.g., points, grid cells, facets etc.) that are assigned values taken from the phenomenon of interest. This work utilizes a two-dimensional tessellation of discrete grid cells. Field models are always approximate because ”we cannot store an infinite number of locations” (Cova and Goodchild, 2002:511) and we cannot collect or store attributes with infinite precision. Any procedure or process used to take measurements of a phenomenon has finite precision, so data obtained from measurements are imperfect. The difference between a measured value and the actual value is error. The amount of error in any measured value is always unknown, however, because the true 176 value exists beyond the limits of equipment or procedural precisions. Because of this inescapable error, the amount of error in a measured value is considered a random variable. An elevation value obtained from measurement is also considered a random variable because it contains a component of fixed truth and a component of random error. A process that generates a random variable is a random process. Accordingly, a process for measuring elevations at locations over a region of interest is considered a spatial random process. Other spatial processes can influence the outcomes of a spatial random process. Consider, for example, an instance of terrain that is flat, easy to access during calm weather and, so, relatively easy to measure given equipment provided (i.e., conditions for taking measurements are optimal). An elevation value taken under such conditions might be expected to contain random error associated only with variation that occurs beyond the limits of equipment or procedural precision. Alternatively, consider an instance of terrain that is rough, difficult to access during stormy weather and, so, very difficult to measure given the equipment provided (i.e., conditions are not optimal). An elevation value taken under such conditions might be expected to contain an error greater than that generated by equipment precision alone because external forces upset the measurement process. Using these two instances as ”straw men,” it is easy to imagine spatial gradients of terrain (from flat to rough) and weather (from calm to stormy) and how they might interact across space and time and, so, differentially perturb the measurement process from one location and one 177 moment to a next. Consequently, a random measurement process can generate a field of random variables that contain noisy errors associated with equipment precision as well as spatial patterns of error introduced by external phenomena that differentially upset the process. Unfortunately, we do not always know where or when a spatial error process (or processes) exists or how such a process might alter our best efforts to collect data in the field or lab. A random field is a conceptual model that is useful for representing the outcomes of a random spatial process. 1 Use of a random field assumes variation over space can be expressed as the sum of three components: a deterministic component akin to either a constant mean or spatial trend (J ournel and Huijbregts, 1978: p40, prefer the term ”drift”); a random and spatially structured error component; and a random but spatially unstructured error component (i.e., noise). In short, the outcome at each location is considered a random variable, the set of all random variables across the region of interest comprises a random field, and a specific set of measured outcomes within the region of interest is 1 Many geostatisticians use the phrases regionalized variable (Journel and Huijbregts, 1978; Burrough, 1991) or random function (see Kyriakidis et al., 1999; Goovaerts ,1997; or Issaks and Srivastava, 1989), but I prefer to use the phrase random field (see Diggle and Ribeiro Jr, 2002; or Ehlschlaeger, 2002) because it provides an explicit link to the field data model (Couclelis, 1992; Cova and Goodchild, 2002), and doing so averts confusion with other uses of the term function. 178 considered a realization of that random field. Following Burrough (1991), a random field is expressed as: Z(u) = m(u) + e’(u) + e”, Vu E R [Eq. 1] where Z(u) is the random field, m(u) is either the constant mean or spatial trend component, e’(u) is the random and spatially structured error component, 8” is the random and spatially unstructured error component, and u is the set of all locations in the region of interest R. Let u..-, i = 1,. . .,n, represent a set of 11 locations where the random field Z(u) has been observed in region R, and let uo represent unsampled locations in region R. The probability density function (pdf) of a random field Z(u) is a model of possible values that 2 can take at locations u and their associated probabilities: f(u; z) = Prob { Z(u) = z }, V 2 [Eq. 2] Because the probability at a single point in a continuous distribution is always zero, a pdf is, for practical purposes, often expressed in terms of the integral between two points in the distribution. The cumulative distribution function (cdf) is a model that defines the probability that outcomes of Z(u) do not exceed a given threshold 2 : F(u; z) = Prob { Z(u) S 2 }, V 2 [Eq. 3] 179 If additional information exists and can be used to condition the random field, then the conditional cumulative distribution function (ccdf) of the random field Z(u) defines the probability that outcomes do not exceed threshold 2 given this information: F(u; z | (info)) = Prob { Z(u) S 2 | (info) }, V 2 [Eq. 4] The expected value or spatial trend m(u) and spatially unstructured error 8” components are models covered in many statistical texts, so another review here would be redundant. I focus instead on the basis geostatistics provides for exploring and modeling the random and spatially structured component £’(u). Spatial structure implies a relationship between proximity and similarity of outcomes: two outcomes located close to each other will tend to have values more similar than two outcomes that are farther apart (after Isaaks and Srivastava, 1989). According to Bailey and Gatrell (1995), if we expect a random field to exhibit spatial structure, then we can anticipate ”observing positive covariance or correlation at short distances and lower covariance or correlation at greater distances” (p 162). Covariance C() is a measure of similarity between outcomes u.- and u; and can be expressed as: C(ui, “1) = El (Z(ui) - m(w» ' (Z(Ui) - m(ul» } [1361- 5] 180 Modeling the spatial component of a random field often requires adopting one or more assumptions when complete knowledge is not available for all locations within the region of interest. For example, it is not always possible to estimate the cdf at a given location because, typically, at most only one observation is available at that location (Shortridge, 2001). So, we must adopt assumptions about the cdf if we want to proceed with sparse data. The assumption of stationarity, for example, holds the expected value and variance of the random field Z(u) are independent of location, or m(u) = m and 02(u) = 02. Adopting this assumption allows the cdf at any location to be characterized by the cdf of the random field. A more commonly used but weaker assumption is intrinsic stationarity, which assumes a constant mean and a constant variance in ”the differences between values at locations separated by a given distance and direction” (Bailey and Gatrell, 1995:162). In geostatistical notation: E [Z(u) - Z(u + h)] = 0 [Eq. 6] Var [Z(u) — Z(u + h)] = E 1 [Z(u) — Z(u + 1012} [Eq. 71 where h is a vector of distances (or distances and directions) that separate pairs of outcomes. Adopting the weaker assumption removes focus from the absolute locations of data and estimating joint cdfs that are supported by few observations and, rather, places focus on the more-numerous distance vectors that separate pairs of random variables. In this sense, spatial structure implies two outcomes 181 located close to each other will tend to have more similar values than two outcomes far apart (after Isaaks and Srivastava (1989:50-51) - regardless of absolute locations. Covariance between outcomes (Eq. 5) can be re-characterized in terms of separating distances: C(11)= E {Z(U) ' Z(u + h)} - E{Z(U)} '1“3{Z(u + h)} [Eq- 8] or, by assuming intrinsic stationarity and substituting a constant mean, m, as: C(h) = E [Z(u) — Z(u + h)] — m2 [Eq. 9] Note, stationarity and intrinsic stationarity are not properties of any phenomenon, but modeling decisions that provide reasonable means of analyzing a random field given sparse data or scant information about data. A variogram 2y(h) is the expected squared difference between data values separated by vector h (Gringarten and Deutsch, 2001). A variogram can be expressed as: 2101) = Var 12(10— Z(u + h)] = E 1 [Z(u) - Z(u + h)] 21 [Eq. 101 The semi-variogram y(h) is one-half the variogram 2y(h) and a measure of variability used commonly to explore spatial structure in data. Semi-variance 182 values increase as outcomes become dissimilar and tend towards zero as outcomes become similar. Again, covariance C() is a statistical measure of similarity and, under the assumption of intrinsic stationarity, inversely related to semi-variance. By definition, the covariance at h = 0, or C(O), is the variance 0?- (Gringarten and Deutsch, 2001). Covariance C(h) is 0 when values separated by h are not correlated. Combining these relations and the concepts embodied in Equations 9 and 10, the link between semi-variance and covariance is: 100 = C0) - C(h) or C(h) = 02-101) 024.111 The link between semi-variance and covariance provides the basis for interpreting a semi-variogram. The ”sill" of a variogram represents the variance 02 of the random field and it corresponds with no spatial correlation. If Z(u) and Z(u + h) are positively correlated, then corresponding semi-variance values are less than the variance, or sill. If Z(u) and Z(u + h) are negatively correlated, however, then corresponding semi-variance values are greater than the variance, or sill. The ”range” of a semi-variogram is the distance h at which observed semi-variance reaches the sill and spatial correlation no longer occurs. Also, a ”nugget” effect might be observed at very short distances. A nugget effect is interpreted as random chaotic or nonlinear variations ” that have no spatial autocorrelation structure” (Gringarten and Deutsch, 2001: 509) 183 For computation, the semi-variogram is expressed as: n(h) 11h > = 2751—1 21201-1 z1u. +1012 1124121 where n(h) is the number of outcomes separated by each distance (or direction) in vector h. Given sparse data, however, values of h typically represent lags or small ranges of distances approximately equal to h (i.e., within some tolerance). When used within a kriging framework, the theory of random fields and spatial structure can be employed to predict outcomes at unsampled locations (i.e., interpolation). According to Bailey and Gatrell (1995), if covariance exists between observed outcomes, then outcomes at unsampled locations ”are not entirely unpredictable" (p183). A kriging framework accomplishes spatial prediction by modeling unsampled outcomes as weighted linear combinations of sampled outcomes. In practice, kriging weights are assigned with respect to: a) a model of covariances observed among sampled outcomes; and b) a model of covariances expected between sampled and unsampled locations. Kriging weights are subject to linear constraints that ensure the expected squared prediction error is zero. In geostatistical notation, a set of kriging weights is: 401,0) = C‘1c1u.-,o) [Eq. 131 184 where Mum) is the set of weights; C —1 is the inverse of a matrix holding covariances between sampled outcomes; and C(ui,0)is the set of expected covariances between sampled and unsampled locations uo. In practice, expected covariances C(h) between sampled and unsampled locations are modeled as a function of distance between u0 and nearby u ,- . Given a sample, distances between sampled and unsampled locations, and information about the spatial structure of covariance over space, the set of predicted outcomes at unsampled locations becomes: 11 Z(uo)=211u.~,o)Z(u.-) 1124141 i=1 A kriging framework can be used to model uncertainty at predicted locations. The mean square prediction error (a.k.a. kriging variance, 0'3 ) is modeled by removing the portion of variance explained by spatial structure from the total variance 02. As the amount of variance explained by spatial structure approaches zero, the kriging variance approaches the total variance 02. This link can be expressed as: at = 02 —cT1u.-,o)C‘1c1u.,o) [Eq. 151 185 or, by substituting kriging weights (Eq. 13) for C ‘1c(u,-,0) , as: 2 2 T 00 =0 'C (111,0)40110) [HQ-16] An expected value Z(u,- )and a kriging variance 03 are sufficient for specifying a cdf at each location (if we also assume the distribution of values at each location is distributed as Normal). Furthermore, information provided by nearby random variables allows a ccdf to be specified. Therefore, alternate possible realizations can be generated at each location in the region of interest by drawing values randomly from each ccdf. Sequential Gaussian simulation (sGs) is a technique for generating alternate realizations of a random field that can be used as inputs to a spatial model. According to Goovaerts (1998:4), ”sequential simulation amounts to modelling the conditional cumulative distribution function (ccdf), ,then sampling it at each of the grid nodes visited along a random sequence. To ensure reproduction of the z-variogram model, each ccdf is made conditional not only to the original It data but also all values simulated at previously visited locations.” The promise of using not a single realization, but multiple alternate realizations of a random field is the ability to reveal uncertainty in results of a spatial model 186 given uncertainty about errors in data. When used within a Monte Carlo simulation framework, sequential Gaussian simulation can reveal uncertainty in results accurately and precisely (Goovaerts, 1998). Monte Carlo simulation is a computer-intensive framework that is useful for revealing contributions of input data errors to results, differences among uncertainties across a geographic region of interest (Burrough and McDonnell, 1998) and uncertainties that lurk in spatial model results. In essence, a typical simulation consists of evaluating a spatial model many times, each time using input data that have been modified randomly via a model that expresses an error distribution or makes assumptions about error behavior. Parameters of the model are derived with respect to known or assumed information about errors in the input data. After a set of 11 input realizations is generated, the spatial model is run 11 times and the set of n spatial model results is summarized to determine the range of possible outcomes at each location. Building the spatial model Schaetzl et al. (2002) built a spatial model for interpolating paleoshorelines along four glacial Lake Algonquin water planes; Figure 5.1 diagrams their workflow. In essence, a differentially upwarped water plane was modeled as a function of location using paleoshoreline position data and a second-order polynomial surface model with parameters estimated via ordinary least squares. The minimum measured paleoshoreline elevation value was subtracted from the 187 upwarped surface and the result represented differential warping relative to the lowest sampled feature. I shoreline e; :1; :~~:,: 1:5: .1 o 1 ‘ 1 . ‘ positions 1?; gala-M; .; :y l {X11112} ,7: '. fit elevation surface via ordinary least squares v 7 v ‘ T! T ... 5 t ‘ r o- ;g; .15.. .1. : iupWarpe‘fl . 7 water plane All “ min {2} 3 mm! subtract min {2} i. ' . . ,'.., .. . . , , . . 5' _ J. _ 1 0' " relatWE‘ ' ’ . ‘ 1 ‘ . ." v .,1- .. .1.1- 4 ..~ 3 ‘1 v . 1 Q - I 1— ‘ \. 31 ¢ . , .i 1 . 11 4 v . .1 ., .,1. 1 ; . . . ' I LLLLL I ' L l L subtract rebound from DEM I ] v1o1o-Iw K I 1 1 Ac . . , .1. . .1.’ ,4 ‘ ’ r ‘ '.‘ "r ‘ 1 ' L... \ 1 A 1 , ., 1. 1' 1 .‘, ...... . . A 1 ‘b' V ‘ ..................................., reclassify Figure 5.1: The spatial model Schaetzl et al. (2002) used to reconstruct landscapes in northern Michigan associated with glacial Lake Algonquin. 188 Next, the relative rebound surface was subtracted from a contemporary digital elevation model (DEM), thereby generating a topographic surface that was down-warped (i.e. differentially depressed) by the immense weight of the Laurentide ice sheet; this is similar to the method used by Krist and Schaetzl (2001). Finally, adjusted elevation values were classified with respect to the aforementioned lowest feature; DEM cells with elevation values greater than the lowest feature were classified above water and the balance was classified not above water. The binary classification scheme served to cartographically flood the depressed landscape, reveal the extent of the modeled water plane as constrained by topography and, hence, interpolate paleoshorelines at locations between measured paleoshoreline positions. Others have built similar models for reconstructing ancient water planes in Scotland (e. g., Smith, Cullingford and Firth, 2000; Firth, Smith and Cullingford, 1993) and interpolating paleoshorelines in Canada (e. g., Leverington and Teller, 2003; Leverington, Teller and Mann, 2002), but each assumed paleoshoreline data do not contain measurement errors (Poole and O’Farrell, 1971), which in my experience is not always the case. Consequently, their results are deterministic and do not convey any imprecision or uncertainty. I use a model similar to the one built by Schaetzl et al. (2002), but instead of generating a single deterministic result, I employ kriging methods and sequential Gaussian simulation (sGs) Within a Monte Carlo framework to take advantage of information I have about random errors in my paleoshoreline data 189 and digital elevation model, respectively. Figure 5.2 diagrams my workflow. The enhanced spatial model allows me to identify locations in my results that are subject to greater or lesser uncertainties. Knowing which locations suffer greater uncertainties is useful for qualifying my results (and prescribing limited uses) and for focusing subsequent fieldwork. Below I present the two major parts of the paleoshoreline interpolation model. The first part acknowledges possible errors that lurk in the paleoshoreline dataset and limitations associated with using ordinary least squares (OLS) to estimate upwarped water plane elevations. The second part acknowledges errors that lurk in the DEM and describes the models I use to generate alternate possible realizations of the contemporary ground surface. Ultimately, the results of these two parts are combined so that a comparison can be made at every DEM grid cell location between the elevation of the contemporary land surface and the elevation of the differentially-upwarped water plane. Unsampled locations with ground surface elevations equal to, or very near to, the co-located water plane elevation are identified as candidates for having unsampled paleoshorelines. 190 ahorelige N68 .1. positioii‘s )— benchmai‘ks {15.1.2} {£11112} __ J i reconstruct water plane Lvia universal kriging J calculate differences at upwarped "1111111111111 classify DEM cells with respect to the water plane classifications ,, {persistence ' ' as V . “359$". ; l L reclassify . 11111-111111 .: representation 91111191911191; l.— 7] control point locations 312:; simulate error surfaces via ordinary kriging and sGs : U100: 1‘99'13311103. _~ 16f DEM-error; I remove errors from DEM '1— 100 alternate f 1 , 2905.5.in . 1 1 . DEM}? Figure 5.2: The spatial model used in this work to digitally reconstruct water planes and landscapes in northern Michigan and western Ontario that were associated with glacial Lake Algonquin. The model was run 100 times for each water plane, which produced a set of 100 results that were summarized to reveal the influences of input data errors. 191 Models of elevation errors in paleoshoreline elevation data Errors in paleoshoreline elevation data are treated as the sum of random and independent and selection errors and measurement errors. I use a statistical model to express these errors and mitigate their influence. Assumed model of selection errors Uncertainty about selection errors in my paleoshoreline data can be addressed by adopting assumptions about selection error and employing a model based on those assumptions (Shortridge, 2001). Based on my experience in the field and via an analysis of the research of others (e. g., Lawson, 1891 ; Goldthwait, 1910; Morton and Speed, 1998), I expect most surveyed elevations are accurate proxies for former water levels. I expect, however, that some measured elevations might overestimate water levels because: a) the base of a wave-cut bluff typically forms at or slightly above water level; b) the bouldery/ gravelly lags that form at the bases of wave-cut bluffs vary in width (pers. obs.) according to local wave conditions and water levels fluctuations (Chapter 2); and c) I began my search for the bouldery/ gravelly lag at the highest point along each base. Accordingly, my search method preferentially selects the highest portion of each the bouldery/ gravelly lag. Therefore, there is a possibility that the highest measurable point on a particular bouldery/ gravelly lag might overestimate average water. 192 I assume most selection errors to be 0m because I expect most surveyed features represent water levels. Wherever errors do occur, I assume most will be small and fall within the interval of 0 and 1 m. I must recognize that, on rare occasions, an extreme error might exist and I assume it will fall within the interval of 1 and 3 m. I do not have reason to expect any selection error will exceed three meters. Selection error at any paleoshoreline location can be modeled as a random variable that is distributed as Exponential. I chose the Exponential distribution for its convenience: it can be adapted easily to describe the expectations and assumptions I have about selection errors (i.e., bounded by zero and decreasing probabilities with increasing error values). The Exponential probability density function is expressed as: f(x)=%e—(x_#)/fl x2#;/5’>0 [liq-17] where: x is a value in the distribution; )1 is the mode of the distribution; and [3 is the mean of the distribution. The specific selection error model employed here assumes an Exponential distribution with parameters )1 = 0.0 m and ,8 = 0.33 m. — (x) / 0.33 f (x) = 3e [Eq. 18] Selection errors are assumed to be spatially unstructured. A histogram (Figure 5.3a) shows one realization of 529 selection errors (one for each paleoshoreline 193 sample I have) that were generated using the Exponential model. Note the distribution is bounded by a value of zero and the frequency of a particular selection error occurring decreases with an increase in value. Empirical model of measurement errors Uncertainty in measurement error can be addressed by characterizing the measurement process and employing a model based on the characterization (Shortridge, 2001). In Chapter 3, I reported elevation errors measured at NGS control stations. From my assessment of these values, measurement error is modeled as a random variable that is distributed as Normal. The Normal probability density function is characterized as: e“(x‘#)2 /<2oz) f(x)= am [Eq. 19] where: x is a value in the distribution; u is the mean of the distribution; and o is the standard deviation. The specific measurement error model (Eq. 20) assumes a Normal distribution with parameters taken from the sample of control station errors: u = 0.16 m and o = 0.77 m. e—(x-0.16)2 /1.1858 f (x) = 1.9301 We 20] 194 A second histogram (Figure 5.3b) shows a realization of 529 measurement errors, which is the same number of ground control points that I have available, that were generated using the Normal model and parameters described above. 250 “i 250 A 200 ' A 200 l 5 5 5* 150- i 5: 150 C I?» C 3 1001- 3 100; ‘D S _'