Variable-resolution building exposure modelling for earthquake and tsunami scenario-based risk assessment. An application case in Lima, Peru

We propose the use of variable resolution boundaries based on Central Voronoi Tessellations (CVT) to spatially aggregate building exposure models for risk assessment to various natural hazards. Such a framework is especially beneficial when 15 the spatial distribution of the considered hazards present intensity measures with contrasting footprints and spatial correlations such as in coastal environments. This proposal avoids the incorrect assumption that a single intensity value from hazards with low spatial correlation (e.g. tsunami) are considered as representative within large sized geocells for physical vulnerability assessment, without, at the same time, increasing the complexity of the overall model. We present decoupled earthquake and tsunami scenariobased risk estimates for the residential building stock of Lima (Peru). We observe that earthquake loss models for far-field 20 subduction sources are practically insensitive to the exposure resolution. Conversely, tsunami loss models and associated uncertainties depend on the spatial correlations of the hazard intensities as well as on the resolution of the exposure models. We observe that for the portfolio located in the coastal area exposed to both perils in Lima, the ground-shaking dominates the losses for lower magnitudes whilst the tsunami does for the larger ones. For the latter, two sets of existing empirical flow-depth fragility models are used, finding large differences in the losses. This study arise awareness about the uncertainties in the selection of 25 fragility models and aggregations entities for exposure modelling and loss mapping.


