Modeling spatial dependencies of natural hazards in coastal regions: a nonstationary approach with barriers

Natural hazards like floods, cyclones, earthquakes, or, tsunamis have deep impacts on the environment and society causing damage to both life and property. These events can cause widespread destruction and can lead to long-term socio-economic disruption often affecting the most vulnerable populations in society. Computational modeling provides an essential tool to estimate the damage by incorporating spatial uncertainties and examining global risk assessments. Classical stationary models in spatial statistics often assume isotropy and stationarity. It causes inappropriate smoothing over features having boundaries, holes, or physical barriers. Despite this, nonstationary models like barrier model have been little explored in the context of natural disasters in complex land structures. The principal objective of the current study is to evaluate the influence of barrier models compared to classical stationary models by analysing the incidence of natural disasters in complex spatial regions like islands and coastal areas. In the current study, we have used tsunami records from the island nation of Maldives. For seven atoll groups considered in our study, we have implemented three distinct categories of stochastic partial differential equation meshes, two for stationary models and one that corresponds to the barrier model concept. The results show that when assessing the spatial variance of tsunami incidence at the atoll scale, the barrier model outperforms the other two models while maintaining the same computational cost as the stationary models. In the broader picture, this research work contributes to the relatively new field of nonstationary barrier models and intends to establish a robust modeling framework to explore spatial phenomena, particularly natural hazards, in complex spatial regions having physical barriers.


