Justification of the Two-Dimensional Model of Electroconductivity for the Earth’s Ionosphere

Conventional two dimensional model for electric fields in the Earth’s ionosphere is analyzed to estimate its error. The main difficulties arise due to asymmetry of the conductivity tensor. We use the energy method and small parameter expansion. To make it possible in spite of asymmetry of the tensor coefficients the problem is reduced to the problem of minimum of proper quadratic energy functional. The variational principle is stated and proved for the 3-D boundary value problem. The error of the 2-D approximation is analyzed in the case, when conductor occupies a flat layer 0 < z < z0 and is homogeneous in z direction, and the vector of magnetic field has only z component. The results of numerical simulation of the electric field penetration from ground to the Earth’s ionosphere with reduction of the 3-D model of the ionospheric conductor to the 2-D model are presented. Precision of such an approach is demonstrated.


Introduction
The tensor of electric conductivity in the Earth's ionosphere is asymmetric, and conductivity in the direction of magnetic field is much larger than other components. So, mathematical models of large scale electric fields and currents in the Earth's ionosphere are usually based on the assumption of infinite conductivity in the direction of magnetic field. Then magnetic field lines are equipotentials, and the model of electric potential distribution is two-dimensional.
A correction to such two-dimensional potential can be found as a solution for the equation that is stated in [12]. Typical values of the terms of this equation are estimated there to demonstrate that the correction is small. Electric fields in low latitude ionosphere are calculated in [18] in frames of 3 − D and 2 − D models, and it is shown that the results are close to each other as far as their 2 − D features are concerned. This paper is devoted to justification of the twodimensional approximation with estimation of the error. A simple case is under analysis. The conductor occupies a flat layer 0 < z < z 0 , it is homogeneous in z direction, vector of magnetic field has only z component.
We use the energy method [15] and expand the solution as a power series of a small parameter ε that is the ratio of transversal and field-aligned conductivities. It works for the problems with asymmetric tensor coefficients because problems of this kind are reduced to the problems of minimum of proper quadratic energy functionals in [2,3]. Since the 3-D boundary value problem under analysis is a specific one, the variational principle is stated and proved in the section 4.
The 2-D model is used for numerical simulation of the electric field penetration from ground to the Earth's Ionosphere in the last section. Its accuracy is shown for real height distributions of the coefficients of conductivity σ which are not constant in contrast with such a restriction in our mathematical proofs.

Statement of the 3-D boundary value problem
In a three-dimensional domain Ω occupied by a conductor, electric field strength e and current density j satisfy the equations div j = q curl e = g j = σe.
The components of the conductivity tensor σ and right-hand sides q, g are given functions of Cartesian coordinates x, y, z.
Here we analyze the domain Ω that is a flat layer 0 < z < z 0 . We suppose periodicity in x, y directions with periods a, b j x | Since we also suppose periodicity of all given functions q, g, σ all components of the vector functions e, j are periodical functions as a result of the conditions (2).
(4) We do not write the limits for integrals when integrate over the whole domain or the whole interval.
Introduce an auxiliary functionẽ For the difference e −ẽ the right-hand side in the second equation (1) becomes (0, 0,g z ), wherẽ Interchange differentiation and integration and use the Newton-Leibnitz formula for the function g z (x, y, z): As the result of the first constraint (4) the integrand equals zero. Hence the new right-hand sidẽ g z (x, y, z) = g z (x, y, 0) does not depend on z. The right-hand side in the first equation (1) changes also. To avoid new designations, we suppose without loss of generality that it is done beforehand and so from the very beginning g x = g y = 0, g z = g z (x, y).
To get a generalized solution for the problem (1-3), suppose that the right-hand sides in (1) are bounded in the norm of the space L 2 (Ω) ∫ q 2 dΩ < q 2 0 , These constraints may be replaced with the less rigid ones that are so called divergency constraints [15]. However, we use more rigid constraints to construct an approximate solution. The function q and its derivatives with respect to x, y are assumed to be bounded Here c 3 , x 0 are positive constants. Since magnetic field is vertical, and the medium without magnetic field is isotropic, the tensor σ is invariant with respect to rotation around z axis. Hence it may be written as Here σ P , σ H , σ ∥ , σ C = (σ 2 P + σ 2 H )/σ P are referred as Pedersen, Hall, field-aligned and Cowling conductivities. The Hall parameter β = σ H /σ P is a measure of asymmetry.
Let us introduce a small parameter ε so that is of the same order of magnitude as σ P . Unbounded increase of σ ∥ is equivalent to ε → 0. Denote the matrix σ withσ ∥ instead of σ ∥ asσ. It is independent of ε. The coefficients of the matrixσ are assumed to be uniformly bounded and the symmetrical part ofσ is assumed to be uniformly positive definite. These properties ofσ may be written as follows: where c 1 is some positive constant. However we use more rigid conditions to construct an approximate solution. We additionally suppose σ to be independent of z ∂σ/∂z = 0, (12) and the second derivatives σ P ,σ ∥ with respect to x, y are assumed to be bounded: where c 2 , c 0 , x 0 are positive constants. The new problem, which solution gives the solution to the problem (1)(2)(3), is stated below in the section 4. The original problem has a solution if the new problem has, but the uniqueness must be proved independently.

Uniqueness of the classical solution
Suppose that there are two smooth solutions for the problem (1-3). Then their difference e, j satisfies these equations with the zero right-hand sides q = 0, g = 0.
We do not take into account a small parameter that means ε = 1 in (9) in this section. We do not use the additional assumptions (7,12,13) either.
In the proof we use the function V : Such a function exists, since the vector field e is irrotational by the second equation (1) and g = 0.
By assumption e is not identically zero. Therefor we can find a number ζ > 0 and a point such that |e| = 2ζ.
In view of continuity of e we can choose a neighborhood of this point in wich |e| > ζ. Denote the volume of this neighborhood by x 3 1 . Consider the following integral over the domain Ω:

