A Modeling and Analysis Framework for Integrated Energy Systems Exposed to Climate Change-Induced NaTech Accidental Scenarios

This paper proposes a novel framework for the analysis of integrated energy systems (IESs) exposed to both stochastic failures and “shock” climate-induced failures, such as those characterizing NaTech accidental scenarios. With such a framework, standard centralized systems (CS), IES with distributed generation (IES-DG) and IES with bidirectional energy conversion (IES+P2G) enabled by power-to-gas (P2G) facilities can be analyzed. The framework embeds the model of each single production plant in an integrated power-flow model and then couples it with a stochastic failures model and a climate-induced failure model, which simulates the occurrence of extreme weather events (e.g., flooding) driven by climate change. To illustrate how to operationalize the analysis in practice, a case study of a realistic IES has been considered that comprises two combined cycle gas turbine plants (CCGT), a nuclear power plant (NPP), two wind farms (WF), a solar photovoltaicS (PV) field and a power-to-gas station (P2G). Results suggest that the IESs are resilient to climate-induced failures.


Introduction
Infrastructures providing water, energy, gas and transportation, etc., to citizens and industries are highly interconnected. Interconnection implies that localized failures can trigger cascade failures leading to service interruption [1]. The analysis of such complex infrastructures typically starts from the identification of independent failure mechanisms and causes and then focuses on dependent and common cause failures (CCF) [2].
External situations of potential dependent failures are related to natural technological (NaTech) events, i.e., natural phenomena that can damage technological infrastructures. The effects of NaTech events on system integrity are widely reported; the flooding that occurred in 1999 at the Blayais (FR) nuclear power plant (NPP), and in 2011 in the Dai-Ichi (JP) NPP and St. Lucie (USA) NPP, are two such examples. NaTech events bring additional stress on components and systems, with effects on their reliability and risk. With respect to climate change, in [3][4][5][6], to name a few, the effects of temperature on the reliability of electric lines, the safety systems of a NPP, the availability of hydropower systems and the cooling water of power systems are analyzed, respectively; in [7,8], the impact of wind and lightning hazards on oil and gas processing plants are evaluated. On the other hand, this work considers "shock" events specifically, such as floods and earthquakes, whose frequency of occurrence and severity are increasing due to climate change [9][10][11][12][13][14] and call for a specific framework of analysis.
In this work, we present a novel modeling framework in which both stochastic and "shock" climate-induced failures (i.e., NaTech events) are jointly considered for the analysis

The Modeling and Analysis Framework
The simulation framework, shown in Figure 1, is operationalized within a double-loop Monte Carlo simulation of n s scenarios, wherein stochastic and climate-induced failures (described in Sections 3.2 and 3.3) are sampled and injected into a power-flow model of the IES (see Section 3.1), where Γ different production plants are integrated. The inputs are the parameters of the system layout (CS, IES-DG, IES+P2G), the sea level projections, the corresponding flooding hazard curves in future (the year y * ) and the types of failures to be considered (stochastic, climate-induced or both).
The pseudocode of the simulation framework is as follows: (1) For i = 0 to n s : (2) i = i + 1 (2.1) For γ = 1 to Γ: 2.1.i Sample the occurrence times t j [yy] of the j-th stochastic failure and climate-induced event (t j > 50 years are neglected, as the longest plant useful life considered (i.e., that of the NPP) is equal to 50 years); 2.1.ii Sample the flooding level and define the fragility curve f γ , γ = 1, 2, . . . , Γ, in relation to the sampled flooding level.
(2.2) Sort t j in increasing order (i.e., a total of n e events are considered for each scenario). (2.3) For j = 0 to n e 2.1.iii If the j-th event is a stochastic failure event: The power-flow model described in Section 3.1 is run to calculate the energy import E i and the energy losses E l (see Equation (11)) during the 24 h, that are stored in E i (i) and E l (i). (2.5) j = j + 1 and return to 2.3.i.
(3) Record the number of outages n (i.e., energy import events (E i (t) > 0)) and the internal losses Int l ; (4) Return to 2. (5) Indicators evaluation: (5.1) SAIFI is computed as in Equation (1): 2) For the y * -th climate change projection, the AFP is calculated as the annual import probability that is given by the number of simulations that result in E i (t) > 0, divided by the total number of simulations: 3) The Zobel resilience metric is computed as: where T is the mean outage duration, i.e., the time between the start of the system failure event and its restoration. (5.4) For the y * -th climate change projection, the LEP is given by the ratio between the number of simulations resulting in Int l (t) > 5.7% (i.e., the Italian average value) during the specific year, and the total number of simulations: (5.5) CVaR is computed with respect to the distribution of energy import E i : where β is the confidence level.

