Transmission Line Modelling of Geomagnetic Induction in the Ocean/Earth Conductivity Structure

David H. Boteler¹, Shibaji Chakraborty², Xueling Shi^2,3, Michael D. Hartinger⁴, Xuan Wang⁵

¹Geomagnetic Laboratory, Natural Resources Canada, Ottawa, Canada.
²Center for Space Science and Engineering Research, Virginia Tech, Blacksburg, USA.
³High Altitude Observatory, National Center for Atmospheric Research, Boulder, USA.
⁴Space Science Institute, Boulder, USA.
⁵Department of Electrical Engineering, Tsinghua University, Beijing, China.

DOI: 10.4236/ijg.2023.148041                

Abstract

During geomagnetic disturbances, electric fields induced in the Earth and in power systems, pipelines and submarine cables can interfere with the operation of these systems. Calculations for submarine cables are complicated by the need to consider not just the induction directly into the cable but also the earth potentials produced at the coast at each end of the cable. To determine the coast potentials, we present a new model of the ocean and earth conductivity structure that spans the whole length of a cable from one coast to another. Calculations are based on the generalised thin sheet approach introduced by Ranganayaki and Madden but converted to a transmission line model that can be solved using standard circuit theory techniques. It is shown how the transmission line model can be used to calculate the earth potential profile from one side of an ocean or sea to the other. Example calculations are presented for a shallow sea, a shallow ocean, and a deep ocean that are simplified approximations to the North Sea, Tasman Sea and Pacific Ocean and show that the peak potentials occur at the coast. An examination is also made of how the width of a shallow sea and the width of the continental shelf affect these coast potentials. The modelling technique and example results provide a guide for more detailed modelling of geomagnetic induction along the routes of specific submarine cables.

Keywords

Geomagnetic Induction, Generalised Thin Sheet, Transmission Line Modelling, Coast Potentials, Submarine Cables

1. Introduction

During geomagnetic disturbances, the magnetic field variations produce electric fields at the surface of the Earth that can cause problems for critical infrastructure on land such as power systems, pipelines, and railways [1] - [6] . Electric fields are also produced in the sea and can affect submarine cables on the seafloor [7] . Modelling the voltages experienced by submarine cables is complicated as it needs to take account of both geomagnetic induction into the cable itself as well as the earth potentials at the coast for the ground connections of the cable power feed equipment [8] . This paper is concerned with the latter problem of how to model geomagnetic induction in the combined ocean and earth conductivity structure to determine the earth potentials at the coasts on either side of a sea or ocean.

Electric fields induced during geomagnetic disturbances are dependent on the amplitude and frequency content of the magnetic field variations and on the conductivity structure of the Earth. Simple one-dimensional (1-D) models of the Earth conductivity have often been used to calculate electric fields [9] [10] . However, such calculations fail to take account of the effect of boundaries between regions of different conductivities. The largest change in conductivity is that between the resistive land and the conducting seawater that occurs at the coast. This change in conductivity results in a build-up of electrical charge at the coast that modifies the electric field in such a way as to produce current continuity across the coast. This “geoelectric coast effect” increases the electric field on the land side of the coast and is a factor that needs to be considered when assessing the geomagnetic hazard to power systems adjacent to the coast and submarine cables [11] [12] [13] .

Electric fields produced in the sea depend on the amplitude and frequency content of the magnetic field variations and on the conductivity structure of the Earth, just as for electric fields on land, but calculation of electric fields on the seafloor also needs to take account of the attenuation of the fields by the conducting seawater. Formulas for a 1-D Earth model and seawater layer have been developed by [14] . The conductivity boundary at the coast also influences the electric field in the sea and is a factor that needs to be considered when examining geomagnetic effects on submarine cables. Ranganayaki and Madden [15] pointed out the importance of the crustal resistance in controlling the leakage of electric currents from the sea through to the conducting mantle and introduced the generalised thin sheet model to represent this. A recent paper [16] introduces a simple method for generalised thin sheet modelling based on transmission line theory. This provides a new approach for modelling geomagnetic induction in the Ocean/Earth system comprising both the seawater and the earth conductivity structure beneath the seafloor and the land at either edge of the ocean.

In this paper we examine geomagnetic induction in the ocean taking account of the land on either side and leakage through the resistive crust to the mantle. This Ocean/Earth system is modelled as a multi-section transmission line with the sections represented by equivalent-pi circuits which are combined to form a nodal admittance network. Analysis of this network leads to a matrix equation that can be solved to give the Earth potentials across the Ocean/Earth System. Example results are presented for simplified models of a deep ocean, a shallow ocean and a shallow sea that roughly correspond to the Pacific Ocean, Tasman Sea, and North Sea. A closer examination is also made of the potentials at the coast and edge of the continental shelf and how they are influenced by the width of the continental shelf.

The purpose of these model calculations is to provide a preliminary understanding of how oceans and seas with varying widths and water depths are expected to respond to geomagnetic induction. These calculations serve as a guide, offering insights into the general behavior of the Ocean/Earth system. However, it is important to note that these results are based on simplified geometry and do not account for the intricate details of depth profiles and conductivity structures specific to certain locations. The intention is to use these model calculations as a starting point for further investigations that consider the detailed depth profiles and conductivity structures of specific areas. By incorporating these specific characteristics into more detailed calculations, a more accurate assessment of the geomagnetic induction response in those locations can be obtained.

The ultimate goal is to apply this refined model to conduct comprehensive analyses for specific locations, enabling a deeper understanding of the Earth’s potential distribution in those areas. This information would be invaluable for various applications, such as assessing the potential impact of geomagnetic induction on submarine cables, optimizing the design of electrical systems in coastal regions, and evaluating the overall electromagnetic response of the marine environment. Therefore, these initial model calculations serve as a valuable stepping stone towards more detailed investigations, allowing researchers to focus on specific locations and refine their understanding of the geomagnetic induction processes occurring in oceans and seas.

2. Ocean/Earth Conductivity Structure