36
Justification of the Two-Dimensional Model of Electroconductivity for the Earth's Ionosphere Using the denotations (14), we can rewrite the integral as ω = We can replace the matrix σ T in this integral with its symmetrical part, since we calculate the quadratic form By positive definiteness of (σ +σ T )/2, (10,11), the integrand is nonnegative. Therefore, the integral over the domain Ω is not less than the integral over the selected neighborhood; hence, Now we turn to our integral written down in the shape (15). Transform identically the integrand: The first term equals V div (σ e) = V div j = 0, since the first equation (1) holds with q = 0. Transform the remaining integral by using of the Gauss -Ostrogradskii theorem: Returning to (14) we obtain In accordance with the last equation (1) it means By the boundary condition (3) the integrand equals zero at the horizontal parts of the boundary z = 0, z = z 0 . By the first periodicity condition (2) the integrals over the surfaces x = 0 and x = a are equal in magnitude and opposite in sign. Hence their sum equals zero. The same about the sum of integrals over the surfaces y = 0 and y = b. We obtain ω = 0 that contradicts (16). Hence, the assumption that e is different from identical zero is false. Thereby the uniqueness of the classical solution for the problem (1-3) is proved.

Variational principle
When the normal component of current density is given at the whole boundary, the equations of electroconductivity (1) are reduced to minimization of the quadratic energy functional in [3]. Here we give analogous statements and proofs for the problem (1-3).
We do not take into account a small parameter that means ε = 1 in (9) in this section. We do not use the additional assumptions (7, 12, 13) either.
One can obtain a new statement that has a symmetrical operator by one of the next two heuristic ideas. It is possible to introduce potentials similar to the usual scalar and vector potentials and to use them not separately but in pair. It is also possible to find the adjoint operator and multiply the original one by the adjoint operator from the right [3]. It is equivalent to introducing of some potentials. Another way is to use the least square method that is equivalent to multiplication of the original operator by the adjoint operator from the left. In contrast with that method we do not differentiate the right-hand sides and the solutions for our problem are more smooth functions. Our way corresponds to the reducing of Cauchy -Rieman equations with nonzero right-hand sides to a pair of the Poisson equations in 2-D space. These equations are separate for the main boundary value problems and the solution of one of them is identical zero because of zero right-hand side. That is why only one potential is usually used, and the same is traditionally done in more general case, when a pair of potentials could be useful.
Omitting details of the presented heuristic ideas which are similar to those in [3], we consider some quadratic functional and prove that its minimization gives pair of potentials that permits to construct the solution for the original problem, which uniqueness is proved above.

