FINDING HIGH GROUND: SIMULATING AN EVACUATION IN A LAHAR RISK ZONE by JOSEPH A. BARD A THESIS Presented to the Department of Geography and the Graduate School of the University of Oregon in partial fulfillment of the requirements for the degree of Master of Science June 2016 !ii THESIS APPROVAL PAGE Student: Joseph A. Bard Title: Finding High Ground: Simulating an Evacuation in a Lahar Risk Zone This thesis has been accepted and approved in partial fulfillment of the requirements for the Master of Science degree in the Department of Geography by: Dr. Christopher Bone Chairperson Dr. Hedda Schmidtke Member and Scott L. Pratt Dean of the Graduate School Original approval signatures are on file with the University of Oregon Graduate School. Degree awarded June 2016 !iii © 2016 Joseph A. Bard !iv THESIS ABSTRACT Joseph A. Bard Master of Science Department of Geography June 2016 Title: Finding High Ground: Simulating an Evacuation in a Lahar Risk Zone Large lahars threaten communities living near volcanoes all over the world. Evacuations are a critical strategy for reducing vulnerability and mitigating a disaster. Hazard perceptions, transportation infrastructure, and transportation mode choice are all important factors in determining the effectiveness of an evacuation. This research explores the effects of population, whether individuals drive or walk, response time, and exit closures on an evacuation in a community threatened by a large lahar originating on Mount Rainier, Washington. An agent-based model employing a co-evolutionary learning algorithm is used to simulate a vehicular evacuation. Clearance times increase when the population is larger and when exits are blocked. Clearance times are reduced when a larger proportion of agents opt out of driving, and as the model learns. Results indicate evacuation times vary greatly due to spatial differences in the transportation network, the initial population distribution, and individual behaviors during the evacuation. !v CURRICULUM VITAE NAME OF AUTHOR: Joseph A. Bard GRADUATE AND UNDERGRADUATE SCHOOLS ATTENDED: University of Oregon, Eugene, Oregon Portland State University, Portland, Oregon DEGREES AWARDED: Master of Science, Geography, 2016, University of Oregon Bachelor of Science, Geography, 2014, Portland State University !vi ACKNOWLEDGMENTS I wish to express sincere appreciation to Dr. Christopher Bone for his assistance in the preparation of this manuscript. In addition, special thanks are due to my other committee member, Dr. Hedda Schmidtke, for offering her time and suggestions to improve this effort. I would also like to acknowledge Dr. Nick Kohler, Megen Brittell, Nick Perdue, Dr. Mike Nelson, and Maya Havasova whose technical support was invaluable to this study. Lastly, thank you to all the UO Geography grads, especially Christine Grummon, Christina Shintani, Jewell Bohlinger, Zach Thill, Anna Moore, and Christine Carolan for your humor, friendship, and support throughout this journey. !vii This thesis is dedicated to my parents, Abby and Matthew Bard, for their universal support of my endeavors (no matter how ridiculous they must have seemed at the time!), to Jeanne Turner for her love, support, and partnership, to all my friends who make my life so great, and lastly, to Peter Menuez and Dr. Doug Foster–great people who left us too soon. !viii TABLE OF CONTENTS Chapter Page I. INTRODUCTION AND LITERATURE REVIEW ................................................ 1 II. METHODS .............................................................................................................. 11 2.1 Study Site ......................................................................................................... 11 2.2 Description of Model ....................................................................................... 13 2.3 Purpose ............................................................................................................. 13 2.4 Entities, State Variables, and Scales ................................................................ 14 2.5 Process Overview and Scheduling ................................................................... 16 2.6 Design Concepts .............................................................................................. 17 2.6.1 Basic Principles ....................................................................................... 17 2.6.2 Emergence ............................................................................................... 18 2.6.3 Objectives ............................................................................................... 18 2.6.4 Learning .................................................................................................. 18 2.6.5 Sensing .................................................................................................... 18 2.6.6 Interactions .............................................................................................. 19 2.6.7 Stochasticity ............................................................................................ 19 2.6.8 Observation ............................................................................................. 19 2.6.9 Initialization ............................................................................................ 20 2.7 Input Data ......................................................................................................... 20 2.7.1 Agents ..................................................................................................... 20 2.7.2 Network ................................................................................................... 21 !ix 2.8 Submodels ........................................................................................................ 21 2.8.1 Creating a Set of Agents for a Model Run .............................................. 21 2.8.2 Routing Strategy ..................................................................................... 22 2.8.3 Re-routing Strategy ................................................................................. 23 2.8.4 Determining Pedestrian Evacuation Potential ........................................ 23 III. RESULTS .............................................................................................................. 25 3.1 ANOVA Results .............................................................................................. 26 3.2 Population Size ................................................................................................ 28 3.3 Proportion of Agents Opting Not to Drive ...................................................... 30 3.4 Response Time to the Call to Evacuate ........................................................... 31 3.5 Exit Point Closure ............................................................................................ 33 3.6 Comparisons of Clearance Times by Evacuation Mode .................................. 34 IV. DISCUSSION ........................................................................................................ 40 4.1 Population: The Effects of More Agents in the System .................................. 41 4.2 Walking Versus Driving: The effect of distance from safety on evacuation potential ............................................................................................... 41 4.3 Response time: The Effect of When Agents Start to Evacuate ....................... 43 4.4 Network Closures: The Effect of Having One Less Way Out ......................... 43 4.5 Model Learning: How Should its Effects Be Interpreted? ............................... 44 4.6 Directions for Future Work .............................................................................. 45 4.7 Conclusion ....................................................................................................... 46 REFERENCES CITED ................................................................................................ 48 !x LIST OF FIGURES Figure Page 1. Map of historic lahar events originating on Mount Rainier ................................... 3 2. Mount Rainier simplified hazard map ................................................................... 4 3. Maps of study area showing number of households and walking time to safe zone ............................................................................................................ 12 4. Overview of the MATSim model run process ....................................................... 16 5. Summary statistics of clearance times ................................................................... 26 6. Cumulative clearance time as a function of population ......................................... 29 7. Cumulative clearance time as a function of population ......................................... 29 8. Cumulative clearance times as a function of proportion of agents opting not to drive ............................................................................................................. 30 9. Clearance times as a function of proportion of agents opting to walk ................... 31 10. Sample Poisson distributions of agent response time parameter ........................... 32 11. Clearance times as a function of response time and population ............................ 32 12. Clearance times in response to exit closures .......................................................... 33 13. Comparing clearance times for walking versus driving, iteration 1 ...................... 35 14. Comparing clearance times for walking versus driving, iteration 20 .................... 36 15. Comparing ideal drive times versus actual drive times, iteration 1 ....................... 38 16. Comparing ideal drive times versus actual drive times, iteration 20 ..................... 39 !xi LIST OF TABLES Table Page 1. Model initialization parameters for parameter sweep ............................................ 20. 2. ANOVA results for iteration 1 ............................................................................... 27 3. ANOVA results for iteration 20 ............................................................................. 28 ! 1 CHAPTER I INTRODUCTION AND LITERATURE REVIEW Volcanogenic mud and debris flows, known by the term of Javanese origin, lahars, are mixtures of water, ash, tephra, and rock stripped from volcanic slopes with a texture of wet cement to thick motor oil (Pierson, Wood, & Driedger, 2014). Large lahars have the power to move boulders up to 10 meters in diameter, entraining wood and other debris as they surge downward through valley drainages. For communities living near volcanoes, lahars pose a significant threat to property and life. In the 20th century alone, lahars were responsible for injuring 5,022 people and the deaths of 29,937 more –12.5 % of the deaths attributed to volcanoes over this time (Witham, 2005). The most catastrophic incident involving a lahar occurred on the evening of November 13, 1985. Following a yearlong period of unrest and anomalous activity, Nevado Del Ruiz, the northernmost volcano in Colombia’s Andes Volcanic Chain, began to erupt. Magma and ash erupting from the volcano formed pyroclastic flows, high density mixtures of ash, fragments of rock, and volcanic gasses pulled down slope by the force of gravity (U.S. Geological Survey, 2016). The combination of extreme heat and physical abrasion scoured and melted the ice and snow on the glaciated summit, sending an immense volume of meltwater into six river valleys draining the mountain (Pierson, Janda, Thouret, & Borrero, 1990). The bulk of the meltwater entered the headwaters of Rio Lagunillas system, eroding sediments deposited by prior eruptions, and transforming into a lahar as it surged downslope (Scott & Vallance, 1995). Two hours after the eruption began, the now massive lahar had overrun the 45 kilometer reach between Volćan Ruiz and the riverside town of Armero. In the minutes following the arrival of the lahar’s first pulse, Armero was inundated: 23,080 people died and 4,470 more were injured (Pierson, Janda, Thouret, & Borrero, 1990; U.S. Geological Survey, 2009). Although Armero was situated upon sizable deposits from at least two previous lahars witnessed in historic times (indicating a precedent for similar events), the call to evacuate was never made. This tragedy was not an unfortunate outcome of best laid plans. Instead, this natural hazard was transformed into a human disaster because of interagency ! 2 bureaucracy, trepidation towards the potential political or financial fallout from calling a false alarm, and the difficulty of evacuating a large population (Voight, 1990; Scott, Macias, Naranjo, Rodriguez, & McGeehin, 2001). The lahar that interred Armero ranks as the 4th deadliest volcano disaster in history. The geography that made Armero vulnerable to this tragedy is not unique. Many cities around the world are built at the foot of a large and unstable volcanoes. In the Pacific Northwest of the United States, towering 4,392 m above sea level, Washington’s iconic Mount Rainier is an active stratovolcano capped by a massive volume of ice and snow. It has a history of producing enormous lahars (Figure 1), including the largest know lahar, the Osceola Mudflow (Crandall, 1971). Today, many communities living in the shadows of the Cascade Range’s tallest peak are vulnerable to a catastrophic lahar exposure (Figure 2). The hazard potential stemming from this volcano evokes comparisons to the Armero disaster because, in the last 10,000 years, at least six large lahars have coursed though the Puget Sound lowlands, recurring on average, every 500 to 1,000 years (Crandall, 1971; Scott & Vallance, 1995). If a similar lahar occurred today, it is almost certain that it would be highly consequential because of the number of vulnerable communities in these lowland areas and the uncertain nature of the hazard (Pierson, Wood, & Driedger, 2014; Wood & Soulard, 2009; Chakraborty, Tobin, & Montz, 2005). Vulnerability is contingent on the attributes or properties of an individual or community and those of the physical hazard. This is expressed conceptually as the product of three interrelated factors: exposure, sensitivity, and adaptive capacity (Turner, 2003). Exposure is simply the physical intersection of people and property with the force of nature (Wood & Soulard, 2009). Sensitivity and adaptive capacity are multi-faceted combinations of the physical and social attributes of a community, an individual, or a system (Wood & Soulard, 2009). Sensitivity is the degree to which similarly exposed entities experience different adverse effects (Turner, 2003). For example, each community built within Mount Rainier’s lahar paths will be impacted differently based on the proportion of residents or vital assets within the hazard zones. Adaptive capacity is the ability to withstand an exposure through planning, or by virtue of the inherent ! 3 physical and social characteristics in place at the time of an event (Wood & Soulard, 2009). Figure 1: Map of historic lahar events originating on Mount Rainier. Three large events occurring in the last 10,000 years are shown (U.S. Geological Survey, 2015). ! 0 Figure 2: Mount Rainier simplified hazard map (USGS, 2014). ! 0 Lahars can be split into two categories based on how they are triggered. Primary lahars result directly from eruptive activity (e.g. hot ash melting snow and ice), and typically follow a period of volcanic unrest (Mothes & Vallance, 2015). Secondary lahars differ in that they are not directly related to an eruption (Mothes & Vallance, 2015). At Mount Rainier, some debris flows may be secondary type lahars produced by a sudden collapse of edifice bedrock (U.S. Geological Survey, 2014). Scott and Valance (1995) noted a “general lack of any association” between eruptive activity and Mount Rainier’s largest lahars. A large, secondary lahar occurring without warning is especially dangerous because the window of opportunity to evacuate is narrowed considerably. Whichever the root cause, evidence in the geologic record shows that large lahars originating on Mount Rainier’s slopes can travel significant distances from their source (Scott, Macias, Naranjo, Rodriguez, & McGeehin, 2001; Diefenbach, Wood, & Ewert, 2015). In the present context, nearly 10 % of Mount Rainier’s historic lahars would cause immense destruction to communities in the Puget Sound lowlands (Crandall, 1971; Scott, Macias, Naranjo, Rodriguez, & McGeehin, 2001). The last momentous event happened on Mount Rainier approximately 550 years ago when an estimated 0.23 km3 of the mountain’s hydrothermally weakened western flank collapsed, resulting in the Electron Mudflow (Scott & Vallance, 1995). Because it cannot be correlated to an eruption, it is suspected to be a secondary lahar that occurred spontaneously before traveling over 64 km and reaching the Puget Sound (U.S. Geological Survey, 2014). If a similar lahar occurred today, there would be no warning until the lahar was already flowing down stream. The city of Orting, 40 km from the mountain and built atop nearly 5 meters of rubble from the Electron Mudflow, would have only about 45 minutes of warning before the lahar’s arrival (Scott, Macias, Naranjo, Rodriguez, & McGeehin, 2001). An areogeophysical survey revealed a mass of hydrothermally altered rock measuring 1.6 km3 perched atop the western flank, corroborating the magnitude of the hazard present at Mount Rainier (Finn, Sisson, & Deszcz-Pan, 2001). Despite the current understanding of this region’s hazard severity, population within the hazard zones is growing (Pierson, Wood, & Driedger, 2014; Strader, Ashley, & ! 1 Walker, 2015). Across the four counties threatened by lahars from Mount Rainier, more than 78,000 people (as of 2009) reside in hazard zones (Wood & Soulard, 2009). Pierce County, which currently has the largest proportion of residents living in lahar hazard zones, is projected to grow by 180,000 residents between 2014-2030 (Diefenbach, Wood, & Ewert, 2015; Washington State Department of Transportation, 2015; Washington State Office of Financial Management, 2015). As the population grows in threatened areas, communities become increasingly vulnerable. Hazard vulnerability is a place-based phenomenon because it is a consequence of the unique combination of physical characteristics, social structures, institutional policies, and differential access to resources at the time of an exposure (Chakraborty, Tobin, & Montz, 2005). Faced by known threat, an individual or community can take action to reduce their vulnerability. Pre-event strategies for reducing exposure include hazard avoidance, hazard modification, and hazard warning systems (Pierson, Wood, & Driedger, 2014). In some instances, structures can be engineered to divert or deflect flows, yet these structures are no match for the sheer size of Mount Rainier’s largest lahars. Hazard avoidance, by limiting types of development in hazard zones, is the most effective strategy. But while conceptually simple, in practice it can be difficult achieve because it requires extended cooperation between the public and government (Pierson, Wood, & Driedger, 2014). For cities like Orting, this tension is highlighted by the fact that the interval between significant lahars can be many lifetimes and the land within the hazards zones is otherwise attractive. After all, no place is without its detracting risks. Even if a long-term avoidance strategy is undertaken, short-term strategies are still needed to address the more immediate threat. Efforts to mitigate sensitivities and improve adaptive capacity through planning and public education can reduce the severity of an exposure. In practice, an evacuation (coupled to a warning system) is an in-the-moment hazard avoidance strategy and a critical short-term action for reducing hazard exposure during volcanic unrest (Marzocchi & Woo, 2007). ! 2 For public officials, issuing an order to evacuate means weighing the risks of a false alarm against the threat to life and property. While choosing to evacuate regardless of the consequences resulting from a false alarm may seem like an easy decision, apprehension over losing institutional credibility and the possibility of litigation stemming from a false alarm can influence both how decision makers act, and where the public will turn for information (Dow & Cutter, 1998; Marzocchi & Woo, 2007). The common denominator underpinning all hazard mitigation strategies is that the extant knowledge about the hazard, and the tactics for minimizing its impacts, is shared amongst official public institutions and vulnerable populations (Pierson, Wood, & Driedger, 2014). To better inform all stakeholders about the dynamics that exert influence over an evacuation, location specific research is needed to explore how hazard perceptions, behavioral responses, and the underlying geography impact hazard zone clearance times. Spatial simulation modeling provides an approach to better understand the dynamics of an evacuation through explicitly defining and representing the rules by which individuals behave and the details of the transportation system hypothesized to regulate the system’s performance. All spatial models can be generalized into two conceptual categories, (1) aggregate and (2) disaggregate, based on how the entities being modeled are represented. Aggregate models are useful for describing the collective effects though general statistical properties, yet processes which govern outcomes cannot be incorporated. In contrast, disaggregate models are useful for learning about how environmental conditions and individual level behaviors influence process (Fothingham & Rogerson, 1993). One disaggregate simulation modeling approach, agent-based modeling (ABM), is advantageous for studying evacuations because each individual entity in the system is represented, which allows for each of their interactions with the environment and other agents to be represented as well. Through explicit representation of process, patterns that would otherwise be difficult to forecast can be studied. Important insights can be gleaned by asking what if questions designed to address specific details of a location or an underlying behavioral assumption. For this reason, ABMs are an important tool for scientists using the lens of complex systems science (CSS) to study system dynamics. ! 3 A core areas of study in CSS is “aggregate complexity” in which agent-to-agent interactions (within their environment) can form structures which may “exhibit learning and emergence” (Manson & O'Sullivan, 2006). Roadway congestion as a result of agent interactions is increasingly being viewed as an emergent outcome of a complex system–a novel, high-level system dynamic that cannot be obtained by simple aggregation of individual activities (Bonabeau, 2002). ABMs are attractive for simulating vehicle-based evacuations because they excel at revealing spatially and temporally explicit details of the transportation network and agent-to-agent interactions that impact network congestion and ultimately clearance times (Bernhardt, 2007). Recently, the transportation simulation package Multi Agent Transportation Simulation (MATSim) has gained favor for studying evacuations by practitioners in the transportation research field as well as in Geography. This software package has been selected for its ability to incorporate behavioral assumptions, agent re-planning capacity, and prior efficacy in researching scenario-based evacuations (Lämmel, Grether, & Nagel, 2010; Durst, Lämmel, & Klüpfel, 2012; Henry & Frazier, 2015). While others have used simulation modeling to explore various evacuation scenarios, there is much diversity in their approaches stemming from differences in environments, scales, hazards, and software packages. Critical features of the hazard or environment that shape the results in one study may not exist in another scenario. And thus, while similarities exist, the nature of each study makes it difficult to apply the lessons learned from one study area to another. For example, many evacuation studies are set in large cities with public transit systems (Shiwakoti, Liu, Hopkins, & Young, 2013). In addition, the hazard may have advanced warning, or the environment facilitates in-place sheltering, like in the Hamburg flooding scenario by Durst et al. (2012). Alternatively, the focus may be on the effect of a particular behavior, for instance, Liu and Murray-Tuite (2014) concentrated on family gathering prior to evacuation under threat from a hazardous waste release. A study may altogether forgo the specificity of a precise hazard scenario for a more general approach, as is the case for the studies of New Orleans, Louisiana (Naghawi & Wolshon, 2012) and Toronto, Ontario (Abdelgawad, Abdulhai, & Wahba, 2010). ! 4 Compared to other study areas, communities like Orting are distinct, and they need to be understood in their own ways. What sets Orting apart is its semi-rural landscape and the important differences between large lahars and other hazards. Outside of the small clusters of compact development, the unstructured open spaces surrounding Orting allow pedestrians to move without the constraints imposed by the urban form. Large lahars are different from phenomena with comparably destructive potential, both spatially and temporally. In contrast to the sudden onset, no-warning scenario, most other catastrophic scale hazards, like hurricanes and some tsunami waves, often can be detected hours to days prior to arrival. Further, the spatial extent of the lahar hazard zones is well known whereas these other hazards are subject to a higher degree of uncertainty surrounding the landfall point and extent of exposure. Evacuees in Orting (and similar cities without public transit options, and where a large proportion lives within walking distance to safety) will rely primarily on private vehicles or walking to evacuate. An additional challenge to the effectiveness of an evacuation is the rural road system which funnels traffic exiting the hazard zones through only a few exit points. A study of hazard perceptions in the Puyallup Valley (where 117 of the 257 respondents were Orting residents) revealed 55 % of the respondents held the opinion that official evacuation routes were inadequate, with the vast majority placing the blame on traffic congestion (Davis, et al., 2006). When asked which mode of transportation they were likely to use, 66.5 % of respondents indicated they planned to evacuate via car whereas 18.7 % intended to evacuate on foot (others indicated an alternative type of vehicle, like an RV, motorcycle, or bicycle, and a single respondent preferred horseback) (Davis, et al., 2006). Traffic congestion during an evacuation situation (and traffic congestion in general) can be thought of as the outcome of demand for roadway space overwhelming the supply of roads (Peeta & Hsu, 2009). Survey responses in Davis, et al. (2006) indicate that an evacuation in the study area may be affected, on the supply side of the equation, by the inherent structural limitations of the transportation network and, on the demand side, by evacuee mode choice. ! 5 The overall objective of this research is to simulate an evacuation of the city of Orting (and the surrounding areas) addressing how the in-place transportation infrastructure and evacuee mode choice impact the overall ability of the population to clear the hazard zone. Specifically, this research evaluates the degree to which four structural and behavioral characteristics impact hazard zone clearance times in the study area, (1) the total number of agents participating in the evacuation, (2) the proportion of the population (for whom walking to safety is feasible) opting to walk instead of drive, (3) the timeliness of agent response to the call to evacuate, and (4) the blockage of various points exiting the hazard zones. Agent-based modeling is used to simulate individual actors, behaving according to their own self-interests, who must negotiate a shared (and possibly overwhelmed) road network to reach safety prior to the arrival of a lahar. This approach will allow us to explore how impactful individual behaviors and the environmental conditions may be to the outcomes of an evacuation–both at the individual level and collectively across the system. This research has potential implications for policy makers and the public alike. The knowledge produced by simulation modeling is useful for focusing dialogue about the efficacy of vulnerability mitigation strategies at a variety of scales, and to inform actions may be useful today and in the future to reduce potential of loss of life and property. From this research, public officials may gain insight on where to focus public education and outreach programs, and which future public works may be most beneficial to the community. For the public, learning about how and why an evacuation is affected by certain actions could help to shift perceptions and motivate individuals to become more educated and prepared for an event. In all, the conversation around vulnerability can be framed in terms of how individuals can help themselves, and others within their communities, to deal with the uncertainty of living in an area where large and destructive hazards are a part of life. ! 6 CHAPTER II METHODS 2.1 Study site The study area for this research is centered around Orting, Washington, a semi-rural city in Pierce County located within the Puyallup Valley lahar hazard zone. The five most active volcanoes in the Cascades are in the state of Washington, and all have the potential to generate large lahars owing to their ice-capped summits and steep slopes. Over 191,000 residents live within these hazard zones, an area containing important economic hubs like the Port of Tacoma (Diefenbach, Wood, & Ewert, 2015). Mount Rainier is responsible for the largest proportion of threatened communities. Each of the twenty- seven communities within Rainier’s lahar paths is partially or entirely built in a hazard zone (Wood & Soulard, 2009). In 1956 Crandall and Waldron (1956) first provided evidence that this region has a history of large lahars when they reinterpreted the origin of a geologic unit covering 549 km2 called the Osceola till. Until this point it was considered be the product of Pleistocene Era glaciation, when in actuality, it was the deposit of an enormous mudflow originating on Mount Rainier (U.S. Geological Survey, 2014). Crandall’s (1971) report is the seminal work responsible for revealing that large lahars in these places are not especially uncommon, and first raising concern about community vulnerability in the Puget Sound lowlands. By 2022, following the regional trend of increasing population, Orting is expected to grow to nearly 8,000 residents, an increase of nearly 13 % over the 2013 total (City of Orting, 2015). During an evacuation each place must contend with its unique geography. Orting, for example is situated on a narrow strip of land between the Puyallup and Carbon rivers measuring between 800-1,800 m wide (Figure 3). Most of the roads leading out of the hazard zone must first cross a bridge over the Puyallup River before ascending the steep valley walls. In the event of a no-notice onset, secondary lahar, the lahar would likely reach Orting with only approximately 45 minutes of warning (Scott, Macias, Naranjo, Rodriguez, & McGeehin, 2001). ! 7 Figure 3: Maps of study area showing number of households and walking time to safe zone. ! 8 The location of the threshold between the hazard and safe zones is known with a high degree of certainty because a lahar traveling down the Puyallup Valley will be constrained within the steep valley walls. Thus, the gradient of exposure is relatively small, meaning that areas just outside the hazard zone will provide refuge for evacuees, especially when compared to a hazard like a hurricane, where the range of exposure severity is distributed over a much larger spatial gradient. In Orting, and many other Puget Sound lowland communities, pedestrian evacuation to high ground is a likely to be a viable option for some proportion of the population, and an important part of a comprehensive evacuation plan. 2.2 Description of model The model used in this research is described here using the Overview, Design, and Details (ODD) protocol developed by Grimm, et al. (2006; 2010). This protocol is commonly employed in agent-based modeling research for providing a consistent manner to communicate agent-based models. The protocol is meant to be a systematic review of the model, at first providing a high-level explanation, with each following section providing more in-depth details than the previous. First, an overview is given with information about the purpose and objectives of the model, details about the study area, model parametrization and scheduling. Next, the design concepts are discussed explaining the conceptual and theoretical underpinnings of the system and entities being studied. Finally, details are given about model initialization, input data, and the submodels. 2.3 Purpose The purpose of this model is to simulate an evacuation of a small semi-rural population threatened by a large lahar to explain how hazard zone clearance times are impacted by evacuation mode choice and the number of agents in the system during an evacuation. The ABM used in this study, MATSim 0.7.0 (2015), is appropriate because individually representing each agent and their behavioral attributes, such as evacuation mode choice and response times, permits a more nuanced assessment of how bottom-up processes, like agent-to-agent interactions, lead to emergent patterns of traffic congestion ! 9 which ultimately impact hazard-zone clearance times. As summarized by Shiwakoti, et al. (2013), other works have employed micro-simulation ABMs to study single and multi- modal evacuations under various hazards contexts. However, in these experiments evacuees are assumed to choose a single common evacuation mode (such as on foot or by automobile), or when more than one mode is is available to the population, the split is between automobiles and mass transit or pedestrians and mass transit. Additionally, these studies are situated in large cities rather than rural areas. Large lahars, too, have unique temporal and spatial characteristics differentiating them from other types of hazards in which ABMs have been used, thus necessitating research focused on this specific hazard and in the context where pedestrian evacuation is a viable alternative vehicular evacuation. This work also incorporates least-cost distance path modeling of pedestrian evacuation potential (Wood & Schmidtlein, 2012) as a method to parameterize the the ABM by identifying the proportion of population for whom walking is an alternative mode of evacuation to driving (see section 2.7.4). 2.4 Entities, state variables, and scales There are two main entities in the model: car agents (referred to hereafter as agents or cars) and the road network (hereafter referred to as roads or network). Each agent is an inclusive entity incorporating all vehicle types into a single representation (e.g. personal vehicles, commercially employed vehicles, etc.). The occupancy of each vehicle is not represented because the research is aimed at understanding how varying levels of vehicular demand put upon the network impact clearance times rather than attempting to predict the specific number of people able to evacuate. At the initiation of the model, each agent has a pre-determined activity plan: simply to exit the hazard zone and travel to a single destination in the safe zone. However, information about which route to take is not part of this original plan. Agents have a memory which can store up to 3 plans holding routes as chains of activities occurring at network locations (e.g. entering or leaving a link is considered an activity) as well as information about the start time and duration of each activity (Horni, Nagel, & Axhausen, 2016, p. 4). Beyond this, cars do not have any state variables that update during model ! 10 runs. All agent activities on the network are recorded in an event file summarizing each run. Further information about the interaction between the event file and the agent’s plan will be discussed in Section 2.5. The network is a vector representation comprised of nodes and links. An important detail of the MATSim representation is that all agent activities take place on the links of the network rather than on the nodes. Each node has one attribute variable: an X,Y coordinate pair defining its location in accordance with the spatial reference system employed by the model (here NAD83 UTM 10N). In this representation nodes exist to simply serve as junction points for the links and do not affect agent travel or network characteristics. Links have three variables affecting the simulation: length, lanes (number of roadway lanes), and a free-speed variable limiting the maximum unimpeded travel speed of the cars. Here, all the variables assigned to the links do not change from their initial values. The network has two types of special links that tie the edges of the network together so all agents can be routed to a single point of attraction. The first type is a single “super- link” serving as the final destination for all agents, this effectively has infinite length and free-speed parameters. The second type of link also has an infinite speed parameter, but a length of 1 meter. These links tie the edges of the network to the super-link so that all network edges attract agents equally. While the special links are not representative of world road networks, they are allowed in the directed-graph computational framework of MATSim. This is a critical factor for calculating routes, and most importantly, the single point of attraction eliminates a top-down, deterministic origin-to-destination routing assignments. The spatial and temporal resolution (the duration of each time-step) of the ABM is one meter and is one second, respectively (MATSim default settings). The network variable length is given in meters and the free-speed rate is given in meters per second. The network is a fixed set of links and nodes; however, five variations of the network are used in the ABM. In the default configuration, all four exit links have typical link parameters. The four network variants each have one of the four exit links effectively ! 11 blocked by setting the free-speed = 0.001 m/s. The exits are roughly located in the northwest, northeast, southwest, and southeast quadrants of the study area. 2.5 Process overview and scheduling The ABM consists of the processes depicted in Figure 4 (Horni, Nagel, & Axhausen, 2016). Each MATSim simulation run is comprised of a pre-determined number of iterations of the initial demand, mobsim, scoring, and replanning cycle–the analysis step takes place after a run is completed, and is not relevant to this section. A set of agents is chosen in the initial demand step (explained in Submodels 2.8.1), agent travel behavior is simulated in the mobsim step, next each agent’s travel is scored, last, during the re- planning phase, agents adjust their behaviors aiming to improve their score in future iterations. Figure 4: Overview of the MATSim model run process (Horni, Nagel, & Axhausen, 2016). Each MATsim iteration is designed to capture agent travel behavior for the duration of the agent’s activities (up to 24 hours) with the resolution of 1 second time steps. By default, the start time of a model iteration is the first scheduled activity. Once an iteration has started, information about the model iteration is recorded in the event file. In turn, the effective end of the iteration is the moment when all agents have completed their chain of activities and the last record for each agent’s activity chain is written to the event file. Although the model is updated every time step, records are not necessarily generated for every time step. Instead, records including time step and event information are written as ! 12 events take place, such as the times that an agent enters or exits a link. During the first iteration all agents attempt to complete their chain of activities via the shortest path as calculated by Dijkstra’s algorithm (Dijkstra, 1959) (explained in Submodels 2.8.2). Following the mobsim step are the scoring and replanning steps. Plan scores, for each agent, are calculated by comparing the measured time intervals needed to complete their chain of activities to the expected time (Charypar & Nagel, 2005). The plan and associated score from the initial iteration is stored in the agent’s memory and serves as the baseline for comparison to plans developed in future iterations. Plans for the next two iterations are added to the agent’s memory until three plans are stored, ranked from highest score to lowest. If plans with higher scores are developed in subsequent iterations, these plans will replace lower scoring plans in the agent’s memory. After scoring, the replanning strategy is employed; 10 % of agents choose a new route (explained in Submodels 2.8.3) while the other 90 % choose a plan with the highest utility score. 2.6 Design concepts 2.6.1 Basic principles The basic principle of this ABM is incorporating behaviors and preferences, such as evacuation mode and response time into a no-notice evacuation simulation. The no-notice scenario is important because the study area has a record of being subjected to catastrophic scale lahars that may occur spontaneously and without warning. The last large lahar in the study area was prior to historic times and thus there is no institutional memory of such an event. In addition, the spatial extent of the hazard is constrained by the river valley such that the gradient between the hazard and safe zones is quite sharp, but within the hazard zone, the degree of devastation is expected to be substantial. The preference to evacuate by car versus to evacuate on foot is evaluated for its impact on system-wide clearance times. Pedestrian movements are not explicitly modeled here, but vehicles that would otherwise participate in the evacuation are, instead, excluded from the simulation in the specific areas where walking to safely is plausible. The hypothesis underpinning this study is that when the travel demand outpaces the finite supply of ! 13 usable roadways, the model will produce an emergent patterns of network latency (traffic congestion) ultimately impacting the the ability of the at-risk population to clear the hazard zone. In contrast, eliminating vehicles in the zones where walking is possible will reduce demand at critical points in space (nearest to the exists) minimizing patterns of congestion, facilitating egress for agents whose sole option is driving. 2.6.2 Emergence The emergent outcomes of this ABM are the patterns of traffic congestion that form when road capacity is overwhelmed by the demand for roadway space. 2.6.3 Objectives The objective of each agent is to complete their chain of activities by selecting the most efficient route. The overall objective is for all agents to clear the hazard zone as quickly as possible, prior to the arrival of the hazard. 2.6.4 Learning Agents learn individually and collectively using a co-evolutionary algorithm. This process is aimed at optimizing the behavior across the system to produce a system-wide state in which no agent can improve their outcome. Each agent maintains a set of plans, each with an associated score, ranking the plan’s fitness compared to an idealized time to complete the plan’s chain of activities. As explained earlier, during the replanning stage, a number of stochastically selected agents modify existing plans during the next iteration. Although this takes place at an individual agent level, the effect is distributed collectively across the set of agents, until no agent can improve their score by acting unilaterally (Horni, Nagel, & Axhausen, 2016, p. 8). 2.6.5 Sensing Agents in this model sense road conditions as well as the presence of other agents sharing the network. Agents have an awareness of the maximum travel speed across each link and the link’s current and maximum capacity (capacity is a function of link length, free speed, and the number of agents currently occupying the link). Agents also sense the occurrence of a disturbance which motivates them to begin evacuating. The disturbance ! 14 in this ABM is the call to evacuate which is prescribed in each agent’s plan. So, while the network conditions are sensed dynamically, the call to evacuate is effectively sensed by proxy. 2.6.6 Interactions Agent-to-agent interactions are indirect. Because each network link has a finite capacity, queues form at the entry points of full links. Agents are excluded from entering full links until the queue has dispersed. Agents can sense the whether a link can be entered, but do not directly sense the presence of other agents. 2.6.7 Stochasticity Stochasticity is present during model iterations and the population selection process. While the model is running, in the replanning step (Figure 4), 10 % of the population uses a rerouting strategy to improve upon their plan score for the next iteration, while the remaining proportion re-uses their highest scoring plan. Selecting the population for each model run is also a stochastic process using the Python random.sample function (Python Software Foundation, 2016). During the population creation process each agent is assigned a time to enter the network after the call to evacuate. The time (in minutes) is chosen by drawing from a set of integers generated by the Python numpy.random.poisson function (SciPy.org, 2015). 2.6.8 Observation Every 10th iteration of the model, an XML event file is created containing records detailing every agent action, with its associated time step, throughout the duration of the iteration. Actions detailed in the event file are the agent’s initial entrance onto the network, each time an agent enters or exits a link, and the agent’s arrival at their final destination. The event file can be parsed to reveal information about individual agents or the population in aggregate. 2.6.9 Initialization The initial state of the ABM at time t = 0 corresponds to the moment an evacuation signal is given. As such, t = 0 is an arbitrary time that doesn’t represent a real time of day. ! 15 It is assumed the entire population senses the signal to evacuate simultaneously, but agents begin evacuating according to the response time parameter. The set of model runs in this study was initialized by running a parameters sweep. Each model run is one of four-hundred and eighty unique combinations of the parameters and their associated values listed in Table 1. At the beginning of each iteration pt(n) agents are distributed randomly throughout each census block in the study area with respect to the proportion of households in each census block to the total number of households in the study area. When an iteration is initialized, MATSim assigns each agent to the nearest network link based on the agent’s X,Y coordinate. Table 1. Model initialization parameters for parameter sweep. Population*number*(pt)* 1000! 2000! 3000! 4000! 5000! 6000! Proportion*of*walkers*(pw)* 0! 0.33! 0.66! 1.0! ! ! Response*time*(rt),*minutes* 0! 3! 6! 9! ! ! Network*(rn),*exit*closed* All!exits!open! Northwest! Northeast! Southwest! Southeast! ! ! 2.7 Input data 2.7.1 Agents The data for parameterizing agents in this study is based on the 2010 US Census count of households per census block. The the data, an ERSI polygon shapefile, was downloaded from the Pierce County (Washington) Open Geospatial Data Portal (2015). The complete set of agents was created using the ESRI ArcGIS 10.3 Random Points tool to allocate 9995 agents to 166 census blocks proportionally according to the number of households per census block to the total households in the study area. Additionally, agents are designated as potential walkers (those that are within a 40-minute walk from ! 16 high ground) and drivers (those located beyond the 40-minute walking threshold). The method for distinguishing these populations is explained below in Submodels 2.8.4. 2.7.2 Network The network used in this ABM was sourced from an ERSI polyline shapefile of roads downloaded from Pierce County, Washington (Pierce County Open GeoSpatial Data Portal, 2016). The total set of roads in the network includes all roads within the study areas hazard zone plus the roads extending 5 km of network distance outside the hazard zone measured from each hazard zone exit point. Service areas were determined using the ESRI ArcGIS Network Toolset’s service area generation function. Where the roads beyond the exit point did not extend a full 5 km, links were added to increase the network length in these areas to ensure the balance of attraction to all network edges. 2.8 Submodels 2.8.1 Creating a set of agents for a model run Creating the set of agents for each model is a two step process, (1) a Bash script iteratively loops through all possible population parameters for determining the agent population for each model run, (2) a combination of the parameters is passed as a set of input variables for a Python script. The Python script takes these inputs and writes the plans.xml file, the population input for a MATSim model run. The Bash script uses a nested for-loop to iterate over the population parameters: population number pt(n) = {1000, 2000, 3000, 4000, 5000, 6000}; the proportion of agents from the population of potential walkers who are assumed to walk rather than evacuate via car pw(proportion) = {0, 0.33, 0.66, 1}; and the Poisson distribution lambda parameter for assigning agent response time rt(minutes) = {0, 3, 6, 9}. Each step in the for-loop passes 1 of 96 unique combinations of these parameters to the Python script. Writing the plans.xml file is a two step routine: First, the population pt(n) and population proportion pw(proportion) parameters are used to select a subset of agents for the model run from the total possible set of agents. Agents are selected from both the driving only population and the potential walker population based on the fixed ratio of ! 17 these populations in the total set of possible agents, 3558:6437. However, as the proportion of potential walker agents who chose to walk increases, fewer agents are included in the subset of agents participating in the simulation run. For example, if the pt(n) = 1000 and pw(proportion) = 0, 355 agents would be selected from the walker/driver population and 644 would be drivers only. If proportion of walkers in increased to pw(proportion) = 0.33, 117 of the 355 agents from the walker/driver population would be excluded from the simulation run. The second step is assigning a response time for each agent by either setting all agents to begin evacuating simultaneously at the call to evacuate rt(minutes) = 0, or by drawing from a random set of Poisson distributed integers with the lambda parameter varied from rt(minutes) = {3, 6, 9} representing the number of minutes it takes for an agent to begin evacuating after the evacuation is initiated. 2.8.2 Routing strategy Routing in MATSim is solved using Dijkstra’s shortest path algorithm which determines the route between two network locations based on the aggregate least-cost path. Cost to traverse a link in a MATSim network is measured in seconds as a function of length/free-speed rather than linear distance. In the MATSim network all links have a fixed length and free-speed and are joined by nodes with fixed locations. These properties mean a route can be calculated using Dijkstra’s algorithm between any two nodes in the network. The basic principle behind this algorithm is that between any two nodes in the network, there is one set of (all possible combinations) of links, that when combined is the least costly of all paths (here measured in seconds). The algorithm starts at a given origin node and iteratively searches the links it is connected to for the link with the shortest length. The shortest link becomes the first segment of the route. In further iterations, this process is repeated again for the links joining the nodes one step from the origin node, and so on for all connected links. At each iterative step, one path is shorter than all others. This path is saved, while the others are discarded until the destination is reached via the shortest path (Dijkstra, 1959; Horni, Nagel, & Axhausen, 2016). ! 18 2.8.3 Re-routing strategy For each iteration, 10 % of the agent population modifies their route from the previous iteration. Modifications are based on Dijkstra’s shortest path algorithm comparing the cost of time; however, link-time values are amplified based on degree of traffic congestion experienced on the link in the previous iteration (Horni, Nagel, & Axhausen, 2016, p. 42; Lefebvre, Balmer, & Axhausen, 2007). 2.8.4 Determining pedestrian evacuation potential To determine which agents have the potential to walk versus drive, a least-cost distance (LCD) anisotropic path distance analysis was performed using ESRI ArcGIS, version 10.3 (Wood & Schmidtlein, 2012). The basic principle of this method is that, for a given walking speed across a hypothetical surface, the type of surface, topography, and whether the direction of travel is up or down slope will modify ideal travel speed (travel speeds are reduced in all cases except for slight downhill travel across paved surfaces). The model returns a raster grid where the value of each cell is the time needed to reach the nearest hazard zone exit point. The inputs for this model are a Lidar derived, digital elevation model (Puget Sound LiDAR Consortium, 2004) and land cover raster from the National Land Cover Database (NLCD, 2011) resampled to match the 2 m resolution elevation data. This method’s advantages are that it can be applied if high-resolution elevation and land cover data are available, and the results are easy to understand. An important methodological limitation is that interactive processes are not integrated into the model and information like the degree of congestion at key egress points can only be estimated (Wood & Schmidtlein, 2012). The authors note that results of LCD models should be considered feasibility baselines rather than expectations for individual evacuees. A lahar large enough to trigger the detection system (USGS, 2014) is estimated to reach the up-river extent of the study in 30-60 minutes and the center of the study area within 60 minutes (Pierce County Open GeoSpatial Data Portal, 2006). The time of 40 minutes was selected as the threshold for walking potential. As such, an agent whose initial location is coincident with a raster cell where t <= 40 minutes is designated as a ! 19 potential walker, whereas agents initially located beyond the 40-minute threshold are driver only agents. ! 20 CHAPTER III RESULTS Overall, two broad trends surfaced from the models results: (1) Clearance times increased with an increase in the number of agents in the system and when network exit points were blocked, and (2) clearance times decreased in response to a larger proportion of agents opting to walk instead of drive and when the response time parameter value was increased. Additionally, clearance times were reduced as a product of the model’s iterative learning process. For each model run, the learning process had a consistent effect of reducing mean clearance times over each subsequent iteration. However, when clearance times are no longer decreasing, any additional iterations are redundant. The graphs of clearance time summary statistics, for all model runs combined (Figure 5), show the overall effect of the iterative learning process (results from iterations 1, 10, and 20 are shown). The greatest proportion of clearance time reduction took place between iterations 1 and 10. Between iterations 10 and 20, further reductions were minimal, indicating that iteration 20 represented the extent of the effect of model learning, and additional iterations were unnecessary (the influence of model learning is explained further in Section 4.5). Thus, for the remainder of this study, results from the iteration 10 are omitted. All analysis hereafter is performed using the results from iteration 1 and 20. To explain the results of this study in more detail this chapter is divided into subsections. The first section reports the results of a sensitivity analysis using the ANOVA test (performed in the R software package) to measure the variance of mean clearances times in response to changes of the input parameters. The following subsections are organized by the factors guiding the research questions. Each factor’s effect on the model’s results will be examined in terms of both its overall importance and whether interactions with other factors influence clearance times. ! 21 Figure 5: Summary statistics of clearance times. Between iterations 1-10 clearance times are reduced through the model learning process. Graphs for iterations 10-20 are virtually identical, indicating no further learning can take place and more iterations are unnecessary. 3.1 ANOVA results To measure the degree of clearance time variation produced by changes in the input parameters, two ANOVA tests of means were performed. The first used results from model iteration 1, and the second with results from iteration 20 because these two conditions represent two theoretical ends of the model learning scenarios. The two ANOVA tests revealed that the variance of mean clearance time produced by changes to each of the individual input parameter values resulted in high F statistics values and p-values below the p <= 0.001 significance threshold. For the model iteration 1 scenario the the input parameters were compared pairwise to evaluate the between- group variance resulting from changes to input parameters. For all combinations of population number pt, proportion of walkers pw, and road network rn, parameters statistically significant differences of were measured below the p = 0.001 level. In contrast, when the parameter for response time rt was included in the pairwise Iteration 1 Iteration 10 Iteration 20 0 50 100 150 200 250 300 350 400 450 500 2000 4000 6000 2000 4000 6000 2000 4000 6000 Population M inu te s maximum mean median sd minimum Summary Statistics of Clearance Times as a Function of Agent Population ! 22 comparisons, no significant differences in the variation between-groups was detected; p- values were greater than the p = 0.1 threshold for all combinations (Table 2). The ANOVA test of outputs from model iteration 20 produced similar results. Again, the effect of all parameters measured individually resulted in differences of within-group variance matching the levels of statistical significance in the iteration 1 scenario. The pairwise comparisons of between-group variance of means yielded similar results for groupings of population number pt, proportion of walkers pw, and road network rn parameters. Similarly, the pairwise comparisons including the response time rt parameter were not statistically significant, with one exception. The pair of population number pt and response time rt returned p-values below the p <= 0.001 level, indicating that the between-group means were significant (Table 3). Table 2. ANOVA results for iteration 1. Model*iteration*1* * Df* Sum*Sq* Mean*Sq* F*value* Pr(>F)* Population*(pt)* 1! 1.54E+09! 1.54E+09! 2.44E+05! <2e?16! Response*time*(rt)* 1! 2.97E+06! 2.97E+06! 4.70E+02! <2e?16! Proportion*of*walkers*(pw)* 1! 4.84E+08! 4.84E+08! 7.64E+04! <2e?16! Road*network*(rn)** 4! 5.34E+08! 1.34E+08! 2.11E+04! <2e?16! Population*(pt)**:**Response*time*(rt)* 1! 2.86E+03! 2.86E+03! 4.52E?01! 0.501! Population*(pt)*:*Proportion*walkers*(pw)* 1! 6.23E+07! 6.23E+07! 9.85E+03! <2e?16! Response*time*(rt)*:*Proportion*walkers*(pw)* 1! 7.07E+03! 7.07E+03! 1.12E+00! 0.291! Population*(pt)*:*Road*network*(rn)* 4! 6.22E+07! 1.55E+07! 2.46E+03! <2e?16! Response*time*(rt)*:*Road*network*(rn)* 4! 4.12E+03! 1.03E+03! 1.63E?01! 0.957! Proportion*walkers*(pw)*:*Road*network*(rn)* 4! 1.34E+07! 3.35E+06! 5.30E+02! <2e?16! ! 23 Table 3. ANOVA results for iteration 20. Model*iteration*20* * Df* Sum*Sq* Mean*Sq* F*value* Pr(>F)* Population*(pt)* 1! 4.76E+08! 4.76E+08! 2.90E+05! <2e?16! Response*time*(rt)* 1! 2.83E+06! 2.83E+06! 1.72E+03! <2e?16! Proportion*of*walkers*(pw)* 1! 2.11E+08! 2.11E+08! 1.28E+05! <2e?16! Road*network*(rn)** 4! 9.60E+07! 2.40E+07! 1.46E+04! <2e?16! Population*(pt)**:**Response*time*(rt)* 1! 2.52E+04! 2.52E+04! 1.53E+01! 9.01E?05! Population*(pt)*:*Proportion*walkers*(pw)* 1! 2.88E+07! 2.88E+07! 1.75E+04! <2e?16! Response*time*(rt)*:*Proportion*walkers*(pw)* 1! 9.17E+02! 9.17E+02! 5.57E?01! 0.455! Population*(pt)*:*Road*network*(rn)* 4! 1.13E+07! 2.82E+06! 1.71E+03! <2e?16! Response*time*(rt)*:*Road*network*(rn)* 4! 4.30E+03! 1.08E+03! 6.53E?01! 0.624! Proportion*walkers*(pw)*:*Road*network*(rn)* 4! 4.32E+06! 1.08E+06! 6.56E+02! <2e?16! 3.2 Population size Across all scenarios, an increase in the agent population resulted in increased clearance times. The cumulative clearance time graphs (Figure 6) show the effects of population increase produces a consistent effect across all model runs. The prevailing pattern (seen in Figure 6) resulting from increasing the number of agents, is the rate of clearance begins to slow sooner, following the exit of a smaller proportion of the total agents. This pattern is much more pronounced for iteration 1 than for iteration 20. The box plots of this same data (Figure 7) show the distribution of clearance times are skewed towards the high end of the data range and for agents above the 50th percentile, maximum exit times can be three times greater than median clearance times. As model runs progress toward the iteration 20, this pattern remains consistent but the overall ranges of these clearance times decrease. ! 24 Figure 6: Cumulative clearance time as a function of population. Figure 7: Cumulative clearance time as a function of population. 1000 2000 3000 4000 5000 6000 1000 2000 3000 4000 5000 6000 Iteration 1 Iteration 20 0 50 100 150 200 250 300 350 400 450 500 Minutes Po pu lat ion Clearance Time as a Function of Agent Population ! 25 3.3 Proportion of agents opting not to drive As the proportion of agents opting not to drive increases, the resulting effect is a strong reduction of clearance times in all scenarios. This produces an effect counter to an increase of the agent population. The cumulative clearance time graph (Figure 8) shows that an increase in the proportion of walkers reduces overall clearance times and also increases the rate of hazard zone clearance. Additionally, the point at which the the curve begins to flatten out (indicating a slowing rate of egress) occurs after a larger proportion of agents have already cleared the hazard zone. Box plots of the clearance times (Figure 9) show a consistent pattern of clearance time reduction as the proportion of walkers increases, but the data remains skewed toward the high end of the range. The proportion of agents opting out of driving has no discernable effect of changing this pattern. And again, as model iterations increase, the range of clearance times shrink and the effect of the proportion of walkers becomes less dramatic. Figure 8: Cumulative clearance times as a function of proportion of agents opting not to drive. ! 26 Figure 9: Clearance times as a function of proportion of agents opting to walk. 3.4 Response time to the call to evacuate Agent response times to the call to evacuate are assigned by making a draw from a Poisson distributed set of numbers according to the agent population and the lambda parameter (Figure 10). Increasing the lambda parameter for a given population size produces a distribution that is less skewed and more bell-shaped. The effect of drawing response time variables from the more bell-shaped distribution reduced mean clearance times to the degree that the effect was found to be statistically significant in the ANOVA test. However, these changes were small. The box plots of clearance times for model iteration 1 (Figure 11) show a slight trend of reduced times toward the center of the distribution. Yet, this doesn’t necessarily result in a lower maximum value. The lower end of the boxes shows, as the response time parameter is increased clearance times are reduced, but the change is slight. This pattern is similar for iteration 20, but less pronounced (not shown). ! 27 Figure 10: Sample Poisson distributions of agent response time parameter. Figure 11: Clearance times as a function of response time and population. ! 28 3.5 Exit point closure The simulations were conducted with five network conditions, one with no blocked exits, and each of the other four having one exit closed. The process of model learning had the largest effect here when compared to the other factors (Figure 12). The patterns of clearance times changed entirely between the first and twentieth iteration of the model, whereas for other factors the underlying patterns shifted or were made more pronounced. Hazard zone clearance times were lowest when all exit point were open for both the iteration 1 and iteration 20, this is as expected because more exit points means the network can accommodate the demand more easily. In relation to model learning, clearance times were affected the least by a closure of the southeast exit, whereas the northwest exit closure scenario had the largest change between iteration 1 to 20. For iteration 1, where agents are choosing the shortest path, the southeast exit has the most effect on clearance times. By iteration 20 though, the southwest exit closure was similar to the closure of the southeast exit. The box plots show that the various exit closure scenarios are responsible for many of the high outlier values seen first iteration plots of other factors. Figure 12: Clearance times in response to exit closures. ! 29 3.6 Comparisons of clearance times by evacuation mode To understand evacuation potential comparatively, each agent’s measured clearance times can be contrasted against the ideal drive times and the time needed to evacuate on foot (grouped into 10 minute intervals). Comparing measured clearance times (driving) versus the time needed to walk to the safe zone (Figure 13 and Figure 14), the first broad observation is that the two plots have a similar overall look. However, there are important differences. Model learning reduces all measured clearance times by nearly one third from model iteration 1 to 20, and at this point the graphed patterns are less distinct as agent route-choice becomes increasing heterogeneous. The most important pattern to recognize is that some agents with the lowest walking times have some of the highest measured evacuation times. The variability of actual drive times becomes greater for agents that are a further walk from safety, especially if the time to walk is greater than 30 minutes. Figure 13 and Figure 14 compare clearance times for walking versus driving (for iteration 1 and iteration 20, respectively). These graphs show the time needed to evacuate by walking (x-axis of individual graphs) compared to the time needed to drive (y-axis of individual graphs), conditioned by both an increase of the population parameter (rows of graphs) and an increasing proportion of agents opting not to drive (columns of graphs). Data is aggregated for all scenarios in which all network exits are open. Each data point represents one agent; blue points represent agents initially located less than a 40-minute walk to safety, pink points represent agents initially located more than a 40-minute walk to safety. In each individual graph, agents are grouped by the time needed to walk to safety (in 10 minute intervals). In the far right column of graphs, all agents with an option to drive vs. walk (agents with less than a 40-minute walk time to safety) choose to walk. Because these agents all opt not to drive, no driving data for these agents is reported. ! 30 Figure 13: Comparing clearance times for walking versus driving, iteration 1. ! 31 Figure 14: Comparing clearance times for walking versus driving, iteration 20. ! 32 Another comparison can be made between the ideal drive time needed by each agent to exit the hazard zone (if the route taken had no traffic congestion) versus the actual measured drive time (Figure 15 and Figure 16). The plots again look similar, but the effect of model learning lowers all actual drive times and makes the plotted patterns less distinct. These plots show that for all evacuees, no idealized route is longer than ten minutes or twenty minutes, for iterations 1 and 20 respectively. Again, short idealized drive times do not necessarily mean short actual clearance times. In fact, actual drive times vary greatly even when agents have identical ideal drive times to safety. When agents are parsed into groups by the time needed to walk to safety, the variability between actual and ideal drive times increases as the agent’s time to walk to safety increases. In all population scenarios, some agents who have the shortest ideal drive times end up with the longest actual drive times. Figure 15 and Figure 16 compare ideal drive times versus actual drive times (for iteration 1 and iteration 20, respectively). These graphs show the ideal time needed to evacuate by driving (x-axis of individual graphs) compared to the actual time needed to drive (y-axis of individual graphs), conditioned by an increase of the population parameter (rows of graphs). Agents are grouped by the time needed to walk to safety (columns of graphs). Data is aggregated for all scenarios in which all network exits are open. Each data point represents one agent; blue points represent agents initially located less than a 40-minute walk to safety, pink points represent agents initially located more than a 40-minute walk to safety. ! 33 Figure 15: Comparing ideal drive times versus actual drive times, iteration 1. ! 34 Figure 16: Comparing ideal drive times versus actual drive times, iteration 20. ! 35 CHAPTER IV DISCUSSION Despite the destructive power of lahars demonstrated by events in recent history, little work has been devoted to evacuation studies that address the specific nature of this hazard and the places where they may strike. The objective of this study is to begin addressing this vacancy in the literature while contributing more broadly to research efforts aiming to understand how individual level behaviors can impact system wide outcomes of an evacuations. A tenet of the complex systems approach is explaining how macro-scale, system wide patterns are generated by micro-scale processes. This approach is guided by suppositions about key behaviors and relationships between agents, which is a fundamental turn from a reductionist methodology employing an aggregated and simplified view of the behavioral components of a system (Manson & O'Sullivan, 2006). MATSim was selected for its ability to represent the processes of many simultaneously interacting agents that manifest as system wide patterns of traffic congestions. This phenomena can alternatively be described as “event driven” (Millington, O'Sullivan, & Perry, 2012) where the patterns generated from cumulative effect of agent-to-agent interactions within the environment cannot be captured unless process is represented. The results of the simulations reveal why a complex systems approach and modeling process is important to evacuation research. Clearance times at the individual level vary widely from agent to agent regardless of their distance from the safe zone when distance is assessed, using simple proximity measures, like network time- distance or path distance (derived from the LCD model). The remainder of this chapter will revisit each topic directing the research questions to explain, in more detail, how each of the topical elements influenced different aspects of the model results. The limitations of this study will be discussed along with how future research could be directed to improve similar efforts in the future. ! 36 4.1 Population: The effects of more agents in the system Varying the number of agents in the system had a predictable effect consistent with the expectations of this study: more agents in the system increased clearance times in all scenarios. This makes sense when placed in terms supply and demand. Traffic jams form when the finite roadway space (supply) is overwhelmed by the pressure of an increased number of agents vying for room to move (demand). A method of the complex systems approach is to vary model input parameter values with the aim of revealing whether the system will respond by exhibiting discernably different behavioral properties when an input parameter is pushed past some threshold value. The questions and analyses in this study were designed to detect outcomes of this nature, i.e. would the system exhibit one behavioral regime while the population was below a certain value, but when population surpassed this threshold value, would different regime of behavior emerge? Increasing the number of agents in the system, however, resulted in consistent and incremental increases in clearance times. The results did not reveal system-wide behavioral changes generated by increasing the number of agents beyond a hypothetical threshold value. This effect is attributed to the first-in-first-out queueing structure of the model (Lämmel, Grether, & Nagel, 2010). This method of modeling traffic flows means that agent-to-agent interactions are indirect. And while traffic congestion does result from agents attempting to occupy the same network space, the formation and diffusion of queues is determined simply by order of arrival at a network junction, and inter-agent behaviors that may exacerbate or relieve queues are ignored. Queue handling is one of the primary controls over system behavior in MATSim. There is no mechanism to incorporate the influence of feedbacks into the overall dynamics of the system. Revisiting this study with a software package incorporating feedbacks into traffic-following behaviors, or other aspects of en route agent-to-agent interactions, may reveal different patterns of clearance times. 4.2 Walking versus driving: The effect of distance from safety on evacuation potential A more nuanced understanding of clearance potential can be gained by comparing actual evacuation times with an agent’s initial location at the start of an evacuation. This ! 37 location can be viewed by two different measures of proximity to the safe zone, (1) the time needed to walk to safety (Figure 13 and Figure 14), and (2) the time needed to drive (Figure 15 and Figure 16). These comparisons show that short evacuation routes (measured in terms of ideal drive times or walking times) do not necessarily correspond to low vehicular evacuation times. The large variations in clearance times that occur despite spatial similarities between agents (in terms of ideal clearance times) points to the degree of uncertainty that is introduced by the uniqueness of this place. This is a finding that would not have been possible to forecast without a methodology explicitly representing agent-to-agent interactions within the constraints of the transportation network. While the ABM approach used in this study made it possible to discover the important findings above, one limitation is that pedestrian agents are not explicitly represented. Pedestrian movements were omitted from this study for two reasons, (1) the focus of this study is understanding the dynamics of a vehicle-based evacuation, and because, (2) in this ABM, movement is constrained by the node-link structure of the network. This representation cannot incorporate the unstructured, open spaces in the study that pedestrians are likely to exploit during an actual evacuation in the study area. In addition, this study recognizes the same bridges (crossing the rivers leading out of the hazard zone) which are bottlenecks to car traffic are also likely to be an impediment to pedestrian travel. Similarly, the bridges are a shared space where pedestrians and vehicles must interact, which is also likely to be important to the outcome of an evacuation. With this in mind, any study making definitive claims about the outcomes of an evacuation in this study area must resolve these tensions. Reconciling in a single model agent travel across open space with agent travel bound to roads is no small challenge. MATSim takes advantage of a directed graph structure to efficiently calculate large numbers of agent activities efficiently, whereas models of movement in unconstrained space typically employ a cellular structure (Lämmel, Grether, & Nagel, 2010). Others have developed evacuation studies where the model space is shared by cars and pedestrians using a cellular representation approach, but the study area was much ! 38 smaller and with fewer agents than were deployed in this study (Mas, Imamura, & Koshimura, 2012). 4.3 Response time: The effect of when agents start to evacuate The results of varying this parameter were statistically significant in the ANOVA test, yet the overall impact if this variable was minimal. The statistically significant result of within-group variance is due to the number of observation in the test and the consistent effect produced for each incremental increase in parameter values. In the pairwise tests, the response time parameter was the only factor that did not register a significant result for between-group variability, meaning there was not enough of an effect produced by the response time parameter to distinguish its effects from that of the other parameters. The one exception was when response time was combined pairwise with population in the iteration 20 scenario. In this scenario, extreme clearance times created by high population and long queues (imposed by the adherence to the absolute shortest path routes) had been relaxed, and the response time parameter had a distinct effect in comparison to population increases in the iteration 1 scenario. In a practical sense these effects were still minimal though. Despite this, that the effect was present provided useful signal of ecological validity: when every agent departs at the same moment the queues form more quickly than when departures times were staggered. Model response time with a Poisson distribution has been used by other studies (Cova & Johnson, 2002), while others have employed sigmoid curves (Mas, Imamura, & Koshimura, 2012), but Lindell and Prater (2007) noted that the sigmoid curve structure (while generally agreed upon as “correct”) becomes more symmetric as the median approaches zero as was first pointed out by Cova and Johnson (2002) in their rational for using the Poisson distribution. The effects of other distributions could be tested to determine whether other methods result in any meaningful differences in a study such as this. 4.4 Network Closures: The effect of having one less way out Exit point closures mean a reduction in the supply of exit routes for a given level of demand, and therefor closures result in higher clearance times across all scenarios. Exit ! 39 closures also restructure the network and shift the availability of road supply in relations to the initial position of all agents in the system. The results showed when the southwest and northwest exits were closed, for iteration 1, clearance times were higher than all other scenarios, but by iteration 20, times were greatly reduced. This shows these two exit points are the greatest attracters in the shortest path scenario and most of the population is closest to these two exits. However, the exit closure scenarios are more difficult to contextualize because, as an artifact of the routing algorithm (if they did not, agents couldn’t find the shortest path), all agents have full knowledge of the network. Thus, agents are aware of these closures when they begin to evacuate. Complete knowledge of network conditions is unrealistic in reality. This software makes a trade off between endowing agents with complete knowledge (which allows for agents to find their own way) and the ability to explicitly direct an agent’s routes, which could permit sending an agent towards an obstacle known to the researcher, for the agent to discover, and thereafter switch to a routing mode where the agent becomes self directed again. This is a software architecture choice, so this could be changed if one were to decide to take on such a challenge. But, doing so is beyond the scope and capability of this study. This is discussed further in section 4.5 with more detail in relation to interpretations of model learning. 4.5 Model learning: How should its effects be interpreted? The discussion above raises questions about the role of agent knowledge and model learning which highlights the inherent methodological challenge of using a routing scheme like the Dijkstra’s shortest path algorithm to move agents in space. If agents are endowed with complete and exact knowledge about the system (that in reality may be impossible have in the best of circumstances) one must proceed with caution and be careful of what inferences are made about an actual population. This is especially true for emergency situations where there is a researcher has an ethical responsibility to interpret results with the utmost concern for stakeholders in a study area. In this research the results are reported as two ends of a continuum of model learning, (1) model iteration 1, where Dijkstra’s shortest path algorithm is employed by each agent, ! 40 and (2) iteration 20, using an adjusted form of the algorithm, where the influence of experience is reflected in the shortest path calculation, a “Nash equilibrium” condition. Lämmel et al. (2010) suggest that these ends can be considered the boundaries of system performance where in one case all agents are absolutely “rational” in their route selection, and a situation where the population is informed by training, and then, agents follow their training regiment to perfection. This study strongly cautions against an interpretation of model results as definitive baseline times for an actual evacuation (for the reasons explained elsewhere in the discussion subsections). However, as noted by Epstein (2008) model results are useful for guiding discussions with stakeholders about how events may take place within a certain system, even if the results need to be qualified as existing within the allowable domain of results for a particular model. 4.6 Directions for future work The population in this study is derived from household counts by census block for the 2010 US Census. For populating a model, it is a very simplistic method, especially for modeling an evacuation scenario, which is a complicated task. Demand modeling can take many forms including estimating if and when agents will decide to leave, assigning where they will go and how they will get there, initial locations, and what exactly is an agent actually meant to represent in terms of the underlying population (Pel, Bliemer, & Hoogendoorn, 2012). This is an area of interest to both transportation and evacuation researchers. As Pel et al. (2012) mention in their review of the multitude of demand models, each method makes assumptions about agent behaviors, perceptions, and what the agent is meant to represent. The choice to use census data at the household level was made because the interest of this research is to uncover the dynamics that constrain an evacuation in this study area and the census derived data provided a way to distribute agents in the system according to one measure of population density in the area. This research acknowledges that this method is too simplistic to make any claims about the number of individuals able to evacuate within a given time frame. Employing a more sophisticated demand model in a similar study is a direction for future work that could provide further insights about an ! 41 evacuation scenario where specific behaviors or agent representations are of interest. It is also important to consider that adding more site-specific detail is not a guarantee that the results will be more useful, nor that the effects of key behaviors will be made more clear. Simplicity has its advantages because the effects generated from each element can be more easily separated from the effects produced by others. The trade off between complexity and simplicity lies within the bigger questions surrounding the purpose of the model and the types of information it is employed to produce. Here, although the demand method is simple, the results shed light on the important effect that space and place have on clearance time variability. 4.7 Conclusion Casting a shadow over all evacuation studies is one big question: What should plans for evacuations and policies for reducing vulnerability look like? In many ways this study confirmed what is already known: traffic congestion is likely to be a problem during an evacuation. This is not a new discovery, in fact, it is widely understood by all stakeholder in the area that a successful evacuation is contingent on what takes place at the river crossings. There are good reasons to strongly consider walking as an alternative to driving during an evacuation in the study area. In many instances, agents closest to the safe zone in terms of walking and ideal drive times were amongst the last to clear the hazard zone by driving (due to bottlenecks at the bridges). Strategies that improve the conditions for walking, like improved wayfinding and an increased number of accessible paths to reach the river crossings would help facilitate a strategy centered on walking. Where impediments to foot travel exist that are common to rural areas, like fences around agricultural fields, the community may be able to agree upon ways to create emergency access points that will allow passage should a lahar occur. A community’s vulnerability is not permanent nor preordained. Although this community (and others like it) are more vulnerable because traffic congestion is likely to hinder an evacuation if a large segment of the population chooses to drive, it doesn’t have to be this way. By addressing the behavioral and environmental conditions known to exacerbate their vulnerability, and promoting those that improve their capacity to adapt, ! 42 communities like this have agency to affect changes that improve their capacity to respond to the threat posed by this region’s large lahars. ! 43 REFERENCES CITED Abdelgawad, H., Abdulhai, B., & Wahba, M. (2010). Multiobjective Optimization for Multimodal Evacuation. Transportation Research Record: Journal of the Transportation Research Board, 2196, 21-33. ! Bernhardt, K. L. (2007, January). Agent-Based Modeling in Transportation. Transportation Research Circular E-C113: Artificial Intelligence in Transportation, pp. 72-80. ! Bonabeau, E. (2002). Agent-based modeling: Methods and techniques for simulating human systems. Colloquium of the National Academy of Sciences, ‘‘Adaptive Agents, Intelligence, and Emergent Human Organization: Capturing Complexity through Agent-Based Modeling.’’. Proceedings of the National Academy of Sciences. ! Chakraborty, J., Tobin, G. A., & Montz, B. E. (2005). Population Evacuation: Assessing Spatial Variability in Geophysical Risk and Social Vulnerability to Natural Hazards. Natural Hazard Review, 6(1), 23. ! Charypar, D., & Nagel, K. (2005). Generating complete all-day activity plans with genetic algorithms. Transportation, 32(4), 369-397. ! City of Orting. (2015). Comprehensive Plan 2015 Draft. ! Cova, T. J., & Johnson, J. P. (2002). Microsimulation of neighborhood evacuations in the urban-wildland interface. Environment and Planning A, 34, 2211- 2229. Crandall, D. R. (1971). Postglacial Lahars From Mount Rainier Volcano, Washington. Geological Survey Professional Paper 677 . ! Crandall, D. R., & Waldron, H. H. (1956, June). A Recent Volcanic Mudflow of Exceptional Dimensions from Mt. Rainier, Washinton. American Journal of Science, 254, 349-362. ! Davis, M., Johnston, D., Becker, J., Leonard, G., Coomer, M., & Gregg, C. (2006). Risk perceptions and preparedness: Mt Rainier 2006 community assessment tabulated results. GNS Science. Institute of Geological and Nuclear Sciences Limited. ! Diefenbach, A. K., Wood, N. J., & Ewert, J. W. (2015). Variations in community exposure to lahar hazards from multiple volcanoes in Washington State. Journal of Applied Volcanology. ! Dijkstra, E. (1959). A note on two problems in connexion with graphs. Numerische Mathematik, 1, 269-271. ! 44 Dow, K., & Cutter, S. L. (1998). Crying wolf: Repeat response to hurricane evacuation orders. Coastal Management, 26(4), 237-252. ! Durst, D., Lämmel, G., & Klüpfel, H. (2012). Large-Scale Multi-modal Evacuation Analysis with an Application to Hamburg. Pedestrian and Evacuation Dynamics. ! Epstein, J. M. (2008). Why Model? Journal of Artificial Societies and Social Simulation, 11(4), 12. ! Finn, C. A., Sisson, T. W., & Deszcz-Pan, M. (2001). Aerogeophysical measurements of collapse-prone hydrothermally altered zones at Mount Rainier volcano. Nature, 409(6820), 600. ! Fothingham, A. S., & Rogerson, P. (1993). GIS and spatial analyitical problems. International Journal of Geographical Infomation Systems, v(1), 3-19. ! Grimm, V., Berger, U., Bastiansen, F., Eliassen, S., Ginot, V., Giske, J., . . . DeAngelis, D. L. (2006). A standard protocol for describing individual-based and agent-based models. Ecological Modeling, 198(1-2), 115-126. ! Grimm, V., Berger, U., DeAngelis, D. L., Polhill, J. G., Giske, J., & Railsback, S. F. (2010). The ODD protocol: A review and first update. Ecological Modeling, 221, 2760-2768. ! Henry, K. D., & Frazier, T. G. (2015). Scenario-Based Modeling of Community Evacuation Vulnerability. Planning, Foresight and Risk Analysis Proceedings of the ISCRAM 2015 Conference - Kristiansand. ! Homer, C., Dewitz, J., Yang, L., Jin, S., Danielson, P., Xian, G., . . . Megown, K. (2015). Completion of the 2011 National Land Cover Database for the conterminous United States-Representing a decade of land cover change information. Photogrammetric Engineering and Remote Sensing, 81(5), 345-354. ! Horni, A., Nagel, K., & Axhausen, K. W. (Eds.). (2016). The Multi-Agent Transport Simulation MATSim. London: Ubiquity. ! Lämmel, G., Grether, D., & Nagel, K. (2010, October 26). The representation and implementation of time-dependent inundation in large-scale microscopic evacuation simulations. Transportation Research Part C: Emerging Technologies. ! Lefebvre, N., Balmer, M., & Axhausen, K. W. (2007). Fast shortest path computation in timedependent traffic networks. 439. Zurich: IVT, ETH Zurich. ! 45 Lindell, M. K., & Prater, C. S. (2007). Critical Behavioral Assumptions in Evacuation Time Estimate Analysis for Private Vehicles: Examples from Hurricane Research and Planning. Journal of Urban Planning and Development, 133(1). ! Liu, S., & Murray-Tuite, P. (2014). Incorporating Household Gathering and Mode Decisions in Large-Scale No-Notice Evacuation Modeling. Computer-Aided Civil and Infrastructure Engineering, 29, 107-122. ! Manson, S., & O'Sullivan, D. (2006). Complexity theory in the study of space and place. Environment and Planning A, 38, 677-692. ! Marzocchi, W., & Woo, G. (2007). Probabilistic eruption forecasting and the call for an evacuation. Geophysical Research Letters. ! Mas, E., Imamura, F., & Koshimura, S. (2012). An Agent Based Model for The Tsunami Evacuation Simulation. A Case Study of The 2011 Great East Japan Tsunami in Arahama Town . 9th International Conference on Urban Earthquake Engineering/ 4th Asia Conference on Earthquake Engineering. Tokyo: Tokyo Institute of Technology. ! MATSim Community. (2015, October 29). MATSim Version 0.7.0. ! Millington, J. D., O'Sullivan, D., & Perry, G. L. (2012). Model histories: Narrative explanation in generative simulation modelling. Geoforum, 43(6), 1025-1034. ! Mothes, P. A., & Vallance, J. W. (2015). Lahars at Cotopaxi and Tungurahua Volcanoes, Ecuador: Highlights from Stratigraphiy and Observational Records and Related Downstream Hazards. In P. Papale (Ed.), Volcanic Hazards, Risks, and Disasters. Elsevere. ! Naghawi, H., & Wolshon, B. (2012). Performance of Traffic Networks during Multimodal Evacuations: Simulation-Based Assessment Read More: http://ascelibrary.org/doi/abs/10.1061/(ASCE)NH.1527-6996.0000065. Natural Hazards Review, 13(2). ! Peeta, S., & Hsu, Y.-T. (2009). Integrating Supply and Demand Aspects of Transportation for Mass Evacuation under Disasters . NEXTRANS Center, Purdue University. U.S. Department of Transportation. ! Pel, A. J., Bliemer, M. C., & Hoogendoorn, S. P. (2012). A review on travel behaviour modelling in dynamic traffic simulation models for evacuations. Transportation, 39, 97–123. ! ! 46 Pierce County Open GeoSpatial Data Portal. (2006, June 1). Retrieved from http://gisdata.piercecowa.opendata.arcgis.com/ ! Pierce County Open Geospatial Data Portal. (2015, April 3). Retrieved from http://gisdata.piercecowa.opendata.arcgis.com/ ! Pierce County Open GeoSpatial Data Portal. (2016, April 8). Retrieved from http://gisdata.piercecowa.opendata.arcgis.com/ Pierson, T. C., Janda, R. J., Thouret, J.-C., & Borrero, C. A. (1990). Perturbation and melting of snow and ice by the 13 November 1985 eruption of Nevado del Ruiz, Colombia, and consequent mobilization, flow and deposition of lahars . Journal of Volcanology and Geothermal Research, 41, 17-66. Pierson, T. C., Wood, N. J., & Driedger, C. L. (2014). Reducing risk from lahar hazards: concepts, case studies, and roles for scientists. Journal of Applied Volcanology 2014, 3:16, 16(3). Pierson, T. C., Wood, N. J., & Driedger, C. L. (2014). Reducing risk from lahar hazards: concepts, case studies, and roles for scientists . Journal of Applied Volcanology, 16(3). Puget Sound LiDAR Consortium. (2004). Pierce County Lowlands 2004 - Bare Earth LiDAR DEM. Retrieved from www.pugetsoundlidar.org Python Software Foundation. (2016). 9.6. random — Generate pseudo-random numbers. Retrieved April 27, 2016, from docs.python.org: https://docs.python.org/3/library/random.html SciPy.org. (2015). numpy.random.poisson. Retrieved April 27, 2016, from http://docs.scipy.org: http://docs.scipy.org/doc/numpy- 1.10.0/reference/generated/numpy.random.poisson.html Scott, K. M., & Vallance, J. W. (1995). Debris Flow, Debris Avalanche, and Flood Hazards at and Downstream From Mount Rainier, Washington . Hydrologic Atlas 729. Scott, K., Macias, J. L., Naranjo, J. A., Rodriguez, S., & McGeehin, J. P. (2001). Catastrophic debris flows transformed from landslides in volcanic terrains: mobility, hazard assessment and mitigation strategies. U.S. Geological Survey Professional Paper 1630. Shiwakoti, N., Liu, Z., Hopkins, T., & Young, W. (2013). An Overview on Multimodal Emergency Evacuation in an Urban Network. Australasian Transport Research Forum 2013 Proceedings. Brisbane: Australasian Transport Research Forum. ! 47 Strader, S. M., Ashley, W., & Walker, J. (2015, June). Changes in volcanic hazard exposure in the Northwest USA from 1940 to 2100. Natural Hazards, 77(2), 1365-1392. Turner, B. K. (2003). A framework for vulnerability analysis in sustainability science. Proceedings of the National Academy of Sciences, 100(14). U.S. Geological Survey. (2009). Retrieved September 1, 2015, from Lessons Learned from the Armero, Colombia Tragedy: http://hvo.wr.usgs.gov/volcanowatch/archive/2009/09_10_29.html U.S. Geological Survey. (2014, 11 10). Retrieved May 5, 2016, from Significant Lahars at Mount Rainier: https://volcanoes.usgs.gov/volcanoes/mount_rainier/geo_hist_lahars.html U.S. Geological Survey. (2014, 11 10). Volcanic Hazards at Mount Rainier. Retrieved from Volcanic Hazards at Mount Rainier U.S. Geological Survey. (2015, September 03). Lahars and Debris Flows at Mount Rainier. Retrieved from https://volcanoes.usgs.gov/volcanoes/mount_rainier/hazard_lahars.html U.S. Geological Survey. (2016, 2 12). Pyroclastic flows move fast and destroy everything in their path. Retrieved from https://volcanoes.usgs.gov/vhp/pyroclastic_flows.html USGS. (2014, November 10). https://volcanoes.usgs.gov. Retrieved from Monitoring Lahars at Mount Rainier: https://volcanoes.usgs.gov/volcanoes/mount_rainier/monitoring_lahars.html Voight, B. (1990). The 1985 Nevado del Ruiz volcano catastrophe: anatomy and retrospection. Journal of Volcanology and Geothermal Research, 42, 151-181. Washington State Department of Transportation. (2015). Population Growth in Relation to the State's Counties. Retrieved September 1, 2015, from http://www.wsdot.wa.gov/planning/wtp/datalibrary/population/ PopGrowthCounty.htm Washington State Office of Financial Management. (2015). 2015 Population Trends. Retrieved September 1, 2015, from http://www.ofm.wa.gov/pop/april1/poptrends.pdf Wood, N. J., & Schmidtlein, M. C. (2012). Anisotropic path modeling to assess pedestrian-evacuation potential from Cascadia-related tsunamis in the US Pacific Northwest. Natural Hazards, 62(2), 275-300. ! 48 Wood, N. J., & Soulard, C. E. (2009). Community Exposure to Lahar Hazards from Mount Rainier, Washington. U.S. Geological Survey Scientific Investigations Report 2009-5211. Wood, N., & Soulard, C. (2009, September 29). Variations in population exposure and sensitivity to lahar hazards from Mount Rainier, Washington. Journal of Volcanology and Geothermal Research, 188, 367–378. ! !