Introduction
The spatial distribution of damage and/or losses expected to be incurred by an extensive building portfolio from a natural hazardous event can be quantified and mapped once a physical vulnerability analysis is performed. For such a purpose, a set of fragility functions per building class is conventionally used. Fragility functions describe the probability of exceeding a certain damage limit 30 state for a given intensity measure (IM) associated with natural hazards, such as spectral acceleration at the yield period (e.g. Fäh et al., 2001) for earthquakes, or tsunami inundation height for tsunamis (Koshimura et al., 2009). These vulnerability calculations are performed at the centroid of every aggregation area of the building exposure model with certain bias (Bazzurro and Luco, 2005), or over weighted points (e.g. Weatherill et al., 2015;Kappos et al., 2008). These aggregation entities can be very diverse, https://doi.org/10.5194/nhess-2021-70 Preprint. Discussion started: 16 March 2021 c Author(s) 2021. CC BY 4.0 License. ranging from administrative units such as district/ communes (e.g. Dunand and Gueguen, 2012); urban blocks (e.g. Papathoma and 35 Dominey-Howes, 2003;Kappos et al., 2008;Figueiredo et al., 2018;Kohrangi et al., 2021); regular grids (e.g. Erdik and Fahjan, 2008;Figueiredo and Martina, 2016) or variable-resolution CVT (Central Voronoi Tessellation) geocells .
Throughout the physical vulnerability assessment, it is implicitly assumed that the intensity observed or estimated at that location (i.e. centroid or weighted points) is representative for the entire aggregation area. Depending on the considered hazard footprint and IM attenuation, this assumption might not be longer valid if the aggregation area is too coarse compared to the correlation 40 length of a highly varying IM. Besides the aggregation of the building exposure itself, the importance of these geographical entities in natural hazard risk assessment is that they are ultimately used to calculate and map the expected damage and loss metrics (e.g. building replacement/repair costs, human casualties). The diverse types of visualisation and interpretations of this kind of geospatial data define the so-called thematic uncertainties (Smith Mason et al., 2017) that can largely impact further decision making processes (Viard et al., 2011). It is then important to find a compromise between the intrinsic resolution of the hazard IM, 45 on one hand, and the cartographic representation of the exposure models and risk metrics on the other hand.
When a geographically distributed hazard intensity measure IM presents no significant spatial variability within distances of the order of tens of kilometres, they are said to be highly spatially correlated (e.g. Gill and Malamud, 2014;Merz et al., 2020). This is the case of hazards with relatively low attenuation and wide-spread footprints such as ash-falls and earthquakes (de Ruiter et al., 2021). For the latter case, when seismic site conditions (e.g. soil amplification) and path effects (e.g. seismic directivity) are 50 insignificant, seismic ground motion correlation lengths between 10 km and 25 km are typical (e.g. Esposito and Iervolino, 2012;Schiappapietra and Douglas, 2020). On the other hand, hazards are called low-spatial correlated if their IM are highly prone to be modified by specific features of the propagation medium. For instance, the modelling of inland IMs from a tsunami (i.e. inundation depth, flow velocity, momentum flux) are highly dependent to be modified by the resolution of the bathymetry and digital elevation models (e.g. Tang et al., 2009); coastal topography (e.g. Goda et al., 2015); coastal morphology (e.g. Song and Goda, 2019), and 55 even built-up areas that heavily interact and modify the inundation footprint and flow velocities (e.g. Kaiser et al., 2011;Lynett, 2016). Moreover, the maximum tsunami IMs also depend on the properties of the triggering earthquake such as its magnitude (e.g. Goda et al., 2014), slip distribution (Miyashita et al., 2020), and directivity of energy radiation (e.g. Kajiura, 1972). Thus, the spatial correlation of inland IMs from tsunamis is very low and remarkably non-linear compared to a much more uniform and highly spatially correlated seismic ground motion. Efforts to visualize uncertainties in the tsunami hazard and risk mapping 60 addressing some of the aforementioned modifiers have been reported in few previous studies (e.g. Goda and Song, 2016;.
Usually, the resolution of exposure models is constrained independently of the hazard, and to a certain extent, also independently of the exposure distribution. That might lead to very poor exposure resolutions where it really matters, i.e. in areas where buildings are densely distributed and/or hazard intensities highly vary with short distances. Or, in contrast, to the 65 unnecessary computation demands for loss assessment in areas with few exposed assets. If aggregation areas of the exposure model are coarser by resolution than the typical correlation lengths of low spatially correlated hazard intensities, then local variations of these intensities would remain hidden in the vulnerability analysis, propagating the associated uncertainties up to the loss estimates.
This dependency between exposure resolution and spatial correlation of hazard intensities has been usually disregarded. Some examples of these practices can be found in soil liquefaction risk. Despite the hazard component can be successfully downscaled 70 (e.g. Bozzoni et al., 2020), thematic uncertainties related to visualisation and interpretation of risk metrics can arise if they are mapped over larger regional administrative units (e.g. Yilmaz et al., 2021) instead of being represented at more hazard-compliant resolutions (e.g. Bozzoni et al., 2021). Similarly, despite that the building exposure models for flood and earthquake vulnerabilities can be aggregated at moderate resolutions (e.g. 4x4 km grid in Dabbeek and Silva, 2019), similar thematic uncertainties can evolve https://doi.org/10.5194/nhess-2021-70 Preprint. Discussion started: 16 March 2021 c Author(s) 2021. CC BY 4.0 License. due to the profound differences between both spatially correlated hazard intensities, and when the calculated losses are mapped 75 over regional administrative units (Dabbeek et al., 2020).
To the best of the authors' knowledge, consideration of different hazard footprints and spatial correlation of their intensity measures for the construction of aggregation entities for exposure modelling has been only seldomly discussed in the literature.
For instance, Chen et al., (2004) described the importance of ensuring a consistent delimitation of the resolution of the exposure model along with the spatial variation of the two considered hazards, earthquakes and hailstorms, that impose damage footprints 80 of very different extents. Complementary, Douglas (2007) and Ordaz et al., (2019) have highlighted the importance of the geographical scale to represent the building exposure models that are affected differently depending on variable hazard footprints.
The study reported in Zuccaro et al., (2018) is perhaps the most advanced framework in the state of the art for the construction of a common aggregation entity for multi-hazard risk assessment, referred to as the minimum reference unit (MRU). This geographical unit coincides with the minimum resolution of analysis of input (i.e. hazard intensities and exposure model) and 85 output elements (i.e. damage and loss estimates) and remarks that despite high-resolution hazard models, one would achieve neither an accurate risk assessment nor meaningful loss mapping if there is no compatibility between the cartographic representation of the building exposure model, the hazard footprints and their IM correlation.
A denser set of geocells in the same area occupied by a coarser regular-sized cell or administrative units provides a denser arrangement of hazard intensity values (when there are local IM variations) to the set of fragility functions per considered hazard. 90 When local IM variations are not sufficiently represented into finer aggregation entities during the vulnerability analysis, thematic uncertainties might appear in the mapping, visualisation and interpretation of the loss estimates. Therefore, besides the conventional epistemic and aleatory uncertainties linked to the hazard, exposure and vulnerability components, thematic uncertainties are also present in the risk chain when the loss metrics are mapped. Awareness of the thematic uncertainties as well as clear and meaningful vulnerability/ loss mappings towards the most relevant hazards a community is exposed to is necessary to 95 improve urban planning, preparedness, mitigation strategies and emergency response actions (e.g. Pang, 2008;Aguirre-Ayerbe et al., 2018).
We can distinguish two types of approaches formerly proposed in the literature to investigate the exposure aggregation for natural hazard risk applications.
(1) To independently represent the building portfolio over a series of aggregation entities such as administrative units, or 100 equidimensional regular grids, and explore their individual contribution to the uncertainty in the losses imposed by certain hazard(s) to ultimately select a representative aggregation model. This option has been explored in Bal et al., (2010), Frolova et al., (2017, Senouci et al., (2018) and Kalakonas et al., (2020) for seismic vulnerability applications, and in Figueiredo and Martina, (2016) for flood vulnerability. These studies discuss the weakness of the physical vulnerability mapping at the individual building scale and over coarse aggregation areas, and highlight the importance of finding an optimal resolution for 105 building exposure modelling while minimizing the uncertainties in the losses estimates. However, these attempts did not explicitly address the spatial correlation or attenuation of the hazard intensity into the predefined aggregation areas, and focused on individual hazards rather than on multi-hazard applications.
(2) Aggregating the exposure models over variable resolution entities that are not necessary administrative boundaries. This has been done in fewer studies. For instance, Muis et al., (2016) assessed the global population exposure to coastal flooding (from 110 storm surges and extreme sea levels) through the application of a hydrodynamic model based on unstructured grids to ensure sufficient resolution in shallow coastal areas. Scheingraber and Käser, (2019) explored the uncertainty in regional building portfolio location for seismic risk through the use of weighted irregular grids. This weighting was carried out as a function of the population density and not using any hazard IM or footprints. Scheingraber and Käser, (2020)  the former procedure in terms of computationally efficiency and treatment and communication of uncertainties in probabilistic 115 seismic risk assessment on a regional scale. Alternatively, aggregating the building portfolio into anisotropic CVT-based geocells (Central Voronoidal Tessellations) is suggested by Pittore et al., (2020).
In this study, we have decided to aggregate the residential building exposure models following the latter aforementioned option, anisotropic CVTs. We are adapting and customizing their derivation to explicitly account for the combination between a lowcorrelated hazard intensity (tsunami inundation) and one exposure proxy (population density) as generating distributions of the 120 centroids of the Voronoi regions. The aggregated building portfolios therein are used for earthquake and tsunami scenario-based risk applications. The use of risk-scenarios has increasingly shown to be useful to estimate physical and cascading damages from natural hazards (e.g. Levi et al., 2018), as well for risk-communication purposes (e.g. Birkmann et al., 2020). We have systematically investigated several risk scenarios for the residential building stock of Metropolitan Lima (Peru) classified in terms of one set of earthquake vulnerability classes and two sets of tsunami vulnerability classes. These building portfolios have been 125 aggregated within six customized CVT models and administrative units at the highest resolution available (i.e. block level).
Through the use of the respective set of fragility functions, we have independently calculated the direct losses from scenario-based physical vulnerability analyses (seismic ground-shaking and tsunami inundation). For such an illustrative application, we have considered six megathrust subduction earthquakes and their respective tsunamis with moment magnitudes ranging between 8.5 and 9.0. We are successively showing that the implementation of this approach is beneficial not only in finding a balance between 130 accuracy and computational demand but also in the efficient representation of the loss estimates while reducing bias generated in the loss mapping. The role of the spatial correlation of both hazard intensities in the efficiency and accuracy of the CVT-based exposure models is also discussed. Since the main scope of this work is to investigate an efficient manner to aggregate the building exposure for risk applications to multiple hazards, we have not investigated the conditional probability relating to cascading events (e.g. , instead, we have assumed that every seismic rupture produces a tsunami. Hence, we are not accounting 135 for cumulative damage on buildings due to consecutive ground shaking and tsunami (e.g. Park et al., 2019;Negulescu et al., 2020;Goda et al., 2021) nor the risk to other seismically induced hazards (i.e. earthquake-triggered landslides, liquefaction, ground failure; e.g. see (Daniell et al., 2017) 2 Methodology to construct variable-resolution exposure model for risk assessment to multiple hazards The proposed methodology is composed of the following steps: 140 1. Simulation of scenario-based hazards (i.e. earthquakes and tsunamis) with the same spatially distributed intensities required by each fragility assessment.
2. Construction of one (or a set of) representative focus map(s). This implies to select and rank (with numerical weights) the hazard intensities or exposure proxies.
3. Generation of CVT-based aggregation entities employing the focus map as an underlying generating distribution and with 145 different numbers of seeding points. 4. Classification of the exposed building stock of interest into vulnerability classes per considered hazard and their aggregation into the CVT-based geographical entities. 5. Scenario-based risk assessment independently per hazard type and discussion of their associated thematic uncertainties in the loss mapping and visualisation. 150 The uncertainties in steps 3 and 5 are explored through the formulation of a condition tree. https://doi.org/10.5194/nhess-2021-70 Preprint. Discussion started: 16 March 2021 c Author(s) 2021. CC BY 4.0 License.

Simulation of scenario-based hazards with spatially distributed intensities
We employ numerical earthquake and tsunami scenarios to replicate historical or hypothetical future events to simulate spatially distributed hazard intensities. For earthquakes, we simulate ground motions from suitable GMPE (ground motion prediction equations). Cross-correlated ground motions are generated for the spectral-periods that serve the fragility functions as intensity 155 measures (IM). For tsunamis, we employ the physical generation and propagation model TsunAWI (Harig et al., 2008) and simulate coastal inundation as IM for tsunamis. The spatially distributed tsunami intensity values (inundation flow depth) are compatible with the IM of the selected set of fragility functions required in the vulnerability analysis.

Construction of focus maps
The focus map drives the construction of a variable-resolution exposure model for aggregating building portfolios. Focus 160 maps were first introduced by Pittore (2015) Through the use of a pooling operator, a focus map highlights areas where a weighted combination of various normalized spatially distributed indicators ( ) jointly assume the larger values. We propose to obtain a focus map that drives the aggregation entities 165 for earthquake and tsunami exposure modelling through the combination of two indicators. (1) Population density ( 0 ) (from aggregated data sources e.g. WorldPop; GPWv4 (CIESIN, 2018)). This indicator is an exposure proxy about the location of residential buildings for which their ground-shaking vulnerability should be addressed. The use of the latter can be a useful indicator when other seismic risk components such as soil amplification conditions are poorly known; come from proxies with coarse resolutions (e.g. topography-based); or when strong seismic site effects are not expected. (2) The tsunami component is 170 constrained through the expected tsunami inundation height ( 1 ) obtained from a "worst-case scenario" approach (i.e. largest feasible intensities) among a series of deterministic scenarios (e.g. Omira et al., 2009;Wronna et al., 2015). For the combination of the two aforementioned normalized input layers, we use a log-linear pooling operator of Eq.
(2). This algorithm assigns a higher sampling probability to spatial locations where both indicators are relevant while penalizing the locations where at least one of the indicators (i.e. tsunami) is negligible. 175 Where Z is a normalizing constant and w i represents the respective weight assigned to score the relevance of each input layer, and ∑ = 1, > 0 ∀ . Thus, the construction of a focus map entails the selection of the weights that rank the importance of every layer, as such carries its own epistemic uncertainties.

Generation of CVT-based exposure models
Selectively increasing the spatial resolution in aggregation areas is beneficial to capture the low spatially correlated hazard intensity 180 values such as tsunami inundation heights. This is achieved by the construction of geocells with variable resolution in the form of Central Voronoi Tessellations (CVT), as proposed by Pittore et al., (2020). During this construction, focus maps are used as underlying generating spatial intensities and are sampled using a Poisson point process (Cox and Isham, 1980) to generate a number of seeding points. These points will be used as centroids of the Voronoi geocells and through an iterative relaxation process will https://doi.org/10.5194/nhess-2021-70 Preprint. Discussion started: 16 March 2021 c Author(s) 2021. CC BY 4.0 License.
converge to the final CVT geocells. The number of seeding points is therefore defining the number of geocells of the resulting 185 tessellation. CVTs are computed in various iteration steps using the simple relaxation method originally proposed by Lloyd (1982) until the distance between the geometrical centroid of the geocell and the weighted mass centroid generated by the raster distribution falls below a defined threshold, or after a given maximal number of iterations achieved. Since the relaxation process is based on the underlying focus maps as generating distribution, CVT cell sizes are inversely proportional to the intensity of the focus map. Each CVT geocell by fact becomes a minimum resolution unit as proposed in Zuccaro et al., (2018), and the resulting 190 tessellation sets the basis for a variable-resolution exposure model.

Condition tree for multi-hazard exposure modelling
The epistemic uncertainties underlying the two former steps are explored using a condition tree with the following hierarchical levels: i.
Selection of a suitable scheme (sets of building classes and their associated fragility functions) to describe the building 195 inventory in the study area for risk assessment.

ii.
Weight arrangement values ( ) that rank every input layer (low spatially correlated hazard intensities or spatial proxies related to the exposure component) in the focus map construction.

iii.
Determination of the number of seeding points sampling the Poisson Point Process driven by a focus map that drives the generation of CVT-based geocells. 200 The condition tree presents a summary of assumptions for uncertainty treatment (Beven et al., 2018). Through the construction of alternative multi-resolution exposure models, the impact of every level of the conditional three is systematically investigated once the vulnerability assessment is performed.

The classification of the building stock into vulnerability classes and aggregation
The building stock of interest is classified in terms of several sets of mutually exclusive, collectively exhaustive vulnerability 205 classes, whose aggrupation describe a set of classes (scheme) specific to the considered hazard (i.e. earthquake and tsunami). A top-down approach is used making use of aggregated census data and ancillary data for the seismic-oriented building classes.
Subsequently, the proportions assigned to each seismic-oriented building classes are reassigned to tsunami oriented classes through the use of inter-scheme compatibility matrices as recently presented in Gomez-Zapata et al., (2021). This method is partly based on the taxonomic disaggregation of every building type within a source and a target scheme as proposed by Pittore et al., (2018) 210 for seismic vulnerability applications, as well as strongly suggested by Charvet et al., (2017) for the definition of tsunami oriented building classes. Then, the classified building stocks are aggregated into every CVT model obtained in the former step.

Scenario-based risk assessment
The fragility of the building portfolio to the considered earthquakes and tsunami scenarios is calculated separately over every aggregation exposure models obtained in Sect. 2.1. This decision is based on the recent findings of Petrone et al., (2020) who found 215 fundamentally different structural responses to both perils. As a consequence, the authors argued that the intensity of the seismic ground motion does not play a significant role in the building's structural tsunami response unless it induces structural yield. The latter is assumed for the vulnerability analysis considering the sole objective of this study of evaluating an optimal exposure model for risk assessment from the considered hazards. The scenario-based risk assessment makes use of a set fragility function associated with each building class and whose IMs are compatible with the hazard intensities modelled in Sect. 2.1. Each of their damage 220 states has an assigned loss ratio to total replacement cost per hazard-dependent building class. According to Petersen et al., (2018), Peru, among all the South American countries, has the largest number of inhabitants, which 225 with a 10% probability of exceedance in 50 years, may experience a ground-shaking greater than VIII (modified Mercalli intensity scale, MMI). This makes Peru the country in which the largest average annual fatalities from earthquakes are expected in South America. In the same study, Lima, with nearly 10 million inhabitants has been identified as the capital city exposed to the highest seismic hazard in the region. This is in line with Salgado-Gálvez et al., (2018). Earthquakes are mostly attributed to the oceanic Nazca Plate subducting beneath the South American plate (Villegas-Lanza et al., 2016). In 1746 a giant megathrust earthquake 230 with an estimated magnitude of Mw 8.8 (Jimenez et al., 2013) ruptured along some 350 km and triggered a tsunami with local runup heights of 15 to 20 m (Dorbath et al., 1990). These events completely destroyed Lima and Callao. The less studied events of 1586 and 1724 triggered tsunami run-ups over 24 m (Kulikov et al., 2005). Moreover, according to Løvholt et al., (2014), Peru has the largest population exposed to tsunami hazard in the American continent, with Lima representing ~ 1/3 of the total country's population. The most important economic activities as well as the main port of the country are located in the study area. 235 Consequently, in Schelske et al., (2014), Lima has been ranked as the second metropolitan area in the world in terms of the value of working days lost relative to the national economy due to earthquakes. This highlights the relevance of integrated vulnerability studies in this study area for the overall economy of Peru.
Local authorities have conducted studies for emergency management and recovery planning considering tsunami and earthquake scenarios (e.g. PREDES, 2009) including qualitative risk estimations. Among others, the Japanese SATREPS project 240 contributed to the improvement of the exposure model of Lima using satellite imagery and census data (Matsuoka et al., 2013).
On the seismic vulnerability side, the project SARA, led by the Global Earthquake Model (GEM), made a large contribution to classify the residential building stock of Peru (Yepes-Estrada et al., 2017) and associate a fragility function per class (Villar-Vega et al., 2017). More specific models for confined masonry have been reported in Lovon et al., (2018). On the tsunami vulnerability side, Adriano et al., (2014) estimated tsunami damage probabilities for two tsunami scenarios over the residential building portfolio 245 classified into four building classes and employing the empirical tsunami fragility functions developed by Suppasri et al., (2013) after the 2011 Japan tsunami. Ordaz et al., (2019) analysed the probabilistic physical risk to both hazards in Callao.

Construction of earthquake and tsunami scenarios for Lima
We have simulated six earthquakes and tsunami scenarios off-shore Peru moment magnitudes Mw 8.5 to 9.0. Finite fault ruptures 250 are modelled using the OpenQuake engine (Pagani et al., 2014) emulating the historical earthquake that occurred in 1746. The basic parameters used in the simulation are; hypocentre location (longitude = -77.93°; latitude = -12.19°; depth = 8 km); strike = 329°; dip = 20°, rake = 90°. Replicating the 1746 earthquake is in line with former studies (e.g. Mas et al., 2014;Pulido et al., 2015;Ceferino et al., 2018a). Spatially distributed ground motion fields are generated using the GMPE proposed in Montalva et al., (2017). Its site term is based on the shear wave velocity in the uppermost 30 meters depth (Vs30) reported in Ceferino et al., 255 (2018b) in which the slope-based Vs30 values (Allen and Wald, 2007) and a seismic microzonation (Aguilar et al., 2013) were compiled and merged to the same resolution (30 arc-seconds ~ 1 km). The aleatory uncertainty in the ground motions is addressed by generating 1,000 realisations per event, as advised in Silva (2016)

Construction of the focus maps
Focus maps have been constructed as inputs to generate CVT-based aggregation boundaries for the building exposure model for 285 seismic and tsunami risk assessment. The spatial population density (PD) in Lima at the block level (INE, 2017) has been combined with a "worst-case" scenario of tsunami inundation height (TI). Distribution of the ground motion has not been used, because, besides its high spatial correlation, negligible nonlinear soil effects in the study area (Pulido et al., 2015) would not allow the visualisation of seismic-intensity-driven hot-spots within a focus map. A Mw 9.0 tsunami scenario was selected among a catalogue of 1,000 offshore scenarios for Lima/Callao that range by magnitude from Mw 8.0 to Mw 9.0. Both map layers have been linearly 290 normalized and combined using the log-linear pooling expressed in Eq. (2) in order to assign a higher probability to the spatial locations where both indicators are relevant. Two sets of weights that rank and combine the layers have been selected to perform a sensitivity analysis at this step. In both sets, tsunami intensities were ranked higher as population density due to their lower spatial correlation. The following weights were accepted for the construction of two focus maps: set (1) PD = 30% -TI = 70%, and set (2) PD = 40% -TI = 60%, and create the basis for the condition tree presented in Figure 3. The construction of the focus map for 295 the first set is shown in Figure 4-a.

Generation of CVT-based exposure aggregation boundaries.
Three seeding sets have been generated by sampling the heterogeneous Poisson point processes defined by the two focus maps including 5,000, 10,000 and 50,000 initial points. We have obtained six CVT aggregation entities for residential building exposure 300 modelling by applying the Lloyd relaxation method as described in Pittore et al., (2020) and recalled in Sect. 2.3. As an example, the resultant CVT-based model obtained from the focus map from set (1) PD = 30% -TI = 70%, and the 5, 000 seeding points (model PD30_TI70_5,000) is depicted in Figure 4-b. The area jointly exposed to the earthquakes and to the largest tsunami footprint (Mw 9.0) is highlighted in pink colour. The distribution of Vs30 values in the study area within this CVT-based example is shown in Figure 5. Due to the contribution of the population density layer (PD), for every Vs30 value at each location, there is a 305 higher density of IM values that are computed where the exposed assets are expected to be concentrated rather than in the locations less densely populated, in contrast with what would occur using a regular grid.

Classification of the building stock into vulnerability classes and aggregation
The residential building stock of Metropolitan Lima (Peru) has been classified in terms of one scheme oriented to seismic 310 vulnerability and into two tsunami schemes with related building classes. They have been constructed following Sect.

Results: scenario-based risk assessment
As suggested by Petrone et al., (2020) due to the fundamentally different structural responses to both perils, the direct economic losses of the aggregated building portfolios for six scenario earthquakes and the corresponding tsunamis have been separately estimated. Hence, all effects related to cumulative damage have been disregarded. The variability of aggregation areas that form every residential building exposure model of the entire Lima/Callao is depicted in Figure 8-a, and listed in Table 1-a. Conversely,  365 if we narrow down the exposed area to the largest tsunami footprint (Mw 9.0), we see that the variability in aggregation areas is very different (Figure 8-b). The CVT-based models with higher resolutions (50,000) geo-cells can reach very fine areas when the focus map considers the weights PD = 30% -TI = 70%. Whilst the block model can reach the largest area values. Moreover, the https://doi.org/10.5194/nhess-2021-70 Preprint.

Seismic risk 380
Seismic losses for the entire study area are initially presented for a selected Mw 8.8 earthquake scenario in order to discuss the implications of the resolution of the exposure model in the economic loss estimates as well as on their associated mapping and visualisation. A comparison for the other five earthquake scenarios is provided in section 3.6.3 for the commonly exposed area to ground-shaking and tsunami inundation. As formerly described, the residential building stock of Lima is classified in terms of the SARA scheme and aggregated considering eight different geographical models (six CVT-based, one block, and one district-based 385 model). Each building class has an associated analytically derived fragility function provided in Villar-Vega et al., (2017)  The seismic vulnerability analysis is performed at every geocell-centroid, where the buildings are aggregated. We consider 390 each IM value resulting from 1,000 realisations of spatially cross-correlated and uncorrelated ground motion fields. The resultant distributions for the Mw 8.8 scenario are displayed in Figure 9. Uncorrelated ground motion fields led to very homogeneous distributions, except at the district level. This finding is aligned with the recent study presented by Scheingraber and Käser, (2020).
Moreover, the latter confirms that if the dimension of the geocells in the exposure model is larger than a typical seismic ground motion correlation length (i.e. 20 km), an artificial bias in the ground motion correlation has to be expected as described in Stafford 395 (2012). We obtain larger median loss values from uncorrelated ground motions. We observe that for the considered scenario, the median loss values are insensitive to the aggregation of the exposure model at varying resolutions. This feature was already described in Bal et al., (2010) for a crustal earthquake damaging a building portfolio in Istanbul while neglecting the crosscorrelation model. We are thus confirming the finding reported in Bal et al., (2010) while further expanding it when a ground motion cross-correlation model is actually considered. Despite these remarkable differences between the areas distributions of the 400 models ( Figure 8, Table 1), we do not observe significant differences in the associated absolute seismic-induced losses, which is explained by the high special correlation of seismic ground motion. The financial loss results that we have obtained are similar to the loss distribution estimated by Markhvida et al., (2017) who reported mean loss values ~ 7 and maxima ~ 35 billion USD. That study investigated the possible losses of the residential 410 building stock of Lima/Callao (aggregated into a regular grid (~ 1 km 2 )) expected from a similar Mw 8.8 scenario. Although the authors employed a different GMPE with respect to the one we have adopted, the ground motion cross-correlation model as well as the set of building classes and fragility functions are the same we have implemented. In Figure 10 the spatial distributions of losses for the considered Mw 8.8 scenario are shown for the "La Punta" sector for two selected ground motion realisations. In accordance with Vamvatsikos et al., (2010) and similar to other studies entailing spatial 415 comparison of risk estimates (e.g. Senouci et al., 2018), the loss values are normalized. Although the absolute value in the loss estimates is not affected by the aggregation entities (Figure 9), large differences arise when the normalised losses are actually mapped. It can be noted as for the same realisation, regardless of the use of correlated or uncorrelated ground motions, the seismically vulnerable areas are still identifiable, albeit with considerable differences. The use of cross-correlated ground motion fields results in smoother mappings. However, the component which imposes the largest impact on the loss estimated from scenario 420 earthquakes is the simulation of the seismic process, as remarked in other studies (e.g. Silva, 2016). This further highlights the importance of using quantile analysis in mapping seismic risk estimates for better visualisation and communication of the uncertainties in an inherently stochastic process (Geller, 2015). They are mapped over six aggregation areas of the building portfolio classified in terms of the SARA vulnerability classes. This is done for two randomly selected realisations with uncorrelated ground motion fields (b, d) and cross-correlated ground motion fields (a,c) using the Markhvida et al., (2018)

Tsunami risk
To constrain the economical consequence model used in tsunami risk assessment, the inter-scheme conversion matrices depicted in Figure 6 are used to obtain the replacement cost values per building class from the correspondent maximum scoring class in SARA. We have assumed loss ratios 5%, 15%, 45%, 65%, 85%, and 100% for each of the six damage states proposed by Suppasri 435 et al., (2013) and similarly, but starting with 15%, for the five ones proposed in De . A similar approach has been recently adopted by Antoncecchi et al., (2020). The impact of using more exhaustive approaches (e.g. (Suppasri et al.,  worth exploring, but out of the scope of this paper. Both tsunami-vulnerability schemes have been associated with a set of empirical fragility functions with tsunami inundation height in meters as IM. They were derived from the same building damage dataset collected after the great 2011 Mw 9.1 Japan earthquake and tsunami with damage state definitions that implicitly accounted for 440 the combined effect of both hazardous events. Thus, despite the extensive use of empirically derived fragility functions from that specific near-field event, care should be taken when using them, not only because they also account for the ground-shaking induced damage, but also because a submarine landslide could have contributed to the tsunami (Tappin et al., 2014). Nevertheless, there is a profound difference in the way the mean intensity values were obtained. Whilst in Suppasri et al., (2013), a linear least squares regression fitting was carried out, in De Risi et al., (2017) a multinomial logistic regression was performed for material-based 445 classes. The latter found very similar regression values to the case when an average simulated flow velocity of 1.84 ms -1 for masonry and wooden buildings classes (predominant in Lima) is integrated into an hybrid fragility model. Making use of these two schemes, we have correspondingly estimated the tsunami-induced losses for the six scenarios and for the seven exposure models. Tsunami loss estimates normalised to the losses by the block-model are presented in Figure 11. Independent of the reference scheme, the two CVT models with the largest number of geocells (50,000) show the closest similarity to the block model 450 (normalized ratio ~1). However, for all the CVT models, this ratio dramatically drops for scenarios with lower magnitudes (8.5, 8.6, and 8.7) which can probably be explained by a smaller tsunami footprint and lower IM spatial correlation. The absolute loss values expected after the six tsunami scenarios are reported in Figure 12 at the block level for the two vulnerability reference schemes. The fragility models of Suppasri et al., (2013) Song et al., (2017). These studies concluded that flow-depth models predict larger probabilities of complete damage for buildings than models which that addressed tsunami velocity in their derivation.
We have computed the discrepancy in the tsunami loss estimations obtained for each CVT model with respect to the block-based model (Figure 13). This is minimised for the larger magnitudes and higher resolution models (50,000 geocells). This analysis shows that the Suppasri et al., (2013)    Estimates at the block level are ~ 3 times less efficient (i.e. larger size, in Table 1) compared to the CVT-models with 470 50,000 geocells for which almost identical results are obtained (Figure 11), and ~ 23 times larger than CVT-models with 5,000 geocells for which overestimations lower than 20% are expected for the three largest magnitudes ( Figure 13). Aggregation entities with lower resolutions lead to overestimations in the tsunami-induced losses. Similar findings were reported in Figueiredo and Martina (2016) in flood risk analysis.
Tsunami loss outcomes for the Mw. 8.8 scenario are mapped and discussed hereafter for the residential building stock in 475 the two areas "La Punta" (Figure 14) and Chorrillos district ( Figure 15). Only geocells with loss values larger than zero are colourmapped. Due to the normalised metric used, no significant differences in the tsunami vulnerability mapping induced by the independent building classification schemas are noticeable. The CVT models at the coarser resolutions (first two rows in every figure) show the largest values, and hence overestimations compared to the block level and other finer CVT-models.
Overestimation of losses decreases with the increase of resolution. Due to the adjacency and compactness of the highest resolution 480 CTV model (fourth row), for "La Punta" we identify at least four zones with a comparatively higher tsunami vulnerability. Considering Figure 15, it can be noted that the overall mapped area is increasingly reduced as the resolution of the CVT models 485 increases. This is due to the lack of residential buildings within the three large parcels, namely "Country Club de Villa", "Reserva Laguna de Villa" and "Refugio de Vida Salvaje Pantanos de Villa" that occupy most of the exposed area in the block-based model.
These zones represent the largest area values in Figure 8-b. This model assigns the largest loss values in the Chorrillos district to these three large blocks due to the assumption of using a single tsunami intensity as representative of the entire enclosing area. Therefore, if the block polygons are too coarse compared with the hazard footprint and IM spatial correlation, biases in the loss 490 assessment are expected. This is highlighting the importance of hazard-driven entities for exposure spatial aggregation.

Comparison between earthquake and tsunami scenario-based induced losses
In Figure 16 we compare the absolute losses induced by each hazard scenario onto the building portfolio exposed to both perils (e.g. Mw 9.0 in Figure 4-b). The CVT-based PD30_TI70_5,000 is used to represent the earthquake-induced losses. The latter was compiled for the cases with and without ground motion cross-correlation model, each sampled with 1,000 realisations. Due to the 500 lack of stochastic realisations in the tsunami case, the respective loss distributions were constructed for both reference schemes with the seven values obtained from the various aggregation entities (6 CVT-and 1 block-based models). Despite the fact that the distributions for seismic and tsunami losses have been obtained independently, the median values are nevertheless illustrative for comparative purposes. Expected losses from earthquake and tsunami scenarios on every commonly exposed area for residential buildings in Lima 505 Figure 16. Comparison of the independent earthquake-induced losses (EQ) for two ground motion field conditions (with and without ground motion cross-correlation model) and the tsunami-induced losses (TS) under two tsunami reference schemes for six magnitude scenarios over every common area exposed to both perils.

510
We observe that in our estimations for the commonly exposed area to both perils in Lima, the earthquake event dominates the median losses at lower magnitudes (Mw 8.5,8.6) whilst the tsunami prevails in the larger ones. The tsunami-induced median losses start to be larger than the earthquake-related ones for the Mw 8.7, although the latter still present high volatility in the extreme values due to the variability in the seismic realisations. From Mw 8.8 on, the tsunami-induced losses are always larger regardless of the tsunami reference scheme implemented. Our findings regarding the role of the earthquake magnitude in the disaggregation 515 of financial loss estimates for every hazard scenario are in line with the results of Goda and De Risi, (2018) obtained for a coastal town in Japan jointly exposed to two decoupled earthquakes and tsunamis risk scenarios (Mw 8.0, 9.0).

Discussion
CVT-based aggregation entities for building exposure modelling can be further customized. For instance, the underlying focus map can be modified in order to integrate other components such as higher resolution seismic microzonation; spatially explicit 520 distributions of the duration of the expected hazard scenario, key aspect in hydraulic risk assessment (Ward et al., 2020); the spatial presence of certain building taxonomic attributes that may drive the physical vulnerability towards a given hazard (e.g. soft-storey in seismic vulnerability and openings or building foundation in tsunami vulnerability (Alam et al., 2018)); high-resolution DEM (digital elevation model). However, caution should be taken not to double count their contribution if the hazard simulations have already been performed using these input data (e.g. DEM in landslide susceptibility and tsunami inundation) as well as a wise 525 selection of their respective weights in the focus map construction. Furthermore, CVT model generation would benefit from further improvements such as outlining an iterative approach that can seek for a minimum geocell size from a convergence criterion imposed by one or several spatially correlated hazard-IM-lengths.
Several of the limitations in this study could be addressed in further studies. For instance, including the computation of an exhaustive set of stochastic tsunami inundation scenarios (per considered magnitude) for loss assessment, as well as including 530 Physics-based tsunami fragility functions based on intensities more relevant to the building failure mechanisms such as momentum flux (e.g. Macabuag et al., 2016;Attary et al., 2017) would benefit future risk simulations. However, this is subject to their actual availability for the typical Peruvian building classes. More comprehensive adaptations of such as "foreign" fragility models (e.g. Suppasri et al., 2019) in Peru, as well as the need of future development of analytical models for the South American context (e.g. 535 (Medina et al., 2019) would benefit future risk assessment studies for Lima. Another area that would benefit from future research is the differential selection of loss ratios with dependencies on the building classes, as for instance recently investigated by Kalakonas et al., (2020) for seismic risk applications. This might be also relevant for tsunami-induced losses that are strongly influenced by the presence and cost of non-structural building elements. Accordingly, more refined financial tsunami consequence models such as the one proposed by Suppasri et al., (2019) and/ or Triantafyllou et al., (2019) are worth exploring when detailed 540 information about prices and built-up areas at the individual building level are available for the study area. In the presented example case, the tsunami economic consequence model was entirely based on the replacement costs of the SARA model obtained from the inter-scheme conversion matrix. On the one hand, the latter makes that the exposed residential assets classified under various schemes have equivalent replacement costs, and thus, the hazard-dependent risk estimates can be comparable with each other. On the other hand, caution should be taken when interpreting the presented results. Neither the damages induced by debris impacts 545 nor scour, relevant for a clearer tsunami vulnerability assessment (Charvet et al., 2015), are included in our modelling. Moreover, it is worth to mention that larger indirect losses can be expected from buildings with other occupancies (e.g. Chen et al., 2018) that we have not considered herein.
CVT-based models can be beneficial to define efficient, multi-hazard aggregation entities for earthquake and tsunami risk assessment, not only in Lima, but also in other coastal cities exposed to similar hazards. Furthermore, it is worth investigating the 550 usefulness of in mapping cumulative damage and losses in hazard sequences i.e. when a first hazardous event modify the fragility of buildings that are then affected by a successive event.

Conclusions
The focus map concept formulated by Pittore (2015) has been used to integrate spatially correlated hazard intensity measures (IMs) with exposure proxies (i.e. population density) to spatially identify hot-spot areas where higher values from both distributions are 555 expected. A focus map can then be sampled by a heterogeneous Poisson point process, as proposed by Pittore et al., (2020) in order to generate variable-resolution aggregation entities in the form of a Central Voronoi Tessellations. The proposed method is subject to epistemic uncertainties, namely the selection of weights for pooling of geospatial layers into the focus map and the selection of seeding points that provide initial seeding set and control the overall number of geocells. Through the application of a condition tree, we have selected different sets of these two components to investigate the impact of customised CVT-based geographical 560 entities used to aggregate building portfolios. This is done through independent risk loss estimates per hazard as well as through the reconnaissance of thematic uncertainties in the loss mapping and visualisation.
Variable-resolution CVT-based exposure models proposed in this work have proved their efficiency in integrating largearea building portfolios for combined earthquake and tsunami loss estimations. Several advantages over conventional models based on administrative aggregation entities are: 565  CVT-based models provide alternative an approach to aggregate an extensive building portfolio constructed from ancillary data (i.e. population) in the case when existing administrative aggregation areas are not suitable (either not publicly available, https://doi.org/10.5194/nhess-2021-70 Preprint. Discussion started: 16 March 2021 c Author(s) 2021. CC BY 4.0 License. or too coarse in resolution) for a certain area of interest, as well as to perform scenario-based risk assessments from various hazards.
 We have observed that they correct some bias in the spatial aggregation of buildings due to the smaller, more compact areas 570 in high-resolution CVT geocells with respect to a coarser, block-based cell. This correction is further propagated to the loss estimates due to the higher density of IM values employed by the respective fragility functions during the loss assessment. This is especially observed in areas of the largest concentration of exposed assets located within the hazard footprint area and where local spatial variations of the IM are expected, leading to more accurate estimates.
 They are computationally more efficient than the block-based models in earthquake and tsunami vulnerability assessment. 575 This is advantageous when thousands of stochastic realisations of hazard scenarios are calculated over the aggregation boundaries that are used to model building portfolios.
 They have shown to be beneficial to map loss estimates in continuous space with adjacent and compact geocells. These features allow the spatial identification of zones with similar vulnerability to the hazards considered and within the area of interest. They contribute to a more intuitive visualisation and interpretation of the loss mapping and hence contribute to 580 raising awareness about epistemic and thematic uncertainties in the loss mapping. That is a significant help for decisionmakers in undertaking urban requalification or disaster risk mitigation activities (e.g. Dolce and Di Bucci, 2015).
For the portfolio exposed to both perils in Lima, we have found that the expected median loss values induced by seismic ground-shaking are insensitive to the representation of the exposure model at varying resolutions. Thus, we are confirming the finding reported in Bal et al., (2010) and expanding it to the case when cross-correlated ground motion fields are considered. 585 However, this is contrasting to tsunami loss models, whose differences with respect to a high-resolution model (i.e. block-based) decrease as the resolution of CVT geocells increases. Similarly, these differences are remarkably minimised for incrementally correlated tsunami-intensities from the large magnitude tsunami scenarios (i.e. Mw 8.8, 8.9, 9.0). According to our observations, the adopted tsunami fragility model based on solely flow-depth as IM and linear square fitting (Suppasri et al., 2013) predicts much larger tsunami-induced losses on the residential buildings portfolios in Lima than the model of De  derived 590 through multinomial logistic regression and with similar values as if flow velocity was accounted. For the residential building portfolio exposed to both perils, we have found that the earthquake scenarios dominate the losses at lower magnitudes (Mw 8.5, 8.6) whilst the contribution of tsunami is dominant for the larger magnitude events. Thus, this study invites to realise the comparative relevance of the tsunami risk in extreme future scenarios in Lima, a coastal mega-city whose safety is vital for the overall economy of the country. 595 Bearing in mind the scope of this study, but also the aforementioned limitations in the discussion section, we are not claiming that the scenario-based economic losses we have presented for the residential building stock of Lima are completely exhaustive. Instead, through the adoption of the condition tree, we have drawn a branched methodological workflow to explore the differential impact of the exposure aggregation models, and the selection of building schemes on the epistemic and thematic uncertainties that are embedded in scenario-based risk applications. As described by Beven et al., (2018), condition trees facilitate 600 the communication of the meaning of the resulting uncertainties while providing a clear audit trail for their analysis that can be reviewed and evaluated by others (e.g. local experts and stakeholders) at a later date. Nevertheless, these aspects are highlighting the importance about keep working on seismic and tsunami fragility models that consider particular construction practices, local hydrodynamics, and remarkably, the relevance of hazard-sound aggregations entities for exposure modelling and loss mapping.
Thus, the continuous understanding of those uncertainty sources is contributing to enhance future risk communication, mitigation, 605 and disaster management activities by local decision-makers.