The energy space
We consider pairs f, p of smooth functions satisfying the following conditions Consider the following symmetrical bilinear form Here u, v and f, p are pairs of smooth functions which satisfy (17)(18)(19), and First, consider the auxiliary integral Transform the integrand: The second summand vanishes identically. Using the Gauss -Ostrogradskii formula, we can transform the rest integral, coming to ∫ f curl n p dΓ.
The integrand equals zero at the boundaries z = 0, z = z 0 by (18). By the first and the third periodicity conditions (17) the integrals over the surfaces x = 0 and x = a are equal in magnitude and opposite in sign. Hence their sum equals zero. The same about the sum of integrals over the surfaces y = 0 and y = b. Therefore the integrals (23) and (22) are equal to zero.
Consider the quadratic form to which the bilinear form (20) is reduced when (u, v) = (f, p). The matrix in (20) is degenerate since its upper blocks result from multiplying its lower blocks by −σ. After adding the duplicated integral (22), the matrix of the quadratic form under the integrating becomes where I is the identity matrix. Since the integral (22) equals zero, the value of the quadratic form is not changed. The matrix K is symmetrical. Hence it has 6 real eigenvalues. Its specific shape can be obtain by substitution of σ (8) and s (21). Without changing its eigenvalues we can put the matrix in a block diagonal form by interchanging rows and columns. The blocks are ( , (24) The eigenvalues of the first and the second blocks are equal. They can be obtained as the roots of the quadratic equation By Viete's theorem λ 1 λ 2 = 1, λ 1 + λ 2 = σ C + 1/σ P . In view of (10) we obtain In view of (11) more rigid restrictions are valid for the rest diagonal elements in (24). Hence all eigenvalues of the matrix K satisfy (25).
Since the integral (22) equals zero, the obtained estimations of the eigenvalues of the matrix K lead to the following estimations from below and above The inequality c 1 ≥ 1 that is a sequence of (10) or (11) was used.
To continue the estimations (27) from below and (26) from above, we consider the functions f and p separately.
Since the average value of the function f equals zero (19), the Poincare's inequality is valid for it. For the parallelepiped Ω it can be written in the shape [15]: For the function p we consider the identity (curl p) 2 + (div p) 2 = (grad p x ) 2 + (grad p y ) 2 + (grad p z ) 2 + div (p div p − (p grad) p), that can be easily verified in Cartesian coordinates. Let integrate it over Ω. In view of the Gauss -Ostrogradskii theorem the integral of the last term is equal to the flux of the following vector through the boundary of the domain p div p − (p grad ) p.
The normal component of the vector equals zero at the boundaries z = 0 and z = z 0 since the vector p has the only nonzero z− component by (18). The x− component equals Only tangential derivatives are involved. They satisfy periodicity conditions as a result of (17) as well as the vector p itself. So the values at the surfaces x = 0 and x = a are equal in magnitude and opposite in sign. Hence their sum equals zero. The same about the sum of the integrals over the surfaces y = 0 and y = b.
Thereby we obtain the identity for the functions p which satisfy (17,18). It permits to transform the right-hand side of (27) Since the average value of the function p z equals zero (19), it satisfies the Poincare's inequality with the same constant as one for the function f (28): The function p x equals zero at z = 0 (18). Hence To square the identity and use the Cauchy -Bunyakovskii inequality By integration over the domain Ω we obtain (32) Since the function p y also equals zero at z = 0 by (18) it satisfies the same inequality

38
Justification of the Two-Dimensional Model of Electroconductivity for the Earth's Ionosphere The inequalities (28, 31, 32, 33) permit to continue the estimation (30) (34) Such an inequality means the quadratic form (20) is positive definite in the space L 2 (Ω). In couple with (30) they mean positive definiteness in the space W (1) 2 (Ω). Since the quadratic form that corresponds to the bilinear form (20) is positive definite, we can use the bilinear form as the inner product in the set of pairs f, p of smooth functions satisfying (17 -19). It is referred to as the energy inner product.
By introducing some inner product we have constructed some Hilbert space. We complete it by adjoining limit elements. The resultant space is referred to as the energy space.
The inequality (26) in view of the identity (29) is the estimation of the quadratic form from above with the sum of norms of the functions f, p x , p y , p z as the elements of the space W

The energy functional
We define the energy functional (35) The linear functionals can be estimated by the Cauchy -Bunyakovskii inequality By inequalities (28, 30 -33) the right-hand sides can be estimated from above by the energy norm with factors q 0 , g 0 (6) which are independent of f, p. Hence linear functionals in (35) are bounded ones.
By the Riesz's theorem each bounded linear functional is representable by means of the inner product into some unique element of the Hilbert space. Therefore we can write the energy functional (35) as where f 0 , p 0 is some element of the energy space. We can transform the expression by distinguishing the square of the difference: The second term is independent of f, p and the first one is positive definite. Therefore the minimum of W (f, p) is attained at f = f 0 , p = p 0 . Since the element f 0 , p 0 belongs to the energy space and is defined uniquely, we have proved that the energy functional has a unique minimum in the energy space.

A weak solution
We use the functions f, p obtained as the result of the energy functional minimization to define the following vector functions Here we use e, j only as a contraction of these expressions, but later we will demonstrate that it is the sought solution for the original problem (1-3).
The minimum condition for the energy functional can be written down, on use made of notations (37), as the following identity that must be hold for all smooth functions u, v satisfying conditions (17)(18)(19): Take where the function V is obtained as the solution for the following boundary value problem for the Poisson equation in the parallelepiped Ω: The function ξ is chosen below. It can be easily verified that the solution for the last problem exists and is unique. By the first boundary condition (39) Therefore the obtained pare of the functions u, v satisfy all the conditions (17 -19) independently of the function ξ, which must be only enough smooth. For such a pair u, v the identity (38) takes the form Transform the second integral by using of the Gauss -Ostrogradskii theorem The function g satisfy the first condition (4). Hence the only integral over the boundary remains. In view of the first boundary condition (39) V = 0 at the boundaries z = 0 and z = z 0 . The sum of the integrals over Computational Research 1(2): 34-51, 2013 39 the boundaries x = 0 and x = a is the integral over y, z of the difference It equals zero because of V and g x are periodic functions. The sum of the integrals over the surfaces y = 0 and y = b equals zero by the same reason. Thereby from (40) we obtain the equality ∫ ξ div p dΩ = 0.
By using the arbitrariness of the function ξ, we can routinely deduce that div p = 0. Assuming the contrary, we take a point Ω at which div p ̸ = 0 and take a neighborhood of this point in which div p preserves sign. Take a function ξ that differs from identical zero, is nonnegative in this neighborhood and equals zero outside the neighborhood. Then the integral (41) is nonzero, that contradicts the identity.
Suppose additionally that all functions are smooth. Then the weak equality div p = 0, that is a result of the minimization of the energy functional, holds pointby-point: div p = 0.
After integration by parts using (42) the identity (38) is omitted because its integrand is identical zero by (18).
When v = 0 and u equals zero at the boundaries the only first integral remains in the identity. By using arbitrariness of the function u we can routinely deduce div j − q = 0.
By the same way we obtain curl e = g.
Then only integrals over boundaries remain in the identity (43). Each of them produces the independent identity.
Arbitrariness of u at the boundaries z = 0 and z = z 0 gives j z = 0 there. Other identities give periodicity conditions of the form (2).
In this way the pare of the functions f, p providing a minimum to the energy functional enable us to construct by formulae (37) the solution to the original boundary value problem (1-3), if we additionally assume their smoothness. In general case we have got the weak solution from the point of view of the identity (38) for arbitrary element u, v of the energy space.
The conditions of a minimum to the energy functional for the functions f, p themselves can be obtained by substitution (37) to the identity (38) that define the weak solution f, p. Under assumption of smoothness substitution of (37) to the equations (1-3) produce the boundary value problem that is completed with the boundary conditions (17 -19).
The boundary conditions (2) for the functions f, p can be simplified taking into account the periodicity of the tangential derivatives that is a result of (17) Since unit determinant of the matrix ( , the first pairs of these conditions can be resolved. It is already done in (44). The boundary conditions (17 -19) are referred to as the main ones because the functions, at which the minimum is searching, must satisfy them. In contrast the boundary conditions in (44) are natural ones since they are satisfied as a result of minimization.
The converse assertion is also easily provable: the solution to the boundary value problem (44, [17][18][19] provides minimum to the energy functional in the energy space. In this way the existence and the uniqueness of the weak solution for the original boundary value problem (1-3) is proved and the minimum principle for the energy functional is justified.
Since the functions f, p belong to the energy space the function f and the components of the vector function p belong to the space W Each functions f, p satisfying (17 -19) can be uniquely represented as a sum of elements of these sets. The function V is obtained as the solution for the Poisson equation with the boundary conditions (45), and the function p 1 is the difference p − grad V . And conversely the sum of arbitrary elements of these sets satisfies (17 -19).
Substitution of the elements f, p 1 and 0, grad V to the bilinear form (20) gives zero integrand because div p 1 = 0 and grad 0 = 0, curl grad V = 0. Hence the inner product equals zero.
We complete these two sets by adjoining limit elements. The resultant spaces are two orthogonal subspaces of the energy space.
The functions p 1 in the elements f, p 1 of the first subspace have zero divergence in weak sense. Hence The second subspace that contains only elements of the shape 0, grad V is in some sense excessive. The vector functions e, j which are of interest in the original problem do not change if elements of this subspace are added. So it is possible to consider the first subspace as the whole energy space slightly changing the statement of the problem. This simplifies the proofs but creates difficulties in constructing numerical algorithms. For this reason we avoid doing so.
The energy functional can be rewritten as the sum of two functionals corresponding to the two orthogonal subspaces. To the second subspace there corresponds only 1 2 since the linear part of the energy functional on the functions (f 2 , p 2 ) = (0, grad V ) is equal to the second integral in (40) that equals zero as it was shown. Thus we can independently minimize the energy functional in the two orthogonal subspaces constructed above.

Thermodynamics
In view of notations (37) with s chosen in (21), we can write down the quadratic form that corresponds to the energy inner product (20) in terms of the original boundary value problem (1)(2)(3): The product e T j is the density of the Joule dissipation, i.e., the thermal energy production that accompanies electric current.
As it was shown in the previous section, we can minimize the energy functional in the subspace which elements have zero divergence ∫ (div p) 2 dΩ = 0.
Then the quadratic form (46) equals the total Joule dissipation in Ω, which justifies usage of the term "energy".
Under the choice instead of (21), we obtain another symmetrical statement for the original problem. Here θ ia an arbitrary symmetrical tensor that leaves the tensor s uniformly positive definite. In this event, In particular, on taking θ to be the product of the identity tensor by the scalar function 1 / T , where T is the absolute temperature, we give this integral the meaning of the total entropy production in the domain Ω.
Thus, the energy norm makes sense from the standpoint of nonequilibrium thermodynamics.
Consider the energy functional in the shape (36). Since the second term (36) is independent of f, p, the minimization of the energy functional implies the minimization of the first term, that is the square of the energy norm of the difference between an arbitrary element f, p and the exact solution. By (46) it is the minimization of the functional or the functional ∫ δe T θδj dΩ, if s is chosen by more general formula (47). Here δ stands for the difference from the exact solution. The quantity δe may be called as the electric field fluctuation. In fact, it is impossible to calculate the integral (48), since the exact solution is unknown, but its minimization is equivalent to the minimization of the energy functional (35), in which all integrals can be calculated for arbitrary f, p.
Thus, from the viewpoint of thermodynamics the square of the energy norm is the entropy production, and the exact solution can be obtained by minimization of the energy of fluctuations.

Small parameter expansion
For the limit value of the small parameter (9) ε = 0, the operator of the boundary value problem (44, [17][18][19] ceases to be positive definite. At the same time the solution becomes two-dimensional. Zero value of ε corresponds to σ ∥ = 0 by (9). To keep finite current density we must put e z = 0 because of the last equality Computational Research 1(2): 34-51, 2013 41 (1) and (8). Then the x− and y− components of the second equation (1) with zero right-hand sides (5) cause e x and e y independence of z. Since e z = 0 the whole vector e is independent of z when ε = 0.
We fulfil the small parameter expansion of the solution under additional restrictions, which are due to the method of proof but not to the essence of the problem. We are to suppose conductivity independece of the height z (12). By the same reason from the very beginning we consider the flat conducting layer and suppose magnetic field to be vertical, when conductivity tensor has shape (8).
Represent the solution in the form Denote the set (u, v, w) as the vector u.
It will be shown that it is enough to put the conditions ∫ r dz = 0, to make this representation unique. Use the main boundary conditions (17 -19) separately for the two-dimensional functions and for the three dimensional correction in the representation (49). Hence all functions are periodic over x, y, and are also satisfied as a result of (50). Substitute the representation (49) to the energy functional (35). The integral that contains the products of the derivatives of two-dimensional functions F, P and those of three-dimensional corrections r, u is the integral over x, y and over z of the following sum Every summand is a product of the functions independent of z because of (12) and two-dimensional definition of F, P (49) and the functions with zero integrals over z by (50) or (53). Hence the integral over z of every summand equals zero and the whole integral can be omitted.
All integrals of the shape (22) can be also omitted because they are equal zero by the same reasons as (22). Therefore after some transformation using independence σ of z (12) the energy functional takes the shape where By (4) and (5)  The integral of rq equal zero by (50) is put into the expression of W q (r) for convenience.
As it can be seen by (55), the functional is factorized, and the two-dimensional functions F, P can be obtained separately. If they are smooth the minimization of W 0 (F, P ) is equivalent to the solution of the equations with conditions (52). Here the operators are twodimensional ones, Curl P = (∂P/∂y, −∂P/∂x).
By the formulae we obtain the periodical solution for the two-dimensional problem Div J =q(x, y), that simulate two-dimensional electric fields and currents. Similar to (57, 52) statements for the first and for the second boundary value problems for the equations (59) are analyzed in [2]. Here we do not include the proof for the periodicity conditions because the differences in the proofs are obvious. State the results only. The operator (57, 52) is symmetrical and positive definite. There exists the unique pair of the functions F, P that satisfies (52) and minimizes the value of the functional W 0 (F, P ). It is the weak solution for the problem (57, 52). By the formulae (58) it gives the periodical solution for the equations (59). The squared functions F, P are summable as well as their first derivatives. The squared vector functions E, J obtained by the formulae (58) are also summable.
The three-dimensional correction r, u to the solution of the two-dimensional problem ought minimize the rest terms of the sum (55).
It will be shown that the minimization of W 3 (r, u) over u results also W 4 (u, v) = 0, that means minimization of this nonnegative term. Since the function u is not involved to other terms of the sum (55), as a result of the minimization of W 3 (r, u) we obtain the expression of u in terms of the function r that is fixed during this minimization over u.

Minimization of the energy functional over u
The function r is fixed in this section. First prove that div u = 0 as a result of minimization of W 3 (r, u) over u. Take into account the set of the functions where the function u minimizes the value of W 3 (r, u), t is an arbitrary number, V is a solution for the problem (39) with an arbitrary function ξ. Then Since the function u minimizes the value of W 3 (r, u), the factor after t equals zero. It means the identity for an arbitrary V from the set under analysis. In view of (39) we obtain the identity with an arbitrary function ξ ∫ ξ div u dΩ = 0.
Hence div u = 0 and it is possible to reject the last term in the functional W 3 (r, u). Now the minimum conditions for W 3 (r, u) over u can be written as (60) and An arbitrary function c(x, y) appears in the last equation (61) because the restriction (50) for the function w. Show that c(x, y) = 0. Integrate the third equation (61) over z between the limits 0 and z 0 taking into account σ P and β independence of z and (60) All terms in the left-hand side equal zero by (53, 50). Hence c(x, y) = 0.
The first and the second equations (61) are the sums of the derivatives over z. Therefore they can be integrated Here ξ, η are arbitrary functions. Prove that they are zero. Integrate these equalities over z between the limits 0 and z 0 taking into account σ P and β independence of z All terms in the left-hand side equal zero by (53, 50). Hence ξ(x, y) = η(x, y) = 0.
Now deduce the equation for the function w from (62, 63). Multiply both equations (62) by the factor σ P . Differentiate the results with respect to x and y respectively and sum them. Add to and subtract from the left-hand side ∂w/∂z The sum of the first three terms equals to −∂ div u/∂z that equals zero by (60). Therefore the equation for w takes shape It defines the unique function w taking into account the last condition (54) and ∂w/∂z| z=0,z0 = 0, that is a result of (60)and (53). However we need the function w that satisfy more strict condition (50). Prove that the obtained function satisfies it. Integrate the equation (64) over z.
The right-hand side equals zero by the first condition (50). The second term in the left-hand side equals zero by (65). Hence it is two-dimensional Laplace equation for the function ∫ w dz. Only a constant may be its periodical solution and it equals zero by (54). Hence the second condition (50) is satisfied.
The function r is fixed during the construction of the functions u, v, w and w is already constructed. It only remains to integrate the equalities (62) over z taking into account (63) and the boundary conditions (53) at z = 0: The conditions (53) at z = z 0 are satisfied, because the integrals of the functions w, r between the limits 0 and z 0 equal zero and σ H is independent of z.
By (62, 63) we obtain the expressions of x−, y− components of the vector curl u, and by (66) its z− component equals zero: This formula provides zero value for W 4 (u, v). The obtained expressions (67, 60) permit to express the minimal over u value of W 3 (r, u) in terms of r: In sum with ε 2 W 2 (r) it equals After the performed minimization (55) over F, P and u it remains to find the function r that satisfies (50) and minimizes the functional εW r (r) = εW 1 (r) − εW q (r) + ε 2 W xy (r).
This way we obtain the minimum of the energy functional W (f, p) (55) over all arguments.

Minimization of the energy functional over r
The minimum conditions for W r (r) (68) over periodic functions r satisfying the first condition (50) under additional suggestion of smoothness is the following boundary value problem: The last condition is added for uniqueness of the solution for the problem. Demonstrate that it can be used instead of the first condition (50). Integrate the first equation over z between the limits 0 and z 0 . Taking into account that σ P is independent of z, ∂r/∂z = 0 at z = 0 and z = z 0 and the definition ofq(x, y) (56) we obtain −Div (σ P Grad where the common factor ε is omitted. The operators are two-dimensional ones. This is a two-dimensional elliptical equation for the function ∫ r dz. By (10) the coefficient σ P is uniformly bounded and strictly positive. The average value of the unknown function equals zero by the last condition (69). Therefore the unique periodical solution is zero ∫ r dz = 0.
Point out that the equation (70) for the function ∫ r dz is a particular case of the equation (57) with σ H = 0. Then the system of the equations (57) is split into independent equations for F and P , the first of which differ from (70) only by nonzero right-hand sides. The obtained equality (71) means that the solution for the problem (69) also satisfies the first condition (50).
The solution for the problem (69) exists and is unique, since for an arbitrary positive ε this is the Neumann boundary value problem for an elliptical equation with uniformly bounded and uniformly positive definite symmetrical coefficient tensor.
By the boundary condition (69) the function r can be formally expanded in Fourier series The Fourier series starts with n = 1 since r 0 (x, y) = 0 by (71).
By this expansion the problem (69) is reduced to the set of two-dimensional problems with a small parameter ε: −ε Div (σ P Grad r n (x, y)) + The periodic over x, y solution must be obtained for each n.

44
Justification of the Two-Dimensional Model of Electroconductivity for the Earth's Ionosphere The shape of the right-hand side is simplified in comparing with the first equation (69) because the only nonzero term of the Fourier series of the functionq(x, y) isq 0 (x, y) =q(x, y), and by (56) q 0 (x, y) −q 0 (x, y) = 0. It corresponds to r 0 (x, y) = 0.
Integrate the equation (73) over x, y. Use the Gauss -Ostrogradskii formula for the first term. Taking into account the periodicity of the involved functions it equals zero. There remains the equalities of the average values for each Fourier coefficient ∫ ∫σ The last condition (69) and the first condition (50) are satisfied in spite of nonzero average values of r n (x, y), because of that the z− dependencies are of the shape cos(nπz/z 0 ) and n ≥ 1.
We construct the solution for each singular pertubed equation (73), following the Lusternik -Vishik method [19]. We expand the unknown function in a power series of the small parameter ε. Since the boundary conditions for the equation (73) are the only periodicity over x, y, the assymptotic expansion does not contain boundary layers. Represent where δ n (x, y) is the remainder term that ought be estimated.
The equations for the functions ρ n (x, y), δ n (x, y) we derive from (73) by equating to zero the sum of the terms which do not contain ε and the sum of the rest terms: δ n (x, y) − ε Div (σ P Grad δ n (x, y)) = Div (σ P Grad ρ n (x, y)) .
By (10) the coefficient σ P is strictly positive and bounded. Hence at the point (x 1 , y 1 ) where the function δ n (x, y) reaches its positive maximum Div (σ P Grad δ n (x, y)) ≤ 0.
By substitution this inequality to the equation (76) we obtain the estimation Using the analogous estimation for the negative minimum of the function δ n (x, y) and the first condition (13) we obtain |δ n (x, y)| ≤ z 2 0 c 2 n 2 π 2 max x,y |Div (σ P Grad ρ n (x, y))| . (77) The expansion in Fourier series corresponds to calculation of the coefficients as the integrals over z and the derivatives of q n (x, y) over x, y can be obtained by differentiating of q(x, y, z) under integrals. Hence the boundedness of the Fourier coefficients and their derivatives over x, y result from (78) and (7) ∂ k q n (x, y) The function 1/σ ∥ are presumed bounded by (13) as well as their first and second derivatives over x, y, and the estimations (79) are proved. Therefore the functions ρ n (x, y) has the same property by virtue of (75) The obtained estimations of the Fourier coefficients give the estimations of the functions themselves The numerical factor is obtained taking into account that the sum of the numbers 1/n 2 from n = 1 till ∞ equals ζ (2) and it can be seen by a table of Riemann's ζ− function that ζ(2) < 1.65.
The boundedness of the first and the second derivatives over x, y of the functions ρ n (x, y) (80) and the estimation (77) cause that the functions δ n (x, y) are bounded and decrease as 1/n 4 when n increases: (82) Therefore the Fourier expansion of the function δ(x, y) absolutely converges and (83) The numerical factor is obtained taking into account that the sum of the numbers 1/n 4 from n = 1 till ∞ equals ζ(4) and ζ(4) < 1.083.
The derivative of the function δ(x, y, z) over z can be calculated by differentiating of its Fourier series of the shape (72) By (82) the coefficients δ n (x, y) decrease as 1/n 4 . Therefore the Fourier expansion of its derivative over z absolutely converges and the derivative is bounded (84) The numerical factor is obtained taking into account that the sum of the numbers 1/n 3 from n = 1 till ∞ equals ζ(3) < 1.25. By the same reasons the second derivative over z is also bounded, but we do not use this property.
Since the Fourier expansions converge the summation of the Fourier coefficients (74) corresponds to the summation of the functions themselves r (x, y, z) = ρ(x, y, z) + εδ(x, y, z). (85) Substitute the representation (85) of the function r(x, y, z) to the problem (69) with (71) instead of the last condition since they are equivalent to obtain the problem for the function ρ(x, y, z) and for the rest term δ(x, y, z): The problem (86) is the set of independent boundary value problems for the ordinary differential equation over z for each fixed pair x, y. The solution for such a problem exists, is unique and can be calculated by integration: where ρ 0 (x, y) is obtained from the last condition (86) as the average value over z of the first summand. By (7,13) we can independently of (81) derive from this formula the boundedness of the first and the second derivatives of the function ρ(x, y, z). In particular, Taking into account the representation (85) and proved boundedness of the functions ρ, δ and their derivatives over z (83, 84, 89) we obtain the boundedness of r and ∂r/∂z: where

The energy norms of the errors
Estimate the derivatives over x, y of the functions r and δ in the energy norm. The derivatives of the function ρ satisfy the estimation (81).
Multiply the equation (69) by r and integrate it over the domain Ω. Transform the integral in the left-hand side by using of the Gauss -Ostrogradskii theorem. In view of the periodicity conditions over x, y and the boundary conditions (69) at z = 0, z 0 the sum of the integrals over all boundaries equals zero. Therefore ∫σ q(x, y)) dΩ.
Use the Cauchy -Bunyakovskii inequality for the right-hand side and omit the nonnegative summand in the left-hand side By the definition ofq(x, y) (56) the last integral is not more than the integral of q 2 .
Using the estimation r (90), boundedness of the functions q (7) and 1/σ P (10) we obtain where c 6 = z 2 0 c 3 c 5 ab/c 1 . To estimate the derivative of the function δ over x, y we multiply the equation (87) by δ and integrate it over the domain Ω. Transform the integrals in the left-hand side and in the right-hand side by using of the Gauss -Ostrogradskii theorem. We omit the integrals over the boundaries since their sum equals zero because of the periodicity conditions over x, y and the boundary conditions for δ at z = 0, z 0 (87). ∫σ Omit the first integral in the left-hand side since it is nonnegative and use the Cauchy -Bunyakovskii inequality for the integrals in the right-hand side Cancel the common nonzero factor ∫ σ P (( ∂δ ∂x Using (81) and (10) we obtain the result estimation where

original problem
The solution for the original problem is represented in terms of f, p by the formulae (37). In view of the presentation (49) it can be represented over the solution of two-dimensional problem F, P and three-dimensional correction r, u. The vector function u is involved only through the vector curl u, that is represented in terms of the derivatives of the function r (67). It permits to simplify the expressions (37) for the solution of the original three-dimensional problem It can be seen by these formulae that the correction to the electric field that is obtained as the solution of the two-dimensional problem (57) is potential. Just the function εr(x, y, z) is the correction to the electric potential and its components ρ(x, y, z) and δ(x, y, z) (85) are estimated in (83, 89).
Hence the error of the electric potential is of order O(ε) if the solution for two-dimensional problem is used as approximation, as it is usually done in practice. The error decreases to O(ε 2 ) if we take the three-dimensional correction into account but approximately take ρ(x, y, z) instead of r(x, y, z).
Since the derivative of the function δ(x, y, z) over z is bounded (84) the errors of e z are j z of order O(ε 2 ) and O(ε) if we approximately use ρ(x, y, z) instead of r(x, y, z). If we use only the two-dimensional solution, i.e., put r = 0 in the formulae (93), then the error is of order O(ε) for e z and O(1) for j z . Both functions e z and j z equal zero in this approximation. Such a simplification is usually done in practice with some remarks about nonzero but small value of e z and finite j z .
The errors of the horizontal components of e, j can be estimated using (91, 92) in the norm of the space L 2 (Ω). As it was mentioned when the two-dimensional problem (57) was analyzed the correspondent parts of e x , e y , j x , j y are summable squared. If we neglect only the correction δ(x, y, z), then the error of e x , e y , j x , j y in the norm of the space L 2 (Ω) is O(ε) in accordance with (92). If we neglect the whole function r(x, y, z) then this error increase to O( √ ε) as it can be seen by (93, 91). The two-dimensional model corresponds to equipotentiality of each magnetic field line and to electric currents between these lines each of which is regarded as an unified object. The correction εr, εu provides currents along a magnetic field line that is necessary because of that the given source of the electric current q depends on z with no correspondence to the two-dimensional model. Naturally, the density of vertical current is not small if q is not specially defined.

Usage of the 2-D model for simulation of the electric field penetration from ground to the Earth's Ionosphere.
Measurements on the Earth's surface show quasi stationary electric field disturbances in the epicentral area before earthquakes. The amplitudes of those disturbances are in the range from several V/m before weak events up to 1 kV/m before strong earthquakes [21]. It is important to know what kind of effect can potentially be observed from space in connection with ionospheric precursors. Variations of the ionospheric electric fields related with earthquakes are reported in many papers [10]. The electric current penetration through the atmospheric conductor is proposed in the model [13] as the main physical mechanism of the lithosphere influence to the ionosphere. Such a model means the equations (1) in the atmosphere and ionosphere and given electric field near ground with boundary condition above the ionosphere that excludes influence of the magnetosphere.
It is shown in our paper [6] that this effect is negligible in contrast with models like [13]. We used the same model of electroconductivity (1) as [13] in the atmosphere, but we also took into account the ionospherical conductor. The 2-D model (59), that is conventional for large scale ionospheric electric field simulation, was used while it is not mathematically proved. The geophysical aspects of such a simulation are presented in our recent paper [5]. Here we use a simplified version of the model to apply the reduction of the boundary value problem presented in previous sections of this paper.

1-D boundary value problems
Here we would like to demonstrate adequacy of the 2-D model to the 3-D model in a very simple case. So we additionally suppose independence of the solution of the horizontal coordinate y and harmonic dependence of x and zero right-hand sides in the equations (1). If magnetic field is vertical and conductivity depends only on height Hall conductivity can be excluded from the equations. Then the problem can be reduced to the problem for the electric potential V (14) that has shape V (x, y, z) = f (z) sin (kx).
The equations (1) means ordinary differential equation for the function f (z) if we suppose the same presentation for the right-hand side q(x, y, z) = q(z) sin (kx). The properties of the boundary value problem differ much if k = 0, but here we are interested in nonzero k only.
where j 0 = σ ∥ (0)E 0 It can be transformed to the condition for the function Currents from the ionosphere to above could be taken into account as the right-hand side in the upper boundary condition that for the function f (z) means The 3-D model is reduced to the 2-D one if we put ε = 0 above some height z up that means σ ∥ = ∞. In accordance with (59) it means the upper boundary condition for the problem below this height where the height integrated horizontal current density J satisfies Ohm's law with the integrated Pedersen and Hall conductivities We did not define these parameters in the section 5, since for a homogeneous layer the integrals (99) differ from local values only by multiplier z ∞ − z up .
The boundary condition (98) is the charge conservation law for the ionosphere: current from the atmosphere j z (z up ) is closed by J in the ionosphere. So the righthand side Q in (98) includes all current sources above z up that we denote as Q 0 sin (kx).
For the function f it takes shape: Since q(z) = 0, j m = 0 the last boundary condition is also homogeneous Q 0 = 0. The numerical method is described in the next sections.

The energy functional
We consider smooth functions f (z) at the interval 0 < z < z up that is the domain Ω and define a norm of W (1) 2 (Ω) space as that is possible in view of nonzero k. Consider the following symmetrical bilinear form that is the energy scalar product. Let us use the Newton-Leibnitz formula for the function f (z): To square the identity and use the Cauchy -Bunyakovskii inequality By integration over the domain Ω we obtain Hence In couple with restrictions (10,11) and universal physical property σ ∥ ≥ σ P that also can be seen in Fig. 1 the inequality (101) proves the equivalence of the energy norm and the norm of the space W (1) 2 (Ω). The energy functional has shape The value f (0) satisfies the same inequality as (101). The first linear functional can be estimated by the Cauchy -Bunyakovskii inequality and so it can be estimated from above by the energy norm with factor L 2 (Ω) norm of the given function q(z) which is independent of f . Hence all linear functionals in (102) are bounded ones. Since the first term is positive definite and the linear one is bounded, the energy functional W (f ) has a unique minimum in the energy space. It is simple to testify that the equations (94, 95, 98) are the conditions of this minimum.

The finite element method
We additionally suppose that the weak solution for the problem (94, 95, 98) f (z) belongs to W (2) 2 (Ω). Then piece-wise linear functions f h can be used to approximate it with error where C is some constant [17]. We use grid with constant step h = z up /N . The grid nodes are numbered by integer 0 < n < N and segments have semi-integer numbers. Denote Conductivity σ ∥ is constant at each grid segment. We simplify approximation of the term σ P k 2 f 2 . Then the the energy functional is a quadratic form of the grid values of f h : Minimum of such a quadratic form means satisfaction to the system of linear algebraic equations The matrix of this system is a symmetrical and positive definite one because of these features of the operator since a piece-wise function f h (z) belongs to the energy space and its norm in the space L 2 (Ω) is equivalent to the Euclidean norm of its grid values.
A simple way to solve this chain of the equations is to construct two auxiliary solutionsf n andf n , which corresponds to Cauchy problems withf N = 1 start value and zero right-hand sides andf N = 0 with original righthand sides. The values in the grid node N −1 can be obtained from the last equation (103) and all other values consequently from the equations number n = N −1, ..., 1. The first equation that corresponds to the node n = 0 is not used. The solution for the system (103) can be constructed asf n + αf n if the value of α is chosen to satisfy the first equation: The coefficient in front of α does not equal zero. Otherwise the functionf n would satisfy the first equation with zero right-hand side and so the system with zero right-hand sides would have nonzero solution that contradicts to positive definiteness of the matrix.
The described finite element method can be used for the original boundary value problem (94, 95, 97) since it is a particular case of (94, 95, 98) when z up = z ∞ and therefore Σ P = 0 in view of definition (99).

Eigenvalues and condition numbers
Condition number of the matrix (103) is its important feature because it defines errors of the solutions and convergence of iterational methods. Since the matrix is a symmetric and positive definite one all its eigenvalues are positive and its condition number equals to the ratio of maximum and minimum eigenvalues.
The simplest way to find the maximum eigenvalue and corresponding eigenvector is to multiply an arbitrary vector with this matrix. The minimum eigenvalue and corresponding eigenvector can be obtained by multiplication an arbitrary vector with inverse matrix. In fact it is not necessary to find the inverse matrix itself. It is enough to solve a few times the system (103) with right-hand side that equals to the solution obtained at the previous step.
For clear demonstration of the properties we use comparison with classical results for 1-D Neumann problem for the Helmholz equation. It is a particular case of our problem (94, 95, 97) with σ P = σ ∥ = 1. Then the system of linear algebraic equations (103) takes shape The equations are written with multiplies 1/h or 2/h as it is usual in the theory of differential schemes [9]. This three-diagonal matrix has N + 1 eigenvectors m = 0, ..., N f m n = cos ) + k 2 .
The minimum and the maximum eigenvalues with eigenvectors f n = 1 and f n = (−1) n correspondingly.
We are interested in grids with hk << 1. So Λ ∆ ≃ (2/h) 2 and condition number We prefer to deal with which are almost independent of grid step h for our problems.
For the Helmholz equatioñ In general case the minimum and the maximum eigenvalues can be found numerically as it is described in previous section. We study the matrix of the equations (103) with the same multiplies 1/h or 2/h as in (105). When z up = z ∞ and therefore Σ P = 0 For the reduced problem with the main part of the ionosphere z > z up = 90 km presented with the boundary condition (100) we obtaiñ Comparison of the values (106-108) shows that condition number is 10 6 times improved by reduction of the problem while it stays 10 6 worse than it could be if the coefficients σ P = σ ∥ = 1.
Of cause, it is simple to solved 1-D problem without such a reduction, but demonstrated improvement of the condition number stays about the same when 3-D model of the ionospheric conductor is reduced to 2-D model. It principally simplifies numerical method [4].

Some geophysical results
Our quasi stationary three dimensional model of electric fields and currents in the conductor that includes the Earth's atmosphere and ionosphere is presented in [5]. Known approaches regarding the ionosphere as a boundary condition at the upper boundary of the atmospheric conductor are analyzed. It is shown that for the investigation of the electric field penetration from ground into the ionosphere it is sufficient to take into account only integral conductivity of the ionosphere. A mathematical simulation has demonstrated that the resulting electric field in the ionosphere is negligible in contrast to the general point of view that such a penetration is a physical process which potentially creates ionospheric precursors of earthquakes. Here we show how the 3-D model of the ionospheric conductor is reduced to 2-D model in [5].
The typical height distributions of the components of the conductivity tensor are shown in Fig. 1 for nighttime conditions during Solar minimum. These conditions provide maximum electric field penetration to the ionosphere. These distributions are obtained in [7] on the base of the ionospherical empirical models IRI, MSISE, IGRF and atmospheric model [20]. The dashed lines in Fig. 1 present the effective Pedersen and Hall conductivities, which describes the ionospheric conductivity averaged over an acceleration period of 3 hours by Ampere's force. A detailed explanation can be found in [7]. Just these distributions we use in calculations as normal night-time ones. The height z ∞ must be above the main conducting layer of the ionosphere. It has small influence on the result if we choose z ∞ = 500 km as we take in fact or 1000 km.
The typical values of vertical electric field perturbation near ground are E 0 = 100 V/m in the domains with horizontal size π/k = 300 km for moderate earthquakes [21,14]. Other right-hand sides equal zero in the boundary value problem (94, 95, 97) in this model.
We solve this boundary value problem numerically as it is described in section 10.3. The obtained height distributions of the components of the electric field are shown in Fig. 2. The horizontal component E x (0, z) and vertical one E z (0.5π/k, z) are plotted at this vertical lines since they have maximums there. It can be mentioned that any grid step h < 1 km does not vary the plots.
In accordance with Fig. 1 it is really small in the main part of the ionosphere: ε ≈ 0.15 at z = 80 km and decreases to 10 −8 with height. In the proofs above we used σ P , σ H , σ ∥ independence of height and here we analyze how the 2-D approach works for the real Earth's ionosphere. Obviously there is no sense to expand the solution as a power series of a small parameter ε in the atmosphere since ε = 1 below z = 50 km.
The upper boundary condition (96) can be modified to take into account adjoint ionosphere and conductor in the magnetospheric tail [1]. Taking these additional conductors into account would decrease the result electric field.
To have large difference from the 3-D solution we choose very low height z up = 70 km. As it can be seen in Fig. 2 this approximate solution differs from the 3-D one only near z up . They differ thrice at z up . If we take z up = 75 km the maximum error is decreased to 25%, z up = 80 km gives 3% error. If we take z up = 85 km when in the domain z up < z < z ∞ the small parameter ε = σ P /σ ∥ < 0.0004 then the error decreases to 0.3%. It ought be mentioned that not this error but accuracy of the horizontal component of the electric field E x in the ionosphere at z > 200 km where satellites can measure it is of value. Even for z up = 70 km such an error is not seen in Fig. 2, since it is about 0.04%. For z up = 80 km it is about 0.003%. Therefore the reduction to the 2-D model of the ionosphere does not vary the principal output parameter of the model. It ought be mentioned that if the event under analysis has smaller horizontal scale the error increases [5]. In the same paper [5] we have analyzed some inadequate simplifications of the ionospheric conductivity which are σ P = 0 above z up = 90 km in the model [13] or σ P = ∞ above z up = 150 km in the model [11].
The electric field E x ≈ 1.5 µV/m that is obtained in our model in fact can not be measured because electric fields in the range of several mV/m produced by other generators always are present in the Earth's ionosphere. So it is necessary to study other physical mechanisms to explain the lithospheric influence on the ionosphere. Some of them are presented in [16].
We expand the solution as a power series of a small parameter ε that is the ratio of transversal and field-aligned conductivities to obtain the two-dimensional problem (57, 52) that also has symmetrical positive definite operator.
The difference between the two-dimensional solution and the exact solution is estimated. The errors of the field-aligned components of electric field e z and current density j z are estimated in the norm of the space C(Ω) by the value O(ε) for e z and by the value O(1) for j z . The errors for the transversal components of the vectors e and j are O( √ ε) in the norm of the space L 2 (Ω). The two-dimensional problem and the problem for the three-dimensional correction are independent in the analyzed case, when conductor occupies a flat layer 0 < z < z 0 and is homogeneous in z direction, and the vector of magnetic field has only z component.
The simplest three-dimensional correction r can be easily obtained if we neglect the correction of the next order of ε. To do so it is enough to construct the function ρ by integrating the given function by the formula (88) and to take it as the approximation for the correction r. When this correction is included, the components of the approximate solution have errors of order ε 2 for e z , of order ε for j z , and the errors of all transversal components e x , e y , j x , j y are of order ε in the norm of the space L 2 (Ω).
In the same way, an approximate solution of the next order of ε can be obtained. For this purpose it is necessary to add terms to the representation (74). The chain of the formulae (75, 76) becomes longer. It is also necessary to require boundedness of higher derivatives of the coefficients and the right-hand sides, than it is done in (7,13).
The 2-D model is used for numerical simulation of the electric field penetration from ground to the Earth's ionosphere. Its accuracy is shown for real height distributions of the coefficients of conductivity tensor σ which are not constant in contrast with such a restriction in our mathematical proofs.
To use the obtained results for more realistic models of ionospheric electric fields [4,8] it is necessary to remove one more simplification that is used in this research. In the ionosphere the conductivity tensor σ has the shape (8) when the local z− axis has its own direction in each point of the domain. If one is not interested in exact justification, it is not a problem to use 2-D model in practice as it is described in [7].