Thermal Assessment of Power Cables and Impacts on Cable Current Rating: An Overview

: The conceptual assessment of the rating conditions of power cables was addressed over one century ago, with theories based on the physical and heat transfer properties of the power cable installed in a given medium. During the years, the evolution of the computational methods and technologies has made more powerful means for executing the calculations available. More detailed conﬁgurations have been analysed, also moving from the steady-state to dynamic rating assessment. The research is in progress, with recent advances obtained on both advanced models, extensive calculations from 2D and 3D ﬁnite element methods, simpliﬁed approaches aimed at reducing the computational burden, and dedicated solutions for speciﬁc types of cables and applications. This paper provides a general overview that links the fundamental concepts of heat transfer for the calculation of cable rating to the advanced solutions that have emerged in the last years.


Introduction
Thermal phenomena are at the basis of power cable current rating. The temperature reached during the time in different points of the cable, especially at the interface between conductor and insulation, determines the cable lifetime, together with other ageing factors due to the characteristics of the materials [1,2]. Power cable design for a given application has to be carried out accurately, because any issues that appear in the cable during its operation could lead to failures. The occurrence of failures requires expensive maintenance, with the creation of joints (or junctions) that break the physical cable continuity and make the cable more vulnerable during successive operation [3].
When a cable is installed, the exact evolution in time of the current that the cable will carry, as well as the characteristics of the external environment, cannot be known in advance. For this purpose, the definition of the current rating of the cable has to come from general hypotheses on the cable operation for a given cable layout and type of installation. The current rating (also denoted as ampacity) refers to the maximum current at which the cable can operate without exceeding the temperature limits for the insulation material. Thereby, the determination of the current rating requires formulating a heat transfer problem, in which the thermal properties of the materials, the heat sources inside and outside the cable, and the mechanisms of heat dissipation are modelled and evaluated [4].
The time-variable nature of the heat sources, together with the possible changes occurring in the outside environment during the time, make the determination of the current rating a time-dependent problem. This problem is also indicated as dynamic cable rating, and consists of the determination of the maximum permissible current loading of the cable during the time. Dynamic cable

Heat Transfer Concepts for Thermal Analysis of Power Cables
This section addresses the heat transfer concepts used in the formulation of the thermal analysis of the power cables. Starting from the basic concepts, the illustration presented contains progressive indications, up to dealing with advanced and detailed aspects that have emerged from recent literature contributions.
The current flowing in a power cable generates heat, which is dissipated through the metallic layers of the cable and its insulation, towards the surrounding environment. Heat is an energy form transferred from one system to another in the presence of a temperature gradient. There are three heat transfer mechanisms involved in this process: conduction, convection, and radiation. The conduction phenomenon occurs through the metallic layers of the cable and its insulation if the cable is buried or is in the air or the water. Convection and radiation occur from the cable surface to the external environment. In the case of convection, the nature of the heat flow classifies this phenomenon into natural convection or forced convection. Natural convection arises when the flow is induced by the buoyancy forces, due to the density differences provoked by the temperature gradients in the air. Forced convection takes place due to the external means such as a with a pump, fan, or the wind. The interactions with the ambient environment makes the heat transfer phenomena more complex for cables installed in free air or water, with respect to the cables installed underground. The radiation refers to electromagnetic waves or photons. Radiation does not necessitate a medium, and its intensity strongly depends on temperature [17]. The heat is transferred by convection and radiation from the cable surface to its surroundings only for the power cables installed in air. The sun is considered an additional energy source if the power cable in the air is exposed to solar radiation [10].

Energy Conservation and the Energy Balance Equation
The energy conservation law is fundamental in the heat transfer assessment of a power cable. The appropriate equation of this law is: where: . Q in heat flow rate that enters in the power cable, and which is generated by the solar radiation for an insulated power cable installed in air, or by the neighbouring cables of a specific power cable buried in the soil; . Q gen heat flow rate generated inside a specific power cable, by Joule, dielectric and ferromagnetic losses; Q stor change of heat flow rate stored inside the power cable; and, . Q out heat flow rate dissipated by heat transfer mechanisms (or heat losses); in the case of the underground installations, the cable system also incorporates the surrounding soil. For an underground cable located in the soil, the conduction phenomenon occurs by all cable layers and from the cable to the soil. The cable length is much bigger than its diameter, the end effects are neglected, and for this reason, the general heat conduction equation is written in two-dimensions only [10]. In the basic model, the cable is assumed to be located in an infinite medium with uniform initial temperature. In this case, the heat conduction equation is written by taking into account the transient conditions that express the variation of the temperature in time: Energies 2020, 13, 5319 4 of 38 where T is the absolute temperature, ∂T ∂x and ∂T ∂y are the temperature gradients in x and y directions, ρ is the thermal resistivity, which is the inverse of the thermal conductivity k = 1 ρ typically used in the heat transfer domain, . q gen is the rate of heat generated per unit volume, c p is the specific heat capacity, ζ is the density, and τ is time.
This equation is a non-homogeneous partial differential equation. Its explicit solution exists only for specific geometries and boundary conditions. The complex geometry of insulated cables makes it impossible to obtain a closed-form solution. Generally, the problem is solved by numerical approaches, such as thermal-electrical analogy or Finite Element Method (FEM) [17].
In a homogeneous environment, the thermal resistivity ρ is constant, thereby the heat conduction equation is: where the thermal diffusivity a = 1 d·ρ·c p = 1 ρ·C v represents how fast heat is transmitted by a material, and C v is the volumetric heat capacity.
If cylindrical symmetry is assumed, the geometric variable considered is the radius r measured from the centre of the cable, and the heat conduction equation becomes: The solution of this equation is written in the form [9,10]: where, by indicating x = r 2 4·a·τ , the exponential integral −Ei(−x) is defined as: For underground cables, the main aspects that make the basic models not applicable depend on various effects, among which: (a) the non-infinite dimension of the soil; (b) the effect of the ambient on the soil properties; (c) the non-homogeneous soil; (d) the finite cable length; (e) the lack of cylindrical symmetry.
Specific indications on the above points are provided below.

