Reliability Evaluation for Integrated Electricity-Gas Systems Considering Hydrogen

Regarded as the cleanest and a versatile energy carrier, green hydrogen generated with renewable energy is receiving increased attention in the transition to a carbon-neutral society. On the other hand, the integration of hydrogen tightens the coupling between power systems and natural gas systems, which highlights the importance of reliability evaluation of the integrated electricity-gas systems (IEGSs) with hydrogen. This paper proposes a reliability evaluation model for the IEGSs considering the effects of hydrogen. The power to hydrogen and methane process is proposed to convert the surplus renewable energy to hydrogen and methane, which are then mixed into natural gas systems. A novel optimal energy shedding model is proposed to explicitly account for the impact of hydrogen on the energy flow. To consider the temporal feature of renewable energy, a sequential Monte Carlo simulation approach is applied to evaluate the reliability of the IEGSs. The numerical simulations are performed on the integrated IEEE 24-bus and 20-node energy systems and the integrated IEEE 72-bus and 40-node energy systems to verify the effectiveness of the proposed model.


I. INTRODUCTION
C ONCERNS about climate issues and energy sustainability trigger the transition to a carbon-neutral society. At the Glasgow Climate Change Conference, over 100 countries have committed to reducing 30 percent of global methane emissions by the end of 2030 [1]. Renewable energy plays a vital role in achieving the ambitious goal [2]. The United States has budgeted more than $62 billion for developing clean energy [3]. In addition, green hydrogen produced using renewable energy through the water electrolysis process is pollution-free and can boost the utilization of renewable energy. Regarded as the cleanest and a versatile energy carrier, green hydrogen is gaining traction in achieving carbon-neutral societies. Meanwhile, the integration of green hydrogen tightens the coupling between power systems and natural gas systems, requiring additional effort in studying the integrated electricity-gas systems (IEGSs) with hydrogen.
Hydrogen has multiple applications (e.g., feedstock for industries and the source for fuel cells) [4]. Of the diverse applications, adding green hydrogen to natural gas provides a prominent way to increase renewable energy utilization and improve energy efficiency. Experiments indicate that slight additions of hydrogen to natural gas systems can improve fuel efficiency while ensuring combustion performance [5]. Adding less than 20% hydrogen to natural gas does not incur security issues [6]. However, hydrogen has different physical properties than natural gas. First, the molecular weight and density of hydrogen are much smaller than natural gas. Therefore, the gas flow dynamics would be changed with the integration of hydrogen [7]. Second, the volumetric energy density of hydrogen is around one-third of natural gas. Blending hydrogen into natural gas would reduce the energy density of the mixed gas. As a result, a higher flow rate is needed to meet the demand required, and the line pack will be reduced, which impairs the short-term operation security [8]. The impact of various hydrogen fractions in natural gas systems is investigated in [5] while the effects of hydrogen on gas flow are not considered. In [9], a two-stage model is proposed to study the operational impacts of various power-to-gas processes on the IEGSs. A transient gas flow model is adopted. The effects of hydrogen on gas flow are not analyzed. An optimal operation model for the IEGSs is proposed with the consideration of gas flow dynamics [10]. The effects of various parameters on gas flow dynamics are analyzed. However, the relation between hydrogen fractions and gas flow equations is not established. Since hydrogen physically differs from natural gas, adding hydrogen to natural gas systems would affect gas flow. Hence, it is urgent to investigate the effects of hydrogen on gas flow dynamics as well as the IEGSs.
On the other hand, reliability evaluation aims at assessing the ability of the systems to provide service without interruption [11]. There are many studies on power system reliability [12], [13], [14]. A novel inverse power system reliability evaluation is proposed to acquire the unknown component reliability parameters using the reliability indices [12]. An analytical model for power system reliability evaluation is proposed based on the sequential Monte Carlo simulation (SMCS) approach [13]. The proposed method enables fast reliability evaluation by partitioning the chronological system state sequence into several mutually exclusive events. Circuit breaker failure is considered in distribution network reliability evaluation based on an analytical approach [14]. Additional efforts have been carried out regarding the reliability evaluation of integrated energy systems [15], [16], [17], [18], [19], [20]. An integrated electricity-gas reliability evaluation model is proposed based on the SMCS approach [15], where steady-state energy flow is adopted. Results show that gas storage can improve the reliability of the IEGSs. A few studies consider dynamic energy flow [16], [17], [18], [19], [20]. Focused on cascading effects, a nodal reliability evaluation of the IEGSs is proposed [16], where the optimal power flow and optimal gas flow are separated and implemented hourly. Case studies demonstrate that considering dynamic cascading effects worsens the reliability of IEGSs. The reliability of integrated electricity-gas-heat systems is investigated in [17] based on the non-sequential Monte Carlo simulation, which shows that considering the dynamic feature of gas and heat systems can enhance the reliability of integrated energy systems. The optimal energy flow is implemented for one day without renewable energy. In [18], [19], [20], dynamic gas flow is modeled through highly nonlinear equations, which are solved independently for each hour using Newton-Raphson methods. Solution optimality and computational efficiency can not be guaranteed due to nonlinearity. The existing studies do not simultaneously address the solution optimality and computational burden issues when adopting dynamic energy flow models. Considering longer dispatch periods (e.g., one week) will exacerbate these drawbacks. Hydrogen storage usually has a large capacity, which can support weekly and even seasonal continuous operation [21]. In addition, the power-to-hydrogen and methane process is not included. Thus, the impact of power-to-hydrogen and methane on system reliability is unclear. Furthermore, the existing studies do not consider the effects of hydrogen on the reliability of the IEGSs. As previously discussed, hydrogen affects the energy flow of the IEGSs because of its different physical properties. Consequently, the reliability of the IEGSs would also be affected. Existing studies only focus on natural gas without hydrogen. Hence, the impact of hydrogen on the reliability of integrated energy systems is unknown.
To this end, this paper studies the reliability evaluation of the IEGSs considering the effects of hydrogen. The power-tohydrogen and methane (P2HM) process is proposed to convert the surplus renewable energy to hydrogen and methane, which are then mixed into natural gas systems. A novel optimal energy shedding model is proposed to explicitly account for the impact of hydrogen on the energy flow dynamics. A two-stage approach is proposed to quantify the effects of hydrogen and complete the linearization for the dynamic energy flow models. To consider the temporal feature of renewable energy, a sequential Monte Carlo simulation approach is applied to evaluate the reliability of the IEGSs. The contributions of this paper are shown as follows: 1) An optimal energy shedding model considering hydrogen is proposed. The P2HM process and its reliability model are proposed. The effects of hydrogen on energy flow dynamics are quantified. 2) A reliability evaluation model of the IEGSs with hydrogen is proposed. The SMCS approach is adopted to consider the temporal feature of renewable energy. A two-stage procedure including different linearization techniques is proposed to ensure computational tractability and accuracy.
3) The effects of the P2HM process and the fractions and physical properties of hydrogen on the reliability of the IEGSs are analyzed. Different energy models (i.e., steadystate and dynamic) and various hydrogen fractions (e.g., 0%, 10%, and 20%) are developed to investigate for the first time the effect of hydrogen on the reliability of IEGSs. The rest of this paper is organized as follows. The proposed optimal load shedding model is introduced in Section II. The SMCS-based reliability evaluation approach for the IEGSs is presented in Section III. The case studies and conclusions are presented in Section IV and Section V, respectively.

