A Unique Electrical Model for the Steady-State Analysis of a Multi-Energy System

: In order to decarbonize the energy sector, the interdependencies between the power and natural gas systems are going to be much stronger in the next period. Thus, it is necessary to have a powerful simulation model that is able to efﬁciently and simultaneously solve all coupled energy carriers in a single simulation environment in only one simulation step. As an answer to the described computational challenges, a unique model for the steady-state analysis of a multi-energy system (MES) using the electrical analogy approach is developed. Detailed electrical equivalent models, developed using the network port theory and the load ﬂow method formulation, of the most important natural gas network elements, as well as of the linking facilities between the power and natural gas systems, are given. The presented models were loaded up into a well-known software for the power system simulation—NEPLAN. In the case studies, the accuracy of the presented models is conﬁrmed by the comparison of the simulation results with the results obtained by SIMONE—a well-known software for natural gas network simulations. Moreover, the applicability of the presented unique model is demonstrated by the MES security of a supply analysis.


Introduction
The main goal of actual energy policies is the decarbonization of the energy sector in the next medium to long term [1][2][3]. According to the joint Ten-Year Network Development Plan 2020 (TYNDP 2020) of ENTSOG and ENTSO-E [4], to comply with the 1.5 • C targets of the Paris Agreement [2], carbon neutrality in the power sector must be achieved by 2040, while in the overall energy sector by 2050. Consequently, the sector coupling between power and natural gas is recognized as the key contributor to achieving the decarbonization targets [5]. Accordingly, renewable energy sources (RES), primarily wind power plants (WPPs) and solar power plants (SPPs), will play an important role in the decarbonization of the power sector. Moreover, the decarbonization of the natural gas system will also have to play a significant role. For that purpose, a relatively new technology called power-togas (P2G) is recognized as one of the most important solutions [4,6]. The role of P2G in the decarbonization process is two-fold. On the one side, P2G technology will convert the surplus generation from intermittent RES into a so-called renewable gas, and help the power sector with variable RES integration [7], while, on the other side, P2G will help the gas sector to achieve the decarbonization targets [8,9]. Moreover, to balance the power system, when additional power injection is needed, e.g., due to a lower RES generation than predicted, the installation of new gas-fired power plants (GFPPs) will take part in the near future [4,10]. Moreover, the natural gas and power sectors are also already linked through the electrically driven gas compressors. In any case, the interdependency between these two energy infrastructures is going to be much stronger in the next period. Thus, the analyses of the coupled and strongly interdependent energy systems, as well as can be easily used in the existing, robust, and well-known electrical software packages. The developed MES model is solved using the extended Newton-Raphson technique that is free of the above-described issues associated with the ordinary Newton-Raphson method. It is worthwhile to mention that the electrical analogy approach allows the setting of only one slack node for the entire MES model.
Finally, this paper is organized as follows: in Section 2, the developed electrical equivalent models of the most important natural gas network elements, as well as the electrical equivalent models of the linking facilities between the power and natural gas systems, are described. Two case studies are conducted in Section 3 to verify and demonstrate the applicability of the presented MES modeling approach, while conclusions are given in Section 4.

