Journal Pre-proof Physics-based simulations of multiple natural hazards for risk-sensitive planning and decision making in expanding urban regions

Rapid urban expansion in many parts of the world is leading to increased exposure to natural hazards, exacerbated by climate change. The use of physics-based models of natural hazards in risk-informed planning and decision-making frameworks may provide an improved understanding of site-specific hazard scenarios, allowing various decision makers to more accurately consider the consequences of their decisions on risk in future development. We present results of physics-based simulations of flood, earthquake, and debris flow scenarios in a virtual urban testbed. The effect of climate change, in terms of increasing rainfall intensity, is also incorporated into some of the considered hazard scenarios. We use our results to highlight the importance of using physics-based models applied to high-resolution urban plans to provide dynamic hazard information at the building level for different development options. Furthermore, our results demonstrate that including building elevations into digital elevation models is crucial for predicting the routing of hazardous flows through future urban landscapes. We show that simulations of multiple, independent hazards can assist with the identification of developing urban regions that are vulnerable to potential multi-hazard events. This information can direct further modelling to provide decision-makers with insights into potential multi-hazard events. Finally, we reflect on how information derived from physics-based hazard models can be effectively used in risk-sensitive planning and decision-making.


Introduction
Urban areas are often hotspots of natural-hazard disaster risks due to their growing populations, ageing and inadequate infrastructure, rising inequality, rapid expansion, and poor spatial planning (e.g., Rosenzweig and Solecki, 2014). For example, an analysis conducted by Gu (2019) assessed exposure and vulnerability to six different natural hazards for 1,860 cities with 300,000 inhabitants or more; 1087 of these cities have a high level of exposure to at least one natural hazard, and 45 cities are highly exposed to three or more natural hazards. Considering multiple hazards is crucial for urban development, as failing to do so can lead to underestimation of the overall impacts and associated risk (e.g., Kappes et al., 2012). Multiple hazards should be considered from the spatial planning perspective when choosing suitable areas for development, as different hazards can affect different areas within an urban region (e.g., Bathrellos et al., 2017). This includes both the expansion and renewal of an urban area, as both can provide an opportunity for reducing disaster risks and introducing new risks. Therefore, risk-informed planning and decision making is considered a prerequisite for safer and more resilient urban futures , Cremen et al., 2022b.
Modelling natural hazards can provide city planners, municipalities, and communities with information on potential hazard scenarios, enabling effective urban design and decision-making that can reduce physical, social, and economic impacts in future urban settings. This information is often presented to stakeholders in the form of hazard maps that display quantitative predictions of susceptibility and/or site-specific intensities of future hazards, for example: landslide susceptibility; peak-ground acceleration or other ground-motion features (e.g., spectral acceleration, significant duration) from earthquakes; maximum flood depth or velocity and flood duration; thickness of volcanic ash deposits; peak wind speeds during storms. Hazard maps are useful for land-use planning and urban design, as they allow urban growth plans to account for threats posed by natural hazards (e.g., Bathrellos et al., 2011;Mesta et al., 2022). However, the use of hazard maps in risk-sensitive urban planning is typically limited to single hazards or multiple independent hazards that do not dynamically interact ). Yet, Global South cities, which account for most global urban growth (e.g., UN Habitat, 2020), will continue to be disproportionally exposed and vulnerable to natural hazards and their interactions. Among other reasons, this is due to a lack of risk-sensitive urban planning (e.g., Johnson et al., 2021;Cremen et al., 2022b). Therefore, integrating models of natural hazards within risk-sensitive urban planning and decision-making frameworks may effectively reduce disaster risk to vulnerable future communities (e.g., Cremen et al., 2022c, This special issue).
Hazard models can broadly be divided into two primary categories: physics-based and non-physical models, including empirical and stochastic methods (e.g., Tilloy et al., 2019). Understanding the advantages and limitations of these model types is crucial, allowing multiple stakeholders to gain the most appropriate information when designing and selecting risk-sensitive urban plans. Non-physical methods apply statistical approaches to observational datasets to estimate broad-scale features of natural hazards. For example, Rickenmann (1999) compiles common empirical relationships between debris flow volume, peak discharge, mean flow velocity and run-out length. These models often have functional forms that are, even if in a simple manner, informed based on physical considerations. Although these and other empirical models are helpful in many applications of hazard assessment, non-physical models are limited when spatiotemporal hazard information on high-resolution urban topography is required. For example, an empirical description of flow inundation calibrated for a particular topography cannot reliably predict inundation where the flow is perturbed by interactions with buildings whose configuration can vary widely in different global cities and future urban layouts. For effective risk-sensitive planning of future urban environments, asset-scale (i.e. at the scale of individual buildings, bridges, roads and other infrastructure) information of hazard intensities and their physical/social impacts is required to reduce disaster risk and increase resilience. In this work, we define 'physics-based' models as methods that simulate the dynamics of natural-hazard events by solving partial differential equations (PDEs) that enforce physical conservation laws (e.g., energy, mass, momentum). Note that physics-based models often incorporate empirical relationships within the system of governing equations. Furthermore, empirical data is often used to prescribe poorly defined aspects of simulation domains. For example, empirical data can be used to define seismic velocity-depth relationships in the absence of a rigorous velocity model (see §3). Physics-based simulations can be applied to high-resolution topography to provide detailed information on hazards in time and space and account for interactions with the built environment on hazard propagation. One advantage of using physics-based models for assisting with future urban planning is that the parameters contained within these models have a physical meaning (e.g., slope angle, basal drag coefficient, rock density, etc.). Hence, changes to these physical parameters in future scenarios can be theoretically accounted for. Conversely, it is difficult to interpret the non-physical parameters that feature in empirical and stochastic models (e.g., Tilloy et al., 2019) and predict how these parameters will evolve with future urban development, climate change, and other physical processes.
Purely stochastic and empirical methods cannot typically provide accurate time-dependent information due to their simplification of the physical processes that control natural hazards. In addition, it is difficult to calibrate empirical and stochastic methods for the largest magnitude events due to their scarcity in historical records. For example, 'giant' earthquakes (moment magnitude, Mw ≥ 9) can have a recurrence interval of hundreds to thousands of years (e.g., Usami et al., 2018), meaning that empirical ground-motion models (GMMs) are often calibrated with more frequent, small-intermediate magnitude earthquakes better represented in empirical databases. This, combined with the simplified empirical representation of physical processes (e.g., fault rupture, seismic wave propagation/attenuation, near-surface soil response), limits the ability of GMMs to provide accurate hazard information that can allow cities to plan for the most hazardous earthquakes (e.g., Freddi et al., 2021). Empirical and stochastic methods often require large quantities of sitespecific data, meaning that they are typically less geographically transferable. This can also be a problem for some physics-based models sensitive to spatial parameter variations. For example, physics-based modelling of earthquake-induced ground motions requires location-specific velocity models of the subsurface. In contrast, empirical GMMs, in some instances, may be more transferable once calibrated for a specific region (e.g., Douglas & Edwards 2016). Furthermore, empirical and stochastic methods are typically less suitable for urban planning, as data in undeveloped areas is often unavailable or poorly resolved in space and time.
Despite the various advantages discussed above, a critical limitation of physics-based models is that a degree of specialist knowledge is required to set up and calibrate physics-based hazard simulations. In addition, physics-based models can be computationally intensive (e.g., Cui et al., 2013;Vasconcellos et al., 2021). As risk-sensitive planning can require iterating through successive urban plans/designs (e.g., Cremen et al., 2022c, this special issue), the computation time of physics-based methods could be a limiting factor relative to other methods. This can be mitigated by using surrogate models or emulators, calibrating against physics-based simulations or real observations to provide computationally efficient, data-driven approximate solutions (e.g., Sarri et al., 2012). While we use this article to outline the importance of physics-based models in advancing effective risk-sensitive planning decisions for future urban development, we also recognise that empirical and stochastic models can also make useful contributions to well-defined components of planning frameworks, such as probabilistic hazard assessment for building code calibration (e.g., Ellingwood 2001). In summary, by understanding the benefits and limitations of the methods outlined above and used below, risk to future communities can be mitigated by using the most appropriate models to simulate natural hazards and their impacts.
the first part of this definition at the urban scale (i.e., multiple hazards that do not occur coincidentally or dynamically interact). Specifically, in this study, we use physics-based models to simulate multiple independent natural hazards in a virtual urban testbed to elucidate the advantages and challenges of using these models for risk-sensitive urban planning. We simulate multiple natural hazard scenarios featuring earthquakes ( §3), fluvial and pluvial flooding ( §4), and debris flows ( §5). All scenarios presented here are simulated in 'Tomorrowville', a virtual urban testbed introduced below in §2. In §6, we synthesise our results and reflect on the identification of potential multi-hazard events. Finally, we reflect on the challenges researchers/modellers, stakeholders, and city planners face in using physics-based models for risk-sensitive urban planning and decision making. We outline how the Tomorrow's Cities Decision Support Environment (TCDSE) design and implementation can limit these challenges.

