On the Integration of Hydrogen Into Integrated Energy Systems: Modeling, Optimal Operation, and Reliability Assessment

The proliferation of renewable energy (RE) brings tremendous challenges to integrated energy systems (IESs). Converting RE into hydrogen, one of the cleanest energy carriers, provides an appealing alternative for decarbonized IESs. Among the various hydrogen applications, blending hydrogen into natural gas systems is already applicable. However, the disparate hydrogen physical properties trigger concerns about hydrogen integration. This paper investigates the integration of hydrogen into the IESs, focusing on blending hydrogen into natural gas systems. The literature on hydrogen modeling, control, operation, planning, and markets are first reviewed. Based on the convex combination methods, a power-to- hydrogen-heat-methane (P2HHM) model with unit commitment is proposed. The steady-state and dynamic gas flow models with the explicit consideration of hydrogen effects are then proposed. Two applications are given based on the proposed hydrogen modeling. A Wasserstein metric-based distributionally robust optimal operation model is proposed first based on the developed P2HHM and gas flow models. Then a sequential Monte Carlo simulation-based reliability assessment model is formulated to analyze the effects of hydrogen physical properties and hydrogen fractions on the optimal operation and reliability of the IESs. Numerical simulations are conducted to analyze the integration of hydrogen and verify the effectiveness of the proposed model.


I. INTRODUCTION
C ONCERNED about global warming and energy sustainability, many countries have proposed ambitious targets for greenhouse gas emission reductions. For instance, by 2030, the U.S. plans to achieve a 50-52 percent reduction in greenhouse gas emissions below 2005 levels [1]. To fulfill the target, decarbonized energy systems are forming where renewable energy (RE) plays a crucial role in the electricity supply. As of 2020, the global RE capacity reaches 2.9 TW, 52% of which is wind and solar energy [2]. Despite the low operation costs and zero carbon dioxide emissions, curtailment of RE frequently occurs because of the intermittence and uncertainty from RE. For instance, the total curtailment of wind and solar generation in China is more than 49 TWh in 2017 [3]. The increasing amount of RE curtailment adversely affects the financial performance and secure operation of the energy systems.
One promising way to address the issues raised by the high penetration of RE is the power-to-hydrogen (P2H) technology. Specifically, the excess RE can be utilized to produce hydrogen through the water electrolysis process. The produced hydrogen is called green hydrogen [4]. Since the only by-product of combusting hydrogen is water and carbon dioxide can be fully utilized in the subsequent methanation process, the production of green hydrogen is free of pollution and greenhouse gas emission. Green hydrogen is gaining attention in the middle of constructing decarbonized energy systems. Hydrogen works as an energy carrier used for various purposes, including feedstocks for industries [5], electricity storage [6], energy resources for electric vehicles [7], and direct injection to natural gas networks [8].
Although green hydrogen is already competitive in a few applications by co-operating the wind farms and P2H facilities (e.g., in Germany and Texas electricity markets) [9], the operating expenses of green hydrogen are still high compared with grey hydrogen and blue hydrogen. Hence, green hydrogen is not readily available for the industrialscale supply [9]. Another green hydrogen application is to generate electricity, which is equivalent to electricity storage. The excess RE is transformed into hydrogen and then stored in hydrogen tanks, which are connected to fuel cell units. When the generation is insufficient, or the electricity price is high, fuel cell units can generate electricity using the stored hydrogen. However, this method is not cost-effective as the round-trip efficiency of power-to-hydrogen-to-power is relatively low [6]. On the other hand, hydrogen-based fuel cell vehicles are seen as a promising way to harness green hydrogen. In comparison with battery electric vehicles, hydrogen-powered electric vehicles can travel long range with zero greenhouse gas emissions [7]. Nonetheless, fuel cell vehicles are premature compared with other electric vehicles. Blending green hydrogen into natural gas systems is receiving attention. Demonstrations in German gas grids show that a certain percentage of hydrogen mixed with natural gas can improve fuel efficiency while guaranteeing combustion performance [8]. Adding no more than 20% hydrogen to natural gas networks is acceptable [10]. Results in [11] show that blending hydrogen into integrated energy systems (IESs) can save operation costs and reduce energy curtailment. Apart from these merits, adding hydrogen would also trigger concerns about the secure operation of the IESs because hydrogen has different physical properties than natural gas [12]. First, hydrogen has a smaller molecular weight and density than natural gas, resulting in changes in the gas flow dynamics [12]. Second, from the perspective of energy content, the volumetric energy density of hydrogen is around one-third of natural gas. The density of the mixed gas will be smaller than the natural 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 [13]. Furthermore, when metals are exposed to hydrogen, hydrogen embrittlement is likely to occur, causing failure, leakage, and even explosion [14]. However, replacing existing natural gas pipelines to accommodate hydrogen is expensive [15]. Hence, it is vital to address the economic and safety issues of the IESs with hydrogen integration. However, there are few studies focused on the effects of green hydrogen on the optimal operation and reliability assessment of the IESs. Regarding the optimal operation, existing studies adopt either energy balance models [27] or conventional gas flow models [22]. Hydrogen effects on energy flow are not investigated. Gas flow with hydrogen is considered in [21], while how hydrogen affects the operation of IESs is not revealed, and renewable uncertainties are not included. However, as demonstrated in case studies, both hydrogen effects and uncertainties affect the optimal operation of integrated energy systems. Regarding the reliability assessment, there are several studies on power systems [53], integrated powergas systems [54], and integrated power-gas-heat systems [55]. However, green hydrogen is not considered. As discussed above, hydrogen affects the energy flow of IESs and thus the reliability of IESs. As a result, the impact of hydrogen on the reliability of IESs is unclear.
To this end, this paper aims to investigate the integration of hydrogen into the IESs, with an emphasis on the impact of blending green hydrogen into natural gas systems. The objective of this paper is to reveal how hydrogen affects the operation and reliability of the IESs with respect to various hydrogen fractions. The proposed models and conclusions can assist system operators to optimally schedule the IESs and enhance reliability while maximizing hydrogen integration. In order to achieve this goal, we first review the literature in the field of hydrogen modeling and control, operation, planning, and markets. Based on convex combination methods, we propose a powerto-hydrogen-heat-methane (P2HHM) model with unit commitment (UC) to ensure model accuracy while reducing the computational burden. We then develop steady-state and dynamic gas flow models with the explicit consideration of hydrogen effects. Finally, targeted at the economic and security issues, two case studies related to the optimal operation and reliability assessment of the IESs with green hydrogen integration are carried out, which have not been reported before. The contributions of this paper are drawn as follows: 1) A P2HHM model with UC. To ensure accuracy while reducing computational complexity, a linear P2HHM model is proposed based on the convex combination methods. UC constraints are then integrated to model the operational characteristics (e.g., loading range, startup, and shutdown limits) of electrolyzers (ELZs). 2) Steady-state and dynamic gas flow models with the explicit consideration of hydrogen effects. The general gas flow partial differential equations (PDEs) are first formulated. Then the PDEs are discretized based on implicit finite difference methods with the explicit consideration of hydrogen effects. Finally, the discretized gas flow models are convexified to mitigate the computational burden. 3) Numerical analysis of hydrogen effects on IES optimal operation and reliability. Wasserstein metric-based distributionally robust optimal operation models are proposed to optimize the operation of the IESs. Sequential Monte Carlo simulation-based approaches are developed to evaluate the reliability of the IESs. These two models are proposed for the first time to consider hydrogen effects. The rest of this paper is organized as follows. Literature on hydrogen modeling and control, operation, planning, and markets is reviewed in Section II. The UC model for the P2HHM process is formulated in Section III. The convexified steady and dynamic gas flow models are derived in Section IV. The case study and conclusion are presented in Section V and Section VI, respectively.