II. PROBLEM FORMULATION
This section first formulates the P2HM process and then develops the gas flow model considering hydrogen impact. Lastly, the optimal load shedding model is proposed.

A. Power to Hydrogen and Methane Model
Generally, hydrogen can supply electricity demand, supply hydrogen demand, or be injected into natural gas systems. Converting hydrogen back to electricity results in additional energy losses due to the relatively low round-trip efficiency of the power-to-hydrogen-to-power process. Since this paper is focused on the coupling between power systems and natural gas systems, hydrogen demand is not considered. Fig. 1 illustrates the P2HM systems, which are composed of electrolyzers, hydrogen storage, and methanation units. Instead of being curtailed, the surplus renewable energy is converted to hydrogen through the water electrolysis process within electrolyzers [5]. The generated hydrogen is then stored in hydrogen storage, which works as a buffer connecting hydrogen production and consumption to improve the security and efficiency of the P2HM systems. One salient feature of hydrogen storage is its enormous storage capacity, which can be up to the order of TWh [21]. By contrast, batteries typically store up to 10 MWh. Hydrogen can be either utilized to synthesize methane through the Sabatier reaction [22] or directly integrated into the natural gas systems [5]. The mathematical model of the P2HM process is as follows. Constraints (1.1)-(1.10) formulate the P2HM process. In the first step, surplus renewable energy is fed into electrolyzers and converted into hydrogen, modeled by (1.1). The input power is limited by the electrolyzer capacity as well as its working status (1.2). In the second step, the generated hydrogen is stored in the hydrogen storage for further utilization (1.3). Hydrogen storage is subject to capacity limits (1.4), (1.5), input limits (1.6), and output limits (1.7). In the third step, hydrogen is transformed into methane (1.8). Lastly, hydrogen and methane are mixed and then blended into the natural gas systems (1.9). Due to the combustion requirements, the volume of hydrogen mixed with methane should be limited (1.10).

