The Corona Phenomenon in Overhead Lines: Critical Overview of Most Common and Reliable Available Models

: Research on corona discharge, shared by physics, chemistry and electrical engineering, has not arrested yet. As a dissipative process, the development of corona increases the resistive losses of transmission lines and enhances the line capacitance locally. Introducing additional losses and propagation delay, along the line, non-linearity and non-uniformity of the line parameters; therefore, corona should not be neglected. The present work is meant to provide the reader with comprehensive information on the corona macroscopic phenomenology and development, referring to the most relevant contributions in the literature on this subject. The models proposed in the literature for the simulation of the corona development are reviewed in detail, and sensitivity curves are provided to highlight their dependence on the input parameters.


Introduction
Corona may be defined as a non-linear phenomenon involved in the initial phase of electrical discharges, resulting in the flow of electric energy from a conductor to the ionized medium [1]. It consists of self-contained localised discharges and of the proximity of electrode surfaces characterized by local high intensity electric fields, stressing the insulating strength of the surrounding dielectric medium [2]. The phenomenon was studied rigorously for the first time by Faraday in 1838 [3].
Some of the very first references to this phenomenon date back to 1838, when Michael Faraday noticed a "quiet phosphorescent continuous glow", and stated that "that form of disruptive discharge which appears as a glow [. . . ] seems to depend on a quick and almost continuous charging of the air close to, and in contact with, the conductor" [3].
Early relevant and extensive studies on corona were performed by Peek in the first decades of the 20th century [4,5]; when applying a threshold voltage to a conductor, a "hissing noise" and a pale luminosity were observed in the proximity of its surface [6]. Indeed, the high frequency spectra typical of discharge processes result in the emission of audible noise and radio interference in the MHz range, research being oriented toward prediction and mitigation of these dissipative side effects [7].
Interaction between charged particles under the accelerating effect of an external electric field, turning into an avalanche phenomenon and, finally, into a process of local discharge by breakdown of the insulating medium, raised the interest of industries. The mechanism of recombination of ionized oxygen atoms, resulting from the discharge process, and oxygen molecules in the air is exploited for the production of ozone (O 3 ) on a utilityscale through ozone generators operated by corona technology [8], which need to trigger efficient discharges in a controlled environment in order to optimize ozone production.
Therefore, since Peek's early works on this topic, a relevant contribution, supported by experimental data, was represented by the book of Gary and Moreau [9], addressing also the computation of the localised losses caused by the corona discharge. Researchers the operation value, are associated with space charge movement under the effect of the alternating electric field [17].

Radio Interference
At higher frequencies, in the range between 0.01 MHz and 30 MHz, radio interference (RI) is a well-known disturbance phenomenon, which may be caused by corona development in proximity of HV d.c. and a.c. lines, affecting the quality of AM radio broadcasting [18]. The International Special Committee on Radio Interference (CISPR), involved in EMC issues and standardization, addresses this topic in detail in [19], with reference to RI due to overhead power lines. To assess the impact of corona on signal transmission at radio frequencies, three aspects are considered: frequency spectrum, lateral profile (i.e., attenuation of the field at ground level in the direction perpendicular to the source), and statistical distribution, to be evaluated following proper measurement procedures [18,20]. Indeed, different weather conditions, features of the conductors surface (e.g., dirt, irregularities, etc.) strongly affect RI levels, up to 20-30 dB [21], and may be taken into account through long-term surveys and statistical analysis. Several empirical and analytical approaches have been proposed in the literature to predict the impact of corona audible noise and RI on the surrounding environment [22] under variable atmospheric conditions, and guide overhead HV lines and substations design to limit emissions.

Corona Inception and Development
Corona is a self-sustained discharge process, developing in regions with high-intensity electric fields, which includes ionization of the surrounding dielectric medium, drift and diffusion of charged species, i.e., electrons, or charged atoms or molecules [23]. Microscopic or macroscopic approaches may be adopted to study and simulate this phenomenon.

Microscopic Level
In order for the corona inception to happen, a ionization process needs to be started by the production of ions and electrons at an exponential rate, also denoted as an avalanche mechanism, which was first addressed by Townsend [24]. With reference to a cylindrical electrode in the air at atmospheric pressure, the main process involved in avalanches development is ionization of neutral molecules by the collision with primary electrons; indeed, these electrons may be available in low concentration in the atmosphere and accelerated by the external electric field, or may be detached from existing negative ions and accelerated until collision. A statistical time lag τ st , in the range of µs, is necessary for the formation of enough primary electrons to start the avalanche process. Experimental evidence showed τ st to depend on the availability of free electrons, and on the rate of rise of the applied voltage; lower mean values of τ st are associated with steeper voltage fronts [2]. Ionization by impact is generally modeled by Townsend's coefficient α, which accounts for the average number of ionizing collisions expected for a single electron along a 1 m path and depends on the air pressure and electric field distribution [25].
When all the primary electrons are involved in collisions, secondary electrons are necessary to act as "seeds" for the following avalanches. They may be produced by several mechanisms, e.g., photoionization for positive corona, which consists in the production of positive ions and electrons resulting from photon-molecule interactions [26], and secondary emission for negative corona, which corresponds to electrons emission by positive ions approaching the electrode surface and being successively neutralized [27].
Outside the ionization region, the only processes of interest are charge movement through drift (depending on the electric field and ion/electron mobility) and diffusion of charged species (linked to the gradient of the space charge density and ion/electron diffusivity), and subtraction of space charge through recombination and attachment. Furthermore, corona experimental evidences have shown different discharge modes with respect to a DC or AC applied voltage source.

DC, AC, Impulse Corona
Depending on the nature of the voltage source applied to the electrode under test, different corona modes develop; the work by Giao and Jordan [28] displays a detailed analysis on this topic, whose fundamental aspects are reported here: • Negative DC corona: when a negative polarity voltage close to the inception value is applied, negative corona starts in the form of current pulses, known as "Trichel pulses"; the pulses' frequency increases with the increasing value of the applied voltage. The intermittent nature of these localised discharges is related to the space charge, which contribute to weaken the electric field in the electrode proximity (holding back the pulses) and to subsequent space charge removal. Further increase in the voltage leads to a pulseless glow, which may turn into negative streamers of growing length if the voltage reaches values close to the insulation breakdown [13]. • Positive DC corona: from experimental data [29], the inception voltage for positive DC corona is slightly lower, compared to negative DC corona. According to [28], this is reasonably due to the electrons being accelerated in the direction of the increasing electric field, i.e., towards the anode, facilitating ionization processes. Mixed opinions may be found in the literature on this subject; for instance, in [13], corona discharge due to negative polarity applied voltages is stated to occur at lower field values, with respect to positive polarity (the difference being enhanced when sharp-shaped electrodes are considered). This is in accordance with experimental values found in [30], and with values for the critical corona onset gradient in [31]. However, it is shared opinion that experimental activities aimed at assessing the corona onset field are strongly dependent on atmospherical conditions, humidity, configuration and available free charges in the vicinity of the conductor under test. Fundamental steps in the study of positive DC corona are found in [32]. It first develops in the form of burst pulses; depending on the electrode radius, gap length and magnitude of the applied voltage, some onset streamers may be produced as well, reaching larger radial distances from the electrode and being choked off by negative ions in the gap. Increasing voltage and higher densities of negative space charge suppress the onset streamers in favor of a diffused glow in the proximity of the anode (Hermstein's glow); streamers finally causing the insulation breakdown are observed at higher voltages, developing from the glow region [28]. • AC corona: when a sinusoidal voltage source triggers corona, combined modes typical of positive and negative DC corona may be observed. In particular, a critical distance can be introduced, depending on the voltage peak value and frequency, accounting for the path covered by the charge particles during the quarters of periods from the voltage peaks to zero crossings; if the critical distance is such that no residual space charge is left by the end of the half period (i.e., suitable time has passed for the ions to migrate toward the electrode of opposite polarity and being neutralized), corona modes in the positive and negative half cycles are expected to be equal to the DC modes with corresponding polarity. Otherwise, the residual space charge will influence the development of the corona modes during the subsequent half-cycle with opposite polarity [28]. • Impulse corona: with fast front voltages, corona mainly develops in the streamer mode, being the role played by the space charge limited by its slower dynamics [33].