II. LITERATURE REVIEW
In recent years, great interest has been shown in hydrogen technology, including hydrogen modeling and control, operation, planning, and markets, due to the crucial role of hydrogen in energy transitions and the construction of carbon-neutral societies. This section gives a brief literature review on these fields.

A. HYDROGEN MODELING AND CONTROL
The hydrogen systems typically consist of electrolyzers, hydrogen tanks, and fuel cells [17]. The modeling [16] and control strategy [17], [18], [19] of these components affect the operation of the hydrogen systems. In [16], the modeling and simulation of alkaline electrolyzers are studied, considering the electrothermal characteristics. Through the modeling and simulation, the relation between cell voltage, current, temperature, and pressure of electrolyzers can be characterized, which facilitates the studies on the electrolyzer operation and planning problems. A control scheme is proposed for the islanded microhydro plants [17], where the control models of electrolyzers, hydrogen tanks, and fuel cells are established. A model predictive control-based control strategy is proposed for hydrogen plants with reversible solid oxide cells [18], which can operate in either fuel cell mode or electrolysis mode. Linear parameter-varying methods are adopted to achieve fast control. In [19], a pressure control strategy for alkaline electrolysis is proposed to expand the loading range of alkaline electrolyzers. The proposed operation curve tracking controllers and model predictive controllers can reduce the minimum load to 10.5% and 10%, respectively. An optimal hierarchical control strategy of the P2H process is proposed in [20] to improve the green hydrogen production efficiency while mitigating the impact of renewable energy on electrolysis current. VOLUME 9, 2022

B. HYDROGEN OPERATION
Several studies focus on the effect of hydrogen injection into natural gas systems [11], [21], and [22]. The optimal operation of IESs considering hydrogen injection into natural gas systems is studied in [11]. The impact of various hydrogen fractions on the optimal operation is investigated. A twostage model is proposed to assess the effects of integrating the power-to-gas process into IESs [21]. 0.1% volume of hydrogen is added to natural gas systems. The effects of hydrogen on voltage security, reactive power, and pressure are studied in [22]. Four gas quality indicators (i.e., Wobbe index, combustion potential, specific gravity, and gross calorific value) are adopted to determine the hydrogen fractions. Additional efforts are made in the operation of IESs with hydrogen [23], [24], [25], [26], [27]. The optimal operation for the distribution networks and district heat networks is studied in [23], where a power-to-hydrogen-heat model is proposed to reduce the output of combined heat and power plants (CHPs) and renewable energy curtailment. The Nash bargaining-based cooperative operation for wind generators and hydrogen fueling stations is proposed in [24]. More than 30% reductions in the operation costs of hydrogen fueling stations are achieved. A two-stage optimal operation model of the integrated electricity-hydrogen energy systems is proposed to coordinate the power systems and the hydrogen generation, transportation, and storage [25]. The scheduling of hydrogen systems is modeled as the vehicle routing problems. A coordinated operation model of power and hydrogen networks considering hydrogen vehicles is proposed in [26]. Sequential convex approximation methods are adopted to handle the nonconvexity stemming from hydrogen hydraulic models. The stochastic optimal operation of hydrogen-based IESs is investigated in [27], considering fuel cell-based CHPs and hydrogen tanks.

