Asymptotic shallow models arising in magnetohydrodynamics

In this paper, we derive a new shallow asymptotic model for the free boundary plasma-vacuum problem governed by the magnetohydrodynamic equation, vital in describing large-scale processes in flows of astrophysical plasma. More precisely, we present the magnetic analogue of the 2D Green-Naghdi equations for water waves in the presence of weakly shared vorticity and magnetic currents. The method is inspired by developed ideas for hydrodynamics flows in by Castro and Lannes (2014) to reduce the $(d+1)$-dimensional dynamics of the problem to a finite cascade of equations which can be closed at the precision of the model.


Introduction
Plasma is an ionized gas consisting of freely moving positively charged ions, electrons and neutrals. It is by far the most common phase of ordinary matter present in the Universe. We are in constant contact with the small amounts that are not, such as the oceans and seas, but electrically charged fluids are everywhere throughout the galaxy, [WD14]. Macroscopic plasma processes are usually described by the so called magnetohydrodynamic equations (MHD), first proposed by the physicist H. Alfvén [Alf42]. Here the behaviour of fluid particles is governed by Newton's second law under the effect of an electromagnetic force described by Maxwell's equations. Furthermore, the movements of present charged particles create an electric field that also affects this magnetic field.
Motivated by the problem of magnetic plasma confinement in laboratory research, the plasma-vacuum interface problem has attracted the interest of the mathematical community in the last decades [Tra12,Tra16 ,GW16,Gu19]. However, studying the full dynamics of the equations is too complex, mainly because the moving surface boundary is part of the solution. This led physicists, oceanographers and mathematicians to derive and replace the original equations by approximate asymptotic systems in specific physical regimes. Those systems are more amenable to numerical simulations and their properties are more transparent.
In the case of the water waves equation describing the motion of an inviscid and incompressible fluid delimited by a free-surface , the most prominent example is the non-linear shallow water equations, also known as Saint-Venant equations, cf. [Ovs76,ASL08)]. In the shallow water regime, when µ := H 2 L 2 1 (where L is the typical horizontal scale and H the typical depth), the non-linear shallow water system is derived from the free-surface Euler equation by averaging and neglecting the O(µ) order terms. As a counterpart of simplifying the model by dropping the O(µ) order terms, we miss completely the dispersive effects, vital for many applications. Keeping them in the equations and just neglecting the O(µ 2 ) order terms one obtains the so called Green-Naghdi equations [GN76] or Serre equations [Ser53]. We refer to [Mak86,ASL08)] for a rigorous derivation 1 and [LB09] for a more recent review. Moreover, the Green-Naghdi system is one of the most used model to perform numerical simulations of coastal flows [CBB06,LM14]. In the equations the flow is assumed to be irrotational or almost irrotational, which breaks down in the presence of rip currents or when underlying currents are present. The complication to describe waves in the presence of non-trivial vorticity is due to its (d + 1)-dimensional nature while for irrotational flows are only d-dimensional.
Besides different approaches to deal with this difficulty, a novel strategy was developed in [CL15], where the additional terms appearing in the momentum equation can be treated without appealing to the resolution of the (d + 1)-dimensional vorticity equation. The strategy followed in [CL15], inspired by the so called turbulence theory, hinges on deriving a cascade of equations which is actually finite at the precision of the model ((O(µ 2 )-order) without any artificial closure. The resulting equations are an extension of the classical Green-Naghdi equations, where the non-hydrostatic pressure terms are affected by the interactions between the horizontal and vertical components of the vorticity. A rigorous justification of the asymptotic expansions and models derived in [CL15] have been studied in a companion paper by same authors in [CL14]. In conclusion, the hydrodynamic shallow water models have been well-studied and extensively used to describe fluid motion in oceans, with direct applications in coastal engineering, [Lan13,Lan20]. Surprisingly enough, most of the electrically conducting fluids appearing in astrophysical plasmas, such as accretion disks, planetary atmospheres or stars is that they are, in some sense, thin. However, the role that magnetic field plays when the fluid is electrically conducting is far from being well understood.
The introduction of magnetic effects into the shallow water system was first proposed by Gilman [Gil00] and is used as a model of the solar tachocline. The tachocline, first coined by Spiegel and Zahn in [SZ92], is a very thin layer in the Sun of thickness of about two and five per-cent of the solar radius which bridged the transition region between the convective zone from the radiative zone, [HRW07]. The approximation of the shallow MHD equation neglects O(µ) order terms and hence can be understood as an magnetic analogue of the non-linear shallow water equations. Since their derivation they have been investigated from a theoretical and modelling point of view. Studies on linear and non-linear waves have been treated in [SBG01] as well as applications of the shear-flow instabilities in [MGH16 ]. The shallow MHD system has also been shown to be hyperbolic and enjoys a Hamiltonian structure [Del02,DeSt01,Ros0213]. Extensions from one-layer to multi-layers models of the shallow MHD equations, experimenting with different stratification settings and enriching the model by choosing different states on each layer (different velocities and different magnetic fields) where also documented in [Hun13,Zei13].
The main purpose of this article is to derive new asymptotic models, in dimension d = 1, 2, for the shallow MHD equation (SMHD) in the presence of vorticity and magnetic currents up to a precision of O(µ 3/2 ). Our strategy is influenced by the work in [CL15] for the water waves equations based on ideas reminiscent of the turbulence theory. The new models take into account the dispersive effects missed in the model by Gilman [Gil00] due to the non-hydrostatic pressure and more important manage to deal with the strongly coupling between the vorticity and the current equations, which are assumed to be weakly sheared plasmas (see §3.2). This is achieved, by getting rid of the (d + 1)-dimensional vorticity-current system, and rather look for an equation involving two d-dimensional quantities, i.e., the shear velocity and magnetic shear which are the key terms to close the cascade equation. Roughly speaking, those quantities represent the contributions of the horizontal vorticity and magnetic current into the horizontal momentum equation. We should mention that the asymptotic models derived in this paper are not rigorously justified which is out of the scope of the present article and is left for future work.
Plan of the paper. The paper is structured as follows. In Section §2 we present the basic equations of the free-surface MHD problem. Section §3 is devoted to the asymptotic analysis of the free-surface MHD equation, where we first recast the MHD problem using the elevation-discharge formulation. Subsection §3.2 deals with non-dimensional averaged-formulation and further remarks on the asymptotic regimes studied on this paper. An asymptotic description of the velocity field and magnetic field is performed in §3.2, which depends on the velocity and magnetic shear term. Section §4 presents the different shallow asymptotic models up to different orders of approximation. In subsection §4.1, we turn to the derivation of the 2D magnetic Green-Naghdi equations, where we first derive an evolution equation for the velocity and magnetic shear §4.1.1 and later a closure equation for the tensors in §4.1.2. Lastly, in §4.2 we also handle the 1D magnetic Green-Naghdi equation, where several terms trivialize. Finally a conclusion and future work perspectives are given in Section §5.
Notation. We will use the following notation throughout the manuscript.
• We denote by d = 1, 2 the horizontal dimension, by X ∈ R d the horizontal coordinate and by z the vertical variable.
• For every vector field F ∈ R 3 we denote by F h its horizontal component and by F v its vertical component. If the vector field F ∈ R 2 we denote by F ⊥ = (F 1 , F 2 ) ⊥ = (−F 2 , F 1 ). Let A ∈ R 2 × R 2 be a matrix, then A t is the transpose of A. Let A, B ∈ R 2 × R 2 two matrices, we denote by A : B the standard matrices product. For any two different vector fields F, G ∈ R 2 , we define the tensor product F ⊗ G = F G t , as the usual outer product.
• We use ∇ to denote the gradient with respect to the horizontal variables and ∇ X,z is the full threedimensional gradient operator. The rotational and divergence operators are given by

The basic equations
In this section we are concerned with a free-boundary problem of ideal incompressible magnetohydrodynamics. The problem consists in finding a variable domain Ω − t occupied by an electrically conducting homogeneous plasma, together with a velocity field U = U (X, z, t), the scalar pressure P = P (X, z, t) and the magnetic field B = B(X, z, t) satisfying the equations of MHD. The elevation of the free surface is parametrized by the graph of a function ξ(·, t), and the non-moving bottom topography is parametrized by a time independent function −H 0 + b(X) where H 0 represents the the depth of the plasma and b the possible variation of the bottom, see Figure 1 . Therefore, the domain occupied by the plasma at time t is The dynamics of the plasma region is governed by the ideal incompressible MHD equations (2.2) where the external forces due to gravity g = −ge z are also taken into account. Above ρ is the density assumed constant, and µ 0 the magnetic permeability constant. It is assumed that the plasma is surrounded by vacuum region, Since, vacuum has no density, velocity or electric current, the pre-Maxwell dynamics apply [GP04,BFKK58)]. In such a case, the magnetic field B is determined by the div-curl ∇ X,z × B = 0, div B = 0, in Ω + t .
(2.4) Figure 1. The plasma-vacuum interface problem setting The plasmavacuum interface is now free to move since the plasma is surrounded by vacuum. The physical quantities of the plasma and vacuum region must satisfied non-trivial jump conditions connecting the fields across the interface (cf. [GP04] for a thorough exposition). The first boundary condition, the so called kinematic boundary condition, is related to the fact that the plasma particles on the surface stay on the surface moving with normal V n = U · N given by where U (X, t) |s = U (X, t, ξ(X, t)) and the N = (−∇ξ, 1) T . The second and third conditions satisfied at the surface, is the pressure balance condition and magnetic jump continuity where [[f ]] = f − f denotes the jump in a quantity across the surface. Hence, Finally, we also impose two boundary conditions at the bottom topography assumed to be a perfect conducting and impermeable material, this is In order to simplify the dynamics of the equations we will assume that the vacuum magnetic field B is identically zero, which of course is a trivial solution of (2.4). In this particular case, the free boundary MHD problem is given by (2.9) with boundary conditions Remark 2.1. In (2.10b) we are neglecting the effects of the surface tensions, which can be also incorporated into the model, cf. [CD19]. Moreover, we also neglect the Coriolis effects induced by planetary rotation.

