Implementation of Unbalanced Thermoelectric Equivalent Circuit for Power Flow and Thermal Rating of Underground LV and MV Cables

This manuscript proposes a time-series temperature-dependent power flow method for unbalanced distribution networks consisting of underground cables. A thermal circuit model for unbalanced three-phase multi-core cables is developed to estimate the conductor temperature and resistance of Medium (MV) and Low Voltage (LV) distribution networks. More specifically, a novel approach is proposed to model and estimate the parameters of the three-phase thermal circuit of 3- or 4-core cables, using the results of Finite Element Method and Particle Swarm Optimization. The proposed approach is generic and can be accurately applied to any kind of 3- or 4-core cables buried in homogeneous or non-homogeneous soil. Furthermore, it is applicable in cases where one or more adjacent cables exist. Using the proposed approach, the conductor temperature of each core can be individually and precisely calculated even under unbalanced loading conditions. The proposed approach is expected to be an important tool for simulating unbalanced distribution networks and estimating the conductor temperatures. The proposed thermal circuit is validated using two 4-core LV and one 3-core MV cables buried in different depths in homogeneous or non-homogeneous soil. Time-series power flow for a whole year is performed in a 25-bus unbalanced LV network consisting of multicore underground cables.


I. INTRODUCTION
P OWER flow is a fundamental tool for several power system applications such as planning, state estimation, real time control and optimization, contingency and stability analysis, etc. Conventional power flow is based on the assumption of a fixed conductor resistance, i.e., the line resistances have pre-specified and invariable values, without considering their dependency on the conductor temperature. Therefore, conventional power flow approaches present inaccuracies, due to the neglect of continuous variation of conductor resistances.
Nowadays, a large portion of distribution networks consists of underground Medium Voltage (MV) and Low Voltage (LV) cables. At present, 55% of all 1 kV distribution lines and 41% of all distribution lines between 1 kV to 100 kV, in Europe, are underground [1]. Looking into the future, more and more EU member states are moving towards underground electricity distribution networks [1]. In addition, a high uptake of Distributed Generators (DGs) and Electric Vehicles (EVs) is also observed, which is possibly going to force distribution networks to their limits [2]. As a result, accurate power flow modeling including the thermal characteristics of the unbalanced underground LV and MV distribution grid is necessary for the planning, analysis, operation, control, and supervision of these networks.

A. Challenges of Modeling Underground MV/LV Cables
Overhead lines present low thermal inertia reaching their steady-state temperature in a few minutes [3], [4]. As a result, the conductor temperature can be estimated based on the conductor current and weather conditions only at the present time instant. Several steady-state temperature-dependent (TD) power flow approaches have been proposed in the literature, so far, for networks with bare overhead conductors. All of them estimate the conductor resistance based on the weather conditions and line loading only at the present time instant, ignoring the line temperatures at the previous time-instants [3], [5]- [9].
On the contrary, underground cables have a much higher thermal inertia, due to the large specific heat constants of the insulation, shield, screen, sheath, and jacket of the cable as well as the soil [4]. As a result, the historical thermal condition of the cable and soil must be considered for accurately estimating the dynamic conductor temperature, increasing the complexity of the calculation process.
Another challenge is the modeling of non-homogeneous soils, which is usually the case in installations where a backfill material exists with different thermal resistivity and specific heat constant than those of native soil [10]. IEC and IEEE standards provide analytical formulas for the modelling of soil, but they are only applicable to a limited number of installation geometries [11], [12]. In practice, several layers of different thermal resistivities and specific heat constants may be present between the cable and the ground surface. In such cases, analytical methods cannot be applied [10].
The complexity of modeling underground MV/LV cables is even further increased in cases where multiple cables are buried next to the examined cable, due to the mutual heating between the cables [13]. More specifically, the heat produced by a buried cable can drastically affect the temperature of an adjacent cable.
Moreover, MV and LV cables are often 3-core or 4-core cables, which are asymmetrically loaded in unbalanced distribution networks. The unbalances are expected to increase in the near future due to the connection of single phase DGs and EVs [14]. As a result, the cores of multicore cables exhibit different temperatures, which should be individually calculated, considering the load unbalances and not simplistically assuming that the three phases are equally loaded, as typically occurs in the existing literature so far [11], [12].