The Case Study
We consider an IES composed of different production plants three possible system layouts, i.e., a standard CS, an IES implem

The Case Study
We consider an IES composed of different production plants that can be allocated in three possible system layouts, i.e., a standard CS, an IES implementing energy hub (IES-DG) and an IES adding bidirectional energy conversion (IES+P2G). In detail, the CS layout considers large-scale plants of either traditional or renewable energy sources, the IES-DG layout contains multiple medium-size prosumers, referred to as energy hubs (EH) [22], whereas the IES+P2G considers the addition of a storage, in the form of a P2G conversion plant, that converts the excess energy in the grid into natural gas and feeds it into the pipelines. For illustrating purposes, the generic IES considered consists of an NPP (for baseload generation), two CCGT plants (for both baseload and peak regulation), a solar PV field, two WFs, a compression station (to overcome losses in the gas pipes) and a P2G station for energy storage (that can be switched off for simulating the CS layout). The system considered is plotted in Figure 2, adapting it from previous works, such as [23,24]. In particular, we assume that the IES mimics the positioning of a number of realistic plants based in central Italy between Lazio and Campania regions ( Figure 3): for the NPP we assume the data of the Garigliano BWR nuclear reactor, CCGT are Napoli-levante and Teverola power plants, WF consists in fleets of Vesta V90 (2000 kW) turbines and solar PVs are fields of 1 kW PV panels with 35 • and 180 • of tilt and azimuth angles, respectively, summing up to 200 MW (Table 1), as a compromise of the results presented in [25]. In each node, energy can be either injected or absorbed into/from the grid: the six production plants (nodes 1 to 6) and the eight user nodes (nodes 7 to 14) are connected in a ring with nominal voltage of 220 kV, where a radial grid (blue in Figure 2) distributes gas and feeds the gas plants and gas customers. In Tables 2 and 3, the physical parameters (length, resistance and reactance) of the electric grid connecting the m-th and the n-th nodes corresponding to the production plants, and of the pipes connecting the ϑ-th and ϕ-th nodes of the pipeline network (length and frictional factor) are listed, respectively. Notice that a virtual additional node (15) is added to the grid (not shown in Figure 1) to account for the "import", when the IES cannot provide enough energy to customers.    Without loss of generality, the power demand of user nodes has been to a typical power demand profile during July 2018 [26], as shown in Fig  been taken as a reference demand because it challenges the IES due to the values that must be supplied while complying, at the same time, with the f nical constraints: minimum load for CCGT plants equal to 40% of the nomin maximum power for the CCGT equal to 110% of nominal power [28] and load only with constant power production at nominal power.  Without loss of generality, the power demand of user nodes has been set according to a typical power demand profile during July 2018 [26], as shown in Figure 4. This has been taken as a reference demand because it challenges the IES due to the high demand values that must be supplied while complying, at the same time, with the following technical constraints: minimum load for CCGT plants equal to 40% of the nominal power [27], maximum power for the CCGT equal to 110% of nominal power [28] and NPP for baseload only with constant power production at nominal power.  The characteristics and models of the EHs that are used in the IES-DG and IES+P2G layouts have been taken from [22]. The EH internal energy demand is satisfied with a small wind turbine and a combined heat and power station (each EH has a total generation capacity of 1200 kW to satisfy its demand and can exchange power with the grid). In our work, this means that EHs can dispatch power to the grid in case of excess of power production to satisfy normal users demand, and the overall power demand of Figure 4 must be discounted by the internal demand and the dispatched power, when the eight EHs, each one connected to one load node, are considered.