The averaged MHD equations and asymptotic analysis
In this section we will recast the MHD equations (2.9) using the so called elevation-discharge formulation (cf. [Lan13] )that proves to be very convenient in order to obtain and understand different asymptotic models. The central idea of the formulation is to get rid of the vertical variable by integrating vertically the horizontal component of the free surface MHD equation. To that purpose, let us denote by U = U (X, z, t) = (U h (X, z, t), U v (X, z, t)) , B = B(X, z, t) = (B h (X, z, t), B v (X, z, t)) the horizontal and vertical component of the plasma velocity and magnetic field, respectively. Moreover, we introduce the horizontal velocity discharge Q, (3.1) and the horizontal magnetic discharge Q B , Integrating vertically the horizontal component of equations (2.9) and using the boundary conditions (2.10a), (2.10c),(2.10d) and (2.10e) gives where P B = P + 1 2µ0 |B| 2 is the magnetic pressure. Next, we can decompose the pressure magnetic field into a hydrostatic magnetic pressure term P h B and non-hydrostatic magnetic pressure term P nh B . It is easy to check that ξ = 0, U = 0, B = 0 is a particular steady state solution of equations (2.9) and therefore the vertical component of the first equation in (2.9) gives the following ordinary differential equation for P B : with boundary condition P B | z=0 = 0, and hence the solution is the hydrostatic magnetic pressure P h B = −ρgz. Similarly when the fluid is not at rest, integrating the vertical component between z and ξ we have that where the non-hydrostatic magnetic pressure is given by Plugging (3.5), we have that the evolution equations are given by In addition, we also decompose the horizontal velocity and magnetic vector field as where for a general function g(·, t) as where h(X, t) = H 0 − b(X) + ξ(X, t) and g (X, t, z) = g(X, t, z) − g(X, t). (3.9) Using (3.7a)-(3.7b) into the equations, we have that We will refer to equations (3.10) as the averaged MHD equations. Although the equations are exact, they are not closed. Indeed, quick inspection of equations (3.10) reveal that the terms R, R b , R m and the nonhydrostatic magnetic pressure term are non explicit in terms of the functions ξ, Q, Q B . Actually, as we will see in §4.1, they will depend on the vorticity ω = ∇ X,z × U and magnetic current j = 1 µ0 ∇ X,z × B equations given by (3.11) 3.1. Non-dimensional free boundary MHD equation. To study the behaviour of the solutions of the full system, is generally too complicated, since it contains all the information of the dynamics. We will therefore simplify the terms which seem less important to us. To determine which terms are more relevant, we will use very well known method based on the principle of dimensioning. We non-dimensionality the equations by using several characteristic length of the problem, namely: the typical amplitude of the waves a s , the typical depth H 0 , the typical horizontal scale L and the order of bottom variations a b . With these characteristic scales we can construct three different dimensionless parameters: where is called the amplitude parameter, β the topography parameter and µ the shallowness parameter.
Remark 3.1. For simplicity we have assumed here that the horizontal scale L is the same in both the longitudinal L x and transversal L y direction, however this could also be considered into the modelling yielding a new parameter γ = Lx Ly called transversality parameter [Lan13]. Since the main goal of this article is to derive shallow models we will assume that µ 1, but not other smallness assumption will be done.
With the before-mentioned parameters, we will define the following dimensionless variables and dimensionless unknowns functions (3.14) Remark 3.2. The nondimensionalization of the time variable, velocity field and magnetic field are based on a linear analysis, similar to the one performed for the water waves equation (cf. [Lan13,CL14,CL15]).
With this new variables at hand, we denote by the free boundary MHD equations (2.9) are given in dimensionless form (omitting tildes) by posed in the dimensionless fluid domain, In the same way, the boundary conditions in dimensionless form are given by In a similar way, the dimensionless version of the averaged MHD equations (3.10) are given by where the dimensionless tensors and non-hydrostatic term is given by Analogously as in the case the free-surface equation Euler equation in the presence of non-trivial vorticity [CL14, CL15], we require a true comprehension of the behaviour of the vorticity ω and the current j. To that purpose, let us do the following digression: due to the incompressibility condition, and the boundary value condition, . Moreover, by a straightforward computation we infer that hence we might rescale them as Hµ 0 gρ j.
The horizontal component, denoted by ω h , j h respectively, is given by . Since we will treat only the case of weakly sheared flows [Tes07, RG12, RG13, CL14, CL15], this is we will assume that ω h , j h of order O(1). To that purpose, we dimensionality the vorticity ω and current j as follows where we have dropped the tidal notation.
Remark 3.3. Deriving different models based on different assumptions on the strength of the vorticity and magnetic current is also possible. However, here we treat the case of weakly shared plasmas, ω µ , j µ are a O(1) quantity with respect to µ. In [GG12, Lan20] models with stronger rotational effects are studied for the water waves problem, however their rigorous well-posedness justification is still an open problem. We recall here that the aim of the present paper is only to derive new asymptotic models for the (SWMHD) equation, but no a rigorous proof.
3.2. Asymptotic expansion of the inner velocity and magnetic field. In this section, we will derive inner asymptotic description of the velocity field U µ and magnetic field B µ . First, let us consider the the following boundary value problem satisfied by the velocity field: The expansion coincides with the one in [CL15] for the water waves equation (where B µ = 0) and is given by where the shear velocity U sh is given by and the operator L acting on a function f is given by Next, we will consider the expansion of the magnetic field B µ , which satisfies: Notice that the boundary value problem for the magnetic field B µ , is completely determined by the magnetic current j µ . In contrast with the velocity field where the irrotational part of the U µ is determined from the tangential component at the surface, the magnetic field does not have the extra degree of freedom (cf. Remark 3 in [CL15]). Taking the horizontal and vertical component of equations we have that We plug the following Ansatz In order to close the asymptotic models up to precision of order O(µ 3/2 ) (cf. §4.2) we need to do the following mild assumption on the vertical component of the current; assuming that j 0 v (x, y) = 0. This can be imposed on the initial data, and check that it is conserved during the evolution of the current equation (3.11), j 0 v (x, y) |t=0 = 0 implies that j 0 v (x, y, t) = 0 ∀t > 0.
(3.30) By the second equation in (3.28), we can write for some potential stream function φ satisfying the the elliptic system with h is smooth bounded function. Therefore, invoking standard elliptic estimates we deduce that B 0 h = 0, [GT83]. 2 In a similar way, we have that for the O( √ µ)-order (3.32) Using the third equation we have that We are implicitly imposing that we want to deal with physically meaningful functions, i.e. functions with finite energy φ ∈ L 2 (R).
for some potential function φ 1 (x, y). Hence, where the magnetic shear B sh is defined as The unknown potential function φ 1 (x, y) satisfies that Equations involving the O(µ)-order terms are given by (3.35) h (x, y, ξ)) = 0. (3.36) By the last equation in (3.36) we have that and hence by the first equation in (3.36) which implies by the argument as in (3.31) that B 2 h = 0. Finally, the system for the O(µ 3/2 )-order satisfies (3.38) The vertical component B 3 v is then given by Therefore, using (3.33) where the operator L is defined as in (3.26). By the third equation in (3.38) φ 3 (x, y) satisfies Therefore, we have that the magnetic field B µ can be describe by the following asymptotic expansion (3.43)