C. HYDROGEN PLANNING
Effective investment in hydrogen-related systems can promote the profitability of the IESs, especially in the context of high hydrogen operating costs. There are several works on the optimal planning of hydrogen systems [28], [29], the IESs [30], [31], and hydrogen supply chains [33], [34]. The capacity of high-temperature electrolyzers is optimized based on a two-stage stochastic programming (SP) model in [28]. Various load conditions and power inputs (i.e., hydro and wind power) to electrolyzers are considered. Results show that the net revenue raises by 29% and 71% with hydro and wind as input power, respectively. The optimal capacity planning for power-to-gas plants including electrolyzers, hydrogen storage (HS), and methanation reactors is studied in [29]. At most 17% of synthetic natural gas (SNG) production costs are saved with the optimal planning decisions. In [30], the capacity of electricity-hydrogen integrated energy systems is optimized, considering the power-to-hydrogen-heat process and hydrogen storage. Two types of hydrogen storage are considered. The daily hydrogen storage operates on an intra-day basis while the seasonal hydrogen storage operates on an inter-day basis. Numerical results show that integrating hydrogen into the IESs can improve economic performance in the case of high renewable energy penetration. The optimal planning of electrolyzers and electric boilers in the IESs is studied in [31], where the renewable uncertainties are considered using distributionally robust chance-constrained models. The total costs are minimized and wind generators benefit the most. A tri-level planning model of IESs is proposed in [32] to optimally invest in the energy hub and power, gas, and hydrogen networks. Studies show that carbon emission of the IESs can be reduced. The optimal sizing of electrolyzers and seasonal hydrogen storage is investigated in [33], considering the linkage between interregional hydrogen supply chains and power transmission networks. The proposed model raises the working hours of electrolyzers and increases the net present value by reducing capital costs and transportation costs. A hydrogen supply chain model is proposed in [34] to minimize the capital costs of hydrogen generation, storage, transmission, and compression components. The scheduling of hydrogen trucks and pipelines is incorporated. A 9% reduction in hydrogen costs is achieved by leveraging the transmission and storage capability of trucks. The model is extended in [35] by further considering the power transmission systems. Uncertainties from load and renewable energy are considered using robust optimization (RO). The proposed co-planning model can reduce the overall costs and renewable curtailment. The co-planning of hydrogen supply chains and hydrogen refueling stations is studied in [36] with the objective of minimizing the least cost of hydrogen. Multiple technology options for each component are considered.

D. HYDROGEN MARKET
Establishing hydrogen markets can exploit the potential of hydrogen supply and consumption and maximize social welfare [37]. However, unlike the electricity markets, hydrogen markets are still at the developing stage [38]. Several studies have been reported on the hydrogen markets [38], hydrogenelectricity markets [37], [39], [40], and hydrogen-electricitygas markets [41]. An equilibrium model for the regional hydrogen markets is proposed, considering renewable energy and transportation costs [38]. Distributionally robust optimization (DRO) is applied to include the uncertainties from renewable energy. A distributed energy sharing model for hydrogen and electricity is proposed to coordinate powerto-gas plants, electric vehicles, and hydrogen vehicles [37]. The total social welfare is maximized. A local energy market for renewable generators, loads, hydrogen vehicles and hydrogen storage systems is established to trade electricity and hydrogen [39]. The merit order principle is applied to settle the market-clearing price. Renewable energy integration, peak load shaving, and utility improvement are boosted through the proposed local market. A dynamic hydrogen pricing strategy and a capacity-based demand response program are developed to schedule the hydrogen storage stations [40], which simultaneously provide energy for the transportation systems and the electricity market. Remarkable profits are made by investing in the hydrogen storage stations and additional profits can be made through demand response programs. A bi-level strategic bidding model for powerto-hydrogen-methane plants is developed, considering multitype markets (i.e., hydrogen, electricity, and natural gas) [41]. The hydrogen market is simulated through a quantity-based competition Cournot model. The hydrogen market becomes dominant with high hydrogen sale prices (e.g., 40 /kg).

III. POWER TO HYDROGEN-HEAT-METHANE MODELING
This section presents the P2HHM process and formulates the P2HHM model with UC based on convex combination methods.