B. Literature Review
Temperature-dependent power flow, known also as Elec-troThermal Coordination (ETC), has been introduced by several research groups to more accurately model distribution and transmission networks, and exploit their maximum power transfer capacity by considering the thermal inertia of conductors and cables.
ETC for cable based transmission grids was introduced in [15], where the authors divided the cable and its environment into a large number of zones, and expressed them into a Thermoelectric Equivalent Circuit (TEE). The TEE interacts with a power flow tool in such a way that the current sources of TEE are updated based on the results of power flow, while the parameters of power flow, i.e., the cable resistances, are updated according to the results of TEE. The authors in [16] modelled underground cables, when real-time measurements of the temperature at the surface of the cable are available. The TEE of [15], [16] was improved by the authors in [17], by reducing the number of thermal zones of the soil, and therefore, the number of components of TEE, without sacrificing the accuracy of the model. In [13], the mutual heating effects caused by neighboring cables were introduced into the TEE, increasing its accuracy in case of multiple buried cables.
Authors in [18] implement Optimal Power Flow (OPF) to minimize the DG curtailment due to the thermal violation of underground cables, transformer and overhead conductors of the network. The power flow is combined with a TEE in a similar sense as in [15] to couple the current calculated from power flow with the conductor temperature calculated by TEE. The TEE is constructed using the approach of [19], in which the cable and soil were divided into a total of 4 zones using formulas from IEC standard [11] and the Van Wormer coefficients [20].
Authors in [4] propose a single-phase time-series Newton-Raphson (NR) power flow method for transmission networks  considering the TEE and the different thermal inertia constants  of transformer, overhead lines and underground cables. Finally, authors in [21] propose an approach to calculate the thermal resistance of the filling layer and inner sheath of three core cables under balanced conditions, using the results of a FEM software. However, they ignore several thermal parameters that play an important role under unbalanced loading conditions, such as the mutual thermal resistances and capacitances between the cores (R m , C m in Fig. 2). Furthermore, the calculation of thermal capacitances between the cores and the surface of the cable (C s1 , C s2 in Fig. 2) as well as the non-homogeneous soil (C cg , C ca , C ag in Fig. 2) is not examined. Therefore, the problem of estimating the core temperatures of unbalanced loaded multicore cables buried inside non-homogeneous soils still remains unsolved.

C. Technical Contribution of This Paper
All the aforementioned works present considerable limitations when implemented in unbalanced LV and MV distribution networks. Firstly, they use the mathematical formulas of IEC standard for calculating the thermal resistance and capacitance of the soil in the TEE. However, these formulas are only applicable to a limited number of installation geometries. When the medium surrounding the cable is composed of materials with different thermal resistivities and specific heat constants, which is usually the case in installations with a non-homogeneous soil or backfill material, standard methods cannot be applied [10]. Furthermore, the aforementioned works assume either singlecore or 3-core cables with balanced three-phase loading. In reality, a large portion of MV and LV cables are 3-core or 4-core cables with a highly unbalanced loading in each conductor, and thus, the temperatures of the cores can significantly deviate from each other.
In this manuscript, we propose a temperature-dependent power flow method for underground distribution networks with the following distinct features: r We propose a novel TEE model for the thermal simulation of underground cables, the parameters of which are estimated and validated based on the results of a Finite Element Method (FEM) software. Therefore, it is generic and applicable to every installation geometry and type of soil. r Our proposed model also considers the mutual heating caused by adjacent cable installations. r Our model accounts for the unbalanced loading of each phase in multicore cables, calculating individually and accurately the temperature of each core.
r The power flow is based on the method of [22], [23], which is applicable to AC and hybrid AC/DC networks operated in both islanded and grid-connected mode as well as in networks with different voltage levels [24].

II. CONVENTIONAL POWER FLOW ALGORITHM
The full configuration of an unbalanced LV network is presented in Fig. 1, including the neutral conductor and the grounding resistances.
Let us define the current and voltage vectors of node i as follows: where I Lir and V iy denote the load current and voltage (in complex form) of node i at phase r = {a, b, c} and conductor y = {a, b, c, n, g}, respectively, as shown in Fig. 1. For a network with m nodes, the load vectors can be expressed as a function of the voltage vectors according to Eq.
where V 0 corresponds to the slack node. Subsequently, for the power flow solution, we remove the first five rows of Eq. (3) that correspond to the slack node, and Eq. (4) is obtained. ⎡ ⎢ ⎣ Then, by transferring all the voltage variables of the left-hand side of Eq. (4) to the right-hand side, Eq. (5) is derived, where I new and Y new are the modified current and admittance matrices.
As a last step, we define the final matrices Y fin and Y fin . The first one consists of the first five columns of Y new , while the second one consists of the remaining columns such that Eq. (6) is then derived from Eq. (5) by subtracting the product Y fin V 0 from both equation sides.
Using Eq. (6), we finally derive Eq. (7), which is iteratively solved until a certain preset tolerance is reached. In Eq. (7), k denotes the iteration number and the vector V fin contains the voltages of all nodes except the slack.