Introduction
A 2019 United Nations report indicates that the world's population is expected to reach 11 billion by 2100.By then, it is estimated that nearly 75 percent of the population will live in urban areas.According to the scientific community, climate conditions would be very different from what we are currently experiencing.Even by 2030, there will be a 33 percent rise in the overall quantity of urban land in highfrequency flood zones compared to 2015 (Gu ¨neralp et al. 2015).Natural hazards pose significant threats to both the environment and human settlements.These hazards can directly impact geomorphological and hydrological processes, as well as biodiversity, with potentially severe consequences for society (Zorn and Komac 2013;Emmer 2018).Scientific literature indicates the significant impact of climate change and natural disasters on populations that are vulnerable to these events (Briere and Elliott 2000;Correa et al. 2011;Benevolenza and DeRigne 2019).Additionally, there have been researches conducted on exposure to natural disasters in the general population and society (Cannon 1994;Cutter 1996; Zhou et al. 2014;Aksha et al. 2019;Raker 2020).It is therefore imperative to thoroughly evaluate the risks associated with natural hazards in order to make informed decisions regarding spatial planning and risk-prevention measures (Spence 2004;Botzen and Van Den Bergh 2009;Cui et al. 2021).This proactive approach is essential for ensuring the sustainability of ecosystems (Sudmeier-Rieux 2006;Gunderson 2010) and minimizing the impact of natural hazards on both the environment and human settlements (Small and Nicholls 2003;Bornstein et al. 2013;Oliveira et al. 2020) over the coming decades.International organizations, such as the World Bank, the European Union (EU), and the United Nations (UN), among others, are aware of the necessity to take into account the long-term effects of natural disasters.Additionally, according to the intergovernmental panel on climate change (IPCC), ''successful risk reduction and adaptation strategies consider the dynamics of vulnerability and exposure and their relationships with socioeconomic processes, sustainable development, and climate change''.
A large number of existing studies in the broader literature have examined the impact of extreme climatic conditions on natural disasters (Sauerborn and Ebi 2012;Phillips et al. 2015;Chaudhuri et al. 2021;Osberghaus and Fugger 2022;Raju et al. 2022).As a result of unpredictable climate change, population growth, and increasingly urbanized societies, a better understanding and prediction of natural disasters has become undoubtedly important (SafarianZengir et al. 2019).Modeling natural disasters is crucial to characterize these phenomena and provide tools to counteract them.Estimating natural hazards, including spatial effects and local conditions, will help in the management and even allow anticipation of events (Cutter and Finch 2008).The 2015-2030 Sendai framework for disaster risk reduction recognizes this need and emphasizes the significance of having strategies in place to mitigate uncontrolled development in hazardous areas in order to better prepare for the disasters that our planet may experience in the future.
Several studies attempt to model natural hazards in order to examine global risk assessments (Morjani et al. 2007;Serra et al. 2013;Calkin and Mentis 2015;Riley et al. 2016;Pittore et al. 2017;Sarkissian et al. 2020).Some of them address the propagation of tsunami and its impact (Sarri et al. 2012;Hayashi et al. 2013;Shao et al. 2019;Rezaldi et al. 2021;Sugawara 2017).In this line, to analyse the spatial distribution of regions affected by natural hazards, statistical inference comes along with Bayesian methodology.Generally, Bayesian statistics are utilized in natural hazards engineering to deal with large-scale problems that involve different types of data inputs and explicitly handle uncertainties (Zheng et al. 2021).For example, studies like (Gaume et al. 2010;Costa and Fernandes 2017;Han and Coulibaly 2017;Barbetta et al. 2018;Bolle et al. 2018) provide a comprehensive review of the applications of Bayesian statistics in flood assessment and monitoring.Besides, Grezio et al. (2009) have presented the challenges in selecting proper models to quantify the uncertainties in the maximum tsunamigenic magnitudes.Literature shows the application of Bayesian inference in analysing probabilistic tsunami hazards (Knighton and Bastidas 2015;Risi and Goda 2017;Smit et al. 2017).
A Bayesian approach with Markov chain Monte Carlo (MCMC) simulation methods has traditionally been used to fit generalized linear mixed models (GLMM) (Wikle et al. 1998).In this context, Shin et al. (2015) have explored the application of the Bayesian MCMC method to estimate the extreme magnitude of tsunamigenic earthquakes.However, MCMC models require considerable computing time for large datasets (Rue et al. 2009;Smedt et al. 2015).Hence, instead of using MCMC, we can use integrated nested Laplace approximation (INLA) methodology, developed by Rue et al. (2009), as it offers faster computation times and is much easier to fit complex models (Ruiz-Ca ´rdenas et al. 2012).As recommended by Lindgren et al. (2011), while processing spatial data we can utilize INLA in conjunction with stochastic partial differential equations (SPDE) to balance the speed and accuracy of the models.However, there are limited contributions using INLA-SPDE approach in the context of natural disasters such as earthquakes or tsunamis.Recently, Wilson (2020) used Bayesian spatial modeling with INLA to analyse earthquake damages from geolocated cluster data.To date, no literature has documented applications of INLA-SPDE to tsunamis in the particular territory of the Maldives which involves a very advanced methodology to handle its irregular and complex land structure (Riyaz and Suppasri 2016).Analyzing tsunami propagation at the island scale is essential to develop well-informed policies for disaster management and to design effective countermeasures.But due to the large domain and high resolution required for modeling, it is also challenging to study tsunami propagation across multiple atolls at the island scale (Rasheed et al. 2022).
The current study is conducted to model and estimate the spatial autocorrelation of tsunami data in the islands of Maldives.The initial point of this work is to explore the application of SPDE with INLA for Maldives tsunami data (HDX 2022), including SPDE triangulations for spatial effect.Maldives consists of 1200 dispersed islands on both sides of the equator.Collection of these islands along with lagoon and reef areas form the complex atoll system of the country as depicted in Fig. 1.It is worthy to mention that, the traditional SPDE method triangulates the entire study area based on continuous geographic boundaries (Krainski et al. 2018).A problem arises while designing the mesh for the entire country's boundary region or, for boundaries of individual atolls.Although the records of tsunami-affected regions are discrete spatial sites situated exclusively on land on reefs for individual islands, the mesh is generated for the entire region inside the geographical boundary including the lagoon and ocean surface.So, it is not realistic and ambiguous especially when we are interested to explore the spatial correlation of tsunami data precisely on the land regions.This leads to the motivation of designing SPDE triangulation precisely on land on reefs for individual islands.However, a stationary model can not be aware of the coastline and the island boundaries and will inappropriately smooth over the features.In spatial modeling, classical models are unrealistic when they smooth over holes or physical barriers.This might result in another unrealistic assumption.In the research work by Bakka et al. ( 2019) a new nonstationary model has been constructed for INLA having syntax very similar to the stationary model.The model, named as barrier model is more realistic with both sparse data and complex barriers and computational cost is the same as for the stationary models (Bakka et al. 2019).Literature shows, there are limited studies that examine how to use barrier models in the field of spatial statistics, despite the fact that the methodology is practical and scientifically robust for analyzing spatial regions with physical barriers or holes.Islands and coastal regions with their unique characteristics and spatial dynamics, present an ideal setting for the application of nonstationary methods.However, the extent to which these approaches have been studied and applied in coastal regions remains relatively limited.Only a few notable works, such as the studies by Jo ´nsdo ´tir et al. (2019), Martinez-Minaya et al. (2019), Bi et al. (2020) and Jaksons et al. (2022), along with the research conducted by Kaurila et al. (2022), have ventured into exploring this domain.These selected works serve as valuable contributions to the understanding of nonstationary methods in coastal spatial analysis.However, given the vast potential and unique characteristics of coastal regions, further research efforts are essential not only to enhance our understanding of coastal spatial processes but also to unlock valuable insights and potential applications for nonstationary approaches in other domains.It is important to mention that in the aforementioned studies, barrier models have been developed with the assumption that water is treated as normal terrain.These models take into account the specific coastlines and boundaries as physical barriers.In the present study, we have explored the barrier model in a converse mode where water bodies (ocean and lagoons) act as barriers for the dispersed islands and natural hazards are the sample events considered precisely on the land area of the islands.
The motivation and aim of this paper are two-fold.On one side we provide a modeling framework to explore and analyze the spatial variation in the incidence of natural disasters such as tsunamis in complex spatial regions like islands and coastal areas.In the current study, we have assessed the impact of the tsunami on coastal communities by quantifying the number of individuals indirectly affected by the event.By examining the population residing in the coastal regions of Maldives who experienced various consequences resulting from the tsunami, such as displacement, loss of livelihoods, and disruption of essential services, we can gauge the wider repercussions on these communities (Fujima et al. 2006).This approach enables a scientific evaluation of the extent to which the effects of tsunami reverberated beyond the immediate areas of impact, providing a comprehensive understanding of the overall societal impact experienced by the coastal communities in Maldives.In particular, the occurrence of tsunami has been analyzed with three spatial modeling scenarios using mesh for the entire geographical boundary of atolls, mesh precisely on land on reefs for islands in individual atolls, and barrier models for atolls.The second aim roots in providing an advanced and realistic computational strategy to design and customize meshes and nonstationary barrier models to examine the spatial dependencies of natural hazards in complex distributed land structures like islands and coastal areas.We present a comparative analysis between the nonstationary barrier model and two traditional stationary approaches.The first approach utilizes SPDE triangulation for the entire region, while the second approach generates SPDE triangulation precisely on the land on reefs.The evaluation of these models relies on quantifying the consequential impact of tsunamis, specifically in terms of the number of indirectly affected individuals on the land on reefs for individual islands.It can be applied in diverse sectors where complex physical barriers are present in road networks (Dawkins et al. 2021), disease control interventions (Cendoya et al. 2022), and categorical areas with different land use.
The study is conducted on 190 records of tsunami-affected islands from the year 2004.R (version R 4.1.2)programming language (R Core Team 2022) has been used for statistical computing and graphical analysis.As part of the data cleaning process and to design some maps, we have used ArcGIS Pro (version 3.0.1)(Redlands CESRI 2022).All computations are conducted on a quad-core Intel i9-4790 (3.60 GHz) processor with 32 GB (DDR3-1333/ 1600) RAM.
The rest of the paper is organized as follows.Section 2 introduces the study area and offers some insights into the dataset used in the current study, accompanied by comprehensive figures and tables.A description of the modeling framework comes in Section 3. The design and mathematical concepts of the barrier model are discussed in Sect.3.1.Model fitting along with model evaluation techniques are reported in Sect.3.2.Section 4 is devoted to present the results and related discussions of the mesh structures and models.Some concluding remarks come in Section 5.