A. POWER TO HYDROGEN-HEAT-METHANE MODELING
There are three types of hydrogen applications (i.e., supply electricity demand, supply hydrogen demand, and inject it into natural gas systems). As converting hydrogen back to electricity is financially ineffective, this paper aims at the hydrogen applications to supply hydrogen demand and inject it into the natural gas systems. Due to the diversities of hydrogen demand, we generally model it as hydrogen consumers.  illustrates the configuration of the P2HHM system, which consists of electrolyzers, water circulation systems, hydrogen storage, methanation reactors, and other auxiliary equipment. The whole system can be divided into a powerto-hydrogen-heat subsystem [23], a hydrogen storage subsystem, and a hydrogen-to-methane subsystem. The whole system takes RE as inputs and outputs hydrogen, heat, and SNG. Specifically, instead of converting RE to hydrogen exclusively, which is common in P2H models, P2HHM can generate hydrogen and heat simultaneously in the powerto-hydrogen-heat subsystem through electrolysis and heat transfer and synthesize methane using methanation reactors in the hydrogen-to-methane subsystem. The hydrogen storage works as a buffer connecting the other two subsystems [42]. In addition, hydrogen storage features an immense storage capacity, which is more than 100 times larger than batteries [30].
Ref. [23] first proposes a linear power-to-hydrogen-heat model by portioning the nonlinear feasible region into several triangles. Binary variables are needed to select the triangles.
Ref. [30] then proposes a simplified power-to-hydrogen-heat model using linear interpolation, which reduces the model accuracy. By contrast, we utilize the convex combination methods to convexify the nonlinear feasible region without binary variables while ensuring model accuracy. In addition, the hydrogen-to-methanation process is included in the proposed P2HHM model.

1) UC MODEL FOR THE ELECTROLYSIS PROCESS
In the electrolysis process, the input RE is divided into hydrogen power and heat power. Hydrogen power is utilized to generate hydrogen in electrolyzers. Heat power is utilized to supply heat demand through heat exchangers. Due to the intermittent nature of RE, the on-off status of electrolyzers may vary frequently. In addition, loading regions restrict the operation of electrolyzers. Therefore, it is necessary to consider the UC of electrolyzers. Based on the convex combination method [43], the electrolysis process with UC can be formulated as follows. Since alkaline water electrolyzers have been used in industrial applications [23], this paper selects alkaline water electrolyzers as the P2H units. However, the proposed P2HHM model is not limited to certain types of electrolyzers because we use the feasible region to formulate the electrolysis process, which is represented by P ELZ e,m , H ELZ e,m and T ELZ e,m .
Given these parameters, various types of electrolyzers can be modeled.

2) MULTIPLE HYDROGEN UTILIZATION
The production and utilization of green hydrogen can be modeled as follows.
f ELZ Equations (1.11)-(1.18) formulate the production and utilization of green hydrogen. In the first step, green hydrogen is produced using RE through electrolyzers, modeled with (1.11). In the second step, to ensure the safety of the P2HHM systems and maintain a stable hydrogen supply, the generated green hydrogen is stored in the hydrogen storage (1.12)-(1.14). Then, the stored hydrogen can be utilized in various ways including supply to hydrogen loads, synthesis of methane, and blending with SNG, which are constrained by (1.15). Methane synthesis is done in methanation reactors using hydrogen and carbon dioxide, modeled by (1.16). Finally, a certain amount of hydrogen is mixed with methane and then injected into natural gas systems, represented by (1.17). To satisfy combustion requirements, hydrogen fractions in the SNG should be limited, which is modeled by (1.18).

3) HEAT TRANSFER PROCESS
In the P2HHM model, a water circulation system is added to control the working temperature of electrolyzers. The thermal energy can be injected into district heating systems through heat exchangers. The input heating power, transferred heating power, and heating losses together determine the working temperature of electrolyzers, which can be formulated as follows [23].

IV. GAS FLOW MODELING A. GENERAL GAS FLOW EQUATION WITH HYDROGEN INTEGRATION 1) GENERAL GAS FLOW PARTIAL DIFFERENTIAL EQUATIONS
State equations, continuity equations, and momentum equations are generally enough to formulate the gas flow inside the gas pipelines [44]. The state equations establish the relation between gas pressure, density, temperature, and other gas and pipe-related physical properties, which are shown as follows.
where p and ρ are the pressure and density, respectively. The molecular weight and compressibility factor depend on the mixed gases (e.g., mole fractions and physical properties). The average molecular weight can be calculated as follows [12].
where y i denotes the mole fractions of gas i. In regard to the compressibility factor, it can be either estimated based on the thermodynamics experimental data or calculated based on the following equation [12].
The simplified continuity equations and momentum equations, neglecting the gravity effects and inertia and kinetic energy, can be written as follows [44].
By defining the volumetric flow rate as f = π·D 2 4 · ρ·v ρ 0 and assuming that T and Z remain constant [44], (2.4) and (2.5) can be rewritten as follows.
Equations (2.6)-(2.7) formulate the gas flow considering hydrogen effects. The effects of hydrogen on gas flow dynamics are quantified by M and Z , which are 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. Most of the existing papers do not consider hydrogen effects [44], [45], or just focus on the natural gas systems with steady-state gas flow [12]. However, as demonstrated later, both hydrogen and gas flow dynamics will affect the integrated electricity-gas energy systems (e.g., operation and reliability). Equations (2.6)-(2.7) are a set of partial differential equations, which can be discretized using the implicit finite difference methods [45].

2) DYNAMIC GAS FLOW MODEL
Based on the backward difference methods, (2.6)-(2.7) can be discretized as follows.
Equation (2.9) can further be linearized using the firstorder Taylor series [45], shown as follows.