The Simulation of the IES Response
A power-flow model is used to calculate voltage V and phase at each bus of a power system, for a specified load, generator power and voltage condition [24]. This entails defining analytical models for each component of the power system, e.g., electric grid, gas pipeline and energy conversion system. In practice, with reference to our case study, the power flow allows computing Vm and in each m-th electric grid node, the active and reactive power and , respectively, generated/absorbed power in the m-th node and, indirectly, pressure in each -th gas distribution node and volume flow rate in each pipeline connecting the -th and the -th nodes of the gas network [29]. The power-flow model provides, at each iteration τ, the steady state solution of � ( ) (hereafter indicated as � ) and of F( � ), correspondingly: where � is the variable matrix composed of voltage phase angle vector ( ) at the m-th electric node, voltage magnitude vector (| ′|), CCGT power generated ( ′ ), P2G and gas compressor power demand ( ′ 2 , ′ ), pressure at the -th pipeline node ( ) and mass flow rate in the -th pipeline edge ( ); ( � ) is the system function matrix composed by active and reactive power (∆ ′, ∆ ′), pressure (∆ ), mass flow rate (∆ ) and The characteristics and models of the EHs that are used in the IES-DG and IES+P2G layouts have been taken from [22]. The EH internal energy demand is satisfied with a small wind turbine and a combined heat and power station (each EH has a total generation capacity of 1200 kW to satisfy its demand and can exchange power with the grid). In our work, this means that EHs can dispatch power to the grid in case of excess of power production to satisfy normal users demand, and the overall power demand of Figure 4 must be discounted by the internal demand and the dispatched power, when the eight EHs, each one connected to one load node, are considered.

The Simulation of the IES Response
A power-flow model is used to calculate voltage V and phase θ at each bus of a power system, for a specified load, generator power and voltage condition [24]. This entails defining analytical models for each component of the power system, e.g., electric grid, gas pipeline and energy conversion system. In practice, with reference to our case study, the power flow allows computing V m and θ m in each m-th electric grid node, the active and reactive power P m and Q m , respectively, generated/absorbed power in the m-th node and, indirectly, pressure p ϑ in each ϑ-th gas distribution node and volume flow rate Q ϑϕ in each pipeline connecting the ϑ-th and the ϕ-th nodes of the gas network [29]. The power-flow model provides, at each iteration τ, the steady state solution of X(τ) (hereafter indicated as X) and of F X , correspondingly: where X is the variable matrix composed of voltage phase angle vector (θ m ) at the m-th electric node, voltage magnitude vector ( V m ), CCGT power generated (P gen ), P2G and gas compressor power demand (P P2G , P compr ), pressure at the ϑ-th pipeline node (π ϑ ) and mass flow rate in the ϑϕ-th pipeline edge (S ϑϕ ); F X is the system function matrix composed by active and reactive power (∆P , ∆Q ), pressure (∆π), mass flow rate (∆S) and compressor power (∆P compr ).
To simulate the IES of Figure 2 working in nominal conditions (in all the considered layouts of CS, IES-DG and IES+P2G) when the power demand is as in Figure 4 [26], the power-flow model is run 24 times, each one considering the system stationary conditions at each hour of the day (τ = 1, 2, . . . , 24). From this, the following quantities can be calculated:

•
The daily energy exchanged by each m-th node: where P m (τ) is the active power of the m-th node at the τ-th hour; • The overall energy supply E s : • The overall energy demand E d : • The system energy losses (%) E l , due to distribution losses along cables: Under normal conditions (i.e., neither stochastic nor climate-induced failures are effecting the operation of the production plants and, therefore, the IES) and when the power demand is as in Figure 4, the energy supplied by each production plant is shown in Figure 5, where the different shades of colors correspond to the energy supplied by each production plant (CCGT, PV, WF and NPP, respectively, in blue, yellow green and red). In all cases, a difference between the demand curve (dashed line) and the supply can be noticed (and quantified with E l as in Equation (11)). In all cases where E i = E s − E d > 0, some energy is imported through the virtual node (15) and, in this work, for simplicity, but without loss of generality, the IES is considered failed, because not capable of fulfilling its function of guaranteeing the regional energy independence and security.
We also notice that: (1) For the IES-DG layout, the difference is smaller than for the CS layout, because the EHs produce (locally) part of the required demanded energy, that is actually discounted to the total demand; (2) In the IES+P2G layout, the positive balance of EHs in the central hours of the day and the PV power generation allows CCGT plants "to follow" the demand curve more gradually, independently from renewable energy production oscillation, thanks to the P2G. Indeed, the P2G, in case of high renewable generation, stores energy in the form of gas and avoids large-power demand oscillations. (3) E l (which increases with the power transmitted by the grid) is the largest for the IES+P2G layout (5.33%), whereas, thanks to local production, E l are the smallest for the IES-DG layout (3.51%) and the CS layout (4.28%) (vlues compared with Italian grid average E l of 5.7% in 2018 [26]). Sustainability 2022, 14, x FOR PEER REVIEW 9 of

The IES Reliability Model
The production plants (γ =1, 2, 3, 4, 5, 6 of Table 1) that compose the IES can fail d to both stochastic and climate-induced failures. Stochastic failures usually originate fro fatigue and strain of components [2], and their uncertain failure (or repair) times are, usual [15], modeled as exponentially distributed: where , ( , * | )is the probability that the γ-th production plant fails (recovers) at ti t (i.e., moves from the nominal state k to a failed state * (or viceversa)), and and are the failure and repair rates, respectively (see Table 4, where and are given each γ-th production plant) [30].

