Georgian Geographical Journal 14 Integrating GIS and Kernel Density Estimation for Multi-Hazard Risk Assessment of Potential Onshore Pipeline Failures Aslan E. Babakhanov1 , Zaur T. Imrani1* 1 Department of Tourism and Recreational Geography, Ministry of Science and Education, Institute of Geography named after academician H.A. Aliyev, Baku, Azerbaijan * Corresponding author: zaur_imrani@mail.ru Introduction Strengthening the durability of Azerbaijan’s oil and gas infrastructure against natural disasters is a critical priority for researchers and industry stakeholders. Pipelines that transport hydrocarbons, while integral to national and transnational energy supply, are especially vulnerable due to the volatile nature of their contents. A failure can lead to irreparable environmental damage, human casualties, and costly socioeconomic disruptions (Lerche et al., 2014). The growing reliance on pipelines for interregional energy exchange, coupled with the apparent increase in natural disaster frequency, has underscored the need for comprehensive risk assessments (Krausmann et al, 2011). These assessments must extend beyond purely technical criteria to include geopolitical and environmental considerations, reflecting the international implications of pipeline failures and the cross-border nature of catastrophic events. Climate change adds further complexity by intensifying extreme weather events and altering hazard patterns, placing critical infrastructure under unprecedented levels of stress (Zio, 2016). As natural disasters become more frequent and severe, safeguarding pipelines requires a flexible and forward- looking framework capable of integrating diverse data sources, advanced modeling, and predictive analytics. This forward-thinking perspective enables risk assessments that anticipate future scenarios rather than merely responding to current conditions. Georgian Geographical Journal, 2025, 5(3) 14-29 © The Author(s) 2025 This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). DOI: https://journals.4science.ge/index.php/GGJ Abstract This paper introduces a newly developed risk-assessment framework that combines Kernel Density Estimation (KDE) with Geographic Information System (GIS) technology to analyze multiple hazards, including earthquakes, floods, landslides, mud volcanoes, and soil erosion, for the Baku-Tbilisi- Ceyhan pipeline in Azerbaijan. By integrating hazard-specific parameters into a unified risk matrix, each hazard’s contribution is weighted, refined, and aggregated to produce a spatially explicit, combined risk map. KDE smooths hazard intensities and reveals overlaps among different risk factors. The resulting high-resolution maps enable more targeted prevention and response measures, guiding planners and stakeholders toward effective pipeline protection strategies. Although the model can demand computational power, it remains scalable and flexible, allowing for adaptation to additional hazards or expanded geographical areas. Furthermore, the methodology underscores the importance of cross-validation in setting KDE bandwidth and in calibrating hazard weights to ensure reliable outputs. Preliminary testing indicates that this integrated model improves the clarity of risk data, highlights areas needing immediate attention, and supports resilience planning across the pipeline corridor. This work can be applied more broadly to critical infrastructure projects in regions where multiple hazards coincide, thereby aiding decision-making processes for disaster risk reduction and sustainable development. Future research will focus on refining statistical models for inter-hazard correlations and incorporating machine learning for predictive analytics. The framework stands as a tool to maintain pipeline integrity in the face of evolving environmental threats. Keywords: GIS, Kernel Density Estimation, Risk Matrix, Pipeline Integrity, Multi-Hazard Risk. Citation: Babakhanov, A.E.; Imrani, Z.T. Integrating GIS and Kernel Density Estimation for Multi-Hazard Risk Assessment of Potential Onshore Pipeline Failures. Georgian Geographical Journal 2025, 5(3), 14-29. https://doi.org/10.52340/ggj.2025.05.03.02 mailto:zaur_imrani@mail.ru https://orcid.org/0000-0002-8790-7945 https://orcid.org/0000-0002-0606-3753 Babakhanov & Imrani 2025 5(3) 15 Developing robust methodologies to quantify and mitigate these risks is thus an imperative, particularly in regions like Azerbaijan where seismic activity, floods, landslides, and mud volcanoes pose substantial threats (Othman et al., 2023). With pipelines traversing extensive distances across varied terrains, even localized disruptions can have extensive consequences. Consequently, the scientific community must invest in innovative tools and interdisciplinary strategies that address physical vulnerabilities, regulatory complexities, and cross-border interdependencies. Adopting a multi-hazard approach allows for a holistic understanding of pipeline resilience, thereby guiding more effective risk management, informed policy-making, and international collaboration. Notably, specific natural disasters recognized in Azerbaijan can directly impact onshore pipelines, highlighting the urgency for targeted, data-driven solutions (Amirova‐Mammadova, 2018; Bagirov et al., 2019): 1. Earthquakes: Azerbaijan experiences substantial seismic activity due to its location in an earthquake-prone zone (Alizadeh et al., 2017). Sudden ground movements or shifts can displace or fracture pipelines. Any rupture in a pipeline carrying hydrocarbons may cause leaks and, in severe instances, lead to explosions or fires. 2. Floods: Flash floods or prolonged heavy rainfall can result in rapid water accumulation and significant soil erosion1. High-velocity water flows may directly damage exposed pipeline segments, while erosion can undermine the structural support of buried pipelines, increasing the risk of bending or rupturing. 3. Landslides: In mountainous or hilly terrains, landslides can occur, particularly following intense precipitation or seismic events2. The downward movement of rock and soil can bury pipelines or exert forces beyond their design limits, leading to material failure. 4. Mud Volcanoes: Azerbaijan is well-known for its numerous mud volcanoes (Kadirov et al., 2005). These geological formations can erupt violently, expelling hot mud and gases capable of overheating, encasing, or fracturing adjacent pipelines, ultimately causing leaks or more severe accidents (Panahi, 2005). 5. Soil Erosion: Continuous erosion, often accelerated by flooding or heavy rainfall, can gradually expose pipelines designed to remain buried (Othman et al., 2023). Without adequate coverage, pipelines become vulnerable to external impacts and stress that can lead to cracks or breaks. 6. Corrosion (Environmental Degradation): While not considered a conventional natural hazard, corrosion is a significant environmental threat to pipelines. Humid or saline conditions accelerate metal deterioration (Cheng, 2015), and if not detected, corrosion can thin pipeline walls, resulting in leaks or complete structural failure. 7. Extreme Weather Conditions: Severe weather events—such as powerful storms, lightning strikes, or extreme temperature fluctuations—can damage pipelines directly or trigger secondary hazards like landslides or floods (Katopodis et al., 2019). While pipelines are built to withstand local climate norms, unusually harsh conditions can exceed their tolerance thresholds. 8. Wildfires: High temperatures from wildfires can damage pipelines and their protective infrastructure (Novacheck et al., 2021). Moreover, the burning of nearby vegetation may contribute to soil erosion, further exposing pipelines to potential mechanical stresses. 9. Ground Subsidence: This phenomenon, caused by natural geological activity or human practices (e.g., mining, groundwater extraction), results in the gradual sinking of the ground surface (Oruji et al., 2022). Uneven subsidence can place bending stress on pipelines, increasing the likelihood of ruptures. Typically, the practice is to represent the impact and risk values for each hazard individually on separate maps. A risk matrix serves as a visual aid, facilitating the evaluation of the cumulative risk level associated with different occurrences, including natural disasters. This matrix is instrumental in pinpointing the probability and potential impact of a given event, thereby allocating it a precise risk rating. However, the variety of natural disasters differs by location, resulting in a diverse set of risk maps. This diversity can introduce complexities in managing and interpreting the array of maps (Samany et al., 2022). 1 Pipeline safety: Potential for damage to pipeline facilities caused by flooding, river scour, and river channel migration (2019). https://www.federalregister.gov/documents/2019/04/11/2019-07132/pipeline-safety-potential-for-damage-to-pipeline- facilities-caused-by-flooding-river-scour-and-river 2 Guidelines for Constructing Natural Gas and Liquid Hydrocarbon Pipelines Through Areas Prone to Landslide and Subsidence Hazards (2009). https://rosap.ntl.bts.gov/view/dot/34640 https://www.federalregister.gov/documents/2019/04/11/2019-07132/pipeline-safety-potential-for-damage-to-pipeline-facilities-caused-by-flooding-river-scour-and-river https://www.federalregister.gov/documents/2019/04/11/2019-07132/pipeline-safety-potential-for-damage-to-pipeline-facilities-caused-by-flooding-river-scour-and-river https://rosap.ntl.bts.gov/view/dot/34640 Babakhanov & Imrani 2025 5(3) 16 Using a combined approach helps streamline this complexity by integrating multiple risk factors into a unified visual representation, simplifying analysis and decision-making (Falcone et al., 2022). To derive a combined risk value encompassing all-natural disasters, one could consider aggregating the distinct risk values (Fotios et al., 2022) for each disaster. However, this method presupposes that all these incidents occur independently, a condition that might not always hold true. The overall risk assessment can be significantly swayed by the interrelations and mutual influences existing between various natural disasters. Recognizing these interdependencies is essential in accurately determining the combined risk value. To maintain the integrity of all assessed risk values, the individual combined risk scores are connected to beforehand computed and proportionately weighted risk values for each type of natural disaster. A critical initial phase involves quantifying the risk (Gemma et al., 2022) tied to each possible natural disaster. This could involve analyzing historical data on the frequency of natural disasters like earthquakes, floods, wildfires, and so on, within a targeted area to gauge their probability. The core objective is to create a robust, multi-dimensional risk assessment model that not only identifies potential risks but also quantifies them in a meaningful and actionable manner (Rasouli & Imrani, 2023). To achieve this, we will employ a risk rating matrix grounded in rigorous mathematical formulations, enabling a detailed and nuanced understanding of risk levels. This matrix will integrate various risk factors listed above to generate a composite risk score for pipeline systems in the context of natural disasters. By combining the spatial analysis capabilities of GIS with the mathematical rigor of risk score calculations, this research endeavors to offer a novel and practical tool for policymakers, engineers, and disaster management professionals. The ultimate goal is to enhance decision-making processes, facilitate proactive risk mitigation strategies, and contribute to the resilience and safety of pipeline infrastructures in the face of increasingly unpredictable natural events. Methods and Materials In this research, we utilize adjusted risk matrices that correspond to each relevant hazard type - encompassing both classical natural disasters (e.g., earthquakes, floods) and environment-based threats (e.g., corrosion). Each hazard is weighted according to its relative impact on pipeline integrity. These weighted factors are then refined using Kernel Density Estimation (KDE) techniques, converting the collected risk scores into a 2D array (Gramacki, 2017). The data assigned to each cell within these matrices can be stored in any relational database (Nasser, 2018) or in straightforward file formats such as comma-separated values (CSV) or JSON. This flexibility enables seamless integration with popular GIS tools like ESRI ArcGIS or QGIS, where the 2D risk arrays can be visualized, analyzed, and used to inform decisions on pipeline protection and maintenance. 2.1 Area of observation This research area (Fig.1) situated within Azerbaijan, known for its susceptibility to various natural disasters such as earthquakes, landslides, floods, and mud volcanoes. This region is identified as the "valley of mud volcanoes" Moreover, satellite imagery revealed the presence of strong wind patterns and noticeable geological faults in the area. 2.2 Risk Matrix When analyzing the effects of natural disasters on pipeline systems, employing a risk matrix is an essential step for illustrating the degree of risk. This study utilizes a standard 5×5 risk matrix (Figure 2), which classifies threats according to their likelihood (ranging from rare to nearly certain) and potential consequences (spanning from negligible to disastrous) (Blokdyk, 2018). By adopting a five- tier system, each hazard can be systematically plotted on the grid, where higher probability and severity correspond to higher risk values. Although summing the individual risk values for all-natural disasters can yield an overall, aggregated risk score, this approach assumes complete independence among events, which is not always accurate. Correlations and interdependencies between different types of natural disasters can significantly alter the true combined risk. Accordingly, such factors must be considered when determining an aggregated risk value. Babakhanov & Imrani 2025 5(3) 17 Figure 1. Area of observation In spatial terms, each risk factor is mapped onto polygonal zones varying from as small as 50×50 meters to areas of several square kilometers (Han et al., 2010). These zones are categorized by risk level according to the matrix. Table 1 outlines a distribution framework for assigning integer or categorical values (e.g., “Low,” “Moderate,” “High”) to each zone based on local hazard intensity. This classification is a cornerstone for subsequent calculations in a GIS environment, ensuring that risk values can be systematically integrated with geospatial data. As a result, risk maps visualize both individual and aggregated hazards, thereby aiding in pipeline safety assessments and informing targeted mitigation strategies. Figure 2. Risk matrix of 5x5 elements Table 1. Distribution of risk levels according to ND factor Natural Disaster Risk Level (1-5, 5 is highest) Description Color Landslides 4 Soil shifts in hilly or mountainous regions that have the potential to disrupt or harm pipelines. Yellow Earthquakes 5 Seismic activity or ruptures leading to significant damages Red Flooding 3 Lead to erosion or the accumulation of sediment Yellow Babakhanov & Imrani 2025 5(3) 18 Corrosion from soil chemistry 3 Specific soil environments may extend the corrosion, resulting in the deterioration of pipeline materials Yellow 2.3 Determining the Risk Values To address numerous hazard events, the contributions from each are consolidated into a unified risk score. We begin by computing the separate risk scores (𝑅𝑖) for every hazard, following the formula: 𝑅𝑖 = 𝐿𝑖 × 𝐶𝑖 where 𝐿𝑖 represents the likelihood and 𝐶𝑖 denotes the impact of the 𝑖𝑡ℎ hazard. After determining the individual risk scores for each hazard (𝑅1, 𝑅2, …, 𝑅𝑛), the average risk score (𝑅𝑎𝑣𝑔) is computed by the subsequent method: 𝑅𝑎𝑣𝑔 = 1 𝑛 ∑ 𝑅𝑖 𝑛 𝑖=1 In this scenario, n represents the total count of hazards or events under evaluation. Averaging these values is one of several methods to consolidate the risk data, and it is particularly effective for the purposes of mapping. It's important to note that the risk scores will be visually differentiated using a color-coding system: scores ranging from 1 to 6 will be marked in Green, indicating low risk; scores between 7 and 12 will be highlighted in Yellow, signifying medium risk; and scores from 13 to 25 will be denoted in Red, representing high risk. Additionally, each specific value within these ranges could be represented by varying shades of the respective color, providing a gradient effect for more nuanced visualization Fig. 2. Indeed, the method of calculating an aggregate risk score, as previously described, presumes that each hazard is independent and of equal importance, an assumption that might not hold true in more complex scenarios. To address this, our approach will integrate a weighting system, which allows for the differentiation in the significance of various hazards or the interdependencies between events. In this system, each hazard is assigned a weight (𝒘𝒊) reflecting its relative importance or impact. These weights are typically values ranging from 0 to 1, with the sum of all weights equaling 1. This ensures that the overall significance of all hazards is proportionately distributed. Following the assignment of weights, the next step involves calculating the weighted risk score (𝑅𝑤𝑖) for each hazard. This is achieved by: 𝑅𝑤𝑖 = 𝑤𝑖 × 𝑅𝑖 To obtain the comprehensive risk score (𝑅𝑡𝑜𝑡𝑎𝑙), the process involves summing all the weighted risk scores for each hazard: 𝑅𝑡𝑜𝑡𝑎𝑙 = ∑ 𝑅wi 𝑛 𝑖=1 This total risk score is then categorized using the designated color-coding system based on its value: 𝐶𝑜𝑙𝑜𝑟(𝑅𝑡𝑜𝑡𝑎𝑙) = { 𝐺𝑟𝑒𝑒𝑛 𝑖𝑓 1 ≤ 𝑅𝑡𝑜𝑡𝑎𝑙 ≤ 6 𝑌𝑒𝑙𝑙𝑜𝑤 𝑖𝑓 7 ≤ 𝑅𝑡𝑜𝑡𝑎𝑙 ≤ 12 𝑅𝑒𝑑 𝑖𝑓 13 ≤ 𝑅𝑡𝑜𝑡𝑎𝑙 ≤ 25 This weighted approach enhances realism by recognizing that certain hazards (e.g., high-magnitude earthquakes) may be far more consequential than others (e.g., minor flooding). However, neither the averaging nor the weighted-sum methods inherently capture potential correlations among hazards— such as how an earthquake might trigger a landslide or how floods may exacerbate soil erosion. Capturing these interdependencies typically requires more advanced statistical or probabilistic models (e.g., Copulas, Bayesian Networks, or multi-variate correlation matrices3). Despite these limitations, weighted risk calculations represent a practical, GIS-friendly means of consolidating disparate hazards into a single metric. By overlaying the final risk values on geospatial layers, decision-makers can pinpoint pipeline segments requiring greater inspection, maintenance, or protective measures. This approach also aligns well with advanced techniques, such as Kernel Density Estimation, which can further smooth risk values spatially and highlight high-risk clusters within the study region, supporting more robust and proactive risk management. 2. 4 Evaluation of 2D risk scores using Gauss kernel method The 2D risk score for natural disasters can be evaluated using geospatial parameters, such as the magnitude of previous events (for earthquakes, landslides, etc.), proximity to fault lines, soil type, slope, 3 Multivariate correlations, https://numericalexpert.com/tutorials/statistics/multivarcorr.php Babakhanov & Imrani 2025 5(3) 19 and rainfall data (for floods and soil erosion). The method also assumes that each type of disaster has its own unique contributing factors. A Gauss kernel function4 is often used in spatial data analysis because it is smooth, symmetric, and its value decreases with distance together with Kernel Density Estimation (KDE). KDE is a non-parametric technique for estimating the probability density function of a random variable. When applied to risk levels, it allows you to smooth the risk measurements and find areas of high and low risk density. The equation (Equation 1.2) provides a smoothed estimate of the risk level distribution. From this, you could estimate probabilities of specific risk levels, find modes of risk (most common risk levels), or perform other analyses based on selection of kernel function for the KDE (Equation 1.1). The Gaussian kernel is a common choice due to its smoothness and nice mathematical properties. K(u) is the Gaussian kernel function, which decreases with the square of the distance from the center of the kernel (u), and f(x) is the kernel density estimate, which is an average of the kernel functions centered at each data point (Xi), with the bandwidth h which controls the amount of smoothing: a large h leads to more smoothing and a smaller h leads to less smoothing. To apply this to measured risk levels, the data points (Xi) would be the measured risk levels, and the estimated density f(x) would give an estimate of the probability density of risk levels at each point in the risk range. Equation 1.1: The equation for a Gauss kernel 𝐾(𝑢) = 1 √2𝜋 𝑒− 1 2 𝑢2 Equation 1.2: KDE equation 𝑓(𝑥) = 1 𝑛ℎ ∑ 𝐾 ( 𝑥 − 𝑋𝑖 ℎ ) 𝑛 𝑖=1 After collecting the data on historical ND, we can then calculate the risk scores per disaster category using minimal set of crucial input parameters: 1. Earthquakes: Variables:  Magnitude (M) of the past earthquakes.  Distance (d) to the epicenter of past earthquakes. The risk score for earthquakes can be calculated as follows: 𝑅(𝑥, 𝑦) = ∑ 𝑀𝑖 𝑖 ⋅ 𝑒 − 𝑑𝑖 2 2σ2 2. Landslides: Variables:  Slope of the terrain (S).  Soil moisture content (M).  Distance (d) to previous landslides. The risk score for landslides can be calculated as follows: 𝑅(𝑥, 𝑦) = ∑ 𝑆𝑖 𝑖 ⋅ 𝑀𝑖 ⋅ 𝑒 − 𝑑𝑖 2 2σ2 3. Flooding: Variables:  Precipitation (P) in the area.  Distance (d) to water bodies such as rivers or lakes. The risk score for flooding can be calculated as follows: 𝑅(𝑥, 𝑦) = ∑ 𝑃𝑖 𝑖 ⋅ 𝑒 − 𝑑𝑖 2 2σ2 In these equations, σ is the standard deviation, controlling the spread of the Gaussian kernel, (x, y) are the coordinates for a specific grid cell in the GIS maps, R(x, y) is the calculated risk score at point (x, 4 The Kernel Cookbook, https://www.cs.toronto.edu/~duvenaud/cookbook/ Babakhanov & Imrani 2025 5(3) 20 y), and the summation i is over all the relevant events or factors within a certain radius around point (x, y). 2.5 Applying the minimal safe distances The minimal safe distance to the object (here the minimal distance to the linear part of pipeline) is a starting point to identify the way we calculate the offset for risk values. Minimal distance is a geometrical distance, a straight line. It varies from disaster factors, so let’s integrate it (Dmin) into our equations: 1. Earthquakes: Variables:  Magnitude (M) of the past earthquakes.  Distance (d) to the epicenter of past earthquakes.  Minimal safe distance to the epicenter of earthquakes (Dmin). The risk score for earthquakes can be calculated as follows: E(𝑥, 𝑦) = ∑ 𝑀𝑖 𝑖 ⋅ 𝑒 − (𝑑𝑖−𝐷𝑚𝑖𝑛)2 2σ2 ⋅ 𝐻(𝑑𝑖 − 𝐷𝑚𝑖𝑛) Here, H (di - Dmin) is the Heaviside step function. It is equal to 0 for di < Dmin (indicating no risk inside the safe distance), and 1 for di >= Dmin. 2. Landslides: Variables:  Slope of the terrain (S).  Soil moisture content (M).  Distance (d) to previous landslides.  Minimal safe distance to previous landslides (Dmin). The risk score for landslides can be calculated as follows: 𝐿(𝑥, 𝑦) = ∑ 𝑆𝑖 𝑖 ⋅ 𝑀𝑖 ⋅ 𝑒 − (𝑑𝑖−𝐷𝑚𝑖𝑛)2 2σ2 ⋅ 𝐻(𝑑𝑖 − 𝐷𝑚𝑖𝑛) 3. Flooding: Variables:  Precipitation (P) in the area.  Distance (d) to water bodies such as rivers or lakes.  Minimal safe distance to water bodies (Dmin). The risk score for flooding can be calculated as follows: 𝐹(𝑥, 𝑦) = ∑ 𝑃𝑖 𝑖 ⋅ 𝑒 − (𝑑𝑖−𝐷𝑚𝑖𝑛)2 2σ2 ⋅ 𝐻(𝑑𝑖 − 𝐷𝑚𝑖𝑛) These equations will yield risk scores that are higher for locations closer to the dangerous object (beyond Dmin), and zero for locations within the safe distance. The more sophisticated approach is considering a gradual decrease of risk within the safe distance, instead of a sudden drop to zero (Section 2.5). The results of evaluated risk scores (RS) of a single contributing factor for each type of disaster (e.g., a single earthquake, a single previous landslide, a single water body) using above equations can be found on Table 2: Table 2. Risk scores per factor, 500 meters, 5 steps Distance (m) Earthquake RS (M=5, σ=500) Landslide RS (S=2, M=0.5, σ=500) Flood RS (P=100, σ=500) 200 0.891 0.577 19.2 700 0.706 0.353 14.12 1200 0.367 0.184 7.34 1700 0.135 0.067 2.7 2200 0.033 0.017 0.66 Here, M is the earthquake magnitude, S is the terrain slope, and M is the soil moisture content for landslides (not to be confused with M for earthquakes). P is the precipitation level. σ is the spread of the Gauss kernel, which is the same for all three types of disasters. Babakhanov & Imrani 2025 5(3) 21 The risk scores were calculated using the updated equations from the “Section 2.4”, with Dmin equal to 200 meters. Note that for the distance of 200 meters, the risk scores are 0 because of the safe distance factor. The risk scores for distances beyond Dmin were calculated by substituting the given values into the equations. The calculations assume a specific value for each parameter, and the actual risk scores may vary greatly depending on these values. The total risk function R(x, y) at location (x, y) would be the sum of these three functions: 𝑅(𝑥, 𝑦) = 𝐸(𝑥, 𝑦) + 𝐿(𝑥, 𝑦) + 𝐹(𝑥, 𝑦) The function R(x, y) will give us a single risk score for each location (x, y) that we can visualize on a 2D map. We can calculate this function for each cell in a 5×5 grid, and then color each cell based on its risk score to create a 2D risk map. This is used as a model for single cell area and actual risk maps would probably use a more sophisticated model and a much larger grid. 2.6 Splitting and scaling the grid In multi-hazard GIS-based risk analyses, subdividing the area of interest from larger cells (Δx,  Δy) into smaller ones can greatly enhance the precision and interpretability of results. Each cell 𝐺𝑖, 𝑗 represents a discrete spatial unit for modeling hazard intensity, vulnerability, and other relevant parameters. Adopting a finer cell size (Δ𝑥small ,  Δ𝑦small) has several advantages: Increased Spatial Resolution. Smaller grid cells capture more localized variations in elevation, land cover, infrastructure density, and other critical attributes. In mathematical terms, the number of cells 𝑁 in a given region of area 𝐴 scales approximately as 𝑁 ∝ 𝐴 Δ𝑥×Δ𝑦 . Hence, decreasing cell size yields a larger 𝑁, facilitating higher- resolution risk modeling. Improved Risk Assessment. By refining each cell’s spatial dimensions, the risk matrix (or any comparable index) can more accurately depict local hazards. For instance, in a flood risk scenario, slight changes in elevation or land cover within a 2×2 km cell may significantly affect water flow and flood extent. A coarser 10×10 km cell (Δ𝑥large,  Δ𝑦large) would mask these localized variations, potentially underestimating or overestimating the actual risk. Localized Analysis. Disasters such as landslides or urban flash floods often affect small areas with disproportionately high severity. A finer grid enables targeted analysis of these hotspots, capturing the nuances of topography and land use that can influence the severity of impacts. This capability is crucial for large-scale disaster management strategies where local conditions can markedly change risk levels. Detailed Mitigation Planning. High-resolution grids reveal specific locations most prone to damage, supporting more efficient mitigation measures. Authorities can use these refined data layers to plan structural defenses (e.g., levees or retaining walls), allocate evacuation routes, and prioritize response resources, thereby reducing both immediate and long-term risk. Enhanced Accuracy of Predictive Models. Many computational models in hydrology, seismology, and other hazard-related fields use grid-based inputs. Smaller cell sizes often improve the fidelity of simulations, though at the cost of greater computational demand. By employing a grid of fine resolution, models can represent spatial heterogeneity more accurately, thereby yielding more reliable predictions of hazard behavior. However, adopting smaller cells increases the computational burden, as the total number of cells N and the associated data complexity grow. Additionally, finer spatial resolution demands correspondingly detailed data inputs (e.g., higher-accuracy digital elevation models, land use surveys). In practical applications, the grid cell size is often a compromise between accuracy and available computational or data resources. Babakhanov & Imrani 2025 5(3) 22 Figure 3. The area of interest with 3x2 cells, each 10x10 km For example, Figure 3 may represent a coarser grid where each cell is approximately 10×10 km, capturing large-scale trends but providing only a coarse view of localized risk. Conversely, Figure 4 might illustrate a finer grid (2×2 km cells), capturing nuanced variations yet increasing both data density and processing time. Calibrating grid size becomes essential for balancing computational feasibility with the desired level of spatial detail. In many multi-hazard pipeline risk assessments, a moderate cell size is often chosen initially, followed by targeted refinement in critical zones where hazards overlap or infrastructure vulnerability is high. This approach ensures that risk maps reflect realistic spatial gradients without overwhelming computational or data-storage capacities. Figure 4. The area of interest with 17x12 cells, each 2x2 km 2.7 Smoothing risk scores by categories When plotting discrete or granular risk values directly onto a map—especially when each value is subdivided into finer sub-scores—visual representations can appear patchy or irregular. To obtain a smoother, more continuous distribution of risk values across a two-dimensional space, Kernel Density Estimation (KDE) offers a robust solution. Under this approach, probability density functions are approximated for each designated risk category (e.g., Green, Yellow, Red), focusing attention on smaller, localized zones that require further scrutiny. In this study, risk values lying outside the specified range of a given category are effectively excluded (or treated as zero) during computation. For instance, when constructing KDE surfaces for Green- category data, only points whose risk scores lie within the “Green” interval are included. Similarly, Babakhanov & Imrani 2025 5(3) 23 Yellow and Red categories are handled in separate density estimations, thus preventing extraneous data points from skewing results. Mathematically, for the c-th category, where 𝑛𝑐 data points (𝑥𝑖 (𝑐) , 𝑦𝑖 (𝑐) ) satisfy that category’s risk range, the two-dimensional Gaussian KDE can be written as: Equation 2: 2D risk score using Gauss kernel 𝑓𝑐(𝑥, 𝑦) = 1 𝑛𝑐ℎ2 ∑ 1 2𝜋 𝑛𝑐 𝑖=1 𝑒 − 1 2 [( 𝑥 − 𝑥𝑖 (𝑐) ℎ ) 2 + ( 𝑦 − 𝑦𝑖 (𝑐) ℎ ) 2 ] Here, categorized risk levels per impact, where (𝑛𝑐) is the number of data points in the (i)th category, (h) is the bandwidth parameter, (e) is the base of the natural logarithm and 𝑥𝑖 (𝑐) and 𝑦𝑖 (𝑐) are the nth pair of 2D data points. This formulation assumes X and Y (the spatial coordinates) are independent variables in the Gaussian kernel, allowing for straightforward, radially symmetric smoothing. The choice of the bandwidth h is pivotal: smaller h values yield highly localized densities at the expense of potential overfitting, whereas larger h values give smoother distributions that may obscure important local variations. Bandwidth selection often leverages cross-validation or rule-of-thumb heuristics (Silverman, 2018). By performing this KDE procedure separately for each risk category, the method generates multiple “layered” density surfaces that can then be visualized independently or composited. The resulting smoothed risk surfaces help highlight areas of greatest concern within each category and avoid the abrupt color transitions that simple, cell-by-cell risk mapping might produce. Consequently, decision- makers and analysts gain a clearer view of spatial risk concentrations, which supports more informed planning and more targeted mitigation strategies. 2.8 Data Processing To streamline the risk analysis within a specified Area of Interest (AIO), a systematic methodology was devised, as summarized in Table 3. The determination of the AIO’s minimum size - which was initially set at 1 km² - drew on a preliminary review of notable natural disaster events within Azerbaijan, ensuring that smaller-scale events would be captured. However, the AIO dimension is not fixed and may be recalibrated for different regions or urban areas where hazard intensity or infrastructure density warrants finer resolution. Adapting AIO dimensions typically involves additional analytical steps. First, detailed hazard datasets (e.g., seismic zonation maps, flood extents, or landslide susceptibility layers) are reviewed to confirm the spatial extent and granularity of the events. Second, the potential impact radius of each hazard is considered - particularly if one hazard (e.g., flooding) tends to spread over a wider area compared to another (e.g., localized ground subsidence). Finally, socioeconomic or administrative factors, such as population centers or critical infrastructure corridors, may also inform the selection of a larger or smaller AIO. Within each AIO, risk levels are determined by referencing previously established matrices for natural hazards (earthquakes, floods, landslides, mud volcanoes, etc.), environment-based hazards (e.g., corrosion), and their associated weights. The resulting risk scores are then aggregated spatially, allowing analysts to pinpoint critical hotspots. By adhering to the structured workflow shown in Table 3, one ensures uniformity across different AIOs, thereby facilitating consistent comparisons - even when hazard profiles differ. Through this stepwise approach, researchers and decision-makers can maintain a clear audit trail of how risk scores are derived, updated, and mapped onto each AIO. This enables rapid re-evaluation if new hazard data emerge - such as revised seismic models or recent flood records -and fosters an adaptable framework suitable for a variety of planning or operational needs. Table 3. Processing Steps for a Single AIO Step Description 1 Determine the AIO: The AIO should be centered around the pipeline to ensure accurate and proportionate GIS visualization of the event area and the pipeline's linear section 2 Identify (ND) and Adjust AIO Size: Identify NDs that could impact the AIO. The AIO's dimensions may be modified if an ND has a wider impact area. Concurrently, determine the risk parameters for each event and their relative weights concerning the pipeline. 3 Calculate Risk Matrices per Event: Develop risk matrices for each identified event, considering the specific risk parameters and weights. Babakhanov & Imrani 2025 5(3) 24 4 Combine Risk Scores: Combine the individual risk scores from each event into a comprehensive risk score matrix. This methodical approach ensures a thorough and systematic assessment of risks within a given AIO, allowing for a nuanced understanding of potential impacts on pipeline infrastructure. 2.9 Aggregated Risks While the aggregated risk score provides a high-level view of the overall risk, it is still tied to the individual risk scores, and these can be used to understand the specific hazards in more detail. It's beneficial to consider both the aggregated and individual risk scores when planning interventions and communicating about risk. In terms of visualization of risk details in an aggregated form, the individual risk being stored as an entity within the relational database. Many GIS frameworks providing the automation facilities (here ArcMap, QGIS and etc.) via scripting languages like Python, JavaScript. So, by adding the referential output per aggregated risk value (here individual colored cell) we can get the detailed information on certain cell of interest in a form of table 4 or other ways, like popup tooltips. Table 4: Aggregated risk score for each location Location Earth-quake Risk Land-slide Risk Flood Risk Mudflow Risk Aggregated Risk A 5.6 4.2 7.3 6.4 23.5 B 7.1 2.3 8.5 5.6 23.5 C 3.4 6.7 6.2 4.5 20.8 D 6.9 4.8 8.1 5.3 25.1 E 4.3 3.5 7.8 6.7 22.3 This table clearly shows both the individual and aggregated risk scores. By comparing the aggregated scores, we can see that location D has the highest overall risk, even though it may not have the highest risk for any individual disaster category. This is a good example of how aggregated risk scores can provide a different perspective than individual risk scores. 2.10 KDE and risk colors The KDE method can be used to generate the smoothed, continuous distribution of the aggregated risk scores. In this context, it would provide a way to understand the distribution of aggregated risk scores across the area of interest, identifying areas of high or low risk density. This could help to inform more granular risk management strategies or interventions, by focusing efforts where the density of high-risk locations is greatest. Figure 5. Histogram of Gaussian aggregated risk score The Fig. 5 is a smoothed histogram (i.e., a probability density function) of the aggregated risk scores using a Gaussian kernel. The resulting plot shows the distribution of risk scores, with the y-axis representing the estimated density of each risk score. Babakhanov & Imrani 2025 5(3) 25 This can help us to visualize how the risk scores are distributed and could potentially inform how we define the ranges for the 'Green', 'Yellow', and 'Red' risk rankings. Then, the distribution of colours on the histogram be colorized based on colour ranges per risk score. On the example of given aggregated values from the table “Table 4”, the risk ranking fall into ‘Red’ coloured risk score and be smoothed as a gradient of ‘Red’ colour. Following the calculation of the probability density function using the Gaussian kernel density estimation (KDE) technique, we visualize the distribution by creating a color-coded histogram plot (Fig. 6). Those gradients then be used to map smoothed aggregated risk values on final GIS maps. This results in a spatial map where each location is coloured based on its normalized risk density, as estimated by the KDE. Darker colours indicate regions with higher risk densities. In practical applications, this spatial risk map can be overlaid with other geographical features such as roads, buildings, and natural landmarks to provide a comprehensive view of the risk landscape. This can aid in identifying high-risk areas and prioritizing interventions. Figure 6. Gradient histogram of Gaussian aggregated risk score 2.11 Assessment challenges In the process of assessing risks to pipeline infrastructure due to natural disasters, several challenges arise, particularly in accurately characterizing the impact of these events. A key aspect of this challenge is the need for a comprehensive definition of risk, which encompasses a clear understanding of the probability and dimensions of potential impacts. The challenges, listed in order of increasing complexity, include: 1. Identification and Evaluation of Risk Parameters: This involves pinpointing specific risk factors and assessing their potential impact. An example is evaluating the effects of earthquakes on pipelines, including subsequent NATECH events (Krausman et al., 2016). 2. Determining Appropriate Risk Weights for Each Parameter: Assigning the correct weight to each risk parameter is crucial for accurate risk assessment. 3. Choosing the Correct Granularity for the Area of Interest: This is essential to avoid oversimplification and ensure detailed risk analysis. 2.12 General steps for risk visualization To visualize the spatial distribution of aggregated risk scores using the colorized histogram, we can employ a GIS framework. This allows us to map the risk scores onto a 2D map representation of the study area. a. Define the study area: Determine the extent and boundaries of the study area where the risk scores will be visualized. b. Obtain the risk scores: Collect or calculate the aggregated risk scores for each location within the study area. c. Create a grid: Divide the study area into a regular grid of cells. The size of each cell depends on the desired resolution for the visualization. Babakhanov & Imrani 2025 5(3) 26 d. Interpolate the risk scores: Assign each cell in the grid a risk score based on interpolation techniques such as inverse distance weighting or kriging. This step ensures that risk scores are assigned to locations that fall between the actual data points. e. Colorize the grid: Assign a color to each cell based on its risk score. Use a color scale or colormap that represents the risk levels of interest. For example, a colormap with shades of red can be used, where lighter shades represent lower risk and darker shades represent higher risk. f.Generate the risk map: Plot the colored grid onto the GIS map, overlaying it on relevant base layers such as roads, topography, or satellite imagery. This provides spatial context and aids in the interpretation of the risk distribution. g. Include a legend: Create a legend that clearly explains the color scale and risk levels associated with each color. This helps viewers interpret the risk map accurately. h. Optional: Add additional features: To enhance the visualization, you can incorporate other geographical features such as administrative boundaries, water bodies, or infrastructure that may be relevant to the risk assessment. By following these steps, the colorized histogram-based visualization of aggregated risk scores can be extended to a GIS environment, allowing for a comprehensive and visually appealing representation of spatial risk distribution. Results The outcome of this study includes the successful implementation of GIS tools to create detailed visualizations of aggregated risk. A significant enhancement in our visualization technique was the incorporation of buffered zones along the pipeline route, which provided a more distinct delineation of risk areas (Petersen, 2020). The series of figures, from Fig. 7 to Fig. 10, illustrate the progression of refined combined risk scores (CRS) derived from various natural hazards within a relatively small region. The smallest unit of analysis, the AIO, covered approximately 2 km², while the entire region under assessment spanned about 100 km². The risk matrices developed were informed by a combination of data on recorded earthquakes, identified geological faults, and mud volcanoes. The tools utilized for these visualizations were ESRI ArcMap 10.8 for mapping, output files generated (CSV files) with Python 3.8. Discussions In this study, we have navigated through the multifaceted challenges of assessing the risk to pipeline infrastructure from natural disasters. Our approach has been to meticulously identify and evaluate risk parameters, determine appropriate risk weights, and select the optimal granularity for the area of interest. Implications of Findings: Our findings underscore the complexity inherent in assessing natural disaster risks to pipelines. The use of a three-dimensional risk assessment plot, as demonstrated in Figure 7 for earthquake impacts, provides a nuanced visualization of risk. This method allows for a more detailed understanding of how different risk parameters interact and affect the overall risk profile. Figure 7-8. Coarse, CRS (left), KDE based, refined CRS (right) Babakhanov & Imrani 2025 5(3) 27 The color-coded risk categories further enhance the interpretability of these risks, making the data more accessible for decision-makers. The application of this methodology to other natural disasters reveals its versatility. However, it also highlights the increasing complexity when more parameters are involved. Each additional parameter requires careful consideration and integration into the overall model, underscoring the need for a robust and flexible risk assessment framework. Challenges and Limitations: One of the primary challenges encountered in this study is the sensitivity of the KDE and color mapping to the choice of parameters, especially the bandwidth h. The selection of h is critical and must be tailored to the specific characteristics of the data and the objectives of the risk assessment. This sensitivity points to a broader challenge in risk assessment: the balance between accuracy and generalizability. Furthermore, our approach assumes the independence of risk factors, which may not always hold true in real-world scenarios. Future Directions: Looking ahead, there are several avenues for further research. First, exploring methods to incorporate the interdependencies between different risk factors would likely yield a more comprehensive and realistic risk assessment model. This could involve the development of more sophisticated statistical models or machine learning algorithms capable of handling complex, interrelated datasets. Second, the adaptation of our model to a wider range of natural disasters would be beneficial. Each type of disaster brings its unique set of challenges and risk factors, and a versatile, adaptable model is essential for broad applicability. Lastly, engaging with experts in geology, meteorology, and disaster management could provide valuable insights and data, enriching the model's accuracy and relevance. Conclusion This study has illuminated several key challenges in determining risk weights, which are crucial for accurate risk assessment in the context of ND and their impact on infrastructure. These challenges are summarized below: 1. Complexity of Risk Interactions: ND often triggers or coincides with other events, complicating the assessment of combined risk impacts. Example: An earthquake increasing the likelihood of a landslide. 2. Data Scarcity: The lack of sufficient data for rare or unprecedented risks introduces significant uncertainty in estimating probabilities and impacts. 3. Risk of Oversimplification: Simplifying risk assessments can lead to underestimating the complexity and interactions of risks, potentially skewing the actual risk levels. 4. Model Dependence: The accuracy of risk weights is directly linked to the efficacy of the risk assessment model. Inaccurate models can result in misleading risk weights. 5. Subjectivity in Risk Weighting: Risk assessments can be influenced by individual or organizational risk tolerance, which varies widely. Balanced risk weighting often requires collective input from professionals or subject matter experts. 6. Evolving Risk Landscape: Risks change over time due to environmental, technological, regulatory, and other factors, necessitating regular updates to risk weights. 7. Perception Bias: Decision-makers' biases in risk perception can lead to overestimating or underestimating certain risks, influenced by factors like visibility, recency, or familiarity. Figure 9-10. KDE based, smaller scale CRS (left), KDE based, fine-tuned CRS (right) Babakhanov & Imrani 2025 5(3) 28 8. Quantifying Qualitative Risks: Some risks are qualitative and difficult to quantify, such as reputational risks or the impact of regulatory changes, posing challenges in numerical risk weighting. Advantages The utilization of aggregated risk scores in assessing various types of risks offers several advantages, as outlined below: 1. Simplicity in Visualization and Interpretation: Aggregated risk scores streamline the process of understanding risks by condensing multiple risk factors into a single, comprehensive score. This simplification aids in easier interpretation and visualization. 2. Enhanced Comparability: These scores enable straightforward comparisons of overall risk across different locations or regions, which is invaluable for decision-makers in resource allocation and prioritization. 3. Holistic Risk Understanding: Aggregated scores offer a broad overview of risk, capturing the cumulative effect of individual risks, especially when they interact or compound. 4. Effective Communication: Communicating a single aggregated risk score is often more straightforward and accessible to stakeholders and the public than explaining multiple individual risk components. Furthermore, the aggregated risk score can be tailored to align with the risk tolerance or appetite of the concerned entity. This adaptability allows for a more relevant and targeted risk assessment. Using a risk matrix for aggregation effectively consolidates risks from various natural disasters into a single value. This comprehensive view is crucial for effective risk management, disaster preparedness, and resource allocation. However, it's essential to periodically update the risk matrix to reflect changes in the risk landscape due to factors like climate change, urban development, and population growth. It's also important to recognize that the aggregated risk score is an indicative tool meant to guide, but not solely dictate, decision-making processes. Factors such as societal values, financial constraints, and political considerations should also be factored into risk management strategies. Despite the challenges, our study successfully visualized risk in an aggregated form, using a single map with finely tuned risk scores. Regular review and updating of risk weights, expert involvement, and careful consideration of risk interactions are key components of this approach. One notable finding is the model's applicability to other risk assessments using the same matrix representation and weight calculations. The flexibility of the model is further enhanced by the potential integration of the Kernel Density Estimation (KDE) method for fine-tuning outputs and the implementation of algorithms like the Moore neighborhood (Ilachinski et al., 2001) for relative AIO region risk score determination. Competing interests The authors declare that they have no competing interests. Authors’ contribution All authors provided critical feedback and helped shape the research, analysis and manuscript. ORCID iD Aslan E. Babakhanov https://orcid.org/0000-0002-8790-7945 Zaur T. Imrani https://orcid.org/0000-0002-0606-3753 Reference Alizadeh, A.A., Guliyev, I.S., Kadirov, F.A., Eppelbaum L.V. (2017). Geosciences of Azerbaijan - Volume II: Economic Geology and Applied Geophysics. Amirova‐Mammadova, S. (2018). Pipeline Politics and Natural Gas Supply from Azerbaijan to Europe: Challenges and Perspectives. Springer. Bagirov, E., Nadirov, R., & Lerche, I. (1996). Earthquakes, mud volcano eruptions, and fracture formation hazards in the South Caspian Basin: Statistical inferences from the historical record. Energy Exploration & Exploitation, 14(6), 585-606. Blokdyk, G. (2018). Risk Matrix a Complete Guide. 5starcooks. https://orcid.org/0000-0002-8790-7945 https://orcid.org/0000-0002-0606-3753 https://orcid.org/0000-0002-0606-3753 Babakhanov & Imrani 2025 5(3) 29 Cheng, Y. F. (2015). Pipeline corrosion. Journal Title, 161-162. https://doi.org/10.1179/1478422X15Z.000000000357 Gemma C., Carmine G., John McCloskey (2022). Modelling and quantifying tomorrow's risks from natural hazards. Science of The Total Environment, Volume 817. https://doi.org/10.1016/j.scitotenv.2021.152552. Gramacki, A. (2017). Nonparametric Kernel Density Estimation and Its Computational Aspects. Springer. Falcone D., G. Di Bona, Forcina (2022). A new method for risk assessment in industrial processes. IFAC- PapersOnLine, 55(19) 1-6. https://doi.org/10.1016/j.ifacol.2022.09.175 Fotios P., Daniele A., Vassilios A., Mohamed Zied Babai, Devon K. Barrow. (2022). Forecasting: theory and practice, International Journal of Forecasting, 38(3), 705-871. https://doi.org/10.1016/j.ijforecast.2021.11.001 Han, Z., & Weng, W. (2010). An integrated quantitative risk analysis method for natural gas pipeline network. Journal of Loss Prevention in the Process Industries, 23(3), 428–436. https://doi.org/10.1016/j.jlp.2010.02.003 Ilachinski, A., & Zane. (2001). Cellular Automata – A Discrete Universe. Kybernetes, 32(4). https://doi.org/10.1108/k.2003.06732dae.007 Kadirov F.A, Lerche I, Guliyev I.S, Kadyrov A.G, Feyzullayev A.A, Mukhtarov A.Sh. (2005). Deep Structure Model and Dynamics of Mud Volcanoes, Southwest Absheron Peninsula (Azerbaijan). Energy Exploration & Exploitation. 23(5):307-332. doi:10.1260/014459805775992717 Katopodis, T., & Sfetsos, A. (2019). A review of climate change impacts to oil sector critical services and suggested recommendations for industry uptake. Infrastructures, 4(4), 74. https://doi.org/10.3390/infrastructures4040074 Krausmann, E., Cruz, A. M., & Salzano, E. (2016). Natech Risk Assessment and Management: Reducing the Risk of Natural-Hazard Impact on Hazardous Installations. Elsevier. Krausmann, E., Renni, E., Campedel, M., & Cozzani, V. (2011). Industrial accidents triggered by earthquakes, floods and lightning: lessons learned from a database analysis. Natural Hazards, 59(1), 285–300. https://doi.org/10.1007/s11069-011-9754-3 Lerche, I., & Bagirov, E. (2014). Impact of Natural Hazards on Oil and Gas Extraction. Springer. Nasser, H. (2014). Learning ArcGIS Geodatabases. Packt Publishing Ltd. Novacheck, J., Sharp, J., Schwarz, M., Donohoo-Vallett, P., Tzavelis, Z., Buster, G., & Rossol, M. (2021). The evolving role of extreme weather events in the U.S. power system with high levels of variable renewable energy. Golden, CO: National Renewable Energy Laboratory. Oruji, S., Ketabdar, M., & Ketabdar, M. (2022). Evaluation of land subsidence hazard on steel natural gas pipelines in California. Upstream Oil and Gas Technology. Othman, A., El-Saoud, W. A., Habeebullah, T., Shaaban, F., & Abotalib, A. Z. (2023). Risk assessment of flash flood and soil erosion impacts on electrical infrastructures in overcrowded mountainous urban areas under climate change. Reliability Engineering & System Safety. https://doi.org/10.1016/j.ress.2023.109302 Panahi, B.M. (2005). Mud Volcanism, Geodynamics and Seismicity of Azerbaijan and the Caspian Sea Region. In: Martinelli, G., Panahi, B. (eds) Mud Volcanoes, Geodynamics and Seismicity. NATO Science Series, vol 51. Springer, Dordrecht. https://doi.org/10.1007/1-4020-3204-8_7 Petersen, K. (2020). Visualizing Risk: Drawing Together and Pushing Apart with Sociotechnical Practices. Rasouli, A., & Imrani, Z. (2023). Object-Based Analysis for Risk and Disasters Assessment: Tutoring on Geo- Environmental Contexts of Azerbaijan. Baku: CN Poliqraf LLC. 430. https://doi.org/10.5281/zenodo.14885894 Samany, N.N., Liu, H., Aghataher, R. et al. (2022). Ten GIS-Based Solutions for Managing and Controlling COVID-19 Pandemic Outbreak. SN COMPUT. SCI. 3, 269. https://doi.org/10.1007/s42979-022-01150-9 Silverman, B.W. (2018). Density Estimation for Statistics and Data Analysis. Routledge. Zio, E. (2016). Critical Infrastructures Vulnerability and Risk Analysis. Eur J Secur Res 1, 97–114. https://doi.org/10.1007/s41125-016-0004-2 https://doi.org/10.1108/k.2003.06732dae.007