1) General Gas Flow Equation:
The gas flow within natural gas pipelines can be formulated by three equations (i.e., state equations, continuity equations, and momentum equations) [23]. We start the analysis with the state equation, shown as follows.
where p and ρ are the pressure and density, respectively. The average molecular weight, M, and compressibility factor, Z, are determined by the composition of the mixed gas. Specifically, the average molecular weight can be calculated as follows [7].
where n i denotes the mole fractions of gas i. The compressibility factor can either be estimated based on the thermodynamics experimental data, or calculated as follows [7].
The value of pseudo-critical pressure and temperature can be determined as follows [7], [26].
The continuity equations and momentum equations are shown as follows.
where v denotes the gas velocity; G and H represent gravity force and pipe height, respectively.

2) Simplified Gas Flow Equation:
Equations (2.6) and (2.7) denote the continuity equations and momentum equations, respectively. The second term in (2.7) models the gravity effects, which can be neglected by assuming horizontal pipelines. The fourth and fifth terms in (2.7) represent inertia and kinetic energy, respectively. Because these two terms account for less than 1% of the solutions, they can also be neglected [23]. As a result, the momentum (2.7) is simplified as follows.
Further, by assuming T and Z remain constant for the pipeline [23], (2.9) can be written as follows.
Equations (2.10), (2.11) formulate the gas flow with the consideration of hydrogen. The effects of hydrogen on gas flow dynamics are quantified by M and Z, which can be calculated using (2.2), (2.3). Changing hydrogen fractions will alter the value of M and Z. Consequently, gas flow and power flow will be affected. The majority of the existing studies do not take into account hydrogen effects [15], [23], or just study natural gas systems with steady-state gas flow [7]. Nonetheless, the reliability of the IEGSs will be affected by hydrogen and gas flow dynamics, which will be verified in the case study section.
Equations (2.10), (2.11) are a set of partial differential equations, which can be discretized using the implicit finite difference methods [24]. To mitigate the computational burden, this paper adopts two diverse strategies to discretize the partial differential equations.
3) Dynamic Gas Flow Model: The partial derivatives with respect to length and time can be discretized through the backward difference methods [24].
The nonlinear item in (2.15) can further be linearized using the first-order Taylor series [24], shown as follows.
wherep P R nm,0 andf nm,t,0 are the initial operating points. Compared with steady-state gas flow models, in dynamic gas flow models, the inflow and outflow can be different such that a certain amount of gas can be stored in the pipes, which is called line pack and is defined as follows.

4) Steady-State Gas Flow Model:
As the inflow equals the outflow, the continuity (2.11) is dropped, and the following equation is adopted to discretize (2.10).
To summarize, the steady-state gas flow is modeled by (2.20), which can also provide the reference operating points for the dynamic gas flow model (2.16), (2.17). It can be seen that injecting hydrogen into natural gas systems will affect both the steady-state and dynamic gas flow as the physical properties of the mixed gas are altered (e.g., Z and M). Based on the above equations, the effects of hydrogen on gas flow dynamics can be quantified.
In general, Z is a variable associated withp P R nm and T based on (2.3). However, treating Z as a variable makes the optimal energy flow shedding a nonlinear programming problem, which can not be efficiently solved. The computational burden of reliability evaluation will also increase significantly as the optimal energy shedding model needs to be implemented thousands of times. Inspired by the fact that Z varies within a small range under similar operating conditions (e.g., pipeline pressures) [26], we first adopt the nonlinear steady-state model (2.20) to calculate Z for each hydrogen fraction based on (2.3)-(2.5). The obtained Z is then regarded as a parameter in the linearized optimal energy shedding model for reliability evaluation. Generally, higher hydrogen fractions lead to smaller M and bigger Z, and the impact of hydrogen fractions on M is greater than Z.

C. Optimal Load Shedding Model
This subsection formulates the optimal load shedding model for reliability evaluation.