The IES Reliability Model
The production plants (γ = 1, 2, 3, 4, 5, 6 of Table 1) that compose the IES can fail due to both stochastic and climate-induced failures. Stochastic failures usually originate from fatigue and strain of components [2], and their uncertain failure (or repair) times are, as usual [15], modeled as exponentially distributed: where P γ,s (t, k * |k) is the probability that the γ-th production plant fails (recovers) at time t (i.e., moves from the nominal state k to a failed state k * (or viceversa)), and λ γ and r γ are the failure and repair rates, respectively (see Table 4, where λ γ and r γ are given for each γ-th production plant) [30].
Wind farm 0.031 4 [29] Climate-induced failures are, instead, originated from a natural event that might affect more than one plant at the same time; therefore, climate-induced failures cannot be considered statistically independent (as we assume for the stochastic failures), calling for a different modeling approach. Failures are, indeed, here modeled as "shock" events affecting all plants at the same time that may or may not fail under the shock received, depending on their fragility. For quantitative evaluation, the probability P γ,c (t|δ) of failure of the γ-th plant due to the δ-th natural event is calculated as: where P(t, γ|δ) is the probability that the γ-th plant fails at time t due to the occurrence of the δ-th event (i.e., the fragility of the γ-th plant to the δ-th event), and P(δ) is the probability of δ-th event occurrence (e.g., a specific δ-th flooding level occurrence). P(t, γ|δ) is calculated using fragility curves that can be obtained by fitting failure databases. In this work, we assume the following fragility curves; for each γ-th power production plant: • Solar PV panels (γ = 1) are assumed to fail with certainty (P(t, γ = 1|δ) = 1) when the flooding level exceeds 1 m and to not fail (P(t, γ = 1|δ) = 0) for lower flooding levels, as plotted in Figure 7 (we conservatively assume that 1 m is the height at which electric equipment is mounted on the PV metal structure and this is the equipment that would be damaged by flooding). • NPP (γ = 2) is assumed to fail with certainty (P(t, γ = 2|δ) = 1) when flooding level exceeds 5.7 m, and to not fail (P(t, γ = 2|δ) = 0) for lower flooding levels, as plotted in Figure 8 [31]. • CCGT power plants (γ = 3, 4) are assumed to fail with P(t, γ = 3, 4|δ) as shown in in Figure 6; different failure probability curves are given for six damage states (negligible, very low, low, medium, relevant and severe) of the concrete walls of CCGT power plants [32]. In this work, the fragility related with the medium damage state is considered (bold line in Figure 8). • WF (γ = 5, 6) are assumed to fail with certainty (P(t,γ = 5, 6|δ) = 1) when flooding level exceeds 2 m, and to not fail (P(t,γ = 5, 6|δ) = 0) for lower flooding levels, as plotted in Figure 9 (we conservatively assume that 2 m is the maximum flooding level withstood by the transformer connecting the plant to the grid, which is the equipment that would be damaged by the flooding).   • WF (γ = 5, 6) are assumed to fail with certainty (P(t,γ = 5, 6|δ) = 1) when floo level exceeds 2 m, and to not fail (P(t,γ = 5, 6|δ) = 0) for lower flooding leve plotted in Figure 9 (we conservatively assume that 2 m is the maximum flooding withstood by the transformer connecting the plant to the grid, which is the eq ment that would be damaged by the flooding).    Flooding occurrence probability P(δ) can be obtained from the outcomes of TSUMAPS-NEAM (http://www.tsumaps-neam.eu/, accessed on June 2019), an international collaborative project aimed at developing tsunami hazard maps for coastal areas in the North-East Atlantic Mediterranean (NEAM) region. As an example, in Figure 10      It is worth pointing out that, as discussed in [11], flooding can originate from dif ent threats, such as heavy rainfall, storm surges or tsunamis that, as a chain of eve often overlap each other. For simplicity, but without loss of generality, in this work assume that floods are modeled independently from each other (i.e., each δ-th even completely resolved before the following one is initiated). It is worth pointing out that, as discussed in [11], flooding can originate from different threats, such as heavy rainfall, storm surges or tsunamis that, as a chain of events, often overlap each other. For simplicity, but without loss of generality, in this work we assume that floods are modeled independently from each other (i.e., each δ-th event is completely resolved before the following one is initiated).

NaTech Events under Climate Change
The climate-induced failures discussed in Section 3.2 deserve particular care when modeled, due to the variability across space (latitude and longitude) and time (shortor long-time projections). This means working with climate data specific to the IES site (rather than on general worldwide projections) and, also, focusing on specific initiating events (i.e., seismic activity), rather than a multitude of generic sets of natural events. In this work, we tailor the analysis on a specific site (latitude 40 • 50 N, and longitude 14 • 15 E, where the IES described in Section 3.1 is located) and consider different time projections (in the y * (year) 2040, 2070 and 2100) for the flooding hazard curves. As flooddriven failures are strongly related to climate change, which induces stronger and more frequent events [9], we model the flooding severity to increase with the sea level rise that is ultimately dependent on the temperature rise [33]. In practice, the 8.  Figure 11), starting from the baseline curve of Equation (14) proposed by [14] as the probability of flood level in 2020:

NaTech Events under Climate Change
The climate-induced failures discussed in Section 3.2 deserve particular care w modeled, due to the variability across space (latitude and longitude) and time (sho long-time projections). This means working with climate data specific to the IES sit ther than on general worldwide projections) and, also, focusing on specific initi events (i.e., seismic activity), rather than a multitude of generic sets of natural even this work, we tailor the analysis on a specific site (latitude 40° 50′ N, and longitude 1 E, where the IES described in Section 3.1 is located) and consider different time projec (in the * (year) 2040, 2070 and 2100) for the flooding hazard curves. As flood-d failures are strongly related to climate change, which induces stronger and more freq events [9], we model the flooding severity to increase with the sea level rise that is mately dependent on the temperature rise [33]. In practice, the 8.  Figure 11), starting the baseline curve of Equation (14) proposed by [14] as the probability of flood lev 2020: 1.05√2 � It is worth noticing that between 2020 and 2100, the flooding occurrence proba almost doubles (at fixed severity) and that the most frequent events consistently inc in severity, endangering the IES integrity, as long as the climate changes.