The basic structure of the Ocean/Earth conductivity is shown in Figure 1. The surface layer is comprised of the deep ocean with shallow seas on either side and all these are underlain by a layer of sediments. On the land on either side of the ocean there will be a surface sedimentary layer that may be very thin or several kilometres thick. Below these are the Earth’s crusts, thinner under the oceans than under the land, and the resistive mantle lithosphere. Deeper in the Earth, in the upper mantle, transition zone and lower mantle the temperature and pressure increase cause partial melting which increases the conductivity. Resistivity values for these layers from the LITHO1.0 model [8] [17] and example seawater and layer thicknesses for the different sections are shown in Table 1.

Figure 1. Basic structure of the Ocean/Earth conductivity for a deep ocean with shallow seas and land on either side.

Table 1. Example parameters for test Ocean-Earth conductivity structure.

3. Transmission Line Modeling

For geomagnetic induction studies, the Ocean/Earth conductivity structure consisting of surface conducting layers (seawater and sediments) over more resistive layers (crust and mantle lithosphere) can be represented as a thin double layer (with the conductive layer on top and the resistive one underneath). Ranganayaki and Madden [15] developed a generalised thin sheet analysis for this structure that is appropriate for frequencies that penetrate through both the conductive and resistive layers. Wang et al. [16] have shown that this generalised thin sheet analysis can be accomplished by representing the double layer as a transmission line with series impedance, Z, given by the resistivity and thickness of the conductive layers and parallel impedance, Y given by the resistance through the resistive layers as shown in Figure 2.

For the surface conductive layers, we need to consider the path for horizontal currents. Thus, the conducting paths through the sea (Layer 1) and the sediments (Layer 2) are in parallel and the conductance, C, of the combined layers is

Figure 2. (a) Layered ocean-earth conductivity model and (b) its equivalent transmission line model. Modified from Wang et al. [16] .

simply the sum of the conductance of the individual layers. Thus

$C={\sigma }_1{d}_{1}w+{\sigma }_2{d}_{2}w$ (1)

where ${\sigma }_1$ and $d_1$ are the conductivity and depth of the seawater layer and ${\sigma }_2$ and $d_2$ are the conductivity and thickness of the sedimentary layer, and w is the width of the layers. These are used to give the series impedance Z for the transmission line model:

$Z=\frac{1}{C}$ (2)

For the resistive crust (Layer 3) and mantle lithosphere (Layer 4) we need to consider the path for vertical currents flowing between the upper conductive layers and the underlying upper mantle layer. Here the resistance of the two layers are in series and the combined resistance is simply the sum of the layer resistances. Thus

$R=\frac{{\rho }_3{d}_3}{w}+\frac{{\rho }_4{d}_4}{w}$ (3)

where ${\rho }_3$ and $d_3$ are the resistivity and thickness of the crust layer and ${\rho }_4$ and $d_4$ are the resistivity and thickness of the mantle lithosphere layer. As for the upper layers, w is the width of the layers. These are used to give the parallel admittance Y of the transmission line model:

$Y=\frac{1}{R}$ (4)

These values are then used to determine the characteristic impedance Z₀ and propagation constant, γ, of the transmission line model given by:

$\gamma =\sqrt{ZY}$ $Z_0=\sqrt{\frac{Z}{Y}}$ (5)

To use the transmission line model to examine the effects of electromagnetic induction we use distributed-source transmission line (DSTL) theory in which the induced electric field is represented by voltages sources distributed along the transmission line [18] [19] . Considering a short element of the transmission line, as shown in Figure 2(b) we can write the transmission line equations:

$\frac{\text{d}V}{\text{d}x}+ZI=E$ (6)

and

$\frac{\text{d}I}{\text{d}x}+YV=0$ (7)

involving the voltage, V, and current, I, produced by an induced electric field, E, in the transmission line. Differentiation and substitution of (6) and (7) lead to the equations

$\frac{{\text{d}}^{2}V}{\text{d}{x}^2}-{\gamma }^{2}V=\frac{\text{d}E}{\text{d}x}$ (8)

and

$\frac{{\text{d}}^{2}I}{\text{d}{x}^2}-{\gamma }^{2}I=-YE$ (9)

If the electric field is assumed to be uniform within a section of transmission line, $\frac{\text{d}E}{\text{d}x}=0$ , and Equation (8) becomes:

$\frac{{\text{d}}^{2}V}{\text{d}{x}^2}-{\gamma }^{2}V=0$ (10)

Equations (9) and (10) have solutions of the form:

$V=A{\text{e}}^{\gamma x}+B{\text{e}}^{-\gamma x}$ (11)

and

$I=\frac{E}{\gamma Z_0}-\frac{A}{Z_0}{\text{e}}^{\gamma x}+\frac{B}{Z_0}{\text{e}}^{-\gamma x}$ (12)

where A and B can be found from conditions at the ends of the transmission line.

At the start of any section, denoted by i^th node, $x=0$ , Equations (11) and (12) give

$V_i=A+B$ (13)

and

$I_i=\frac{E}{\gamma Z_0}-\frac{A}{Z_0}+\frac{B}{Z_0}$ (14)

while at the end of that section with length L, denoted by k^th node, $x=L$ , Equations (11) and (12) give

$V_k=A{\text{e}}^{\gamma L}+B{\text{e}}^{-\gamma L}$ (15)

and

$I_k=-\frac{A}{Z_0}{\text{e}}^{\gamma L}+\frac{B}{Z_0}{\text{e}}^{-\gamma L}+\frac{E}{\gamma Z_0}$ (16)

Equations (13) and (15) can be combined to solve for A and B

$A=\frac{V_k-V_i{\text{e}}^{-\gamma L}}{{\text{e}}^{\gamma L}-{\text{e}}^{-\gamma L}}$ (17)

and

$B=\frac{V_i{\text{e}}^{\gamma L}-V_k}{{\text{e}}^{\gamma L}-{\text{e}}^{-\gamma L}}$ (18)

Substituting for A and B in Equation (11) for voltage gives

$V=\left(\frac{V_k{\text{e}}^{\gamma L}-V_i}{{\text{e}}^{\gamma L}-{\text{e}}^{-\gamma L}}\right){\text{e}}^{-\gamma \left(L-x\right)}+\left(\frac{V_i{\text{e}}^{\gamma L}-V_k}{{\text{e}}^{\gamma L}-{\text{e}}^{-\gamma L}}\right){\text{e}}^{-\gamma x}$ (19)