1) Transmission System Constraints:
A DC power flow model is adopted to formulate the operation of power systems.
Constraint (3.1) denotes the nodal power balance. Constraint (3.2) restricts the curtailment of power demand. Constraint (3.3) indicates that the hydrogen can only be generated using available renewable energy, which is called green hydrogen. Constraints (3.4), (3.5) limit the power flow through transmission lines. Constraints (3.6)-(3.8) represent the feasible region of conventional power generators.
2) Natural Gas Systems Constraints: The dynamic gas flow model is shown as follows.
f P,out Constraint (4.1) indicates the nodal gas balance. Constraint (4.2) limits the gas shedding. Constraints (4.3)-(4.5) denote the feasible region of gas well output, node pressure, and flow through pipes. Constraint (4.6) restricts the node pressure of operating compressors. Constraints (4.7), (4.8) are the relaxed version of the continuity equations. Constraints (4.9), (4.10) are the relaxed version of the moment equations. In the case of pipe outages, the related continuity equations and moment equations are relaxed; otherwise, they equal (2.16), (2.17), respectively. For steady-state gas flow, the continuity (4.7), (4.8) should be dropped, and the moment (4.9), (4.10) are replaced by (2.20) with the same relaxed form.
3) Optimal Load Shedding Model for the IEGSs: To calculate reliability indices, the system adequacy should be evaluated for each system state, which can be done using an optimal load shedding model [11]. Therefore, we propose an optimal load shedding model for the IEGSs. The objective is to minimize the amount of power and gas load shedding and renewable energy curtailment, subject to various operational constraints.  (5) for each dispatch block (e.g., one week), the amount of load and renewable curtailment can be determined as well as the reliability indices. The specific value of load shedding costs depends on the cost of energy interruptions and the priority of the electric buses and gas nodes, which can be estimated based on the interrupted energy assessment [30].
The time step for discretizing the partial differential equations affects the accuracy and complexity of the proposed reliability evaluation approach. As discussed and demonstrated in [23], [24], the hourly time scales can balance accuracy and efficiency. In addition, reliability evaluation focuses on energy adequacy assessing whether the demand can be supplied in case of outages, which is typically conducted on an hourly basis [11]. Hence, the time resolution in this paper is set as one hour to balance accuracy and computational burden. Note that short time scales are still applicable to the proposed approach at the cost of increased computation time.
Unlike unit commitment which centers on the economic aspect, reliability evaluation cares more about the security aspect. Specifically, reliability evaluation aims to assess whether and how the demand can be satisfied in the case of outages. In the contingency situation, the foremost objective is to meet the demand rather than minimize the operating costs. Hence, the shutdown cost, start-up cost, and minimum up and down time constraints, which are widely considered in the unit commitment problems, are not modeled in the proposed reliability evaluation process. This assumption is widely adopted in reliability evaluation [11] and also in line with other similar studies [13], [15].

III. RELIABILITY EVALUATION
The section presents the reliability evaluation procedure for the IEGSs based on the proposed optimal load shedding model.

A. Component State Generation
The first step to conduct reliability evaluation is to build component outage models and sample component states. This paper adopts the two-state outage model [11]. That is, the component i is either down or up.
In the sequential Monte Carlo simulation approach, given the current component state and reliability parameters, the state duration can be sampled as follows [11].
where U i is uniformly distributed random variables in (0, 1). The sequential Monte Carlo is based on the well-known Monte Carlo method. Based on the principle of the Monte Carlo method, a random variable in (01) is first sampled, and then the state duration of the component is determined using (6.2) or (6.3), depending on the current status of the sampled component. For the operating components (e.g., y i,t equals 1), (6.2) is used to sample the up time; otherwise, (6.3) is used to calculate the failure time. The up and down states will appear alternately during the sampling year. Once all the components are sampled, the corresponding system state duration can be determined. The above formulations are applied to all the components in the IEGSs, except for the wind farms.

B. Reliability Indices
Reliability indices are proposed in this subsection to evaluate the reliability level of power systems, renewable energy, and natural gas systems. Provided the sample size is large enough, we use the mean values of the samples to calculate reliability indices. On the power system side, the loss of load probability (LOLP) and the expected demand not supplied (EDNS) are defined as follows [11].
where s refers to the sampling year and 1(•) is the indicator function having the following form: Accordingly, these indices can be extended to each bus.
Similarly, the loss of renewable energy probability (LORP), the loss of gas probability (LOGP), the expected gas not supplied (EGNS), and the expected renewable energy not supplied (ERNS) are defined as follows. The coefficient of variation (COV) is then used to judge if the reliability evaluation algorithm converges [11]. As the proposed reliability indices cover electricity, renewable energy, and gas, the following criterion is adopted.