3) STEADY-STATE GAS FLOW MODEL
In this case, the continuity equations (2.6) are dropped as the inflow equals the outflow. Based on the backward difference methods, (2.7) can be rewritten as follows. In (2.12), the square items at the right side can be replaced by new variables and the left side items can be linearized using piecewise linearization approaches [44].
To summarize, the steady-state gas flow model is formulated by (2.12) and the dynamic gas flow model is formulated by (2.8), (2.10), and (2.11). Based on the above equations, the effects of hydrogen on gas flow dynamics can be quantified.
Remark 4.1: The dynamic gas flow model can also be linearized using piecewise linearization approaches [44]. However, as demonstrated in [57], the computation time for the dynamic gas flow models increases exponentially with the piecewise segments and the system scale. On the other hand, the first-order Taylor series relies on initial operating points, which vary with system states. Therefore, we adopt the piecewise linearization approaches to calculate the initial operating points and then the first-order Taylor series to ensure computational tractability.

V. CASE STUDY
This section verifies the effectiveness of the proposed P2HMM model and hydrogen-integrated gas flow model. Two common problems in energy systems are taken as examples (i.e., optimal operation and reliability assessment with renewable energy). Specifically, for the operation problem, we formulate a Wasserstein metric-based distributionally robust optimal operation model for the integrated electricitygas-heat energy systems. For the reliability assessment problem, we adopt sequential Monte Carlo simulation (SMCS) to evaluate the reliability of the integrated electricity-gas energy systems.

A. APPLICATION I: INTEGRATED ENERGY SYSTEM OPTIMAL OPERATION
Since green hydrogen is produced by renewable energy, renewable uncertainties would affect not only power systems but also hydrogen systems. First, spinning reserves should be optimized to handle renewable uncertainties. In light of the quick response ability of electrolyzers, part of spinning reserves can be provided by the hydrogen systems. Second, the on-off status of ELZs should be optimized to avoid unexpected and frequent alteration, which may harm the safety of hydrogen systems. Therefore, it is necessary to consider the renewable uncertainties when hydrogen is integrated into the IESs. SP and RO are the other two most popular methods to characterize uncertainties. SP relies on a large number of scenarios or precise probability distributions to achieve acceptable accuracy [25]. However, numerous scenarios make the model computationally intractable, and precise probability distributions are difficult to access. Although only limited information is required by RO (e.g., budgets and boundaries of uncertain variables), RO focuses on the worst-case scenario, yielding conservative solutions [35]. By contrast, DRO surpasses SP and RO in computational tractability and solution quality by exploiting historical data. Hence, this paper adopts DRO to consider the renewable uncertainties.
In this example, we propose a Wasserstein metric-based distributionally robust optimal operation (DRO) model for the integrated electricity-gas-heat energy systems (IEGHESs). Due to space limitations, we omit the detailed modeling of distribution systems, natural gas systems, and district heat systems. Readers can refer to [46], [47], and [48] for more details. Here we focus on the modeling of affine rules for tackling renewable uncertainties and the reformulation based on the Wasserstein metric.

1) AFFINE ADJUSTMENT FOR UNCERTAINTIES
At the time of making decisions, the RE output is uncertain to the operator. To handle renewable uncertainties at the real-time stage, the spinning reserve is utilized, which allows power adjustment after the uncertainties are revealed. Affine rules are adopted to allocate the spinning reserve to different units (i.e., conventional generators (GENs), CHPs, and ELZs) and penalties are introduced to punish the reserve violations [49], as follows.
Since we focus on green hydrogen produced using renewable energy, this paper only considers renewable uncertainties, but not load uncertainties. On the other hand, compared with renewable output, loads can be predicted with high accuracy [58]. However, load uncertainties can be considered by redefining the uncertain variable w t as the net load.

2) WASSERSTEIN METRIC-BASED MODEL REFORMULATION
The objective function of the proposed optimal operation model consists of UC costs, operation costs, power adjustment costs, and reserve penalties. The deterministic version of the four items is shown as follows.
The proposed model can be regarded as a two-stage problem. In the day-ahead stage, UC and economic dispatch decisions are made before the revealing of uncertainties. Thus, f UC and f O are deterministic and independent of uncertainties. After the realization of uncertainties, power is adjusted to handle uncertainties in the real-time stage, which is constrained by spinning reserves defined in the day-ahead. Hence, f A and f R are associated with uncertainties w. Let x be the variables independent of uncertainties w and y be the variables dependent of uncertainties w. Correspondingly, let H(x) equal the sum of f UC and f O and Q(y, w) equal the sum of f A and f R . Considering the worst-case realization of uncertain variables, the proposed DRO model is reformulated as follows.
s.t. distribution network constraints (4.6) district heat constraints (4.7) natural gas network constraints (4.8) P2HHM constraints (1) (4.9) Spinning reserve constraints (3.2)-(3.4) (4.10) The above model determines the optimal operating points of various power sources to minimize the sum of day-ahead operation costs and the worst-case expected real-time power adjustment costs plus reserve penalties considering the uncertainties from renewable energy. Note that we do not consider carbon emissions in the current formulation, which can be added to the objective functions or constraints.
Remark 5.1: At distribution levels, the IESs can be operated by one entity, while at transmission levels, different energy systems (i.e., power systems, natural gas systems, and heating systems) usually belong to different entities with various interests, and they only share limited operational information (e.g., the gas requirement of gas-fired generators, the gas demand of power-to-gas units, and the electricity demand of electric boilers) to facilitate the coupling. In this paper, we aim to maximize social welfare and thus follow the widely-adopted assumption that there is a centralized coordinator [11], [44] who operates the integrated energy systems and has access to all the operational information. However, this assumption can be dropped by adopting a distributed framework [25], [48], where only the boundary information is shared.
The worst-case expected cost E B [Q(y, w)] poses challenges in solving the proposed model because of the expectation operator and two different types of objective functions inside. Firstly, to mitigate computational burden, the following approximation can be made [50].
Then, based on the Wasserstein metric, the two items on the right side can be reformulated as two linear programming problems [49]. Interested readers can refer to [49], [51], and [52] for more details. As a result, the original model (4.5)-(4.10) can be reformulated as a mixed-integer linear programming problem that is readily solved using commercial solvers.