Data
The Republic of Maldives is located in the southwestern region off the coast of India in the Indian Ocean.This unique island nation is one of the smallest countries in Asia having a chain of coral islands across an archipelago more than 800 kilometers long and 130 kilometers wide.The archipelago consists of about 1190 coral islands grouped into 20 natural atolls.Out of which 189 islands are inhabited (Isles, The Presidents Office 2022).Figure 1 illustrates the geographical location of the Maldives and the complex island structure for the atolls of the country.
The natural hazards dataset of Maldives for the year 2004, contains 190 records of tsunami affected islands, all being inhabited islands, and provides the number of direct and indirect affected people for individual islands.The dataset is published by an open data sharing platform, Humanitarian Data Exchange (HDX) managed by the United Nations Office for the Coordination of Humanitarian Affairs (OCHA) under a creative commons attribution 4.0 international license (HDX 2022).We note that the shape files for atolls and islands of Maldives are accessed from the open data portal1 of the Maldives Land and Survey Authority, Republic of Maldives.These shape files possess a resolution of 1:4000 m.A spatial dataset of natural hazards like cyclones, typhoons, storms, floods, and water shortages of Maldives can be accessed from the same open portal of the Humanitarian Data Exchange (HDX) disaster dataset of Maldives (HDX 2022).We have used tsunami data as a showcase for the current study.The dataset for 190 islands provides detailed information about deaths, injuries, destroyed houses, and also about people who were directly and indirectly affected by the tsunami.In order to simplify the modeling process, we have used the number of indirectly affected people as the response variable.The specific coordinates of tsunami-affected islands are derived from the Humanitarian Data Exchange (HDX) disaster dataset of Maldives (HDX 2022), with the coordinates provided in the World Geodetic System 1984 (WGS84) coordinate reference system.It is worth noting that the shape files utilized in the study also employ the same reference system.
According to an independent evaluation report (2012) by the Asian Development Bank, the 2004 tsunami had a devastating impact on the island nation of Maldives.The country experienced a disaster of national proportions, with 39 islands severely damaged (Asian Development Bank 2012).Figure 2 (left) depicts the locations of tsunami affected islands for individual atolls of the Maldives which include enclosed lagoon or basin, fore reef, sub tidal reef, pass reef flat, and land on reefs.An interactive map of the Maldivian island structure and tsunami affected locations is published in an ArcGIS online map, which can be accessed through this link2 .The present study focuses on assessing the impact of the tsunami on coastal communities in Maldives through a comprehensive analysis of 12 distinct atolls spanning the entire country from north to south.These atolls include Haa Alifu, Haa Dhaalu, Shaviyani, Noonu, Raa, Baa, Kaafu, Meemu, Laamu, Gaafu Alifu, Gaafu Dhaal and Seenu as highlighted in Fig. 2 (right).By strategically selecting these atolls, we ensure that the study covers the entire geographical range of the country, encompassing all 39 islands highly affected by the tsunami.Additionally, the chosen atolls cover key urban areas, such as the capital city of Male ´, as well as other significant cities in Maldives.Atolls with a considerably lower number of tsunami-affected islands are intentionally excluded from this study to maintain a focused and representative analysis.As discussed in Section 1 the atolls of Maldives are a collection of disjoint islands, similarly, almost all the atolls are disjoint land surfaces.But some atolls such as Haa Alifu, Haa Dhaalu, Shaviyani, and Noonu atolls in the north, in the north-west Raa and Baa atolls, and in the south Gaafu Alifu and Gaafu Dhaal, these three regional groups share common boundaries.In the current study, these eight atolls, based on their common geographical boundaries are considered as three distinct combined spatial regions.Details about the regional integration are depicted in Figs. 6, 7, and 8 in the Appendix.In the rest of the study, we have referred to these three atolls groups as, Shaviyani group, Baa group, and Gaafu Dhaal group respectively.Considering these 3 merged atoll groups, the initial 12 unique atolls are now reduced to 7 distinct atoll groups, which are the focal point of investigation in the present study.
These integrations allow us to explore a substantial number of tsunami-hit islands in continuous spatial regions.We report that the present study examines 135 islands from 7 different atolls or groups of atolls as shown in Fig. 2 (right).This is approximately 71 percent of the total 190 tsunami-affected islands of Maldives.Detail records of the number of affected islands and indirectly affected people for the selected atolls are reported in Table 1.
Figure 3 (left) shows the locations of tsunami affected regions of Baa group atolls.Figure 3 (right) depicts detailed island distribution on reef areas for the same.In both cases, tsunami affected regions are highlighted in red.We report that all datasets used in the current study are collected from sources without restrictions and that have open access.