C. Reliability Evaluation Approach
This paper adopts the SMCS approach to evaluate the reliability of the IEGSs. The SMCS approach utilizes chronological data to evaluate the reliability, which covers one year on an hourly basis [11]. The state duration for each component is sampled using (6.2)-(6.3), and then the system state can be analyzed using optimal load shedding models. Since the reference points used for the first-order Taylor series are closely related to the system state, which is randomly generated, it is necessary to determine the reference points for each system state instead of using the common values [24]. Therefore, a two-stage optimal energy shedding approach is developed based on (5). Specifically, given the system outage state, the steady-state gas flow is conducted first to calculate the reference points and fix the gas flow directions. Then, the dynamic gas flow is performed to optimize the amount of load shedding and renewable curtailment. For each sampling year, the reliability indices are calculated using (6.4)-(6.16). After obtaining the reliability indices, (6.17) is used to check the convergent status. The detailed procedure for the IEGS reliability evaluation is given in Algorithm 1 and the flowchart is shown in Fig. 2.

IV. CASE STUDY
This section verifies the effectiveness of the proposed reliability evaluation approach for the IEGSs. The IEGSs are composed of a modified IEEE 24-bus power system [11] and a 20-node natural gas system [23], as shown in Fig. 3 for each block do 7 Solve the steady-state model based on (5) by dropping the continuity equations and using the linearized version of (2.20) 8 Given the reference points and gas flow directions, solve the dynamic model based on (5) 9 end for 11 Calculate reliability indices based on (6.4)-(6.16). 12 Check stopping conditions using (6.17). 13 end for 14 Output the final reliability indices. 5, 8, and 12 of the natural gas systems, respectively. Gas-fired generators are placed at buses 1, 13, and 15 of the power systems, which are connected to nodes 5, 14, and 2 of the natural gas systems, respectively. The 20-node natural gas system consists of two gas wells, 19 pipes, and two compressors. The reliability parameters are based on [15]. The yearly normalized electricity demand, gas demand, and wind farm output are depicted in Fig. 4. To consider the renewable uncertainty, representative scenarios are adopted. In this paper, a longer scheduling period for each dispatch block (i.e., 168 hours) is adopted to fully utilize the large capacity of hydrogen storage, which is equivalent to one week. As a result, a year with 8736 hours, rather than 8760 hours, is considered, which can be exactly divided into 52 weeks. The annual load curve with 8736 hourly load points is obtained from the IEEE RTS [11]. The maximum number of sampling years is set as 50. The minimum coefficient of variation for stopping the algorithm is set as 0.04. All the detailed parameters used for case  studies can be found in the online warehouse [25]. The proposed reliability evaluation approach is implemented in Python, and the optimal load shedding model is solved using Gurobi 9.1 with default settings. All the simulations are implemented on a PC with 16GB RAM and a 3.6 GHz Intel Core i7 processor.
Remark 4.1: Unlike the economic dispatch, the proposed optimal energy shedding model is to minimize the energy shedding in the case of outages rather than minimizing the generation costs and allocating the profits. Hence, the clearing mechanism of electricity markets is not mandated, and the duration of each dispatch block can be set as 168 hours, considering the large size of hydrogen storage.

A. Reliability Evaluation of the IEGSs
This subsection reports the reliability evaluation results of the IEGSs. Due to security concerns, the maximum hydrogen fractions are set as 20%. To demonstrate the merit of the proposed model, the following three cases with 20% hydrogen fractions are developed.
C0: dynamic gas flow with P2HM and hydrogen effects. C1: steady-state gas flow with P2HM, but no hydrogen effects.
This case works as the baseline. C2: dynamic energy flow with hydrogen effects, but no P2HM. Further, the reliability indices of the three cases are compared in Table I. It can be seen that the best reliability performance is observed in the proposed model. The worst LOLP, LOGP, GWh/yr. to 0.0004 and 0.3592 GWh/yr., respectively. In addition, natural gas systems receive more enhancement than power systems. These results demonstrate that the proposed model can enhance the system reliability of the IEGSs and maximize the utilization of renewable energy.
Further, the nodal EDNS, ERNS, and EGNS are compared. Table II reports the nodal EDNS of all the load buses in the three cases. It can be observed that the nodal EDNS is related to certain cases. For the baseline case, the maximum and minimum nodal EDNS appear in buses B18 and B1 with 32.5142 GWh/yr. and 0.0805 GWh/yr., respectively. By contrast, the maximum and minimum nodal EDNS in the proposed model reduce to 31.2064 GWh/yr. and 0.0652 GWh/yr., respectively. The bus with the maximum nodal ENDS is shifted to B19. Although the maximum nodal EDNS in C2 is smaller than C0 and C1, its minimum nodal EDNS is also larger.    Table IV compares the nodal EGNS in the three cases. Compared with C1 and C2, the proposed model improves the EGNS of all the load nodes in the natural gas systems, which is different from the power system side. The maximum and minimum nodal EGNS are reduced from 123.1283 MSm 3 /yr. and 0.8139 MSm 3 /yr. to 117.3800 MSm 3 /yr. and 0.4487 MSm 3 /yr., respectively. In addition, the nodes close to P2HM (e.g., 6, 11, 15, and 16) achieve higher reliability improvement.