A Introduction of the Proposed Thermoelectric-Equivalent Circuit
In this section, we propose a TEE circuit model for modeling a 4-core underground cable, as illustrated in Fig. 2. It consists of the examined 4-core underground cable, buried near to an adjacent thermal body. The soil is modeled through a set of a thermal resistance and capacitances between the surface of the investigated cable and the surface of the ground (R cg , C cg1 , C cg2 , C cg3 , C cg4 , C cg5 , C cg6 ), a set of a thermal resistance and capacitances between the surface of the investigated cable and the surface of the adjacent cable (R ca , C ca1 , C ca2 , C ca3 , C ca4 , C ca5 , C ca6 ) and a set of a thermal resistance and capacitances between the surface of the adjacent cable and the surface of the ground (R ag , C ag1 , C ag2 , C ag3 , C ag4 , C ag5 , C ag6 ), as depicted in Inside the cable, four heat sources exist (Q a , Q b , Q c , Q n ), at the position of the conductors, which denote their Joule losses due to current flows. The heat path between the conductors is denoted with a set of thermal resistance and capacitance (R m , C m ). Finally, the heat path between the conductor and the surface of the cable is denoted with a set of thermal resistances and capacitances (R s , C s1 , C s2 ), which represent the various layers of the cable. It is noted that the parameters C con , R m , C m , R s , C s1 , C s2 are the same in all conductors, due to the geometrical symmetry of the cable, which implies that the paths of the heat of the four conductors inside the cable will have equal thermal properties.
In order to calculate the parameters of the TEE circuit, we execute some simulations in a FEM software by modeling the examined underground installation. The thermal resistances of the TEE circuit of Fig. 2 are calculated based only on the results of FEM simulation, while the capacitances are calculated using the results of FEM simulation and Particle Swarm Optimization (PSO), as explained in the next sub-sections. TEE circuit modelling for underground cables is not a new concept and is widely applied in the literature, e.g., in [15]- [19]. In all existing TEE circuits so far, the estimation of the thermal parameters (thermal resistances and capacitances) is performed using analytical formulas, which are applicable only in specific types of underground installations under balanced loading conditions. For instance, the existing TEE circuit models found in the literature do not model mutual heating between the cores of multicore cables (expressed with the parameters R m , C m in the TEE circuit of Fig. 2). This is mainly because they either assume single-core cables or they have been based on the assumption that all the cores of multicore cables are equally loaded, and therefore no mutual heat flow exists between them. However, this assumption and simplification is not correct in the case of unbalanced loaded cores in multicore cables. Moreover, the existing models are restricted only to specific types of installations, for instance, only in homogeneous soils without backfill materials. The proposed model, on the other hand, has generic applicability in all types of soils under both balanced and unbalanced loading conditions, since its parameters are not estimated by analytical formulas but from FEM simulations and PSO.

