A Method for Finding Exact Solutions to the 2D and 3D Euler–Boussinesq Equations in Lagrangian Coordinates

We study the Boussinesq approximation for the incompressible Euler equations using Lagrangian description. The conditions for the Lagrangian fluid map are derived in this setting, and a general method is presented to find exact fluid flows in both the two-dimensional and the three-dimensional case. There is a vast amount of solutions obtainable with this method and we can only showcase a handful of interesting examples here, including a Gerstner type solution to the two-dimensional Euler–Boussinesq equations. In two earlier papers we used the same method to find exact Lagrangian solutions to the homogeneous Euler equations, and this paper serves as an example of how these same ideas can be extended to provide solutions also to related, more involved models.


Introduction
We continue our investigations of finding explicit solutions to the incompressible Euler equations in the Lagrangian framework.In [12,13] we considered the homogeneous case while here we treat the Boussinesq approximation of the heterogeneous case.The idea behind the Boussinesq approximation is that the density of the fluid does not fluctuate very much around the mean value, so that one supposes that the density is constant, except in the term that represents the gravity.There are many physical situations where this simplification is justified.The general introduction to this type of models can be found in [11], while in regard to the Lagrangian description in general we refer to [3].
One of the first explicit solutions to the homogeneous Euler equations in Lagrangian formulation is due to Gerstner [9], and it turns out that solutions of the Gerstner type show up in many related models.For example the Gerstner solution plays a prominent role also in the heterogeneous Euler equations, even without using the Boussinesq approximation.In [8] it was shown that the Gerstner wave can also be a barotropic flow, i.e. a flow where pressure and density are functions of each other, see also a modern exposition by Stuhlmeier [14].Gerstner type solutions are also relevant in many physical models that are somehow modifications of the Euler equations.For example Constantin [4] found a barotropic Gerstner type exact solution to the equatorial water wave equations with the beta-plane approximation of the Coriolis effect.For more papers on the applicability of the Gerstner waves see for example [10,7,6,5] and especially the survey articles [2,15].In [16] there is also a Gerstner type solution to the first and second terms of the asymptotic expansion to the heterogeneous Euler equations.One remarkable property of Gerstner type solutions that makes them popular, is that they can in many cases be interpreted as free surface solutions, or in other words one can model the interface of two different fluids.Apparently there are no other known explicit solutions with this property.
Many explicit solutions to stratified fluid flows are shear flows, which are flows with two-dimensional horizontal motion that only depends on height.Most of the studies use the Eulerian framework.As for the Lagrangian description, Yakubovich and Shrira [17] found solutions with columnar motion using the Boussinesq approximation, but it seems that elsewhere the Eulerian framework is used.
As stated above there are numerous models for geophysical fluid flows, depending on the particular context.We have chosen to analyze the Euler-Boussinesq model in the present article.For a thorough discussion of this model and its applicability we refer to [11].As in [12,13], we use a separation of variables method to compute our solutions, that is, we search for fluid particle maps that can be expressed as the product of a time-only dependent matrix and a spatial-only dependent vector.This approach leads for example to a two-dimensional Gerstner type flow (4.2), as well as a plethora of other solutions both in the 2D and the 3D setting.Solution (4.2) does not satisfy the free boundary condition like most Gerstner type flows, but it still gives valid internal flows, which are an important application of the Euler-Boussinesq equations.Note in particular that our approach is not restricted to the Euler equations and Euler-Boussinesq equations, and we are confident that it could be used profitably to analyze other related models.For example Abrashkin [1] derived the governing Lagrangian equations for equatorial beta-plane flows, allowing anyone to readily attempt the separation of variables strategy that we use below.Also the Gerstner type solutions in all those different models discussed above could have been found by our method.This paper is organized as follows.In Section 2 we give basic definitions and specify the model which will be analyzed.In Section 3 we discuss columnar flows and generalize one family of solutions which was given in [17].In Section 4 we find exact solutions to the two-dimensional Euler-Boussinesq equations while utilizing the similarities to the homogeneous situation [12].In Section 5 we consider the three-dimensional case and use the framework outlined in [13] to find solutions to the Euler-Boussinesq equations.
In some cases we are able to find explicit periodic and nonperiodic solutions.In other cases one can show that a solution exists for all times for convenient parameter values.Finally there are cases where all one can say that a local solution is well-defined.We allow for both stably and unstably stratified fluids as both types of situations can often be covered by a single formula, though when we give specific examples we concentrate on the stably stratified case.Some solutions that we find are neither stably nor unstably stratified.We also note that it has not been possible to explore all different cases in the present article.