B. Impact of Gas Flow Dynamics on the IEGS Reliability
This subsection analyzes the impact of gas flow dynamics on the IEGS reliability. The following two cases are developed.
C3: dynamic gas flow with hydrogen effects; C4: steady-state gas flow with hydrogen effects; The hydrogen fractions increase by 5% from 0% up to 20% for each case. Fig. 6 illustrates the effects of gas flow models on the reliability of the IEGSs at different hydrogen fractions. Since LORP and ERNS are not sensitive to the gas flow models, they are not depicted. It can be seen that considering gas flow dynamics reduces LOLP, LOGP, EGNS, and EDNS, which indicates better reliability. Another important observation is that hydrogen fractions have diverse impacts on reliability in these two cases. In the case of dynamic gas flow models (i.e., C3), increasing hydrogen fractions improves LOLP and EDNS but exacerbates LOGP and EGNS. By contrast, the opposite trends in LOLP, LOGP, and EGNS are observed when steady-state gas flow models are adopted (i.e., C4). As a result, the reliability gap between C3 and C4 narrows with the increase of hydrogen fractions.

C. Impact of Hydrogen Properties on the IEGS Reliability
This subsection analyzes the impact of hydrogen physical properties on the IEGS reliability. The following case is added.
C5: dynamic gas flow models without hydrogen effects; hydrogen fractions increase by 5% from 0% up to 20%. Fig. 7 illustrates how the hydrogen physical properties affect reliability at different hydrogen fractions. It can be seen that considering hydrogen effects escalates LOLP, LOGP, EDNS, and EGNS, which means the reliability of IEGSs is impaired. The reliability gap between C3 and C5 expands with the increase of hydrogen fractions. Specifically, when hydrogen physical properties are not included, increasing hydrogen fractions lowers LOLP, LOGP, EDNS, and EGNS. The trends in LOGP and EGNS are opposite to C3.