B Calculation of Thermal Resistances
For the calculation of thermal resistances, the following simulations are executed in FEM: Simulation 1: Calculation of R s , R cg , R ag . We load the four conductors with a current of equal magnitude, so that the heat flowing through the R m thermal resistances is zero. More specifically, the four conductors are equally loaded with a current I∠0 o , I∠90 o , I∠ − 90 o , I∠180 o , respectively, so that they generate an equal heat Q con . Furthermore, the adjacent cable is loaded with a suitable current so that the temperatures at the surfaces of the two cables is equal and no heat is flowing through the resistances R ca . The total thermal resistance from the conductor to the surface of the cable (3 ·R s ) is calculated by Eq. (8): where T con is the temperature of the conductor, while T sur_cab is the temperature at the surface of the investigated cable, both calculated and measured in the FEM software. The total thermal resistance from the surface of the investigated cable to the surface of the ground (6·R cg ) is calculated by Eq. (9): where T sur_gr is the temperature at the surface of the ground, while 4 · Q con denotes the total heat produced by the equally loaded four conductors of the cable. The total thermal resistance from the surface of the adjacent cable to the surface of the ground (7 · R ag ) is calculated by Eq. (10): where T sur_adj is the temperature at the surface of the adjacent cable, while Q adj denotes the total heat produced by the adjacent cable, both measured by the FEM simulation software. Simulation 2: Calculation of R m . We load the conductors a , b, and c with I∠0 o , 2I∠180 o and I∠0 o , respectively. Moreover, the adjacent cable is loaded with a suitable current so that the temperatures at the surfaces of the two cables are equal and no heat is flowing through the resistances R ca . The heat produced by the conductor of phase b flows through the thermal resistances R m towards the conductors a and c as well as through the resistances R s , towards the surface of the cable according to Eq. (11).
where Q b , T a con T b con , T c con denote the heat produced by the conductor b, the temperatures of conductors a, b, c, respectively, all measured in the FEM simulation software, while R s is known from Simulation 1 above. R m is easily calculated from Eq. (11).
Simulation 3: Calculation of R ca . The four conductors of the investigated cable are equally loaded with a current I∠0 o , I∠90 o , I∠ − 90 o , I∠180 o , respectively, so that they generate an equal heat Q con . The adjacent cable is not loaded. The total heat produced by the investigated cable (4 · Q con ) flows through the thermal resistances R cg , R ca according to Eq. (12).
where Q con , T sur_cab , T sur_adj , T sur_gr are all measured in the FEM simulation software, while R cg has been calculated in Simulation 1. Therefore, R ca is easily calculated from Eq. .

C. Calculation of Thermal Capacitances
As shown in the TEE of Fig. 2, there is a total of 22 unknown thermal capacitances. Tο calculate their values, we perform time-series simulations in FEM software by setting specified current values for the conductors of the cable, in specified time periods, as shown in Tables I, II, III, and IV. A time-varying waveform for the conductor temperatures is obtained such as the ones shown in Figs. 8, 9, 10, and 11. The thermal capacitances are estimated so that the response of the proposed TEE circuit of Fig. 2 fits, as good as possible, with the response of the FEM simulation software, for the same conductor current values.
The estimates are realized using PSO, which is a heuristic optimization approach with global optimum searchability [25] and implementation simplicity. Although PSO presents an increased execution time to find the global optimum, this is not a problem in our application, because the optimization algorithm is executed only once (offline) for every cable installation and not in every load variation. Therefore, in our application, the main interest is the ability of the optimization algorithm to find the global optimum as opposed to a lower computation time. For this reason, PSO is a suitable optimization tool for this application. Of course, instead of PSO, other heuristic optimization algorithms with global optimum searchability could be also applied such as Genetic Algorithms, Ant Colony optimization, Artificial Bee Colony, etc.
The thermal circuit of Fig. 2 is mathematically expressed in a form of state-space equations, as follows: where denote the heat produced by the four conductors of the investigated cable, Q adj is the heat produced by the adjacent cable, while T sur_gr is the temperature at the surface of the ground. The vector u is the input of the thermal circuit of Fig. 2, which determines its response. The thermal resistances and capacitances are included in the matrices A, B.
Assuming that the states x 1 (t), x 2 (t), x 3 (t), x 4 (t) of the vector x denote the temperatures of conductors a, b, c, n, respectively, the following equation expresses the objective function that is minimized with PSO.
Where t n is the total simulation time by the FEM software. T r con (t) for r = {a, b, c, n} is the temperature of conductor r of the examined cable at time instant t calculated by the FEM software. The optimization variables are the 22 unknown thermal capacitances of Fig. 2 namely C con , C s1 , C s2 , C m , C cgi , C cai , C agi for i = {1, …6}. They should be constrained to have a positive value throughout the execution of PSO. In practice, Eq. (14) searches for the optimal thermal capacitance values in order to minimize the deviation between the dynamic response of the conductor temperatures obtained from FEM simulation and the proposed TEE circuit of Fig. 2.