The Equivalent Electrical Analogy Models
In order to simulate the natural gas network, it is only important to know gas pressures, volumetric gas flows, and gas temperatures at the element suction (inlet) and discharge (outlet) side (if applicable). In other words, it is not relevant what happens inside the network element in detail, but what goes in/out of it. Hence, in this paper, the so-called network port theory and the load flow method formulation, known in the power systems analysis, were applied to define the electrical equivalent models of the gas network elements.
In the network port theory, the most used were the one-and two-port network models. The two-port network model, shown in Figure 1, has two-terminal pairs: input (port 1) and output (port 2), and it is represented by four external variables: input and output phase voltages (Vel1f and Vel2f, respectively), as well as input and output currents (Iel1 and Iel2, respectively) [34]. It should be noted that subscript el in this paper was used to distinguish electrical and natural gas variables. Furthermore, the two-port network can be treated as a black box modeled by the relationships between the four external variables. Similarly, the one-port network has only an input port (one-terminal pair) and two external variables: Vel1f and Iel1 [35]. The currents flowing into the model were assumed to be positive, and the voltages had positive reference polarities [36].
In this paper, we differentiated two types of variables: internal and external. The internal variables referred to the element ports themselves, and the mandatory ones were, according to the load flow method [37], voltage magnitude Vel and angle Ael (in this case, the line voltage took place, so the index f was omitted), as well as active Pel and reactive Qel power. Additional variables may also be specified, such as the number of transformer taps, while some of them, such as the current flow, were derived. The external variables can refer to any other element in the model. Only the internal variables were mandatory. It should be noted that, as it has a strong influence on the form of equations, the model's load flow equations together with their first partial derivatives by internal and external variables (elements of the Jacobian matrix) had to be specified in the load flow calculation process. Furthermore, it is common in the load flow formulation that the parameter values used in the load flow equations are expressed in the p.u. system (marked in bold) with In this paper, we differentiated two types of variables: internal and external. The internal variables referred to the element ports themselves, and the mandatory ones were, according to the load flow method [37], voltage magnitude V el and angle A el (in this case, the line voltage took place, so the index f was omitted), as well as active P el and reactive Q el power. Additional variables may also be specified, such as the number of transformer taps, while some of them, such as the current flow, were derived. The external variables can refer to any other element in the model. Only the internal variables were mandatory. It should be noted that, as it has a strong influence on the form of equations, the model's load flow equations together with their first partial derivatives by internal and external variables (elements of the Jacobian matrix) had to be specified in the load flow calculation process. Furthermore, it is common in the load flow formulation that the parameter values used in the load flow equations are expressed in the p.u. system (marked in bold) with the base power of 100 MVA. Moreover, these equations were expressed in an implicit form where the right part (=0) was omitted. Since the reactive power does not exist in natural gas networks, the active power was equal to the apparent power, i.e., the power factor was equal to 1. Consequently, the three-phase active power can be calculated as: where V el represents the line voltage.
It should be noted that in the p.u. system of the balanced power network, the phase and line voltages were equal [38,39].
The volumetric gas flow Q can be calculated as [32]: where z, R u , T, M w , p, and . m are the compressibility factor, the universal gas constant, the absolute gas temperature, the molecular weight of the natural gas composition, the gas pressure, and the gas mass flow, respectively.
The compressibility factor z can be calculated using several well-known equations. In this paper, Papay's equation [40] was used: where p r and T r are the reduced gas pressure and temperature, respectively. Since the reduced gas pressure p r was determined as a relation of the actual pressure to the critical value, to make the power flow equations shorter and easier to read, (3) was segmented into two constant parts. These parts were further, in the load flow equations, multiplied by the gas pressures (port voltages). Thus, for the two-port network element, the constant parts were: where k 11 and k 21 , as well as k 12 and k 22 , are the first and the second constant parts of the first and second port compressibility factor equation, respectively. k pr and k Tr are constants calculated as the sum of all relations of gas portions in the natural gas mixture to its critical values (gas pressures and temperatures, respectively). T 1 and T 2 represent gas temperatures at both element sides. In the following subsections, the load flow equations of the most important natural gas network elements, as well as P2G, and GFPP, are described.

Gas Compressor Station
The gas compressor station in this paper was modeled as the two-port network element. A detailed electrical analogy model of the gas compressor station can be found in [32].
In the load flow formulation, the gas compressor station was to be considered as the PV type, i.e., the active power and the voltage had to be specified/known. Moreover, as the reactive power does not exist in the natural gas networks, the reactive power had to be set to zero.
The next set of equations depended on the compressor's operation state, and they were based on the gas mass flow conservation through the compressor as well as on the gas pressures relation according to the compression ratio [32].
If the gas compressor was in operation, then the inlet gas mass flow minus the gas driver consumption had to be equal to the outlet gas mass flow, i.e.,: where k 0 represents the gas mass flow consumption by the driver, expressed in % of the inlet gas mass flow. By applying the electrical analogy approach to (8) and the described segmentation techniques to (3), and expressing the volumetric gas flow (i.e., electric current) using (1), the gas mass flow through the compressor station in p.u., i.e., the power equation, can be expressed as follows: where kVb 1 and kVb 2 are inlet and outlet base pressures (voltages) in bar that should be applied to convert p.u. voltages into natural units as demanded by (3). Moreover, the gas mass flow at the discharge side was flowing out of the two-port network model, so the positive sign came in front of the P el1 part in (9).
If the compressor driver was electric-powered, then k 0 was equal to zero. The discharge gas pressure, i.e., the voltage equation, depended on the compression ratio applied. In this paper, three types of compression ratios could be set. The first one was the user-specified fix compression ratio r c , so the compressor station inlet and outlet gas pressures (voltages) related as: The second compression ratio option was to set the outlet pressure to the nominal value: The third compression ratio option was to copy the p.u. pressure from the inlet to the outlet side: If the gas compressor was bypassed, then the inlet and outlet sides had the same gas pressures, temperatures, and compressibility factors, while the driver consumption k 0 was equal to zero. The application of the former to (9) resulted in the inlet and outlet active powers equality (with opposite signs since discharge flow was out of the two-port network model): and in gas pressure equality (p.u. pressures had to be converted to the natural units by multiplying with their base values): A gas compressor was also bypassed if it was electric-powered and the electric voltage of the bus from which the gas compressor driver was supplied (external variable) became an outage.