D. Case Study in the IEEE 72-bus Power Systems and 40-Node Natural Gas Systems
In this subsection, the proposed reliability evaluation approach is applied to the IEEE 72-bus power systems and 40-node gas systems. This system is constructed by interconnecting three IEEE 24-bus power systems and two 20-node gas systems, as shown in Fig. 8. On the power system side, the generators and electrolyzers of zone 1 and 3 are identical to the IEEE 24-bus test system used in the previous subsection. For zone 2, three wind farms are placed at bus 221, 222, and 223, respectively, and all the conventional generators are coal-fired. Zone 1 is interconnected with zone 2 through tie-lines 23-217, 13-215, and 7-203, and zone 3 through tie lines 21-325, 325-212. Zone 2 is interconnected with zone 3 through tie-line 223-318. On the natural gas system side, the settings of zone 1 and zone 3 keep the same as the 20-node test systems used in the previous subsection. These two zones are interconnected through pipe N9-N309. The maximum number of sampling years is set as 40.
The minimum COV is set as 0.03.
First, the following three cases with 20% hydrogen fractions are developed to compare the reliability performance.
D0: dynamic gas flow with P2HM. D1: steady-state gas flow with P2HM, but no hydrogen effects.
This case is the baseline. D2: dynamic energy flow without P2HM. Table V compares the reliability indices of the three cases. It can be seen that the proposed model achieves the best reliability performance. Compared with the baseline D1, LOLP, LOGP, EDNS, and EGNS in the proposed model are improved by 1.39%, 15.67%, 2.75%, and 24.63%, respectively. Compared with D2, where P2HM is not considered, LOLP, LOGP, EDNS, and EGNS in the proposed model are improved by 1.56%, 4.61%, 2.05%, and 14.42%, respectively. The most remarkable improvement is LORP and ERNS, which are improved by 95.61% and 99.33%, respectively. In addition, natural gas systems gain more improvement than power systems.
Table VI compares the nodal EDNS of several representative buses in various cases. Specifically, the buses with the best and worst EDNS in each case and the buses connected to natural gas systems are selected as the representative buses. It can be seen that nodal EDNS depends on the cases. In D1, the minimum and maximum EDNS are 0.1147 GWh/yr. and 31.0833 GWh/yr., respectively, which are observed at buses 201 and 118. By contrast, buses 216 and 106 in D0 have the minimum and maximum EDNS, respectively. D2 has the same trends as D0. The maximum and minimum EDNS obtained by the proposed model are smaller than the other two cases.      Table VI, the maximum and minimum EGNS obtained by the proposed model are smaller than the other two cases. In addition, the nodes close to P2HM (e.g., 6, 16, 315, and 316) achieve higher reliability improvements.
Further, the following cases are developed to study the effects of gas dynamics and hydrogen physical properties on the reliability of the IEGSs at various hydrogen fractions. The hydrogen fractions increase by 10% from 0% up to 20%. D3: dynamic gas flow with hydrogen effects. D4: steady-state gas flow with hydrogen effects. D5: dynamic gas flow without hydrogen effects. Fig. 9 illustrates how the gas flow dynamics affect the reliability of the IEGSs at various hydrogen fractions. It can be observed that the proposed model, D3, obtains the smallest LOLP, LOGP, EDNS, and EGNS. In addition, the hydrogen factions have various effects on the reliability indices. In D3, as the hydrogen fractions grow, the LOLP and EDNS reduce while LOGP and EGNS grow. In D4, increasing hydrogen fractions results in smaller LOLP, EDNS, LOGP, and EGNS. Fig. 10 illustrates how the hydrogen physical properties affect the reliability of the IEGSs at various hydrogen fractions. In this case, it can be seen that considering hydrogen physical properties lowers LOLP but elevates LOGP, EDNS, and EGNS. As the hydrogen fractions grow, LOLP, LOGP, EDNS, and EGNS in D5 decline, and the reliability gap between D3 and D5 increases.

V. DISCUSSION
This section discusses the accuracy of the linearization methods in dynamic gas flows, concerns on practical system applications, and the impact of hydrogen storage size on system reliability.

A. Accuracy of First-Order Taylor Series-Based Linearization Method
To numerically demonstrate the accuracy of the proposed approach, in the modified IEEE 24-bus power systems and 20-node gas systems, we compare the reliability indices obtained by the following approaches.
As SOCP can guarantee global optimality, M2 works as the benchmark. We consider a one-year simulation and then calculate the reliability indices. The time horizon for each dispatch block is changed to 24 hours as the SOCP-based approach can not obtain results within 4 hours when the scheduling period is set as 168 hours. To make sure that the hourly system state is the same for the different approaches, the random seed is fixed. Below reports the comparison.
It can be seen that the proposed approach, M1, is very close to M2 while greatly reducing the computation time. By contrast, the conventional Taylor series linearization approach, M3, which is based on fixed reference points, leads to enormous errors, although it has the least computation time. Note that the computation time is only for the one-year simulation (i.e., one iteration). In practice, reliability evaluation requires dozens of iterations to converge. It can be concluded that by combining different linearization strategies, the proposed two-stage approach can ensure the accuracy of reliability evaluation while improving computational tractability.

B. Applications to Practical Systems
In this part, we present the solution time of the proposed reliability evaluation approach for the two test systems and also discuss the scalability and computational performance of the proposed approach when applied to practical systems. Table X compares the solution time of various cases in the two test systems. S1 refers to the modified IEEE 24-bus power systems and 20-node gas systems, and S2 refers to the modified IEEE 72-bus power systems and 40-node gas systems. It can be seen from Table X that time spent on evaluation accounts for the largest portion and the total solution time increases with the system size. However, the time on each iteration is acceptable. Note that, unlike economic dispatch or unit commitment, reliability evaluation is not integrated into the energy management system and thus does not need online implementation or has strict requirements on calculation speed.
As can be seen from Table X, the computational burden in reliability evaluation primarily lies in the evaluation process, which repeatedly implements the optimal energy shedding model for various system states. For each iteration, this procedure will be implemented dozens of times. A number of iterations are needed to ensure convergence. Hence, there are two possible ways to enhance the scalability of the proposed method. The first way is to reduce the time spent on calculation for each iteration. This way can be done by reducing the system size by clustering the sub-regions that have similar consequences or are out of interest as one equivalent node [27]. Although the accuracy will be slightly reduced, computational tractability is ensured. The second way is to reduce iterations. This way aims to improve the convergence rate by adopting advanced sampling strategies [28], [29]. For instance, combining cross-entropy optimization and importance sampling techniques can greatly improve the sampling efficiency and thus reduce the required iterations to converge [28]. Note that the mentioned strategies can be directly integrated into the proposed approach and the evaluation process.