D. Discussion About the Proposed TEE Circuit
The proposed TEE circuit of Fig. 2 does not include the diagonal thermal parameters between the diagonal cores b-n and a-c. In Fig. S3 of the supplementary material, we quote the form of the extended TEE circuit along with some guidance on how to calculate the mutual and diagonal thermal resistances of that circuit. However, based on our observation from the simulations, the neglect of the diagonal thermal parameters affects negligibly the accuracy of the TEE circuit. This occurs because, the diagonal parameters are indirectly considered in the TEE of Fig. 2 through the reduction of the R m value (please see the supplementary for more details). Therefore, we recommend the usage of the TEE of Fig. 2 for the modeling of unbalanced 4-core cables, because it has a simpler form than that of Fig.  S3, with almost the same accuracy (the difference is less than 0.5°C for the two circuits).
The proposed TEE circuit of Fig. 2 consists of 34 thermal capacitors to model the thermal inertia of soil and cable. It is pointed out that a lower number of capacitors could be applied to form the proposed TEE circuit, but in that case, the accuracy of the model could be sacrificed. On the opposite, a higher number of capacitors would only increase the complexity of the circuit without offering a substantial improvement in accuracy. Based on our understanding from the investigated topologies, Fig. 2 includes the ideal number of capacitors since it combines low modelling complexity with high accuracy.

IV. VALIDATION OF THE PROPOSED TEE CIRCUIT MODEL
The proposed TEE circuit model of Section III is validated here, using four different underground cable installations. The validation is performed by comparing the results of the proposed TEE circuit model against FEM simulation.

A. Investigated Installations
Four different underground cable installations are considered here, as explained below: r Installation 1: A 4-core LV cable is buried at a depth of 0.5 m inside backfill material, as presented in Fig. 3. The thermal resistivity and specific heat constant of the backfill are 0.45 Km/W and 850 J/(kgK), while for the native soil they are 1 Km/W and 1500 J/(kgK). The density of the backfill and native soil are 1700 kg/m 3 and 2600 kg/m 3 , respectively [16]. The 4-core cable was taken from [27, page 26] with a nominal cross-sectional area of 35 mm 2 . All conductor cores have the same cross-sectional area. An adjacent cable is installed at a distance of 25 cm.

B. Validation Results
The conductors of the cables of the aforementioned 4 installations are loaded over the time, as shown in Tables I-IV. It is pointed out that the simulation time length of Installation 3 is longer than the other three installations, due to the high thermal inertia of Installation 3, as a result of the deep depth, high specific heat and high density of the native soil.  In most time periods the three-phases are unbalanced loaded, which is usually the case in MV and LV networks. The conductor temperatures calculated by the FEM software and the proposed TEE circuit, for the four installations, are depicted in Figs. [8][9][10][11]. As shown, the response of the proposed TEE circuit is in good agreement with the simulations of the FEM software, over the whole investigated time period, confirming the accuracy of the proposed approach. More specifically for Installation 1, the maximum difference between the two models is 0.5°C at the 2 nd hour. In Installation 2, the maximum difference is 0.4°C at the 7 th hour. In Installation 3, it is 0.6°C at the 45 th hour, while in Installation 4, it is 0.6°C at the 3 rd hour.   Another important thing to note is the large deviation between the temperatures of the conductors, in periods where they are unbalanced loaded. For instance, in Installation 1, the temperature of the conductors a, b, c, n at the 7 th hour is 80.52°C, 90.73°C , 79.57°C, 71.1°C, respectively. It is evident that there is a temperature deviation between the conductors b and c, which is higher than 10°C. As a result, the assumption of a balanced loading and equal temperatures between the conductors, which is typically considered so far in the thermal modeling of multicore cables, could be significantly misleading in power system analysis.
It is pointed out that in Fig. 11, the neutral conductor has a similar temperature with phase c (74.44°C and 74.24°C, respectively at 3 rd hour), although the current through the neutral is significantly lower. This occurs due to the reduced crosssectional area of the neutral in Installation 4.