Overview of Tomorrowville
Tomorrowville is a virtual urban testbed which has been designed to represent typical physical and socio-economic aspects of evolving cities in the Global South ( Figure 1). An active river channel is located at the base of the valley floor. During intense rainfall, ephemeral channels develop in incised valleys. Settlements are located on the eastern plateau and clustered in the central valley region. For the purpose of this analysis, we assume that the region is affected by three types of hazards: flooding, debris flows, and earthquakes.
The topographic domain of Tomorrowville is derived from a 2 m resolution digital surface model (DSM) of a ~6 km 2 section of the Kathmandu Valley (KTMV), derived from tri-stereo Pleiades satellite imagery. This area of the KTMV was selected to represent Tomorrowville due to several key considerations: principally, community engagement identified three natural hazards (earthquakes, flooding, and debris flows) that impact this area; this area is covered by a high-resolution digital elevation model (DEM) that can be used for high-resolution urban flooding and debris flow simulations; we have access to lithological/sedimentological data which can be used to determine model parameters.   In the scenarios presented here, several assumptions were made in developing the physics-based models. For example, in a real landscape, model parameters would have to be calibrated and validated against observed historical events to ensure appropriate accuracy. In the absence of reliable historical data commonly used to calibrate the models, it would be essential to explicitly define the bands of uncertainty associated with uncalibrated models. Because of the nature of this research, where the models have been developed for a virtual urban testbed for demonstration purposes of the TCDSE concept, the calibration of model input parameters and source conditions and validation of model output against observations has not been carried out for Tomorrowville. Instead, this study presents comparative scenarios for each hazard, with scenario intensities selected to represent typical hazards that affect Global South cities (see below). These scenarios are subsequently used to demonstrate the importance of modelling multiple hazards and incorporating climate change projections when designing risk-sensitive future urban plans. For the TCDSE, various hazard scenarios would be simulated for each hazard to produce a database of dynamic intensity maps. These can be, in turn, used within the TCDSE's 'Physical Infrastructure Impact' module to assess impacts from natural hazards on the built environment, for instance, through appropriate fragility (i.e., the probability of various damage levels as a function of a hazard intensity measure) and consequence models (e.g., Gentile et al., 2022; Cremen et al., 2022c, this special issue).