Macroscopic Level
Ionization and interactions among charges at microscopic level provoke a non-linear and hysteretic relation between the total charge and the voltage at the electrode surface at macroscopic level, as sketched in Figure 1. When a monotonically increasing voltage is applied to a reference conductor, a proportionality factor, i.e., the conductor capacitance, governs the q-v relation until corona develops, in correspondence to the inception voltage V inc ; a further voltage increase results in the space charge growing at a rate faster than linear, albeit delayed of time τ sp in the range between 25 µs and 100 µs in the hemisphere-plane gap configuration, depending on the gap length, and on the voltage rate of the rise and peak value [2,34] necessary for its formation [26]. Due to this delay time τ sp , the peak value of the charge may be observed when the time derivative of the applied voltage has already reversed in sign; at later times, for decreasing values of voltage, the linearity of the q-v relation is restored. The last branch of the loop (D and E in Figure 1) shows a larger slope with respect to the one of branches C and D. At point E, the dynamic and the magnitude of the electric field produced by the applied source are no longer predominant, compared to the electric field produced by the space charge, which contributes to the mechanism of space charge subtraction, resulting in a value of the derivative dq/dv larger than C geo [35]. Numerous experimental activities [12,36], supported by theoretical considerations [2], have put in evidence that the electric field inception value E inc depends on the steepness of the applied voltage. As a result, different formulae for the computation of E inc have been proposed in the literature.

Inception Criterion
In this section, the approaches from Peek [5], Olsen and Yamazaki [11], and Mikropoulos et Zagkanas [37] are reviewed, their applicability is clarified, and the main ideas and considerations underlying the formulation are discussed. The majority of the approaches do not consider the influence of the steepness of the applied voltage on the value of the inception field, their validity being limited to DC or power frequency applications. Criteria accounting for this aspect, as well as for non-coaxial arrangements for the direct application to transmission lines (TLs) studies, are found to be lacking in the literature.
Since these criteria refer to solid conductors, a correcting factor, i.e., the irregularity factor, is often exploited to compute the inception voltage associated with stranded conductors and bundles. Its role is also discussed in this section.

Inception Field by Peek
Early experimental research activities by Peek [4][5][6] led to the well-known expression for the critical inception field E inc of the insulating air, i.e., the value of the electric field at which visual effects of corona, such as luminosity of the air surrounding the conductor, may be first observed: where K P = 0.3 in the original formulation by Peek. As expected, the inception field in (1), expressed in kV/cm, only depends on the conductor radius r 0 , and not on any other geometrical feature of the configuration under test (which, on the contrary, affects the onset voltage); δ is the air density factor, and f p is the polarity factor, which takes into account the effect of voltage polarity on E inc since experimental evidence has shown that the inception electric field associated with positive corona may hold different values with respect to those related to negative corona (with reference to the coaxial configuration, [38]); m is a constant to take into account the irregularity of the conductor surface, distinguishing between perfectly smooth conductors and rough surfaces due to weather conditions and aging under operation. From the manipulation of (1), assuming δ = 1 and axial-symmetric configuration, Peek's criterion may be alternatively interpreted as the value of the electric field at the conductor surface allowing to reach a field equal to E 0 at distancer = r 0 + 0.301 √ r 0 from the conductor axis, E 0 being independent from the configuration under analysis. Indeed, according to (1), the inception field holds higher values for conductors, presenting a smaller radius r 0 , as the proposed formula, although expressed in terms of r 0 , relies on a requirement for the electric field atr.
In the literature, many empirical expressions similar to (1) may be found, computing E inc through the following general relation: where coefficients A, B, and C (ranging between 23-35, 0.15-1, and 0.3-0.5, respectively) hold different values for the different expressions proposed in the literature, which are classified in [39] for d.c. corona.

Inception Field by Olsen and Yamazaki
When dealing with the evaluation of corona inception field at power frequencies, the onset criterion by Olsen and Yamazaki [11] takes into account the predominant interactions among free charges in proximity of the conductor; in particular, differently from the approaches in Section 3.3.1, the proposed corona onset criterion is associated not only with the maximum electric field at r 0 , but also with the electric field trend in the area surrounding the conductor surface, still not accounting for any steepness effect by the applied voltage. Two main phenomena are considered: charge formation through impact ionization, and space charge reduction through attachment. The former phenomenon is associated with Townsend's first ionization coefficient α, the latter with the attachment coefficient η. The inception criterion, derived from corona onset experimental data, is expressed by the following analytical condition: In (3), K(r) represents the number of free electrons at distance r from the conductor axis. Coefficients α = α(p, E) and η = η(p, E) have to be intended as functions of the air pressure p and of the instantaneous electric field E(r), assumed radial; δ 0 denotes the thickness of the corona layer, to be intended as the distance from the conductor surface at which free charges produced by impact ionization and those canceled by attachment are equal, i.e., the distance where the following condition is verified: Expressions for α and η can be found in [40], and read as follows: where E is expressed in kV/cm and p in mmHg. Both the coefficients α and η are computed considering that the space charge influence on the electric field is negligible before inception occurs. This assumption made, starting from the configuration under study, i.e., coaxial arrangement, or conductor above a perfect electric conductor (PEC) surface, and the conductor features (i.e., single conductor of radius r 0 , or stranded conductor), the electric field E(r) to be plugged into (5) and (6) may be found from the known applied voltage. Hence, relations (5) and (6), along with the trend of the electric field surrounding the conductor, have to be taken into account simultaneously to compute the inception value of E inc through the solution of Equation (3).

Inception Field by Mikropoulos and Zagkanas
A modification to the general formula (2) for the inception electric field (in kV/cm) is proposed by Mikropoulos and Zagkanas in a coaxial arrangement [37]: where K s = 0.42 · sign dv dt 10/90 dv dt 10/90 0.345 (8) where the voltage time derivative at the right-hand side of (8) is the steepness of the applied voltage in kV/µs, assuming its waveform to increase linearly between 10% and 90% of its peak value. In fact, experimental data from [37] are in good agreement with Peek's inception field only in the case of slow-front voltage sources (up to 10 kV/µs); when different impulse voltages are applied to a reference conductor, electric field onset values show a dependence on the impulse steepness, which is neglected in expressions based on Peek's formula. Hence, the steepness correction factor (8) is introduced to take into account the increase in E inc with steeper fronts of the applied voltage.