Compressor's Electric Driver Consumption
If the gas compressor was electric-powered, then it was necessary to connect an appropriate load to the supplying node in the power network. In this paper, the electric load of the compressor driver was modeled as the one-port network element and as the PQ type in the load flow formulation (active and reactive powers had to be specified/known). The electric load of the compressor driver needed to track the compressor operation state. Hence, the external variables needed to be defined as follows: • V el-ext1 -compressor's voltage at the suction side; • P el-ext1 -compressor's active power at the suction side; • P el-ext2 -compressor's active power at the discharge side.
In this paper, if the gas compressor was in operation, the reactive power was set to the fixed value: where kq and kVAb are the reactive power constant settings by the user in kvar and base power of 10 5 kVA, respectively. Moreover, the driver's active power consumption could be determined in two ways. The first one was to calculate the active power consumption as a percentage of the gas mass flow through the compressor. Accordingly, to determine the power consumption, gas energy had to be converted into kWh using the gas heating value. In some countries, a low heating value (LHV) is used, while in others, a high heating value (HHV) takes place. So, in this paper, the heating value was defined as a constant that could be set by the user. The volumetric gas flow at the compressor's suction side could be determined using the defined external variables V el-ext1 , P el-ext1 , and (1). Since the gas heating value was given per standard gas volume, the calculated actual volumetric gas flow in A, i.e., in m 3 /s, needed to be converted to the standard gas flow (gas flow at standard conditions: 15 • C and 101.325 kPa, in Nm 3 /s) as follows: where subscripts st and 1 denote standard conditions and the compressor's station suction side, respectively. The suction gas pressure p 1 was represented with V el-ext1 , and the compressibility factor z 1 was calculated using the segmentation of (3), i.e., multiplying (4) and (5) with V el-ext1 . Finally, the compressor driver power consumption P el1 was defined as: (1−k11·Vel-ext1·kVb+k21·V 2 el-ext1 ·kVb 2 ) It should be noted that the gas pressure p 1 in (16) was in Pa, while in the load flow equation the p.u. voltage multiplied by the base value gave pressure in bar. Because of that, factor 10 5 exists in (17). Moreover, factor 3600 in (17) was a result of the conversion of kWh/h into kWh/s. The second possibility was to set the active power manually to the fixed value: where kp is the active power constant setting by the user in kW. If the gas compressor station inlet and outlet active powers in the unique electrical MES model (P el-ext1 and P el-ext2 ) were equal, then the compressor was bypassed, and, consequently, the driver's electric power consumption (P el1 and Q el1 ) equaled zero.