Earthquake hazard
Earthquakes produce damaging shaking, which globally has caused billions of dollars of economic losses and almost one million fatalities so far in this century (e.g., EM-DAT 2021). In earthquake-prone regions, understanding the potential for strong ground shaking is a critical constraint for the design of new structures and the performance assessment of existing structures and infrastructure components within the built environment. This has led to the development of a range of methods to assess seismic hazard (e.g., Cornell 1968; Kramer 1996; Baker, Bradley and Stafford 2021). Seismic hazard analysis relies on two main components: (i) modelling potential rupture scenarios in the vicinity of the region of interest; and (ii) estimating the ground-motion properties (or local intensities) caused by a given scenario rupture. Compilations of ground-motion records from past earthquakes of different sizes across an area have led to a series of empirical GMMs -based on regression analysis -estimating the probability distributions of ground-motion intensity measures at a site, given an earthquake of a certain magnitude occurring at a nearby location. Specifically, these models relate the magnitude and style of a seismic event, the distance from the source, and other site-specific parameters (e.g., local soil properties) to the amplitude and other features (e.g., significant duration) of the expected shaking and its variability at a given target site (e.g., Douglas and Aochi 2008; Douglas and Edwards 2016; Douglas 2017). While providing a systematic and effective method for estimating ground-motion parameters, such ground-motion modelling suffers from several critical shortcomings. First, the paucity of high-magnitude/high-intensity shaking recordings at small distances requires extrapolation of the models into the most important parts of the domain. The scarcity of instrumental recordings, particularly in the Global South where urban expansion is concentrated, often requires the import of GMMs into areas for which their efficacy is questionable and where validation is not always possible. Finally, nonlinear time-history analysis of engineered systems arguably represents the most advanced procedure for structural seismic performance assessment and fragility analysis (e.g., Gentile and Galasso, 2021; Silva et al., 2019). Such analysis procedures require reliable ground-motion time series to investigate structural response to dynamic loading. Generally, the input ground-motion time series for nonlinear structural analysis are selected (and eventually scaled) from a database of existing records to represent target seismic characteristics at a given location (e.g., Iervolino et al., 2010). As discussed above, the inherent scarcity or total absence of suitable recorded ground motions for some specific scenarios (e.g., large-magnitude strike-slip events recorded at close source-to-site distances) makes the use of alternative options unavoidable.
In particular, physics-based ground-motion simulation has emerged as a complementary and possibly alternative approach to empirical ground-motion modelling, addressing the abovementioned challenges ( . These physics-based simulated ground motions can capture complex source features (such as spatially variable slip distributions, rise-time, and rupture velocities); path effects (geometric spreading and crustal damping); and site effects (wave propagation through basins and shallow site response), providing a valuable supplement to recorded ground motions. Hence, in the specific context of high-resolution hazard modelling for urban planning, physics-based ground-motion simulations can provide estimates of full ground-motion time series, enabling an improved assessment of the complete range of ground-motion features. Moreover, these waveforms can be coupled with advanced nonlinear structural models to perform nonlinear time-history analysis of buildings and infrastructure components and derive sophisticated, numerical models of structural fragility (e.g., Gentile et al., 2022, this special issue). In this way, physics-based ground-motion simulations provide risk modellers and decision-makers with a scientifically rigorous basis for estimating the likely impacts of strong shaking at the scale of individual assets.
Below, we present simulations from physics-based and empirical GMMs to compare the advantages and limitations of these methods and to provide insight into how physics-based simulation can assist in urban planning and decision-making. We examine the potential of a physics-based set of earthquake simulations to expose some invariants of the shake distribution for a series of idealised events around Tomorrowville. We demonstrate that the relative amplitude of local shaking is strongly constrained by the near-surface velocity structure and identify robust patterns of relative amplitudes of ground shaking that permit detailed predictions of spatial distributions of relative shaking amplitude regardless of the fault mechanism, locations, directivity or even magnitude of the event.
3.1 Model setup

Geological setting
We first construct a simple but reasonable velocity structure for the crust below our testbed area. 3D crustal models (often referred to as 'velocity models') provide the 3D domain over which the wave equation is solved. The crustal model provides the 3-D variation of geophysical and geotechnical parameters required for wave propagation calculations. The principal parameters to describe the model are the P-and S-wave velocities, density, and anelastic attenuation -additional parameters are needed if nonlinearities are considered.
Tomorrowville, similarly to many areas slated for urban development, is located around a river channel, under which a deeper sedimentary layer is assumed. Soft sedimentary layers influence earthquake-induced ground motions (e.g., Roten Figure  2 presents the extent of the river channel within the wider basin, the location of Tomorrowville, and the shear-wave speed (V s ) for the entire simulation domain.
The depth-dependent velocity structure of the simulation domain, including density ( ), the shear wave speed, V s , primary-wave speed, V p , and anelastic attenuation factors (Q s and Q p ), is generated based on the values assumed for V s at the surface for the river channel, basin interior, and basin exterior (see Table 1). We use Brocher (2005; equation 9) to relate V p to V s for the basin exterior region. Note that in the absence of site data, the use of an empirical velocity model is essential to enable the estimation of ground motion using physics-based simulations. We assume that anelastic attenuation is given by = /20 and = 2 . In addition, we introduce some stochastic variability in the underlying crust by perturbing the velocity at each point, producing a fractal spatial correlation in the velocity structure. Currently, this perturbation has an arbitrary amplitude proportional to the velocity at any point but provides realistic spatial correlations consistent with observations in deep boreholes (e.g., Bean and McCloskey, 1993). These choices are unlikely to significantly alter the main conclusions of this paper.  We now explore the shaking modelled across Tomorrowville due to different earthquake events. Here, we compute the shaking at every node of the surface of the computational domain so that we can compute the peak spectral accelerations at each point for each event. We simulate 14 Mw = 6 events, illustrated in Figure 3a; and 10 Mw = 5 events, illustrated in figure 3b. As shown in Figure 4, the slip distribution remains the same for each magnitude. Still, the strike, dip, rake and hypocentral depth vary, producing a wide range of expected source directivity for any location. The aim of considering these two sets was to scrutinise the effect of the underlying geological structure on the spatial groundmotion intensity measures for a range of modelled events. Two simulation cases (scenario EQ1 from the Mw 5.0 set and scenario EQ2 from the Mw 6.0 set) are used to assess the utility of physics-based ground-motion simulations for risk-informed urban planning (see §3.3). Kinematic rupture models for the considered scenarios in Figure 4 are generated based on the model developed by Liu et al. (2006) and Schmedes et al. (2013), in which the correlation between the slip, rise time, peak time and rupture velocity in the neighbouring sub-faults are considered based on a large ensemble of dynamic simulations. Figure 4 presents the moment release across the rupture for the Mw 6.0 and Mw 5.0 scenario events.

Results
The results are presented in Figure 5 for 5%-damped pseudo-spectral accelerations (pSAs) at the vibration period of 1.0 and 2.5 s. These periods can capture the medium-and long-period groundmotion amplitudes relevant for the physical-damage assessment of any mid-and high-rise buildings in Tomorrowville. The geometrical mean of the two horizontal ground-motion intensity measures is shown only for the Tomorrowville region. Figure 5 shows that the soft sedimentary layer near the river channel and the basin interior amplify the ground shaking significantly for both Mw 6.0 (i.e., Figure   These results are presented in the Tomorrowville region as well as a 100x100 km extent to reveal the general trend in estimates from empirical models. It is noted that the empirical results are produced for a V s 30 = 300 m/s which is the V s 30 value in the basin and river channel region (with no variation due to the specific depth sampling considered in the adopted velocity model). The intensity measure estimates follow a general race-track pattern. These results, which are representative outputs of empirical models, do not reproduce the small-scale, site-specific spatial variation that can be simulated using physics-based models (Figures 5 and 7).   All flood scenarios discussed here are assumed to occur during a monsoon season. It is assumed that the soil is fully saturated by pre-monsoon rainfall so that, during the monsoon season, all rainfall becomes surface runoff. Although some runoff may be diverted to stormwater drainage systems, it is reasonable to assume that during the intense rainfall events presented in this study, typical urban drainage systems would reach capacity quickly, and the majority of rainfall would contribute to surface runoff. The TOPMODEL m parameter (Beven and Kirkby, 1979), required to calculate rainfall runoff in the Caesar-Lisflood model, is chosen to be small (m = 0.003) to represent this reactive rainfall-runoff system. In the absence of detailed information about land cover, sediment grain size, and bedforms, a uniform Manning's coefficient (which characterises any surface roughness that impacts water flow) of 0.04 (m 1/3 s -1 ) is used in all flood scenarios presented here (Chow, 1959). In a real urban landscape, a spatially distributed Manning's coefficient could be calibrated to improve the accuracy of flowrouting through a domain with multiple land cover and channel types. . Therefore, using these two discharge hydrographs and rainfall time-series data allows us to explore the potential impact of climate change on future flooding hazard in Tomorrowville.
The initial conditions for all three scenarios are identical, with an initial water depth in the river reach only and no water on the floodplains. The discharge hydrographs shown in Figure 8 are applied at the upstream boundary for all scenarios. For scenarios F2 and F3, as well as a discharge hydrograph, an hourly rainfall time series (mm/hr) of an estimated 1-in-25 year and 1-in-100 year rainfall event, respectively, is applied to every cell in the model domain ( Figure 8). Because of the size of the catchment (~6 km 2 ), we use a uniform rainfall distribution across the catchment. An open boundary condition is used at the outlet of the river (Figure 9).

Results
Maps of maximum flood inundation for each of the three scenarios are shown in Figure 9. For all scenarios, overbank flooding occurs along certain parts of the river, inundating a small number of buildings with up to 1-2 m of water. In Figure 9a, where intense rainfall is neglected (F1), it can be seen that the flood inundation is confined to the region close to the river bank. Comparing the results of Figure 9a-c, the flood inundation across the catchment increases when localised rainfall is considered. This highlights the need to use a combined fluvial/pluvial model to generate appropriate flood maps in urban regions prone to flash flooding. In addition, because of the terraced nature of the landscape, in Figure 9b-c, the ephemeral rivers channel rainfall runoff to certain low-lying parts of the To test the effect of building footprints on the flow routing, we ran one additional scenario for the 1in-25 year combined flood event using the 2 m DEM, which did not have TV0 building heights embedded. The results are shown in Figure 10. For the overall catchment, the inundation pattern is similar. However, when we zoom in on the built-up areas (e.g., Figure 10), we can see that the flood pattern is different, and the flood depths are greatly reduced when buildings are not considered (less than 0.4 m compared to depths more than 1 m to 1.5 m when TV0 buildings are included). The presence of the building footprints in the DEM acts as an obstacle to the flow, causing it to be routed down certain streets and lanes between buildings. Because the flow cannot 'spread out' as it does in Figure 10b,

Debris flow hazard
Rainfall-driven debris flows and flash floods are energetic flows that often feature strong erosion and deposition (collectively referred to as morphodynamics). Their flow dynamics are inherently coupled to morphodynamics by: increasing or decreasing the mass of the flow through erosion or deposition of sediment, respectively; modifying the topography over which present and future flows are routed; and modifying the magnitude and form of basal drag. Debris flows damage buildings via inundation and by exerting dynamic forces on structures. When hydrodynamics forces are neglected, these forces are approximated by the hydrostatic pressure in the flow, allowing for the flow depth to be used as a proxy to estimate building damage resulting from debris flows, lahars, flash floods and other shallow overland flows (see , in this special issue, for a detailed discussion on appropriate fragility models).

Model setup
We simulate rainfall-driven debris flows using the LaharFlow dynamic hazard model Langham et al., 2021). Previously LaharFlow has been used to simulate two-phase, overland flows that feature strong morphodynamics, such as lahars, flash floods, and debris flows in arid-urban environments (Tierz et al., 2017;Hogg, 2018;. LaharFlow, as with many other models that simulate surface flows, utilises a shallow-layer formulation which assumes that the pressure distribution is hydrostatic and flow properties are vertically averaged. The model comprises a coupled system of partial differential equations that enforce conservation of mass and momentum in both the fluid and solid phases (Langham et al., 2021). These equations are solved numerically over a 2 m resolution DEM of the target area in which building elevations are embedded. Note that in conserving mass in the fluid and solid phases separately, the model explicitly solves for the solids concentration of the flow, which is allowed to vary with morphodynamics. The model adopts wellestablished closures for morphodynamics: the erosive flux is proportional to the excess Shields stress (Meyer-Peter and Muller, 1948) and an erosion rate coefficient that is a function of the flow concentration. This smooth, monotonic function (hyperbolic tangent) accounts for the phenomenology that, for a given momentum flux, granular flows are more erosive than dilute flows. The depositional flux is calculated using a hindered settling law that accounts for the suppression of deposition as sediment concentration increases (Soulsby, 1997;Richardson and Zaki, 1954).
LaharFlow utilises a novel drag formulation to capture different basal stress modes present in flows across the entire range of sediment concentrations, from dilute, fluid-dominated flows to concentrated granular flows. Each of these limits utilises a widely adopted drag law: dilute flows use a Chézy drag formulation (e.g., Shalaby, 2018), whereas concentrated flows adopt a Coulomb law with a velocity-dependent friction coefficient (Forterre and Pouliquen, 2008). A smooth, monotonic weighting function is used to determine the magnitude of basal drag for intermediate values of sediment concentration. Note that in the absence of morphodynamics at the dilute limit, the system of equations reduces to the familiar shallow water equations.
The source condition must be carefully chosen to ensure that the resulting debris flow is realistic. The domain is split into broadly three sub-domains that are approximately north-south oriented: a central river and low-relief zone separate two high-relief plateaux. Debris flows and other rainfall-driven overland flows require topographic focusing of surface runoff during high-intensity rainfall. The eastern plateau contains several valleys that are eroded by ephemeral channels that transport water and sediment to lower relief areas. We select a single valley (catchment area ~ 0.2 km 2 ) as a location for the source of scenarios because the selected valley, and associated micro-catchment, is sufficiently large to focus a hazardous volume of surface runoff during intense rainfall, potentially inundating an area that could be deemed suitable for urban development. In contrast, smaller valleys/catchments are less likely to produce hazardous flows.
Debris flows are initiated by a volumetric flux distributed over a circular region with a radius of 50 m (location indicated in Figure 12). This source area is chosen to be sufficiently large to ensure that flow in the source region respects the shallow-layer assumptions by preventing flow velocities and erosion in the source region from being unrealistically strong. We also ensure that flow is only introduced in the intended valley by selecting a source area smaller than the valley's width. The volumetric flux is calculated by multiplying the catchment area (0.2 km 2 ) by a time series hydrograph of rainfall intensity (m.s -1 ). In our assumptions, we neglect infiltration and other water sinks during strong surface runoff promoted by the steep slopes in the source region.
We use peak rainfall intensities of 10, 50, and 150 mm hr -1 for low-, medium-, and high-intensity debris flow scenarios, respectively. Note that the rainfall used to initiate debris flow scenarios is more intense (i.e., higher hourly rate) and of shorter duration than used for pluvial flooding scenarios. We find that, for the topography present in the target valley, short-duration, extreme intensity rainfall is required to erode and entrain non-negligible amounts of sediment necessary for these flows to be classified as debris flows. We base our source conditions on empirical intensity-duration thresholds for rainfall events that triggered recorded debris flows and landslides (Guzzetti et al., 2008). Despite the different hourly rainfall rates used for pluvial flooding and debris flow scenarios, the cumulative rainfall used to initiate these flow hazards is similar. As climate change will likely increase the frequency and intensity of extreme rainfall events in the case study area (e.g., Sheikh et al., 2015), considering a broad range of source rainfall intensities also allows us to account for a range of potential future climatic conditions in Tomorrowville. The source fluxes for these scenarios are summarised in Figure 11. In each scenario, we impose the following source conditions: 1. The source flux linearly increases from 0 (m 3 /s) to the maximum value over the first 500 seconds. The maximum flux is maintained for 50 minutes and then linearly removed over a period of 100 seconds. Linearly ramping the source up and down assists with model stability and limits the formation of artificial fronts in the flow, which sharp changes in the source flux can introduce. Note that in the absence of historical observations or flow data, the simple structure of an idealised hydrograph allows for a simpler interpretation of the resulting flow relative to more complex hydrographs. 2. The debris flow is initiated by setting the sediment concentration in the source flux to 10 % solids (by mass). The addition of sediment allows for equilibrium between erosion and deposition to be established at earlier times. The addition of sediment to the source also accounts for sediment that is likely to be entrained whilst being routed from catchment to source. The sediment concentration is linearly ramped up to and down from its maximum value of 10 %, as described for the volumetric source flux. 3. We continue running the simulation for an hour after the source is removed. This ensures that hazard maps include locations where the topography has focused the flow after removal of the source (typically furthest downstream).
The source condition is varied in all scenarios, and all other parameters are fixed (summarised in Table  2). Given that Tomorrowville is a virtual urban testbed, model parameters are selected to be consistent with values in the literature for similar magnitude debris flow events. The erosion parameters listed in Table 2 were tuned to yield erosion depths (tens of cm) consistent with direct measurements of erosion from debris flows (Berger et al., 2011). Although stronger erosion has been observed from debris flows (> 1 m), we limit erosion to a maximum depth of 1 m to ensure we are not capturing extreme behaviour. As discussed above, LaharFlow utilises a novel drag law that is a smooth function of solids concentration. Table 2 summarises the end member drag limits: a dimensionless Chézy drag coefficient of 0.04 lies within the range of calibrated field events studies (e.g., Rickenmann et al., 2006). The minimum and maximum Pouliquen slope values correspond to the tangents of the angles at which a granular layer will begin and cease flowing, respectively (e.g., Pouliquen and Forterre, 2002). The Voellmy switch value can be interpreted as the sediment concentration (20 %) at which granular effects begin dominating basal shear stress. Finally, the switch rate is tuned to control the sharpness of this transition. J o u r n a l P r e -p r o o f

Results
We represent the hazard from debris flow scenarios using maps of maximum flow depth (see Figure  12). In all scenarios, debris flows are routed through a single channel in the valley downstream of the source before spreading on the valley floor; partially inundate the settlement located on the valley floor; transport water and sediment into the river after following a channel to the north of the settlement. When simulating debris flow scenarios, we do not actively model fluvial transport of water/sediment that enters the river channel. Therefore, when debris flows enter the river channel, conservation of mass and momentum leads to the artificial development of expanding lakes (see Figure 12). In reality, fluvial transport would likely prevent lake formation by advecting material downstream. Given that, in all simulations, these lakes never inundate any settlements or individual buildings, this assumption is unlikely to significantly impact debris flow routing through the urban settlement. Furthermore, we expect the effects from backward (i.e., upstream) propagation of the flow from these lakes to be limited, as these lakes mostly develop after the occurrence of maximum flow depth in the urban settlement. The strength of the source condition determines the severity of inundation in the settlement. Figure 12 demonstrates that as the source intensity increases, the extent (i.e., number of assets affected) and depth of inundation increase. For higher intensity scenarios, an additional channel, located towards the southern extent of the settlement, becomes active, propagating the flow through the domain (Fig 12b&c). In moderate to low scenarios, some ephemeral channels never become active. Given the large volume of material that can be deposited into the river for high-intensity scenarios, there is scope for considering multi-hazard interactions between fluvial flooding and debris flows (see below).

Discussion and conclusions
Our results demonstrate that physics-based models can be used to effectively simulate multiple independent hazards that impact an urban area. We now discuss the utility of these methods when developing risk-sensitive urban plans that consider multiple natural hazards. Physics-based models can produce high-resolution spatially distributed maps of hazard intensity for a range of scenarios. These can help decision makers identify regions of the landscape exposed to one or multiple hazards, from district level to individual buildings. When considering future urban development, the dynamic nature of the landscape and the hazards and the dynamic interactions between hazards and landscapes must be accounted for. For example, as the urban landscape evolves, the construction of new assets or the demolition of others can change the routing of flows (e.g., flood, debris flows, fire, etc.) through the landscape, thus changing the spatial intensity of flow hazards. This can be achieved using physics-based simulations only. We demonstrated this in the flood hazard modelling scenarios. When building footprints were incorporated into the DEM directly, the flood inundation increased in certain regions of Tomorrowville ( Figure 10). These features associated with the high-spatialresolution mapping of building layouts are unlikely to be captured in empirical models that are not calibrated at this resolution. Physics-based models can also account for changes to flow behaviour and source conditions that arise from changes to surface runoff characteristics (e.g., from the construction of paved roads and new buildings). In addition, when considering future urban development, it is essential to account for the effects of climate change. Physics-based models can account for perturbations to parameters and source conditions likely to be caused by climate change. We demonstrated this in §4, where we incorporated rainfall projections from Shrestha et al. (2022) into our flooding scenarios to demonstrate how climate change may alter flood hazard in Tomorrowville in the near future. This highlights the need to base future risk-sensitive urban plans on predicted rainfall events or projected climate scenarios and not only on return intervals produced from historical events. Physics-based models can be used to develop hazard maps for combinations of different physical urban layouts and hazard intensities when exploring potential future urban spaces.
The spatial correlations of hazards occur due to their dynamic behaviour and the physical characteristics of the domain. Most notably, settlements in the vicinity of the river channel are most vulnerable to future hazards due to: (a) overbanking of the river (fluvial flooding); (b) damage from debris flows and pluvial floods that are routed onto the valley floor; (c) maximum shaking from earthquakes occurs due to the soft sediments that surround the river channel and due to the basin geometry. Our results demonstrate that physics-based models can identify regions within Tomorrowville that are vulnerable to multiple independent hazards. Furthermore, using empirical GMMs as an example (Figure 6), we demonstrate that empirical methods cannot resolve local-scale variations that contribute to the spatial correlation of hazards. It is important to identify multi-hazard interrelationships, as the cumulative impact of multiple interacting hazards can be greater than the sum of their individual impacts (e.g., Kappes  . We have shown that independent physics-based hazard simulations, in combination with the literature and stakeholder knowledge, can be used to assist with the identification and location of multi-hazard prone regions. For Tomorrowville, we identify -at least qualitativelypotential multi-hazard interactions from simulations of individual hazards. In addition, for applying the TCDSE in an actual location, we also envisage a process of compiling and analysing a database of previous hazard and multi-hazard events from community knowledge, municipal government, risk management and media records. The resulting simulations should utilise DEMs that account for future urban designs (i.e., that represent future distributions of assets and populations). Several challenges exist in developing physics-based models to simulate natural hazards in an urban environment, including the identification of major hazards and their potential interactions, the choice of hazard intensity scenarios, the choice of model parameters, and the uncertainty associated with all of the above (e.g., Tilloy et al., 2019). However, using state-of-the-art physics-based models of multiple hazards can help to identify areas in a rapidly expanding urban landscape which are highly sensitive to multiple hazards. Understanding how different hazards interact with a landscape can help modellers, urban planners, policymakers and communities to identify areas more suitable for residential, industrial, leisure, and other land uses and for planning/designing future infrastructure systems.
Validation of physics-based models is essential to understand their accuracy and capability. In the case of a physics-based model for earthquake-induced ground motions, the appropriateness of input models describing the source rupture, 3D crustal structure, and shallow site conditions are complex and regionally varying. Hence, the accuracy of ground-motion simulations is region-and even sitespecific. Furthermore, as ground-motion time series are complex transient signals, the power of simulations to adequately capture the relevant features of these signals varies depending on the specific focus of the analysis. When modelling natural hazards for risk-informed planning purposes, numerous sources of uncertainty must be understood and accounted for. Such modelling uncertainty is due to inaccuracies and idealisations made in the physics-based model formulations and the choices of source conditions and input parameters. While a detailed discussion of uncertainties in physics-based hazard simulations is beyond the scope of this work (and is also hazard-and method-dependent), we provide some general discussion and a few illustrative examples here. Specifically, three principal sources of modelling uncertainty in physics-based hazard models must be communicated to and considered by decision makers, including: errors in the physical domain; errors in the parameters contained within the governing equations; errors in the source/input conditions used to propagate hazards throughout the domain. Some of these issues can be addressed by thoroughly calibrating and validating the model using past events. However, it is impossible to calibrate and validate future projected hazards that depend on the urban layout (e.g. floods, mass movements, fires) that occur in regions that have yet to be developed. This is especially important for modelling the routing of hazardous flows since climate change, and alterations to the physical environment caused by urbanisation will significantly affect modelling outputs, as discussed above. Physics-based models can simulate these physical changes by carefully selecting model parameters and source conditions based on future land use and climate projections. Moreover, different hazard models (and different modelling approaches/methods) will have different levels of uncertainty. Therefore, decision makers must be aware of, and adequately account for, those uncertainties when attempting to reduce risk in future urban environments. When using simulation-based risk-reduction frameworks, uncertaintiesincluding those related to hazard intensities-should be propagated through to the relevant impact metrics. For instance, the simulation-based framework for earthquake risk-informed and peoplecentred decision making in future urban planning proposed by Cremen et al. (2022b) explicitly accounts for uncertainties in the future projections of underlying variables (e.g., asset location and structural or non-structural features, building fragility, age and income profile of inhabitants). Monte Carlo sampling is extensively used to capture uncertainties in end-to-end calculations, from hazards to impacts and risk.
To illustrate some of these issues further, let us now consider uncertainties related to the domain of flow hazard simulations. A large source of error when simulating present-day flow hazards in urban environments is caused by artefacts in the DEM, such as: trees and vehicles, which provide obstacles to flow, and bridges over river channels, which can act as dams. These errors can be accounted for in present-day simulations by carefully processing the DEM to identify and remove these artefacts. However, the eventual layout of future urban environments is typically unknown precisely. Given that we have demonstrated the importance of incorporating future building layouts into the domains of flow simulations (see §4), decision-makers must recognise that perturbations to future urban plans can effectively alter future flooding hazard and therefore acts as a source of uncertainty in the riskinformed modelling process.
Some uncertainties associated with the limitations of individual physics-based models can be effectively communicated to stakeholders and decision-makers by leveraging a combination (ensemble) of different models for a given risk assessment. However, this approach requires an accepted protocol for ranking and aggregating models. The uncertainty associated with input parameters and source conditions can be accounted for in the planning process by simulating a range of scenarios co-produced with local stakeholders and decision-makers. For example, in the absence of historical flood maps from satellite imagery or maps developed as part of a post-disaster survey, communities can work with modellers to identify areas prone to flooding (e.g., Mulligan et al., 2019) as communities experiencing frequent natural hazards are known to have a wealth of local knowledge which can inform the modelling process (Sakic Trogrlić et al., 2021). Indeed, a key aspect of the TCDSE is to reduce risk in future cities by engaging with a broad spectrum of stakeholders throughout the planning stage through a process of co-production .
In summary, we use a case study of a virtual urban testbed to highlight the importance of using physics-based models to simulate multiple natural hazards when designing risk-sensitive urban plans. We present simulations of independent earthquake, flood, and debris flow scenarios which we use to reflect on the process of integrating hazard models into risk-sensitive planning frameworks. We demonstrate that physics-based methods can be used to assist with the identification of urban areas that are vulnerable to multi-hazard events. We also show that including building elevations in highresolution DEMs of future urban topography is necessary for planning that accounts for the routing of hazardous flows. We conclude by recommending that future risk-sensitive planning frameworks include physics-based models that can dynamically simulate coupled multi-hazard events that are identified through advanced modelling and co-production with local stakeholders.