Results
The framework presented in Section 2 has been applied for the analysis of th described in Section 3. The following analysis considers the IES+P2G layout (i.e., the complete layout among CS, IES-DG and IES+P2G) subject to stochastic and climat duced failures under a climate change projection over the horizon to 2100 (i.e., the It is worth noticing that between 2020 and 2100, the flooding occurrence probability almost doubles (at fixed severity) and that the most frequent events consistently increase in severity, endangering the IES integrity, as long as the climate changes.

Results
The framework presented in Section 2 has been applied for the analysis of the IES described in Section 3. The following analysis considers the IES+P2G layout (i.e., the most complete layout among CS, IES-DG and IES+P2G) subject to stochastic and climateinduced failures under a climate change projection over the horizon to 2100 (i.e., the most conservative and threatening scenario), to analyze which failure mechanism (stochastic or climate-induced) plays the major role.
We consider n s = 10,000 alternative stochastic scenarios, which are run on MATLAB with a Dell XPS: Intel Core i7-8550U CPU @1.8 Ghz and 8 GB RAM.
In Figure 12, the SAIFI boxes are shown for comparing the results when stochastic, climate-induced failures and (their) combined effects are considered. It is worth mentioning that, as already stated, state-of-the art framework of analysis entails accounting only for stochastic failures. With this proposed approach, instead, it can be quantified the joint effect of the large number of import events due to the stochastic failures (due to the large failure rate of the CCGT), and the lower number of climate-induced failures (due to low-frequency and high-severity events). a matter of fact, all the CVaR values, when considering climate-induced failures, are larger than those of stochastic failures. AFP and LEP, shown in Figures 13 and 14, respectively, are close to 1 when stochastic (only) and (both) stochastic and climate-induced failures are considered, due to the large share of power produced by the CCGT; any stochastic failure of any single plant may lead to an energy import, whereas climate-induced failures, which more likely affect small power production plants, can be easily covered by the remaining operating plants.
The effect of climate-induced failures on the performance of the IES+P2G layout is noticeable and quantifiable, with the proposed modeling and analysis framework: for example, the resilience of the system due to climate-induced shock events (blue line in Figure 15) demonstrates less resourcefulness for recovering from natural hazards events with respect to only stochastic failures (red line). This is due to the fact that, in general terms, low severity (high frequency) climatic events affect PV and WF (characterized by low production capacity and long repair times), whereas the frequent stochastic failures that affect the CCGT plant are more promptly recovered, as they are likely to occur and maintenance readiness is expected to recover the systems.  The R and CVaR (for β = 90%, 95% and 99% confidence levels) are provided in Table 5, and again for comparing the results of stochastic, climate-induced failures and (their) combined effects; climate-induced failures are shown to have the largest impact on E i . As a matter of fact, all the CVaR values, when considering climate-induced failures, are larger than those of stochastic failures. AFP and LEP, shown in Figures 13 and 14, respectively, are close to 1 when stochastic (only) and (both) stochastic and climate-induced failures are considered, due to the large share of power produced by the CCGT; any stochastic failure of any single plant may lead to an energy import, whereas climate-induced failures, which more likely affect small power production plants, can be easily covered by the remaining operating plants.      The effect of climate-induced failures on the performance of the IES+P2G layout is noticeable and quantifiable, with the proposed modeling and analysis framework: for example, the resilience of the system due to climate-induced shock events (blue line in Figure 15) demonstrates less resourcefulness for recovering from natural hazards events with respect to only stochastic failures (red line). This is due to the fact that, in general terms, low severity (high frequency) climatic events affect PV and WF (characterized by low production capacity and long repair times), whereas the frequent stochastic failures that affect the CCGT plant are more promptly recovered, as they are likely to occur and maintenance readiness is expected to recover the systems.  In conclusion, the analysis of the case study here provided shows how the proposed framework can be operationalized to jointly consider and model stochastic and climateinduced failures. The outcomes obtained warn against the fact that neglecting the contribution of the latter could expose an IES to additional risks: the identification and quantification of these by the proposed framework can also inform the proper planning and designing IES, and support decision making and policy-making when adaptation to climate change is to be duly considered.

Conclusions
In this work, we have presented a modeling framework for the analysis of IES that considers both stochastic failures and climate-induced events. The framework, which is operationalized with a double-loop Montecarlo simulation, has been applied to a typical IES+P2G layout of power production exposed to flooding scenarios, to evaluate a pool of representative indicators, such as SAIFI, AFP, R index, LEP and CVaR. Results suggest that power systems withstand better climate change when energy production systems are integrated, as in the IES+P2G layout. This is because, as expected, the failures on certain plants can be mitigated by the capability offered by other plants and EHs connected to the network. It is worth mentioning that, if the IES+P2G suffers (from the resilience point-ofview) the P2G electric demand that makes the system more vulnerable to climate-induced failure events, it is the preferred solution when the infrastructure is suffering from a gas shortage. As a last remark, it is important to point out that the presented results, due to the large uncertainty that affects the climate change models, need a rigorous uncertainty quantification and analysis to allow for decision making with confidence. Future work In conclusion, the analysis of the case study here provided shows how the proposed framework can be operationalized to jointly consider and model stochastic and climateinduced failures. The outcomes obtained warn against the fact that neglecting the contribution of the latter could expose an IES to additional risks: the identification and quantification of these by the proposed framework can also inform the proper planning and designing IES, and support decision making and policy-making when adaptation to climate change is to be duly considered.

Conclusions
In this work, we have presented a modeling framework for the analysis of IES that considers both stochastic failures and climate-induced events. The framework, which is operationalized with a double-loop Montecarlo simulation, has been applied to a typical IES+P2G layout of power production exposed to flooding scenarios, to evaluate a pool of representative indicators, such as SAIFI, AFP, R index, LEP and CVaR. Results suggest that power systems withstand better climate change when energy production systems are integrated, as in the IES+P2G layout. This is because, as expected, the failures on certain plants can be mitigated by the capability offered by other plants and EHs connected to the network. It is worth mentioning that, if the IES+P2G suffers (from the resilience point-ofview) the P2G electric demand that makes the system more vulnerable to climate-induced failure events, it is the preferred solution when the infrastructure is suffering from a gas shortage. As a last remark, it is important to point out that the presented results, due to the large uncertainty that affects the climate change models, need a rigorous uncertainty quantification and analysis to allow for decision making with confidence. Future work will, therefore, focus on the development of uncertainty quantification approaches that can capitalize the Monte Carlo simulations results that the framework presented here can provide.

Conflicts of Interest:
The authors declare no conflict of interest. Probability of failure of the γ-th plant due to the δ-th climate induced natural event P γ,s (t, k * |k)

List of Symbols
Probability that the γ-th production plant fails (recover) at time t (i.e., moves from the nominal state k to a failed state k * (or viceversa)) due to stochastic events P m Active power P(δ) Probability of δ-th event occurrence P( f l. level, y * ) Probability of flood level in y * P(t, γ|δ) Probability that the γ-th plant fails at time t due to the occurrence of the δ-th event P compr Gas compressor power demand P gen CCGT power generated P P2G P2G power demand p ϑ Pressure in the ϑ-th gas distribution node Q ϑϕ Volume flow rate in each pipeline connecting the ϑ-th and the ϕ-th nodes of the gas network