Similarly substituting for A and B in Equation (12) for current gives

$I=-\frac{1}{Z_0}\left(\frac{V_k{\text{e}}^{\gamma L}-V_i}{{\text{e}}^{\gamma L}-{\text{e}}^{-\gamma L}}\right){\text{e}}^{-\gamma \left(L-x\right)}+\frac{1}{Z_0}\left(\frac{V_i{\text{e}}^{\gamma L}-V_k}{{\text{e}}^{\gamma L}-{\text{e}}^{-\gamma L}}\right){\text{e}}^{-\gamma x}+\frac{E}{\gamma Z_0}$ (20)

Equation (19) shows that the voltage on the transmission line is comprised of two components that fall off exponentially with distance from either end of the transmission line. The exponential fall-off has a characteristic length $\frac{1}{\gamma }$ referred to as the “adjustment distance”. For short and medium length transmission lines these two components overlap and the effects from both ends have to be taken into account to calculate the voltages. However, if the transmission line is sufficiently long then these contributions do not overlap. This “non-overlapping” criterion requires that the transmission line length, L, is greater than four adjustment distances, i.e. $L>4\frac{1}{\gamma }$ which then means ${\text{e}}^{\gamma L}\gg 1\gg {\text{e}}^{-\gamma L}$ and the transmission line is referred to as “electrically-long”. For such an “electrically-long” transmission line Equations (19) and (20) reduce to

$V=V_i{\text{e}}^{-\gamma x}+V_k{\text{e}}^{-\gamma \left(L-x\right)}$ (21)

and

$I=\frac{1}{Z_0}\left(V_i{\text{e}}^{-\gamma x}-V_k{\text{e}}^{-\gamma \left(L-x\right)}+\frac{E}{\gamma }\right)$ (22)

3.1. Thevenin and Norton Equivalent Circuits

The appearance of this electrically-long transmission line to another circuit can be represented by a Thevenin equivalent circuit comprising a voltage source, V_Th, in series with an impedance, Z_Th, as shown in Figure 3(a).

These can be calculated from the open-circuit voltage and short-circuit current [18] . Consider the view looking into end (node) “i” of the transmission line;

Figure 3. (a) Thevenin equivalent circuit, (b) Norton equivalent circuit.

in this case the “electrically-long” condition means that the contribution from the voltage at the other end (node) “k”, V_k is negligible and can be ignored. When the end (node) “i” of the transmission line is open circuit the current at $x=0$ is zero and equation (22) gives

$V_{oc}=-\frac{E}{\gamma }$ (23)

When the end (node) “i” of the transmission line is short-circuited the voltage at that end, $V_i=0$ and equation (22) gives

$I_{sc}=\frac{1}{Z_0}\frac{E}{\gamma }$ (24)

Thus, the Thevenin components are

$V_{Th}=V_{oc}=-\frac{E}{\gamma }$ (25)

$Z_{Th}=\frac{V_{oc}}{I_{sc}}=Z_0$ (26)

A similar derivation can be done for the other end (node) “k” of an electrically-long transmission line and gives the same Thevenin impedance but a Thevenin voltage $V_{Th}=\frac{E}{\gamma }$ . The Thevenin voltage represents an electric field directed towards the end of the transmission line (as at end “k”) so gives a voltage with a minus sign where the electric field is directed away from the end (as at end “i”).

The Thevenin equivalent circuit can be converted to a Norton equivalent circuit (Figure 3(b)), comprising a current source, J_N, in parallel with an admittance, Y_N, given by

$J_N=\frac{V_{Th}}{Z_{Th}}=\frac{E}{\gamma Z_0}=\frac{E}{Z}$ (27)

$Y_N=\frac{1}{Z_{Th}}=\frac{1}{Z_0}$ (28)

3.2. Equivalent-Pi Circuit

To combine multiple transmission line sections, it is useful to convert each section into equivalent-pi circuits [19] as shown in Figure 4.

The equivalent-pi circuit can be represented with a voltage source, E', and series impedance, Z', admittances, $\frac{Y^{\prime }}{2}$ , to ground at each end as shown in Figure 4(a) [19] with circuit components given by

$Z^{\prime }=Z_0\sinh \gamma L$ (29)

$\frac{Y^{\prime }}{2}=\left(\cosh \gamma L-1\right)\frac{1}{Z_0\sinh \gamma L}$ (30)

$E^{\prime }=\frac{E}{\gamma }\sinh \gamma L$ (31)

Alternatively, the equivalent-pi circuit can be represented with an equivalent

Figure 4. Equivalent-pi circuit for a transmission line section. (a) with a voltage source, E', and series impedance, Z'; (b) with a current source, I_E and parallel admittance, Y_E.

current source, I_E, in parallel with a series admittance, Y_E, as shown in Figure 4(b). In this case the equivalent-pi components are given by:

$Y_E=\frac{1}{Z_0\sinh \gamma L}$ (32)

$\frac{Y^{\prime }}{2}=\left(\cosh \gamma L-1\right)\frac{1}{Z_0\sinh \gamma L}$ (33)

$I_E=\frac{E}{Z}$ (34)

3.3. Transmission Line with “Active” Terminations

For electromagnetic induction into a transmission line with multiple sections we can consider a single section and represent the sections on either side by their Thevenin equivalent circuits as shown in Figure 5.

When the transmission line length is considerably shorter than the adjustment distance it is referred to as “electrically-short” and the equivalent-pi components reduce to [20] :

$Z^{\prime }=ZL$ (35)

$\frac{Y^{\prime }}{2}=0$ (36)

$E^{\prime }=EL$ (37)

In this case, electromagnetic induction in the transmission line in Figure 5 can be represented by the circuit shown in Figure 6.

The current along the transmission line is then given by:

Figure 5. Transmission Line with distributed series impedance, Z, parallel admittance, Y, and voltage sources representing the electric field, E, with “active” terminations at each end represented by Thevenin equivalent circuits.

$I=\frac{V_1+EL+V_2}{Z_1+ZL+Z_2}$ (38)

The potential at the left end of the transmission line is then:

$U_1=V_1-I{Z}_1$ (39)

