Practical Model Applicable to Investigating the Coast Effect on the Geoelectric Field in Connection with Studies of Geomagnetically Induced Currents

The horizontal geoelectric field induced at the Earth's surface during temporal variations of the geomagnetic field drives Geomagnetically Induced Currents (GIC) in ground-based conductor networks, such as electric power transmission grids. Usually studies of the geoelectric field in connection with GIC research employ the assumption of a one-dimensional (1D) Earth conductivity structure. At ocean-land interfaces the conductivity contrast in the horizontal direction is large making the 1D assumption invalid in coastal areas. Geoelectromagnetic induction literature contains many publications dealing with complicated 2D or 3D Earth models. However, these studies mainly aim at the investigation of the Earth's structure and often do not provide an easy theoretical discussion and clear numbers to be used in practical power engineering investigations of GIC in networks located at coasts. Regarding GIC applications, James Gilbert introduces a two-dimensional (2D) Earth conductivity model for studying the " coast effect " near an 10 R. Pirjola ocean-land interface and derives an integral equation for the surface geoelectric field [4]. In this paper, we show that a model similar to that discussed by Gilbert can also be handled by forward calculations without an integral equation. Thus the treatment considered in this paper is simple enough to be easily understood physically and mathematically but complicated enough to provide realistic information about the coast effect on the geoelectric field for practical GIC purposes. Preliminary numerical computations indicate a satisfactory agreement with Gilbert's results about the enhancement of the horizontal geoelectric field on the land near a coastline.


Introduction
At the Earth's surface, Space Weather, which is driven by the activity of the Sun, manifests itself as Geomagnetically Induced Currents (GIC) in ground-based conductor networks, such as electric power transmission grids, oil and gas pipelines, telecommunication cables and railway circuits, e.g.[2].In general, GIC constitute a hazard to the particular system.In power networks, this is due to saturation of transformers, which can lead to blackouts and permanent damage of transformers, e.g.[7].In pipelines, space weather may cause problems associated with corrosion and its control, e.g.[5].So far, GIC effects on railway systems have not been much investigated, but clear indications of them exist, e.g.[3,13,14].The history of GIC dates back to early telegraph equipment in the 1800's, e.g.[2,9].Due to the use of optical fibre cables, today's telecommunication systems are probably rather insensitive to GIC.Power networks constitute the most important systems regarding GIC issues nowadays.
The horizontal geoelectric field, which is induced during a temporal variation of the geomagnetic field, as stated by Faraday's law of induction, is the actual cause of GIC in a network.Thus GIC research much concentrates on the determination of the geoelectric field.This field is primarily created by currents and charges in the ionosphere and magnetosphere and secondarily affected by currents and charges produced in the conducting Earth.A common assumption in model studies of the geoelectric field is that the conductivity of the Earth only varies with depth, i.e. the Earth structure is one-dimensional (1D).This assumption simplifies the modelling, makes numerical computations faster and is usually acceptable in connection with GIC studies (at least as the first approximation).At this point, it should be noted that the geoelectric field is integrated along paths determined by the transmission lines to calculate GIC in a power network [10].Thus small spatial details are unimportant in the geoelectric field and thus need not be known in GIC studies.
It is well known that the horizontal geoelectric field is amplified in coastal areas, which is due to the large lateral conductivity contrast at an ocean-land interface.This phenomenon is called the "coast effect".The amplification, which concerns the electric component perpendicular to the coastline, is physically due to charges accumulating at the ocean-land boundary, thus ensuring the continuity of the current flow.The enhancement of the geoelectric field means that GIC also tend to be large in conductor networks located at and near a coast.In spite of the qualitative understanding of the amplification of the geoelectric field in coastal areas, a quantitative description is more difficult.Because of the lateral change of the conductivity, a 1D Earth model is naturally insufficient to describe the coast effect.In other words, precise knowledge of the magnitude of the amplification and of the distance how far away from the coastline the amplification extends requires a better conductivity model.On the other hand, as mentioned above, spatial details of the geoelectric field are insignificant.Taking also into account that GIC computations usually have to be fast enough to be of any practical help, it is understandable that an Earth conductivity model appropriate to investigating the coast effect in connection with GIC studies should not be too complicated and detailed.
It is necessary to clearly emphasise that the published literature on geoelectromagnetic induction contains many papers that deal with sophisticated 2D and 3D Earth (and ocean) conductivity models accompanied by elaborate numerical computations.However, the application of those models is generally focussed on investigations of the ground structure including geological interpretations.Therefore, the models are obviously unnecessarily complicated and thus impractical to be utilised and considered, for example, by the power engineering community concerned about GIC in power networks located in coastal areas.Furthermore, it seems to be difficult to extract concrete numerical values from the literature that easily provide an estimate of the enhancement of the geoelectric field as a function of the distance from the coast.A principal objective of this paper is to remove this shortcoming by introducing a model to be applied in connection with GIC studies that is simple enough to be easily understood physically and mathematically but complicated enough to provide realistic information R. Pirjola about the coast effect on the geoelectric field.The work discussed in this paper thus especially has practical significance for the research of ground effects of space weather, and so it should not be regarded as a conventional geoelectromagnetic or magnetotelluric investigation of the Earth's conductivity structure.
In his impressive and valuable paper [4], James Gilbert presents a two-dimensional (2D) model suitable for studying the coast effect on the geoelectric field.Gilbert particularly stresses GIC applications of his work.In his model, the depth of the ocean is a linearly increasing function of the distance from the coast, and the total magnetic field is considered known at the surface of the land and at the ocean surface.By assuming that the (linear) spatial change of the ocean depth is sufficiently slow and that the fields do not vary too much in the horizontal direction, the plane wave and layered-Earth formulas presented in [1] can be used approximately.Gilbert investigates both the "H polarisation" and the "E polarisation" case, which refer to the magnetic field parallel and perpendicular to the coastline, respectively.The former case, where the electric field is perpendicular to the coastline, is more interesting and important regarding GIC studies since it exhibits the pronounced enhancement of the electric field.For the H polarisation, Gilbert ends up at an integral equation for the electric field at the ocean floor and at the land surface, which together form the xy plane of his Cartesian coordinate system with the y axis lying at the coastline (equations (10) and (12) in [4]).In the integral equation, the electric field is a function of x, and the kernel is a modified Bessel function K 0 thanks to the simplifying assumption that the land and the Earth below the ocean constitute a uniform half-space.
In this paper, we make use of Gilbert's study [4] but the geophysical model applied here has a small, though evidently insignificant, difference.The theory is also formulated in a more general way as concerns the ocean depth.Moreover, the general equations do not presume that the Earth is uniform but any layered structure is acceptable.We consider the H polarisation case and derive a forward expression for the perpendicular x component of the geoelectric field, thus making the mathematical treatment simpler by avoiding the integral equation.The electric field can be calculated by utilising an inverse Fourier transform or a convolution integral.Numerical computations by using the formulas are also discussed including comparisons with Gilbert's results.

