A 3-d Rbf-fd Elliptic Solver for Irregular Boundaries Printer-friendly Version Interactive Discussion a 3-d Rbf-fd Elliptic Solver for Irregular Boundaries: Modeling the Atmospheric Global Electric Circuit with Topography a 3-d Rbf-fd Elliptic Solver for Irregular Boundaries Printer-friendly Version

A numerical model based on Radial Basis Function-generated Finite Differences (RBF-FD) is developed for simulating the Global Electric Circuit (GEC) within the Earth's atmosphere , represented by a 3-D variable coefficient linear elliptic PDE in a spherically-shaped volume with the lower boundary being the Earth's topography and the upper 5 boundary a sphere at 60 km. To our knowledge, this is (1) the first 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 first RBF-FD method to use irregular 3-D stencils for discretization to handle the topography. It benefits from the mesh-free nature of RBF-FD, which is especially suitable for modeling high-dimensional problems 10 with irregular boundaries. The RBF-FD elliptic solver proposed here makes no limiting assumptions on the spatial variability of the coefficients in the PDE (i.e. the conductivity profile), the right hand side forcing term of the PDE (i.e. distribution of current sources) or the geometry of the lower boundary.


Introduction
The global electric circuit (GEC) is a system of currents within Earth's atmosphere.The system is defined by the volume between two highly conductive shells, one the surface of the Earth and the other the lower ionosphere.These two highly conductive shells can be thought of as a leaky capacitor.The currents in the system are driven by electrified clouds that produce a source current, which then holds the ionosphere at a fixed 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 field mills on the ground to estimate its strength as well as global distribution of thunderstorms.Introduction

Conclusions References
Tables Figures

Back Close
Full The first 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 influence 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 GEC model proposed here solves the full three dimensional problem with the Earth's real topography as the bottom boundary, without making any limiting assumptions on the conductivity (the coefficients of the 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 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 profile (i.e.coefficients of the PDE).Section 6 builds on Sect.5, using the same conductivity profile 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 coefficients and the right hand side forcing term.Lastly, Sect.8 gives some timing results, followed by conclusions.Introduction

Conclusions References
Tables Figures

Back Close
Full 2 Global electric circuit model

Formulation
The 3-D electric potential for a given conductivity distribution and electrified cloud current sources can be determined by the equation 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 defined 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 electrified 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,

Conclusions References
Tables Figures

Back Close
Full and u s the source potential,

Integrated quantities
In this paper, we are also interested in integrated quantities derived from the fairweather 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 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 can be computed the net current at the top boundary, and the net charge within the domain,

Conclusions References
Tables Figures

Back Close
Full where both quantities must be conserved.

RBF-FD method
Radial Basis Function-generated Finite Differences (RBF-FD) can be considered a natural generalization of classical Finite 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 the 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 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).RBF interpolants can be augmented with polynomials to increase accuracy.In this work, MQ-RBF interpolants augmented with a constant are used, so that the constraint n i =1 w i = 0 can be satisfied, allowing the solution to exactly reproduce a constant which increases accuracy (Lehto, 2012;Flyer et al., 2012Flyer et al., , 2015;;Fornberg and Flyer, 2015b, a).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 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., , 2015;;Fornberg and Flyer, 2015b, a) have proven to be very beneficial in modeling the GEC.For instance, RBF-FD is a mesh-less local method that only depends on the Euclidean distance between neighboring nodes.This feature makes the method independent of the number of dimensions, and as a result, it is straightforward to program even for threedimensional 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 ill-conditioned due to the highly variable and exponential nature of the conductivity σ.To overcome this issue, it is pos-Introduction

Conclusions References
Tables Figures