V. PROPOSED TEMPERATURE-DEPENDENT POWER FLOW ALGORITHM
Proceeding from the validation of the proposed TEE circuit model, we present in this section the proposed temperaturedependent power flow algorithm, as follows: Step 1: Firstly, the TEE circuit parameters of Fig. 2 are calculated offline for each cable installation, using the PSO optimization procedure proposed in Section III.
Step 2: After obtaining the parameters of the TEE circuit in Step 1, we formulate Eq. (13) to mathematically express the TEE circuit.
Step 3: We discretize Eq. (13) so that it is expressed in the following form: In this manuscript, the discretization was executed in MAT-LAB with the c2dm command, because of MATLAB fast and efficient inbuilt algorithms.
Step 4: Then, the iterative process begins by calculating the three-phase power flow using the algorithm described in Section II.
Step 5: From the currents obtained via the power flow solution in Step 4, we calculate the temperature of each conductor through Eq. (15). The vector u is updated through the Joule's law, e.g., Q a = I a 2 · R a ac , where I a is the current of conductor α calculated by the power flow in Step 4, while R a ac is the AC resistance of conductor a.
Step 6: From the conductor temperatures calculated in Step 5, we update the conductor resistances, as follows: where R θ is the resistance value at θ°C, R 20 is the resistance value at 20 o C, α is the temperature coefficient of resistivity, T θ is the conductor temperature at θ°C, and T 20 is the conductor temperature at 20°C. This iterative process is updated from Step 4 until convergence.

VI. TIME-SERIES TEMPERATURE-DEPENDENT POWER FLOW SIMULATION
A time-series temperature-dependent power flow is executed here, for a whole year, in a 25-bus unbalanced LV network, using the proposed algorithm.

A. Investigated Network
The examined underground 25-Bus network is depicted in Fig. 12 [3]. Every bus (except bus 1) supplies three single-phase residential loads, one in each phase. Thus, the network supplies in total 72 single-phase customers. The network was assumed to operate in grid-connected mode although the proposed power flow algorithm can be applied in islanded operation, as well [3].
For the loads, we used real measured data from [28], which were normalized, as shown in Fig. 13. Each load connected to phases a, b, c is composed from the product of the normalized load and the base load S base A = 1500 + j1125, S base B = 1800 +   j1350, S base C = 1600 + j1200, respectively, in order to form an unbalanced LV network. All loads were assumed to have a power factor equal to 0.8.
A single-phase photovoltaic (PV) is connected to phase a of buses 13,19,25. Each PV has a base power of 5 kW multiplied by the normalized solar irradiance of Fig. 14. The solar irradiance is measured in the region of Sindos (a suburb of Thessaloniki, Greece) [29], for the whole year of 2018. Fig. 15. Temperature of the ground's surface over the time using the model of [30].
The distance between the buses was assumed to be 70 m, which is a typical distance in LV networks. The buses are connected with underground LV cables. All the cables are assumed to be buried according to the Installation 4 (see Fig. 7). The adjacent cable is neglected in this study.
A grounding resistance (Z gri in Fig. 1) of 2 Ohm was assumed for all the buses of the network, which is a typical resistance value for all new electrical installations in Greece.

B. Modelling of the Ground Surface Temperature
The temperature of the ground surface significantly affects the thermal behavior of underground cables. This temperature varies both in a seasonal and in a daily mode. For the calculation of the ground's surface temperature over the whole year, in this paper, we applied the simplified model of [30] (see Eq. (11) of [30]) using the parameters of [30, Table 1]. These parameters have been estimated based on experimental observations of bare soil in the region of Athens, Greece. The temperature of the ground's surface is depicted, over a whole year, in the small sub-figure of Fig. 15. The large figure of Fig. 15 depicts the ground's surface temperature from 4000 th to 5000 th hour of the year. As shown, both the seasonal and the daily temperature variations are considered. In the summer season, the temperature of the ground's surface varies between 20°C and 45°C, while in winter from 3°C to 11°C.

C. Results About the Conductor Temperature and Resistance
The temperature of each conductor (three phases and neutral) for the cable connecting the buses 1-2 is depicted in Fig. 16. The smaller sub-figure on the top-left depicts the conductor temperature for the entire year. As expected, the temperature is proportional to the loading condition (see sub- Fig. 13) and ground surface temperature (see sub- Fig. 15), which take their maximum value in the summer period. The maximum conductor temperatures are observed on July 3. The corresponding conductor currents for this day are shown in Fig. 17.
As shown, the maximum peak temperatures of phases a, b, c, and neutral n are 77.81°C, 85.56°C, 82.72°C and 76.24°C , respectively. A temperature deviation of around 8°C is observed between the conductors of phase a and b. This large difference highlights the importance of simulating appropriately the unbalanced thermal loading of multicore cables, instead of   simplistically assuming that all cores have the same temperature, as occurs so far in the literature.
Finally, it should be noted from Fig. 17 that the conductor loading on July 3 is around 250 A, which is well above the static current rating of 158 A referred in [27, page 24]. Therefore, the proposed algorithm can constitute an accurate analysis tool for studying and utilizing the maximum dynamic current rating of LV and MV multicore cables.
The resistance of each conductor of the cable connecting the buses 1-2 is depicted in Fig. 18, during the whole year. It is clarified that in this example, we used a resistivity temperature coefficient equal to 3.39·10 −3 K -1 (see α in Eq. (15)), which is a typical value for copper. The neutral conductor has a higher resistance due to its reduced cross-sectional area. From this example, it is indicated that the conductor resistance varies significantly not only seasonally, but also hourly. Therefore, the assumption of a constant resistance adopted by the conventional  non-temperature-dependent (NTD) power flow (with fixed resistances) can be significantly misleading.