Role of the Irregularity Factor: Stranded Conductors and Bundles
In the current section, some more in-depth information will be given about the irregularity factor m, when dealing with stranded conductors and bundles.
This factor allows to evaluate the corona inception field of any conductor, starting from V inc computed for a reference conductor with the same outer diameter and smooth cylindrical surface. Indeed, m may include also several physical or geometrical features of the conductor generating irregularity, i.e., grease, dirtiness or dust, and scratches. Several tests were conducted by Baker et al. [41], showing the influence of surface irregularities on corona development in heavy rain conditions, or for wet surface/foggy weather; differences were observed between new (as from manufacturer) and aged conductors, the former showing more pronounced corona activity with heavy rain [41]. This is attributed to the greased surface, which causes water to lay in droplets, enhancing locally the electric field. Aged conductors, instead, show hydrophilic behavior, letting the water flow toward their underside. Hence, factor m might take a different value, depending on the conductor aging process. Instead, the effect of the aging process of conductors operating in polluted areas, i.e., the increased roughness of their surface, is investigated in [42].

Stranded Conductors
When dealing with stranded conductors, the diameter of the external strands, and their arrangement affect the value of the inception voltage; it is usually taken into account through the stranding factor, defined as the ratio between V inc computed for the stranded conductor, and V inc of a smooth conductor sharing the same outer diameter [43].
Experimental data on the stranding factors with stranding ratios (i.e., ratio between the diameter of any external strand and the outer diameter of the conductor) in the range 0.05 ÷ 0.25 and different conductors diameters were collected by Stone [43]. As expected, the found values converge to 1 for the smallest stranding ratios; for the growing stranding ratio, the stranding factor is found to decay faster when dealing with larger conductors, suggesting that the inception of corona may not be related only to the maximum value of the electric field on the surface of the strands.
Hence, Yamazaki and Olsen [11] elaborated the inset criterion, reviewed in Section 3.3.2, based on the electric field distribution in the proximity of the conductor.
This criterion justifies the values assumed by the stranding factor, i.e., the reduced V inc at which the corona inception is observed; the stranded surface results in a different electric field distribution, and in the enhancement of its maximum value with respect to the case of a smooth conductor for a given voltage.
Even though we have referred to the potential, since it is often used in numerical codes as a threshold value for the simulation of the corona inception, the electric field distribution is the primary cause of corona. Bousiou et al. [38] showed the trend of the electric field magnitude at increasing distance from the conductor outer radius r 0 , for different stranding ratios, showing that, for r → r 0 , the electric field for stranded conductors tends to a value approximately equal to 1.4 times E inc for the corresponding smooth conductor, independently from the number of external strands; however, as expected, the trend of E(r) approximates that of smooth conductors as the number of external strands increases. Furthermore, at equal applied voltage, the electric field holds higher values with a decreasing number of strands, also at larger distance from the conductor axis [44]; it is reasonable that the inception depends on the distribution of the electric field and not only on its maximum value; thus, lower values of E inc found for higher stranding ratios [38] are justified.
When accounting for the stranding, misleading interpretation may occur about the role played by the surface regularity factor m in (1). In fact, only V inc , as a function of the stranding factor, effectively decreases, due to stranding. Hence, the value of E inc , with the irregularity factor m, corresponds to the inception field of an equivalent smooth conductor sharing the same inception voltage with the original stranded one, but not the same inception field at the surface.

Bundles
The technical solution of bundle conductors, for high-voltage transmission applications, is well-known for increasing lines loadability and reducing the occurrence of corona; this latter effect has been largely investigated in the literature, especially with the increasing spread of EHV and UHV lines.
Peek formula (1) does not apply to bundle configurations; research has been conducted to establish a suitable corona inception criterion. Waiving the approximation of uniform radial distribution of the electric field, it is necessary to refer to the maximum value E max attained by the electric field in the bundle proximity, which can be computed through the charge simulation method [45], or approximately for the coaxial configuration [46] and for bundles over a conductive plane [47] through the general formula as follows: where n is the number of subconductors of radii r 0 of the bundle; r b is the bundle radius; v is the voltage applied to the subconductors of the bundle; d denotes the radius of the external conductive cylinder in the coaxial configuration, or the height of the bundle over the conductive plane; and r eq = n r 0 r n−1 1/n is the equivalent radius of the bundle [47]. Indeed, the exact solution of the electrostatic problem of charged cylinders over a conductive plane, also accounting for proximity effects, was solved through conformal mapping in [48].
Data presented by Maruvada et al. in [36] show a slight increase in the corona inception field for bundles with respect to single conductors; indeed, in the former case, charged ions and electrons might need an increased electric field to properly start an avalanche, due to the asymmetric radial distribution. Nevertheless, the authors of [46] observed that V inc for bundle configuration may be estimated satisfactorily by the applied voltage, causing E max to be equal to the inception field of any of the subconductors, as if it was alone on the ground.
The influence of the number of subconductors and bundle radius was assessed in [46]. In particular, while increasing the number of subconductors always results in a reduced maximum electric field, and hence, in a higher inception voltage, the dependence on the bundles radius is more complex. For a fixed number of subconductors, an optimum bundle radius can be identified, which maximizes the inception voltage (associated with the configuration generating the weakest E max ).

Models for Corona Simulation
In the literature, several models are available for the simulation of the corona phenomenon, in order to reproduce the multi-valued hysteretic relation between the charge and the voltage. The most relevant models included here were classified as in Table 1, according to the different approaches adopted for their formulation. It should be noted that the following models are the preferred choice when the effect of corona discharge needs to be considered in TL studies with respect to multi-species discharge models. In fact, despite being more accurate and adherent to the physical reality, the latter models require larger computation times, usually by means of the finite element method (FEM) with a suitable iterative solver for the solution of non-linear charge transport equations. The diagram in Figure 2 shows, at a glance, the models here reviewed. We refer to the first subset as "physics-based" models. These models use a simplified description of the physical mechanisms underlying the corona phenomenon, which would require a coupled electro-fluid dynamic plasma model to be fully determined [49,50]. Complexity arising is the main drawback of such detailed models, which conflict with the aim of formulating simple tools suitable for engineering purposes. Nonetheless, the main advantage of their application consists in their general applicability, through non-linear relations accounting for all the key physical aspects from a macroscopic point of view. In these models, corona losses derive directly from the area described by the q-v hysteretic curve and hence, no additional conductance is introduced.
Starting from measured q-v curves in known configurations, "empirical" models attempt to define q or C (p.u.l. quantities) as a function of the applied voltage in order to approximate the measured data. The most remarkable advantage is the easiness of their implementation; as a downside, they all require a careful choice of the values of a set of input parameters in order to obtain q-v curves close to those resulting from measurements.
Finally, availability of transient-programs (e.g., ATP/EMTP, PSCAD) has fostered the development of "circuit-based" corona models to be easily implemented in conjunction with widely tested models of multiconductor TLs and other equipment. Circuit-oriented models introduce shunt circuital elements to reproduce the propagation effects of corona phenomenon, i.e., attenuation and time delay. The additional lumped components may be generally expressed as the product of some p.u.l. voltage-dependent electrical parameters (i.e., C , G ) and the spacing ∆ adopted for the insertion of these equivalent additional shunt circuits. Their circuit structure reveals an advantage for the implementation in commercial transients programs, through simple construction of equivalent circuits to be plugged to the line under study. As a drawback, these models compute attenuation and delay due to corona only at chosen voltage nodes, simulating a phenomenon which is actually distributed as a lumped equivalent one.
We note that both the empirical and the circuit-based models may lead to the definition/computation of a voltage-dependent, i.e., time-varying, capacitance; strictly speaking, if the conductor capacitance is assumed to be non-constant, the associated p.u.l. capacitive current i c should be computed as follows: In Equation (10), C dyn is known as the p.u.l. dynamic capacitance of the conductor and is defined as the p.u.l. charge derivative over the voltage [13]. In the absence of corona, C dyn is constant and equal to the conductor geometric capacitance. The dynamic capacitance approach is frequently used in the literature to account for the effect of distributed corona on propagation along high-voltage transmission lines, deriving the expression of C dyn from one of the available corona models. However, the general validity of this method is questionable, especially if the adopted corona model assumes a delay between the charge and voltage. In fact, when dealing with transmission line models, the assumption of quasi-transverse electromagnetic propagation implies that static-like expressions for the fields are valid in the transverse plane with respect to the direction of propagation; hence, the p.u.l. line capacitance should account for the static, i.e., instantaneous, relation between p.u.l. charge and voltage. Including the delayed dependence of the charge over the voltage in C dyn would be in contrast with the assumptions underlying the transmission line theory, leading to erroneous results. Alternative approaches with distributed voltage-controlled current generators may be also employed [51].