3) NUMERICAL ANALYSIS
The proposed optimal operation model is applied to the IEGHESs consisting of the IEEE 33-bus distribution networks, the 20-node natural gas networks, and the 32-node district heat networks, as shown in Fig. 2. In the distribution networks, three CHP units are installed at buses 18, 22, and 33, which are connected to heat source nodes 1, 31, and 32 of the district heat networks, respectively. Three conventional generators are installed at buses 2, 10, and 25, respectively. Generators at buses 2 and 25 are gas-fired generators. Three wind generators equipped with P2HHM systems are placed at buses 1, 13, and 28, respectively. The P2HHM systems are connected to nodes 5, 13, and 19 of the natural gas networks and nodes 1, 31, and 32 of the district heat networks. In the natural gas networks, two gas wells are placed at nodes 1 and 8, respectively. In the district heat networks, three electrical boilers are installed at nodes 1, 31, and 32, respectively. The total capacity of GENs, CHPs,  and renewable generators is 2.4 MW, 2.4 MW, and 7.5 MW, respectively. The total capacity of electrolyzers and electrical boilers is 2 MW and 3 MW, respectively. The peak power, heat, and gas load are set to 3.715 MW, 3 MW, and 0.5 km 3 , respectively. The penalty for the unserved spinning reserve is 1000 $/MWh. The confidence level is set to 0.95. All the detailed parameters used for case studies can be found in the online warehouse [56].
To analyze the effects of hydrogen, the following cases are developed. S1-S3: steady-state gas flow without hydrogen effects. SH1-SH3: steady-state gas flow with hydrogen effects. T1-T3: dynamic gas flow without hydrogen effects. TH1-TH3: dynamic gas flow with hydrogen effects. For each group, hydrogen fractions increase by 10% from 0% up to 20%.
The operation costs of all the cases are compared in Table 1. Regarding the effects of hydrogen fractions, it can be seen that the total operation costs decrease with the increase of hydrogen fractions because of the reduced gas well production. Meanwhile, the operation costs of CHP and ELZs rise first and then decline. On the one hand, the volume of SNG injected into natural gas systems rises with the hydrogen fractions as hydrogen has a lower heating value, which can be seen from the reduction in operation costs of gas wells. Hence, when the hydrogen fraction increases from 0% to 10%, the power consumption of ELZs increases as more hydrogen is produced. CHPs need to increase their output to maintain power balance. Therefore, the operating costs of CHPs and ELZs rise. On the other hand, as the amount of injected SNG is subject to the operational constraints of natural gas systems, the increment of SNG is less than that of hydrogen when the hydrogen fraction increases from 10% to 20%. Hence, the volume of generated hydrogen from ELZs decreases. Consequently, the power consumption of ELZs drops, and so do the CHPs. Hence, the operating costs of CHPs and ELZs fall. Note that the power consumption of ELZs is related to the available renewable energy. Regarding the effects of gas flow models, lower operation costs are observed when the dynamic gas flow models are adopted. This is because the line pack is considered and part of the gas stored in pipes is leveraged to supply gas loads, which decreases the gas well output. In addition, GEN, CHP, and ELZs have smaller operation costs. Regarding the effects of hydrogen physical properties, when the steady-state gas flow is adopted, hydrogen physical properties have no impact on the total operation costs in all the hydrogen fractions. However, when the dynamic gas flow model is adopted, considering hydrogen physical properties raises the operation costs compared with the cases without considering hydrogen physical properties. And the gap becomes larger as the hydrogen fraction increases. It indicates that neglecting the physical properties of hydrogen leads to an incorrect assessment of the scheduling plan.
To illustrate the performance of the proposed DRO approach, three methods are compared (i.e., deterministic approach (DET), SP, and DRO). The in-sample size increases from 5 to 100. Another 5000 samples are used to test the out-of-sample performance. The in-sample analysis can be regarded as the day-ahead stage performance since the uncertainties are not realized at the time of making decisions. In contrast, the out-of-sample analysis can be regarded as the real-time stage performance.
The results are reported in Table 2. It can be seen that in both cases the lowest in-sample cost and highest out-ofsample cost are observed in DET. It indicates that although DET can achieve the least costs in the day-ahead stage, the severest penalties would occur in the real-time stage. As a result, the solutions obtained from DET are not reliable. By contrast, although DRO has the highest in-sample cost, the least out-of-sample cost is achieved in DRO. SP has a moderate performance between DET and DRO. Concretely, the in-sample cost of SP is always lower than that of DRO, while the opposite trend is observed in the out-of-sample scenarios. With the growth in the sample size, the in-sample cost of SP increases while that of DRO decreases. Because more distribution information is revealed, the ambiguity level between the estimated distribution and true distribution is reduced. Meanwhile, the bias in SP is reduced such that the operation cost gradually converges to the true value. Therefore, the gap between SP and DRO shrinks. On the other hand, increasing scenarios contributes to reducing outof-sample costs. Regarding the computational performance, DET takes the least time, followed by SP and then DRO. And the computation time increases with the growth of sample size in SP and DRO. However, even with 100 samples, the computation time does not hike too much. The results indicate the scalability of the proposed approach.
Similar to the results reported in Table 1, considering the physical properties of hydrogen increase the in-sample and out-of-sample costs, regardless of the algorithms.
In a word, the proposed DRO approach has the least outof-sample costs, which can be regarded as real-time operating costs. Both the in-sample costs and out-of-sample costs decrease with the increase in the sample size. In addition, the computation time does not increase significantly with the increase in sample size.