Non-Infinite Dimension of the Soil
In Equations (4) and (5), the soil is considered homogeneous and infinite. In this case, with a single cable, the thermal field is composed of concentric lines. In practice, the cables are buried at a given depth, and the shape of the thermal field depends on the presence of the ground surface. In turn, the ground surface can be considered at the same temperature or different temperatures. The traditional solution considers Kennelly's hypothesis of having an isothermal interface between soil and air [18], further discussed in [19,20]. Kennelly's hypothesis leads to the definition of the image method, which is useful to calculate the temperature growth at any point of the soil. The cable system is considered an infinitely long cylindrical heat source buried in a uniform medium, and it is not possible to consider the convective boundary conditions at the soil surface [9,21]. The heat source + . q l (with thermal flux . q l in W/m) is buried at a given depth in uniform soil. Heat is transferred from Energies 2020, 13, 5319 5 of 38 the heat source to all points having lower temperature by conduction. If the heat source is enclosed in a duct or pipe containing air, then the heat is also transferred by convection. The soil surface is considered as the symmetry axis. The heat sink, emitting the heat − . q l , is the reflected image of the heat source. The heat sink has the same distance L and magnitude above the soil surface as the heat source has below the soil surface. In the case of multiple cables, or cables with multiple cores in which the electric current flows, the superposition effect is applied [10].
The heat conduction equation due to the heat sink is: The heat conduction equation due to the heat source is: The temperature increase at any point in the soil (in Figure 1, the point indicated is N) is the sum of two temperature increments (i.e., the temperature increment due to the heat source in the ground, and the temperature increment due to its fictive image above the soil): Energies 2020, 10, x FOR PEER REVIEW 5 of 36 surface is considered as the symmetry axis. The heat sink, emitting the heat − , is the reflected image of the heat source. The heat sink has the same distance L and magnitude above the soil surface as the heat source has below the soil surface. In the case of multiple cables, or cables with multiple cores in which the electric current flows, the superposition effect is applied [10]. The heat conduction equation due to the heat sink is: The heat conduction equation due to the heat source is: The temperature increase at any point in the soil (in Figure 1, the point indicated is N) is the sum of two temperature increments (i.e., the temperature increment due to the heat source in the ground, and the temperature increment due to its fictive image above the soil): Consider the hypothesis that the outside diameter of the cable is much lower than the distance from the surface of the ground to the cable centre. Therefore, the temperature growth at the outside cable surface is [10]:  Consider the hypothesis that the outside diameter D of the cable is much lower than the distance L from the surface of the ground to the cable centre. Therefore, the temperature growth at the outside cable surface is [10]: The effect of a non-isothermal ground surface is considered in [21] by adding a fictitious layer that moves the isothermal surface in the direction opposite to the cable location, as an application of the additional wall method shown in [22]. The adaptation proposed in the paper is based on applying the Fourier transform for converting the heat transfer problem from two dimensions to one dimension. In the transformed Fourier domain, the heat transfer coefficient is no longer dependent on soil thermal resistivity, installation depth, and cable dimensions; it depends only on the physical properties of the air and the heat that is dissipated by the cable. To model the non-isothermal soil surface, [21] proposed an accurate method in which the heat transfer coefficient may be computed in the Fourier domain. The application of this procedure was possible because, after the transformation in the Fourier domain, the heat transfer coefficient depends only on the air physical properties and the heat dissipated by the cable, and does not depend on cable dimensions, cable depth, and soil thermal resistivity ρ. This method has been shown to be compatible with standardized methods (IEC and IEEE).
The equation that describes the conversion of the two-dimensional heat transfer problem into a one-dimensional problem is: where x and y, respectively, are the horizontal and vertical coordinates, F{.} denotes the Fourier transformation, and s is the variable of the transformed Fourier domain. Furthermore, the cable (the heat source) is defined as a Dirac function f (x) that is applied at a depth L with respect to the soil surface, in which . q gen only appears at the position of the directly buried cable, and . q l = 0 anywhere else. In this case, the expression F T(x, y) = 1 points out that a buried cable in the space domain is expressed in the Fourier domain by a straight line with . q l = constant ( Figure 2). In the spatial domain, the soil temperature depends on the coordinate x, with the maximum temperature T s at x = 0 (i.e., above the cable), progressively reduced when the coordinate x increases. In the Fourier domain, the cable is considered a constant line source, and the soil surface temperature T s is constant. This temperature is iteratively calculated to verify the heat equation at the soil surface: where h is the heat transfer coefficient in the Fourier domain. The heat conduction equation (Laplace equation in steady-state, without . q gen ) after the Fourier transformation is: This equation satisfies the following boundary conditions: at  In this case, the heat conduction equation after the Fourier transform is Equation (13). Consider, the same line source. The boundary conditions are Equation (15) and the following equation: The temperature gradient is constant at y= : The solution of the equation with these boundary conditions is: By comparing Equation (16) with Equation (19), the thickness of the fictitious layer becomes: Figure 3 shows the scenario of when the fictitious layer is added. The mirror line is coincident with the horizontal isothermal line. In addition, different temperatures are found at the soil surface, which is crossed by other isothermal lines. The solution of Equation (13) satisfying the boundary conditions (14) and (15) is: The non-isothermal soil surface is modelled by introducing a fictitious layer with thickness d. The image of the cable line is then found at a distance L = L + d ( Figure 2).
In this case, the heat conduction equation after the Fourier transform is Equation (13). Consider, the same line source. The boundary conditions are Equation (15) and the following equation: The temperature gradient is constant at y = L : The solution of the equation with these boundary conditions is: By comparing Equation (16) with Equation (19), the thickness of the fictitious layer becomes: Figure 3 shows the scenario of when the fictitious layer is added. The mirror line is coincident with the horizontal isothermal line. In addition, different temperatures are found at the soil surface, which is crossed by other isothermal lines.

The Effect of the Ambient on the Soil Properties
Underground cables are buried at a given depth in the soil. The cable temperature can be affected by the phenomena occurring in the ambient, depending on the depth and on the presence of layers that reduce the effect of the external ambient, such as pavements or conduits. When the latter elements are absent, the variation of the temperature on the soil surface affects the internal part of the soil and the cable temperature. The effects can be seen on different time horizons.
In uniform soils, the heat flow is conductive. However, as indicated in [23,24], this assumption is not generally valid if there is a strong moisture gradient. The moisture gradient determines an isothermal vapour flux. In this case, part of the total soil evaporation results from subsurface evaporation that occurs under the soil surface. This evaporation affects the heat flux and the temperature profile.
The total soil heat flux is given by the following equation: where the isothermal latent heat flux depends on the presence of a moisture gradient, and represents the latent heat carried from evaporating subsurface layers by the isothermal vapour flux [24], while is the conductive soil heat flux, which incorporates only the thermal latent heat flux, and considering the soil depth is given by: During daytime, when a drying soil is heated, the thermal flux is positive, while the thermal flux is negative, and reduces the sum at the soil surface. The variation of ( , )

The Effect of the Ambient on the Soil Properties
Underground cables are buried at a given depth in the soil. The cable temperature can be affected by the phenomena occurring in the ambient, depending on the depth and on the presence of layers that reduce the effect of the external ambient, such as pavements or conduits. When the latter elements are absent, the variation of the temperature on the soil surface affects the internal part of the soil and the cable temperature. The effects can be seen on different time horizons.
In uniform soils, the heat flow is conductive. However, as indicated in [23,24], this assumption is not generally valid if there is a strong moisture gradient. The moisture gradient determines an isothermal vapour flux. In this case, part of the total soil evaporation results from subsurface evaporation that occurs under the soil surface. This evaporation affects the heat flux and the temperature profile.
The total soil heat flux . q tot is given by the following equation: where the isothermal latent heat flux . q latp depends on the presence of a moisture gradient, and represents the latent heat carried from evaporating subsurface layers by the isothermal vapour flux [24], while . q ds is the conductive soil heat flux, which incorporates only the thermal latent heat flux, and considering the soil depth z is given by: Energies 2020, 13, 5319 9 of 38 During daytime, when a drying soil is heated, the thermal flux . q ds is positive, while the thermal flux . q latp is negative, and reduces the sum . q tot at the soil surface. The variation of . q ds (z, τ) with the soil depth, depends on both temperature changes and the subsurface phase change (evaporation or condensation at the soil depth). The variation of . q ds (z, τ) leads to changes in the sensible and latent heats, to satisfy the energy conservation law: where . q gen is the rate of heat generated per unit volume for the local heat sink or source (e.g., water phase changes) in W m 3 . Neglecting . q gen and the spatial variations in k, Equations (22) and (23) give the uncoupled heat diffusion equation: In cylindrical coordinates, the uncoupled heat diffusion equation becomes: If heat and vapour flows are strongly coupled in the soil, the total soil heat flux is: .
where . q c is the conductive heat flux, and . q vT is the thermal latent heat flux (both thermal fluxes are proportional to − dT dz and can be combined as . q ds ) [24]. The energy conservation law in this case becomes the coupled equation [25]: The solution of the coupled equation needs numerical methods to be addressed. In fact, in general cases, many aspects cannot be solved analytically. With numerical methods, the soil is partitioned into adjacent volumes separated by boundaries. In this context, it is possible to model some particular cases with non-uniform soils, irregular boundary conditions, multi-dimensional flows, and non-linear equations.
Analytical solutions are possible under appropriate simplifications. In particular, these simplifications refer to considering uniform soil (or knowing the analytic representation of the variation of the soil thermal properties with depth) and assuming that the thermal properties of the soil do not change with the temperature. Further requirements to set up analytical solutions include the analytic representation of the boundary conditions and the initial conditions. Under these simplifications, it is possible to find the analytical solutions of two types of variation in time [24]: For periodic variations, starting from a time series with N values (where N is an even number), it is possible to calculate a number of harmonics not higher than H = N/2. The soil surface temperature is then expressed as: where ω 1 = 2·π p is the fundamental angular frequency, calculated by considering a period of 24 h when the daily wave is considered, or a period of 12 months for the annual wave; T s is the mean temperature at the soil surface; A h is the amplitude of the h-th harmonic; and ϕ h is the phase angle of the h-th harmonic.
The results of the Fourier analysis are rigorously valid when the wave is repeated without changes in successive time intervals, namely, each day for the daily wave, or each year for the annual wave. In addition, if the soil thermal properties do not vary with the depth, and the temperatures at all depths are assumed to vary around the same T s , the depth penetration of soil surface temperature T(z, τ) is the sum of the penetrations of each harmonic [26]: where D 1 = 2a ω 1 is the damping depth of the fundamental (h = 1). More generally, changes in T s with the soil depth can be addressed by studying the evolution of T s,z for different depths z, separately with respect to the sinusoidal terms.
From Equation (29), the conductive soil heat flux becomes: where For non-periodic variations, the Laplace transform of T(z, τ) is a function of z and the Laplace parameter s only [27]: The Laplace transform of Equations (24) and (25) is the ordinary differential equation [27]: Considering the difference T (z, τ) between T(z, τ) and an initial isothermal value of the soil, the boundary condition is expressed as: and for a semi-infinite soil a solution to Equation (31) is: where B is a constant, which is a function of s and depends on the boundary conditions used.