Notation
Let A : R → R n×m and v : R n → R m .We denote the columns of A by A i and its entries by a ij .Now the minors of A are denoted by Similarly the the minors of dv are 2) It will also be useful to define where A ′ i is the (time) derivative of A i and G ij is useful only when n = 3.For derivatives of v we use multiindices so that for example At times we will also say that functions v1 and v 2 are an anti-Cauchy-Riemann pair, or an anti-CR pair, if

Model
The n-dimensional heterogeneous incompressible Euler equations are given by the system Here u is the velocity field, ρ is the density, g is the acceleration due to gravity, and e n is the vertical unit vector parallel to gravity.Let us write these equations in the Lagrangian framework.Let D ⊂ R n be a domain and let us consider a family of diffeomorphisms φ t : D → Ω t = φ t (D).The coordinates in D are denoted by z and in Ω t by x.We can also define Now given such φ we can define the associated vector field u by the formula Then (u, ρ, p) solves (2.4) if det(dφ) ̸ = 0 and where p = p • φ and ρ = ρ • φ.Hence ρ is a function of spatial variables only, which is one great advantage of using the Lagrangian description. 1 Typically one cannot explicitly recover u from φ, since it requires computing the inverse of φ t .One exception to this is obtained when φ = A(t)z for some square matrix A, yielding u = A ′ A −1 x.
The standard way to apply the Boussinesq approximation is to assume that ρ is constant in every term except when it is multiplied by g.Thus the Boussinesq approximation takes into account how density variations affect buoyancy.This would mean that Newton's second law (2.6b) would be replaced by the equation where ρ is the average density.However, as was shown in [17], we can make the model slightly more accurate without making the equations any more difficult to study.Supposing that n = 3 and taking the curl of (2.6b), we obtain Now assuming only ∇ρ × (dφ T φ ′′ ) = 0 in (2.8), dividing by ρ and defining ρ = g ln ρ we obtain If we had used (2.7), we would have obtained the same equation but with ρ defined as ρ = g ρ/ρ instead of ρ = g ln ρ.Hence it should be kept in mind that ρ is not really density, but the density ρ can be recovered from ρ using either ρ = g ρ/ρ or ρ = g ln ρ.
Integrating the left hand side of (2.9) with respect to t we obtain the equivalents of what are the Cauchy invariants for the Euler equations: This integrated form is often useful since it removes the second-order dependence of φ 1 and φ 2 .But calling the components of h the "Cauchy invariants of the Euler-Boussinesq equations" is perhaps a bit of a stretch since the time integral of φ 3 produces an arbitrary function of z and so there is no canonical way to choose h.
In the two-dimensional case h is just a scalar and it is convenient to write it in the following form:

Separation of variables
As in our previous papers [12,13], we try to find solutions of the form φ(z, t) = A(t)v(z), where A : R → R n×m and v : D ⊂ R n → R m .To this end we need convenient formulas for det(dφ) and h.Using (2.1), (2.2) and the Cauchy-Binet formula we obtain (2.12) Then using (2.3) we compute that where y i is a function such that y ′ i = a ni .Some properties of φ and ρ are gathered in the following Lemma.Lemma 2.2.Let (φ, ρ) be a solution to the Euler-Boussinesq equations.