B. APPLICATION II: INTEGRATED ENERGY SYSTEM RELIABILITY ASSESSMENT
In this example, we evaluate the reliability of the integrated electricity-gas energy systems (IEGESs) based on the SMCS approach. State enumeration and non-sequential Monte Carlo simulation are the other two widely adopted methods for the reliability evaluation of composite systems. However, these two methods are not able to consider the temporal feature of renewable energy. Since green hydrogen is produced using renewable energy, neglecting this feature will lead to inaccurate reliability results. Hence, we adopt SMCS to evaluate the reliability of the IEGESs. The two-state outage models are adopted. The component models and optimal load shedding models of transmission systems for reliability assessment are widely studied. Interested readers can refer to [53] for more details. In the rest of this example, we illustrate the reliability assessment indices and procedure and the case studies.

Remark 5.2:
Since the chronological renewable energy, power load, and gas load are considered, other methods (e.g., state enumeration and non-sequential Monte Carlo) can not be applied. In addition, we do not contribute to improving the accuracy or speed of reliability evaluation approaches but use the well-established SMCS approach. Hence, we do not compare various methods in this application as we aim at analyzing the impact of hydrogen and gas flow models on the reliability of the IEGESs with the modeling and formulation proposed in this paper.

1) RELIABILITY ASSESSMENT INDICES AND PROCEDURE
In reliability evaluation, the outage probability and the unserved energy are the two common indicators of interest. Therefore, we develop three groups of reliability indices for the power systems, renewable generators, and natural gas systems, respectively.
The loss of load probability (LOLP) and the expected demand not supplied (EDNS) for the power systems, the loss of renewable energy probability (LORP) and the expected renewable energy not supplied (ERNS) for the renewable generators, and the loss of gas probability (LOGP) and the expected gas not supplied (EGNS) for the natural gas systems are defined as follows [55].
where 1(·) is the indicator function. The coefficient of variation (COV) is used to check the stopping conditions [54]. As the proposed reliability indices cover electricity, renewable energy, and gas, the following criterion is adopted.
where V (·) is the variance.
460 VOLUME 9, 2022 The detailed procedure for the IEGES reliability evaluation is given in Algorithm 1 and the flowchart is shown in Fig. 3.  Solve the optimal load shedding model to determine the amount of load curtailment. 7 end for 8 Calculate reliability indices based on (5.1)-(5.6) 9 Check the stop conditions using (5.7). 10 end for 11 Output the final reliability indices.