And the potential at the right end of the transmission line is:

$U_2=I{Z}_2-V_2$ (40)

For the special case of a boundary between two regions that both extend far enough that they can be represented by Thevenin equivalent circuits for “electrically-long” transmission lines as shown in Figure 7.

The current across the boundary is given by

$I=\frac{V_{Th1}+V_{Th2}}{Z_{Th1}+Z_{Th2}}$ (41)

And the potential at the boundary is:

$U_{boundary}=V_{Th1}-I{Z}_{Th1}=I{Z}_{Th2}-V_{Th2}$ (42)

4. Network Modelling

Geomagnetic induction affects all the sections of the Ocean/Earth conductivity structure shown in Figure 1 and Table 1. Thus, in the transmission line model of this structure the voltages produced in one section are not just the result of

Figure 6. Equivalent circuit for an “electrically-short” transmission line with “active” terminations.

Figure 7. Boundary between Region 1 and Region 2 with each region represented by its Thevenin equivalent circuit.

electric fields in the neighbouring sections. Thus, we need to consider all the induced electric fields in that section but are also affected by the induced transmission line sections together and their influence on each other. This is done by using the equivalent-pi circuits for the ocean and shallow sea sections and Norton equivalent circuits for the land at either end to create a network model of the Ocean/Earth conductivity structure.

The circuit components for the basic Ocean/Earth structure are shown in Figure 8(a). This considers an east-west cross-section across an ocean which comprises 3 sections: middle one for the deep ocean and one at each side for the continental shelf. Each section is represented by its equivalent-pi circuit. The land at each end is represented by its Norton equivalent circuit with components J_W and Y_W at the western end, and J_E and Y_E at the eastern end. The admittances to ground from adjacent equivalent-pi circuits and from the Norton equivalent circuits are combined to give the nodal admittance network as shown in Figure 8(b).

To calculate the voltages produced in this transmission line model, the first step is to combine admittances from neighbouring sections in Figure 8(a) to give the admittance to ground from each node in Figure 8(b). Nodes 1 and 4 involve the admittance from the Norton equivalent circuit and the admittance from the equivalent-pi sections for the shallow sea:

Figure 8. Circuit models for the Ocean/Earth conductivity structure. (a) Equivalent-pi circuits for the transmission line models of each seawater section and Norton equivalent circuits for the land on each side. (b) Nodal admittance network.

$y_1=\frac{{y^{\prime }}_A}{2}+y_W$ $y_4=\frac{{y^{\prime }}_C}{2}+y_E$ (43)

Nodes 2 and 3 involve the admittances from the equivalent-pi circuits for the deep ocean and shallow seas:

$y_2=\frac{{y^{\prime }}_A}{2}+\frac{{y^{\prime }}_B}{2}$ $y_3=\frac{{y^{\prime }}_B}{2}+\frac{{y^{\prime }}_C}{2}$ (44)

The series admittance of the equivalent-pi sections becomes the series admittance of the corresponding section in the nodal network:

$y_{12}=y_A$ $y_{23}=y_B$ $y_{34}=y_C$ (45)

Similarly, the current sources from the equivalent-pi sections become the current sources in the nodal network

$j_{12}=j_A$ $j_{23}=j_B$ $j_{34}=j_C$ (46)

and the current sources for the Norton equivalent circuits at the ends are also connected to nodes 1 and 4.

Applying Kirchoff's current law that the algebraic sum of the currents entering any node is zero, i.e. the sum of the currents entering the nodes from the lines equals the current flowing to ground, we can write equations for each node:

$i_W-i_{12}=i_1$ (47a)

$i_{12}-i_{23}=i_2$ (47b)

$i_{23}-i_{34}=i_3$ (47c)

$i_{34}-i_E=i_4$ (47d)

The current in each line is determined by the potential difference between the nodes at the ends of the line, the admittance of the line and the equivalent current source for the line

$i_{12}=\left(v_1-v_2\right)y_{12}+j_{12}$ (48a)

$i_{23}=\left(v_2-v_3\right)y_{23}+j_{23}$ (48b)

$i_{34}=\left(v_3-v_4\right)y_{34}+j_{34}$ (48c)

While the currents from the Norton equivalent circuits are given by the current source values

$i_w=j_W$ (49a)

$i_E=j_E$ (49b)

The currents to ground at each node i_i are given by the node voltages times the admittances to ground, $i_i=v_i{y}_i$ . Making these substitutions in Equation (47) above gives for each node:

$j_W-\left(v_1-v_2\right)y_{12}-j_{12}=v_1{y}_1$ (50a)

$\left(v_1-v_2\right)y_{12}+j_{12}-\left(v_2-v_3\right)y_{23}-j_{23}=v_2{y}_2$ (50b)

$\left(v_2-v_3\right)y_{23}+j_{23}-\left(v_3-v_4\right)y_{34}-j_{34}=v_3{y}_3$ (50c)

$\left(v_3-v_4\right)y_{34}+j_{34}-j_E=v_4{y}_4$ (50d)

Rearranging gives

$j_W-j_{12}=\left(y_1+y_{12}\right)v_1-y_{12}{v}_2$ (51a)

$j_{12}-j_{23}=-y_{12}{v}_1+\left(y_{12}+y_2+y_{23}\right)v_2-y_{23}{v}_3$ (51b)

$j_{23}-j_{34}=-y_{23}{v}_2+\left(y_{23}+y_3+y_{34}\right)v_3-y_{34}{v}_4$ (51c)

$j_{34}-j_E=-y_{34}{v}_3+\left(y_{34}+y_4\right)v_4$ (51d)

The currents on the left hand side represent the total of the equivalent current sources directed into each node, $J_i=j_{i-1,i}+j_{i,i+1}$ . Thus, we can rewrite the equations in the form

$J_1=\left(y_1+y_{12}\right)v_1-y_{12}{v}_2$ (52a)

$J_2=-y_{12}{v}_1+\left(y_{12}+y_{23}+y_2\right)v_2-y_{23}{v}_3$ (52b)

$J_3=-y_{23}{v}_2+\left(y_{23}+y_{34}+y_3\right)v_3-y_{34}{v}_4$ (52c)