Let ψ : D → D be an arbitrary diffeomorphism and let
If H is a regular m × m matrix with constant entries and ṽ = Hv, Ã = AH −1 , then ( Ãṽ, ρ) is a solution.
3. In the 3D case, if R is a constant rotation matrix such that (Rφ) 3 = φ 3 , then (Rφ, ρ) is also a solution.In the 2D case there is no nontrivial rotation R for which Rφ is always a solution.
The change of coordinates in the first part of this Lemma typically allows us to assume without loss of generality that (v 1 , v 2 ) = (z 1 , z 2 ) in the 2D case and (v 1 , v 2 , v 3 ) = (z 1 , z 2 , z 3 ) in the 3D case.In the present paper we always use this simplification, but in practice the general form allows for more flexibility since the inverse map needed for this transformation cannot always be explicitly computed.The second part of the Lemma can be used to bring the problems to simpler form without loss of generality.From this part it also follows that if A i or ∇v i are linearly dependent over R, the solution reduces to a case of lower m.On the other hand, ∇ρ being an R-linear combination of ∇v i does not imply that the solution reduces in this way.
The formulas (2.12) and (2.13) allow us to deduce what kind of constraints we should set for the spatial functions v i and ρ.We want the time derivatives of det(dφ) and h to vanish for all t, so, for example in the two-dimensional case, fixing any t produces constraints of the form 1≤i<j≤m where β ij , γ ij and γ i are constants.Thus the spatial functions need to satisfy a number of constraints like these.By substituting these constraint equations back to the formulas of det(dφ) and h we also immediately obtain the conditions that the time functions a ij are required to satisfy.In [12] and [13] we derived the spatial constraints from the formulas of det(dφ) and h for the homogeneous Euler equations, and especially in [12] we went into greater detail in the 2D case to show what were all possible constraint sets that are essentially different for the cases that we studied.In the case of the Euler-Boussinesq equations this analysis, which is done using the second part of Lemma 2.2, yields very similar results, so in the present paper we mainly consider cases that are similar to those in [12] and [13].Here the presence of ρ in the formula of h makes the analysis slightly different, and while we consider several different possibilities for ρ as well, we omit large amount of situations that where the interplay between ρ and v is more intricate.
The equations for h are second-order differential equations for a nj .This makes the problem very hard in general, but we can still find some exact solutions, though we often need to restrict to special cases where there are not many terms in φ n .Even so, the number of different solutions we can find is so vast that we have to restrict to the cases that seem interesting as well as sufficiently simple.Many times, when we have an underdetermined system and are able to choose some functions arbitrarily, we choose the ones that have these second-order derivatives.Thus the system becomes a first-order system in terms of the unknown functions and is more easily solvable.However, the resulting solution formulas are sometimes quite ugly, and it is possible to analyze the systems with other approaches as well.

Columnar flows
Before studying any specific cases, we would like to show a general method of obtaining solutions to the Euler-Boussinesq equations from a specific type of flows that satisfy the homogeneous Euler equations.We consider the 3D case first, the 2D case is similar.
In many cases that we considered in [13], we could find solutions to the 3D Euler equations that were of the form ) where a is always nonzero.Flows like this feature columnar motion, where the vertical columns can stretch and contract but otherwise remain intact as the flow evolves.These are also solutions to the Euler-Boussinesq equations when density is an arbitrary function of z 3 , as ∇ρ × ∇φ 3 vanishes and condition (2.10) reduces to the usual Cauchy invariants of the Euler equations.But in this case we can look for more solutions to the Euler-Boussinesq equations by choosing ρ = c 0 z 3 and In [17] a particular solution of this form was presented in the stably stratified special case a(t) = 1, c 0 < 0. In this case the general solution of φ 3 is where the constant N , given by N 2 = −c 0 , is the Brunt-Väisälä frequency, f i are arbitrary, and φ 1 and φ 2 have to satisfy the usual 2D Euler equations.
Actually one can describe the solutions for arbitrary a.First we note that f has no effect on det(dφ) and h 3 , whereas for the time derivatives of h 1 and h 2 we have Both (h 1 ) ′ and (h 2 ) ′ vanish if and only if This is a linear ODE where z appears only as a parameter so the general solution is of the form Substituting this back to (3.3) we see that a i are two linearly independent solutions of y ′′ + qy = 0 , where q = −(a ′′ + c 0 )/a .(3.4) Note that q and hence a i are well-defined since we must anyway choose a such that a ̸ = 0 for all t.
So the solution set is parametrized by a via q, and now by classical theorems if we choose a such that q > 0 then a j are oscillating and bounded solutions.Note that there is no condition whatsoever for f 1 and f 2 , other than that they do not depend on z 3 .
Of course we can also choose c 0 = 0 to obtain new solutions to the usual Euler equations, although in that case a, a 1 , and a 2 are linearly dependent and we may assume a 2 = 0.
In the 2D situation choosing ρ = c 0 z 2 the analogous φ is given by The incompressibility condition (2.6a) already says that det(dφ) = (φ 1 ) 10 a(t) is independent of time, so we may assume by a coordinate transformation that φ 1 = z 1 /a(t).Then the formula of h gives us the same condition as in the three-dimensional case so the solution can be written as where a i are the solutions of (3.4) and f i are arbitrary.This is the only 2D solution of this form, so below we will only use the generalization method of this Section in 3-dimensional situations.
Here det(dφ) = 1 can be assumed by scaling A. By a linear transformation described in part 2 of Lemma 2.2 we may also assume that either a 21 = 0, or a 21 and a 22 are linearly independent.
If a 21 = 0, the incompressibility condition requires a 22 ̸ = 0, and in this case we must have where f is arbitrary.Then the equations are where c is a constant.If a 22 is chosen arbitrarily, then a 11 = 1/a 22 and a 12 = 1 a22 (c 0 y 2 − c)a 2 22 dt.We can also derive the Eulerian description for this solution by computing x .
If a 21 and a 22 are both nonzero and linearly independent, then ρ must be linear.In this case we may assume by a linear transformation that ρ = c 0 z 2 .If we write A as then the incompressibility condition is satisfied and the equation for h gives where now y ′ 1 = a 21 = b sin(θ).Thus we can choose b and θ arbitrarily and solve for ℓ to obtain all the solutions.In Eulerian coordinates the solution is x .

