A Gravity Study of Holocene Active Structures in the Puget Lowland of Washington State by James Paul Taylor A thesis submitted to the Graduate Faculty of Auburn University in partial fulfillment of the requirements for the Degree of Master of Science Auburn, Alabama May 5, 2013 Copyright 2013 by James Paul Taylor Approved by Lorraine W. Wolf, Chair, Professor of Geology and Geography Mark G. Steltenpohl, Professor of Geology and Geography James A. Saunders, Professor of Geology and Geography ii Abstract This study presents models based on existing and new gravity and magnetic data for two regions in the Puget Sound, Washington, area: the Bellingham basin and the Muckleshoot basin. The principal goals of the project are to identify basin boundaries and determine whether and how faults mapped outside of these basins continue beneath their thick sedimentary cover. In the Bellingham basin, new gravity maps reveal that the north basin boundary is located further south than previously thought. Cross-sectional models presented in this study focus on two previously identified structures, the Drayton-Harbor magnetic lineament and the proposed Birch Bay fault. Late Holocene displacements along these structures have been observed by other workers on the western Washington coast, and analyses of magnetic and LiDAR data suggest that these faults extend eastward into the basin. Models are consistent with the inland continuation of the Birch Bay fault towards the city of Bellingham and also suggest that a previously unidentified fault or fold could exist nearby. Results indicate that the Bellingham basin may consist of smaller subbasins, rather than one large basin, as previously mapped. In the Muckleshoot basin, gravity maps support the previously defined basin boundaries. Magnetic data and new gravity data support the existence of at least two subbasins within the Muckleshoot that are likely separated by a fault. A 2-D cross-sectional model focuses on the possible connection of the east-trending Tacoma fault (west of the basin) and the northwest-trending White River fault, mapped to the east of the Muckleshoot iii basin. A connection between these two fault systems would have significant implications for hazard estimates in terms of the length and size of these seismogenic structures and the maximum magnitudes that could be generated by faulting. iv Acknowledgments I would like to thank my faculty advisor, Dr. Lorraine Wolf, for her enormous support and guidance throughout this entire project. Without her motivation and encouragement the completion of this study would not have been possible. I would also like to thank my committee members Dr. Mark Steltenpohl and Dr. James Saunders for their patience, insight, comments, and overall contribution. I would like to especially thank Dr. Richard Blakely for supplying compilations of previously acquired gravity and magnetic data. His help and insight was critical to the completion of data processing. I would like to thank Dr. Brian Sherrod for supplying LiDAR digital elevation models as well as his recommendations for field mapping and problem solving. I would like to thank Dr. Michael Brown, Dr. Paul Bodin, and the University of Washington for loaning their gravimeter for my field component. I would like to thank Dr. Frank Danes and the University of Puget Sound for their assistance with establishing the Covington gravity base station. I would like to thank the United States Geological Survey for partially supporting this work through the National Earthquake Hazards Reduction Program through a grant to L. Wolf (USGS-G11AP20042). Thanks to the Geological Society of America for partially funding this project through a graduate research grant. The support of the Auburn University faculty and staff was felt through the duration of this task. I especially thank the faculty and friends that I have made in the Geology and Geography department. Thanks to Dr. Luke Marzen for allowing me to use v the Trimble XH GPS unit for establishing my new field stations. The constant support of my Auburn friends motivated me and helped me work towards this milestone in my life. Finally, I would like to thank my parents, family and friends for their sacrifice, patience, support, and encouragement throughout my life. I would not be where I am today without their love and selflessness. vi Table of Contents Abstract ............................................................................................................................... ii Acknowledgments ............................................................................................................. iv List of Tables ................................................................................................................... viii List of Figures .................................................................................................................... ix List of Abbreviations ....................................................................................................... xiv Introduction ..........................................................................................................................1 Geologic Background and Previous Work ..........................................................................3 Tectonic Setting and Regional Geology ..................................................................3 Bellingham basin .....................................................................................................6 Muckleshoot basin .................................................................................................16 Methodology ......................................................................................................................23 Gravity Data ...........................................................................................................23 Data Collection ......................................................................................................23 Magnetic Data ........................................................................................................27 Data Processing ......................................................................................................27 Data Modeling .......................................................................................................29 Cross-sectional Profiles .........................................................................................35 Results ................................................................................................................................38 Bellingham .............................................................................................................38 vii Bellingham Gravity Grid Maps .......................................................................38 Bellingham Magnetic Grid Maps.....................................................................41 Bellingham Cross-sections...............................................................................44 Gravity and Magnetic Models: A-A? Profile .......................................48 Gravity and Magnetic Model: B-B? Profile .........................................51 Muckleshoot ...........................................................................................................53 Muckleshoot Gravity Grid Maps .....................................................................53 Muckleshoot Magnetic Grid Maps ..................................................................57 Muckleshoot Cross-section ..............................................................................61 Gravity and Magnetic Model: C-C? Profile .........................................61 Discussion ..........................................................................................................................65 Bellingham basin ...................................................................................................65 Birch Bay fault .....................................................................................68 Drayton Harbor fault ............................................................................69 Bellingham Bay feature .......................................................................69 Muckleshoot basin .................................................................................................70 Disclaimer ..............................................................................................................74 Conclusions ........................................................................................................................76 References ..........................................................................................................................79 Appendix ............................................................................................................................86 viii List of Tables Table 1. Stratigraphy of the Bellingham basin (Brown and Dragovich, 2003; Brown and Gehrels, 2007) ..........................................................................................................7 Table 2. Stratigraphy of the Muckleshoot Basin (Johnson et al., 2004; Jones, 1996; Tre?hu et al., 1994; Finn, 1990; Lees and Crosson, 1990; Yount et al., 1985) ......17 Table 3. Blocks/polygons and their depths for A-A?, with densities, susceptibilities, and lithologies listed for each body (Brown and Dragovich, 2003; Brown and Gehrels, 2007) ......................................................................................................................50 Table 4. Blocks/polygons and their depths for B-B?, with densities, susceptibilities and lithologies listed for each body (Brown and Dragovich, 2003; Brown and Gehrels, 2007) ......................................................................................................................52 Table 5. Blocks/polygons and their depths for C-C?, with densities, susceptibilities and lithologies listed for each body (Johnson et al., 2004; Jones, 1996; Tre?hu et al., 1994; Finn, 1990; Lees and Crosson, 1990; Yount et al., 1985) ...........................64 ix List of Figures Figure 1. Regional tectonic map from Kelsey et al. (2012) modified from Wells at al. (1998). Arrows show motion of the Juan de Fuca plate relative to the North American plate (after McCaffrey et al., 2007). The Oregon coast block is moving northward relative to the North American plate toward the deforming Northern Cascadia forearc .......................................................................................................4 Figure 2. Generalized tectonic map of the Puget Lowland and surrounding regions modified from Sherrod et al. (2008). The Muckleshoot study area is shown with a green box. The Bellingham Study area is shown with a red box. Abbreviations: BB=Bellingham basin, EB=Everett basin, SB=Seattle basin, SU=Seattle uplift, TB= Tacoma basin, BCF=Boulder Creek fault, SF=Seattle fault, TF=Tacoma fault, WRF=White River fault, CRBF=Coast Range boundary fault, RMF=Rattlesnake Mountain fault, B=Bellingham, E=Everett, O=Olympia, S=Seattle, T=Tacoma. Redrawn from Brocher et al. (2001). Other abbreviations represent structures as listed in Sherrod et al. (2008) ..............................................5 Figure 3. Geologic map of the Bellingham study area modified from Brown and Dragovich (2003). Line B-B? shows the location of a geologic cross-section. Solid black lines indicate contacts (depositional or intrusive). Dotted lines indicate concealed faults. Magenta lines with arrows represent anticlines (arrows away) and synclines (arrows towards). Other abbreviations represent units and structures as listed in Brown and Dragovich (2003) ................................................8 Figure 4. Geologic cross-section (west to east) of line B-B? (Figure 3) modified from Brown and Dragovich (2003). The yellow region, bounded by vertical arrows, represents the location of the Bellingham study area. Abbreviations: CN=Chuckanut Formation, EA=Easton Metamorphic Suite, OC=Orcas Chert, BP=Bell Pass M?lange, ES=East Sound Group, CN=Chilliwack Group, NK=Nooksack Formation. Other abbreviations represent units and structures as listed in Brown and Dragovich (2003).....................................................................9 Figure 5. (a) Map from Kelsey et al. (2012) of northwestern Washington and southern British Columbia showing the western portion of the Bellingham study area. Line X-X? marks the location of the south to north seismic line shown in Figure 5b. (b) Seismic reflection profile X-X? reproduced by Kelsey et al. (2012) from Hurst (1991). This south-to-north seismic reflection profile shows multiple faults and the Birch Bay anticline. This on-land seismic reflection line was shot by American Hunter Exploration, Limited .................................................................11 x Figure 6. Map of Bellingham study area modified from Kelsey et al. (2012), showing major faults, scarps, anomalies, and interpreted uplift, along with geologic units. Line X-X? marks the location of the south to north seismic line shown in Figure 5b (Hurst, 1991). Solid red lines show probable faults .........................................13 Figure 7. LiDAR DEM (acquired summer 2006 by U.S. Geological Survey) from Kelsey et al. (2012) showing the western portion of the Bellingham study area and three proposed Holocene faults (solid red lines). Anticlines axes are marked by red dashed lines. Line X-X? shows the location of the seismic reflection profile reproduced in Figure 5b .........................................................................................14 Figure 8. (a) LiDAR image showing a north-west trending scarp near Drayton Harbor in the Bellingham study area. Both white boxes are bisected diagonally by a scarp (Kelsey et al., 2012; B. Sherrod, pers. comm., 2011); (b) 2-km south-side-up scarp southeast of Blaine, Washington; (c) Two south-side-up scarps located southeast of Figure 8b ............................................................................................15 Figure 9. (a) Isostatic residual gravity map including the majority of the Muckleshoot study area (black box) modified from Johnson et al. (2004). (b) Magnetic anomaly map including the Muckleshoot study area. Abbreviations: MB=Muckleshoot basin, PS=Puget Sound, SFZ=Seattle fault zone, SU=Seattle uplift, T=Tacoma, TB=Tacoma basin. Other abbreviations represent features as listed in Johnson et al. (2004) ................................................................................................................18 Figure 10. Seismic tomographic image (top) from 2.5 km depth, isostatic residual gravity map (middle), and magnetic anomaly map (bottom) for the Tacoma and Muckleshoot basins (Liberty, 2007). Red lines are Boise State high-resolution land-based seismic reflection profiles. Lines A-A?, B-B?, and C-C? coincide with profiles in Figure 11. WRF= White River fault. Basins are outlined with dashed black lines. Other lines represent features as listed in Liberty (2007) ...................21 Figure 11. Cross-sections of the Muckleshoot basin from profiles in Figure 10 (Liberty, 2007). (a) A-A? is a west to east tomographic cross-section; (b) B-B? is a south to north cross-section that shows seismic data, borehole data, and estimated Quaternary strata depth; (c) C-C? is a south to north cross-section that combines seismic data, borehole data, and velocity data .......................................................22 Figure 12. Locations of gravity sites/stations in the Bellingham study area. The top map shows the existing measurements before field work began (red dots) superimposed on an area made up of 15 (1:24000) USGS topographic maps. The bottom map shows both the existing locations and new locations (blue dots) acquired during this investigation ..........................................................................25 Figure 13. Locations of gravity sites/stations in the Muckleshoot study area. The top map shows the existing measurements before field work began (red dots) xi superimposed on an area made up of 15 (1:24000) USGS topographic maps. The bottom map shows both the existing locations and new locations (blue dots) acquired during this investigation ..........................................................................26 Figure 14. Illustration of how a spherical object produces (a) a slanted anomaly before filtering, and (b) a symmetric anomaly after applying reduction to pole filter (modified by Bajgain, 2011; after Blakely, 1995) .................................................29 Figure 15. Diagram showing how the minimum curvature gridding process achieves the smoothest possible surface. Iterations continue until the observed data are closest to the grid nodes (green dots) (Geosoft: Topics in Gridding) ................................30 Figure 16. Two minimum curvature grid maps showing the effects of a vertical derivative filter. (a) Top map shows total field anomaly (TFA). (b) Bottom map enhances shallow features as a result of applying the vertical derivative filter ....................31 Figure 17. Grid maps showing the effects of poorly chosen grid size. (a) Top map shows vertical rippling in TFA due to poorly chosen grid size. (b) Bottom map shows how a vertical derivative filter exaggerates the effects of an inappropriate choice of grid size. White stripes on the right of both maps result from spacing changes in sampling protocols combined with incorrect grid size ......................................33 Figure 18. (a) TFA map including an artifact caused by a nearby smelting plant. (b) Map showing a vertical derivative filter of Figure 18a. Filtering creates an unwanted bullseye artifact. (c) TFA map where anomaly has been masked. (d) Map showing a vertical derivative filter of Figure 18c. Masking removes the bullseye artifact .35 Figure 19. Complete Bouguer anomaly (CBA) gravity map of the Bellingham study area. Black dots=station locations. A-A? and B-B?=profile locations. Solid blue lines=faults. Dashed blue lines=inferred faults. Abbreviations: BB=Birch Bay, BH=Bellingham, BHB=Bellingham Bay, BL=Blaine, DH=Drayton Harbor .......39 Figure 20. Gravity (CBA) map of the Bellingham study area with vertical derivative filter applied. A-A? and B-B?=profile locations. Solid blue lines=faults. Dashed blue lines=inferred faults. White arrows=previously unseen high anomalies. Abbreviations: BB=Birch Bay, BH=Bellingham, BHB=Bellingham Bay, BL=Blaine, DH=Drayton Harbor, DHF=Drayton Harbor fault, BBF=Birch Bay fault ........................................................................................................................40 Figure 21. Magnetic (TFA) map of the Bellingham study area. Black dots=station locations. A-A? and B-B?=profile locations. Solid blue lines=faults. Dashed blue lines=inferred faults. Abbreviations: BB=Birch Bay, BH=Bellingham, BHB=Bellingham Bay, BL=Blaine, DH=Drayton Harbor ....................................42 Figure 22. Magnetic (TFA) map of the Bellingham study area with vertical derivative filter applied. A-A? and B-B?=profile locations. Solid blue lines=faults. Dashed xii blue lines=inferred faults. Abbreviations: BB=Birch Bay, BH=Bellingham, BHB=Bellingham Bay, BL=Blaine, DH=Drayton Harbor ....................................43 Figure 23. Magnetic ternary map of the Bellingham study area generated by combining the grids from Figures 21 and 22 with an analytic signal grid. A-A? and B- B?=profile locations. Solid blue lines=faults. Dashed white lines=inferred faults. Abbreviations: BB=Birch Bay, BH=Bellingham, BHB=Bellingham Bay, BL=Blaine, DH=Drayton Harbor, T= Total field, DZ=vertical derivative, AS=analytic signal .................................................................................................45 Figure 24. Gravity and magnetic profiles with geologic cross-section for A-A?. Abbreviations: QS= Quaternary Sediments, CN= Chuckanut Formation, EA= Easton Metamorphic Suite, BP= Bell Pass M?lange, CH= Chiliwack Group, NK= Nooksack Formation, BBF=Birch Bay Fault, DHF=Drayton Harbor Fault ..........46 Figure 25. Gravity and magnetic profiles with geologic cross-section for B-B?. Abbreviations: QS= Quaternary Sediments, CN= Chuckanut Formation, EA= Easton Metamorphic Suite, BP= Bell Pass M?lange, CH= Chiliwack Group, NK= Nooksack Formation, BBF=Birch Bay Fault ........................................................47 Figure 26. Gravity (CBA) map of the Muckleshoot study area. Black dots=station locations. C-C?=profile location. Solid blue lines=faults. Abbreviations: CV=Covington, K=Kent, PS=Puget Sound, R=Renton, T=Tacoma, WRF=White River fault ..............................................................................................................54 Figure 27. Gravity (CBA) map of the Muckleshoot study area with vertical derivative filter applied. C-C?=profile location. Solid blue lines=faults. Arrows=aligned gravity highs. Abbreviations: CV=Covington, K=Kent, PS=Puget Sound, R=Renton, T=Tacoma, WRF=White River fault...................................................56 Figure 28. Magnetic (TFA) map of the Muckleshoot study area. C-C?=profile location. Solid blue lines=faults. Abbreviations: CV=Covington, K=Kent, PS=Puget Sound, R=Renton, T=Tacoma, WRF=White River fault ......................................58 Figure 29. Magnetic (TFA) map of the Muckleshoot study area with vertical derivative filter applied. C-C?=profile location. Solid blue lines=faults. Abbreviations: CV=Covington, K=Kent, PS=Puget Sound, R=Renton, T=Tacoma, WRF=White River fault ..............................................................................................................59 Figure 30. Three grid composite of the Muckleshoot study area generated by combining the grids from Figures 28 and 29 with an analytic signal grid. C-C?=profile location. Solid blue lines=faults. White arrows= highlighted sharp anomaly. Abbreviations: CV=Covington, K=Kent, PS=Puget Sound, R=Renton, T=Tacoma, WRF=White River fault .....................................................................60 xiii Figure 31. Gravity and magnetic profiles with geologic cross-section for C-C?. Abbreviations: G1=Glacial and post glacial strata 1, G2=Glacial and post glacial strata 2, OP=Ohanapecosh Formation, P1= Puget Group 1, P2= Puget Group 2, CR= Crescent Formation, TF/WRF=Tacoma fault/ White River fault, vertical purple polygon=igneous dike .................................................................................63 Figure 32. Vertical derivative maps of gravity (a) and magnetic (b) data for the Bellingham basin with locations of structures and boundaries as determined by this study. White dashed line=basin boundary. Black dotted line=possible unidentified northwest striking structure. Red dashed lines=location of faults proposed by Kelsey et al. (2012). (a) Map showing gravity. (b) Map showing magnetics. DHF=Drayton Harbor fault. BBF=Birch Bay fault .............................67 Figure 33. Maps with vertical derivative filter, showing locations of structures and boundaries as determined by this study for the Muckleshoot basin. White dashed line=basin boundary. Black dotted line=possible unidentified northwest striking structures. (a) Map showing gravity. (b) Map showing magnetics. TF/WRF=Tacoma fault/White River fault seen in cross-section (Figure 31) .......73 xiv List of Abbreviations LiDAR Light detection and ranging DEM Digital elevation model GPS Global positioning system TFA Total field anomaly CBA Complete Bouguer anomaly DZ Vertical derivative AS Analytic signal CMY Cyan-magenta-yellow mGal Miligal nT nanotesla NAD North American Datum UTM Universal Transverse Mercator QS Quaternary glacial sediments CN Chuckanut Formation EA Easton Metamorphic Suite BP Bell Pass M?lange CH Chiliwack Group NK Nooksack Formation BB Birch Bay xv BBF Birch Bay fault DHF Drayton Harbor fault TF Tacoma fault WRF White River fault 1 Introduction Knowing where faults are and how they behave is critical information for evaluating earthquake hazard. The region surrounding Seattle, Washington, is situated in an active subduction zone, which has created a complex network of faults. Subduction zones have the potential to generate the largest earthquakes, yet the last major earthquake affecting the area occured around 1700 A.D. (Atwater, 2003). The Puget Sound region is home to millions of people and major cities that could potentially be affected in a major earthquake event, making it important to gain an understanding of the nature of these faults and how they interconnect. Recent paleoseismic studies, light detection and ranging (LiDAR) data, potential field studies, and geologic mapping have uncovered new evidence of seismic activity associated with some faults in Holocene time. What these studies on individual fault segments imply or how the larger fault systems interact remains unknown, particularly because many of the faults are buried or otherwise obscured by sediments deposited in this forearc basin formed in response to plate subduction. This project focuses on structures that have been previously recognized in two main regions of the Puget Sound area: the Bellingham basin and the Muckleshoot basin. It uses detailed gravity data to study the extent and relation of recently discovered faults and other structures in and surrounding both basins. The project (1) increases our understanding of fault locations and orientations in the site areas, (2) yields new data that 2 could possibly constrain fault models of the region, and (3) provides insight into the behavior of other fault systems and basins nearby. These results help us understand the regional earthquake hazard, which is critical for estimating ground motions of future earthquakes. 3 Geologic Background and Previous Work Tectonic Setting and Regional Geology For such a dynamic plate tectonic environment and locality, surprisingly little detail is known about the subsurface of Western Washington. This area is located at an active plate margin between the North American and Juan de Fuca plates (Figure 1). Continuous seafloor spreading drives the Juan de Fuca plate toward the North American plate, where it is subsequently subducted. This subduction zone boundary extends along the west coast from northern California to British Columbia and is the driving force behind the actively deforming Cascadian forearc (Kelsey et al., 2012) (Figure 1). The Sierra Nevada block is forced northward against the fixed Canadian buttress (Wells et al., 1998; McCaffrey et al., 2000, 2007). Stress in the Pacific Northwest region is characterized by north-south to north- northeast-south-southwest compression. These long-term stresses have shaped western Washington since the Tertiary (Zoback, 1992; Wells et al., 1998). This sandwiching causes shortening of approximately 4.5 mm/yr in the forearc. Networks of faults and folds form in the northern extents of the Cascadian forearc in response to variable distribution of the regional north-south strain (Mazzotti et al., 2002; McCaffrey et al., 2007) (Figure 2). Features formed in this dynamic environment are exposed at the surface in some localities and can be easily mapped. In the forearc basins, however, basement 4 rock units that could potentially contain evidence of fracture and folding are buried under younger sediments, making it necessary to use different methods to map those structures. Figure 1. Regional tectonic map from Kelsey et al. (2012) modified from Wells at al. (1998). Arrows show motion of the Juan de Fuca plate relative to the North American plate (after McCaffrey et al., 2007). The Oregon coast block is moving northward relative to the North American plate toward the deforming Northern Cascadia forearc. 5 Figure 2. Generalized tectonic map of the Puget Lowland and surrounding regions modified from Sherrod et al. (2008). The Muckleshoot study area is shown with a green box. The Bellingham Study area is shown with a red box. Abbreviations: BB=Bellingham basin, EB=Everett basin, SB=Seattle basin, SU=Seattle uplift, TB= Tacoma basin, BCF=Boulder Creek fault, SF=Seattle fault, TF=Tacoma fault, WRF=White River fault, CRBF=Coast Range boundary fault, RMF=Rattlesnake Mountain fault, B=Bellingham, E=Everett, O=Olympia, S=Seattle, T=Tacoma. Redrawn from Brocher et al. (2001). Other abbreviations represent structures as listed in Sherrod et al. (2008). 6 Bellingham basin The Bellingham basin is a sedimentary basin located in the extreme northwest corner of Washington named after the town of Bellingham (Figure 2). The study area is bounded to the north by Canada and to the west by the Puget Sound (Figures 1 and 2). The study area centers on key structures and does not encompass the entire basin. The Bellingham basin is the northernmost basin in the forearc and has existed since the mid- Tertiary. The precise basin extent is unknown but it has been suggested to extend to the west into the Puget Sound and at least 15 km northward into Canada. The Eastern extent is believed to be the same as the modern broader Georgia basin (Gallup, 1957). No faults or folds have been definitively mapped at Earth?s surface in the sedimentary basin itself, although there is a large amount of evidence supporting the presence of these structures, especially in the Birch Bay area. Any structure present is easily hidden by Pleistocene and younger sediments that blanket the basin (Easterbrook, 1963; Easterbrook, 1976). The stratigraphy of the Bellingham basin can be subdivided into six main layers (Table 1). Several units described here have been classified differently by previous workers. The youngest layer is composed of unconsolidated Quaternary glacial deposits. This overlies the Chuckanut Formation, Easton Metamorphic Suite, Bell Pass M?lange, Chilliwack Group, and the Nooksack Formation. These units are very diverse though most are of marine origin. Brown and Dragovich (2003) published a state map with descriptions of these units (Figure 3). This map contained several cross-sections, one of which extends very close to the Bellingham study area. It is labeled as line B-B? (Figure 4). 7 Table 1. Stratigraphy of the Bellingham basin (Brown and Dragovich, 2003; Brown and Gehrels, 2007). Unit Time Interval Tectonic Setting Protolith Glacial Sediments Quaternary Sedimentary basin, fluvial Previously deposited units of all types (Igneous, sedimentary, and metamorphic). Chuckanut Formation Tertiary (Eocene- Oligocene) Intracontinental rift zones, likely strike- slip pull-apart basins Arkosic and lithic sandstone, siltstone, conglomerate, and coal. Mostly fluvial deposition with some marine. Easton Metamorphic Suite Jurassic Ocean ridge, trench, transitional area to island arc Metabasite, siliceous carbonaceous phyllite, metagraywacke, metabasalt, gabbro-tonalite plutons, and metatuffs. Bell Pass M?lange/ Orcas Chert Triassic- Middle Jurassic Ocean floor and trench with some exotic accreted rocks Ribbon chert, mudstone, pillow basalt, limestone, volcanic lithic clastic rocks. Chilliwack Group/ East Sound Group Devonian- Permian Island arc Volcanic lithic sandstones and siltstones, basaltic to dacitic volcanic rocks, fossiliferous to reefoidal limestones. Andesitic to dacitic rocks, with lesser limestone lenses, shale, argillite, graywacke, and conglomerate. Nooksack Formation Middle Jurassic- Early Cretaceous Island arc Marine volcanic-lithic sandstone, siltstone, and local conglomerate. Local andesitic-dacitic volcanic member. 8 Figure 3. Geologic map of the Bellingham study area modified from Brown and Dragovich (2003). Line B-B? shows the location of a geologic cross-section. Solid black lines indicate contacts (depositional or intrusive). Dotted lines indicate concealed faults. Magenta lines with arrows represent anticlines (arrows away) and synclines (arrows towards). Other abbreviations represent units and structures as listed in Brown and Dragovich (2003). 9 Figure 4. Geologic cross-section (west to east) of line B-B? (Figure 3) modified from Brown and Dragovich (2003). The yellow region, bounded by vertical arrows, represents the location of the Bellingham study area. Abbreviations: CN=Chuckanut Formation, EA=Easton Metamorphic Suite, OC=Orcas Chert, BP=Bell Pass M?lange, ES=East Sound Group, CN=Chilliwack Group, NK=Nooksack Formation. Other abbreviations represent units and structures as listed in Brown and Dragovich (2003). 10 The Brown and Dragovich (2003) B-B? cross -section shows large faults to the east and west of the study area and some folding throughout. In the basin, the units are approximately horizontal with some notable folding occurring to the west. Further evidence of these structures can be seen in a seismic line shot by American Hunter Exploration, Limited. Unfortunately, these data are only available in the form of a reflection profile published by Hurst (1991) (Figure 5). In this profile, a well-defined antiform is shown as the Birch Bay anticline. A network of multiple meter-scale faults is also visible. These structures occur at shallow depths and suggest that the Bellingham basin has experienced northward shortening in the Quaternary (Kelsey et al., 2012). Interestingly, magnetic anomalies are present in the area of the seismic line. Kelsey et al. (2012) filtered the magnetic data in order to highlight and emphasize the shallow magnetic sources. 11 Figure 5. (a) Map from Kelsey et al. (2012) of northwestern Washington and southern British Columbia showing the western portion of the Bellingham study area. Line X-X? marks the location of the south to north seismic line shown in Figure 5b. (b) Seismic reflection profile X-X? reproduced by Kelsey et al. (2012) from Hurst (1991). This south-to-north seismic reflection profile shows multiple faults and the Birch Bay anticline. This on-land seismic reflection line was shot by American Hunter Exploration, Limited. 12 Using an upward continuation, Kelsey et al. (2012) found that a network of linear magnetic anomalies exist in the study area (Figure 6). One of these magnetic lineaments intersects the seismic profile shown in Figure 5. This particular lineament trends northwest and can be traced westward to the shoreline of Birch Bay. Kelsey et al. (2012) found evidence of Holocene motion in the immediate Birch Bay area through careful examination of LiDAR digital elevation models (DEMs) (Figure 7) and subsequent shallow core sampling. By filtering the LiDAR DEMs, vegetation was removed. This allowed uplifted shoreline and beach berms in the area to be identified. The northern Puget Sound region is thought to have experienced gradual uplift since the late Holocene (Clague and James, 2002; James et al., 2009). Detrital wood and seed samples were picked from the cores for radiocarbon dating. Dates returned from the organic samples suggest that this regional trend is not uniform locally. Kelsey et al. (2012) show that their northern samples were taken from areas that experienced an abrupt uplift, while samples taken from the southern location showed evidence of abrupt subsidence. From these results, Kelsey et al. (2012) postulate that the trace of the Birch Bay fault lies somewhere between these two field areas. Other LiDAR images show multiple scarps trending west-northwest through the study area, which are coincident with magnetic lineaments formed by high magnetic gradients (B. Sherrod, pers. comm., 2010). One major scarp can be seen south of the town of Blaine and seems to parallel the Drayton Harbor magnetic lineament (Figure 8). The scarp terminates in the basin, but the magnetic anomaly continues eastward towards the town of Sumas. Similar magnetic lineaments have been associated with scarps directly related to Holocene fault activity. Near Sumas, the Boulder Creek, Sumas, and 13 Figure 6. Map of Bellingham study area modified from Kelsey et al. (2012), showing major faults, scarps, anomalies, and interpreted uplift, along with geologic units. Line X-X? marks the location of the south to north seismic line shown in Figure 5b (Hurst, 1991). Solid red lines show probable faults. 14 Figure 7. LiDAR DEM (acquired summer 2006 by U.S. Geological Survey) from Kelsey et al. (2012) showing the western portion of the Bellingham study area and three proposed Holocene faults (solid red lines). Anticlines axes are marked by red dashed lines. Line X-X? shows the location of the seismic reflection profile reproduced in Figure 5b. 15 Figure 8. (a) LiDAR image showing a north-west trending scarp near Drayton Harbor in the Bellingham study area. Both white boxes are bisected diagonally by a scarp (Kelsey et al., 2012; B. Sherrod, pers. comm., 2011); (b) 2-km south-side-up scarp southeast of Blaine, Washington; (c) Two south-side-up scarps located southeast of Figure 8b. 16 Vedder Mountain faults trend west-southwest into the basin (B. Sherrod, pers. comm., 2010). Some of these faults coincide with magnetic lineaments that extend westward in to the basin. Muckleshoot basin The Muckleshoot basin is located south of Seattle and east of Tacoma (Figure 2). It was first described by vanWagoner et al. (2002) as being 7 to 9 km deep. No faults or folds have been definitely mapped in the sedimentary basin itself, though there is evidence that large structures identified on the basin boundaries possibly extend beneath it. Any structure present could easily be hidden by the Pleistocene and younger sediments that cover the basin. Most units described in the basin and surrounding areas have been given multiple names, depending on the publication source and the location of the unit. However, the stratigraphy of this area can generally be described as consisting of four main units. The glacial and post-glacial strata constitute the youngest layer in the basin. This overlies the Puget Group and the Crescent Formation. The Ohanapecosh Formation partially enters the study area in the east. It is a volcaniclastic flood basalt deposit that is not laterally continuous in the area. All of these layers contain multiple rock types that have varying properties. For this reason the glacial and post-glacial strata and the Puget Group have been subdivided (Table 2). Potential field maps (Finn, 1991; Blakely, 1999) along the western portions of the study area reveal anomalies associated with basement structures (Figure 9). These maps 17 Table 2. Stratigraphy of the Muckleshoot Basin (Johnson et al., 2004; Jones, 1996; Tre?hu et al., 1994; Finn, 1990; Lees and Crosson, 1990; Yount et al., 1985). Unit Period Time Interval Tectonic Setting Protolith Glacial and Post Glacial Strata 1 Qua ter na ry Pleistocene-Holocene Sedimentary basin, fluvial Clay, silt, and sand derived from rivers and erosion of Pleistocene deposits. 2 Pliocene- Pleistocene- Sedimentary basin. Rivers, streams, glaciation Clay, silt, and sand derived from rivers and erosion of Pleistocene deposits. transitional and interglacial deposits Te rtiar y Ohanapecosh Formation Oligocene Flood basalts Volcaniclastic deposits and flood basalts Puget Group 1 Oligocene Non-marine to marginal marine Volcanics and non-marine sedimentary rocks. Correlatives to the Blakely Formation (Oligocene) 2 Eocene Non-marine to marginal marine Volcanics and non-marine sedimentary rocks. Correlatives to the Tukwilia Formation (Eocene) Crescent Formation Eocene Marine Marine basaltic rocks that are dense and magnetic. Sedimentary rocks are inter- bedded at the top of the formation. 18 Figure 9. (a) Isostatic residual gravity map including the majority of the Muckleshoot study area (black box) modified from Johnson et al. (2004). (b) Magnetic anomaly map including the Muckleshoot study area. Abbreviations: MB=Muckleshoot basin, PS=Puget Sound, SFZ=Seattle fault zone, SU=Seattle uplift, T=Tacoma, TB=Tacoma basin. Other abbreviations represent features as listed in Johnson et al. (2004). 19 show gravity and magnetic highs occurring in areas experiencing local uplift near sedimentary basins (Johnson et al., 2004). Gravity data clearly reveal the Seattle uplift to the north and the Tacoma basin to the south. The Muckleshoot basin can be seen along the eastern edge of the map and appears to be made up of two depressions (Figure 9a). This subdivision could be caused by an extension of the White River and/or Tacoma faults into the basin. Although there are various other deformational structures present in the area, the Tacoma fault system to the west and the White River fault system to the east are of particular interest (Figure 2). The White River fault was interpreted by Tabor et al. (2000) as being a normal fault with a subvertical dip to the southwest. Newer data suggest that it is a reverse fault that has a steep dip and verges to the south (Box et al., 2003). It is thought to extend from the Puget Lowland to east of the Cascades. The Tacoma fault is an active south-verging reverse fault that extends from the western edge of the Muckleshoot basin westward through the Puget Sound (Blakely et al., 2007) (Figure 2). It has been clearly defined in offshore seismic reflection profiles (Pratt et al., 1997; Brocher et al., 2001; Johnson et al., 2004), LiDAR data, paleoseismic data (Haugerud et al., 2003), and aeromagnetic data (Blakely et al., 2002). The Tacoma fault is oriented on strike with the western extent of the White River fault. Some evidence supports activity on the Tacoma fault during Quaternary (Sherrod et al., 2004). Aeromagnetic data show a west-northwest trending lineament that appears to follow the Tacoma fault; however, this anomaly continues eastward into the Muckleshoot basin, where the Tacoma fault disappears (Blakely et al., 2007). 20 If these two faults are joined beneath the basin, then the earthquake hazard in the region could be significantly higher than previously thought. Hazard increases with fault size because bigger faults can potentially have larger rupture areas. Liberty (2007) modeled three transects through the Muckleshoot basin (Figure 10). The locations of these profiles have been plotted over a seismic tomography image of 2.5 km depth, isostatic residual gravity map, and a magnetic anomaly map (Liberty, 2007). The Muckleshoot basin is outlined in the eastern portion of each map and is characterized by velocity, gravity, and magnetic lows. An apparent subdivision of the basin is evident in the gravity map. Plotting the location of the White River fault reveals that the basin anomaly has a similar trend. Profile A-A? is a west-east tomographic cross-section that shows the seismic velocities of basin rocks to a depth of 10 km (Figure 11a). The B-B? cross-section extends from south to north and incorporates three seismic lines in conjunction with borehole data (Figure 11b). This cross-section is 2 km deep and shows the thickness of Quaternary fill. C-C? is also a south-north cross-section through the Muckleshoot basin (Figure 11c). It combines a tomographic velocity profile, a published industry seismic line, and borehole data to constrain the depth of Quaternary sediments, along with various zones of the Tertiary Puget Group, down to a depth of 5 km (Liberty, 2007). 21 Figure 10. Seismic tomographic image (top) from 2.5 km depth, isostatic residual gravity map (middle), and magnetic anomaly map (bottom) for the Tacoma and Muckleshoot basins (Liberty, 2007). Red lines are high-resolution land-based seismic reflection profiles. Lines A-A?, B-B?, and C-C? coincide with profiles in Figure 11. WRF= White River fault. Basins are outlined with dashed black lines. Other lines represent features as listed in Liberty (2007). 22 Figure 11. Cross-sections of the Muckleshoot basin from profiles in Figure 10 (Liberty, 2007). (a) A-A? is a west to east tomographic cross-section; (b) B-B? is a south to north cross-section that shows seismic data, borehole data, and estimated Quaternary strata depth; (c) C-C? is a south to north cross-section that combines seismic data, borehole data, and velocity data. 23 Methodology Gravity Data Gravity contrasts can be directly related to changes in density. High or low Bouguer gravity anomalies occur when there are changes in the crustal rock density from one location to another (Phillips et al., 1993). For this reason, high gravity gradients can be used to determine contacts and boundaries between rock units. In this study, gravity will mainly be utilized to define lithologic contacts and structural features, such as basin boundaries. For the gravity study, all new stations were established on land and values were measured from the ground surface. These raw measurements were processed with a standard series of corrections to adjust for variability in Earth parameters, such as shape, rotation, topography, etc., and instrument drift. This chapter describes the details involved in the field component data collection and subsequent data processing for this study. New gravity data used in this study were collected in the summer of 2011. Existing gravity and magnetic data, from previous campaigns, were also obtained, along with some well-log data. Data were then processed and analyzed. Data Collection In the summer of 2011, 320 new gravity measurements were acquired and added to the existing gravity database (Dater, 1999; R. Blakely, pers. comm., 2011) for areas of 24 interest to this study (Figures 12 and 13). Gravity measurements from previous investigations were mainly located at road intersections. Intersection locations are desirable because they are easy to locate, often have safer places to make measurements, and, in some areas, are set up on a one-mile grid. New measurements locations followed this same strategy and were chosen to achieve minimal coverage of approximately one measurement per square mile. Several transects of half-mile intervals were also taken across target areas in order to get data for detailed modeling of smaller features. The latitude, longitude, and elevation of each location were recorded using a Trimble? GeoXH Global Positioning System (GPS) unit, with an initial horizontal positioning accuracy ranging from 6 to 9 m. Local GPS base stations were used to help adjust for environmental conditions present at the time of each measurement during post- processing. Post-processing of the data using Trimble? GPS Pathfinder? Office increased horizontal accuracy to less than 20 cm and vertical accuracy to less than 40 cm. Gravity values were obtained with a LaCoste and Romberg model G geodetic gravity meter, borrowed from the University of Washington?s Earth and Space Sciences department. Gravity base stations established in the 1960s and 1970s were used to merge new data with the existing database, as well as to correct for unavoidable instrument drift over the collection period. The first and last measurement each day was taken at a base station, following a standard looping procedure. The Muckleshoot site did not have any feasible base stations closeby, so a local base station was established in Covington, Washington, from an established base station at the University of Puget Sound through a repeated looping technique. The Covington base station was located at the main entrance of the City of Covington building (47.3590427322?, -122.118273863?) (Appendix). 25 Figure 12. Locations of gravity sites/stations in the Bellingham study area. The top map shows the existing measurements before field work began (red dots) superimposed on an area made up of 15 (1:24000) USGS topographic maps. The bottom map shows both the existing locations and new locations (blue dots) acquired during this investigation. 26 Figure 13. Locations of gravity sites/stations in the Muckleshoot study area. The top map shows the existing measurements before field work began (red dots) superimposed on an area made up of 15 (1:24000) USGS topographic maps. The bottom map shows both the existing locations and new locations (blue dots) acquired during this investigation. 27 Magnetic Data Magnetic minerals are usually associated with igneous and metamorphic rocks. Detrital sedimentary rocks are generally nonmagnetic but if magnetic mineral grains were present in their source rocks, they can be magnetic. In western Washington, forearc basins have been filled with sediments derived from a variety of source rocks. Strongly magnetic igneous and metamorphic rocks are present nearby (Figure 6). They undoubtedly have been weathered and incorporated into the basin fill over time. In some circumstances, fault zones can be correlated to areas of high magnetic gradient. Data Processing Gravity and magnetic data were gathered and organized separately in preparation for modeling. Both existing and new data were used in this study. Geosoft?s Oasis Montaj? software was used for processing the data. This software offers many tools for data manipulation. It assists in data correction, filtering, mapping, and gridding. Oasis Montaj? is designed to work with Geosoft?s GM-SYS? profile modeling software. GM- SYS? assists in building and constraining two-dimensional models at depth. The program uses variables such as depth, thickness, density, and susceptibility to define rock units and structures. GM-SYS? and Oasis Montaj? have a data linking feature which allows curser tacking between all open models. This feature is extremely useful for precisely correlating cross-sectional features of the model with anomalies and locations expressed in grid base maps. For the new gravity data acquisition, station locations were recorded digitally using a Trimble? GeoXH GPS in the field. All station readings from the gravimeter were 28 manually recorded. The gravity values being used in this study were first corrected for instrument drift. This correction is made by beginning and ending each day?s measurements at a base station and using these values to calculate the rate of drift during the loop. Drift was assumed to be linear over the time of the loop. Subsequent data processing included the free-air, Bouguer, latitude, tide and terrain corrections. Free?air correction adjusts for the variable station elevations throughout the survey. Bouguer correction adjusts for the mass present between sea-level and the station. Latitude, tide, and terrain correct for variations in the Earth?s shape and rotation, tidal effects, and local topography, respectively. Gravity values from prior surveys were already reduced and new data were merged with these existing data. Several stations established in previous surveys were re-occupied to check agreement between data sets. Average percent error for Complete Bouguer anomaly (CBA) values were 1.37 mGal. This study used only preexisting magnetic data received in the form of a spreadsheet that contained total field anomaly (TFA) values and location coordinates for each measurement. This database was prepared for import by formatting the layout to allow the Oasis Montaj? software to correctly recognize and identify data channels. These preparations created properly arranged databases to be utilized for creating magnetic maps and models. These models were used to locate basin features and to aid in anomaly identification. No data corrections were necessary for magnetic values prior to gridding. However, magnetic curves can be affected if the ambient field and magnetization direction is not vertical (Blakely, 1995). The magnetic data underwent a reduction to the magnetic pole to help prevent this from affecting the models. In this process, the vertical field component is treated as if the source were at the north magnetic 29 pole (90? inclination). This migrates the field from the observed magnetic inclination and declination to a vertically polarized field (Bajgain, 2011). Once the data have been reduced to pole, symmetrical features generate symmetrical anomalies (Figure 14). This aids in the data interpretation as anomalies can more definitively be associated with geometry and or magnetic properties of the rock units (Geosoft technical note). Figure 14. Illustration of how a spherical object produces (a) a slanted anomaly before filtering, and (b) a symmetric anomaly after applying reduction to pole filter (modified by Bajgain, 2011; after Blakely, 1995) Data Modeling Minimum curvature gridding was utilized to generate both gravity and magnetic grid maps. This method fits the smoothest possible surface to the given data values by estimating grid values at points on a coarse grid, while using the inverse distance average of the actual data within a defined search radius (Figure 15). A value equal to the average of all data points is assigned to areas where no data points fall within the assigned radius. An iterative method is used fit the grid to the observed data closest to the grid nodes. The program continues until the best fit is obtained. It then halves the grid cell size and (A) Before filtering (B) After applying filter 30 repeats the process until a best fit surface is achieved under the given parameters (Geosoft: Topics in Gridding) (Figure 15). Figure 15. Diagram showing how the minimum curvature gridding process achieves the smoothest possible surface. Iterations continue until the observed data are closest to the grid nodes (green dots) (Geosoft: Topics in Gridding). Vertical derivative convolution was also used to make gravity and magnetic grid maps. This method generates new maps from previously made minimum curvature grids (Figure 16). Derivatives are effective tools for observing rate of change and can be used to more clearly see edges of geologic contacts/anomalies (Telford et al., 1990). Shallow features are also highlighted in the process. This is achieved by enhancing short wavelength features and filtering long wavelength features in the study areas (Oasis Montaj: Geophysics, 2009). Deeper basement features are generally associated with longer wavelengths, so filtering them out allows shallow anomalies (higher frequency) to become more visible and pronounced (Kinabo, 2007) (Figure 16b). 31 Figure 16. Two minimum curvature grid maps showing the effects of a vertical derivative filter. (a) Top map shows total field anomaly (TFA). (b) Bottom map enhances shallow features as a result of applying the vertical derivative filter. A B 32 Analytic signal grid maps were generated for both gravity and magnetic data sets. These maps can be used to help locate source bodies in the subsurface. This is especially useful when trying to locate magnetic bodies. The analytic signal is calculated by taking the square root of the sum of squares of derivatives (?T/?x) in three orthogonal directions (x, y, and z) where T is total field anomaly (Kinabo, 2007) (Equation 1). The gravity values for the complete Bouguer anomaly were selected for gravity modeling in this study. The initial grid is a map of the values from all merged stations. This grid was then masked to only show values within the two study areas. Auxiliary data were collected and imported into the program to assist in model building. These data sets included shape files, geologic maps, LiDAR DEM?s, and county boundaries. Vertical derivative convolution filters were applied and analytic signal maps were made to aid the interpretation. Magnetic data were processed into multiple grid maps using the minimum curvature method. The TFA channel was selected as the data to grid. The resulting map covers the entire magnetic survey area. As with the gravity data, masking was used to remove all values outside of the two basins of interest. Masking visually allows more detail to be seen as smaller value differences are represented by greater color changes. For both gravity and magnetic data, a poor choice of grid cell size can result in artifacts that can be exaggerated in the filtering process (Figure 17b). An example of this can be seen in Figure 17. A close look reveals a fine vertical rippling pattern that looks smeared, (1) 33 Figure 17. Grid maps showing the effects of poorly chosen grid size. (a) Top map shows vertical rippling in TFA due to poorly chosen grid size. (b) Bottom map shows how a vertical derivative filter exaggerates the effects of an inappropriate choice of grid size. White stripes on the right of both maps result from spacing changes in sampling protocols combined with incorrect grid size. A B 34 especially apparent in areas of high contrast (Figure 17a). The white vertical bars on the right of the map are artifacts resulting from differences in sample spacing combined with the improper cell size. Regriding to a finer cell size removes these ripples (Figure 16a). Manmade features can show up in magnetic data and it is important to identify them. For example, a metal smelting plant forms a sharp, high/low artifact near the center of (Figures 18a and 18b). Masking prior to filtering provides a solution to this problem (Figures 18c and 18d). Ternary maps were created in addition to the other grid maps. Ternary maps combine three grids into one map. The three grids used in the creation of the magnetic ternary map were the original magnetic total field anomaly grid (T), magnetic vertical derivative grid (DZ), and the magnetic analytic signal grid (AS). The ternary map plotted the grids into a cyan-magenta-yellow (CMY) color model. Plotting is done by assigning each of the separate data grids to control one of the desired colors. The coloring method was selected to use histogram equalization (Oasis Montaj: advanced imaging, 2012). When combined, a new map is created that can assist in making geological interpretations about the magnetic anomalies present in the study area by defining boundaries with noticeably different magnetic properties (Kinabo, 2007). 35 Figure 18. (a) TFA map including an artifact caused by a nearby smelting plant. (b) Map showing a vertical derivative filter of Figure 18a. Filtering creates an unwanted bullseye artifact. (c) TFA map where anomaly has been masked. (d) Map showing a vertical derivative filter of Figure 18c. Masking removes the bullseye artifact. Cross-sectional Profiles Multiple cross-sectional models were created in the study areas. These models were constructed from transects of the Bouguer gravity anomaly and the total field anomaly using Geosoft?s GM-SYS? software. GM-SYS? is modeling software that uses A B C D 36 various constraints such as density, depth, and magnetic susceptibility to generate 2D gravity and magnetic profiles. To properly calculate the model?s magnetic response, magnitude and orientation of the Earth?s field are required. One profile in each basin was chosen to be parallel to the potential field gradient to minimize three-dimensional effects. Polygons were constructed in the model to represent major lithologies. Each layer was assigned density and susceptibility properties based on the major rock types assumed to be present. Initial models were constructed by utilizing previously existing geologic cross-sections, seismic data interpretations (if available), and seismic velocity profiles as depth and shape guides for polygon boundaries. Block properties were assigned based on lithologic descriptions. Once the basic geologic cross-section had been constructed and properties were entered, polygons were manipulated to alter the calculated gravity and magnetic curves. The objective was to make the calculated curves match the observed curves by creating geologically plausible models for the area. A ?DC? level is the average level of anomaly response as calculated by GM- SYS?. This level is largely dependent on the total model thickness. GM-SYS? gives several options for dealing with the ?DC? level through applying a ?DC?-shift. For this study, the option to automatically calculate the DC shift was chosen. With this option, GM-SYS? measures the average difference between observed and calculated profiles. Then the program adds the difference to the calculated curve, adjusting the average difference to zero (Geosoft: Knowledge base DC-shift, 2012). 37 Gravity curves were matched first by manually manipulating polygon orientations and properties. Once the observed and calculated gravity curves were in agreement, the total field magnetic curves resulting were adjusted by altering magnetic susceptibilities assigned to the polygons. Magnetic susceptibilities can vary within a given layer so blocks created in the gravity modeling were sometimes subdivided to allow values to range. Each new block maintained the same density as its corresponding layer. The inversion tool was used in GM-SYS? to automatically assign susceptibilities to multiple blocks at the same time within predefined ranges for each layer. This yields a calculated magnetic curve that more closely matches the observed curve than if done using a forward modeling approach. 38 Results Bellingham Bellingham Gravity Grid Maps A gravity map of the Bellingham site was generated to locate anomalies and better define basin boundaries (Figure 19). This map represents gridded complete Bouguer values from the combined gravity data sets. All CBA values are negative. High gravity values are represented by warm colors, and cool colors denote areas of low gravity values. The highest value (-39 mGal) is seen in the northcentral and southwestern regions of the map. The lowest gravity value (- 98 mGal) is seen in the central and northwestern regions. A gravity high dominates the northern border with Canada. There is also a sizable high gravity anomaly in the southeastern corner of the map. Multiple high gravity anomalies occur at and near the intersection of profiles A-A? and B-B? (Figure 19). There is a large low gravity anomaly in the center of the study area. It transitions into a second low anomaly feature to the northeast. The white area to the west is the result of poor data coverage. A vertical derivative convolution filter was used to generate a second grid map (Figure 20). This filter enhances shallow features by removing much of the regional, deeper sources that contribute to the observed gravity values. Highs and lows in the second grid occur in the same general areas as the previous model; however, the filtered map provides greater detail associated with shallow 39 Figure 19. Complete Bouguer anomaly (CBA) gravity map of the Bellingham study area. Black dots=station locations. A-A? and B-B?=profile locations. Solid blue lines=faults. Dashed blue lines=inferred faults. Abbreviations: BB=Birch Bay, BH=Bellingham, BHB=Bellingham Bay, BL=Blaine, DH=Drayton Harbor. 40 Figure 20. Gravity (CBA) map of the Bellingham study area with vertical derivative filter applied. A-A? and B-B?=profile locations. Solid blue lines=faults. Dashed blue lines=inferred faults. White arrows=previously unseen high gravity anomalies. Abbreviations: BB=Birch Bay, BH=Bellingham, BHB=Bellingham Bay, BL=Blaine, DH=Drayton Harbor, DHF=Drayton Harbor fault, BBF=Birch Bay fault. 41 structures. Gravity values on the second map range from -0.0122 mGal/km (low) to 0.0150 mGal/km (high). A gravity low is seen on the shoreline to the west that was not evident in the unfiltered model. The highs at and near the intersection of the profiles are sharper and more prominent. Bellingham Magnetic Grid Maps A magnetic map of the Bellingham study area was generated to locate anomalies to assist in identifying structures (Figure 21). Most of the magnetic data used to create this grid map were previously acquired through an aerial survey. Warm colors represent magnetic highs and cool colors show magnetic lows. The lowest TFA values (-510 nT) are observed to the south and east of the map. The highest observed TFA (919 nT) are shown in pink and occur in several areas in the northwestern half of the map. Data coverage allows for the same 3358 km? area to be gridded. Gaps in data coverage result in artifacts in the eastern part of the study area. A slight ripple effect is observed in the extreme southeast corner of the map. This is an artifact that occurs due to poor data coverage to the east. The whited-out area along the western shoreline is a mask that covers an extremely high anomaly caused by a metal smelting facility. A vertical derivative convolution filter was used to generate a second TFA map (Figure 22). This grid map enhances shallow, shorter wavelength features. Values range from -1.102 nT/km (low) to1.943 nT/km (high). Highs and lows occur throughout the map and are no longer confined to just a few areas. Masking the metal smelting facility enhanced lineations nearby that would have otherwise not been visible. Rippling to the east has been amplified as a result of the filtering process and spacing differences. 42 Figure 21. Magnetic (TFA) map of the Bellingham study area. A-A? and B-B?=profile locations. Solid blue lines=faults. Dashed blue lines=inferred faults. Abbreviations: BB=Birch Bay, BH=Bellingham, BHB=Bellingham Bay, BL=Blaine, DH=Drayton Harbor. 43 Figure 22. Magnetic (TFA) map of the Bellingham study area with vertical derivative filter applied. A-A? and B-B?=profile locations. Solid blue lines=faults. Dashed blue lines=inferred faults. Abbreviations: BB=Birch Bay, BH=Bellingham, BHB=Bellingham Bay, BL=Blaine, DH=Drayton Harbor. 44 A ternary map was generated by combining the previous two grids with an analytic signal grid (Figure 23). Total field (T) is shown in cyan, vertical derivative (DZ) is yellow, and analytic signal (AS) is magenta. This map emphasizes areas where anomalies are common between grids. Bellingham Cross-sections Two gravity and magnetic profiles (A-A? and B-B?) were constructed for the Bellingham study area (Figures 24 and 25). A-A? crosses the Bellingham basin from south to north and is 28.2 km long. Coordinates are defined using the NAD_83 datum. The study area follows the UTM zone 10N projection method. The coordinates for the starting point of A-A? are (531177.7, 5401249) m, and the coordinates of the end point are (530497, 5429384) m. B-B? runs from southwest to northeast diagonally though the basin and is 41 km long. The coordinates for the starting point of B-B? are (525202.8, 5402913) m, and the coordinates for its endpoint are (557119.4, 5428703) m. The decision was made to place both profiles in areas where both gravity and magnetic data coverage was most dense. Figures 24 and 25 show the potential field properties of the rocks in the area. Positive or negative contrasts are shown as relative to an assumed background value. The properties of specific rock units are assigned to various polygons or ?blocks? generated in the cross-sectional modeling process. Both cross-sections extend 12 km below the surface. Surface contacts of rock units in and outside of the modeled cross-sections were used to constrain boundaries between polygons and the location of structures, such as faults and folds. 45 Figure 23. Magnetic ternary map of the Bellingham study area generated by combining the grids from Figures 21 and 22 with an analytic signal grid. A-A? and B-B?=profile locations. Solid blue lines=faults. Dashed white lines=inferred faults. Abbreviations: BB=Birch Bay, BH=Bellingham, BHB=Bellingham Bay, BL=Blaine, DH=Drayton Harbor, T= Total field, DZ=vertical derivative, AS=analytic signal. 46 Figure 24. Gravity and magnetic profiles with geologic cross-section for A-A?. Abbreviations: QS= Quaternary Sediments, CN= Chuckanut Formation, EA= Easton Metamorphic Suite, BP= Bell Pass M?lange, CH= Chiliwack Group, NK= Nooksack Formation, BBF=Birch Bay Fault, DHF=Drayton Harbor Fault. 47 Figure 25. Gravity and magnetic profiles with geologic cross-section for B-B?. Abbreviations: QS= Quaternary Sediments, CN= Chuckanut Formation, EA= Easton Metamorphic Suite, BP= Bell Pass M?lange, CH= Chiliwack Group, NK= Nooksack Formation, BBF=Birch Bay Fault. 48 Units (polygons) within the cross-section were assigned estimated values based on lithologies described in existing maps and publications. While these sources provided a guide to constrain values, the lack of publicly available data for the Bellingham basin at a depth of 12 km limited the number of model constraints. None of the thirteen 7.5? geologic quadrangle maps present in the area have been mapped and/or published on a 1:24,000 scale. Only a 1:100,000 scale Bellingham map was available for surface geology (Lapen, 2000). Major geologic contacts are still visible at this scale though greater detail would have been helpful for constructing polygons for the cross-sections. Gravity and Magnetic Models: A-A? Profile Gravity and magnetic data for the south to north transect (A-A?) are modeled by six, roughly horizontally bodies (Figure 24). Very little well data are available in the immediate transect areas. All logs that are on or near this transect are from drill holes less than one kilometer in depth. The deepest available well log in the basin is three kilometers deep and does not reach basement rocks and thus provide only a minimum estimate of depth to basement. The observed gravity is the lowest (-76.1 mGal) at the southern starting point of A-A?, and is the highest (-52.4 mGal) near the northern endpoint. All the values contained within A-A? are negative. The gravity curve is dominated by a gradual increase with northward progression. One notable exception to this trend occurs in the southern quarter of the transect, where a notable gravity high occurs. This gravity spike coincides with the noticeable magnetic jump of 20 nT over one kilometer seen in the observed magnetics curve. A-A? contains no abrupt magnetic peaks but it does show several other locations where TFA steps up or down. One of these steps changes occurs on the northern end of the transect. It coincides with a linear scarp visible 49 in the LiDAR digital elevation models (Figure 8). Magnetic susceptibilities for the six bodies were estimated to fall within certain ranges. These ranges were constrained by accepted values for the general rock types represented by each body. Densities were also assigned to polygons based on acceptable values for the lithologies represented. The uppermost layer, QS, consists of Quaternary glacial sediments (Table 3). It is the youngest unit in the area and overlies the Chuckanut Formation (CN). CN?s protoliths are arkosic and lithic sandstone, siltstone, conglomerate, and coal. Deposition for these units is considered to have been mostly fluvial with some marine. The Easton metamorphic suite (EA) underlies CN and is made up of metabasite, siliceous carbonaceous phyllite, metagraywacke, metabasalt, gabbro-tonalite plutons, and metatuffs. The Bell Pass m?lange (BP) is the next deepest block (Brown and Dragovich, 2003). The Bell Pass m?lange is thought to grade into the Orcas Chert somewhere beneath the study area, so constituents from both units were considered in the model. Their protoliths include ribbon chert, mudstone, pillow basalt, limestone, and volcanic lithic clastic rocks. The Chiliwack Group (CH) underlies BP and it grades into the East Sound Group to the west. Rock types from both units where included, since they also meet beneath the study area. As a result, this block represents a diverse assortment of lithologies. These include volcanic lithic sandstones and siltstones, basaltic to dacitic volcanic rocks, fossiliferous to reefoidal limestones, andesitic to dacitic rocks, with lesser limestone lenses, shale, argillite, graywacke, and conglomerate. The bottom block in the cross-section represents the Nooksack Formation (NK). Its protoliths are largely made up of marine volcanic-lithic sandstone, siltstone, local conglomerate, and a local andesitic- dacitic volcanic member (Brown and Dragovich, 2003). 50 Table 3. Blocks/polygons and their depths for A-A?, with densities, susceptibilities, and lithologies listed for each body (Brown and Dragovich, 2003; Brown and Gehrels, 2007). Body symbol Unit name Depth to top (km) Maximum thickness (km) Modeled rock density (g/cm?) Modeled rock susceptibility (cgs) Rock description QS Quaternary Sediments 0 3.31 2.0 0.000016- 0.00098 Glacial Sediments CN Chuckanut Formation 0.66- 3.31 3.24 2.5 0.000098- 0.001102 Arkosic and lithic sandstone, siltstone, conglomerate, and coal. Mostly fluvial deposition with some marine. EA Easton Metamorphic Suite 3.74- 5.42 1.93 2.8 0.000232 Metabasite, siliceous carbonaceous phyllite, metagraywacke, metabasalt, gabbro-tonalite plutons, and metatuffs. BP Bell Pass M?lange/ Orcas Chert 5.26- 7.33 1.97 2.7 0.004002 Ribbon chert, mudstone, pillow basalt, limestone, volcanic lithic clastic rocks. CH Chiliwack Group/ East Sound Group 6.72- 9.08 1.93 2.8 0.00765 Volcanic lithic sandstones and siltstones, basaltic to dacitic volcanic rocks, fossiliferous to reefoidal limestones. Andesitic to dacitic rocks, with lesser limestone lenses, shale, argillite, graywacke, and conglomerate. NK Nooksack Formation 8.26- 10.67 N/A 2.67 0.000801 Marine volcanic-lithic sandstone, siltstone, and local conglomerate. Local andesitic-dacitic volcanic member. 51 Gravity and Magnetic Model: B-B? Profile Gravity and magnetic data for the B-B? profile were also modeled using six roughly horizontal bodies (Table 4 and Figure 25). B-B? runs from southwest to northeast across the basin. This profile intersects A-A?, adding the constraint that both models must tie in this area. Very little well log data are available along this transect. All logs that are on or near this transect are less than one kilometer in depth. The deepest available well log in the area remains the three kilometer deep well mentioned previously and does not encounter basement. The observed gravity is the lowest (-77.2 mGal) at the southwestern starting point of B-B?, and highest (-61.4 mGal) at 5 km from the northeastern endpoint. The observed gravity curve for B-B? has a similar regional trend as the curve in A-A?. Gravity is lowest to the southwest and gradually increases towards the northeast. A notable deviation from the region trend occurs at 4 km from the southwestern starting point. There is a 6.5 mGal increase over a distance of 1.9 km at 6 km in the model. A second gravity peak occurs 23 km from the southwestern endpoint. It is defined by a drop in gravity of -1.48 mGal over a 3.18 km and is less prominent than the previously discussed feature. The highest gravity value of -61.4 mGal is observed at the peak of a third notable gravity anomaly. It is located at 34.6 km from southwest endpoint in the model and is defined by a slight increase in gravity followed by a decrease. Minor fluctuations, relative to these three spikes in observed values occur between the second and third peaks. All three gravity peaks coincide with changes along the observed magnetic profile. Observed magnetic data are lowest (-25.3 nT) in the center of transect 52 Table 4. Blocks/polygons and their depths for B-B?, with densities, susceptibilities and lithologies listed for each body (Brown and Dragovich, 2003; Brown and Gehrels, 2007). Body symbol Unit name Depth to top (km) Maximum thickness (km) Modeled rock density (g/cm?) Modeled rock susceptibility (cgs) Rock description QS Quaternary Sediments 0 1.93 2.0 0.000052- 0.000513 Glacial Sediments CN Chuckanut Formation 0.63- 1.93 4.85 2.5 0.000035- 0.00058 Arkosic and lithic sandstone, siltstone, conglomerate, and coal. Mostly fluvial deposition with some marine. EA Easton Metamorphic Suite 3.6- 5.77 2.42 2.8 0.000221- 0.000462 Metabasite, siliceous carbonaceous phyllite, metagraywacke, metabasalt, gabbro-tonalite plutons, and metatuffs. BP Bell Pass M?lange/ Orcas Chert 5.52- 7.53 2.55 2.7 0.000353- 0.007348 Ribbon chert, mudstone, pillow basalt, limestone, volcanic lithic clastic rocks. CH Chiliwack Group/ East Sound Group 6.9-9.7 2.3 2.8 0.000238 Volcanic lithic sandstones and siltstones, basaltic to dacitic volcanic rocks, fossiliferous to reefoidal limestones. Andesitic to dacitic rocks, with lesser limestone lenses, shale, argillite, graywacke, and conglomerate. NK Nooksack Formation 8.61- 11.87 N/A 2.67 0.000796 Marine volcanic-lithic sandstone, siltstone, and local conglomerate. Local andesitic-dacitic volcanic member. 53 B-B?, and highest (302.68 nT) at the northeastern endpoint. All but eighteen of the 98 values are positive. The magnetic curve for B-B? is more complex regionally than the one generated for A-A?. The observed anomaly is low at the southwestern edge and then jumps 74.2 nT over 2 km. A regional trough occurs at the center of the profile. It is formed from a gradual decrease in TFA followed by a sharp increase of 88.7 nT over 3.57 km. This magnetic low coincides with the second observed gravity peak. The curve gets steeper over the last 6 km. Over this distance, there is an increase in susceptibility of 200.45 nT. Numerous smaller fluctuations occur between these three points. The first gravity anomaly peak coincides with one of these. It is defined by a relatively slight increase in magnetic TFA of 8.12 nT followed by a decrease of 20.98 nT. Muckleshoot Muckleshoot Gravity Grid Maps Gravity grid maps were constructed for the Muckleshoot study area to help identify basin extents and subsurface features. Figure 26 is a color shaded grid map of all new and old complete bouguer gravity measurements. This study area is approximately 50 km by 44 km (2200 km?). All observed gravity measurements are negative. Values range from -28 mGal (high) to -104 mGal (low). The lowest values show up in the center of the map as well as the southeast corner. A large high is seen along the western edge of the map. Another high is visible wrapping along the eastern edge of the previously mentioned low. It enters from the north and ends before exiting the southern edge of the map. A second grid map was generated form the first using a vertical derivative 54 Figure 26. Gravity (CBA) map of the Muckleshoot study area. Black dots=station locations. C-C?=profile location. Solid blue lines=faults. Abbreviations: CV=Covington, K=Kent, PS=Puget Sound, R=Renton, T=Tacoma, WRF=White River fault. 55 convolution filter (Figure 27). This filter enhances the shallow gravity signals. Highs and lows occur in roughly the same areas as the prior map. Gravity values on this new map range from -0.0140 mGal/km (low) to 0.0109 mGal/km (high). A notable low shows up along the northern edge of the map. It enters from the west and trends east-west before deviating southeast into the central low. The southwest corner of the area also has a low that is not in the previous map. A shaper contrast can be seen in some areas of the map. This is most notable in the southeast corner along the White River fault. Several other isolated highs occur in the middle of the central low. The southern edge of the map shows a high that appears to be a boundary of the central low. 56 Figure 27. Gravity (CBA) map of the Muckleshoot study area with vertical derivative filter applied. C-C?=profile location. Solid blue lines=faults. Arrows=aligned gravity highs. Abbreviations: CV=Covington, K=Kent, PS=Puget Sound, R=Renton, T=Tacoma, WRF=White River fault. 57 Muckleshoot Magnetic Grid Maps A magnetic grid was constructed to aid in locating and identifying magnetic anomalies that could be associated with subsurface features (Figure 28). The lowest magnetic TFA values (-2316 nT) are observed as blue areas to the south and east of the map. The highest TFA values (520 nT) are shown in pink and occur along the western edge and in several areas in the southeastern corner of the map. Several white pixels and lines in the southeast corner are artifacts data coverage gaps. Another map was constructed using a vertical derivative filter (Figure 29). Differences in data coverage cause the appearance of the same white anomalies along with slight rippling in the immediate surrounding areas. The filter enhanced TFA contrasts caused by shallow source magnetic anomalies. Values range from -9.420 nT/km (low) to 1.897 nT/km (high). Highs and lows occur throughout the map. End values on the scale bar reflect the highest and lowest values in the area. Numerous high contrast anomalies dominate the borders of the map while lower contrasts features occur more centrally. A three grid composite map was generated by combining vertical derivative, total field, and analytical signal grids (Figure 30). This third map is similar in appearance to the vertical derivative grid. Vivid reds and blues represent highs and lows that occur in all three grids. Muted muddy greens and browns mark areas in which anomalies are dissimilar between grids. 58 Figure 28. Magnetic (TFA) map of the Muckleshoot study area. C-C?=profile location. Solid blue lines=faults. Abbreviations: CV=Covington, K=Kent, PS=Puget Sound, R=Renton, T=Tacoma, WRF=White River fault. 59 Figure 29. Magnetic (TFA) map of the Muckleshoot study area with vertical derivative filter applied. C-C?=profile location. Solid blue lines=faults. Abbreviations: CV=Covington, K=Kent, PS=Puget Sound, R=Renton, T=Tacoma, WRF=White River fault. 60 Figure 30. Three grid composite of the Muckleshoot study area generated by combining the grids from Figures 28 and 29 with an analytic signal grid. C-C?=profile location. Solid blue lines=faults. White arrows= highlighted sharp anomaly. Abbreviations: CV=Covington, K=Kent, PS=Puget Sound, R=Renton, T=Tacoma, WRF=White River fault. 61 Muckleshoot Cross-section A cross-section (C-C?) was constructed along a 42 km west-east transect through the Muckleshoot study area (Figure 31). The coordinates for the starting point are (547248.8, 5238997) m, and the coordinates for its endpoint are (591509.3, 5238997) m. The transect was oriented parallel to the magnetic gradient and crosses the basin where gravity and magnetic coverage was most dense. The profile reaches depths of 10 km. Surface geology, rock descriptions, and velocity profiles were used in assigning estimated physical properties for cross-section layers. The Muckleshoot gravity and magnetic data was modeled by six main bodies (Table 5). Gravity and Magnetic Model: C-C? Profile Observed gravity and magnetic values were profiled along the transect (Figure 31). Observed gravity is lowest at the center of the profile (-95 mGal) and highest (-45 mGal) at the starting point. From the starting point, gravity decreases with eastward progression until the midpoint. After the midpoint, gravity increases gradually for 13 km and then starts falling again. A notable trough occurs 13 km from the start where gravity decreases 9 mGals over a distance of 1.2 km. The magnetic profile begins 11 km from the start. It has been intentionally shortened in the west to avoid magnetic sources along the Interstate 5 corridor. This curve begins with a steep decrease in TFA followed by a gradual increase until 32.7 km from the start. This is where another sharp drop occurs followed by dramatic fluctuations with an increasing trend. The highest TFA (-31.65nT) occurs 29.5 km from the start and coincides with a vertical igneous dike (Purple feature 62 shown in Figure 31). This peak marks a 46.9 nT increase over 1 km. The lowest observed TFA is -171 nT and is 36.75 km from the start. 63 Figure 31. Gravity and magnetic profiles with geologic cross-section for C-C?. Abbreviations: G1=Glacial and post glacial strata 1, G2=Glacial and post glacial strata 2, OP=Ohanapecosh Formation, P1= Puget Group 1, P2= Puget Group 2, CR= Crescent Formation, TF/WRF=Tacoma fault/ White River fault, vertical purple polygon=igneous dike. 64 Table 5. Blocks/polygons and their depths for C-C?, with densities, susceptibilities and lithologies listed for each body (Johnson et al., 2004; Jones, 1996; Tre?hu et al., 1994; Finn, 1990; Lees and Crosson, 1990; Yount et al., 1985). Body symbol Unit name Depth to top (km) Maximum thickness (km) Modeled rock density (g/cm?) Modeled rock susceptibility (cgs) Rock description G1 Glacial and Post Glacial Strata -0.7- 0 2.34 1.9 0.00001- 0.00007 Clay, silt, and sand derived from rivers and erosion of Pleistocene deposits. G2 -1.16- 2.17 2.82 2.0 0.00002- 0.00028 Clay, silt, and sand derived from rivers and erosion of Pleistocene deposits. transitional and interglacial deposits OP Ohanapecosh Formation -1.39- -1 1.17 3.0 0.0005 Volcanoclastic deposits and flood basalts P1 Puget Group -2.7- 3.3 2.53 2.22 0.00006- 0.00078 Volcanics and non-marine sedimentary rocks. Correlatives to the Blakely Formation (Oligocene) P2 -1.05- 6 3.37 2.7 0.0001- 0.00045 Volcanics and non-marine sedimentary rocks. Correlatives to the Tukwilia Formation (Eocene) CR Crescent Formation -2.59- 8.59 N/A 3.0 0.00127 Marine basaltic rocks that are dense and magnetic. Sedimentary rocks are inter- bedded at the top of the formation. 65 Discussion Bellingham basin New gravity data added in this study help to refine the geometry of the Bellingham basin and elucidate several previously identified structures within and adjacent to the basin (Figure 19). These structures include the Birch Bay fault and the Drayton Harbor fault. Results from this study suggest a basin that is much smaller than that which was previously proposed (Figure 32). In previous studies, the Bellingham basin is depicted as widening towards the coast and extending northward into Canada (Figure 6) (Kelsey et al., 2012). Figure 32a shows that the basin has a broad arcuate shape that likely extends into Canada on its western and eastern flanks, dropping to approximately 13 km south of the Canadian border in its central area. The gravity low defining the basin, although large, is not continuous, suggesting that the Bellingham basin is comprised of at least two, possibly three, subbasins. Station density along the southern edge of the basin is low and therefore the geometry of the basin in this section is not well constrained. Similarly, the western boundary of the basin is not well defined by the available data. Given that the majority of stations measured on the western coastline show lower gravity values relative to adjacent landward stations to the east, it is possible that the basin extends westward into Puget Sound. Recent marine based studies focused in these offshore areas will better 66 A 67 Figure 32. Vertical derivative maps of gravity (a) and magnetic (b) data for the Bellingham basin with locations of structures and boundaries as determined by this study. White dashed line=basin boundary. Black dotted line=possible unidentified northwest striking structure. Red dashed lines=location of faults proposed by Kelsey et al. (2012). (a) Map showing gravity. (b) Map showing magnetics. DHF=Drayton Harbor fault. BBF=Birch Bay fault. B 68 constrain the basin?s western boundary (R. Blakely, pers. comm., 2012). Although there is some debate over possible westward extension of the Bellingham basin, the gravity maps presented here do not support the location of the previously proposed northern basin boundary. Birch Bay fault The newly acquired gravity data reveal several high gravity anomalies (white arrows in Figure 20) north of the Bellingham Bay that had not been seen before. These seem to be associated with the Birch Bay fault (BBF) proposed by Kelsey et al. (2012). This fault trends southeast from Birch Bay and appears to offset the two gravity highs (Figure 20). Magnetic data show a linear feature that also appears bisected and offset along the proposed BBF (Figures 21 and 22). The trace of the BBF is well defined on the coast in the vertical derivative maps but becomes less obvious as it progresses inland (Figure 22). The magnetic ternary map shows a sharp contrast in the vicinity of Birch Bay (BB) further supporting the BBF?s existence and orientation near the coast, but is not distinctive further inland to the east (Figure 23). The proposed BBF can be projected to the intersection of profiles A-A? and B-B?. At the projected location of the BBF on the cross-sectional profiles, there is a pronounced high in the observed gravity curves. This high is coincident with the crest of an anticline modeled with a down-to-the-south offset along the BBF (Figures 24 and 25). In paleoseismic studies, conducted by near the coast of BB, Kelsey et al. (2012) also suggest a north-side-up orientation of the BBF. A steeply dipping south vergent reverse fault is consistent with north-northeast-south-southwest regional compression (Zoeback, 1992; Wells et al., 1998). 69 Hurst (1991) published an industry seismic line in which the Birch Bay anticline and several faults are visible. These faults are shown to cut through the anticline at high angles and are consistent with the proposed model. It is also possible, however, to match the observed gravity with a series of very steep non-faulted anticlines and synclines. Given the extreme dip and size of these folds and evidence of offsets along the coast (Kelsey et al., 2012) the faulted fold model was preferred. Drayton Harbor fault Figure 20 shows a bisected gravity high near Drayton Harbor that directly correlates with a LiDAR scarp (Figure 8) and the postulated location of the DHF for a distance of 15 km (Kelsey et al., 2012). At this same location, a pronounced magnetic lineament is also visible in data plotted using a vertical derivative convolution filter (Figure 22). Figure 24 shows the DHF projected onto the cross-section along A-A?. The fault?s location and orientation is largely based on the LiDAR DEM?s and the filtered magnetic maps (Figures 8 and 22). Observed gravity data along the profile do not show an obvious deviation at this location. However, observed magnetic data do show a slight deviation along a generally northward decreasing trend at the projected location of the DHF (Figure 24). Attempts to model this last section of the profile without a south- dipping reverse fault were unsuccessful. Bellingham Bay feature Closer observation of the vertical derivative maps suggests that a third previously unidentified structure may be present (black dotted line in Figure 32). This feature trends north-northwest and incorporates both Birch Bay gravity highs (Figure 32a). It can be 70 traced southward through the intersection point of both profiles. Linear lows occur along the west side of these gravity highs. Similarly, a south trending magnetic lineament extends though the Bellingham Bay (Figure 32b). This feature could be an extension of one of the known folds southeast of Bellingham Bay (Figure 3). A fault with an east-side up, west-side down orientation could also explain the gravity and magnetic data. It is possible that the feature represents the southern extent of the BBF. It is worth noting that the anomaly locally coincides with the Nooksack River delta and appears to follow the river for a short distance. Sediments containing significant amounts of magnetic mineral grains could be present. Magnetic sedimentary deposits could explain some of these smaller magnetic anomalies but not the observed gravity data or the north-south trending feature as a whole. The river could also be following a fault. Muckleshoot basin New gravity data added in this study area help to gain a better understanding of the Muckleshoot basin and illuminate several previously identified structures within and adjacent to the basin (Figure 2). These structures include the Tacoma and White River faults. Previous studies show the Muckleshoot basin as an inland isolated gravity low east of the Puget Sound (Johnson et al., 2004) (Figure 9a). The gravity low appears to be subdivided into two subbasins by a northwest-trending high (Figures 9 and 10). The Muckleshoot gravity maps prepared for this study include the basin boundaries as well as the White River fault (WRF) and a possible extension of the Tacoma fault. In the CBA gravity map, the Muckleshoot basin is clearly defined by a gravity low seen in the center of Figure 26. This low changes orientation north of profile 71 C-C?. This subdivision could be the northwest extension of the WRF, the southeast extension of the Tacoma fault system (Figure 2), or a basement high that has not been previously identified. The vertical derivative map of the gravity data further supports the existence of at least two subbasins within the area currently recognized as the Muckleshoot basin (Figure 27). Numerous localized gravity highs align through the basin, with a distinct change in trend near the location of profile C-C? and on strike with the WRF (Figures 27 and 33). Several magnetic lineaments are visible within the basin, although the trend of these lineations is not uniform. The filtered magnetic TFA map most notably reveals a northwest trending lineament between Renton and Kent, Washington, which is also a prominent feature visible in the composite grid map (white arrows in Figures 29 and 30). A magnetic lineament possibly corresponding to the northwest extension of the WRF crosses the profile C-C? near its midpoint and appears as a fault (TF/WRF) in the geologic cross-sectional model (Figures 31 and 33b). Three seismic profiles from Liberty (2007) provide a guide to model subsurface geology (Figure 10). Liberty?s A-A? profile is at the same latitude as C-C? profile, though C-C? is shorter and only includes the Muckleshoot basin. Gravity and magnetic profiles along C-C? reveal that the basin is more complex than the Bellingham basin (Figure 31). The gravity profile roughly matches the observed values but reveals that several unknown features or orientations must exist. Few publications were available in this area to provide subsurface constraints on the models. In addition, models are two-dimensional and are assumed to be approximately continuous in the third dimension. Nearby features that are 72 A 73 Figure 33. Maps with vertical derivative filter, showing locations of structures and boundaries as determined by this study for the Muckleshoot basin. White dashed line=basin boundary. Black dotted line=possible unidentified northwest striking structures. (a) Map showing gravity. (b) Map showing magnetics. TF/WRF=Tacoma fault/White River fault seen in cross-section (Figure 31). B 74 not present along a profile can still affect the observed TFA and could be one source of error for the cross-section. The TF/WRF was projected onto the western half of profile C- C? and it is interpreted to be a possible link between the Tacoma fault and WRF systems (Figure 33). It is modeled as thrust fault that dips steeply to the north, consistent with previous interpretations (Figure 31) (Liberty, 2005; Johnson et al., 2004). Although the cross-sectional model fits the observed data fairly well, it fails to match the observed gravity data in some areas. This misfit is especially apparent in the northern portion of the profile where geologic maps indicate the presence of mafic dikes and plutons (Figure 31). More work needs to been done in this area to better determine the nature of this basin at depth. Disclaimer Gravity and magnetic profiles are used to build geologic cross-sections that reflect the observed changes in density and magnetic susceptibility. There are some limitations when interpreting gravity and magnetic data. The source of a given anomaly is unknown along with its geometry. Models are also limited by constraints at depth. Lateral variations in geometry, density, and susceptibility contribute to the values observed. It is important to realize that these models only express one explanation for subsurface geometries. Multiple combinations of features and geometries can produce the same gravity and magnetic model. More subsurface constraints are needed to further support model geometry and probability. 75 Work on this thesis was partially supported by the U.S. Geological Survey (USGS), Department of the Interior, under USGS award number (Wolf, USGS- G11AP20042). The views and conclusions contained in this document are those of the authors and should not be interpreted as necessarily representing the official policies, either expressed or implied, of the U.S. Government. 76 Conclusions Gravity and magnetic modeling of both the Bellingham and Muckleshoot basins in western Washington yield new insight into the configuration of these basins and important faults that traverse them. Boundaries of the Bellingham basin were determined through gravity models and discovered to be different than previously thought. The Bellingham basin appears to be smaller and possibly subdivided. Some boundaries in the east coincide with faults, suggesting that faulting had a role in determining basin geometry. Results from the Bellingham study area support the existence of at least two faults beneath the Bellingham basin, as proposed previously by Kelsey et al. (2012). This study clarifies the location of the BBF and DHF beneath basin sediments. The BBF is modeled as a reverse fault that dips steeply to the north. It deeply penetrates an anticlinal structure that trends northwest through the study area, and is consistent with a compressional tectonic environment. The DHF was modeled as a reverse fault that dips steeply to the south. Although its expression is clearly visible in filtered magnetic data and LiDAR DEM?s, it is not apparent in the gravity data. Lack of a prominent gravity anomaly could suggest that this is a smaller structure. 77 Gravity and magnetic grid maps also reveal evidence for a structure that has not yet been identified. This feature trends north-northwest and cuts through the Bellingham Bay. The anomaly could be a buried continuation of known folds that exist to the southeast. Both the filtered magnetic grid and the ternary map show northwest-trending magnetic lineaments that appear to be offset by this new feature, suggesting that it could be a fault. If this is a fault, its proximity to the town of Bellingham might pose a significant hazard. Gravity data from the Muckleshoot study area suggest basin boundaries that agree with previous interpretations. Gravity models support the existence of at least two subbasins within the Muckleshoot. The location and orientation of these subbasins suggests that faulting may have segmented a larger basin. This segmentation could be an expression of the westward continuation of the WRF, the eastward continuation of the Tacoma fault, or a connecting structure between these two fault systems. Both study areas would benefit from further study. The addition of marine gravity data, especially west of the Bellingham site area, would provide information on the geometry of the basin to the west. More subsurface constraints, in the form of deep wells and detailed geologic mapping, would be beneficial and assist with future modeling and interpretations, as this study was limited by the lack of data from deep wells and by having access to only one industry seismic line. 78 Knowing the location of faults and how they behave is critical information for evaluating earthquake hazard. The basins surrounding the Puget Sound, Washington, are situated in an active subduction zone, which has created a complex network of faults. The Puget Sound region is home to millions of people and major cities that would be affected in a major earthquake event, making it important to gain an understanding of the nature of these faults and how they interconnect. 79 References Atwater, B.F., Tuttle, M.P., Schweig, E.S., Rubin, C.M., Yamaguchi, D.K., Hemphill- Haley, E., 2003, Earthquake recurrence inferred from paleoseismology, Developments in Quaternary Sciences, v. 1, p. 331-332. Bajgain, S.K., 2011, Gravity and magnetic modeling of basement beneath Alabama Gulf Coastal Plain [MS thesis]: Auburn University, p. 24-74. Blakely, R.J., 1995, Potential theory in gravity and magnetic applications: New York, Cambridge Univ. Press, p.1-441. Blakely, R.J., John, D.A., Box, S.E., Berger, B.R., Fleck, R.J., Ashley, R.P., Newport, G.R., and Heinemeyer, G.R., 2007, Crustal controls on magmatic-hydrothermal systems: A geophysical comparison of White River, Washington, with Goldfield, Nevada, Geosphere; April 2007, v. 3, no. 2, p. 91?107. Blakely, R.J., Wells, R.E., and Weaver, C.S., 1999, Puget Sound aeromagnetic maps and data, U.S. Geological Survey Open File Report, 99-514. Blakely, R.J., Wells, R.E., Weaver, C.S., and Johnson, S.Y., 2002, Location, structure, and seismicity of the Seattle fault zone, Washington; evidence from aeromagnetic anomalies, geologic mapping, and seismic-reflection data, Geological Society of America Bulletin, v. 114, no. 2, p. 169-177. Box, S.E., John, D.A., Ashley, R.P., and Rytuba, J.J., 2003, Kinematic analysis of fault- slip data from mining districts in the Tertiary Cascades, OR and WA: Geological Society of America Abstracts with Programs, v. 35, no. 6, p. 512. 80 Brocher, T.M., Parsons, T., Blakely, R.J., Christensen, N.I., Fisher, M.A., Wells, R.E., and SHIPS Working Group, 2001, Upper crustal structure in Puget Lowland, Washington: Results from 1998 Seismic Hazards Investigation in Puget Sound, Journal of Geophysical Research, v. 106, p. 13,541?13,564. Brown, E.H., and Dragovich, J.D., 2003, Tectonic elements and evolution of northwest Washington, Washington Division of Geology and Earth Resources, Geologic map GM-52. Brown, E.H., and Gehrels, G.E., 2007, Detrital zircon constraints on terrane ages and affinities and timing of orogenic events in the San Juan Islands and North Cascades, Washington, Canadian Journal of Earth Sciences, v. 44, p. 1375-1396. Clague, J.J., and James, T.S., 2002, History and isostatic effects of the last ice sheet in southern British Columbia, Quaternary Science Review, v. 21, p. 71?87, doi:10.1016/S0277-3791(01)00070-1. Dater, D., Metzger, D., and Hittelman, A., 1999, Land and Marine Gravity CD-ROMs, U.S. Department of Commerce, National Oceanic and Atmospheric Administration, National Geophysical Data Center, Boulder, CO. Easterbrook, D.J., 1963, Late Pleistocene glacial events and relative sea level changes in the northern Puget lowland, Geologic Society of America Bulletin, v. 74, no. 12, p. 1465?1484, doi:10.1130/0016-7606(1963)74[1465:LPGEAR]2.0.CO;2. Easterbrook, D.J., 1976, Geologic map of western Whatcom County, Washington, U.S. Geological Survey Miscellaneous Investigation Map, I-854-B, scale 1:62,500. Finn, C., 1990, Geophysical constraints on Washington convergent margin structure, Journal of Geophysical Research, v. 95, p. 19,533 ? 19,546. 81 Finn, C., Phillips, W.M., and Williams, D.L., 1991, Gravity anomaly and terrain maps of Washington, scale 1:500000 and 1:1,000,000, U.S. Geological Survey Geophysical Investigation Map, GP-988. Gallup, W.B., 1957, Will these sedimentary basins [Pacific Coast] prove oil bearing?, Oil Gas Journal, v. 55, p. 142. Geosoft, 2012, Knowledge base: DC-Shift in GM-SYS profile modeling: http://www.geosoft.com/support/knowledge-base/generalinfo/DC-Shift-in-GM- SYS-Profile-Modelling (accessed January 2013). Geosoft technical note: Reduction to the Magnetic Pole for High Frequency Amplitude Anomalies: http://www.geosoft.com/media/uploads/resources/technical- notes/Reduct_to_Mag_Pole_High_Freq_High_Amp_Anomalies.pdf (accessed January 2013). Geosoft: Topics in Gridding: Technical Workshop manual, Geosoft: http://www.geosoft.com/media/uploads/resources/technical- papers/topicsingriddingworkshop.pdf (accessed January 2013). Haugerud, R.A., Harding, D.J., Johnson, S.Y., Harless, J.L., and Weaver, C.S., 2003, High- Resolution LiDAR Topography of the Puget Lowland, Washington, Geological Society of America Today, v. 13, no. 6, p. 4-10. Hurst, P.D., 1991, Petroleum geology of the Bellingham basin, Washington, and evaluation of the AHEL and Partners Birch Bay No. 1 well, Washington Geology, v. 19, p. 16?18. James, T., Gowan, E.J., Hutchinson, I., Clague, J.J., Barrie, J.V., and Conway, K.W., 2009, Sea-level change and paleogeographic reconstructions, southern Vancouver Island, British Columbia, Canada, Quaternary Science Review, v. 28, p. 1200? 1216, doi:10.1016/j.quascirev.2008.12.022. 82 Johnson, S.Y., Nelson, A.R., Personius, S.F., Wells, R.E., Kelsey, H.M., Sherrod, B.L., Okumura, K., Koehler, R., Witter, R.C., Bradley, L.A., and Harding, D.J., 2004, Evidence for late Holocene earthquakes on the Utsalady Point Fault, northern Puget Lowland, Washington, Bulletin of the Seismological Society of America, v. 94, no. 6, p. 2299? 2316. Jones, M.A., 1996, Thickness of unconsolidated deposits in the Puget Sound Lowland, Washington and British Columbia, U.S. Geological Survey Water Resource Investigative Report, 94-4133. Kelsey, H.M., Sherrod, B.L., Blakely, R.J., and Haugerud, R.A., 2012, Holocene faulting in the Bellingham forearc basin: Upper-plate deformation at the northern end of the Cascadia subduction zone, Journal of Geophysical Research, v. 117, B03409, doi:10.1029/2011JB008816. Kinabo, B.D., 2007, Incipient continental rifting: insights from the Okavango Rift Zone [Ph.D. dissertation]: Rolla, University of Missouri, p. 11-68. Lapen, T.J., 2000, Geologic map of the Bellingham 1:100,000 quadrangle, Washington, Washington Division of Geology and Earth Resources, Open File Report 2000-5. Lees, J.M., and Crosson, R.S., 1990, Tomographic imaging of local earthquake delay times for 3-D velocity variation in western Washington, Journal of Geophysical Research, v. 95, p. 4763 ? 4776. Liberty, M.L., 2007, Seismic reflection imaging across the eastern portions of the Tacoma fault zone, U.S. Geological Survey Project Award Number: # 07HQGR0088. Mazzotti, S., Dragert, H., Hyndman, R., Miller, M.M., and Henton, J., 2002, GPS deformation in a region of high crustal seismicity: N. Cascadia forearc, Earth 83 Planetary Science Letters, v. 198, p. 41?48, doi:10.1016/S0012-821X(02)00520- 4. McCaffrey, R., Long, M.D., Goldfinger, C., Zwick, P.C., Nabelek, J.L., Johnson, C.K., and Smith, C., 2000, Rotation and plate locking at the southern Cascadia subduction zone, Geophysical Research Letters, v. 27, p. 3117?3120, doi:10.1029/2000GL011768. McCaffrey, R., Qamar, A.I., King, R.W., Wells, R., Khazaradze, G., Williams, C.A., Stevens, C.W., Vollick, J.J., and Zwick, P.C., 2007, Fault locking, block rotation and crustal deformation in the Pacific Northwest, Geophysical Journal International, v. 169, no. 3, p. 1315?1340, doi:10.1111/j.1365- 246X.2007.03371.x. Oasis Montaj, 2012, Oasis montaj how-to guide: Oasis montaj advanced imaging: http://updates.geosoft.com/downloads/files/how-to- guides/Oasis_montaj_Advanced_Imaging.pdf (accessed January 2013). Oasis Montaj, 2009, Geophysics: montaj Extension developed by Geosoft: http://www.geosoft.com/media/uploads/resources/brochures/OM_G_fs_2009_10_ web.pdf (accessed January 2013). Phillips, J.D., Duval, J.S., and Ambroziak, R.A., 1993, National geophysical data grids; gamma ray, gravity, magnetic, and topographic data for the conterminous United States, U.S. Geological Survey Digital Data Series DDS-09, 1 CD-ROM. Pratt, T.L., Johnson, S., Potter, C., Stephenson, W., and Finn, C., 1997, Seismic reflection images beneath Puget Sound, western Washington State: The Puget lowland thrust sheet hypothesis, Journal of Geophysical Research, v. 102, no. B12, p. 27,469? 27,489. 84 Sherrod, B.L., Brocher, T.M., Weaver, C.S., Bucknam, R.C., Blakely, R.J., Kelsey, H.M., Nelson, A.R., and Haugerud, R.A., 2004, Holocene fault scarps near Tacoma, Washington, Geology, v. 32, p. 9 ? 12. Sherrod, B.L., Blakely, R.J., Weaver, C.S., Kelsey, H.M., Barnett, E., Liberty, L., Meagher, K.L., and Pape, K., 2008, Finding concealed active faults: Extending the southern Whidbey Island fault across the Puget Lowland, Washington, Journal of Geophysical Research, v. 113, B05313. Tabor, R.W., Frizzell, V.A., Jr., Booth, D.B., and Waitt, R.B., 2000, Geologic map of the Snoqualmie Pass 30 ? 60? quadrangle, Washington, U.S. Geological Survey Geologic Investigations Series Map I2538, scale1:100,000. Telford, W.M., Geldart, L.P., and Sheriff, R.E., 1990, Applied Geophysics: Second Edition, New York, Cambridge Univ. Press, p.105-109. Tre?hu, A.M., Asudeh, I., Brocher, T.M., Luetgert, J., Mooney, W.D., Nabelek, J.L., and Nakumura, Y., 1994, Crustal architecture of the Cascadia forearc, Science, v. 266, p. 237 ? 243. van Wagoner, T.M., Crosson, R.S., Creager, K.C., Medema, G.F., Preston, L.A., Symons, N.P., and Brocher, T.M., 2002, Crustal structure and relocated earthquakes in the Puget Lowland, Washington from high resolution seismic tomography, Journal of Geophysical Research, v. 107 (B12), 2381, doi 10.1029/2001JB000710. Wells, R.E., Weaver, C.S., and Blakely, R.J., 1998, Forearc migration in Cascadia and its neotectonic significance, Geology, v. 26, p.759-762. Yount, J.C., Dembroff, G.R., and Barats, G.M, 1985, Map showing depth to bedrock in the Seattle 300 by 600 quadrangle, Washington, scale 1:100000, U.S. Geological Survey Miscellaneous Field Study Map, MF-1692. 85 Zoeback, M.L., 1992, First- and second-order patterns of stress in the lithosphere: The world stress map project, Journal of Geophysical Research, v. 97, p. 11,703- 11,728. 86 Appendix Covington gravity base station Covington base station (red dot) located at the main (south) entrance of the City of Covington building near the flag poles. State: Washington City: Covington Building: City of Covington building Elevation: 117.746 m Coordinates: 47.3590427322?, -122.118273863? Gravity: 980723.561