Gas Pressure Reducing and Metering Station
The gas pressure reducing and metering station (GPRMS) in this paper was also modeled (such as the gas compressor station) as the two-port network element and the PV type in the load flow formulation. A detailed GPRMS electrical analogy model can be found in [41].
As for the gas compressor, the GPRMS reactive power also had to be set to zero. Regarding the gas mass flow balance through the GPRMS, two possibilities exist: GPRMS with or without gas preheating prior to entering the control valve [41]. In the first case, when the preheating process took place, the gas mass flow in p.u., i.e., the power equation, through GPRMS was similar as for the gas compressor, but without the gas compressor driver's consumption, i.e.,: In the second case, when the gas was not preheated prior to entering the control valve, the Joule-Thomson effect had to be taken into account. Accordingly, the discharge gas temperature can be calculated as follows [41]: (20) where JT is the Joule-Thomson coefficient in • F/100 psi. In this paper, JT of 7 • F/100 psi, as suggested by the American Gas Association (AGA), was used, but the user could set any value to the JT variable in the model.
Since the discharge gas temperature Equation (20) contains suction and discharge gas pressures, the discharge compressibility factor equation had to be segmented into six constant terms k 1v2 -k 6v2 : Applying (20)-(26) into (19), the gas mass balance through GPRMS without a preheating process resulted in: In both cases, the discharge pressure was set to the nominal value, so (11) could be applied as the voltage equation.
If the suction gas pressure was lower than the nominal discharge value, then the GPRMS was bypassed, so (13) and (14) could be applied.

Gas Pipeline
The electrical analogy model of the gas pipeline in a steady state is given in [33], and it is represented by a serial connection of electrical resistance and the constant electric potential source. The former represents the resistance to the gas flow through the pipeline, while the latter imports the effect of the pipeline inclination, i.e., potential energy into the electrical analogy. If the pipeline is horizontal, the latter part of the model does not exist.
However, in [33], the standard volumetric gas flow is analog to the electric current, while in this paper, the actual volumetric gas flow took place since it was directly proportional to the gas velocity upon which the gas flow was constrained through the pipeline. Thus, it was necessary to derive the resistance equation in terms of the actual gas flow and p.u. values. The gas pressure drop along the pipeline can be calculated as [33]: where E p and R p are the constant electrical potential source, i.e., potential energy due to the pipeline inclination, and pipeline resistance, respectively. The potential energy E p can be determined as follows [33]: where H 1 and H 2 are heights of the pipe upstream and downstream side, respectively. Subscript avg represents the average value between the upstream and downstream variables. Moreover, in [33], several equations for the gas pipeline resistance in a fully turbulent flow are given, according to the Darcy friction factor expression applied. In this paper, the AGA fully turbulent gas pipeline resistance equation was used, as a tradeoff between the accuracy of the calculation results and the simplicity of expression: where L, G, η p , D, and ε are the pipeline length, gas specific gravity, pipeline efficiency factor, pipeline diameter, and surface roughness, respectively.
The gas specific gravity G can be substituted with molecular weight M w as follows [33]: The pipeline efficiency factor η p represents additional friction imposed by fittings (e.g., valves, bends, tees, etc.), corrosion, dust deposition, etc. [33].
To convert the standard to actual volumetric gas flow (28), (30) and (16) should be applied.
The gas pipeline in this paper was, such as the compressor station and GPRMS, modeled as the two-port network element and the PV-type in the load flow formulation. The reactive power of both pipeline sides also has to be set to zero.
After applying the segmentation to the compressibility factors using (4)-(7) and the substitution of the current with active power and voltage using (1) in (30), and applying (30) to (28), the pressure equation is as follows: where R p const and E p const represent constants calculated as follows: Moreover, the absolute value of P el1 in (32) was needed as the gas pipeline model had to be reciprocal, i.e., upstream and downstream sides could be freely switched in the Energies 2021, 14, 5753 9 of 24 network model. The same situation could also be observed when the gas flow changed its direction through the pipeline.
The gas mass flow had to be equal at both pipeline sides, i.e., the mass flow equality must be satisfied. Hence, applying the segmentation of z and (1) into (2), the power equation is as follows: 2.5. Gas Load, Gas Input, Underground Gas Storage In this paper, gas load, gas input, and underground gas storage were modeled as the one-port network elements and the PQ type in the load flow formulation.
The below-derived load flow models were similar for all three gas network elements. The difference was only in the sign of the desired volumetric gas flow that is usually defined at standard conditions Q st and given in 1000·Nm 3 /h. For the gas load, a positive value was used, for the gas input a negative value was used, while for the underground gas storage both signs were possible, depending on the operation state (positive for the gas injection process and negative for the gas withdrawal process).
Having defined Q st , the gas mass flow . m in kg/s can be calculated as follows [32]: As in all the natural gas network elements described above, the reactive power had to be set to zero.

