A 3-D RBF-FD solver for modeling the atmospheric global electric circuit with topography (GEC-RBFFD v1.0)

. A numerical model based on radial basis function-generated finit differences (RBF-FD) is developed for simulating the global electric circuit (GEC) within the Earth’s atmosphere, represented by a 3-D variable coefficien linear elliptic partial differential equation (PDE) in a spherically shaped volume with the lower boundary being the Earth’s to-pography and the upper boundary a sphere at 60 km. To our knowledge, this is (1) the firs numerical model of the GEC to combine the Earth’s topography with directly approximating the differential operators in 3-D space and, related to this, (2) the firs RBF-FD method to use irregular 3-D stencils for discretization to handle the topography. It benefit from the mesh-free nature of RBF-FD, which is especially suitable for modeling high-dimensional problems with irregular boundaries. The RBF-FD elliptic solver proposed here makes no limiting assumptions on the spatial variability of the coefficient in the PDE (i.e., the conductivity profile) the right hand side forcing term of the


Introduction
The global electric circuit (GEC) is a system of currents within Earth's atmosphere.The system is define by the volume between two highly conductive shells, one for the surface of the Earth and the other for the lower ionosphere.These two highly conductive shells can be thought of as a leaky capacitor.The currents in the system are driven by elec-trifie clouds that produce a source current, which then holds the ionosphere at a fi ed potential relative to the Earth.Far away from storm clouds, into the so-called fair weather region, this potential difference between the ionosphere and ground produces an electric current that is ∼ 2 pA m −2 globally.This global return current can be measured by current probes and electric fiel mills on the ground to estimate its strength as well as global distribution of thunderstorms.
The firs modeling efforts focused on decomposing the domain into separate regions and solving the problem analytically with a spherical harmonics decomposition (Hays and Roble, 1979;Roble and Hays, 1979).However, in order to obtain solutions, these models needed to impose constraints on the source and conductivity distributions.Further advancements have focused on modeling the system with an electrical engineering approach of resistors and capacitors aligned in series and parallel (Rycroft et al., 2008;Odzimek et al., 2010).Other models have focused on how individual aspects of the system change, such as how aerosols and clouds influenc the resistivity within the domain and what effect that has globally on the solution (Tinsley and Zhou, 2006).All of these previous modeling efforts have had to either make assumptions on the solution or simplify the domain, omitting topography, to obtain a feasible solution.For an excellent overview of the GEC and recent progress made we refer readers to Williams and Mareev (2014).
The RBF-FD (radial basis function-generated finit differences) GEC (global electric circuit) model proposed here solves the full three-dimensional problem with the Earth's Published by Copernicus Publications on behalf of the European Geosciences Union.
V. Bayona et al.: GEC-RBFFD v1.0 real topography as the bottom boundary, without making any limiting assumptions on the conductivity (the coefficient of the partial differential equation (PDE)) or source distribution (the right hand side forcing term of the PDE).The structure of the paper is as follows.Section 2 introduces the PDE with its corresponding boundary conditions that will be discretized and solved.Section 3 gives a brief introduction to the RBF-FD discretization of differential operators with references for more in-depth study.Section 4 is the core of the paper, describing the numerical implementation.Section 5 gives a test case with a known solution for method validation, using an analytic conductivity profil (i.e., coefficient of the PDE).Section 6 builds on Sect.5, using the same conductivity profil but changing the forcing term to actual observational data.Section 7 is the hardest case in which discrete data is used for all inputs into the PDE, i.e., both coefficient and the right hand side forcing term.Lastly, Sect.8 gives some timing results, followed by conclusions.