m = 4
We turn immediately to case m = 4, skipping the analysis of case m = 3 completely.There are solutions in case m = 3 for several choices of ρ but there is no room for us to consider them here.
The same is true for the Euler-Boussinesq equations with an identical proof since the incompressibility condition (2.12), which is the same for both systems, is all that is required to prove this claim.
We will consider all four cases in this Section.The equations for the entries of A are not underdetermined in the first three cases so the second-order ODEs will make them difficult to study.We still find one special solution to all these cases.The fourth case, however, is underdetermined and we are able to present a general solution for one choice of ρ.

Case 1
Let v satisfy g 13 + g 24 = g 14 − g 23 = 0, in which case v = (z 1 , z 2 , f 1 , f 2 ), where f 1 and f 2 are an anti-CR pair.It is shown in [12] how we can use linear transformations to bring A to a simpler form without loss of generality.In this case, [12, Lemma 5.3] implies that we may assume that In this form A satisfies the incompressibility condition (2.12) with which has to be nonzero in D.
Here ρ could be chosen to be for example an arbitrary linear combination of v i but the equations would become too difficult for us to find any exact solutions.Instead we choose the more simple ρ = c 0 z 2 , in which case (2.13) yields h =(Q 12 − c 0 y 1 )g 12 + Q 34 g 34 + (Q 13 − Q 24 − c 0 y 4 )g 13 + (Q 14 + Q 23 + c 0 y 3 )g 14 =c 12 g 12 + c 34 g 34 + c 13 g 13 + c 14 g 14 .
We have chosen the subscripts of the constants c ij according to the corresponding g ij .In later cases we often use this notation without further notice.
This looks like a usual Gerstner type solution from the case of Euler equations [12, Theorem 5.1], except that φ 1 and φ 2 have been scaled.In case c 0 < 0 and c 1 < 1, which describes the stably stratified situation, the particle trajectories are ellipses stretched in the vertical direction.The curves of constant density are formed by the particles whose centers of trajectory are on the same horizontal line, see Figure 4.1.
Unfortunately this equilibrium point is not hyperbolic so the linearization is inconclusive in regard to the stability of the equilibrium.However, the Jacobian has one zero eigenvalue and four purely imaginary eigenvalues if c 0 < 0 and α < c 1 < 1, where α ≈ 0.71 is the smaller positive root of q = y 8 − 12y 4 + 3.So one might suspect that the flow is "stable" for α < c 1 < 1 and "unstable" for c 1 < α.Indeed the numerical computations seem to suggest this, see Figures 4.2 It is of interest to check whether the Gerstner type solution (4.2) satisfies the free boundary condition for some choice of f 1 and f 2 .Assuming that β(s), s ∈ I ⊂ R, is a regular curve in the z plane such that ∂D = β(I), pressure should be constant along β; in other words, we need to satisfy for all t and s.Unfortunately it turns out that this condition cannot be met as we will show.The type of Boussinesq approximation we use does not matter so let us calculate the partial derivatives of p from the where 0 .This implies that f 1 and f 2 are constant in β(I) and thus constant everywhere since they are an anti-CR pair.This leads to a trivial solution as we further obtain f 1 = f 2 = 0.