Power-to-Gas
The power-to-gas facility in this paper was also modeled as the two-port network element and the PQ type in the load flow formulation. The first port was dedicated to the power network side, while the second port was connected to the natural gas network side. The input values were power consumption kp in MW, the overall process efficiency η P2G in %, as well as the temperature of the produced gas T 2 in K.
Since the P2G facility consumed only active power and bearing in mind that reactive power does not exist in natural gas networks, the reactive power of both ports had to be set to zero.
In the power equation of the first port, the active power was set to the defined fixed value: In the power equation of the second port, the active power was calculated using (1), (3), (6), (7) and (16) as follows:

Gas-Fired Power Plant
As power-to-gas, the gas-fired power plant (GFPP) was also modeled as the two-port network element where the first node was connected to the power network, while the second node was connected to the natural gas network.
The gas temperature T 2 in K, the overall power plant efficiency η GFPP in %, and the power generation kp in MW represented the input values.
Furthermore, there were two possibilities in the GFPP operation regime. In the first one, the voltage of the power bus at which the GFPP was connected could be regulated, while in the second one, a fixed value of the reactive power generated kq in Mvar could be set. In the load flow formulation, the former was represented as the PV type, while the latter was represented as the PQ type.
Moreover, the GFPP model defined the minimum allowed gas pressure of the gas node at which GFPP was connected. If the threshold setting was violated, the GFPP was disconnected, so the active and reactive powers of both GFPP's sides had to be set to zero.
If the GFPP regulated the voltage of the connecting bus to the desired value kv expressed in % of the nominal value, the voltage equation was applied: Moreover, if the fix reactive power generation kq in Mvar was set, the reactive power equation was executed: Since the reactive power does not exist in the natural gas network, it had to be set to zero at the second port.
The active power equation of the first port was similar as for power-to-gas, but with the opposite sign: Finally, the active power equation of the second port was the same as (39), but the variable η GFPP instead of η P2G should be used.

Case Studies
In this section, two case studies were conducted to verify the presented equivalent electrical analogy models and applied modeling approach described in the previous section.
The presented equivalent electrical analogy models were written in the C/C++ scripts and compiled to the Dynamic Link Library (DLL) files, which could be further assigned in the NEPLAN (a well-known software package for power network simulations) graphical editor as the User-Defined Models (UDM) [42]. The C/C++ user's code must contain several routines and follow certain instructions that are described in the NEPLAN manual [42].
In the first case study, an extensive gas network was modeled in the NEPLAN electric module, using the presented equivalent electrical models of the gas network elements, as well as in the well-known commercial software for the natural gas network analyses-SIMONE. The verification of the results obtained by the NEPLAN was conducted by a comparison with the results obtained by SIMONE.
After the verification of the accuracy of the applied methodology in the first case study, in the second case study, an electrical power network was added to the same model together with the main coupling facilities between the power and gas networks: the GFPP, the P2G, and the electric-driven natural gas compressor station. The goal of the second case study was to demonstrate the applicability of the presented modeling approach as well as to emphasize the importance of having the unique MES simulation model.
The same gas mixture from Table 1 was used in all simulations within both case studies. Additionally, in all calculations, the gas heating value of 9.260625 kWh/Nm 3 was used. Reference conditions were 15 • C and 101.325 kPa.