Non-Homogeneous Soil
A detailed analysis of the soil thermal properties is needed for the design and installation of underground cables and pipelines, to avoid premature damages. The main thermal property to assess the soil is its thermal resistivity. The main aspects that affect the thermal resistivity are the soil type, the geometrical layout of the soil components, and the moisture content. The composition of the soil includes different particles, which form aggregates of sand, silt, colloids, and pore spaces. In addition, water may be present in the soil into different extents [19].
The soil thermal resistivity under loading conditions is found in the range of 0.5 to 1.2 m·K/W [28]. This value could become higher because of the heat dissipated from the underground cables. The heat dissipation phenomenon provokes the soil instability, leading to thermal damage of the cable [29]. The soil thermal resistivity is affected by the soil temperature, soil porosity, and water-vapour transport. The corresponding parameters have to be modelled, to avoid errors in the temperature distribution compared with real conditions [30].
The standard IEC 60287 [31] indicates how the effective thermal resistivity and thermal resistance of the soil can be calculated. The ability of the soil to keep its thermal resistivity constant in the presence of a heat source is named thermal stability of the soil [32]. The soil moisture content has a significant effect on buried cables capacity and soil thermal resistivity. The major concern is that the heat transfer from the cable into the soil leads to a considerable moisture migration far from the buried cable in the case of unfavourable conditions. Furthermore, around the cable, a dry zone characterized by uniform high thermal resistivity could form. This phenomenon occurs because air, considered as a weak conductor, divides the solid particles of the soil. If the soil moisture content increases, the soil thermal resistivity decreases because water is a good conductor. In the dry zone, the thermal resistivity in the dry state is considered. The wet zone is characterized by a thermal resistivity in the saturated state [33]. As such, the dry soil has a thermal resistivity higher than the wet soil. If in the beginning, the soil thermal resistivity decreases quickly as the moisture content is raised after a while, the decrease rate comes to be lower [34].
Gouda et al. [35] specified that in the presence of the dry zone, the capacity of the buried cable is decreased by applying a derating factor that depends on the type of soil. Outside the dry zone, the soil thermal resistivity is also uniform but corresponds to the soil moisture content [9]. The dry zone provokes an increment of temperature in the cable sheath. This phenomenon leads to a deterioration of the cable insulation and an eventual formation of hot spots in the cable [36]. A practical assumption that can be considered for practical purposes is that the dry zone does not modify the profile of the isotherms compared with their profile when the soil was moist. Only the numerical values of some isotherms are changed [10]. When the temperature difference between the external surface of the cable and the ambient temperature exceeds a critical limit (which depends on temperature, type of soil, and moisture content), the drying of the soil forms a zone in which the thermal resistivity increases. The isotherm corresponding to the critical limit gives the boundary of the dry zone. In the dry zone, the uniform (high) thermal resistivity is assumed. Outside the dry zone, a uniform thermal resistivity is also considered. The only changes from the uniform conditions of the soil depend on the non-uniformity is caused by drying. The assumption considered is no longer valid in the presence of backfills with different characteristics from the uniform soil around the cable. The probability of the soil drying out increases when the direction of the cable is crossed by another heat source [10]. The moisture migration induced by thermal gradient changes the thermal environment and needs appropriate modelling of the temperature response of the cable [37].
Donazzi et al. [38] specified that some backfill materials such as sand, cement, silt, as well as water, can be used to improve the thermal conditions of the cables. For example, the dry zone phenomenon in the backfill, initiated at various temperatures and velocities, depends on the soil type and the quantity of mud [35]. Donazzi et al. [38] also highlighted the significance of the critical water content, supposing that this is independent of the environment temperature in practice (when the soil temperatures around buried cables are less than 80 • C). Groeneveld et al. [39] experimentally demonstrated the effect of temperature on the critical water content in the soil and specified that a higher temperature reduces the capacity of water-keeping of soil. The effect of ambient temperature and ambient saturation on the temperature of the cable conductor is addressed in [40] and [30]. Moreover, the critical temperature for the formation of the dry zone, and the ratio of the dry to wet thermal resistivity, depending on the soil components but not on the loading on the cable. An essential factor to obtain the time required for the dry zone formation is the heat flux [41]. Both the critical temperature and the ratio of the dry to wet thermal resistivity do not depend on the heat flux transferred from the buried cables to the soil [36]. The heat flux at the cable surface is useful to obtain the time required for soil to get [41]. Some mathematical models have been developed, assessing the dry zone phenomenon around the buried cable [42][43][44][45][46][47].
The hot spots formed in the buried cable can be mitigated by using some solutions such as: − the addition of a corrective backfill with low thermal resistivity; − the heat sources around the buried cables must be insulated; − the forced convection for the fluid around the buried cable; − an insulating fluid for the inner cooling of the cable; − installing in the hot spot zone a forced cooling system [48].
One of the most used solutions to avoid soil drying around the cable is based on water cooling [48,49]. The most adopted practice is to install, in parallel with the buried cables, the pipes in which the cooling water passes. In this case, the heat transfer modelling must take into account the soil parameters, depth of cable burial, location of the water pipes, and other factors varying along the cable route. This cooling method is not proper for the extruded insulation cables-except for the case in which a tight sheath for water is used. Tobin et al. [48] proposed a prototype of a chilled-water heat removal system, applied to underground urban distribution systems. The proposed system raised the possibility of carrying current with about 60% for 13 to 115 kV cables. More recently, Klimenta et al. [50] proposed a solution based on hydronic asphalt pavements. In the case in which the water cooling cannot be used to reduce the hot spots, Brakelmann et al. [47] propose gravitational water cooling as a solution.