Practical model applicable to investigating the coast effect
13 2. Theory

General model
Let the surfaces of the land and the ocean define the xy plane of a Cartesian coordinate system with the z axis pointing downwards.(Ignoring the sphericity of the Earth naturally means that the model is applicable to local and regional studies.)We assume that the ocean depth denoted by d is a function of x, but independent of y.In a model similar to that in [4], the function d = d(x) would have the expression where η gives the slope of the linearly increasing depth and Θ(x) is the Heaviside step function being zero and one for negative and positive values of x, respectively.Thus the coastline coincides with the y axis, and the land and the ocean lie at x < 0 and x > 0, respectively.It should be noted that the present assumptions combined with equation (1) do not exactly correspond to Gilbert's model [4] because his xy plane coincides with the plane consisting of the land surface at x < 0 and of the ocean floor at x > 0. So, in his model, the land and ocean surfaces make an angle and the coordinate system is tilted (see Fig. 1 in [4]).Naturally in practice, a linear increase of the ocean depth is realistic only until a certain distance x, after which d is more or less constant.
In the present theoretical discussion, we do not restrict the treatment to a linear variation of the depth as given by equation ( 1).So we keep the function d(x) more general, but the change of d with x is assumed to be slow enough to enable the use of the layered-ocean/Earth approximation.The (total) horizontal magnetic (H) field is considered known at the xy plane.Let us assume that it is parallel to the y axis (H polarisation) and varies sufficiently slowly with x and y, so that the plane wave approximation is valid.We denote this H y component by H 0 (= H 0 (x,y)).In this paper, the dependence of all fields on the time t is assumed to be sinusoidal (e iωt ) with the angular frequency ω, whereas the term e st is used in [4] referring to the Laplace transform.
By definition, the surface impedance is a quantity that provides relations between electric and magnetic fields at a surface, e.g.[8].The Earth underlying the ocean, i.e. located at z > d(x), is assumed to be layered and characterised by the plane wave surface R. Pirjola impedance Z d .In principle, Z d is local depending on x due to variation of d, i.e. the z coordinate of the ocean floor varies with x.However, we assume the variation of d to be small enough compared to other length scales of the ocean and of the Earth structures, so that Z d can be regarded as independent of x.
Based on the above assumptions and on equation ( 13) in [1], we obtain a formula for the magnetic field at the ocean floor at z = d(x).It is parallel to the y axis and denoted by H d (= H d (x,y)): where the propagation constant k and characteristic impedance Z of the ocean water in the quasi-static approximation are defined by and The conductivity of the ocean water is σ, and µ 0 denotes the permeability assumed to equal the vacuum value (= 4π⋅10 -7 VsA -1 m -1 ) everywhere.(Note that the propagation constant is often defined as −iωµ 0 σ , i.e. differently from equation (3).)It must still be emphasised that equation ( 2) is exact only in the strict plane wave and layeredocean/Earth case, so now when d depends on x and H 0 on x and y, equation ( 2) is only approximate but it is assumed to be acceptable.
The following definitions of the Fourier transform and the inverse Fourier transform of any function f from the space domain (x) to the wavenumber domain (b) are used throughout this paper: Using equation ( 6) we obtain E d in the space domain: The possible dependence of E d on y is implicit in formula (7).In order to calculate the electric field at the surface (z = 0), we make use of equation (15) in [1].This electric field denoted by E 0 (= E 0 (x,y)) is parallel to the x axis and has the expression Equations ( 2), (5) (with f(x) = H d (x,y)), ( 7) and ( 8) enable the calculation of E 0 (x,y) from d(x) and H 0 (x,y), provided that the Earth's conductivity structure at z > d(x) and the ocean water conductivity are known.This procedure involves a Fourier transform (equation ( 5)) and an inverse Fourier transform (equation ( 7)) but no integral equation needs to be solved.
The convolution theorem can be used to replace equations (5) (with f(x) = H d (x,y)) and (7) with where Z u (x) is the inverse Fourier transform of Z u (b) (equation ( 6)) and, as above, the possible dependence of E d and H d on y is implicit.Consequently, instead of using equations ( 2), ( 5), ( 7) and ( 8), it may sometimes be more practical to calculate Z u (x) and then apply equations ( 2), ( 9) and ( 8) (in this order) to the determination of E 0 (x,y) for different [H 0 (x,y), d(x)] combinations.
Although d = d(x) giving the ocean depth is an arbitrary function of x in the theoretical treatment presented, in practical situations d should be equal to zero for a large range of x values, which give the locations of land areas.For these x values, E 0 and E d are naturally equal as also directly seen from equation ( 8).An appropriate form of d(x) to be used in studies of the coast effect would seem to be that given by equation ( 1), but we have to remember that ηx may also be replaced with another function of x.