Methodology
Random spatial events generate irregularly scattered point patterns over areas of interest.In these cases, spatial point process models are useful tools to perform precise statistical analysis (Juan et al. 2012;Loo et al. 2011;Karaganis and Mimis 2006).Moreover, we can find recent studies (Verdoy 2019;Opitz et al. 2020) on spatial point processes that are able to identify spatial auto-correlations and interactions between points in the pattern.From Fig. 3 (left) it seems that atolls with enclosed lagoon and land on reefs are continuous land structures.But each atoll is originally a collection of the number of distributed islands as depicted in Fig. 3 (right).By considering the total number of indirectly affected people from individual tsunami-hit islands, we open the door to consider Poisson regression models in combination with a Bayesian framework.Instead of using MCMC, we have used computationally faster solutions for latent Gaussian models by using a Laplace approximation for the integrals with the INLA method (Rue et al. 2009).It focuses on models that can be expressed as latent Gaussian Markov random fields (GMRF) (Rue and Held 2005).
Our approach combines a spatial Poisson regression method with an INLA Bayesian framework.In particular, let Y i and E i be the observed and expected number of indirectly affected victims of the tsunami on the i-th island.We assume that conditional on the relative risk, q i , the number of observed events follows a Poisson distribution: where the log-risk is modeled as Here, n i accounts for the spatially structured random effects and i stands for an unstructured zero mean Gaussian random effect and logGamma precision parameters 0.5 and 0.01, defined as penalized complexity (PC) priors (Simpson et al. 2017).Z i represents the covariates.It is worth noting that, no covariates are included in our study.But it is possible to incorporate relevant covariates into similar models in future studies (Rue et al. 2009).We assigned a vague prior to the vector of coefficients b ¼ ðb 0 ; . ..; b p Þ which is a zero mean Gaussian distribution ).As we run several cases with different priors, we find that our results are robust against other alternative priors.
To bypass the problem of inefficiency in the estimation under a general INLA approximation, we have used another computationally tractable approach based on SPDE models (Lindgren et al. 2011).On one hand, we used SPDE to transform the initial Gaussian Field (GF) with Mate ´rn covariance function to a GMRF.In particular, the spatial random process n, here represented by nð:Þ explicitly denotes dependence on the spatial field, follows a zeromean Gaussian process with Mate ´rn covariance function represented as is the Laplacian operator, a ¼ ðm þ d=2Þ is the smoothness parameter, s is inversely proportional to r, W(x) is a spatial white noise and j [ 0 is the scale parameter, related to range r, defined as the distance at which the spatial correlation becomes small.For each m, empirically derived definition r ¼ ffiffiffiffiffi 8m p =j with r corresponding to the distance where the spatial correlation is close to 0.1 (Blangiardo and Cameletti 2015).Note that we have d ¼ 2 for a two-dimensional process, and we fix m ¼ 1 so that a ¼ 2 in our case.
INLA-SPDE requires a triangulation or mesh structure to interpolate discrete event locations to estimate a continuous process in space (Krainski et al. 2018).In the current study, the spatial coordinates of each tsunami-affected island are employed as the target sites over which we have constructed the meshes using Delaunay triangulation technique.It is important to highlight that the SPDE triangulation techniques utilized in the study have been classified into two distinct categories: stationary and nonstationary approaches.These approaches are applied to the selected seven atoll groups as part of the analysis.Additionally, within the stationary category, two subcategories have been defined.The first subcategory involved the creation of an SPDE mesh that encompasses the entire study area, including both islands and water bodies.The  second subcategory focused on the design of SPDE mesh precisely on the land surface, specifically targeting the land on reefs.While the second category employed a nonstationary approach, which primarily revolved around the development of barrier models.Section 3.1 provides a detailed discussion of the theoretical background and application techniques associated with the nonstationary barrier models.
A comprehensive understanding of individual mesh configurations for all seven atoll groups is reported in Appendix Section B. This includes the mesh structures used for both stationary models and the non-stationary model.Specific details of these mesh structures can be found under the heading ''Mesh of Atolls'' in Appendix Section B. However, for the purpose of demonstrating the applications of both stationary and non-stationary models in a clear and understandable manner, we have selected Seenu atoll as an illustrative example.By focusing on Seenu Atoll, we aim to provide a comprehensive explanation and facilitate a better understanding of the methodology and findings of our study.
Because of the highly distributed nature of the island structure in each atoll, a continuous spatial structure is initially chosen for modeling, and triangulation is performed on the entire study area.According to Verdoy (2019), the best fitting mesh should have enough vertices for effective prediction, but the number should be within a limit to have control over the computational time.Based on this principle, a series of meshes with varying numbers of vertices are constructed for each of the three modeling scenarios mentioned in Section 1. Finally from the battery of meshes, the best-fitting mesh is selected.In each case, the deviance information criterion (DIC), Watanabe-Akaike information criterion (WAIC), and computation time are used to assess the performance of the models and to select the best fitting model by balancing model accuracy against complexity (Spiegelhalter et al. 2002).
After comparing the performance of the models with different meshes, the best-performing models have been selected, and the meshes used in those models are referred to as the best-fitted meshes.These best-fitted meshes strike a balance between model accuracy and computational efficiency, ensuring that the models can provide accurate and reliable results in a reasonable amount of time.The left panel of Fig. 4 depicts the selected mesh for the entire Seenu atoll, which includes both the land on reefs (highlighted in blue) and the surrounding lagoons and water bodies.The grey representation within the figure illustrates the SPDE mesh, which consists of 1189 vertices.Additionally, the locations of the islands affected by the tsunami are marked in red, providing a visual reference for their positioning within the context of the SPDE mesh and the Seenu atoll.
From Fig. 4 (left) it is obvious that the SPDE mesh is generated for the entire geographic boundary of the atoll, including the water bodies like the lagoon and ocean surface.However, the records of the tsunami-affected areas are limited to discrete spatial sites that are located on land on reefs for individual islands.Thus, the traditional stationary method of SPDE triangulation using the boundary of the entire region is not realistic.As a result, we need to design triangulation only using boundaries of the land on reefs for individual islands in the atoll.Recent studies (Chaudhuri et al. 2022(Chaudhuri et al. , 2023) ) explore the application of explicit network triangulation where SPDE mesh is restricted to linear networks rather than the entire study area.Using a similar methodology, we have designed SPDE mesh precisely on land on reefs instead of defining the entire atoll boundary.As mentioned earlier, balancing between precision and computation time we have finetuned the number of vertices to identify the best-fitted mesh.The grey representations on the right panel of Fig. 4 depict the best fitting mesh having 427 vertices projected only on the land on reefs.Locations of tsunami-affected islands are marked in red.
According to Wood et al. (2008) and Bakka et al. (2019), the default stationary method of designing triangulation only on the land on reefs has a serious drawback in that the Neumann boundary condition is often unrealistic and severely impacts the results.A comprehensive discussion of the challenges associated with boundary conditions and other limitations has been discussed in Section 4.

Barrier model
The models designed using SPDE triangulations discussed in the previous subsection assume stationarity and isotropy that is, the autocorrelation between two locations depends solely on the Euclidean distance.But coastal areas and regions with physical barriers often exhibit complex spatial patterns due to the presence of features such as coastlines, islands, or channels.These patterns cannot be adequately represented by stationary models that assume spatial homogeneity.Thus, while modeling events on dispersed island structures stationarity is an unrealistic assumption (Bakka et al. 2018).Similar coastline problems are reported by Ramsay (2002), Wood et al. (2008), andScott-Hayward et al. (2014).Moreover, stopping the triangulation at the coastline imposes the Neumann boundary conditions, also leading to unrealistic models (Bakka et al. 2018).This issue is common while exploring coastline and complex island problems.Other examples of physical barriers include road networks, power lines, categorical health sectors, and areas with different land use.To deal with the coastline problem, several studies proposed solutions by computing the shortest distance in water (Wang and Ranalli 2007;Scott-Hayward et al. 2014;Miller and Wood 2014).Ramsay (2002)  The comprehensive literature review reveals a clear requirement for a generalized model that effectively retains the primary strengths of previous approaches while simultaneously addressing and overcoming their main weaknesses.The generalized model should exhibit robustness in relation to the process of selecting the boundary polygon.Moreover, it should maintain a computational cost comparable to the stationary alternatives.It should also incorporate a smoothness parameter and a range parameter, which collectively determine the distance of spatial similarity (Bakka et al. 2019).Additionally, it is essential to emphasize that the ease of use for researchers in practical applications should not significantly exceed that of the stationary models, thus maintaining usability parity.In this line, Bakka et al. (2019) proposed an approach to handle nonstationary and anisotropic spatial processes with emphasis to handle complex archipelago structures in which the coastline is used as a physical barrier (Bakka et al. 2019).In their proposal, Bakka et al. (2019) approximated them using a finite element method based on the SPDE method.A system of two SPDEs is presented in this case, one for the barrier area, and the other for the remaining area.The following system of stochastic differential equations has a solution that is specifically a nonstationary spatial effect, denoted by u(s).
where u(s) is the spatial effect, X b is the barrier area and X n is the remaining area and their disjoint union gives the whole study area X. Ranges for the barrier and remaining areas are represented by r and r b respectively.r u is the marginal standard deviation.r is equal to ð o ox ; o oy Þ and W(s) stands for white noise.
It is worth noting that the barrier model is based on viewing the Mate ´rn correlation as a collection of paths through a simultaneous autoregressive (SAR) model, rather than as a correlation function on the shortest distance between two points.The local dependencies are manipulated to cut off paths crossing the physical barriers.In the next step, the new SAR model is formulated to SPDE format to represent the Gaussian field, with a sparse precision matrix that is automatically positive definite (Bakka et al. 2019).
In the study by Bakka et al. (2019), the water body has been considered as normal terrain, and distinct coastlines and boundaries are used as physical barriers.In contrast, in the current study, we have defined boundary polygons of individual land on reefs for each island as our study area and the water body acts as the physical barrier.Besides, tsunami-hit locations are the sample events considered precisely on the land area of the islands.Figure 5 depicts the barrier object (left panel) and barrier object with triangulations (right panel) generated using inla.barrier.polygonfunction from the R-INLA package.However, it includes a modification where the triangulation is designed using a barrier model where the water bodies (ocean and lagoons) act as physical barriers.Within the left  5, the grey region signifies the physical barrier, while the white area represents the land on reefs where the analysis of spatial dependency is conducted.Meanwhile, the right panel showcases the triangulation, including the physical barrier in grey.In both panels, the points marked in red indicate the locations of areas affected by tsunamis, which are utilized as event locations in the model.Furthermore, this section serves as an opportunity to present Seenu atoll as an example case for our discussion.For additional details and mesh structures of barrier objects of other atoll groups, please refer to Appendix Section B under the heading ''Mesh of Atolls''.

Model fitting
Based on the discussions in the previous subsections we designed a set of hierarchical Bayesian models with Poisson likelihood and PC priors (see Eq. 2).The steps for modeling the application include the spatial effect created with the mesh using the spatial locations by Mate ´rn covariance and then implemented using individual mesh structures.We have also considered independent and identically distributed Gaussian random effect, represented as iid in the modeling process.In total, we have fitted 7 different models for each atoll group.Meshes used in both stationary and nonstationary models for individual atoll groups are reported in Appendix Section B. The seven models represent different combinations of mesh types, incorporating independent and identically distributed (iid) random effects in both stationary and nonstationary scenarios.In order to assess the significance of spatial effects, one of the models employed in the study intentionally omits any consideration of spatial effects.It is worth mentioning that, no covariates are used in designing the models for the current study.Details of each model are shown in Table 2.The dependent variable across all models is the number of individuals indirectly affected by tsunami.
As we have a battery of competing models, we compare them using the deviance information criterion (DIC) (Spiegelhalter et al. 2002), which is a Bayesian model comparison criterion, represented as where DðhÞ is the deviance evaluated at the posterior mean of the parameters, and p D denotes the effective number of parameters, which measures the complexity of the model (Spiegelhalter et al. 2002).An alternative is the Watanabe Akaike information criterion (WAIC) which follows a more strict Bayesian approach to construct a criterion (Watanabe 2010).pWAIC is similar to p D in the original DIC (Gelman et al. 2014).The lowest values of DIC and WAIC suggest the best-fitted model.

Results and discussion
As outlined in Section 3, our study employs a range of methodological approaches to analyze the spatial variability in the occurrence of natural disasters, specifically focusing on tsunamis in complex spatial regions having physical barriers.In this study, we have selected seven atoll groups located across the northernmost to the southernmost regions of Maldives, serving as a representative showcase.Each atoll group possesses unique characteristics, distinguished by factors such as shape, land area, and number and size of component islands, as well as the number of islands affected by tsunami.Additionally, the distributions of tsunami-affected islands within these atoll groups display unique patterns.In the two northern atoll groups (Shaviyani and Baa), the locations of affected islands are observed both along the coastline and within the enclosed boundaries of the reefs (refer to Figs. 9, 10).Conversely, in the central and southern atoll groups, the majority of tsunami-affected regions are concentrated along the periphery of the land on reefs (see Figs. 11,12,13,14,15).
In our analysis, we utilized two different scenarios when creating the SPDE triangulations: one for the entire region and another only focused on the land on reefs.Consequently, there exists a wide variation in the number of nodes (vertices) for each atoll group between these two cases.Table 4 in the Appendix shows a comparative overview of the node counts for each atoll group in both scenarios.It is observed that while considering the entire region, the number of nodes ranges from 1189 in Seenu atoll to 3001 in Gaffu Dhaal atoll group.Conversely, when considering only the land on reefs, the node counts are considerably smaller, ranging from 427 in Seenu to 1561 in Meemu atoll.Notably, the largest disparity between these two approaches is observed in the case of the Gaffu Dhaal atoll group, with a variation from 3001 to 1283 nodes.Using these meshes, we applied two models to analyze the spatial dependencies of natural hazards in each atoll group.
In addition to compare the performance of these two stationary models, we have implemented a nonstationary barrier model.The nonstationary model differentiates between component islands with ocean or lagoon areas in between, enabling the analysis of the atoll as a whole while considering the boundaries of land on reefs for individual islands.These characteristics directly influence the goodness-of-fit of the models, allowing for a comprehensive assessment of their effectiveness in capturing the spatial dynamics and patterns of the studied phenomena.By considering these characteristics, we gain valuable insights into the unique features and complexities of the spatial region, contributing to a better understanding of the performance and suitability of the models for analyzing the given context.Specifically, we conducted five stationary models and two nonstationary models.Table 2 provides a comprehensive list of the models considered, and in the subsequent sections, we present and discuss the results of the model fittings with a focus on specific aspects of the atolls.
Table 3 shows the DIC and WAIC values, which serve as widely utilized diagnostics for evaluating the quality of models in a Bayesian context (Rue et al. 2009).We report that the values in bold signify the best-fit models for each atoll group.Lower values indicate better model performance.Our observations indicate that, based on the DIC criteria, the nonstationary model incorporating the iid (Model 6) outperforms other models for nearly all atolls.However, it is important to note that both Model 6 and Model 7 exhibit better performance than Model 5, while Model 2 and Model 3 outperform Model 1 for each atoll.These findings suggest that the choice of mesh, specifically considering the entire region of atolls that includes physical barriers such as lagoons and the ocean, has a substantial impact on model performance.Consequently, the stationary model considering only land on reefs and the nonstationary barrier model, which accounts for the exclusion of physical barriers, outperform the other stationary model and provide a more realistic representation of the spatial characteristics.
Another interesting observation is related to the use of iid random effects.When comparing pairs of models with and without the use of iid for both stationary and nonstationary models (M1-M5, M2-M7, and M3-M6), it becomes evident that models incorporating iid consistently outperform their respective contrasting models.This suggests that the inclusion of independent and identically distributed random effects has a significant impact on model performance when combined with spatial effects.However, we do observe some disparities in certain atoll groups.For instance, the Shaviyani and Baa atoll groups exhibit exceptionally high DIC values for Model 2. These two atolls are characterized because the events occur both, on the coastline and also in the pieces of land inside the boundary reef (see Figs. 9,10).This distribution could explain why the model which applies SPDE, considering only the land on reef region, does not work properly in these two cases as the mesh used does not take into account the space between the component islands.In contrast, when applying the models to the Kaafu, Laamu, and Gaafu Dhaal atolls, lead to negligible differences between the stationary model only on land on reefs and nonstationary barrier model.This suggests that both types of models perform similarly in these atoll groups.It is possible that the structures of these atolls play an important role in this similarity.The tsunami-affected islands in these atolls are primarily concentrated on the boundary of the reefs, with only one instance in the case of Gaafu Dhal atoll where the event is located in the middle of the atoll (as depicted in Fig. 14).This particular distribution of events on the boundaries can explain the similarities observed between the two models used.In this case, the mesh structures do not encounter major challenges in accurately modeling the events, as they are primarily focused on the boundaries.Therefore, distinguishing between land and water (as done in the barrier model) may have less importance in these atolls, leading to comparable performance between the stationary and nonstationary models.
In the case of the Meemu atoll, all models show almost the same DIC values.This can be attributed to the characteristics of the atoll, which consists of a single land region and where all the tsunami events are located around the boundary reef (see Fig. 12).Consequently, in this scenario, the distinction between boundaries and water using a mesh may not be relevant, as there are only events within a single land region.As a result, the model's fit remains consistent irrespective of whether the applied mesh considers the atoll as a whole, treats the component islands individually, or employs the barrier model.Furthermore, when examining the atoll located further south (Seenu), noticeable distinctions arise when comparing the stationary models with the nonstationary model, particularly when iids are not included.These results are surprising because the structure of the territory of this atoll and the locations of the events are similar to that observed in the previously mentioned atoll groups.One notable distinction is that the Seenu atoll exhibits a relatively lower number of locations affected by tsunamis compared to the other atolls considered.Consequently, the number of events emerges as an important factor influencing the model fit in this particular context.
The comparison between stationary and nonstationary models at the atoll level provides valuable insights into the spatial variation of tsunamis in complex spatial regions like the islands of Maldives.The results suggest that a barrier model is an effective and realistic option, particularly in situations with regions having physical barriers.Notably, the findings also highlight the significance of the locations of Regarding the limitations and challenges, the most important one is related to boundary effects in the SPDE-INLA approximation.This approach creates artifact spatial dependencies on the boundary.In a standard mesh, as long as it is well constructed, the boundaries are in the outer limits of the spatial domain of interest and, therefore, those dependencies can be identified and eliminated.However, in a more complex mesh, such as barrier models, the boundaries lie within the spatial domain of interest.This fact makes it difficult, sometimes excessively, to identify and subsequently eliminate artifact spatial dependencies.Another important constraint lies in the choice of the range fraction, represented as r b r (refer to Eq. 5).This parameter introduces a degree of arbitrariness into the model results.The choice of an inappropriate range fraction can significantly impact the model outcomes and their interpretation.It is essential to select an appropriate range fraction that effectively captures the spatial dynamics of the barrier and aligns with the specific characteristics of the studied region.Another limitation arises when dealing with physically thin barriers.The barrier model may yield ambiguous results in such cases, as its ability to accurately capture the effects of these thin barriers is inherently limited.This introduces uncertainties in the model outcomes, and caution should be exercised when interpreting the results in such scenarios.Conversely, when modeling artificial thicker barriers, there is a risk of overlapping segments, creating fictitious spatial structures that do not correspond to the actual physical characteristics of the region.These overlaps can distort the model's output and compromise the reliability of the results.Therefore, when dealing with thin barriers that require manual adjustment, it becomes crucial to exercise careful consideration in selecting the appropriate barrier thickness parameter.Additionally, the investigation of tsunami propagation across multiple atolls at the island scale presents considerable difficulties, necessitating the use of high-resolution shape files.However, the availability of enhanced spatial resolution maps of the Maldives enables the development of detailed island-specific structures and the creation of more efficient nonstationary models.

Conclusions
The barrier model exhibits robustness in relation to the process of selecting the boundary polygon in a more realistic environment.Furthermore, it maintains a computational cost comparable to the stationary alternative.The model also incorporates a smoothness parameter and a range parameter, which collectively determine the spatial similarity distance.
This study exhibits certain limitations, such as the presence of boundary effects, the requisite for precise finetuning of the range fraction, and the optimal thickness of the barrier.Despite these difficulties, in this work, we have managed to identify them.However, a different approximation is needed in which the SPDE-INLA approximation does not cause these fictitious spatial dependencies.Moreover, these fictitious dependencies cause a very low predictive capacity of the model.Consequently, it is essential for future research to prioritize the exploration and development of an alternative approach that eliminates these artificial dependencies, thereby enhancing the model's predictive performance and enabling more reliable and robust spatial analysis.
However, this study demonstrates notable advantages.Primarily, the proposed nonstationary technique outperforms the traditional stationary spatial models to analyse the spatial dependencies of natural hazards in complex distributed land structures like islands and coastal areas.It entails equivalent computational expenses as the stationary spatial models, thereby streamlining their analysis.In general, the proposed model is easy to use and can deal with both sparse data and complex land structures having physical barriers.It provides a more accurate representation of the spatial dynamics of natural hazards that facilitates the identification of vulnerable areas, enabling targeted mitigation measures to be implemented.Additionally, understanding the spatial variability of natural hazards aids in post-effect disaster management.It allows for the assessment of the extent and severity of the impact on coastal communities, enabling efficient allocation of resources for response and recovery efforts.This knowledge of spatial dependency could be incorporated to the development of coastal communities that possess resilience and are better equipped to effectively mitigate and respond to future natural hazards.9,10,11,12,13,14,15 and Table 4.

Fig. 1
Fig. 1 Republic of Maldives geographical location and island structure

Fig. 2
Fig. 2 Left: Locations of tsunami affected islands for individual atolls of Maldives.Right: Study area (integrated 7 groups of atolls highlighted) where K m ð:Þ is the modified Bessel function of second order, and m [ 0 and j [ 0 are the smoothness and scaling parameters, respectively.INLA approach constructs a Mate ´rn SPDE model, with spatial range r and standard deviation parameter r.The parameterized model we follow is of the form ðj 2 À DÞ ða=2Þ ðsSðxÞÞ ¼ WðxÞ ð4Þ

Fig. 3
Fig. 3 Tsunami affected regions of Baa and Raa atolls.Left: Boundaries of atolls.Right: Boundaries of land on reefs for component islands proposed a methodology defining boundary conditions that use a smoothing penalty together with the Neumann boundary condition.Furthermore, Wood et al. (2008) and Sangalli et al. (2013) demonstrate alternative solutions based on the Dirichlet boundary condition where a known value or, function is used along the boundary.
; for s 2 X n ;

Fig. 4
Fig. 4 Left: SPDE triangulation with tsunami-affected regions for the entire atoll boundary.Right: Triangulation generated only for land on reefs

Table 1
Records of tsunami effects in the selected 12 Maldivian atolls

Table 2
Competing models with choice of SPDE mesh and iid

Table 3
DIC and WAIC values according to the seven models and for the considered atolls , as observed in the case of Shaviyani and Baa atoll groups, where it has a substantial impact on the model fit.It can be stated that, through the utilization of nonstationary barrier models, we can capture the intricate spatial variations and complexities present in coastal and island regions. events

Table 4
Number of vertices in SPDE mesh for seven atoll groups Atoll group SPDE mesh (entire region) SPDE mesh (only land)