Case 2
Suppose that we have the constraints g 13 = g 24 = 0, or that v = z 1 , z 2 , f 1 (z 1 ), f 2 (z 2 ) .Reducing as in [ This satisfies the incompressibility condition (2.12) with Again for density we choose the simplest situation ρ = c 0 z 2 and moreover we choose θ = 0. Inspecting condition (2.13), we first find that b 12 and b 11 /ℓ are linearly dependent, and by a linear transformation we may assume that b 12 = 0. Then we obtain for ℓ and b 11 Further let b 11 = y ′ .This can be transformed to the system 10y .
If k 0 = 0 there is an explicit solution Here k 1 does not appear in A, k 2 = 0 can be assumed by translation of t and by scaling we may take k 3 = c 5 0 .Thus this solution can be written as This solution is unstably stratified regardless of how c 0 is chosen.

Case 4
Studying the homogeneous situation in [12], we dismissed the fourth case as it only provided solutions that reduced to case m = 3.Here this does not happen; in fact solution (3.6) is already an example of this case, albeit with spatial variables named differently.
In this case we can rewrite the spatial constraints as g 23 = g 24 = 0 by changing the order of the components of v.In this case we have v = z 1 , z 2 , f 1 (z 2 ), f 2 (z 2 ) , and due to the incompressibility condition (2.12)A can be assumed to be of the form . We can choose for example b 1 and θ arbitrarily and easily solve these equations for b 2 , b 3 and b 4 , yielding We may assume c 12 = 0 without loss of generality. where As in the first 2D case (4.1), we may use linear transformations to assume without loss of generality that either ρ = f (z 3 ) + c 0 z 2 , in which case a 31 = a 32 = 0, or ρ = c 0 z 3 .Suppose first that ρ = f (z 3 ) + c 0 z 2 .Since a 31 = a 32 = 0, the rotation matrix R must be of the form for some function θ(t), and y ′ 3 = b 33 .In this case w = (0, 0, 2θ ′ ) and h can be written as h There are different ways to find partial solutions to this system, but one fairly general solution can be found by assuming only that c 12 = 0. We may then take b

m = 5
Much like for the 2D cases with m = 4, we will only find a rather small subset of the solutions in the 3D case when m ≥ 5.The three cases considered in [13,Section 5] readily yield solutions to the Euler-Boussinesq equations as well.For example in [13, Section 5.1] there was the solution φ = Av with where θ(t) is arbitrary and f 1 and f 2 are an anti-CR pair with respect to z 1 and z 2 .Recall that a solution to the Euler equations is also a solution to the Euler-Boussinesq equations if ∇ρ × ∇φ 3 = 0.Here φ 3 = θ ′ z 3 , so choosing ρ to be a function of z 3 we obtain a solution to the Euler-Boussinesq equations.The other two cases in [13,Section 5] are where b 2 ℓ ′ = 1; and , where b 2 ℓ ′ = 1.These are also solutions to the Euler-Boussinesq equations with ρ = ρ(z 3 ).
We can also find the similar type of solutions in the other three cases mentioned at the start of [13, Section 5] if we require that φ 3 = a 33 z 3 .The derivation of the following formulas is very similar to the above three cases, which were derived in [13], so we present them without proofs.In the following formulas θ is an arbitrary function with θ ′ > 0. If G 14 + G 25 = G 15 − G 24 = 0 and φ 3 = a 33 z 3 , then the solutions of the Euler equations, and the solutions of the Euler-Boussinesq equations with ρ = ρ(z 3 ), are where f 1 and f 2 are an anti-CR pair.
If G 14 = G 25 = 0, we have And if G 14 + G 25 = G 15 = 0, then we have (5.5) In the above three cases φ is of the form (3.1) so we can extend these solutions as shown in Section 3; for example the solution (5.3) can be extended as follows: φ = Av, where , where f 1 and f 2 are an anti-CR pair, and f 3 and f 4 are arbitrary.Here a i are linearly independent solutions of (3.4) with a = θ ′ .
We can add the same expression to (5.4) and (5.5) as well.
5.3 m = 6, case 1 , where f 1 is arbitrary and f 2 and f 3 are an anti-CR pair.We studied this case for the homogeneous Euler equations in [13, Section 7] and a large part of the calculations we need in the present case was already shown there.We choose ρ = i c i v i and look for solutions of the form In this case det(dφ) = f 2 100 + f 2 010 f 1 001 ̸ = 0 .We collect the rest of the constraints from the expressions of h j , which are All the equations that were used to prove [13,Lemma 7.4] are also present here.In that Lemma we showed that from these equations it follows that where 8(θ ′ ) 3 + c 5 cos(2θ) + c 6 sin(2θ) + c 56 = 0 .One particular periodic particle path in this case is shown in Figure 5.1.This solution is neither stably nor unstably stratified in particular; instead the denser particles are periodically either at the top or at the bottom of the fluid.

and 4. 3 .
We have varied s in these examples; varying other initial values gives similar results.

where again b 11 ,
b 22 , and θ are arbitrary, b 12 = b 11 w 3 b 22 b 11 dt and a i are linearly independent solutions of (3.4) with a = 1/(b 11 b 22 ).
4.1 m = 2Let us start by studying the two-dimensional case, which is the easier case.In case m = 2 we always have v = (z 1 , z 2 ) by part 1 of Lemma 2.2, and by (2.12) and (2.13)A satisfies det(dφ) = a 11 a 22 − a 12 a 21 = 1 h = a