Special case
We now ignore the possible y dependence.If, similarly to [4], the Earth underlying the ocean (as well as at land areas) is uniform with a conductivity σ 1 , Z u (b) has the following (quasi-static) expression corresponding to the two-dimensional H polarisation case: where the propagation constant k 1 is defined by formula (3) replacing σ with σ 1 .The validity of equation ( 10) can be seen, for example, from equation ( 6) in [4], but equation (10) can also be easily derived based directly on Maxwell's equations in the quasi-static approximation with the assumptions that the fields are independent of y and that the magnetic field is parallel to the y axis.The computation of Z u (x) by using equation ( 6) with Z u (b) from equation ( 10) is clearly impossible.(Note that, for the corresponding E polarisation surface impedance, the computation of the inverse Fourier transform would be possible.)The difficulty in the present H polarisation situation arises from the fact that Z u (b) in equation (10) does not go to zero but actually approaches infinity, when |b| increases, so the integral in equation ( 6) does not converge.However, referring to equations ( 15) and ( 16) in [12], we see that the inverse Fourier transform of 1 Z u (b) is proportional to the Hankel function of the second kind and of the zero order H 0 (2) (−ik 1 x).By noting further that H 0 (2) (−ik 1 x) is proportional to the modified Bessel function K 0 (k 1 x), the present formulation is similar to that introduced in [4], as expected.