C. Impact of Hydrogen Storage Size
To analyze the impact of hydrogen storage size on system reliability, the size factor of hydrogen storage size is set as 0.25, 0.50, 1.0, 1.5, and 2.0. The third case is the same as C3 where the dynamic gas flow with hydrogen effects is adopted. For each case, the hydrogen fraction increases by 10% from 0% to 20%. Other parameters remain the same as the modified IEEE 24-bus power systems and 20-node gas systems. Fig. 11 reports the effects of hydrogen storage size on system reliability with respect to different hydrogen fractions in the modified IEEE 24-bus power systems and 20-node gas systems. It can be seen that increasing hydrogen storage size can improve system reliability. The largest enhancement is observed when the size factor increases from 0.25 to 0.5 and from 0.5 to 1.0. However, when further increasing the hydrogen storage size to 1.5 times or 2.0 times of the original size, the enhancement becomes marginal. In addition, for a specific scale factor, increasing hydrogen fractions improves the reliability of power systems but worsens the reliability of natural gas systems, which is consistent with the previous case studies. Fig. 11. The effects of hydrogen storage size on system reliability with respect to different hydrogen fractions.

VI. CONCLUSION
This paper proposes a reliability evaluation model for the IEGSs considering the effects of hydrogen. The power to hydrogen and methane process is proposed to convert the surplus renewable energy to hydrogen and methane, which are then blended into natural gas systems. A novel optimal energy shedding model is proposed to explicitly account for the impact of hydrogen on the energy flow. Then a two-stage approach including different linearization techniques is proposed to ensure computational tractability and accuracy. To consider the temporal feature of renewable energy, the sequential Monte Carlo simulation approach is applied to evaluate the reliability of the IEGSs.
Based on the numerical simulations in the two IEGSs with different scales, the following conclusions are drawn: 1) P2HM can enhance the system reliability of the IEGSs and maximize the utilization of renewable energy. In the first test system, P2HM reduces LORP and ERNS by 99.72% and 99.93%, respectively. In the second test system, P2HM reduces LORP and ERNS by 95.61% and 99.33%, respectively. Regarding the nodal reliability, the nodal EGNS of all the load nodes can be improved, especially for those close to P2HM. For EDNS, the influence varies with cases. 2) Gas flow dynamics can enhance the reliability of the IEGSs. In the first test system, dynamic gas flow models reduce LOLP, LOGP, EDNS, and EGNS by 0.53%, 22.27%, 0.87%, and 9.25%, respectively. In the second test system, dynamic gas flow models reduce LOLP, LOGP, EDNS, and EGNS by 1.39%, 15.67%, 2.75%, and 24.63%, respectively. 3) Adding hydrogen worsens the reliability of natural gas systems due to the hydrogen physical properties. The detrimental influence increases with growing hydrogen fractions. When hydrogen fractions rise from 0% to 20%, LOGP and EGNS increase by about 2.74% and 0.58% in the first test system, respectively, and by 2.45% and 3.03% in the second test system. By contrast, the influence of hydrogen on power systems are case dependent. LOLP increases by 0.01% in the first system but decreases by 0.05% in the second system. EDNS decreases by 0.07% and 0.06% in the two systems, respectively. Generally, hydrogen has a more significant impact on natural gas systems than power systems. The following directions can be future work: 1) Reliability evaluation of integrated power-gas-hydrogen systems. Since electricity, natural gas, and hydrogen can be converted to each other, reliability evaluation of integrated power-gas-hydrogen systems can help determine the optimal penetration level of hydrogen.
2) The coordinated operation and planning of the integrated energy systems considering multiple hydrogen pathways [32]. Hydrogen supply chains contain multiple choices for production, transmission, and distribution. Incorporating various technologies into the integrated energy systems is critical to constructing a low-carbon society. 3) Energy trading with distributed energy resources, responsive demand and green hydrogen [33]. Since production costs of green hydrogen are still high, incorporating green hydrogen into energy trading with distributed energy resources and responsive demand can enhance the profitability of green hydrogen.