Formulation
The 3-D electric potential for a given conductivity distribution and electrifie cloud current sources can be determined by the equation −∇ • (σ (r, θ, λ)∇u) = S(r, θ, λ), (1 where σ is the conductivity, u is the electric potential and S is the source distribution.This equation is derived by applying Ohm's law to the steady-state current continuity equation.The 3-D domain is define as where k(θ, λ) is the Earth's surface (i.e., topography) and r b is the altitude from sea level where the top boundary is enforced.In this paper, the mean radius of the Earth is set to r earth = 6400 km and r b = 60 km.As boundary conditions, zero electrostatic potential is enforced along the Earth's surface, and zero net current at the upper boundary, which leads to the potential where R is the global resistance and I s is the upward current at the top boundary generated by the electrifie clouds.Since Eq. ( 1) is linear, the electrostatic potential u can be split as u = u f + u s , where u f is the fair weather potential, and u s the source potential,

Integrated quantities
In this paper, we are also interested in integrated quantities derived from the fair-weather and source potential.Scaling the fair-weather potential in Eq. ( 4), as ûf = u f /RI s , leads to solving it with the boundary condition ûf r=r b = 1.Then, the atmospheric resistance R is computed from ûf as The upward current I s is computed from the solution of Eq. (5) as In both Eqs. ( 6) and ( 7), is the surface of the sphere that encloses the domain at the top boundary.As a result, the GEC solution is equal to u = RI s ûf +u s , from which we can compute the net current at the top boundary, and the net charge within the domain, where both quantities must be conserved.

RBF-FD method
The radial basis function-generated finit differences (RBF-FD) method can be considered a natural generalization of classical finit differences (FD) (Shu et al., 2003;Tolstykh and Shirobokov, 2003;Wright and Fornberg, 2006).As in FD, RBF-FD approximates a linear differential operator Lu at node x k ∈ R d as a linear combination of the values of the function at the n closest nodes, The main difference lies in how the differentiation weights w i are computed.While FD enforces Eq. ( 10) to be exact for

Infinitel smooth RBFs
polynomials evaluated at the node x k , RBF-FD enforces it for RBF interpolants: where φ(r) is a radial basis function, • is the Euclidean distance, and λ i are the RBF coefficients Some examples of smooth RBFs are listed in Table 1.Unlike FD, in which the interpolation problem is not guaranteed to be non-singular for scattered nodes in n dimensions (n ≥ 2), RBF-FD is guaranteed to be non-singular no matter how the n nodes (assumed distinct) are scattered in any number of dimensions (Fasshauer, 2007;Fornberg and Flyer, 2015b).
Augmenting RBF interpolants with polynomials can be beneficial In this work, MQ-RBF interpolants augmented with a constant are used, so that the constraint n i=1 w i = 0 can be satisfie and the solution exactly reproduces a constant (Lehto, 2012;Flyer et al., 2012Flyer et al., , 2015b;;Fornberg and Flyer, 2015b, a).This results in a less oscillatory interpolant and thus more accurate derivative approximations.Further augmentation with more polynomials is currently being studied in Flyer et al. (2015a) and Bayona et al. (2015).As a result, the system of equations that determines the RBF-FD differentiation weight w i to approximate Lu is The weight w n+1 is discarded after the system is solved.
Some of the main features of RBF-FD (Bayona et al., 2010;Bayona and Kindelan, 2013;Stevens et al., 2009;Lehto, 2012;Flyer et al., 2012Flyer et al., , 2015b;;Fornberg and Flyer, 2015b, a) have proven to be very beneficia in modeling the GEC.For instance, RBF-FD is a meshless local method that only depends on the Euclidean distance between neighboring nodes.From the implementation point of view, this feature makes the method independent of the number of dimensions and, as a result, it is straightforward to program even for three-dimensional domains such as the one considered in this work.In addition, RBF-FD approximations can achieve high-order accuracy, at the same time yielding highly sparse differentiation matrices.This is specially important when applied to this kind of elliptic problem where there are millions of unknowns and a large linear system of equations must be solved.

Numerical implementation
As described in Sect.2, the solution of the GEC model Eq. ( 1) is given by u = u f +u s , where u f = ûf RI s and u s are the solutions of Eqs. ( 4) and ( 5), respectively.The differential operator Lu = −∇ •(σ ∇u) might be numerically illconditioned due to the highly variable and exponential nature of the conductivity σ .To overcome this issue, it is possible to take advantage of the fact that σ > 0 and improve the conditioning by rewriting the PDE Eqs. ( 4) and (5) as and respectively.In the following subsections, the numerical implementation is explained in detail.

Change of variable
Figure 1 shows the Earth's topography used in the numerical model.The left side of the figur shows the height above sea level averaged for a 1.9 • ×2.5 • grid in latitude and longitude, with the actual scaling of the problem r ∼ r earth .Since the highest averaged region is 5 km, as seen in Fig. 1a, compared to r earth ≈ 6400 km this scaling misses all the topographical features with the Earth appearing flat In order to increase the topographical resolution of the model, a change of variable is considered: where A and B are constants determined by enforcing the conditions www.geosci-model-dev.net/8/3007/2015/Geosci.Model Dev., 8, 3007-3020, 2015 In RBF-FD modeling, there is a well-known tradeoff between accuracy and ill-conditioning, unless what is known as a stable algorithm is used (Larsson et al., 2013;Fornberg et al., 2013); however, these algorithms increase the computational cost by a factor of about 10.The method itself suffers from numerical ill-conditioning for small values of εh, where ε is the shape parameter and h the internodal distance.In order to achieve the best accuracy and avoid ill-conditioning, the RBF shape parameter ε must be selected for every resolution h.
To make the method attainable to the scientifi community and overcome the necessity of selecting ε for varying h, we have used an alternative approach in this work.We propose to take advantage of the change of variable and select the computational domain for every resolution such that where N H is the number of nodes in the latitude-longitude angular direction, N r is the number of nodes in the radial direction and h H and h r are the angular and radial internodal distances, respectively.As a result, the extent of the computational domain changes for every N H and N r , but εh r and εh H are fi ed and independent of the resolution.The condition number is also fi ed and thus the problem of selecting the shape parameter is bypassed.It can be selected once and used for any resolution.Thereby, the accuracy of the solver is also fi ed, in this work at slightly greater than fourth-order for the resolutions considered, assuming the variable coefficien σ of the PDE is analytic.However, in more realistic applications σ comes from discrete data that is never more than C 1 and thus the accuracy of the solver in such cases is a moot point.

Spatial discretization
Spatial discretization is similar to a nested shell model (Wright et al., 2010;Flyer and Fornberg, 2011).The majority of the domain is discretized horizontally by using a spherical shell formed by N H icosahedral nodes and radially by repeating this spherical shell from sea level to the top boundary at every spacing of h r .This results in N r radially aligned spherical shells of N H nodes.However, to incorporate topography, the following alterations need to be incorporated.
1.An algorithm has been developed to distribute nodes along the topography, with approximately twice as many nodes on land as over the oceans to accommodate the steeper gradients of the orography (see Fig. 2a).
2. When part of a spherical shell intersects land, the nodes that fall under the Earth's topography are discarded.For example, the firs spherical shell above sea level is at an altitude of 500 m.In the top left panel of Fig. 4, corresponding to r = 500 m in the test case, the white areas show the topography and thus where the nodes have been discarded.
3. The last item to be done is to rearrange the nodes where a shell intersects land in order to have quasi-uniformly distribution in that region.Thus, nodes on each shell lying above the surface are repelled in the latitudelongitude direction (using a charge-type repulsion algorithm) while holding the nodes on the boundary fi ed; this allows the nodes near the Earth's surface to follow the topography more closely, yet keeping the radial distance between nodes fi ed, and preserving conditioning of the matrix system to be solved.
Consequently, there are two different regions in terms of the structure of the node layouts and thus the shape of the stencils used to approximate the differential operator on the left hand side of Eq. ( 15).A near-surface region formed by the  nodes close to the topography (< 8 km), where the differential operator Eq. ( 15) at any node is approximated by Eq. ( 13) using the closest 56 nodes in 3-D space (found via a k-d tree search) and forming a true 3-D stencil.In contrast, above all topography (i.e., > 8 km), the nodes retain their nested shell formation, resulting in a 2-D + 1-D stencil formation.A hybrid FD/RBF-FD approach is implemented, where classical 5-node FD approximations in the radial direction (1-D) are combined with 21-node RBF-FD approximations for the angular derivatives (2-D).
To compute the RBF-FD differentiation weights for the 3-D Laplacian and gradient that appear in the PDE, the system of Eq. ( 13) must be solved for at each node in the domain.By using the chain rule and taking the derivatives with respect to the square of the Euclidean distance (the argument of the RBF), it is possible to write the action of the differential operator on the RBF in a very elegant way.If d = x − x j 2 is the square of the Euclidean distance between an RBF centered at the node x j = {r j , θ j , λ j } and evaluated at x = {r, θ, λ}, the three-dimensional Laplacian and gradient can be written in the scattered-node region as and In the structured nested-shell region above 8 km, the angular terms of the differential operators can be written following the procedure described in Wright et al. (2010), where the author noticed that the approximation is invariant under rotations.In this case, the square of the Euclidean distance takes the form d = 2 (1 − sin θ) and the surface Laplacian s on a spherical shell can be written as Since the same spherical node layout is repeated in the structured region, the angular derivatives are computed once on a unitary sphere and scaled by 1/r(ξ ) 2 for every layer, where r(ξ ) is the radius of each layer.

Handling topography: ghost nodes
The stencils that incorporate boundary nodes will be more one-sided and might have skew shapes due to terrain.The weights that approximate the differential operators on those stencils might lose some properties, such as the positiveness in the case of the Laplacian.As a consequence, the stability of the numerical solver may be affected.The spectrum of the eigenvalues can behave oddly as the shape parameter decreases, with some eigenvalues crossing the imaginary axis and the differentiation matrix becoming unstable.Naturally, these eigenvalues do not have physical meaning and are only a numerical artifact.In order to avoid this issue, the concept of "ghost nodes" is implemented.The name comes from the fact that these nodes are used in Eq. ( 13) to approximate the differential operator on near-boundary/boundary nodes, making the stencils more symmetric, but no equations are ever enforced at these nodes as they are outside the domain.For most boundary nodes, a ghost node is introduced directly outside the domain, i.e., under the topography or directly above 60 km (the only caveat to this is when the terrain becomes to steep, as in the Andes or Himalayas, making the ghost nodes close to overlapping; in these cases a smaller one-sided stencil is used to maintain stability of the solver).
To determine the value of the function at the ghost nodes, the PDE is enforced on the boundary, in addition to the Dirichlet boundary conditions.Hence, the resulting system of equations has as many unknowns as equations and preserves uni-solvency.In addition, the interior stencils near the irregular boundary recover a more symmetrical shape and the stability of the solver improves.This procedure is also used at the top boundary to enable the use of 5-node stencils in the radial direction.

The elliptic solver
Once the differentiation weights are computed according to Eq. ( 13), they are assembled into a matrix that approximates the left hand side of the PDE.Each row of the matrix represents the discretized PDE at a single node.The left panel of Fig. 3 shows the sparsity pattern of the assembled matrix for a 4 • × 1 km resolution (a total of ≈ 154 000 nodes).One immediately notices that two different stencils have been used.The rows in the upper left corner where the pattern is much denser and unstructured corresponds to the 3-D 56node stencils below 8 km.The rest of the matrix, with its pentadiagonal-type pattern, corresponds to the 21 + 5-node stencil in the structured region above 8 km.Drastic changes in the pattern of a matrix generally impede iterative solvers, making for much poorer and slower convergence.The panel on the right of Fig. 3 shows the same matrix but after applying reverse Cuthill-McKee reordering, giving a consistent sparsity pattern with a nicer bandwidth for the iterative solver.Even though only 0.016 % of the entries are nonzero, the bandwidth of the matrix is about 5N H ≈ 12 000 nodes, much too large to use Gaussian elimination (i.e., the "\" operator in MATLAB).However, due to the high sparsity of the resulting differentiation matrix, the iterative method GMRES(20) (generalized minimal residual method) proved ideal for solving the linear system of equations, where the restarting parameter is set to 20.In order to greatly decrease the number of iterations necessary for convergence, GMRES must be preconditioned with a solution that results from a simplifie version of the PDE yet captures the main features.Since the conductivity increases exponentially with altitude and varies by orders of magnitude less in the angular directions, a simple exponential conductivity profil σ (r) = σ 0 e r/c is a good firs approximation and leads to a spatial operator, that can be solved extremely fast.As a result, Eq. ( 22) leads to a good preconditioner which is (1) repeatedly called a function fil by the GMRES solver for the original problem and (2) itself solved with GMRES(20) using incomplete LU factorization as preconditioner.The residual tolerance is set to 10 −9 and the maximum number of outer iterations is set to 10. Runtimes for solving the problem at different resolutions are listed in Sect.8.

Conductivity
The spatially varying conductivity σ appears as the variable coefficien in the PDE.It can be an analytic function or a discrete data set output from a different numerical model, such as the Whole Atmosphere Community Climate Model (WACCM; for further details about how the conductivity is computed see Baumgaertner et al., 2013Baumgaertner et al., , 2014)).In the latter case, the conductivity is interpolated to the node distribution used in this paper.

Sources
The storm counts based on the 2-D TRMM (Tropical Rainfall Measuring Mission) satellite data (Liu et al., 2008) provide the radial current density averaged over a 12-year period at 20 km altitude on a 1 • × 1 • grid between −38 and 38 • in latitude.To approximate the 3-D source term of the PDE from 2-D data, two approaches are proposed.

Dipole approach
the firs approach is to distribute 3-D dipoles over the Earth according to the spatial distribution of the data, each one with charge centers at altitudes r p i and r n i from the surface, where I ± is the dipole's current strength, ρ (θ − θ i , λ − λ i ) is the orthodromic distance from the charge center (θ i , λ i ), and a and b determine the widths of the dipoles in the radial and angular directions, respectively.The source term is then represented as the sum of these dipoles: The previous approach is commonly used in the literature to charge the ionosphere (Tzur and Roble, 1985;Mallios and Pasko, 2012).However, questions concerning the actual dipole's width, the manner in which the dipoles are distributed with respect to the TRMM satellite data, and on how to assign to the current strength I ± may arise.For the purpose of assessing the model with inputs obtained directly from such 2-D satellite data, it might be more natural to use the following approach, where the sources are enforced as a boundary condition which is always one dimension less than the PDE.

TRMM-BC approach
in this alternative approach, problem ( 5) is replaced by where J r (θ, λ) is the radial current density at 20 km provided by the TRMM satellite data and the domain of the problem is define as This approach is independent of dipole parameters within the model.
In Sects.6 and 7, both the dipole and TRMM-BC approaches are considered and compared in Tables 2 and 3.For the fair weather potential, u f = RI s ûf , the two approaches for source treatment yield field that only differ by the scaling factor I s .Thus, to observe the relative spatial variations in the fair weather field only one approach needs to be considered.

A test case for method validation
Before considering more realistic cases in terms of conductivity and source terms, a simplifie test case with a known analytic solution is proposed to validate the numerical scheme.A simple exponential conductivity profil that varies only in the radial direction (a good firs approximation to real atmospheric conductivity) is considered,   where c = 6 km and σ 0 = 5 × 10 −14 s m −1 .In this case, Eq. ( 15) then reduces to the 3-D problem where the domain = {(r, θ, λ) : −π/2 ≤ θ ≤ π/2, −π ≤ λ < π, k (θ, λ) ≤ r ≤ r b } and the function k (θ, λ) is the Earth's topography shown in Fig. 2b under the change of variables (Eq.16) with β = 1.This is an extreme stretching of the topography; however, it will allow for sharp gradients and more skewed stencils under 8 km that will test the robustness of the solver.The functions f and g are computed by assuming the exact solution Due to the exponential stretching direction, the resolution in the physical domain is variable in the radial direction, from 500 m close to the topography to 2 km near the top boundary, with N r = 60.In the angular direction the resolution is approximately 4 The solution is given in the left column of Fig. 4 at 500 m, 5 km, and 20 km above sea level and the error in the right column.At 500 m, the firs layer of nodes above the boundary, all RBF-FD stencils that approximate the differentiation operator on the right hand side of Eq. ( 27) involve boundary nodes and about 30 % of those involve nodes that lie directly on land surfaces.The importance of this is that the 3-D nearboundary stencils are more skewed and irregular, leading to a degradation of diagonal dominance in the matrix in Fig. 3.This in turn would have the expectation of decreasing accuracy.However, as can be seen in the error for 500 m the solution is accurate over most of the domain to O(10 −5 ).In fact, the error at 5 and 20 km is almost identical (the latter being slightly larger due to a coarser radial resolution from the exponential stretching), showing that the numerical treatment of the boundary has not impacted the accuracy of the solver.The highest errors are at the poles due to the solution having the steepest gradient as the poles are approached.

Forcing the PDE with observational data: analytic conductivity profil
Given the validation of the numerical scheme in the previous section, a natural progression for model performance would be to now force the PDE with observational data sources for more realistic modeling, yet using the same exponential conductivity profil as in Sect. 5.The forcing or source term corresponds with TRMM satellite data for the month of April at 12:00 UT.In the dipole approach, dipoles are spatially distributed according to the TRMM satellite data shown in Fig. 5; however, the strength of the 3-D current sources cannot be accurately represented since there is no information in the radial direction from the TRMM data.In this case, r n = 8 km and r p = 15 km, with I ± = ±4.2A has been chosen.This gives a Wilson current of 1 A at 20 km.The radial and horizontal widths are 1.5 and 150 km, respectively.In order to resolve them numerically, a numerical resolution smaller than the dipole's width is required.As a result, a resolution of 0.5 km × 0.75 • is used, which results in 9 108 837 nodes (N H = 73 962 and N r = 120).In the second approach (TRMM-BC), the data is directly implemented as the lower boundary condition, previously discussed in Sect.4.6.

Dipole approach
When the dipole approach is used to spatially approximate the current density distribution from the TRMM data, the integrated quantities obtained are listed in Table 2. Notice that the global total resistance for the case without topography is R = 233 , which is 10 larger than the case with topography.This is expected as the column of air above sea level is decreased when topography is included.The source current at the top boundary is I s = 1325 A and is independent of the topography.It charges the ionosphere and generates a potential difference equal to 308.7 and 295.7 kV for the cases without and with topography, respectively.In both cases, the net current at the top boundary I top and the net charge Q within the domain are numerically conserved, as shown in Table 2.

TRMM-BC approach
For the purpose of assessing the model with inputs obtained directly from TRMM satellite data, the alternative TRMM-BC approach proposed in Sect.4.6 is used, where the TRMM satellite data is directly enforced as a boundary condition.The corresponding integrated quantities are listed also in Table 2.The global resistance does not change from that calculated with the dipole approach because it is independent of the source term.However, the source current I s at the top boundary is 368 A larger than in the case based on dipoles.To alleviate this discrepancy, one could use the TRMM-BC approach to scale appropriately I ± in the dipole approach.The TRMM-BC approach also numerically conserves I top .Since there is a lack of knowledge of the charge centers, as explained in Sect.4.6, Q cannot be computed when using the TRMM-BC approach.
Figure 6 displays the convergence rate when computing the integrated quantities R and I s for both the analytic conductivity profil (Eq.26) and the discrete model data profil that will be used in the next section.Since the conductivity of the atmosphere naturally varies by orders of magnitude more in the radial direction than in the angular directions, convergence is considered with respect to N r , where the reference solution is set to N r = 320.As can be seen, when the variable coefficien of the PDE is an analytic function, a slightly greater than fourth-order convergence of the elliptic solver is achieved, as expected (see Sect. 4.1).

Fair weather field
Figure 7 shows the fair weather potential distribution u f at 1.5 km (top figure and 6 km (bottom figure above sea level.The effect of the topography on the GEC modifie the column resistance over the higher elevations, as noted earlier. Figure 8 shows the fair-weather current density (J r , J θ , J λ ) at 20 km above sea level.In the top figure notice that the larger radial current density is localized over the higher elevations.Furthermore, the topography also modifie the horizontal current density, especially at higher elevations, as can be appreciated in the middle and bottom panels of Fig. 8. Notice that in the horizontal components, (J θ , J λ ), the positive fl w of current is immediately neighbored by a negative fl w.However, the strength of the horizontal components are 2-3 orders of magnitude smaller than the radial one, with much of the Earth close to or at sea level having near-zero current density in these directions.

Forcing the PDE with observational data: model data conductivity profil
The last step to illustrate the robustness of the RBF-FD model is to consider model data with steep gradients as opposed to a smooth analytic function for the conductivity, using the two approaches for treating the TRMM source data.In lieu of an analytical exponential conductivity profile we now consider a conductivity profil from model data computed with WACCM (https://www2.cesm.ucar.edu/working-groups/wawg), as described in Baumgaertner et al. (2014).This discrete conductivity profil includes aerosols and fair-weather clouds as well as topography, varying not only in the radial direction but also in latitude and longitude.An example of its radial profil at 0 • latitude and 180 • longitude versus the exponential profil can be seen in Fig. 9, noting how cloud layers centered at 2 and 13 km cause steep gradients in the conductivity.As in the previous section, the results of the RBF-FD solver will be evaluated by noting whether they are consistent with physical expectations, both with regard to integrated quantities as well as to the fairweather fields In the fina subsection we will show a full solution of the GEC.

Integrated quantities
The resistance R is of course independent of the sources and thus will be the same for the approaches of source treatment.Although clouds and aerosols are known to increase the atmospheric column resistance, as seen in Fig. 10, Table 3 shows that the total integrated resistance is lower than what resulted when the exponential conductivity profil was used (see Table 2).However, it is important to note that the color bars in Fig. 10 use the same color map, with shades of blue being a lower column resistance than green.Thus, from the mid-latitudes to the polar regions, the column resistance is lower with the WACCM conductivity profile with the lowest values at 16.4 m 2 as opposed to 16.8 m 2 and thus resulting in a lower total integrated resistance.
In contrast, the source current I s at the top boundary differs dramatically between the two approaches in Table 2 due to the fact that the TRMM-BC approach does not incorporate any of the conductivity profil below 20 km, the region in which the conductivity is severely altered by cloud layers as was seen in Fig. 9.In fact, using the exponential or WACCM conductivity with the TRMM-BC approach makes little difference in I s (1693 A as opposed to 1712 A).With regard to the discrepancy between the two approaches for calculating the potential difference at the top boundary, the TRMM-BC approach can be used to scale the dipole source approximation to achieve the same values for u top , as noted in Sect.6.1.The net current at the top boundary I top for both approaches and the net charge Q within the domain for the dipole approach are numerically conserved, as shown in Table 3.
As noted earlier, Fig. 6 also displays the convergence rate when computing the integrated quantities R and I s for the discrete WACCM model data profile Since the conductivity data is only C 0 , one cannot expect greater first-orde convergence from any numerical method.As can be seen in Fig. 6, this is achieved.

Fair weather field
With regard to the fair weather fields the best way to examine the output of the RBF-FD solver is to see if the results are consistent with our expectation of what the physics should be.The two cloud layers at 2 and 13 km directly modify, at those altitudes, the radial electric fiel and current density as shown in the top row of Fig. 11, causing a jump in the field which would be expected.In the angular directions, one would expect the field to bend around the cloud layers, increasing the divergence of those field between the cloud layers (Baumgaertner et al., 2014).This can in fact be seen by the increased bulge between the cloud layers in the plots of the latitudinal and longitudinal components of the electric fiel and current density, shown by the middle and bottom rows of Fig. 11.

An example of the full solution for the GEC
Until now, we have only illustrated results calculated from 3-D fair weather potentials, u f , or those calculated from the 3-D source potential u s .To view a full solution (u = u f + u s ) of the GEC model based on the 3-D RBF-FD solver, we plot in Fig. 12 3-D isosurfaces of the radial current density J r = σ ∇u • r corresponding to ±3.5 pA m −2 .This isosurface was chosen since it shows the clearest image of the structure of currents within the GEC, specially with regard to topography.All panels in the figur are centered on the North Pole (NP). Figure 12a shows both the positive (red) currents fl wing upward as thin column currents from the Earth to the cloud layer as well as from the cloud layer to the ionosphere.The empty "ring" region corresponds with negative currents inside the cloud layer as shown in Fig. 12b.This fig ure also shows the downward currents to the high elevations such as the Rocky mountains, Greenland, and the Himalayas.Figure 12c shows the combined isosurfaces 3.5 pA m −2 (red) and −3.5 pA m −2 (blue) of the previous two plots.

Timing results of the model
In order to give the reader a feel for how long it takes to solve the GEC model with the RBF-FD elliptic solver, Table 4 shows some run-times at different resolutions for the WACCM conductivity with clouds, which is the computationally most intense.All test cases were conducted on a MacBook Pro 2.7 GHz Intel Core i7.The code was written  in MATLAB and run under version 2013a and used a peak of 5 GB of memory for the calculations.The results under the GMRES(20) column show the number of iterations and computing time that it takes to solve the resulting system of equations.Speeding up the algorithm is possible through parallel implementation of the method.This would require scaling GMRES across multiple CPUs or GPUs (Li and Saad, 2013), taking into consideration careful partitioning of the nodes to ensure proper load balancing across processors.The latter can be done using a domain partitioning library such as ParMETIS (http://glaros.dtc.umn.edu/gkhome/metis/parmetis/overview).However, implementation of RBF-FD for scattered node layouts on parallel computing architectures is a novel topic, only having been addressed since 2012 (Bollig et al., 2012;Erlebacher et al., 2014;Tillenius et al., 2015).

Conclusions
This paper advances the research front in two different fields First, it presents a novel numerical elliptic solver based on RBF-FD that can handle irregular boundaries, as the Earth's topography.This required novel developments, such as 1. an algorithm for node distribution on the Earth's surface; 2. a repelling algorithm to maintain quasi-uniformity of nodes where stencils intersect the boundary; 3. a novel spatial discretization scheme that consists of two types of stencils, one to handle the irregular near topography regime below 8 km and the other the regular regime above 8 km; 4. strategies to combat loss of accuracy near boundaries and maintain stability of the solver; 5. a preconditioner especially designed to aid the elliptic solver due to the drastic change in the sparsity pattern of the matrix from the use of two completely different types of stencils.
On the atmospheric science front, the new solver is the firs to make no limiting assumptions on the inputs to the PDE, including geometry.For instance, in previous numerical models, the surface of the Earth has been assumed spherical.By ignoring topography within the domain, the total resistance was off by 10 .This modifie resistance effects where the currents in the domain fl w.The higher elevation regions have a lower column resistance and therefore more current was able to return in the fair weather GEC at these locations.With complete fl xibility of model inputs, two different approaches for the treatment of current sources, given 2-D satellite data at 20 km, were also developed to solve for currents and electric field within the Earth's atmosphere.The fir t approach involved placing 3-D dipoles globally to represent individual thunderstorms, where the satellite data provided only the spatial distribution but not the strength of the current charges.In contrast, the alternative method implemented the current density strength from the 2-D satellite data as a boundary condition.It was shown that the latter method, although giving better integrated quantities, as the upward current at the top boundary (I s ), might suffer from the lack of knowledge of the conductivity distribution near the sources which are below 20 km.Therefore, to determine the effect that different conductivity distributions have on the GEC one should utilize dipole sources within the model, and then scale the current and potential at the upper boundary by the approach that directly uses the satellite data.
To conclude, this novel solver allows for complete fl xibility of model inputs and thus will further investigations of the currents and electric field arising through different physical perturbations to the GEC.With higher fidelit data sets being produced by global climate models and even real data, one needs to utilize a solver that will couple these parameters without any limitations or assumptions.

Figure 1 .
Figure 1.Earth's topography: (a) actual scale of the problem (the color bar is the grid-averaged altitude above sea level in kilometers).(b) Result of the change of variable (Eq.16) selecting ξ 0 and ξ b as in Eq. (18), where N H = 18 566, N r = 70, h H = h r = 0.05 and β = 0.05.

Figure 2 .
Figure 2. Example of the discretization of the Earth's topography under the change of variable (Eq.16): (a) β = 0.05 and ∼ 150 km resolution at sea level.(b) β = 1 and ∼ 400 km resolution at sea level.

Figure 4 .
Figure 4. Left: numerical solution at different distances from the origin.Right: corresponding error.

Figure 5 .
Figure 5. Radial current density (pA m −2 ) at 20 km from sea level obtained through TRMM satellite data for the month of April at 12:00 UT.

Figure 6 .
Figure 6.Convergence rate when computing the integrated quantities R and I s using the analytical conductivity (Eq.26) and the conductivity computed from discrete model data used in Sect.7. The dash-dot lines represent the order of convergence.

Figure 7 .
Figure 7. Fair weather electric potential (kV) along the 1.5 km (top) and 6 km (bottom) constant height surface above sea level.In the top figure the regions in white are the intersection of the 1.5 km constant height surface with the Earth's topography.

Figure 9 .
Figure 9. Analytical-exponential conductivity used in Sect.6 compared to an example of the conductivity from WACCM model output with clouds and aerosols at 0 • latitude and 180 • longitude.

Figure 12 .
Figure 12.North Pole view of the Earth showing the isosurfaces 3.5 pA m −2 (red) and −3.5 pA m −2 (blue) for the global solution of the GEC model.

Table 2 .
Integrated quantities for the exponential conductivity (Eq.26) using both methods in Sect.4.6 with and without topography.

Table 3 .
Integrated quantities for WACCM conductivity with clouds, aerosols, and topography.

Table 4 .
Runtime results for the WACCM conductivity with clouds on a MacBook Pro 2.7 GHz Intel Core i7.