Theoretical comments on the approximations
The Earth underlying the ocean is handled by using the plane wave surface impedance Z d in equations ( 2) and ( 8) whereas the wavenumber-dependent surface impedance Z u (b) is considered in equation ( 7) and the inverse Fourier transform Z u (x) of Z u (b) in equation ( 9).This is, of course, in principle inconsistent.However, we should notice that the application of Z d in equations ( 2) and ( 8) concerns the transfer through the highly-conducting ocean water whereas Z u is used for relating E d to H d at the ocean floor, and generally large conductivities are more favourable for the plane wave assumption, e.g.[12] and references therein.On the other hand, it must be remembered that the derivation of equations ( 13) and (15) in [1], which are utilised in equations ( 2) and ( 8) above, applies the formula Z d = E d /H d at the ocean floor.All this means that approximations included in the technique described above for calculating E 0 (x,y) in terms of H 0 (x,y) are theoretically not consistent and fully justified, though probably acceptable.A similar approximation also seems to be made in [4] where equation (2) clearly corresponds to the plane wave assumption but equation ( 6) includes the wavenumber dependence.
It was mentioned in Section 1 that the coast effect is physically a result from charge accumulation at the ocean-land boundary, which implies the continuity of the current flow by enhancing the electric field in the land and decreasing the electric field in the ocean.In Section 2.1, we first say that the electric field E d is parallel to the ocean floor and later regard E d as parallel to the x axis.Both are approximations.When E d is parallel to the floor the electric current in the ocean actually does not carry charges to the oceanland boundary surface (which is in fact the case in the treatment in [4] as well).An electric field parallel to the x axis leads to charge accumulation at the boundary but such an electric field obviously does not correspond to the real situation either.This means that, while the formulas derived in Section 2.1 are applicable to calculations of the coast R. Pirjola effect on the geoelectric field, physical details of the processes involved cannot be concluded from them.

Numerical results
The formulas derived in Section 2 for calculating E 0 (x,y) involve Fourier integrals.This would indicate the use of the Fast Fourier Transform (FFT) in numerical computations.An alternative is the Fast Hankel Transform (FHT) [6] successfully utilised, for example, in another geophysical model study that also included Fourier integrals [11].However, in this paper, whose objective is to go through and introduce the theory for a handy method for investigating the coast effect in association with GIC research that also corresponds to the technique presented in [4] but avoids the integral equation, we do not go deep into considering and evaluating numerical computation tools.Thus, we now just compute the Fourier integrals by approximating them simply with sums and performing appropriate smoothing of the results to remove unnatural oscillations due to numerical inaccuracies.
We assume that the known surface magnetic field H 0 is uniform all over the xy plane.Its value is set equal to 0.0796 Am -1 , which corresponds to a magnetic B field of 100 nT (B = µ 0 H).The depth of the ocean follows equation (1) with a slope η of 1/30 for ≤ 150 km, and for x > 150 km the depth is constant (= 5 km).The ocean water conductivity is assumed to be σ = 4 Ω -1 m -1 .The Earth is uniform corresponding to the special case discussed in Section 2.2.The perpendicular horizontal electric field component 0 now has to be solved from equations (2), ( 5), (7) and ( 8) because the function Z u (x) needed in the convolution (9) cannot be computed (see Section 2.2).Note that E 0 is a function of x but does not depend on y because H 0 is uniform.
We assume different values for the period of oscillation T = 2π/ω and for the Earth's conductivity σ 1 .The solid curve in Fig. 1 depicts the absolute value of E 0 , which is parallel to the x axis, at the land surface when T and σ 1 equal 60 s and 10 -3 Ω -1 m -1 , respectively.For comparison, the small circles give the corresponding results from Fig. 2 in [4], which is a plot in scaled units on both axes.For Fig. 1, data from [4] are obtained simply by making measurements with a ruler from Fig. 2 in [4] and by performing the appropriate scalings to absolute units.As pointed out in Section 2.1, the time factor e st referring to the Laplace transform is used in R. Pirjola practical applications to GIC research, the field at the land is of importance.Therefore, and because computations of the geoelectric on the ocean side of the interface seem to suffer from larger numerical problems, we restrict the study to the land side in this paper.