Effective Soil Thermal Conductivity
The effective thermal conductivity (k eff ) is a thermal property of a multi-phase (air, water, solid) soil. It represents the soil capacity to transfer heat by conduction under unit temperature gradient, as a function of the volumetric fractions of the soil microstructure, soil phases, and the phases connectivity [51,52].
The main factors that influence the effective thermal conductivity of soil under isothermal conditions are: soil mineralogy, solid particle shape and size, gradation, cementation, porosity, packing geometry, water content (or saturation degree), soil temperature, and stress level [53]. A synthesis of these factors is provided below.
In soil mineralogy, the soil is considered a multi-phase system (a three-or four-phase system) because its composition includes solid (mineral) particles, water, air, and sometimes ice (in some cold zones) [15]. The solid particles contain soil minerals (e.g., one of them is the quartz with the biggest thermal conductivity 8.4 W/(m·K) compared to other minerals, see [54]), which are surrounded by water and air [55].
Furthermore, the thermal conductivity of the dry soil at a temperature of 10 • C is less than 0.5 W/(m·K) and depends on the packing density and mineral composition [56]. In the dry zone, the air impedes heat conduction, and this heat transfer phenomenon takes place through the contact points of the solid particles. If air is replaced by water, a significant enhancement of the heat conduction is observed. Therefore, the order of the thermal conductivities is k air < k dry soil < k water < k water−saturated soil < k mineral [15,[56][57][58].
Solid particle size and shape have a significant influence on the positioning of primary and secondary solid particles. In natural soils, the soft particles are included into bigger particles of different sizes and shapes. In addition, the number of contact points of the solid particles has a significant influence on the soil thermal conductivity. As it is known, the heat transfer in soils relies on the solid phase and occurs across the contact points, especially in the dry zone, because the air thermal conductivity is reduced compared to the solid particles of soil (the air thermal conductivity is 0.0026 W/(m· K)) [16]. Moreover, fewer contact points and bigger solid particles lead to increased soil thermal conductivity [59]. Gradation represents the distribution of various sizes of individual solid particles inside a soil zone. The soil with a good gradation presents a good heat transfer because the little solid particles fill the interstitial space of pores and raise the coordinates among the solid particles [15].
Cementation also influences soil thermal conductivity if the solid particles of soil are cemented together by binders or clay (the thermal conductivity is 1.28 W/(m·K) [54] the contact area will increase. The soil thermal conductivity will majorly increase [60,61].
Porosity influences soil thermal conductivity. The void ratio is a parameter to assess the compactness soil and represents the ratio between the voids volume to the solid volume. Based on the sketch of the volumetric ratios of soils shown in Figure 4, the void ratio is calculated as follows: Energies 2020, 10, x FOR PEER REVIEW 13 of 36 Packing geometry highlights that good coordination among the solid particles increases the soil thermal conductivity [60]. The packing density is the ratio between the solid volume to the total volume, The effect of the soil density on soil thermal conductivity is relatively low. The rise of the soil density leads to the significant growth of the contact point's number, but not at a substantial increment of the soil thermal conductivity [62].
Water content plays an important role in obtaining soil thermal conductivity. In unsaturated soils, the growth of the soil thermal conductivity with the water content increase highlights the significant contribution of the pore conduction [15,34]. Water movement also influences soil thermal conductivity. At temperatures less than 0 °C the water frozen in the soil and the soil thermal conductivity are modified. Conversely, at high temperatures, water is changed into water vapour molecules, and the soil thermal conductivity increases.
The influence of the soil temperature on the soil thermal conductivity is analysed in [63]. The thermal conductivity of nine soil samples at the soil temperature ranging from 30 to 90 °C was measured. The results showed that in moist soil, the thermal conductivity had a significant increment with the soil temperature, obtaining values three to five times the 30 °C value when the sample soil temperature was 90 °C. Hiraiwa and Kasubuchi [64] noted that in the case of sandy soils, the soil thermal conductivity increased as the temperature of the soil increased. Tarnawski and Gori [65] measured the thermal conductivity of the soil for soil temperatures ranging from 5 to 90 °C in the case of four soil moisture content domains. They demonstrated that soil thermal conductivity has a low variation with the soil temperature and water content at reduced moisture contents. Smits et al. [66] measured thermal conductivity for two grains of sand with different solid particle sizes for a temperature range from 30 to 70 °C and variable saturation. They observed that the increase of the soil thermal conductivity with the temperature is for temperatures higher than 50 °C at low to intermediate saturation. When the soil is close to saturation, and at the lowest saturations, the temperature did not have a measurable effect on the thermal conductivity. At the temperatures ranging from 30 to 50 °C, the thermal conductivity has a small variation with the temperature.
Stress level also plays an important role, in the sense that higher stress leads to higher contact radii resulting in an increment of the soil thermal conductivity. In addition, under higher stress, the granular chains enhance the heat transfer in soil [15,67]. Many thermal conductivity models have been proposed to obtain accurate predictions of the effective thermal conductivity. The models have The lower the void ratio, the greater the thermal conductivity [57]. Porosity is the ratio between the voids volume and total or bulk volume of soil.
Packing geometry highlights that good coordination among the solid particles increases the soil thermal conductivity [60]. The packing density is the ratio between the solid volume to the total volume, The effect of the soil density on soil thermal conductivity is relatively low. The rise of the soil density leads to the significant growth of the contact point's number, but not at a substantial increment of the soil thermal conductivity [62].
Water content plays an important role in obtaining soil thermal conductivity. In unsaturated soils, the growth of the soil thermal conductivity with the water content increase highlights the significant contribution of the pore conduction [15,34]. Water movement also influences soil thermal conductivity. At temperatures less than 0 • C the water frozen in the soil and the soil thermal conductivity are modified. Conversely, at high temperatures, water is changed into water vapour molecules, and the soil thermal conductivity increases.
The influence of the soil temperature on the soil thermal conductivity is analysed in [63]. The thermal conductivity of nine soil samples at the soil temperature ranging from 30 to 90 • C was measured. The results showed that in moist soil, the thermal conductivity had a significant increment with the soil temperature, obtaining values three to five times the 30 • C value when the sample soil temperature was 90 • C. Hiraiwa and Kasubuchi [64] noted that in the case of sandy soils, the soil thermal conductivity increased as the temperature of the soil increased. Tarnawski and Gori [65] measured the thermal conductivity of the soil for soil temperatures ranging from 5 to 90 • C in the case of four soil moisture content domains. They demonstrated that soil thermal conductivity has a low variation with the soil temperature and water content at reduced moisture contents. Smits et al. [66] Energies 2020, 13, 5319 14 of 38 measured thermal conductivity for two grains of sand with different solid particle sizes for a temperature range from 30 to 70 • C and variable saturation. They observed that the increase of the soil thermal conductivity with the temperature is for temperatures higher than 50 • C at low to intermediate saturation. When the soil is close to saturation, and at the lowest saturations, the temperature did not have a measurable effect on the thermal conductivity. At the temperatures ranging from 30 to 50 • C, the thermal conductivity has a small variation with the temperature.
Stress level also plays an important role, in the sense that higher stress leads to higher contact radii resulting in an increment of the soil thermal conductivity. In addition, under higher stress, the granular chains enhance the heat transfer in soil [15,67]. Many thermal conductivity models have been proposed to obtain accurate predictions of the effective thermal conductivity. The models have been divided into three model types: mixing models, empirical models, and mathematical models.
The mixing models, called theoretical/physical models, consider the soil as a three-phase system (solid, water and air) in which the phases are represented as a particular combination of series and parallel in the soil sample [15]. The mixing models are based on the mixing laws (arithmetic, geometric, and harmonic mean) of the series model and parallel model. The series and parallel models refer to Wiener bounds (or upper and lower bounds) of thermal conductivity and do not depend on the pore structure of porous medium [68]. Combinations of such series and parallel models are extensively presented in [15].
The series model considers a constant heat flux through each layer ( Figure 5a). The phases have different thermal conductivities and develop different temperature gradients [15].  In this case, the effective thermal conductivity of soil is: where is the number of phases (M = 3 for air, water, and solid), n is the volume fraction of each phase, and n is the thermal conductivity of each phase.
The parallel model considers a different heat flux through each phase that depends on the thermal conductivity of each phase (Figure 5b). The phases develop the same temperature gradients.
In this case, the effective thermal conductivity of soil is: In this case, the effective thermal conductivity of soil is: where M is the number of phases (M = 3 for air, water, and solid), µ n is the volume fraction of each phase, and k n is the thermal conductivity of each phase. The parallel model considers a different heat flux through each phase that depends on the thermal conductivity of each phase (Figure 5b). The phases develop the same temperature gradients.
In this case, the effective thermal conductivity of soil is: Another theoretical model is the De Vries model [23] and its simplified form presented in [69]. The De Vries model considers the soil as a mixture of ellipsoidal particles aleatory placed in the wet soil (continuous water) or the dry soil (ambient). In the wet soil, the solids and air represent the dispersed phase (the phase that is present in the particle shape), and in the dry soil, the solids and water represent the dispersed phase. The De Vries model requires many shape factors of ellipsoidal particles that are difficult to obtain. In this case, the effective thermal conductivity is very influenced by the shape factors and water content [70].
The effective thermal conductivity is: where M is the number of phases (air, water, solid); G i is the percentage of the mean temperature gradient of each phase. It is influenced by the shape factors of the ellipsoidal particles and soil components. The case i = 0 represents the continuous phase, θ i is the volume fraction of each phase; G i is the percentage of the mean temperature gradient along with each phase; and k i is the thermal conductivity of each phase. The model developed by Gori [71] was based on a soil cubic mixing cell and consists of a comparison between the analytical predictions and the experimental results of the effective thermal conductivity. The hypothesis to obtain the effective thermal conductivity was to consider parallel and horizontal isotherms or vertical heat flux lines. The experimental data obtained in [71] on the unsaturated frozen soils showed that the hypothesis with parallel and horizontal isotherms better predicts the effective thermal conductivity.
In the case of saturated frozen or dry soils, a cubic cell representing the soil was placed inside a cubic space (Figure 6a). The unfrozen water was distributed along the cubic cell length L t . Furthermore, the ratio between the cubic cell length L t and the solid length L s depends on the soil porosity.
In the case of the parallel and horizontal isotherms, the analytical expression of the effective thermal conductivity is: where k eff is the effective thermal conductivity for parallel and horizontals isotherms, ρ T = 1 k T is the effective thermal resistivity, and k c is the thermal conductivity of continuous phase (for the frozen soils k cp = k i and for the dry soils k cp = k a ). The other thermal conductivities refer to the ice (k i ), the air (k a ), the water (k w ), and the solid (k s ). Furthermore, u w = V uw V tot cell is the ratio between the unfrozen water volume and the total volume of the cubic cell, V uw is the unfrozen water volume, V tot cell is the cubic cell volume. In the Gori model, the solid phase is placed in the centre, with the bridges increasing around the solid particle or water films. The remaining space is occupied with the air. The change of solid, water and air in the dry condition, low moisture condition and unsaturated condition with bridges around is shown in Figure 6. Other theoretical models are widely presented in [16] and [51]. The theoretical models depend on a number of correlated factors, such as the ones recalled at the beginning of this subsection. Hence, the formulation of these models is quite challenging [52].
The empirical models are based on the mathematical and numerical assessment of experimental values between effective thermal conductivity and the measured soil properties (such as the degree of saturation or temperatures) [51]. The empirical models determine the expression between the relative thermal conductivity and the water content or saturation degree by normalizing the effective thermal conductivity [72]. The normalized thermal conductivity affirms that the effective thermal conductivity is estimated with the Kersten number, that is a linear combination of the saturated thermal conductivity sat and the dry soil conductivity dry [73]: The saturated thermal conductivity sat is determined with the geometric mean method: where is the porosity. The dry thermal conductivity dry can be expressed in function of the porosity [73]: The upper bound and the lower bound that appear in Equation (44) are not considered for the sake of facilitating the calculation. New empirical models have been formulated by various researchers by modifying this model [16,51,52]. The mathematical models are adapted from predictive models of physical properties. These properties are electric and hydraulic conductivities, dielectric permittivity, and magnetic permeability. These models are calculated by using the mathematical methods that take into account the volume fractions and thermal conductivity of each phase [51]. The details about these models are extensively presented in [16,51].

Finite Cable Length
The general hypothesis for carrying out studies of cable rating is that the cable length is virtually infinite. However, in some cases, there are situations in which this hypothesis is no longer In the case of the unsaturated frozen soils, the analytical expressions of the effective thermal conductivity are more complex and are computed for two cases: the water content is less than the ratio V wa V v (Figure 6a); the water content is higher than the ratio V wa V v (Figure 6b). In the Gori model, the solid phase is placed in the centre, with the bridges increasing around the solid particle or water films. The remaining space is occupied with the air. The change of solid, water and air in the dry condition, low moisture condition and unsaturated condition with bridges around is shown in Figure 6.
Other theoretical models are widely presented in [16] and [51]. The theoretical models depend on a number of correlated factors, such as the ones recalled at the beginning of this subsection. Hence, the formulation of these models is quite challenging [52].
The empirical models are based on the mathematical and numerical assessment of experimental values between effective thermal conductivity and the measured soil properties (such as the degree of saturation or temperatures) [51]. The empirical models determine the expression between the relative thermal conductivity and the water content or saturation degree by normalizing the effective thermal conductivity [72]. The normalized thermal conductivity affirms that the effective thermal conductivity is estimated with the Kersten number, that is a linear combination of the saturated thermal conductivity k sat and the dry soil conductivity k dry [73]: The saturated thermal conductivity k sat is determined with the geometric mean method: where σ is the porosity. The dry thermal conductivity k dry can be expressed in function of the porosity [73]: The upper bound and the lower bound that appear in Equation (44) are not considered for the sake of facilitating the calculation. New empirical models have been formulated by various researchers by modifying this model [16,51,52]. The mathematical models are adapted from predictive models of physical properties. These properties are electric and hydraulic conductivities, dielectric permittivity, and magnetic permeability. These models are calculated by using the mathematical methods that take into account the volume fractions and thermal conductivity of each phase [51]. The details about these models are extensively presented in [16,51].

Finite Cable Length
The general hypothesis for carrying out studies of cable rating is that the cable length is virtually infinite. However, in some cases, there are situations in which this hypothesis is no longer true. An example is the short-conduit, typically used for providing additional protection to buried cables that cross the streets or pass near other pipelines [74]. The length of the short-conduit and the conditions of heat dissipation depend on the specific case. The short-conduit cable is composed of a buried section and a conduit section (Figure 7). The classical methods, as well as the method indicated in the standard IEC 60287 [75], cannot be applied directly to this situation. An approximation is possible only if the heat transfer along the axis is not modelled. However, in actual conditions, the heat dissipation is better in the buried section. Hence, part of the heat generated in the conduit section reaches the buried section in the axial direction, creating an axial temperature gradient [73]. For this reason, the calculations carried out by ignoring the axial phenomena lead to incorrect results for the assessment of the thermal cable phenomena, with possible underestimation of the maximum temperature in the buried section and overestimation of the maximum temperature in the conduit section. The latter would create a disadvantage for the full usage of conduit sections [73]. approximation is possible only if the heat transfer along the axis is not modelled. However, in actual conditions, the heat dissipation is better in the buried section. Hence, part of the heat generated in the conduit section reaches the buried section in the axial direction, creating an axial temperature gradient [73]. For this reason, the calculations carried out by ignoring the axial phenomena lead to incorrect results for the assessment of the thermal cable phenomena, with possible underestimation of the maximum temperature in the buried section and overestimation of the maximum temperature in the conduit section. The latter would create a disadvantage for the full usage of conduit sections [73]. The effects of the length of the short-conduit cable, buried depth and soil resistivity have been investigated in [76] through 3D thermal simulations. In general, the initial point for the computation of the cable radial temperature has been based on a known thermal parameter of the environment [77,78]. In practice, however, the soil thermal conductivity is variable. For short-conduit cables, the accurate determination of the thermal environmental parameters in real-time is a challenging task.  The effects of the length of the short-conduit cable, buried depth and soil resistivity have been investigated in [76] through 3D thermal simulations. In general, the initial point for the computation of the cable radial temperature has been based on a known thermal parameter of the environment [77,78]. In practice, however, the soil thermal conductivity is variable. For short-conduit cables, the accurate determination of the thermal environmental parameters in real-time is a challenging task. For this purpose, a specific method for determining the current real-time capacity of short-conduit cables has been proposed in [73]. A simplified quasi-3D thermal model is set up, and an iterative procedure is formulated starting from the real-time temperature at the conduit surface, measured in the conduit section, the temperature at the cable surface in the buried part, as well as the current that flows in the cable. The axial heat flow is updated during the iterations, together with the other variables.

Lack of Geometrical Symmetry
Even when the cable is modelled as a system with cylindrical symmetry (also including possible junctions), the overall system around the cable could not be geometrically symmetric. For example, for cables in conduits, a simple geometrically symmetric structure could adopt a concentric configuration. However, this configuration is not physical and does not correspond to the symmetric distribution of the temperatures, mainly because of convection. In [79], physical non-symmetric cases with "eccentric" and "cradle" configurations have been analysed, indicating that one of the main shortcomings in setting up a thermal model could be the over-simplification of the cable location in unfilled conduits. The effect is that the average temperature at the conductor-insulation interface in the case of concentric formation can be higher than for the eccentric and cradle configurations.
Other causes of lack of symmetry appear in helicoidal cables [80], in which the location of the cores changes at different horizontal cross-sections. Modelling such a system requires a software tool able to support 3D simulations, and the preparation of the geometry would be better supported by the definition of the symmetric structure, followed by the application of a torsional operator.
Non-symmetry also appears when a failure occurs in a given point of the cable. In this case, the properties of the materials change around the failure point and can create more non-symmetry when the failure evolves in time. This thing makes the cable modelling in transient conditions more complex, but is also relevant to the identification of failure paths, leading to hot spots and progressive deterioration of the materials. In this case, lack of symmetry also appears when there are multiple cables close to each other and the effect of the cable failure during time causes thermal effects to other cables, with temperature variations due to the mutual thermal coupling.

Non-Buried Cables
The underground cables mainly addressed in this paper cover most of the existing installations. However, there are other installations for which the thermal model has to be adapted. Some cases are described in chapter 10 of [9], including cables located within protective walls (such as covered trays, protective risers, or tunnels), and cables placed in uncovered trays. When the cables are installed in air, the main heat transfer mode becomes the convection, either natural, or forced (if an airflow along the cable is present), and radiation. In contrast, conduction in the air is typically neglected [9].
Submarine cables are used in high voltage direct current electric transmission systems, generally with single-core cables. For submarine cables, in any case, part of the path close to the cable terminals is buried in the soil, and is typically the bottleneck on the thermal side, with the possible occurrence of high temperatures. Among the various solutions, using a cooling pipe with circulating water would improve the thermal conditions [81]. In this case, the thermal model becomes more elaborated, to include the effects of the solutions for cooling the cable.
In recent years, with the progressive integration of offshore wind farms, the number of three-core submarine cables has increased. For these cables, specific thermal analysis is needed, taking into account the internal and external characteristics. Submarine cables are subject to thermo-mechanical stresses, which can be analysed in a multi-physics environment [82]. Further issues that appeared, owing to the increased usage of submarine cables for offshore applications, include the crossings among submarine cables, which can cause temperature increase at the crossing points. These cases cannot be handled by using the classical image method, and their analytical formulation would also depend on the model of the subsea thermal environment [83]. In particular, in this thermal environment, the effect of sediments on the heat transfer characteristics is not negligible [84]. The sediments are a mixture of solid and liquid components, which impact on conductive heat transfer. Moreover, if the permeability of the sediments is high, the convective heat transfer can be significant [85]. The thermal properties of the sediments (e.g., thermal conductivity and volumetric heat capacity) have to be considered in the overall thermal model.

Thermal Models of Power Cables and Electrothermal Analogy
The thermal resistance R t is the ratio of the temperature difference between the two faces of a material to the rate of heat flow per unit area. The thermal capacitance C t represents the ability of a material to absorb and store the heat for using it later. In the electrothermal analogy, the thermal circuit is modelled by representing an equivalent electrical circuit, in which currents are equivalent to heat flows, and voltages are equivalent to temperatures. If the thermal parameters are independent of temperature, the equivalent circuit is linear. In this case, to solve a heat transfer problem, the superposition principle is applied [10].
The conductive thermal resistance of a cylindrical wall (e.g., cable insulation) having the inner diameter d i and the outer diameter d e is: where ρ is the thermal resistivity. The conductive thermal resistance of a plane wall is: where δ is the thickness of the wall, and S cond is the cross-section of the wall. The convective thermal resistance is: where S conv is the convective heat transfer surface, and h conv is the convective heat transfer coefficient. The radiative thermal resistance is: where S rad is the radiative heat transfer surface, and h rad is the radiative heat transfer coefficient. The total heat transfer coefficient for a cable installed in air is: Energies 2020, 13, 5319 20 of 38 The thermal capacity C t for a coaxial configuration (e.g., cylindrical insulation of a cable) having the inner diameter d i and the outer diameter d e , is: where V is the volume of the cylindrical configuration. The thermal conductance is a measure of the rate of heat flow through a body, and it is the reciprocal of the thermal resistance: The typical thermal models of the cables have been constructed by resorting to the electrothermal analogy (Table 1). Mostly, the resistors represent the thermal resistances; the capacitors represent the thermal capacities, which are essential to model the thermal transients in real operating conditions. The generators represent the heat sources due to different types of losses (e.g., Joule losses in the conductor, dielectric losses, and losses in the sheath and armour) [11,86].

Electrical Parameters Symbol and Unit Thermal Parameters Symbol and Unit
Electric conductivity In Table 2, the values of thermal resistivity and thermal capacity for some materials and soils in different conditions are reported. Typical values of other coefficients for the materials used in cables and backfill materials can be found in [10,31].
Specific representations of the equivalent circuit considered to model the cable in operational conditions have been presented for the cable and the outer part (the soil and the ambient). The soil can be divided around a buried cable into many concentric layers. Compatibility with the IEC standards is obtained by representing each layer with an equivalent thermal circuit of T-type, composed of the thermal capacitance of the layer, as well as the thermal resistance of the layer divided by two (see Figure 8) [87]. Table 2. Values of thermal resistivity and thermal capacity for some materials [10]. Specific representations of the equivalent circuit considered to model the cable in operational conditions have been presented for the cable and the outer part (the soil and the ambient). The soil can be divided around a buried cable into many concentric layers. Compatibility with the IEC standards is obtained by representing each layer with an equivalent thermal circuit of T-type, composed of the thermal capacitance of the layer, as well as the thermal resistance of the layer divided by two (see Figure 8) [87].

Insulating materials
The thermal resistance st for each layer of the soil is given by: where s is the thermal resistivity of the soil, s is the soil thickness, and is is the inner diameter of the soil layer.  The thermal resistance R st for each layer of the soil is given by: where ρ s is the thermal resistivity of the soil, δ s is the soil thickness, and d is is the inner diameter of the soil layer.
The thermal capacity C st for each layer of the soil is given by: where d es is the outer diameter of the soil layer, and C t soil is the thermal capacity of the soil. The RC ladder of the soil is obtained by linking together all the RC ladder layers of the buried cable ( Figure 9). The electrothermal circuit is obtained by linking together the ladder model of the cable with the ladder model of the soil [88]. The equivalent thermal resistances of the soil are calculated as: Energies 2020, 10, x FOR PEER REVIEW 21 of 36 The thermal capacity st for each layer of the soil is given by: where es is the outer diameter of the soil layer, and t soil is the thermal capacity of the soil. The RC ladder of the soil is obtained by linking together all the RC ladder layers of the buried cable ( Figure 9). The electrothermal circuit is obtained by linking together the ladder model of the cable with the ladder model of the soil [88]. The equivalent thermal resistances of the soil are calculated as: (54) Figure 9. The lumped thermal parameter network of a buried cable and the soil.
To represent fast transients, an RC circuit with a low RC time constant is used. Conversely, the long transients, which are useful for determining the heat to be transferred at far soil layers, are studied with an RC circuit with a high RC time constant [87].
The relevant temperature rise e ( ) in transient conditions is the one determined at the external surface of the hottest cable. The IEC Standard 60853 establishes how to calculate this temperature rise. Let us consider the distance , from the centre of the hottest cable p to the centre of the generic cable w in a group of cables (for single-core cables, the number of cables is = 1). Let us further consider the distance , between the centre of the hottest cable p and the centre of the image of the generic cable w.
For long transients, the temperature rise is: For fast transients, the effect of the images is considered negligible, so that: To represent fast transients, an RC circuit with a low RC time constant is used. Conversely, the long transients, which are useful for determining the heat to be transferred at far soil layers, are studied with an RC circuit with a high RC time constant [87].

Methods for Thermal Analysis of Power Cables
The relevant temperature rise T e (τ) in transient conditions is the one determined at the external surface of the hottest cable. The IEC Standard 60853 establishes how to calculate this temperature rise. Let us consider the distance d p,w from the centre of the hottest cable p to the centre of the generic cable w in a group of N C cables (for single-core cables, the number of cables is N C = 1). Let us further consider the distanced p,w between the centre of the hottest cable p and the centre of the image of the generic cable w.
For long transients, the temperature rise is: For fast transients, the effect of the images is considered negligible, so that:

Methods for Thermal Analysis of Power Cables
The thermal analysis of power cables and overhead lines is typically conducted with analytical and numerical methods [78]. In the sections below, the main methods for the two categories are described. Table 3 summarizes some differences among the considered methods. The characteristics taken into account are: the computational burden, which represents the resources required by a computing machine to solve the problem; versatility, which describes the possibility to model complex scenarios Energies 2020, 13, 5319 23 of 38 characterized for example by non-homogeneous geometries or materials; geometrical dimension, which defines the number of geometrical dimensions of the simulated domain; multi-physics approach, which pinpoints the possibility to carry out a simulation combining several physical domains. The temperature distribution inside and outside the cable is computed by analytically solving the heat diffusion equation.
The increase of the transient temperature of the external surface of a cable with respect to the temperature of the soil can be determined by considering the cable as a heat source immersed in a homogeneous medium in which the initial temperature is uniform. In this case, at any point in the soil the transient temperature T(τ) is expressed by using Equation (5) having the solution indicated in Equation (6). The principal issue of this analytical approach is that inhomogeneity of materials and cable structures cannot be taken into account [89]. Moreover, it is possible to solve the heat flow equation in different coordinate systems by separating the variables only when the internal heat generation is null. Therefore, only a limited number of geometries and boundary conditions can be studied with this approach [90]. Some techniques to partially overcome this problem have recently been proposed. For instance, an extension of the approach provided by IEC 60287 is presented in [87]. In this method, which stays simple-to-use, the soil is divided into concentric layers, each one characterized by a thermal resistance and capacitance forming a lumped circuit. Unfortunately, at finite depth, the isotherms around a buried power cable are highly asymmetric, and cannot be easily modelled with concentric layers of the soil. In [91], the soil is thus analytically modelled through non-concentric layers.

Numerical Methods
Numerical methods can be adopted to solve the heat transfer diffusion equation. The main numerical methods used in literature to determine the temperature distribution inside the cable and in the external environment are: Finite Difference Method (FDM) [45,92], Finite Element Method (FEM) [93][94][95], and thermal-electrical analogy [17]. FDM and FEM require making a mesh of discrete points in which the temperature is computed. In FDM, mesh grids are generally in cylindrical and rectangular coordinates, while in FEM, different mesh shapes can be chosen, taking into account the geometry of the objects under analysis. Techniques to increase the number of mesh elements have been developed. It is considered that the temperature gradient is higher or closer to the point of interest in the analysis [96]. Both for FDM and FEM, the speed of solution increases as the number of points studied decreases, but at the expense of the solution accuracy [45].
For both methods, inhomogeneity of the materials can be taken into account. However, only FEM allows the simulation of complex scenarios. In particular, FEM is adopted to investigate the impacts of the trench geometry and of the backfill material type and formation, such as for example the presence of multiple circuits [97], ground surface heat [89], cable trench profile [98], concrete and asphalt cover [76,99,100], and mixtures for bedding [101].
From the mathematical formulation point of view, FDM requires to solve the finite difference temperature equations that have the form of Equation (57) to compute the temperature distribution [92]: where T i,j is the temperature at the node of the mesh (identified by the indexes i, j), and α, β, γ, ξ are constant values, in the function of the location of the central node, of the thermal characteristics of the soil, and the heat transfer propagation type (conduction or convection). Following the FEM method, the transient temperature equation for a cable partitioned into N f finite elements can be written in the general form reported in Equations (58) and (59) [93]: .
where C is N f by N f heat capacitance matrix, Θ is an N by 1 column vector of node temperatures, K C and K h are N f by N f thermal conductance and convection matrices, and . q gen (τ) and . q h (τ) are N f by 1 column vectors of heat fluxes arising from internal heat generation and surface convection.
Some software-based FEM (e.g., COMSOL Multiphysics and ANSYS) allows multi-physics simulations that are particularly convenient in determining the dynamic thermal rating of electric cables [102,103]. For example, if the electromagnetic and temperature fields are computed together, simulation results are expected to be more realistic since they are not based on the simplified hypothesis of field distribution [104,105]. In this case, the currents in the three-phase current system, in the screen and into the soil can be computed considering mutual coupling and skin effects. Therefore, the temperature distribution can be calculated more accurately [104]. Moreover, thanks to the multi-physics approach, the molecular properties of the conductors and insulating materials that have an impact on the temperature distribution field can be taken into account. For instance, the interface between metal and polymer, which is influenced by different parameters such as the resonance or mismatch of phonon vibrational mode frequencies and the morphology of the insulation material (crystalline, amorphous and lamellae), can be evaluated [106,107].
Unfortunately, multi-physics analysis usually increases the complexity of the simulation. In fact, considering multi-physics domains often implicates multiscaling, which means that the combined physical models have significant differences in space or time scales [104]. For example, in the cable current rating evaluation, an increment of the complexity could be due to the different mesh density requirements of the combined physical models: the eddy current calculation in the thin solid screen requires several layers of finite elements across the screen thickness, whereas for heat transfer calculation a single layer is enough. Researchers usually adopt this approach only if the usage of the same mesh for both the domains implicates an acceptable increment of the computational time [104]. Especially for cable line ordinary constructions and depending on the purpose of the simulation, simplifying hypotheses can be assumed without having a significant impact on the results [104,108,109]. In [109], a comparison between the results obtained through the simplified approach proposed by IEC 60287 [31] and through multi-physics FEM simulations is provided for various cable line layouts. For an underground three-phase line formed by three single-core cables (cross-section of 630 mm 2 , XLPE insulation) in flat formation, the conductor temperatures differ with about 4 • C, which means a percentage deviation of around 5% if the FEM simulation are assumed as by reference [109]. Vice-versa, for complex scenarios such as the cable duct described in the paper (10 rows, 4 ducts per row, in which both 185 mm 2 and 240 mm 2 single core cables are installed, XLPE insulation), simplified IEC 60287 model provides lower temperatures of the cable (difference up to 27 • C) with a percentage deviation of around 20% [109].
When simulation purposes are required to consider large volumes/surfaces, and therefore geometrical objects characterized by significant different dimensions are involved in the simulation (e.g., a cable line and a large soil volume), analytical models are usually more appropriate than numerical methods that need to mesh the domains. The mesh creation process could result in a complex task. The critical point of FDM and FEM consists of the high computational resources needed to make the mesh and solve the heat diffusion equations.
The thermal-electrical analogy, also known as the matrix approach, is the third numerical method described in this paper and presents similarities with the FDM [88,110]. As already mentioned in Section 3.1, this method exploits the analogies between heat transmission and electrical equations. The cables and the surrounding environment are modelled through electrical components such as resistors, capacitors and generators. The state variable temperature T is represented as the electric potential at different nodes, and the heat flows are solved through an RC circuit. Generally, it is assumed that the length of the cable is much longer than its diameter, that no axial variations occur and that the thermal flux distributes only in the radial direction [17]. Moreover, the thermal resistance of the conductor is not considered, when the heat source is located in the outside section of the conductor. The system of first-order differential equations that represents the circuit can be written in a matrix form reported in Equation (60): where χ(τ) is the vector of the state variables, . q gen (τ) is the input vector determined by the Joule losses in the conductor core, T a (τ) is the ambient temperature, and A, B and c are the dynamic matrix, the input coefficient matrix, and the disturbance vector, respectively.

Steady-State and Dynamic Cable Rating
The previous sections have been dedicated to the thermal quantities, and to the ways to determine the temperatures in the cable and the surrounding environment. The next paragraphs introduce the electrical quantities and present the methodologies to compute the cable current rating. In particular, the steady-state and dynamic cable rating approaches are described. Furthermore, methods to consider the impact of harmonic currents and the uncertainty of the input parameters (e.g., thermal resistivity of the soil, environmental temperature, cable loading) are presented.

Steady-State Cable Rating Calculations
The current rating of the cable can be calculated as the continuous current carried by the cable, such as the continuous conductor temperature will be equal to the maximum allowable conductor temperature (this value depends on the insulation material). With these assumptions, steady-state conditions are assumed for the useful life of the cable. Starting from the thermal model of the cable and of the surrounding, the IEC Standard 60287 [31] provides the current rating equations with a constant load (i.e., 100% load factor) taking into account all losses arising in the cable (Joule losses, dielectric losses, armour and screen losses, etc.). The insulation surrounding the conductor is represented with a Π equivalent model. The permissible current rating I r for a cable with n load-carrying conductors can be calculated as: where: ∆T is the allowable conductor temperature rise above the ambient temperature, given by the difference between the permissible maximum conductor temperature and the ambient temperature; W d are the dielectric losses for the insulation surrounding the conductor; R t1 is the thermal resistance between one conductor and the sheath; R t2 is the thermal resistance between the sheath and the armour; R t3 is the thermal resistance of the external serving of the cable; R t4 is the thermal resistance between the cable surface and the surrounding medium; λ 1 is the ratio of losses in the metal sheath to total losses in all conductors; λ 2 is the ratio of losses in the armouring to total losses in all conductors; R e is the electric resistance of the conductor evaluated at the maximum allowable conductor temperature.
In particular, Equation (61) is applied for buried cables where drying out of the soil does not occur, or for cables placed in air. The IEC Standard 60287 [110] provides equations that can be applied for buried cables in the presence of partial drying-out of the soil and where drying-out of the soil is to be avoided.
On the basis of Equation (61), the cable current rating depends on the electrical and thermal parameters of the cable and on the thermal parameters of the soil.
When a higher value of the current rating is required, the analysis of Equation (61) can help understanding how to vary the parameters to increase the current ratings. A review on this problem is provided in [11], and the main results of [11] are summarized here. The cable current rating can be increased by reducing the thermal resistances, that is, by reducing the burial depth, by increasing the cable spacing, by using thermally controlled backfill with very low thermal resistivity or with natural or forced cooling; or by reducing the electrical resistance, by using special configurations for the conductors such as insulated wire Milliken-type conductor. Of course, the increase of the admissible temperature, obtainable when insulation materials exhibiting better thermal performances are used, will result in a greater current rating.

Dynamic Cable Rating
The application of the steady-state current rating of cables could resultingly be conservative in some cases due to the variation of the conductor temperature during the useful life of the cable according to the inputs of the thermal models. Hence, during the years, the assessment of the thermal stress of the lines has passed from the steady-state thermal rating to a more general Dynamic Line Rating (DLR) or Dynamic Thermal Rating (DTR), able to characterize the thermal transients and their consequences better. The DLR is aimed at determining the actual current rating of the line on the basis of continuous measurements or solution of the thermal model of the cable. Conservative assumptions on input variables, such as those assumed in the steady-state cable rating calculations, are no longer considered. Considering dynamic current ratings instead of static ratings allows increasing the estimated capacity of the cable.
While many contributions appeared in the past with respect to overhead lines, underground cables have been recently taken into consideration to apply the concepts of dynamic rating [87,111].

Effect of Harmonics on Cable Rating
The thermal-electric model initially proposed for the definition of cable current rating was extended to account for the presence of harmonic distortion. Indeed, Joule losses, as well as dielectric losses, are affected by harmonic distortion. The problem is complex, and some proposals are available in the relevant literature.
Harmonic currents in distribution systems are increasing due to the increased usage of power electronics-based appliances, converters and, recently, of inverter interfaced distributed generation and electric vehicle chargers. The presence of current components at harmonic frequencies imposes to evaluate the behaviour of the cables at those frequencies, and several aspects have to be considered. First, the harmonic components of the current contribute to the heating of the cable. Moreover, the resistance of the cables varies with the frequency and, in particular, at frequencies higher than the fundamental power frequency, the resistance increases due to the skin effect and proximity effect. Specific attention is needed for the neutral conductor (when present) because of the non-cancellation of zero sequence harmonic currents. The additional heating due to harmonic distortion may lead to a higher cable temperature. Then, it has to be considered in the determination of cable current rating when the cable carries a distorted current.
In several papers addressing these issues, the effect of harmonic currents supplied by the cable on its ampacity is handled by defining proper derating factors of the cable [112][113][114]. The derating factor is defined as "the ratio of the RMS value of a distorted current with a specific harmonic signature to the RMS value of a current of the fundamental frequency that produces the same losses in the cable as the distorted one" [99]. It can be calculated only when the model of the cable at harmonic frequencies is available, and the harmonic signature of the current is given.
Meliopoulos and Martin [112] address the problem of the evaluation of the cable ampacity when the cable supplies highly distorted current. In particular, low voltage supply systems are considered. The proposal is based on the extension of the Neher-McGrath equations to account for the additional losses due to harmonics.
The Joule losses of the cable for distorted current P loss are calculated in [112] as: where R eh is the conductor resistance at the hth harmonic, I h is the RMS value of the hth harmonic current, and H is the maximum harmonic order. The evaluation of the resistance at the hth harmonic requires the application of specific relationship able to account for the dependence of the skin effect and the proximity effect on the frequency. Models for evaluating the resistance of a cable conductor at a generic frequency are available starting from classical books. Additionally, in [112], the most relevant formulas are reported.
To calculate the losses, the spectrum of the current carried by the cable has to be known. A derating factor is proposed in [112] to determine the distorted current that produces the same losses of an undistorted current. Given the harmonic signature {I B , α 1 , . . . , α H }, being I B the base RMS current, and α h per unit value of the hth harmonic with respect to the base value I B, the derating factor is defined as: where R e1 is the resistance at the fundamental frequency, and I B is the base RMS current. In [113] an expression is proposed for the derating factor depending on the "harmonic signature" and on the ratios between the conductor resistances at harmonics and the conductor resistance at the fundamental power frequency.
A finite-element analysis is applied in [46] to analyse the effect of harmonic currents on PVC-insulated, low-voltage (0.6/1.0-kV) power cables symmetrically loaded and placed in free air; four-conductor cables (three phases and neutral) are considered and cables with a cross-section of the neutral conductor equal to or less than that of the phase conductors are taken into consideration. The results reported in [46] indicate that the derating factor depends on the cable configuration and on the type of non-linear loads the cable will supply. In [99], the influence of the metallic tray on the ampacity derating factor is evaluated; in particular, it is demonstrated that the derating factor increases with the cable cross-section.
An application of the harmonic derating factor to pipe-type cables is provided in [114]. The particular case of concentric neutral cables used in North America for power distribution systems, is addressed in [115], while cables with impregnated paper insulation are considered in [116].
A recent contribution [117] is aimed at including the presence of harmonic currents generated by electric vehicles (EV) chargers in the evaluation of temperatures and ampacity of medium voltage (MV) cables. The model used in [117] is an extension of the IEC formula for the assessment of the ampacity. Specifically, the variation of the resistance with the frequency, as well as the variation of the losses in the sheath layer and in the steel armour layer, is considered. While in non-distorted conditions Equation (60) has to applied, in distorted conditions (and neglecting the dielectric losses) the current rating I r,distorted conditions of MV cables is evaluated as in [117]: where T max is the permissible maximum conductor temperature, T a is the ambient temperature, and λ 1h and λ 2h are the ratios of the losses in the metal sheath layer and in the steel armour layer, respectively, at the hth harmonic, to the total conductor losses.
In [117], several cases have been considered with different electric vehicle chargers, and the ampacity of a three-core XLPE Medium Voltage cable is evaluated with respect to typical charging profiles. The results demonstrate that the derating factor decreases as the cross-section of the cable increases. Moreover, the derating factor depends on the harmonic distortion of the current and, therefore, depends on the typology of electric vehicle chargers that are supplied. The values of the derating factor evaluated in [117] range from 89% (in the worst scenario, for single-phase uncontrolled rectifiers rated at 6.6 kW) to about 99.7% (for uncontrolled rectifier topology with power factor correction).

Probabilistic Models and Risk Analysis for Calculation of Current Rating
As demonstrated in Section 4.1, the current rating of a buried cable depends on several factors; among these, the soil thermal resistivity and the ambient temperature are recognized to be random in nature [118], therefore, to account for this randomness, some authors proposed a probabilistic approach.
In [119] the temperature of the cable is evaluated considering the random changes of the thermal resistivity of the native soil and of the backfill, the ambient temperature, and the cable loading. In particular, the thermal resistivity is linked to the soil moisture content, and the relationship is investigated. A Monte Carlo simulation is performed to select random values for the uncertain parameters and, then, for each set of values of moisture (and, then, of the thermal resistivity), of ambient temperature, and of load current, a finite element analysis allows to evaluate the temperature of the conductor to obtain accurate predictions of the cable current rating.
In [120] a method to calculate underground cable current rating based on accurate evaluation of soil thermal resistivity is proposed. Experimental thermal resistivity probability distributions at a selected site are obtained from monitoring of soil thermal resistivity and rains. The results reported in [120] demonstrate that the thermal resistivity of the soil can vary over a month and can also maintain a low value for regular rainfall patterns. Starting from the thermal resistivity probability distributions, the authors derive conservative assumptions.
Shabani and Vahidi [121] propose a procedure aimed at optimizing the current rating of underground cables allocated in backfill considering the uncertainty of parameters such as soil thermal resistivity, ambient temperature, and load current. A Monte Carlo procedure is set to simulate the random variables and, for each set of random numbers, an optimization problem is solved. The objective function includes the cost of the backfill and the deviation of load current and cable ampacity. The result of the procedure is the probability density function of the cable ampacity. Zarchi and Vahidi [122] apply the Hong Point Estimate Method to characterize the random temperature of underground cables in duct banks. The uncertain input variables are the ambient temperature, the soil thermal resistivity, the backfill thermal resistivity and the burial depth; for each of them, a proper probability density function is assumed. An optimization of the cable configuration, based on this method, is solved and the total ampacity of underground cables in duct banks has been related to the chosen confidence level.

Conclusions
This paper has addressed the thermal models of underground cables, starting from basic models with general hypotheses towards the adoption of more detailed specifications to address practical cases. Heat transfer concepts needed for deriving the thermal model of an underground cable have been summarized, and applications to particular cases have been defined (e.g., the non-infinite dimension of the soil, cable with a finite length). Moreover, the electrothermal analogy has been applied to the cable thermal model, since it has been widely used in the domain of power systems. The methods used to simulate the heat transfer in the cables and in the external medium have been summarized to provide an overview of the main contributions. Finally, the cable thermal models have been applied to the determination of the cable current rating.
Underground power cables have been extensively studied and modelled for decades; therefore, many contributions are available in the relevant literature. This review paper has summarized the historical aspects and has delineated the recent evolutions.
From the literature review, it emerges that there are various developments in progress, which need further improvements. With respect to the methods, detailed FEM representations, also in 3D, of cables with the surrounding soil and ambient also in non-uniform conditions are becoming viable, thanks to the computational speed available. However, there is also an interesting development of improved models and simplified models that can provide results comparable with FEM. A topic that is becoming more attractive is the application of accurate thermal models of cables in the field of the dynamic line rating, recognized as an action that allows better utilization of the cable lines in many conditions, and is also useful for postponing the investments to upgrade the installed cables.
Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflict of interest.