$J_4=-y_{34}{v}_3+\left(y_{34}+y_4\right)v_4$ (52d)

This can be written in matrix form

$\left[\begin{array}{c}\begin{array}{c}{J}_1\\ J_2\end{array}\\ \begin{array}{c}{J}_3\\ J_4\end{array}\end{array}\right]=\left[\begin{array}{cccc}{y}_1+y_{12}& -y_{12}& 0& 0\\ -y_{12}& y_{12}+y_{23}+y_2& -y_{23}& 0\\ 0& -y_{23}& y_{23}+y_{34}+y_3& -y_{34}\\ 0& 0& -y_{34}& y_{34}+y_4\end{array}\right]\left[\begin{array}{c}\begin{array}{c}{v}_1\\ v_2\end{array}\\ \begin{array}{c}{v}_3\\ v_4\end{array}\end{array}\right]$ (53)

The first matrix on the right-hand side of Equation (53) is termed the admittance matrix. This, multiplied by the column matrix of nodal voltages equals the column matrix of nodal current source values

$\left[J\right]=\left[Y\right]\left[V\right]$ (54)

Taking the inverse of the admittance matrix and multiplying by the nodal current sources then gives the voltages at the nodes.

$\left[V\right]={\left[Y\right]}^{-1}\left[J\right]$ (55)

These nodal voltages represent the voltages at the ends of each transmission line section and can be substituted into Equations (19) and (20) to give the voltage and current profile along the length of each section.

5. Example Modelling Results

To illustrate the application of the modelling technique and to examine the voltages produced in typical Ocean/Earth conductivity structures we present results for three models representing 1) a shallow sea, 2) a shallow ocean, and 3) a deep ocean.

Model 1. Shallow Sea

This model consists of a shallow sea 600 km wide, with a depth of 100 m, with land at either end as shown in Figure 9(a). The seawater section has the same properties as the continental shelf in Table 1. This is a simplified model approximation to the North Sea.

Model 2. Shallow Ocean

This model consists of a deep ocean section 1800 km wide, with a depth of 1 km, with continental shelf 100 km wide and 100 m deep on either side as shown in Figure 9(b). The total width from coast to coast is 2000 km. The parameters for each section are the same as those shown in Table 1, except the ocean section has a seawater depth of 1 km. This is a simplified model approximation of the Tasman Sea.

Figure 9. Example models for (a) shallow sea, (b) shallow ocean, (c) deep ocean.

Model 3. Deep Ocean

This model consists of a deep ocean section 7800 km wide, with a depth of 4 km, with continental shelf 100 km wide and 100 m deep on either side as shown in Figure 9(c). The total width from coast to coast is 8000 km. The parameters for each section are as shown in Table 1. This is a simplified model approximation of the Pacific Ocean.

Transmission Line Parameters

For each model the layer thicknesses and resistivities from Table 1 are used in Equations (1)-(4) to calculate the transmission line series impedance, Z, and parallel admittance, Y. These values are then used in Equations (5) and (6) to calculate the transmission line propagation constant, γ, characteristic impedance, Z₀, and adjustment distance, $\frac{1}{\gamma }$ for each section. The values obtained are shown in Table 2. The “Shallow Sea” parameters are used for the continental shelf sections in models 2 and 3 as well as for the seawater section in model 1. For all of the example models the transmission line parameters for the land section are used in Equations (27) and (28) to calculate the Norton equivalent circuit components.

The voltage profiles for the three example models are shown in Figure 10. These show a transition from the linear profile characteristic of an “electrically-short” transmission line, seen for the Shallow Sea model to the S-shaped profile characteristic of an “electrically-long” transmission line, seen for the Deep Ocean model.

6. Potentials at Boundaries

The above results show that the earth potentials peak at the coast which is the boundary between the land and sea regions with their different conductivity structures. To examine the coast potential more closely and to examine how this is affected by the width of a shallow sea or continental shelf we make calculations for different widths of the Shallow Sea model and for different widths of the continental shelf in the Deep Ocean model.

6.1. Dependence on Shallow Sea Width

The calculations for the Shallow Sea model (Figure 9(a)) were repeated for different widths of the sea. The land on either side was assumed to extend far enough from the coast that it can be represented by the equivalent circuits for an

Table 2. Transmission line parameters for the sections of the example models.

Figure 10. Voltage profiles for the three example models of the Ocean/Earth conductivity structure, calculated with electric field values of 0.1 V/km for the deep ocean and 0.3 V/km for the continental shelf, shallow sea, and land sections. (a) Shallow sea, 600 km wide and 100 m deep; (b) Shallow Ocean, 2000 km wide and 1 km deep; (c) Deep Ocean, 8000 km wide and 4 km deep. Note the different scales on the axes of these plots.

“electrically-long” transmission line. We assume values for the shallow sea and land electric fields of 0.3 V/km. The width of the continental shelf is varied from 0 to 2000 km and calculations made for the potential at the coast, U_C (Figure 11).

A shallow sea with land on either side can also be represented by the circuit of Figure 5, where the transmission line represents the shallow sea and the Thevenin equivalent circuits represent the land on either side. Where the width of the

Figure 11. Dependence of the coast potential on the width of a shallow sea for an electric field in the sea and the land of 0.3 V/km.

shallow sea is less than its adjustment distance then the transmission line model reduces to the circuit shown in Figure 6 and the coast potentials U_C are given by U₁ and U₂ from Equations (39) and (40).

If the width of the shallow sea is extended to (unreasonably) large values, it can be represented by an “electrically-long” transmission line and the coast potential can be calculated using the circuit in Figure 7 where region 1 represents the shallow sea and region 2 represents the land. The coast potential is then given by Equation (42). For an electric field of 0.3 V/km in both the shallow sea and the land and inserting the appropriate values from Table 2, U_C = 77.46 V. This represents a limiting value for the coast potential for the specified electric field and transmission line parameters.

6.2. Dependence on Width of the Continental Shelf