Concluding remarks
The horizontal geoelectric field occurring at the Earth's surface during a space weather storm is the key quantity regarding Geomagnetically Induced Currents (GIC) in ground-based conductor networks.In connection with GIC studies, the geoelectric field is commonly modelled by assuming that the Earth has a one-dimensional (1D) structure.Due to the fact that GIC are controlled by spatial line integrals of the geoelectric field, small lateral inhomogeneities of the Earth's conductivity do not play a significant role making the 1D assumption usually acceptable.However, in coastal areas the vicinity of the ocean-land interface implies that a 1D model is not correct.
Regarding practical power engineering work on GIC issues, the many versatile 2D and 3D Earth/ocean models, which are published in geoelectromagnetic literature and which generally aim at investigating the ground structure and geology, are not the most appropriate.With an emphasis on GIC applications an approximate model and technique are introduced in [4] for a quantitative investigation of the coast effect on the geoelectric field, which is obtained from an integral equation.In this paper, we make use of [4] but formulate and solve the geoelectric field in terms of a forward expression, and so the analysis becomes simpler as the integral equation is avoided.
Equations ( 2), ( 5), ( 7) and (8) or equations ( 2), ( 9) and ( 8) of this paper enable the calculation of the horizontal electric field E 0 (x,y) perpendicular to the coastline at the ocean and land surfaces (= the xy plane) when the magnetic field H 0 (x,y) parallel to the coastline is known at this plane.The depth of the ocean d, which depends on the x coordinate, the ocean water conductivity and the layered Earth structure are also assumed to be known.A useful form of d = d(x) is expressed by equation (1).It states that the ocean-land interface lies at x = 0, and a linearly deepening ocean occupies the region x > 0. But the linear dependence ηx may also be replaced by any other function of x.Utilising equations (2), ( 9) and ( 8) in practice, it is necessary to calculate the kernel function Z u (x) only once, as it remains the same for different H d functions in the integrals in equation (9).However, if the Earth is uniform (as in the numerical examples in this paper), Z u (x) is not possible to be computed.
Results from numerical calculations performed in this paper are not the same as those in [4] but the agreement can be considered satisfactory, which would seem to justify the use of the equations of this paper for quantitative estimations of the effect of an ocean-land interface on the horizontal geoelectric field in coastal areas.However, numerical computations would definitely require additional studies later to make them more accurate and efficient.Future efforts should also include comparisons between geoelectric field values obtained by applying the model presented in this paper and data measured near coastline.

Table 1:
Real and imaginary parts of the electric field values at the land point corresponding to the left-hand ends of the solid curves in Figures 1-4.This point, which lies at about x = -300 km, represents the most distant site from the ocean-land interface (located at x = 0) included in the calculations of this paper.The Figure and Plane Wave columns show the values corresponding to Figures 1-4 and the values that are obtained in the case of no ocean, respectively.

R. Pirjola
Fig. 1: The solid curve depicts the absolute value of the horizontal geoelectric field at the land surface calculated by using equations ( 2), ( 5), ( 7) and ( 8).The geoelectric field, which is shown as a function of the distance from the ocean-land interface, is perpendicular to the interface located at x = 0.The period T and the Earth's uniform conductivity σ 1 equal 60 s and 10 -3 Ω -1 m -1 , respectively.The values of the other parameters are given in the text in Section 3. The small circles give the corresponding results from Fig. 2 in [4].

R. Pirjola
x dependencies of d and H 0 make H d given by equation (2) a function of x.Substituting H d (x) for f(x) in equation (5) results in H d (b), where the possible y dependence is implicit.Multiplying H d (b) by the b dependent surface impedance Z u (b) of the Earth occupying the region z > d(x) gives the Fourier transform of the perpendicular electric field component, denoted by E d , at the ocean floor and parallel to the floor.Based on the above assumption that the variation of d with x is small, Z u (b) corresponds to a layered Earth with a flat surface, and we regard E d as parallel to the x axis.Note that the plane wave value (b = 0) of Z u (b) equals Z d introduced above.