D. Power Flow Results of the Proposed and the Non-Temperature-Dependent (NTD) Approach
The proposed temperature-dependent power flow method is compared here against the NTD power flow method of [22]. For the NTD power flow, the resistances of the lines are assumed constant over the whole year at three different temperatures: at a) 40°C, b) 60°C, and c) 90°C. It is noted that 90°C is the maximum allowable conductor temperature for XLPE insulated cables, according to [27, page 18]. The voltage profile of the most loaded phase a, at two different time instants, is investigated in Figs. 19 and 20. More specifically, we investigated the following time instants: 1. the 4409 th hour, which corresponds to 3 rd July at 5 p.m. It is a hot hour with high loading. 2. the 860 th hour, which corresponds to 4 th February at 8 p.m.
It is a cold hour with high loading. As shown, a resistance value at 40°C leads to optimistic results, which are on the unsafe side, and therefore, studying power networks with NTD algorithms at this conductor temperature is not recommended. On the contrary, a resistance value at 90°C estimates larger overvoltages than actually occur, and therefore, it is a conservative consideration. A resistance value at 60°C seems to produce accurate results since its voltage profile almost coincides with the TD algorithm at the 4409 th hour. However, the results of Fig. 20 indicate that 60°C is a very and 20 that the resistance of the cables varies significantly over the year, and therefore, the assumption of a fixed resistance can produce misleading results, when studying the voltage profile of the network.
The total annual power losses of the network are depicted in Fig. 21 for all approaches. As shown, the proposed TD power flow algorithm estimates the lowest active power losses at 15.09 MWh. The NTD power flow calculates the losses at 15.67 MWh, 16.9 MWh and 18.7 MWh for the cases the resistances are fixed at 40°C, 60°C and 90°C, respectively. Therefore, significant deviations of the active power losses are observed between the proposed TD and the NTD approaches, highlighting the importance of accurately modeling the temperature effects on the resistance of underground cables. It is clarified that the difference of the reactive power losses is negligible between the examined approaches, due to the equal line inductances that are considered.

E. Computation Time of the Proposed Approach
The proposed approach needs an average of 0.05 sec to solve the power flow for one time instant. In this example, we executed the power flow with one hour resolution for a whole year, i.e., 8760 time instants. Thus, the total computation time of the algorithm over the whole year is around 7.3 minutes. The NTD power flow requires 0.03 seconds to solve the power flow for each time instant and 4.4 minutes, in total, for the whole year. The increased computation time of the proposed method is mainly due to the inversion of the Y BUS matrix of the network (the Y fin matrix in Eq. (7)) at each time instant, as a result of the variation of the conductor's resistances. On the opposite, in NTD power flow, the Y BUS matrix is factorized only once at the first time instant and remains constant over the whole year. Nevertheless, by using several techniques of numerical analysis, e.g., matrix decomposition [31], the computation time of the proposed method can be significantly reduced. A PC with a 64-bit Intel Core i7, 3.4GHz CPU and 16GB RAM was used for all simulations, for the sake of fairness.