To examine how the coast potential is affected by the width of the continental shelf we repeated the calculations for the deep ocean model but with different continental shelf widths. The ocean width was set to 10,000 km. This is unreasonably large but makes sure that the deep ocean section is “electrically-long” so that the potentials at the coast under study are not affected by the coast at the other side of the ocean. We assume values for the deep ocean electric field of 0.1 V/km and continental shelf and land electric fields of 0.3 V/km. The width of the continental shelf is varied from 0 to 2000 km and the calculations made for both the potential at the coast, U_C and the potential at the edge of the continental shelf, U_E (Figure 12).

To understand why the coast potential (black curve in Figure 12) initially increases and then decreases as the width of the continental shelf increases, it is useful to consider the transmission line model shown in Figure 5. Here the Continental Shelf is represented by the transmission line in the centre of the model, with the Deep Ocean and Land represented by the Thevenin equivalent circuits at either end.

Figure 12. Dependence of the potentials at the coast, U_C, and at the edge of the continental shelf, U_E, on the width of the continental shelf.

Table 3. Thevenin circuits.

The Thevenin equivalent circuit components for the different regions are calculated using the adjustment distances and characteristic impedances from Table 2 and the electric field values. For the electric field values of 0.1 V/km in the deep ocean and 0.3 V/km in the continental shelf and land sections as used in the modelling for Figure 8, we obtain the parameters for the Thevenin equivalent circuits shown in Table 3.

The circuit in Figure 5 is used to examine two limiting cases: 1) when the continental shelf width is small so that it can be represented by an “electrically-short” transmission line, and 2) when the continental shelf is exceedingly wide so that it can be represented by an “electrically-long” transmission line.

For an “electrically-short” transmission line its equivalent-pi components are given by Equations (35) to (37) and the circuit of Figure 5 reduces to that shown in Figure 6. The current along the transmission line is then given by Equation (38) and the potential at the edge of the continental shelf (boundary between ocean and continental shelf) is given by Equation (39) and the potential at the coast (boundary between continental shelf and land) is given by Equation (40).

As the width of continental shelf is increased the modelling results in Figure 12 shows that U_C increase to a width of about 200 km and then decreases, while U_E shows a steady decrease. Both U_C and U_E approach limiting values for very large widths of the continental shelf. To determine these limiting values, we can use a model comprising just two regions where both extend far enough away from the boundary that they can be represented by Thevenin equivalent circuits for “electrically-long” transmission lines as shown in Figure 7.

Boundary between Continental Shelf and Land

For the boundary between a very wide continental shelf and the land, substituting the appropriate Thevenin circuit components from Table 3 into Equation (41) gives:

$I=\frac{154.92+77.46}{387.3+774.6}=0.2\text{A}$ (56)

Substituting this into Equation (42) gives

$U_C=154.92-0.2\times 387.3=77.46\text{V}$ (57)

Boundary between Deep Ocean and Continental Shelf

For the boundary between the deep ocean and a very wide continental shelf, substituting the appropriate Thevenin circuit components from Table 3 into Equation (41) gives:

$I=\frac{118.32+154.92}{84.5+387.3}=0.579\text{A}$ (58)

Substituting this into Equation (42) gives

$U_E=118.32-0.579\times 84.5=69.38\text{V}$ (59)

Now, with the circuit models of Figures 5-7, and the limiting cases considered above, we are able to explain the dependence of U_C and U_E on the width of the continental shelf shown in Figure 12. For the hypothetical case with a continental shelf of zero width the boundary between the deep ocean and the land has an earth potential of 99 V. For continental shelves of increasing widths (but still short enough to be considered “electrically-short”), the circuit of Figure 6 shows that the induced voltage in the continental shelf section adds to the voltages from the deep ocean and land sections to increase the current across the boundaries. Equation (40) shows that this increase in current increases the potential at the coast and, because of the large value of Z₂ (774.6 Ω), this has a noticeable effect on U_C. Equation (39) shows that the increase in current decreases the potential at the edge of the continental shelf but, because the value of Z₁ (84.5 Ω) is smaller than Z₂, this produces a comparatively smaller change in the value of U_E as shown in Figure 12.

The circuit of Figure 6 only applies for continental shelf widths that are small enough to satisfy the “electrically-short” condition for the transmission line modelling. In this condition the transmission line is short enough that there is no significant leakage of current from the conductive surface layers through the resistive layers underneath. When the continental shelf width is increased this condition no longer applies and some of the current in the surface layers leaks out through the resistive layer underneath so the current across the boundary starts to decrease and the potential at the coast starts to go down. Figure 12 shows that this decrease in U_C starts for continental shelf widths greater than 250 km.

When the continental shelf is much wider so that its transmission line model can be considered “electrically-long”, Equations (56) and (57) show that the current across the continental shelf/land boundary drops to a limiting value of 0.2 A and the potential at the coast goes to the limiting value of 77.46 V. At the boundary between the ocean and the continental shelf Equations (58) and (59) show that the current drops to a limiting value of 0.579 A and the potential at the edge of the continental shelf goes to a limiting value of 69.38 V.

7. Discussion

This transmission line modelling of geomagnetic induction in the Ocean/Earth conductivity structure was motivated by the need to calculate the variations in potential at the coasts at either end of a submarine cable. The power feed equipment at each end of a submarine cable will be connected to local earth which will have the coast potential, U_A, at one end of the cable, A, and coast potential, U_B, at the other end of the cable. During geomagnetic disturbances, a submarine cable will experience a voltage, V_C, comprised of two parts: the induced electromotive force, ${\mathcal{E}}_C$ , induced by the magnetic field variations directly in the cable and the potential difference between the ends of the cable

$V_C={\mathcal{E}}_C+U_A-U_B$ (60)

where: ${\mathcal{E}}_C={\displaystyle {\oint }_C\stackrel{\to }{E}\cdot \text{d}\hat{l}}$ with the integration taken along the path of the cable,

U_A is the Earth potential at the end “A” of the cable,

U_B is the Earth potential at the end “B” of the cable.

The calculations for a specific submarine cable will depend on the bathymetry of the sea or ocean and the conductivity structure beneath the seafloor and on the land at either end of the cable. Here we have described the calculation method and presented examples of simplified Ocean/Earth conductivity structures to illustrate how the modelling is used and the results that will be obtained.