Back Close
Full sible to take advantage of the fact that σ > 0 and improve the conditioning by rewriting the PDEs (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.On the left side of the figure, it is shown 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 and β is a parameter which controls the topography stretching.Under this change of variable, the Earth is mapped over a sphere of radius ξ 0 and the radial coordinate is exponentially stretched, as shown in the right side of Fig. 1.
In RBF-FD modeling, there is a well-known trade-off between accuracy and illconditioning, unless what is known as a stable algorithm is used (Larsson et al., 2013;Introduction Conclusions References Tables Figures

Back Close
Full   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 scientific 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 fixed and independent of the resolution.The condition number is also fixed 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 fixed, in this work at slightly greater than fourth-order for the resolutions considered, assuming the variable coefficient σ 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 mute 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 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

GMDD Introduction Conclusions References
Tables Figures

Back Close
Full 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 first 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 latitude-longitude direction (using a charge-type repulsion algorithm) while holding the nodes on the boundary fixed; this allows the nodes near the Earth's surface to follow the topography more closely, yet keeping the radial distance between nodes fixed, 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 Introduction

Conclusions References
Tables Figures

Back Close
Full 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

GMDD Introduction Conclusions References
Tables Figures

Back Close
Full 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 singular.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 unisolvency.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 56 node stencils below 8 km.The rest of the matrix, with its pentadiagonal-type pattern, correspond to the 21 + 5 node stencil in the structured re-Introduction

Conclusions References
Tables Figures

Back Close
Full gion 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 re-ordering, 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.\ operator in MATLAB).However, due to the high sparsity of the resulting differentiation matrix, the iterative method GMRES(20) 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 simplified 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 profile σ(r) = σ 0 e r/c is a good first 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 as a function file by the GMRES solver for the original problem and (2) itself solved with GMRES(20) using incomplete LU factorization as preconditioner.Runtimes for solving the problem at different resolutions are listed in Sect.8.

Conductivity
The spatially-varying conductivity σ appears as the variable coefficient 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 (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.

GMDD Introduction Conclusions References
Tables Figures

Back Close
Full

Sources
The storm counts based on the 2-D TRMM 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 first 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 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 (Eq.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 defined 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 the fair weather potential, u f = RI s ûf , the two approaches for source treatment yield fields that only differ by the scaling factor I s .Thus to observe the relative spatial variations in the fair weather fields 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 simplified test case with a known analytic solution is proposed to validate the numerical scheme.A simple exponential conductivity profile that varies only in the radial direction (a good first 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- where the domain 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 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 profile
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 profile 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 can not 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) is the resolution.In the second approach (TRMM-BC), the data is directly implemented as the lower boundary condition, previously discussed in Sect.4.6.

GMDD Introduction Conclusions References
Tables Figures

Back Close
Full

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 profile (Eq.26) and the discrete model data profile that will be used in the next section.Since the conductivity of the atmosphere

Conclusions References
Tables Figures

Back Close
Full naturally varies by orders of magnitude more in the radial direction then 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 coefficient of the PDE is an analytic function slightly greater than fourth-order convergence of the elliptic solver is achieved, as expected (see Sect. 4.1).

Fair weather fields
Figure 7 shows the fair weather potential distribution u f at 1.5 km (top figure) and at 6 km a.s.l.(bottom figure).The effect of the topography on the GEC modifies 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 a.s.l.In the top figure, notice that the larger radial current density is localized over the higher elevations.Furthermore, the topography also modifies the horizontal current density, especially at higher elevations, as it can be appreciated in both the middle and bottom panels of Fig. 8. Notice that in the horizontal components, (J θ , J λ ), the positive flow of current is immediately neighbored by a negative flow.However, the strength of the horizontal components are two to three 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 profile
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 profile from model data computed with the Whole Atmosphere Community Climate Model (WACCM, https://www2.cesm.ucar.edu/working-groups/wawg),as described in (Baumgaertner et al., 2014).This discrete conductivity profile includes aerosols and fair-weather clouds Figures

Back Close
Full as well as topography, varying not only in the radial direction but also in latitude and longitude.An example of its radial profile at 0 • latitude and 180 • longitude vs. the exponential profile 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 the fair-weather fields.
In the final 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 profile 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 profile 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 can not expect greater first-order convergence from any numerical method.As can be seen in Fig. 6, this is achieved.

Fair weather fields
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 field and current density as shown in the top row of Fig. 11, causing a jump in the fields which would be expected.In the angular directions, one would expect the fields to bend around the cloud layers, increasing the divergence of those fields 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 field and current density, shown by the middle and bottom rows of Fig. 11. 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 type of stencils.

An
On the atmospheric science front, the new solver is the first 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 modified resistance affects where the currents in the domain flow.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 flexibility of model inputs, two different approaches for the treatment of current sources, given 2-D satellite data at 20 km, was also developed to solve for currents and electric fields within the Earth's atmosphere.The first 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 flexibility of model inputs and thus will further investigations of the currents and electric fields arising through different physical perturbations to the GEC.With higher fidelity 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.from radial basis functions, J. Comput. Phys., 212, 99-123, doi:10.1016/j.jcp.2005.05.030, 2006. 3528 Wright, G. B., Flyer, N., and Yuen, D. A.: A hybrid  Full Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Fornberg
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

)
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 • with N H = 2562 nodes per shell.The 3-D domain contains a total of 168 601 nodes (149 856 interior domain nodes, 6966 nodes on the bottom boundary (topography), 6655 ghost nodes under the topography, 2562 nodes for the top boundary and another 2562 nodes above the top boundary).The solution is given in the left column of Fig. 4 at 500 m, 5, and 20 km a.s.l. and the error in the right column.At 500 m, the first 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 near boundary 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 decreasing accuracy.However, Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 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 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 figure are centered on the north pole (NP).Figure 12a shows both the positive (red) currents flowing 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 figure also shows the downward currents Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Figure 1 .
Figure 1.Earth's topography: (a) actual scale of the problem (the colorbar 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.