Case Study 1
The gas network used to verify the applied approach represented an imaginary network consisting of 35 nodes, 27 pipelines, 4 gas compressor stations (CS), 6 GPRMSs, 7 gas exits/loads, 2 gas inputs, and 1 underground gas storage (UGS). It was designed to cover a wide range of input data: gas inputs and loads, pipeline diameters, roughness, lengths, and elevations. The model of the gas network, constructed in SIMONE, is shown  Figure 2. Two valves VA1 and VA2 served to switch off the pipelines N2-N7 and N7-N8 out of operation in situations when the gas flow velocities through them fell below the critical value. In those situations, all downstream natural gas flows through the pipeline N2-N8. The gas network used to verify the applied approach represented an imaginary network consisting of 35 nodes, 27 pipelines, 4 gas compressor stations (CS), 6 GPRMSs, 7 gas exits/loads, 2 gas inputs, and 1 underground gas storage (UGS). It was designed to cover a wide range of input data: gas inputs and loads, pipeline diameters, roughness, lengths, and elevations. The model of the gas network, constructed in SIMONE, is shown in Figure 2. Two valves VA1 and VA2 served to switch off the pipelines N2-N7 and N7-N8 out of operation in situations when the gas flow velocities through them fell below the critical value. In those situations, all downstream natural gas flows through the pipeline N2-N8.  Within case study one, six simulations were performed, as indicated in Table 2. The same gas pipeline parameters, presented in Table 3, were used in all defined scenarios. The pipeline efficiency was taken as 1, and the total length of the gas network was 681.5 km. The simulation parameters, shown in Table 4, differ for the UGS injection and the withdrawal process. When the UGS is in the gas injection phase, the gas reduction station GPRMS5 and the compressor station CS4 were in operation. Otherwise, when natural gas was withdrawn from the UGS, the GPRMS6 and the CS3 were in operation. During the withdrawing process, the second gas source (i.e., connection to another gas network) INPUT2 changed the role and became a gas consumer. For all compressors, the driver's gas consumptions were neglected; therefore, in this case, they could be considered as electric-driven compressors. The main goal was to determine node gas pressures and branch volumetric gas flows. The deviations between the obtained results were calculated as: where NEPLAN and SIMONE denote results obtained by models developed in NEPLAN and SIMONE, respectively. The results of the node gas pressures and branch volumetric gas flows are shown in Tables 5 and 6, respectively. It should be noted that for each node, Kirchhoff's current law had to be satisfied. This meant that the sum of all volumetric gas flows going into a node had to be equal to the sum of all volumetric gas flows going out of the same node. Since the gas temperature was taken as a constant along the pipeline and the pipe cross-section remained unchanged, according to (2) the volumetric gas flow at the pipe exit was greater than at the pipe inlet side, as a result of the gas pressure drop. For this reason, both the pipe inlet and the outlet volumetric gas flows were compared to the SIMONE results. If some node only had two connected branches, then the volumetric gas flows of both branches at the sides connected to that node were the same. In that case, only one value was presented in Table 6. However, if some node had three or more connected branches, then Table 6 contains the volumetric gas flows of all connected branches. In those situations, the name of the second node of the corresponding branch was specified in parentheses in order to provide information about the volumetric gas flow branch affiliation.  Due to different equations used for the calculation of the gas pressure drop along the gas pipeline in SIMONE and the proposed equivalent electrical model in NEPLAN, it was evident that deviations between the obtained results would exist. The node gas pressure deviations, as well as the branch volumetric gas flow deviations, were calculated and presented in Figure 3. In normal situations (scenarios S1 and S4), the absolute node gas pressure deviations were up to 2.49% (standard deviation 0.56%), while the absolute branch volumetric gas flow deviations were up to 3.50% (standard deviation 0.82%). Larger deviations were recorded in the N-1 situations when CS1 or CS2 were out of operation, i.e., bypassed. However, the standard deviations were slightly higher and were still very acceptable (0.81% and 1.18% for the node gas pressure and branch volumetric gas flow deviations, respectively). Hence, it can be concluded that the developed equivalent electrical models could efficiently simulate natural gas network steady-state load flows. bypassed. However, the standard deviations were slightly higher and were still very acceptable (0.81% and 1.18% for the node gas pressure and branch volumetric gas flow deviations, respectively). Hence, it can be concluded that the developed equivalent electrical models could efficiently simulate natural gas network steady-state load flows.

Case Study 2
In this case study, the electrical network, together with the coupling facilities, was added, so a unique electrical model of the MES, as shown in Figure 4, was constructed. The power network consisted of 14 nodes, 13 loads, 12 lines, 2 GFPPs, 1 WPP, and 2 transformers that couples 2 voltage levels of 110 kV and 20 kV. The coupling facilities between the gas and the power network were two GFPPs, the P2G and the electric driven gas compressor CS2, marked with bold lines in Figure 4. It should be noted that the GFPP1el and the GFPP2el were numerically equal to the GFPP1 and the GFPP2, respectively, and should have been used only when any of the networks were going to be analyzed separately, such as in case study one. In those situations, GFPP1 and GFPP2 should be disconnected. The frequency of the power system was taken as 50 Hz.