The results show that the maximum potentials (positive or negative) occur at the ends of the transmission line, corresponding to the coastline on either side of the ocean. A key parameter is shown to be the “adjustment distance” that is the inverse of the propagation constant of the transmission line. For ocean widths that are more than four times the adjustment distance the end potentials are the product of the electric field and the adjustment distance. For a wide ocean, such as in example 3, the “electrically-long” criterion is met so the end potentials, in the absence of any conducting path through the land, would be the electric field value of 0.1 V/km multiplied by the adjustment distance of 1183 km, giving a value of 118.3 V. The model result (99 V) is slightly less because of the effect of the continental shelf and the current flow into the land.

Shallow seas tend to be “electrically-short” which means that the earth potential produced at one coast is influenced by the distance to, and conductivity structure, at the opposite coast. Therefore, modelling of geomagnetic induction in a shallow sea needs to include the whole sea from coast to coast, as well as the conductivity structure of the land on either side. Modelling of geomagnetic induction in a deep ocean needs to take account of the continental shelf on either side of the ocean, and it is shown that the width of the continental shelf has a significant influence on the potentials produced at the coast.

The results obtained are valid for models that meet the “thin sheet” condition. This is typically met for induction by longer period variations which penetrate through the seawater and crustal layers. Short period variations have a smaller skin depth so do not penetrate as far into the Earth. Other methods, such as finite element modelling (FEM) need to be used in such cases. However, the results presented here can provide a set of useful test cases for testing such modelling software.

Calculations of geomagnetic induction in submarine cables under the oceans involve a number of differences compared to calculations for long conductors, such as cables, pipeline and power systems, on land. The principal difference is that the geomagnetic disturbance is attenuated by the conducting seawater so produces smaller electric fields on the seafloor of a deep ocean than on land. The coast potentials are significant for all systems but on land, the potential gradients add to the induced electric field to increase the geomagnetically induced currents (GIC) in systems, while for submarine cables the coast potentials act against the induced emf in the cable and reduce the impact on the power feed for the cable.

Submarine cables are the longest conductors on Earth and extend for much larger distances than systems on land, thus the spatial characteristics of a geomagnetic disturbance have greater significance. For systems on land, initial calculations are made by assuming that the magnetic field variations are uniform across the region of the system under study, although this is recognised as an approximation and more work is being done to include the spatial structure of the disturbance in GIC calculations. However, the length of submarine cables means that the geomagnetic disturbances at one end of a cable may be significantly different from those at the other end. The example calculations presented in this paper have used a uniform disturbance to simplify the calculations. Future calculations for realistic disturbances need to take account of the differences in the disturbances along the whole length of the cable, however this will be challenging because knowledge of the geomagnetic field variation over the ocean is less constrained by measurements than over the land.

8. Conclusions

Geomagnetic induction in the Ocean/Earth conductivity structure can be modelled using the generalised thin sheet approach of Ranganayaki and Madden [15] which represents the surface layers of the Earth by a double layer with a conductive top part and a resistive bottom part. The conductive top part is given by the combined conductance of the seawater and sediment layers, while the resistive bottom part is given by the resistance through the crust and lithosphere mantle.

A generalised thin sheet can be modelled as a transmission line with series impedance, Z, given by the conductance of the top layer, and parallel admittance, Y, determined by the resistance of the bottom part [16] . This allows distributed source transmission line (DSTL) theory and network theory to be applied to calculate the Earth potentials produced by an induced electric field in the Ocean/Earth conductivity structure.

Transmission line theory allows the potentials within a transmission line section to be modelled in terms of the line’s propagation constant, γ, and characteristic impedance, Z₀, and the potentials at the ends of the line. A critical parameter is the inverse of the propagation constant referred to as the “adjustment distance”. Simplified calculations are possible for “electrically-short” transmission lines with length much less than the adjustment distance and for “electrically-long” transmission lines with length greater than four adjustment distances.

Complex structures can be modelled using network theory where transmission line sections are represented by their Norton equivalent-pi circuit involving a current source in parallel with an admittance which is combined to form a nodal admittance network. Using Thevenin’s current law for each node a matrix equation is formed relating the current sources [J] to the admittance matrix [Y] and the nodal voltages [U]. Matrix inversion of [Y] and multiplication by the current sources [J] then gives the nodal voltages [U].

Example calculations for a shallow sea (100 m depth), shallow ocean (1 km depth) and deep ocean (4 km depth) show the Earth potential variation along the seafloor and that they are largest at the coasts. For a shallow sea, the coast potentials increase with increasing width of the sea but asymptotically approach a limiting value for large widths. For a deep ocean, the coast potential varies with the width of the continental shelf, increasing to a maximum value at a width, for the example considered, of 250 km and then decreases to a limiting value for large widths of the continental shelf.

Open Science