Model
Ref.

FEM Models
Although marginal to the present review, it is appropriate to briefly revise the different electro-hydrodynamic models [72] proposed for simulating corona discharge, which are usually solved by means of the FEM [73,74]. The models are based on the classical engineering approach, which replaces the microscopic point of view with a continuum modeling. They lead to a system of fully coupled nonlinear partial differential equations that can be conveniently crafted to account for the specific transport properties and constitutive relations of the corona discharge processes. The simplest and most idealized model is the single-species approach, which considers only one species of ionic charges, whose density is n [1/m 3 ], moving in the electric field with a constant mobility µ. The electrostatic potential ϕ = −∇e, where e is the electric field, and the transport current of the charges n are governed by the following system of equations, comprising the Poisson's equation and the continuity equation under the so-called drift-diffusion approximation: where D is the diffusion coefficient (usually the diffusion is neglected) and U = G − R is the net difference between the generation G and the recombination R rates for the charge carriers (q is the absolute value of the electron charge). More realistically, the two-species model [75] assumes interactions among positive p and negative n ions, and the threespecies [76,77] model assumes also the electron attachment to neutral molecules m. In the latter case, the equations read as follows: Approaches including several species (seven and more) have been proposed in the literature [78]. Additionally, in [75] the effect of the wind velocity w has been observed to be of the same magnitude of the drift term.
Often, these approaches reduce ionization phenomena, localised in a limited layer close to the conductor surface [38], to simplified boundary conditions ensured by a suitable space charge distribution. For instance, while values of the electrodes voltages are assumed to be known (Dirichlet's conditions), additional space charge may be computed and injected in the surrounding area in order for the electric field at the conductor surface to hold a constant value equal to E inc for the entire corona development. This frequently adopted hypothesis is also known as Kaptzov's assumption [79].
Since the potential influences the space charge density and vice versa, the system of differential equations is non-linear, and algorithms exploited for the simulation of corona rely on iterative solvers to compute the sought solutions. The system of differential equations is usually solved by means of the FEM [80] under DC and AC steady states [75,80] or transient conditions [81,82]. Well-known multi-physics commercial software is available on the market for this purpose. Some assumptions are commonly made as follows: the potential on the electrodes surface is known (Dirichlet-type condition); the electric field on the surface of the conductor experiencing corona is approximately constant to its known onset value, in accordance with Kaptzov's assumption (Neumann-type condition); diffusion is neglected; the ionization area is neglected due to its limited thickness [73][74][75][76][77][78][79][80][81][82][83]; and the surface of the coronating conductor is considered the only source of new space charge.
Referring the interested reader to the appropriate literature [74], FEM is briefly summarized by four fundamental steps: discretizing the computational domain into a specific number of sub-regions called elements; deriving matrix equations for each element approximating the physical unknowns with scalar (most often polynomials) and/or vector basis functions over each element; assembling the element equations into a single matrix (usually of sparse nature) and, finally, solving the obtained system of equations with direct or iterative (parallel) solvers. If the system of equations is non-linear, Newton-Raphson-like methods are used. FEM on triangular and tetrahedral meshes is attractive for several reasons: it exhibits second order accuracy for the diffusion equation and allows a high degree of geometrical flexibility with a relatively small number of spatial unknowns. Yet, attention is needed on possible spurious solutions, which arise in the diffusion problem and derive from the geometrical scale and parameter differences. Hence, stabilization schemes are mandatory.

Correia de Barros' Model
An accurate analysis of the corona effect was proposed by Correia de Barros, originally on a single cylindrical conductor [53], and later extended to multiconductor systems [52]. Further details may be found in some earlier works, such as [54].
The model assumes that a time delay τ exists between the time t inc at which the critical inception field E inc is reached and the time at which the discharge is visibly triggered in the vicinity of the conductor surface. This time delay is composed of two different time lags, i.e., τ = τ st + τ sp . The first term τ st is the statistical time necessary for the formation of seed electrons by detachment (the concentration of free electrons in a standard atmospheric condition, related to the electronegativity of oxygen in the air, would be too low to start an avalanche); the second term τ sp is the critical time required for the formation of the space charge. The former τ st depends on the gap size (i.e., on the conductor height), the preexisting amount of available charged particles and the magnitude of the overvoltage with respect to V inc [26]; lower values of τ st were observed for larger applied voltages, tending toward a minimum value of about 0.25 µs [84].
The dynamics of formation and radial spread of the space charge around the conductor is modeled by considering a growing radius r inj (t) of the injection layer, as well as an exponential growth of the generated space-charge density. The key parameter r inj (t) permits here, as in other models, to separate two distinct phenomenological regions, i.e., the injection discharge region (also denoted as the glow region or ionization region) nearby the emitter conductor, where there is a source term for charge production caused by electronic impact and attachment, and the drift region away from it, where charge concentrations are low and charges are mainly electro-convected by the ruling electric field. The space surrounding the conductor is subdivided in M cylindrical layers with radial thickness s i and outer boundary equal to r i , with i = 1, 2 . . . M as depicted in Figure 3. The configuration is studied through an equivalent coaxial one, in which the conductor holds the same capacitance C geo ; hence, r M = r 0 exp 2π 0 /C geo = 2h [35]. Assuming a constant streamers average velocity v st [85], the instantaneous radius r inj (t) is given at any time t > t inc + τ by the following: where r s is the outer limit for the injection area. The larger is r s ; the wider is the resulting q-v loop; and the higher is the energy dissipated per cycle. From experimental evidences on transients events, the discharge area associated with negative polarity corona shows minor extension and typically displays a narrower loop [62]; this can be taken into account through the adoption of a lower value of r s under negative surges. The algorithm consists of two main steps. In the first step, we evaluate the electric field e 0 (t) at r 0 , employing the currently applied voltage v(t) and the space charge of the previous time step. Applying the electrostatics Green's reciprocity theorem [86], under the assumptions of linearity and ground permittivity much greater than that of the medium surrounding the conductor, the following relations hold: where and where q e = 1602 × 10 −19 C is the elementary electron charge in absolute value and q ind (t) is the p.u.l. charge induced in the conductor by the total space-charge q sp (t); k i is a geometric coefficient with the dimensions of an area. In (15), p i and n i are the volumetric densities of positive and negative charged particles in the i-th layer (expressed in 1/m 3 ).
As discussed in detail in [87], the surface electric field e 0 (t) is found to remain constant during discharge and close to the onset field strength E inc : this is known as Kaptzov's assumption [79], widely used in simulations to obtain quantitative estimates of the glow corona with moderate computational effort; however, it is strictly valid for a non-stationary corona only when the background electric field changes slowly [88].
Depending on the polarity of the electric field e 0 , when | e 0 (t) |> E inc , the ionization phenomenon is simulated through the injection of a new generated p.u.l. space charge of the same polarity, whose volumetric density dρ(t) (charged particles injected per cubic meter) is assumed to be uniform in the injection area r 0 < r < r inj (t) and given by the following: where τ 0 in the exponential factor accounts for the dynamics of space charge injection.
In (17), S inj (t) = π r 2 inj (t) − r 2 0 is the injection area and k inj (t) corresponds to the integral in (16) performed from r 0 to r inj (t).
The electric field e i (t) at the separation surfaces between the M layers is found recursively from e i−1 (t) as the following: where s i and r i = r i +r i−1 2 are the thickness and the average radius of the i-th layer, respectively. In the second step of the algorithm, the space-charge densities p i (t) and n i (t) are computed from the electric fields of the previous step through the solution of the following system of 2M non-linear differential drift equations (the time dependence is dropped for clarity): dp i dt with i = 1, 2, . . . M. These relations assume a constant mobility µ p and µ n of the positive and negative charges, respectively, as well as their recombination through a constant recombination coefficient R. Typical values for µ p , µ n , and R are reported in Table 2.
In (19), α = i and β = i − 1 under positive electric field e i , whereas α = i + 1 and β = i under negative electric field e i . Additionally, Dirichlet boundary conditions are assumed, i.e., p 0 = p M+1 = n 0 = n M+1 = 0. Once the space-charge densities p i (t) and n i (t) are known, the total space-charge q sp (t) in (15) is computed as the following: Finally, the total p.u.l. charge q (t) needed for plotting the q-v hysteretic loop is computed as q (t) = C geo v(t) + q ind (t) + q sp (t). From the numerical derivative dq /dv, the dynamic capacitance C dyn is computed.
The main parameters required for the implementation of the model are synthesized in Table 2; the range of variability given for τ s , v st , and r s allows to select these parameters to fit the experimental data. In Figure 4, the trend of the space-charge density ρ sp is shown, when a voltage source is applied to a conductor of radius r 0 = 1.32 cm at height h = 7.5 m over a perfectly conducting ground. The waveform of the voltage source is given by the following expression: where V max = 2V inc , V inc = 263.25 kV, η ≈ 0.78, T 1 = 167 µs and T 2 = 2648 µs, time to peak T p = 493 µs, rise time (between 10% and 90% of the voltage peak value) T r = 250 µs, and time to half value T h = 2500 µs. It should be noted that the parameters T 1 and T 2 in (21) differ from the time to peak and the time to half value of the resulting waveform; hence, values of T 1 and T 2 were chosen to compute results in Figures 4 and 5 corresponding to a standard 250/2500 µs switching impulse.

Malik's Model
In this model [56], as long as the voltage v(t) is greater than V inc , the dynamic behavior of C dyn is simulated through an apparent increase in the conductor radius r 0 , which is replaced by the radius r c (t) ≥ r 0 . This latter corresponds to the external boundary of the whole space-charge around the conductor and it must not be confused with r inj of the previous model, defining only the ionization area.
The model is developed for a single conductor above a conducting ground plane, at height h, and assumes a time delay τ (ranging between 0.1 µs and 0.5 µs) in the formation of corona charge with respect to the instantaneous value of the voltage applied to the conductor. The time-dependent p.u.l. charge q (t) is given by the following: where the radius r c (t) in the presence of corona can be determined by solving the following non-linear equation at each time instant t: In (22) and (23), α < 1 is a multiplicative factor taking into account the reduction in the electric field in the corona area after the inception. The main assumption of the model is to consider the electric field inside the corona sheath (r 0 ≤ r ≤ r c ) constant and equal to αE inc ; consequently, the surface electric field e 0 (t) (oppositely to the Kaptzov's assumption) undergoes an abrupt discontinuity when the corona discharge starts, which can result in a discontinuity in the q-v relation. Comparison with (23) shows that the parameter α introduces a step discontinuity between r 0 and r c (t) at t inc , where the lower is α, the larger is the step discontinuity, and the wider is the resulting corona hysteretic loop. Along with τ, this is the main parameter to be tuned for optimal fitting of measured results (Table 3). Table 3. Typical values of the input parameters required by Malik's model.

Parameter
Value From (22), the dynamic capacitance C dyn is computed as dq/dv; when the maximum charge is reached and the voltage v(t) decreases, the dynamic capacitance is considered constant and equal to the geometric value C geo .

Cooray's Model
The theory formulated by Cooray [59] is based on a coaxial geometry. Herein, the original formulation is extended to the practical configuration of a cylindrical conductor over a PEC plane in order to obtain expressions suitable for comparison with other models, but still assuming a radially symmetric electric field.
Cooray describes the physics of the corona phenomenon, identifying four different stages. In the first stage, the voltage v(t) applied to the conductor increases progressively; hence, a proportional relation is assumed between q (t) and v(t), with C geo being the proportionality coefficient. As the voltage reaches the inception threshold V inc , the discharge starts as well as the second stage of the model. For a positive (negative) voltage surge, a positive (negative) spatial charge progressively settles around the conductor; this phenomenon is taken into account through a time-dependent increasing radius r c (t). The model is based on some relevant assumptions on the physical development of the discharge phenomenon. Denoting with r c (t) the external radius of the corona charge distribution, the corresponding electric field e(t, r c ) is forced to the value E c , which depends on the atmospheric conditions, the conductor characteristics, and the applied voltage polarity. A time-dependent expression is adopted for the electric field at the conductor surface e 0 (t): it decays exponentially from its inception value E inc to E c , which is assumed to be the minimum field value to guarantee streamers propagation. The expression reads as follows: In (24), τ d is the time constant defining the electric field decay, and t inc denotes the corona inception time. The distribution of the corona space-charge q sp + (without loss of generality, we refer hereafter to positive corona) is assumed to be dependent on the inverse of the radial distance r from the conductor central axis (i.e., of the type ρ sp + r −1 ).
It is worth noting that, unlike previous models, any time delay in the spatial charge formation is neglected. Based on these assumptions, the following relations can be derived from trivial electrostatics considerations, given h r 0 and angular symmetry: where the total charge q (t) is the sum of the charge on the conductor surface q 0 (t) and the corona spatial charge q sp + (t): where 2h r 0 , r c (t) (as in the case of transmission lines conductors), and ρ sp + (t)r −1 represents the radial distribution assumed for the p.u.l. space charge density in C/m 3 .
From (26) and (27), we readily obtain the following: Considering both the conductor and its perfect image located at depth h below the ground surface, the instantaneous voltage is given by the following: where the time-dependence of r c (t) is dropped for conciseness. Inserting (28) into (29), the non-linear Equation (29) can be solved by an iterative method in the unknown r c (t); then, the total charge q (t) is obtained by (25) and (26). Table 4. Values of the input parameters required by Cooray's model.

Parameter Value
After the inception, the ionized area is supposed to expand as long as the sign of the voltage derivative is positive. The third stage of the model begins when the maximum voltage is reached and the voltage derivative changes sign. The model assumes that the corona sheath radius and the charge density are fixed to their maximum values r M c and ρ M sp + , respectively, due to the slow mobility of the space charges. The voltage decrease is associated with an initial progressive reduction and a subsequent change in the sign of q 0 (t), until e 0 (t) = E ib , i.e., the electric field causing the inception of back-corona and the development of negative space-charges.
The back-corona phenomenon is the fourth stage of the model. A second ionizing process begins in the radial direction from r 0 , neutralizing progressively the previously settled positive charge and setting new negative space-charge q sp − (t). This new negative charge is assumed to be distributed in the area between r 0 and an increasing radius r cb (t), with the same distribution as in (27). The electric field is assumed constant and equal to E ib in the back-corona area, i.e., r 0 ≤ r ≤ r cb (t). Hence, the negative charge density can be readily obtained as ρ sp − = 0 E ib . Finally, the voltage of the conductor may be expressed as follows: Equation (30) can be solved to find r cb (t); once the radius is known, the total charge is computed as q (t) = q 0 (t) + q sp + (t) + q sp − (t). Values of the input parameters required by Cooray's model are summarised in table 4.
A last remark concerns the transition between the first and the second stage. Let us focus exactly on the corona inception at time t inc ; if we assume a continuous function r c (t) at time t inc , i.e., r c (t) → r 0 for t → t inc , from (26), we may write the following: Therefore, we should conclude that a non-zero negative space-charge would exist in an ionized corona of infinitesimal thickness. This limit demonstrates how the model produces a step discontinuity between r 0 and r c (t) when the corona discharge begins, as in Malik's model.

Inoue's Model
The expression proposed for the dynamic capacitance [61] is the following: with the ratio [v(t) − V inc ] α 1 −1 /v(t), considered to be dimensionless, and the following: In (32), α ranges in the interval 2 ÷ 2.1, and in (33), the value of σ κ is in the interval 100 ÷ 450 pF/m. Both the values have to be optimized to fit the experimental data; to this aim, the branch of the measured q-v loop ranging between V inc and the maximum voltage value V max may be fitted through two curves with different values of κ and α for [61].
The model does not make any distinction between the cases of positive and negative corona discharge. The capacitance C dyn (t) is continuous with C geo at the inception instant, when v(t) = V inc ; no discontinuity due to the formation of the space charge is predicted.

Gary's Model
According to the model presented by Gary [62], the p.u.l. charge after the corona inception and for dv/dt > 0 is given by the following: and the dynamic capacitance of a single conductor under corona can be computed as follows: where the coefficient B is given by the following experimental formula, which distinguishes the cases of impulses of different polarities: The model predicts a step discontinuity between C dyn (t) and C geo at the inception of corona discharge (C dyn (t) → BC geo for v(t) → V inc ), which depends on the parameter B in (36), i.e., on the surge polarity and the conductor radius r 0 .

Podporkin and Sivaev's Model
Sivaev and Podporkin originally proposed in [63] an expression for the p.u.l. charge q (t) under corona as follows: From the definition given in (10) for the dynamic capacitance, C dyn is computed as follows: where κ 1 = 1.17, κ 2 = 0.87, and a is equal to 0.08 and 0.036 for positive and negative impulses, respectively. The expression in (38) is valid for a single conductor, while a modified expression suitable for bundled conductors can be found in [63]. At the inception voltage, the method predicts a discontinuity in the transition from geometric to dynamic capacitance: in fact,

Suliciu's Model
The corona phenomenon is simulated assuming that the current i c drained to ground by any elemental transverse section in A/m, defined in (10), may be expressed as the sum of two contributes [64]: The p.u.l. current i sp , which corresponds to the time derivative of the space charge only, can be computed as follows: and C 2 > C 1 > C geo , K 1 > 0 and K 2 > 0; a common choice forṽ 1 > 0 andṽ 2 > 0 is v 1 =ṽ 2 = V inc . Parameters C 2 , C 1 are on the order of C geo in pF/m, while K 1 and K 2 range between less than 1 kHz and 5 MHz; however, the values of these parameters should be adapted in order to properly fit measured data. Integrating (40) over time to obtain q sp (t), the total charge may be computed through the following relation: Finally, the dynamic capacitance C dyn can be found by performing the numeric derivative of the p.u.l. charge in (42) with respect to the applied voltage v(t). Huang et al. extended the original model in order to take into account an applied oscillating voltage source [65].

CIGRÉ Model
An alternative approach for corona modeling may be found in [66] in which a linear charge-voltage relation is assumed at voltages lower than V inc , and for the descending branch of the q-v loop (i.e., for dv/dt < 0), while the following expression is considered for the total charge during corona development: In (43), K is a constant to be determined from the fitting of experimental data, or graphically as a function of the conductor diameter [66]; C I may be intended as a step discontinuity assumed for the capacitance at the corona inception (to be derived from fitting of measured data); in fact, for v(t) → V inc , the right-hand side of (43) reduces to the product C geo + C I v(t). However, the value of C I should not exceed 1 pF/m [66].
From the definition given for the dynamic capacitance, C dyn can be derived as follows:

Hara and Umoto's Model
Several models are based on the assumption that the p.u.l. corona losses p can be computed through the following quadratic frequency-dependent relation originally proposed by Peek, and expressed in W/m [4]: where f is the frequency in Hz, v(t) and V inc are in V, δ is the air relative density, and κ is a corona loss constant, expressed in s/(Ωm), found from the analysis of experimental data.
Starting from (45), Hara and Umoto proposed to include the effect of corona in a finite-difference time-domain (FDTD) updating scheme by connecting in parallel two lumped non-linear voltage-dependent components along the line when the electric field e 0 (t) reaches the inception value E inc . The equivalent circuit is represented in Figure 6a; additional lumped shunt branches include the capacitance ∆C, modeling corona space-charge effect, and a conductance G C , accounting for additional corona losses. The expression for the currents drained by the branches, i ∆C and i G respectively, were originally proposed in [68], and then recalled by Lee [69]: with The factors κ C and κ G may be expressed as follows: where σ C is approximately in the range 15 ÷ 35 F/m, and σ G (corona loss constant) is in the interval 5 × 10 6 ÷ 20 × 10 6 S/m, depending on the configuration under study. For sufficient accuracy, the distance between voltage nodes equipped with additional lumped corona branches should be less than ∆ = 70 m; indeed, a further increment in ∆ would result in noticeable deviation from the benchmark experimental data [69].

Motoyama and Ametani's Model
In order to avoid the non-linear components present in Hara and Umoto's model, Motoyama and Ametani [70] proposed a further development. The expressions given in (47) for the lumped equivalent capacitance and conductance are replaced by three-levels piece-wise constant functions of the voltage v(t), which are valid for dv/dt > 0: and where The parameters κ C and κ G are given by expressions (48), in which σ C ranges approximately between 2 ÷ 33 F/m, and σ G is in the interval 0.5 × 10 6 ÷ 13 × 10 6 S/m. The proposed model may be easily embedded in any transients program, plugging the equivalent circuit in Figure 6b to selected voltage nodes; the model is claimed to be more practical with respect to Lee's model: it should allow to plug lumped circuits for corona at larger intervals, even just at the poles for ∆ in the range 350 ÷ 450 m, still holding accuracy [70].

Maccioni-Araneo et al. Model
The model refers to a single conductor, and introduces a couple of lumped additional shunt branches [71] as shown in Figure 6c: the first branch consists in the voltage-controlled current generator i ∆C (t) simulating the increment in the capacitive current; the second branch is represented by a non-linear conductance G ∆C (t), accounting for corona losses. As the voltage v(t) reaches V inc , the voltage-controlled current generator is turned on, fictitiously reproducing the increase in capacitance ∆C(t) due to corona. As the voltage starts to decrease, it is turned off so that ∆C(t) = 0 on the descending branch of the q-v curve, and the conductor charge depends linearly on the voltage through C geo ∆ .
Following the theory of the physics-based models, when dv/dt > 0, ∆C(t) is computed as the following: where the time dependent radius r c (t) is computed through the following non-linear equation: The current generator i ∆C (t) = d[∆C(t)v(t)]/dt is given by the following: Equations (52)- (54) can be easily inserted into any FDTD scheme. The value of the lumped equivalent conductance G C (v) is computed according to (47b) [69] for v(t) ≥ V inc and dv/dt > 0. Once the voltage applied to the conductor has reached its maximum value V max and begins to decrease, the conductance G C is kept connected to the line, according to [89]; however, the value of G C is no longer voltage-dependent, but is assumed to be constant and equal to the value attained by (47b) when the corona radius has reached its maximum value r = r max .
With reference to (53), the main simplifying hypothesis adopted by this model consists in assuming the corona sheath to be conductive enough so that the voltage drop across the ionized area may be disregarded; hence, the potential v(t) is intended to be measured at the outer surface of the space-charge region, not at r 0 on the conductor surface. Table 5 is provided to compare the models reviewed here in terms of the number of required input parameters, and applicability with respect to the voltage source applied to the conductor under test. As to their validity with respect to the voltage source, all the models should be exploited when a unipolar waveform is applied to the conductor, except for the model by Correia de Barros, which is suitable for studies involving any voltage waveform, both impulsive and AC steady-state (due to the modeling of the drift of negative and positive charges). In this section, curves are provided to address the influence of the parameters required by each model on the shape of the final q-v curves. The reference voltage impulse displayed in Figure 7 is applied to a reference conductor (whose characteristics are summarized in Table 6):

Results: Comparison of the Models Responses to the Variability of Input Parameters
where τ 1 = 10 µs, τ 2 = 75 µs, and η = 0.636 is the amplitude correction factor. To this aim, the voltage maximum value is chosen to be equal to V max = 3V inc , V inc for the reference conductor being computed, according to Peek's formula (1) since the validity of the formula proposed by Mikropoulos and Zagkanas [37] is limited to coaxial configurations, and the criterion proposed by Olsen and Yamazaki applies to power frequencies.  Table 6. Parameters of the reference configuration in the inset of Figure 7. In the following sections, we address the variability of the hysteretic q-v loops computed by the different approaches with respect to the input parameters required by the models. In particular, for a model depending on n p input parameters, an n p -dimensional finite region of space is defined, where the n p parameters may vary in well-defined ranges of admissible values; q-v curves are displayed corresponding to one hundred points of this space, i.e., to one hundred different combinations of the input parameters, the values of which are randomly extracted from a uniform distribution in their variability range.
The results of Malik's model are in Figure 8b. The input parameters required by the model, i.e., τ and α, range between 0.1 ÷ 0.5 µs and 0.5 ÷ 0.9, respectively. The model does not show significant dependence on the input parameters; their impact may be noticed at the inception. In fact, the different time delays τ result in slight variations in the location of the discontinuity in the q-v loops at the inception. Furthermore, lower values of α result in a faster growth of the equivalent external radius; hence, the slope of the rising branch in Figure 8b is more pronounced, enlarging the total loop area. Charge-voltage curves computed by Cooray's model are shown in Figure 8c. The variability ranges assumed for the input parameters are: The slope of the ascending branch of the loops, after the inception, is related to E c : values of E c closer to (much smaller than) E inc result in a slow (fast) growth of the charge with respect to the voltage. If the electric field on the conductor surface reaches the value E ib at the voltage tail, the charge decay will be faster, with lower values of residual p.u.l. charge. In fact, the descending branch will deviate from the expected linear decay with slope equal to C geo (see Figure 8d).
As to the shape of the hysteretic curves, the models result in different shapes; with reference to Malik's model, the loops are tight and of definite shape, predicting the maximum p.u.l. charge values to be lower with respect to results by Cooray and Correia de Barros, which, on the contrary, are in good agreement.
Comparison of the descending branches of the curves reveals some differences. The constant value of the p.u.l. capacitance assumed by Malik, equal to C geo , results in the constant slope of the corresponding branch. Instead, Cooray's curves show more pronounced slopes for dv/dt < 0, due to the progressive accumulation of negative p.u.l. charge on the conductor surface, and to the subsequent inception of back corona (depending on the value of E ib ).
Correia de Barros' results are computed for a larger time-window, up to 1 ms, in order to properly show the curvature of the descending branches for v(t) → 0. If the other models present a residual p.u.l. charge caused by the applied voltage amplitude and time derivative approaching zero, q (t) vanishes in Figure 8a for v(t) → 0 and t → ∞, due to the gradual drift of residual charged particles surrounding the conductor. Figure 9 collects the responses to the applied voltage source computed by means of the empirical models in Section 4.3. Sensitivity studies are conducted by choosing some variability ranges for the parameters required by the models (which should be set by fitting experimental data), from analogy with values adopted in previous applications in the literature (if available), or arbitrarily. Inoue's model gives rise to the widest range of variability of the output curves; hence, its predictive capability is questionable, and comparison with measured data is essential for the suitable tuning of the input parameters and computation of reliable results. The curves in Figure 9a are obtained for α ∈ [2, 2.1] and σ κ ∈ [100, 450] pF/m (similar values are adopted in [61] for the first branch of the q-v loop after the inception). The product σ κ · α determines the slope of the non-linear sections.

Sensitivity Curves for Empirical Models
Gary's approach requires the choice of the parameter B to fit experimental data (see Equation (36) for original values by Gary). In order to address the model dependence on B, curves in Figure 9b were computed with B ∈ 1.23 ÷ 1.55, i.e., in the range defined by values in Equation (36), chosen as limiting values of the interval. Wider loops correspond to larger values of B. Figure 9c presents the curves computed through the approach by Sivaev and Podporkin. Since κ 1 = 1.17 and κ 2 = 0.87 should be of practical use for 0.5 < r 0 < 3 cm and 10 < h < 30 m [63], here, we adopt limited (arbitrary and heuristic) ranges of variability to account for the different height of the conductor h = 7.5 m < 10 m): κ 1 ∈ [1.1, 1.2], κ 2 ∈ [0.8, 0.9], while a is in the interval 0.036 ÷ 0.08.
As for Suliciu's model, the following set of input parameters was chosen: C 1 ∈ C geo , 2C geo , C 2 ∈ C geo , 4C geo ensuring C 2 > C 1 , K 1 ∈ [0.5, 5] MHz, and K 2 ∈ [0. 5,5] kHz. The output curves show a characteristic shape; in particular, the maximum p.u.l. charge value may not correspond to the maximum applied voltage. This is due to the fact that the capacitive current derived toward the ground depends not only on the voltage time-derivative, but also on the voltage instantaneous value (Equation (41)). Sensitivity curves computed by the Cigré model with respect to C I , K are shown in Figure 9e. Considering the values of K in [66] corresponding to the diameter of the reference conductor, q-v loops are displayed for K ∈ 2.5 × 10 −6 , 5 × 10 −6 pF/(V·m), while the value of C I < 1 pF/m [66] is chosen between 0.1÷1 pF/m. Growing values of C I introduce harsher step-like discontinuities at the inception voltage, i.e., for v ≈ 300 kV; instead, larger values of K enhance the rate of rise of the non-linear branch of the curves.

Sensitivity Curves for Circuit-Based Models
In Sections 4.4.1 and 4.4.2, the parameter denoted with σ C is responsible for variation of the capacitance computed by the corresponding models. In Figure 10b,c, we present the curves computed by the models of Hara-Umoto and Motoyama-Ametani, respectively, with σ C ranging in the interval 15 ÷ 35 F/m (from values adopted in [69], in good agreement with [70]). Both models are strongly dependent on the input parameter, as may be deduced by the noticeable variability in the displayed curves. Greater values of σ C correspond to larger q-v loops. The hysteretic curves in Figure 10b are wider than the curves in Figure 10a. This is due to the simplified approach proposed by Motoyama and Ametani (see Section 4.4.2). In fact, for the voltage intervals V inc ≤ v(t) < 2V inc and 2V inc ≤ v(t) < 3V inc , the piece-wise constant capacitance defined by Motoyama and Ametani exceeds the values computed by Hara and Umoto's approach.
The single parameter required by the model of Maccioni-Araneo et al. is σ G , which defines the conductance accounting for the corona losses in its equivalent circuit representation. Instead, the capacitive current, denoted with i ∆C in Section 4.4.3, only depends on the geometrical features of the configuration under study and on the conductor voltage. Hence, the total charge corresponding to the applied voltage (which is derived by integrating the additional corona capacitive current i ∆C with respect to time and accounting for the instantaneous electrostatic charge) is displayed in Figure 10c (the independence of any input parameters results in a single q-v loop).

Discussion
The available methods for the computation of the corona onset voltage for DC and AC applications show good agreement with the measured values [90] and produce results with differences up to 4% [91], which are acceptable for the estimation of conductors' surface voltage gradients. Inception voltages computed using Peek's equation agree well, provided that an appropriate value of the irregularity factor m is included [92]. Despite its simplicity, Peek's equation, which is largely used for engineering purpose, strongly depends on the correct choice of m, and this choice is not trivial since m varies significantly [38] with the conductors geometry, surface state, and weather conditions. Hence, it is necessary to derive a suitable value from past experience and/or available measurements.
Moving to stranded conductors-being that the electrical field at the tip of each strand is about 14% higher than the field on the surface of a smooth cylindrical conductor of the same diameter [91]-and bundles, the situation is even less trivial. Usually, engineers refer to an equivalent conductor, which approximates the bundle [93]. However, the surface electric field of the equivalent is lower than that of the bundle, so further corrections are needed when calculating the onset voltage. Despite the large literature, difficulties to correctly quantify the corona onset voltage still arise in these cases.
Successively, we have discussed and compared the reviewed models in terms of their advantageous and disadvantageous characteristics. A practical overview on the considerations to be made when choosing a specific model for implementation is presented in Table 7. Only Correia de Barros' approach models the drift of the charged particles. Hence, it allows to consider the dynamic of the p.u.l. charge approaching zero for large simulation times and decaying values of the applied voltage, at the cost of a heavier computational burden. On the contrary, Malik and Cooray's models, being based on a macroscopic approximation of the physical phenomenon, are easier to implement. In particular, since the shape of the q-v loops computed through Malik's model is narrower and hardly dependent on the input parameters (see Figure 8b), the model may not be suitable to fit differently shaped measured q-v curves.
Only Inoue's model does not predict discontinuities in the p.u.l. capacitance value at corona inception; with reference to Suliciu's approach, a discontinuity may be introduced depending on the values attributed toṽ j (j = 1, 2). The introduction of a harsh discontinuity of the p.u.l. dynamic capacitance (i.e., an abrupt change in the propagation velocity of traveling waves) requires careful implementation in FDTD approaches in order to avoid numerical instability [51].
Finally, the implementation of circuit-based models requires a compromise to be made between the simplicity of plugging a limited number of shunt branches to the TL under study and the accuracy in modeling the corona effect as a non-lumped phenomenon. Table 7. Features of the reviewed models.

Advantages Disadvantages
Correia de Barros · Accurate physical modelling of the phenomenon; · Drift of the charges is simulated.
· Structured coding is necessary; · Requires longer computation times; · Several input parameters are required.
· Predicts narrow q-v curves; · Limited versatility when fitting experimental data.

Cooray
· Simple implementation; · Implementation of back corona in the descending branch of q-v loops.
· Predicts an instantaneous q-v relation (unlike other physics-based models).

Inoue
· Simple implementation; · Versatility when fitting experimental data.
· Strong dependence on the input parameters: questionable predicting capability; · No distinction between positive and negative corona.
Gary et al.
· Simple implementation · Distinction between positive and negative corona.
· Predicts the discontinuity of the dynamic capacitance at the corona inception.
Podporkin and Sivaev · Simple implementation. · Predicts a step discontinuity of the dynamic capacitance at the corona inception.
Suliciu · Allows to model q-v loops presenting non-simultaneous maximum values of charge and voltage.
· Several input parameters are required.
CIGRÉ · Simple implementation. · Predicts the discontinuity of the dynamic capacitance at the corona inception.

Hara and Umoto
· Trivial implementation in programs for electro-magnetic transient studies.
· Less accurate representation of corona as a distributedphenomenon.
· Excessive simplification of the non-linear phenomenon, through constant-valued additional ∆C; · Less accurate representation of corona as a distributed phenomenon.
Maccioni, Araneo et al. · Computed q-v curves have fixed shape: they cannot be adjusted to better match measured data; · Less accurate representation of corona as a distributed phenomenon.

Conclusions
This paper reviews and summarizes the main phenomena involved with corona discharge: microscopic and macroscopic dynamics, electric field distribution at the inception. Special attention is devoted to the influence of corona discharge on the p.u.l. parameters of overhead transmission lines, referring to the additional power losses and the enhanced capacitive behavior, which are addressed in the second part of the paper. This point is ad-dressed through a comprehensive and critical review of a wide selection of models available in the literature for the simulation of the relation between p.u.l. charge and voltage.
With respect to a chosen voltage source applied to a reference conductor, the corona models are implemented to compute the charge-voltage hysteretic loops and evaluate the models' sensitivity to variations in the required set of input parameters.
Most of the models display a remarkable dependence on the input parameters, which may be properly adjusted in order to reproduce charge-voltage measured loops. However, due to this tunability feature, the models are suitable for reproducing experimental data rather than for predictive studies of the electrical behavior of transmission lines.
Further updates in research are necessary to study approaches of general applicability with limited dependence on additional parameters, which may be exploited to assess the influence of corona discharge on propagation phenomena along transmission lines regardless the availability of pre-existing measured values of charge and voltage in the analyzed configurations.