Case Study 2
In this case study, the electrical network, together with the coupling facilities, was added, so a unique electrical model of the MES, as shown in Figure 4, was constructed. The power network consisted of 14 nodes, 13 loads, 12 lines, 2 GFPPs, 1 WPP, and 2 transformers that couples 2 voltage levels of 110 kV and 20 kV. The coupling facilities between the gas and the power network were two GFPPs, the P2G and the electric driven gas compressor CS2, marked with bold lines in Figure 4. It should be noted that the GFPP1el and the GFPP2el were numerically equal to the GFPP1 and the GFPP2, respectively, and should have been used only when any of the networks were going to be analyzed separately, such as in case study one. In those situations, GFPP1 and GFPP2 should be disconnected. The frequency of the power system was taken as 50 Hz.
The parameters of electric lines and loads were given in Tables 7 and 8, respectively. The total length of the power network was 285 km, and all lines were overhead types.
Each of the two identical 110/20 kV transformers, with parameters presented in Table 9, on the primary side had an on-load tap changer with a total of 21 steps of 1.5% of nominal voltage. The secondary voltage was regulated to 103% of the nominal voltage. Moreover, one of the advantages of a unique electrical MES model was the opportunity to have only one slack node for both networks. For this reason, the third (imaginary) transformer, called TR-slack, was added into the unique model to couple busbars INPUT1 and B1 having the network feeders (slack nodes) connected. Parameters for the TR-slack are also given in Table 9, and it had the same on-load tap changer as TRF1 and TRF2.  The parameters of electric lines and loads were given in Tables 7 and 8, respectively. The total length of the power network was 285 km, and all lines were overhead types.   The WPP park consists of 25 units of 5 MW, so the total generation was set to 125 MW and 0 Mvar.
The GFPP1 has an efficiency of 60%, and the operating point of active power was set to 10 MW, while the reactive power controlled the voltage of node B7 to 1.03 p.u.. The minimum allowed gas pressure was set to 85% of the nominal gas pressure at node EXIT1.
The GFPP2 also has an efficiency of 60%, and the minimum allowed gas pressure was 85% of nominal gas pressure at node EXIT6. The operating point was set to 70 MW and 2 Mvar.
Dissimilar to case study one, in case study two, the compressor driver consumption was set to 1.3% of the inlet gas mass flow to all compressors. For compressor CS2, which was the only electric-driven compressor, in this case, the driver consumption of 1.3% of the inlet gas mass flow was recalculated to the electric side via the above-described gas heating value at reference conditions. Pipe efficiency was not changed in this case, i.e., it was left as 100% (no additional gas pressure loss along the gas pipeline was introduced).
Within case study two, three scenarios (S7, S8, and S9) were simulated as described below. All simulations were performed for both UGS operation states: the injection and the withdrawal process.
The goal of scenario S7 was to identify the influence of the operation state of GFPP1, GFPP2, and P2G on the natural gas network node pressures and volumetric gas flows. The performed simulations showed that the same results were obtained for both the UGS operation processes, and they can be summarized as follows: • When GFPP1 and GFPP2 were in operation, the most influenced gas nodes were N6 and N25, respectively, (gas pressures were decreased for 0.01 p.u. and 0.088 p.u., respectively). At the same time, the volumetric gas flows at nodes EXIT1 and EXIT6 were increased by 1.1% and 16.4%, respectively, both compared with situations when GFPP1 and GFPP2 were out of operation. The decreased gas pressures of N6 and N25 were still higher than the nominal gas pressures of nodes EXIT1 and EXIT6; therefore, the influences of GFPP1 and GFPP2 on the gas network could be neglected.