2) NUMERICAL ANALYSIS
In this example, the proposed model is applied to the IEGESs consisting of the modified IEEE 24-bus power systems and the 20-node natural gas systems [54], as shown in Fig. 4. The modified IEEE 24-bus power system is comprised of 38 transmission lines, 22 conventional generators, and three wind farms. The three wind farms and the associated powerto-hydrogen-methane (P2HM) systems are installed at buses 21, 22, and 23, respectively. The P2HM systems are connected to nodes 5, 8, and 12 of the natural gas systems, respectively. Generators at buses 1, 13, and 15 are gas-fired and connected to nodes 5, 14, and 2, respectively. The 20-node natural gas system consists of three gas wells, 19 pipes, and three compressors. The total hours of each year, NT, is 8736 h, which is divided into 52 weeks to exploit the large capacity of hydrogen storage. Therefore, the number of time slots, T , in each dispatch block is 168 h. The maximum  number of samples 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 [56].
To verify the merit of the proposed model, the following three cases with 20% hydrogen fractions are developed.
C0: dynamic gas flow with P2HM. C1: steady-state gas flow with P2HM, but no hydrogen effects. This case works as the baseline.
C2: dynamic energy flow without P2HM. Table 3 compares the reliability indices of C0, C1, and C2. It can be seen that the proposed model can improve the reliability of the IEGESs. Specifically, compared with C1, C0 improves LOLP, LOGP, EDNS, and EGNS by 0.53%, 22.27%, 0.87%, and 9.25%, respectively. The natural gas systems receive the biggest boost. Regarding the integration of P2HM, compared with C2, C0 improves the LOLP, LOGP, EDNS, and EGNS by 1.47%, 9.35%, 2.14%, and 4.37%, respectively. The most remarkable improvements are the LORP and ERNS that plummet from 0.1431 and 493.4201 GWh/yr. to 0.0004 and 0.3592 GWh/yr., respectively. The above results indicate that the proposed model can improve the reliability of the IEGESs.
Further, to investigate how the gas flow dynamics and hydrogen physical properties affect the reliability of the IEGESs, the following cases are developed.
C3: dynamic gas flow with hydrogen effects C4: steady-state gas flow with hydrogen effects C5: dynamic gas flow without hydrogen effects When not considering hydrogen effects, the molecular weight and compressor factor are identical to the natural gas. The hydrogen fractions increase by 10% from 0% up to 20% in the above cases.   5. illustrates the effects of gas flow dynamics on the IEGES reliability at various hydrogen fractions. LORP and ERNS are excluded because they are unsusceptible to the gas flow dynamics. It can be seen that considering gas flow dynamics lowers LOLP, EDNS, LOGP, and EGNS, which means the reliability of the IEGESs is improved. The reliability gap between C3 and C4 narrows as the hydrogen fraction rises. In addition, hydrogen fractions have diverse effects on the system reliability in these two cases. In dynamic gas flow models, higher hydrogen fractions lead to lower LOLP and EDNS but larger LOGP and EGNS. The opposite trends in LOLP, LOGP, and EGNS are observed in the steady-state gas flow models. Fig. 6. illustrates the effects of hydrogen physical properties on the IEGES reliability at various hydrogen fractions. Fig. 6 shows that reliability indices are increased after considering hydrogen, indicating hydrogen physical properties have a detrimental influence on reliability. Specifically, considering hydrogen effects on gas flow in C3 exacerbates LOLP, LOGP, EDNS, and EGNS, compared with C5. As the hydrogen fraction grows, LOLP, LOGP, and EGNS increase while EDNS decreases. In contrast, all the reliability indices are improved with the increase of hydrogen fractions when not considering hydrogen physical properties. The difference in reliability between C3 and C5 expands as hydrogen fractions rise.
The improvement in reliability comes from the increased volume of injected SNG and the line pack. The volume of SNG injected into natural gas systems rises with the hydrogen fractions, which can improve the reliability of the natural gas systems. In the dynamic gas flow models, the additional improvement in reliability results from the line pack. Gases stored in the pipes can be used to mitigate load shedding in emergency states. According to the definition of the line pack in (2.11), increasing hydrogen fractions reduces K 1 and thus reduces the volume of the line pack given the same pressure level when considering hydrogen effects. The adverse impact resulting from hydrogen's physical properties outweighs the benefits of increasing hydrogen fractions. Consequently, the improvement in reliability diminishes.

VI. CONCLUSION
This paper investigates the integration of hydrogen into the integrated energy systems, focused on blending green hydrogen into natural gas systems. First, the literature on hydrogen modeling and control, operation, planning, and markets are reviewed. A power-to-hydrogen-heat-methane model with unit commitment is proposed based on convex combination methods. The steady and dynamic gas flow models with the explicit consideration of hydrogen effects are then proposed. Based on the developed P2HHM and gas flow models, a Wasserstein metric-based distributionally robust optimal operation model and a sequential Monte Carlo simulationbased reliability assessment model are proposed to analyze the effects of hydrogen physical properties and various hydrogen fractions on the optimal operation and reliability of the integrated energy systems. Based on the defined case study settings, the following conclusions can be drawn through the numerical simulations: 1) In regard to the optimal operation, increasing hydrogen factions can save the total costs, and considering gas flow dynamics can also reduce the total costs. The total costs decrease from $3527.90 to $2504.08. Hydrogen physical properties have impacts on the dynamic gas flow model. When the dynamic gas flow model is adopted, considering hydrogen physical properties raises the total costs and the gap in the total costs becomes larger as the hydrogen fraction increases. In the case of 20% hydrogen fractions, considering hydrogen effects in the dynamic gas flow model increases the total costs by $24.27, or around 0.97% of the total costs. By contrast, there is no increment in the steadystate model. The proposed Wasserstein metric-based DRO model can achieve lower out-of-sample costs while ensuring computational tractability. In the case of 100 samples, the proposed model yields $3684.63 in-sample cost using 47.59s, and the corresponding out-of-sample cost is $2959. 10.
2) In regard to the reliability evaluation, gas flow dynamics and power-to-hydrogen can enhance the reliability of the IEGESs and minimize renewable energy curtailment. Considering dynamic gas flow reduces LOLP, LOGP, EDNS, and EGNS by 0.53%, 22.27%, 0.87%, and 9.25%, respectively. Considering power-to-hydrogen reduces LORP and ERNS by 99.72% and 99.93%, respectively. However, injecting hydrogen worsens the reliability of the natural gas systems and the adverse influence expands with the increasing hydrogen fractions. When the hydrogen fractions increase to 20%, LOGP and EGNS increase by around 2.74% and 0.58%, respectively. For power systems, the influence on reliability indices is different. In this test system, ENDS is reduced by around 0.07% while LOLP is increased by around 0.01%.