The computational model, SCUBAS: Submarine Cables Upset by Auroral Streams is developed using Python to compute geomagnetic induction effect on submarine cables [8] . A pip installer is also available to install the model in any python environment. Example codes are also available to describe the parameters and illustrate the use of SCUBAS [8] . The majority of analysis and visualization was completed with the help of free, open-source software tools such as matplotlib [21] , IPython [22] , pandas [23] , and others (e.g., [24] . Our code is published in Zenodo repository [25] .

Acknowledgements

We thank Dr. Lidia Nikitina for useful comments on the manuscript. DHB was supported by Natural Resources Canada. MDH and SC were supported by NASA grant 80NSSC19K0907. XS is supported by NASA award 80NSSC19K0907 and the NASA DRIVE Science Center for Geospace Storms (CGS) under award 80NSSC22M0163. NRCan contribution number: 20220638.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Boteler, D.H., Pirjola, R. and Nevanlinna, H. (1998) The Effects of Geomagnetic Disturbances on Electrical Systems at the Earth’s Surface. Advances in Space Research, 22, 17-27. https://doi.org/10.1016/S0273-1177(97)01096-X
[2]	Gummow, R.A. (2002) GIC Effects on Pipeline Corrosion and Corrosion Control Systems. Journal of Atmospheric and Solar-Terrestrial Physics, 64, 1755-1764. https://doi.org/10.1016/S1364-6826(02)00125-6
[3]	Pulkkinen, A., Bernabeu, E., Thomson, A., Viljanen, A., Pirjola, R., Boteler, D., Eichner, J., Cilliers, P.J., Welling, D., Savani, N.P., Weigel, R.S., Love, J.J., Balch, C., Ngwira, M., Crowley, G., Schultz, A., Kataoka, R., Anderson, B., Fugate, D., Simpson, J.J. and MacAlester, M. (2017) Geomagnetically Induced Currents: Science, Engineering and Applications Readiness. Space Weather, 15, 828-856. https://doi.org/10.1002/2016SW001501
[4]	Rajput, V.N., Boteler, D.H., Rana, N., Saiyed, M., Anjana, S. and Shah, M. (2021) Insight into Impact of Geomagnetically Induced Currents on Power Systems: Overview, Challenges and Mitigation. Electric Power Systems Research, 192, Article ID: 106927. https://doi.org/10.1016/j.epsr.2020.106927
[5]	Ingham, M., Divett, T., Rodger, C.J. and Sigley, M. (2022) Impacts of GIC on the New Zealand Gas Pipeline Network. Space Weather, 20, e2022SW003298. https://doi.org/10.1029/2022SW003298
[6]	Patterson, C.J., Wild, J.A. and Boteler, D.H. (2023) Modeling the Impact of Geomagnetically Induced Currents on Electrified Railway Signaling Systems in the United Kingdom. Space Weather, 21, e2022SW003385. https://doi.org/10.1029/2022SW003385
[7]	Meloni, A., Lanzerotti, L.J. and Gregori, G.P. (1983) Induction of Currents in Long Submarine Cables by Natural Phenomena. Reviews of Geophysics, 21, 795-803. https://doi.org/10.1029/RG021i004p00795
[8]	Chakraborty, S., Boteler, D.H., Shi, X., Murphy, B.S., Hartinger, M.D., Wang, X., Lucas, G. and Baker, J.B.H. (2022) Modeling Geomagnetic Induction in Submarine Cables. Frontiers in Physics, 10, Article 1022475. https://doi.org/10.3389/fphy.2022.1022475
[9]	Marti, L., Yiu, C., Rezaei-Zare, A. and Boteler, D. (2014) Simulation of Geomagnetically Induced Currents with Piecewise Layered-Earth Models. IEEE Transactions on Power Delivery, 29, 1886-1893. https://doi.org/10.1109/TPWRD.2014.2317851
[10]	Boteler, D.H. and Pirjola, R.J. (2019) Numerical Calculation of Geoelectric Fields that Affect Critical Infrastructure. International Journal of Geosciences, 10, 930-949. https://doi.org/10.4236/ijg.2019.1010053
[11]	Gilbert, J.L. (2005) Modeling the Effect of the Ocean-Land Interface on Induced Electric Fields during Geomagnetic Storms. Space Weather, 3, S04A03. https://doi.org/10.1029/2004SW000120
[12]	Pirjola, R. (2013) Practical Model Applicable to Investigating the Coast Effect of the Geoelectric Field in Connection with Studies of Geomagnetically Induced Currents. Advances in Applied Physics, 1, 9-28. https://doi.org/10.12988/aap.2013.13002
[13]	Goto, T. (2015) Numerical Studies of Geomagnetically Induced Electric Field on Seafloor and Near Coastal Zones Incorporated with Heterogeneous Conductivity Distributions. Earth, Planets and Space, 67, Article No. 193. https://doi.org/10.1186/s40623-015-0356-2
[14]	Boteler, D.H. and Pirjola, R.J. (2003) The Magnetic and Electric Fields Produced in the Sea during Geomagnetic Disturbances. Pure and Applied Geophysics, 160, 1695-1716. https://doi.org/10.1007/s00024-003-2372-6
[15]	Ranganayaki, R.P. and Madden, T.R. (1980) Generalized Thin Sheet Analysis in Magnetotellurics: An Extension of Price’s Analysis. Geophysical Journal International, 60, 445-457. https://doi.org/10.1111/j.1365-246X.1980.tb04820.x
[16]	Wang, X., Boteler, D.H. and Pirjola, R. (2023) Distributed-Source Transmission Line Theory for Modeling the Coast Effect on Geoelectric Fields. IEEE Transactions on Power Delivery. https://doi.org/10.1109/TPWRD.2023.3279462
[17]	Pasyanos, M.E., Masters, T.G., Laske, G. and Ma, Z. (2014) Litho1.0: An Updated Crust and Lithospheric Model of the Earth. Journal of Geophysical Research: Solid Earth, 119, 2153-2173. https://doi.org/10.1002/2013JB010626
[18]	Boteler, D.H. (1997) Distributed-Source Transmission Line Theory for Electromagnetic Induction Studies. Proceedings of 1997 Zurich EMC Symposium, Zurich, 18-20 February 1997, 401-408.
[19]	Boteler, D.H. (2013) A New Versatile Model of Geomagnetic Induction of Telluric Currents in Pipelines. Geophysical Journal International, 193, 98–109. https://doi.org/10.1093/gji/ggs113
[20]	Boteler, D.H., Trichtchenko, L. and Edwall, H.E. (2013) Telluric Effects on Pipelines. Proceedings of the CEOCOR Symposium, Florence, 6-7 June 2013, 1-16.
[21]	Hunter, J.D. (2007) Matplotlib: A 2D Graphics Environment. Computing in Science & Engineering, 9, 90-95. https://doi.org/10.1109/MCSE.2007.55
[22]	Perez, F. and Granger, B.E. (2007) IPython: A System for Interactive Scientific Computing. Computing in Science & Engineering, 9, 21-29. https://doi.org/10.1109/MCSE.2007.53
[23]	McKinney, W. (2010) Data Structures for Statistical Computing in Python. Proceedings of the 9th Python in Science Conference, Austin, 28 June-3 July 2010, 56-61. https://doi.org/10.25080/Majora-92bf1922-00a
[24]	Millman, K.J. and Aivazis, M. (2011) Python for Scientists and Engineers. Computing in Science & Engineering, 13, 9-12. https://doi.org/10.1109/MCSE.2011.36
[25]	Chakraborty, S. (2023) shibaji7/submarine_cable_modeling: Geomagnetic Induction in Submarine Cables (V1.0). Zenodo. https://doi.org/10.5281/zenodo.7055321