F. Proposed TEE Versus Existing TEE Models
The thermal behavior of multicore underground cables depends on the geometry and size of the cable as well as on the geometry and thermal properties of the surrounding soil. The external thermal resistance (T 4 in IEC standard) is computed in both the IEC [11] and IEEE [12] standards using the Neher-McGrath method [32]. The method of Neher-McGrath is also applied by all recent papers focusing on the thermal rating and power flow of underground cables, e.g., [4], [13], [15]- [19]. However, the Neher-McGrath method is not applicable in cases, where the surrounding soil is composed of materials with different thermal properties like the Installation 1, 2 and 4 of Figs. 3, 4 and 7, respectively. According to the simulation results of [10, section III.B], errors up to 11% can be resulted in the standard methods, due to the presence of the backfill material.
Moreover, all the earlier works studying the power flow in underground installations, e.g., [4], [13], [15]- [19], have considered only single core cables, in which the layers of the cable, e.g., insulation, jacket etc., are placed concentrically around the conductor. In single-core cables, the thermal parameters (thermal resistances and capacitances) of the cable itself can be easily computed, using standard formulas. On the other hand, in our applications, we have multicore cables which are asymmetrically loaded. To the best of our knowledge, the thermal parameters between the cores of the cable (see R m , C m in Fig. 2) as well as between the cores and the surface of the cable (see R s , C s1 , C s2 in Fig. 2) can not be computed with analytical formulas.
In order to highlight the high accuracy of the proposed TEE circuit, under unbalanced loading, against the IEC standard and other existing TEE circuits, we performed simulations in Installation 3 (see Figs. 5, 6) using two different loading conditions: a) moderately unbalanced case of Table V (it is the same with  Table III), and b) highly unbalanced case of Table VI. For each loading, the positive-and negative-sequence currents are depicted in the last two columns of the tables. The zero-sequence current is zero due to the absence of neutral conductor, as it usually occurs in European medium voltage networks. The following TEE approaches are compared: 1) the proposed TEE approach,  Fig. 2], where the mutual resistances are neglected (practically, in Fig. 2, it means R m →Ý), 3) the commercial software CYME, which uses the IEC thermal equations and simulates only balanced loaded cables.
Since the currents of Tables V and VI are unbalanced, we employed CYME by equally loading the three cores with the positive-sequence current of Tables V and VI. For the "Moderately unbalanced case", the results of the proposed TEE circuit are shown in Fig. 10, while the results of the TEE circuit of [ [21], Fig. 2] and CYME are shown in Figs. 22 and 23, respectively. As shown in Fig. 10, the temperature responses of the proposed TEE circuit coincide almost completely with the FEM software for the loading of Table V.   On the opposite, the TEE of [ [21], Fig.2] presents considerable deviations from the FEM software, which reach up to 5°C at 45 th hour for phase a. Regarding CYME software, it calculates the same temperature for all phases ignoring the temperature unbalances between the cores. Nevertheless, CYME estimates a maximum temperature to around 89°C at 45 th hour, which is 5°C higher than the maximum temperature estimation of FEM software (at 45 th hour for phase a).
For the "Highly unbalanced case", the response of the proposed TEE circuit, the TEE circuit of [ [21], Fig.2] and CYME are depicted in Figs. 24, 25 and 26, respectively. Again, in Fig. 24, the response of the proposed TEE circuit coincides almost completely with the FEM. On the other hand, the TEE of [ [21], Fig. 2] deviates by around 7°C from FEM at 45 th hour. Finally, the maximum temperature of CYME in Fig. 26 is 95.3°C at 45 th hour, which fits well with the maximum temperature of phase a calculated by FEM.
From the aforementioned cases, it is concluded that the proposed TEE circuit calculates more accurate temperature responses for unbalanced loaded multicore cables, in comparison to the existing TEE approaches.

VII. CONCLUSION
This paper proposes a time-series temperature-dependent three-phase power flow method for unbalanced distribution networks, consisting of underground cables. A novel approach is proposed to estimate the parameters of the three-phase TEE circuit of multicore cables, using the results of a Finite Element Method (FEM) software and Particle Swarm Optimization (PSO). The proposed approach is generic and can be accurately applied to any kind of 3-core or 4-core cables buried in homogeneous or non-homogeneous soil. It is also applicable in cases, where one or more adjacent heating bodies exist. Using the proposed approach, the conductor temperature of each core of multicore cables can be separately and precisely calculated, even in networks with highly unbalanced loads. The proposed approach can constitute an accurate tool for power flow analysis of unbalanced distribution networks and estimation of conductor temperatures.
Additional validation results of the proposed TEE circuit of Fig. 2 are provided in the supplementary material, under extreme loading conditions. Moreover, in the supplementary material, explanations are provided about the diagonal thermal parameters between the diagonal cores of 4-core cables.