The 2D magnetic Green-Naghdi equation
We deal here with the derivation of the magnetic Green-Naghdi type 3 of equation in the two dimensional case (d = 2), however we will only treat with models up to precision O(µ 3/2 ). Structural complications arise in order to close the two dimensioanl cascade of equations if we want to push the expansion further to full order O(µ 3/2 ). Those difficulties are due to the strong coupling between the magnetic field and velocity field, as we will notice in the computations.
Let us first present the, nonlinear shallow MHD equations, which are an approximation of order O(µ) of the full MHD equations (3.18), in the sense that we drop all the terms of order O(µ). At this order, we do not need to assume that the vertical component j 0 v is zero as we did when deriving the asymptotic expansions in §3.2. Actually, the asymptotic expansion would be given by 1 (x, y)). Plugging the asymptotic expansions into the non-hydrostatic magnetic pressure (3.19) and the tensors (3.20)-(3.21) we have that The non-hydrostatic and turbulent effects coming from the tensors are not present in this first approximation. Therefore, neglecting terms of order O(µ) in (3.18), one obtains the following system of equations: The system (4.3) can be understood as the magnetic MHD version of the well-known nonlinear shallow water equations, [Ovs76].
Remark 4.1. The non-linear shallow MHD equations presented above where first derived by [Gil00], and studied since their derivation intensively, see [Hun13,Mak13] for recent reviews. In [Del02], the authors studied the hyperbolic character of the non-linear shallow MHD equations. General theory for symmetric hyperbolic systems (cf. [BS07)]) assures the local well-posedness for initial data (ξ(x, 0), In recent paper, Trakhinin [Tra19 ] investigated the structural stability of shock waves and current-vortex sheets in the non-linear shallow MHD equation (4.3). 3 We have coined the equation 2D magnetic Green-Naghdi equation due to the analogy with the Green-Naghdi equations for water waves. However, since we only deal with precision up to O(µ 3/2 ), this could be understood as a Green-Naghdi equation with medium amplitude wave assumption = √ µ.

4.1.
The magnetic 2D Green-Naghdi equation. Next, we will compute the contributions of the nonhydrostatic terms (3.19) and the tensors (3.20)-(3.21) up to a precision of O(µ 3/2 ). Using the asymptotic expansion (3.23) and (3.42) into (3.20) and (3.21) yields Similarly, we can compute the non-hydrostatic magnetic pressure contributions. Dropping O(µ 3/2 ) terms we get where the operators F and D , are defined as and We have derived the non-hydrostatic magnetic pressure contributions using the operators F and D which we borrowed from the formulation in [BCLMT11]. Therefore, from inserting the expressions (4.4)-(4.6) and the pressure (4.7) in the equations (3.18) we infer that (4.8) where Noticing that ∂ t h + ∇ · (hU h ) = 0 and denoting by B h = √ µ 1 h ∇ ⊥ φ we can rewrite equation (4.8) in a more compact way, namely (4.9) Hence, to close the equations we have find closure equations for the quantities R, R b , R m . To that purpose we have first to derive an evolution equation U sh in (3.25) and magnetic shear B sh defined in (3.34). With this equations at hand, we will be able to find the closure equations and close the system. 4.1.1. Evolution equation for U sh and B sh . The dimensionless the vorticity and the current are given by Hence the dimensionless vorticity-current equations of (3.11) (4.10) Let us first derive an evolution equation for the shear velocity U sh . The computations are an adaptation of §2.3 in [CL15] .We recall the main steps, highlighting the new modifications. The horizontal component of the vorticity equation (4.10) is given by Using the expansion (3.23)-(3.24) and (3.42)-(3.43), we have that (4.12) Using (4.13) Taking the perpendicular operator ⊥, integrating the equation between z and ξ and using boundary condition ∂ t ξ + ∇ · Q = 0 at the surface we have that (4.14) Applying operator 1 h ξ −1+βb to (4.14) and subtracting the resulting equation from (4.14), we infer that Remark 4.2. Due to the assumption j v 0 ≡ 0 in the asymptotic description of the magnetic field described in Subsection 3.2 we notice that equation (4.15) coincides with the equation (2.32) in [CL15] for water waves. Without the j v 0 ≡ 0 assumption, we could also derive an evolution equation for U sh involving non-trivial contributions of the magnetic field, however we would not be able to close the system.
In a similar way, the horizontal component of the current density equation (4.10), is given by Plugging in the asymptotic expansion (3.23)-(3.24) and (3.42)-(3.43) and dropping terms of order O( √ µ), (4.16) Taking the perpendicular operator ⊥, integrating in z and using boundary condition B µ | z=ξ · N µ = 0 we have that Noticing that ∇ ⊥ U h ⊥ : B sh = (B sh · ∇ ⊥ )U h ⊥ and using the vectorial identity Hence, the equation can be rewritten as Taking the average to (4.18) and subtracting it from (4.18), Therefore, taking the time derivative on the tensors R, we have that Using the fact that ∂ t ξ + ∇ · (hU h ) = 0, we obtain that Mimicking the computations but for the magnetic shear equation for B sh given by we arrive to Finally we compute the closure equation for R m , which differentiating the tensor R m gives and To last term is given by Collecting all the computations, and recalling that ∂ t ξ + ∇ · (hU h ) = 0, we have that 4.1.3. Full 2D magnetic Green-Naghdi equation. We can now express the two-dimensional magnetic Green-Naghdi equation, dropping O(µ 3/2 ) terms by where the operators F, D are defined as 2h (∇(h 1 ∇b · g) − h 2 ∇bg · g) + ∇b∇b · g, D(g) = −2R 1 (∂ x g · ∂ y g ⊥ + (∇ · g) 2 ) + R 2 (g · (g · ∇)∇b), and while the tensors R, R b , R m and R S b stand for Remark 4.3. If we compare the equations with the one derived for the rotational water-waves in [CL15], we observe new two new phenomena. First we have a self interaction of the magnetic shear B sh coded in R b and the interactions between the shear velocity and magnetic shear gathered in the term R m and R t m . 4.2. The 1D magnetic Green-Naghdi equation. For the sake of completeness we also derive the equations in the one-dimensional case. We consider velocity and magnetic fields given by Therefore the scalar vorticity and magnetic current are given by Moreover we will have that U sh = (u sh , 0) t and B sh = (b sh , 0) t , and hence Therefore, straightforward modifications of the two-dimensional case, yield the following one-dimensional magnetic Green-Naghdi equations with the one-dimensional versions of F, D are given by x b∂ x b while we recall the one dimensional tensors R and R b are given by Remark 4.4. We notice that the main difference between the one-dimensional setting and the two-dimension version of the magnetic Green-Naghdi equations is that the one dimensional version of the averaged MHD equations (3.10) trivialize to (4.22) Hence the magnetic elevation discharge Q B = 0, if we want to Q B has finite energy, therefore excluding the interaction of the velocity and magnetic shear (present in (4.20) as R m ) and the coupling between the magnetic field B h in the momentum equation, also present in (4.20).

Conclusion and future work
In this paper, we have derived a new shallow water models for the free-surface MHD equation in the presence of vorticity and magnetic currents. The most essential ingredient of the model is that the resulting equations are d-dimensional which avoid the need to solve the (d + 1)-dimensional nature of the vorticitycurrent equations, reducing the complexity of the system from a mathematical and numerical point of view. The strategy follows closely the ideas developed in [CL15]. It is shown that the additional terms appearing in the momentum equation due to the presence of the vorticity and current effects, satisfy additional two dimensional advection type equations which couple to the system. The advected quantities describe the self-interactions of the shear velocity and magnetic shear induced by the vorticity-current system, and the coupled interactions between the shear velocity and magnetic shear.
Different techniques were developed to perform numerical simulations for the non-linear shallow MHD equation first derived in [Gil00], for example, the constrained transport approach [DeSt02,Ros03], projection method [Ros0213] or central upwind scheme [KT00,ZAQ14]. Therefore a natural perspective is to take into account the vorticity-current effects in the numerical simulations, allowing the modelling of underlying currents in plasmas or sheared plasmas, using the models derived in this article. Another future research direction, is related to the rigorous justification or the well-posedness of the derived models, which to the best of the authors knowledge is an open problem even for the non-linear shallow MHD equations derived by [Gil00].