•
When the P2G was in operation, the most influenced gas node was N12 (gas pressure was decreased by 0.014 p.u.). Interestingly, the volumetric gas flow produced by P2G (0.026 m 3 /s) was a little higher than the volumetric gas flow at node N14 demanded by EXIT3 (0.022 m 3 /s). This meant that the losses in the natural gas network could be reduced since natural gas did not need to flow anymore from node N9 to node N14 (130 km). Its flow reached velocities below the minimum allowed (less than 1.5 m/s). Thus, it can be concluded that P2G had a positive influence on the gas network.
In scenario S8, the influences of the CS1 operation state on the security of the power network supply were analyzed. In a normal operation, the CS1 regulates the output pressure to the nominal value. In the UGS injection process, when the CS1 was bypassed, the gas pressure of node N6 fell to 29.3 bar, as shown in Figure 5. Since the nominal pressure of node EXIT1 was 35 bar, GPRMS1 was bypassed, so gas pressures at nodes N6 and EXIT1 were equal. The gas pressure of 29.3 bar at node EXIT1 was lower than the minimum gas pressure required by GFPP1; thus, GFPP1 fell out of operation. This was reflected in the power network in such a way that lines B12-B13 and B9-B10 were overloaded (112.5% and 135.9%, respectively). The overloading of overhead lines in contingency situations is usually allowed up to 120%, so the overloading of line B12-B13 was not critical. The overloading of line B9-B10 was outside the upper limit; therefore, the load curtailment was expected in the 20 kV power network. This situation could be overcome if the UGS was in the withdrawal process. In that case, the gas pressure at node EXIT1 would be 0.89 p.u. and GPRMS1 would not fall out of operation.   In scenario S9, the influences of the CS2 operation state on the security of the power network supply, as well as the influences of a contingency in the power network on the CS2 driver power supply, were analyzed. The latter represents the influence of the power network on the security of the natural gas network supply. In a normal operation, CS2 regulates the output pressure to the nominal value. If CS2 was bypassed, the gas pressure at node N23 would become equal to the gas pressure at node N22, as shown in Figure 6. Consequently, in both UGS operation processes, the gas pressure at node EXIT6 would fall below the minimum required by GFPP2. Thus, GFPP2 would fall out of operation, which would cause voltage drops at B5, B6, B12, and B14 below the minimum allowed of 0.9 p.u., as well as the overloading of lines B1-B3, B3-B6, B9-B11, and B9-B10. Only the latter had an overload above the allowed 120%. Consequently, load curtailment in the power network could be expected to overcome this issue.
The same situation, as described in scenario S9, could occur if line B1-B2 fell out of operation. As a result, the CS2 driver would stay without a power supply, so the CS2 would have to be bypassed. Thus, the operation problems in the power network could be reflected in the natural gas network and returned to the power network, worsening the situation in the power network. Without the implementation of an integrated approach to modeling multi-energy systems, the problem mentioned above would not be observed. The proposed unique electrical model of a multi-energy system that uses the electrical analogy approach easily copes with this problem.

Conclusions
In this paper, the steady-state electrical equivalent models of the most important gas network elements (gas compressor station, gas pressure reducing and metering station, pipeline, gas load, gas source, and underground gas storage), as well as the steady-state electrical equivalent models of the linking facilities between the power and natural gas systems (gas-fired power plant and power-to-gas facility) were given. The equivalent models were developed using the network port theory and the load flow method formulation, known in the power systems analysis. In that way, the overall multi-energy system could be simultaneously solved as one electrical network and in one simulation step, i.e., without high resource-consuming iterative inter-model calculations associated with the separate network modeling approach. Moreover, the use of the extended Newton-Raphson method exceeded the calculation prerequisites associated with the ordinary Newton-Raphson method.
Two case studies were conducted. The developed models of network elements were used for modeling MES in the electrical software package-NEPLAN. In the first case study, the high accuracy of the presented models was confirmed by comparing the simulation results of one gas network modeled in the electrical analogy with the results obtained by SIMONE (a well-known software package for natural gas network simulations). The applicability of the presented approach was demonstrated in the second case study, where the electrical network was added to the gas network model together with the linking facilities. Finally, the case studies showed that the presented approach could be easily applied to the security of supply analyses (scenarios S8 and S9), as well as to determine the interdependencies between the coupled network infrastructures (scenario S7). Moreover, by conducting analyses on a multi-energy system, it was possible to detect some system vulnerabilities that could not be observed if the networks were modeled separately. For example, scenario S8 showed that operation problems in the electrical power network caused by outages in the natural gas network could be avoided if the operating regime of the underground gas storage was changed. Additionally, scenario S9 showed that the outages in the power network could affect the natural gas network, causing problems there, and, consequently, go back and worsen the situation in the power network.
The application of the presented methodology was limited to steady-state networks. Thus, to be able to use the model on an hourly basis, it is necessary to take into account natural gas stored in the gas pipelines, so-called linepack. The development of such a quasi-stationary model of a multi-energy system is the direction of our future research.