Towards an Analytical Understanding of Internal Wave Attractors

Time harmonic inviscid internal wave motions constrained to fully closed domains generically lead to singular velocity fields. In spite of this difficulty, several techniques exist to solve such internal wave boundary value problems. Recently it has been shown that for a domain with the shape of a trapezium, solutions can be written in terms of a double sine Fourier series. However, the solutions were found by a numerical technique and thus not all coefficients of the series are available. Unfortunately, for questions related e.g. to regularization of the inviscid singular solutions, the knowledge of the asymptotic behavior of the spectrum for large wave numbers is essential. Here we discuss solutions of internal wave boundary value problems for which the spectra are known, at least asymptotically. We further describe shortcomings of the found solutions that need to be overcome in the future. Finally, we sketch applications of the solutions in the context of viscous energy dissipation.


Introduction
In the oceans, internal waves are generated mainly due to tides, moving water across bottom topography (Vlasenko et al., 2005).Such internal waves can travel large distances but they also might become trapped between lateral and vertical boundaries of the oceans.Extending the pioneering work by Magaard (1968) to fully closed basins, Maas and Lam (1995) showed that the characteristic curves of the hyperbolic internal wave equation converge towards a limit cycle, referred to as "wave attractor".Such a wave attractor marks a singularity in the velocity field.From a mathematical point of view, the non-existence of regular solutions for inviscid internal wave boundary value problems (BVPs) reflects the Correspondence to: U. Harlander (uwe.harlander@tu-cottbus.de)fact that such problems are ill-posed: a hyperbolic equation for the streamfunction ψ has to be combined with the elliptic condition, ψ=0, along the boundary.
Several techniques exist that can solve internal wave BVPs (Swart et al., 2007).However, the solutions obtained by most of these methods cannot easily be written down in closed form.Recently, Harlander and Maas (2007) showed that for an ocean with a boundary geometry like a trapezium, solutions in the form of Fourier series can be found.The coefficients of the series were found by applying a numerical collocation technique.Subsequently, the authors showed that such solutions are useful to understand dissipative damping.However, for a full understanding, the asymptotic behavior of the spectrum for large wave numbers is essential.From numerical solutions we obtain important but incomplete information.So the question arises whether analytical wave attractor solutions can be constructed, or at least solutions that give information about the asymptotic behavior of the coefficients.Ogilvie (2005) has shown for a tilted rectangle geometry, that the total rate of energy dissipation is asymptotically independent of the viscosity.We think that such important findings might be investigated further with the help of analytical wave attractor solutions.
The present paper is organized as follows.In Sect. 2 we briefly pose the mathematical problem we want to solve.In Sect. 3 we discuss the internal wave attractor solution for a trapezium, obtained recently via a collocation technique.The description of the technique is taken from Harlander and Maas (2007).Subsequently, in Sect. 4 we construct a wave attractor solution by using boundary data that (i) show focusing, and for which (ii) the Fourier transform is known.With this information we are able to find the correct Fourier coefficients of wave attractor solutions in the limit of large wave numbers.Finally, in Sect.(2).The figure is adopted from Harlander and Maas (2007).
In practice, the series has to be truncated, then representing an approximate solution to (1) where N = N + J 1 + J 2 is the total number of discrete points where data are specified (see dotted lines along the lower part of the slope and above the two fundamental intervals in Fig. 1).Along the slope it is demanded that ψ = 0, in the two fundamental intervals we prescribe ψ z .In discrete form we obtain where k = 1, 2, y s (i) are the grid points along the slope (with z s (i) = αy s (i) + c, α and c constants), y f 1 (i) the grid points along fundamental interval 1, and y f 2 (i) the grid points along fundamental interval 2 (see Fig. 1).The vector b stands for a prescribed ψ z at the grid points along the fundamental intervals.In matrix form (4)-(5) simply read where A is a known N × N matrix, a is an unknown coefficient vector of length N , and b a known fundamental data vector of length N .To find the coefficient vector a we have to invert A.

The linearized 2-D internal wave equation
The partial differential equation governing inviscid, linear, time-harmonic, two dimensional internal wave motion is well known (Vlasenko et al., 2005).By using the transformation y=((N 2 −ω 2 )/ω 2 ) −1/2 ŷ, where ŷ ist the horizontal coordinate, N is the buoyancy frequency and ω the internal wave frequency, it reads where ψ is the streamfunction and z the vertical coordinate.M denotes the vertical lake/ocean cross section and ∂M its boundary.Note that by using the transformed coordinate y, the stratification as well as the wave's frequency are part of the horizontal coordinate.A derivation of this boundary value problem (BVP) from the Boussinesq equations can be found in standard textbooks of physical oceanography, or in the recent paper by Harlander and Maas (2007).We do not repeat the derivation here but in the following section we sketch briefly the collocation method that has been used by Harlander and Maas (2007) to solve the BVP (1).

The boundary collocation method
In contrast to other methods to solve hyperbolic BVPs, the boundary collocation method gives spectral information about the solution.This can help to better understand typical features of hyperbolic BVPs.In the following brief description of the boundary collocation method we follow closely Sect.4.1 of Harlander and Maas (2007).
The domain M we begin with is a square with width π (see Fig. 1).For this geometry, the BVP (1) is solved by where the a n are constant coefficients and n∈N .Next, M is changed by replacing one vertical boundary by a sloping wall.Fundamental intervals corresponding to a square with a sloping sidewall are shown in Fig. 1 by two bold line segments along the upper boundary (see Fig. 1).In these intervals, ψ z can be prescribed arbitrarily.This follows from the fact that webs of characteristics (i.e., characteristics connected at boundary reflections), launched from fundamental intervals have no reflection point within fundamental intervals.In other words, any characteristic can uniquely be mapped into one of the fundamental intervals.They can be found by following the characteristic starting at the lower left corner, and the characteristic starting at the bottom of the slope.The latter characteristic shows one reflection at the left boundary.
Obviously, the ψ n of (2) for a square domain do no longer satisfy the boundary conditions when a slope is introduced.Nevertheless, we can still write the solution in the form (2). In practice, the series has to be truncated, then representing an approximate solution to (1) where N=N+J 1 +J 2 is the total number of discrete points where data are specified (see dotted lines along the lower part of the slope and above the two fundamental intervals in Fig. 1).Along the slope it is demanded that ψ=0, in the two fundamental intervals we prescribe ψ z .In discrete form we obtain where k=1, 2, y s (i) are the grid points along the slope (with z s (i)=αy s (i)+c, α and c constants), y f 1 (i) the grid points along fundamental interval 1, and y f 2 (i) the grid points along fundamental interval 2 (see Fig. 1).The vector b stands for a prescribed ψ z at the grid points along the fundamental intervals.In matrix form (4)-( 5) simply read where A is a known N× N matrix, a is an unknown coefficient vector of length N, and b a known fundamental data vector of length N. To find the coefficient vector a we have to invert A.
Figure 2 displays an approximate solution of (1) for the trapezium by using ( 3) and ( 6).The solution, shown in Fig. 2a, agrees with the solution presented by Maas et al. (1997) (their Fig. 1b).The location of the square shaped wave attractor is clearly visible.We note in passing that a zero level set exists at y=π.Consequently, the solution shown in Fig. 2a is also valid mathematically for a square with length π, although, from a physical viewpoint, an oblique sidewall with focusing reflections is necessary for the existence of a wave attractor.
The spectrum, that is the a n in (3), shown in the upper part of Fig. 2b, clearly reflects the anti-symmetry of the solution.All coefficients with an odd index are zero which guarantees ψ=0 along y=z=π, i.e. the anti-symmetry of ψ with respect to the diagonal axis of the wave attractor.The first part of the spectrum, determining the largest scales of the solution, shows rather irregular oscillations of large amplitude.For larger wave numbers (n>30 say), the coefficients are organized as wave packets.They determine the scales resolved close to the wave attractor.By removing single wave packets from the spectrum, we remove some part of the fine structure of the wave attractor but do not alter the gross structure of the solution.
In the lower part of Fig. 2b we plotted a 2 m , m=2, 4, 6, • • •, together with m −2 in a log-log-scale.It appears that for larger m, the coefficients converge faster to zero than m −2 .However, as we will discuss below, this sudden increase in convergency is very likely a numerical artefact.

Towards analytical solutions
Let us now try to construct an exact solution in the form that is a solution where all coefficients a n are known and not just N. As was said above, such a solution automatically fulfills the boundary condition ψ=0 for a square with length π.

General idea
In contrast to the previous section we now give ψ z along the whole upper boundary and not just in the fundamental intervals.Let us assume where ỹ=y−L p , and L p is the position of focusing.To find the yet unknown coefficients a n , we multiply ( 12 Fig. 2 displays an approximate solution of (1) for the trapezium by using (3) and ( 6).The solution, shown in Fig. 2a, agrees with the solution presented by Maas et al. (1997) (their Fig. 1B).The location of the square shaped wave attractor is clearly visible.We note in passing that a zero level set exists at y = π.Consequently, the solution shown in Fig. 2a is also valid mathematically for a square with length π, although, from a physical viewpoint, an oblique sidewall with focusing reflections is necessary for the existence of a wave attractor.
The spectrum, that is the a n in (3), shown in the upper part of Fig. 2b, clearly reflects the anti-symmetry of the solution.All coefficients with an odd index are zero which guarantees ψ = 0 along y = z = π, i.e. the anti-symmetry of ψ with respect to the diagonal axis of the wave attractor.The first part of the spectrum, determining the largest scales of the solution, shows rather irregular oscillations of large amplitude.For larger wave numbers (n > 30 say), the coefficients are Compared to the first example, the convergency of ( 19) is even faster than 1/k.The spectrum for n → ∞ reads It is instructive to compute the solution (7) by using all coefficients (20) and not just the one with large index n.Note that we still get an analytical solution of (1) for a square domain with length π.However, away from the attractor, the data deviate from f (ỹ) given in (17).Fig. 4a shows the solution for a = π/ √ 2, obtained by inserting (20) into (7).The solution as well as the spectrum shown in Fig. 4b shows a striking similarity to the solution for the trapezium (Fig. 2).
The spectrum decays roughly with 1/n 2 and again we find a series of wave packets, responsible for the wave attractor's fine structure.Note however that the center region of the attractor cannot be visualized properly.Due to the truncation of the spectrum, the fine structure is removed from the inner part of the attractor.The same feature was already observed in the solutions found via the boundary collocation method (see last paragraph of section 3).
The attractor can visually be separated from the slowly varying 'background field' by plotting the integrand | ∇ψ | 2 of the dissipation rate that reads for Rayleigh friction where ν is the kinematic viscosity coefficient.This scale sep-Fig.3. Difference between the left-hand-and right-hand-side of (15) as a function of the wavenumber for l=1, 2, 3 (thick, dashed, thin line), where a 2 =lπ 2 /2.

First example
The perhaps most simple data useful for our purpose are given by with L p =π/2.The constant a 2 =lπ 2 /2, l=1, 2, • • •, follows from the boundary condition ψ z (y) | z=π =f ( ỹ)=0 at y=0, π.Note that for any internal wave frequency (which is hidden in the y-coordinate, see Sect. 2) we find an infinite number of possible constants a.In some sense, this reflects the known fact that the spectra of internal wave BVPs are degenerate: for any frequency, an infinite number of "modes" exists (in contrast to elliptic BVPs, where for any eigenfrequency we have just one mode).
The anti-symmetry of f ( ỹ) immediately implies To find the a n with the even indices we use the sine Fourier transform of f ( ỹ), reading (Erdélyi et al., 1954) where J 1 (•) is a member of the Bessel functions of the first kind (J ±ν (•), ν=1).It can be shown that in the limit k→∞ Inserting ( 10), ( 11), ( 14), and ( 15) into (9), and using the anti-symmetry (13), we find for n→∞ With formula ( 16) we have found the asymptotically correct spectrum of a focusing solution.Note that its decay is much faster than the decay of the spectrum found numerically for the trapezium (Fig. 2b), that decays approximately with n −2 .Before we discuss a second example we note that (15) holds approximately even for rather small k. Figure 3 shows the difference between the left-hand-and right-hand-side of (15) for l=1, 2, 3 (thick, dashed, thin line).For l=1, e.g., the error decays faster than 1/k.

Second example
Let us now consider data that imply a similar spectral decay than the trapezium shown in Fig. 2. A proper choice is For this function, the sine Fourier transform reads (Erdélyi et al., 1954) Again we use the fact that for ) Compared to the first example, the convergency of ( 19) is even faster than 1/k.The spectrum for n→∞ reads It is instructive to compute the solution ( 7) by using all coefficients (20) and not just the one with large index n.Note that we still get an analytical solution of (1) for a square domain with length π.However, away from the attractor, the data deviate from f ( ỹ) given in (17).Figure 4a shows the solution for a=π/ √ 2, obtained by inserting ( 20) into (7).The solution as well as the spectrum shown in Fig. 4b shows a striking similarity to the solution for the trapezium (Fig. 2).The spectrum decays roughly with 1/n 2 and again we find a series of wave packets, responsible for the wave attractor's fine structure.Note however that the center region of the attractor cannot be visualized properly.Due to the truncation of the spectrum, the fine structure is removed from the inner part of the attractor.The same feature was already observed in the solutions found via the boundary collocation method (see last paragraph of Sect.3).
The attractor can visually be separated from the slowly varying "background field" by plotting the integrand | ∇ψ | 2 of the dissipation rate that reads for Rayleigh friction where ν is the kinematic viscosity coefficient.This scale separation is demonstrated in Fig. 4c.It is obvious that the inner Adv.Geosci., 15, 3-9, 2008 www.adv-geosci.net/15/3/2008/region of the attractor, where dissipation becomes large (and even singular), is removed by truncating the series (7).Seen from a more physical point of view, Fig. 4c demonstrates that in a viscous fluid, for which the fine structure close to the internal wave attractor is damped out, energy dissipation and mixing will take place at a certain distance away from the attractor.This distance depends not only on the size of ν, but also on the geometry of the basin.E.g., for the trapezium, the distance is smaller than for the case shown in Fig. 4.

Conclusions
Singular flow components can be observed in freely propagating wave beams of locally excited internal waves (Chashechkin and Prihodko, 2007), as well as along the limit cycles of trapped internal waves that are focused due to multiple reflections from oblique boundaries (Maas et al., 1997).
The limit cycles are called wave attractors and the solutions of the corresponding internal wave boundary value problem are referred to as wave attractor solutions.In the present paper we have demonstrated that wave attractor solutions can be constructed for a square with length π by prescribing von Neumann boundary conditions ψ x =f (x) (x is a substitute for y or z) at one side of the square.The solutions are written in the form of a double sine Fourier series.For n→∞, the coefficients a n , n=1, 2, • • • of the solution converge towards the coefficients of the sine Fourier transform of f .In other words, the "focusing" of the solution (i.e., the wave attractor) corresponds with the focusing prescribed by f .It was pointed out by Harlander and Maas (2007) that for an understanding of the effect of viscosity on the structure of the wave attractor, an analytic expression of the spectrum of the attractor would be helpful.With ( 16) and (20) two examples of such expressions have come in.Under the assumption that the time scale of the damping is much larger than the time scale of the internal wave period, Harlander and Maas (2007) found that the coefficients are damped by the factor exp(−ν( √ 2n) 2p t), where p=1, 2, • • • specifies the damping characteristic (p=0, Rayleigh friction, p=1 momentum diffusion, p>1, hyper-diffusion).With the expressions for the coefficients the effect of dissipation can now be computed for any ν, p, and n.
With the second example, we have discussed a solution that shows a similar spectral decay than the solution for the trapezium.It would be interesting to find data f that lead to a weaker damping.Such a solution would raise questions about the convergency of the time mean dissipation rate.For Rayleigh friction, e.g., the spectrum has to decay at least with n −2 to obtain a finite time mean dissipation rate (Harlander and Maas, 2007).
The major shortcoming of the solutions presented is that we do not have an expression for the focusing part of the boundary.Mathematically, the solutions are correct, however, from a physical point of view internal waves need to   reflect from an oblique or curved boundary to become focused.The solution shown in Fig. 2a holds for a square with length π but there exists also a zero level of ψ that defines . 5. Enlargement of the wave attractor reflection region of Fig.
pe of the boundary segment that leads to focusing is imitly given by f but explicitly unknown.Fig. 5 shows a gnification of the region of reflection at the wall y = π.highlighted two zero levels that seem to be connected at z) ≈ (2.6, 0.53).The red and green curve correspond to tened circles defined by parameter α defines the flatness of the circle.For red curve we use the plus sign and (α, s, r, y 0 , z 0 ) = 15, 3, π/4, π, 0), for the green curve the minus sign and s, r, y 0 , z 0 ) = (0.4, 1, π/2, π/2, π/2).Closer to the attor, the zero level (shown by blue lines) can accurately be cribed by hyperbolas the upper part we use the positive sign and b, y 0 , z 0 ) = (0.35, 1/2, −π/2, −π −a), for the lower part negative sign and (a, b, y 0 , z 0 ) = (1/2, 0.35, a, −π/2).w we can approximately follow the sloping boundary, that esponsible for the focusing: starting on the red curve at 0 and following the curve until it is connected with the en curve.Then switch to the green curve and follow it uprd towards the blue curve.Finally, follow the blue curve makes it difficult to resolve its neighboring region.Not ther that attractor solutions do not depend much on the cise shape of the focusing boundary as long as the frequ is not at the edge of an attractor interval.Attractor solu at such critical frequencies have a large or small aspe tion (i.e. the parts of the attractor with positive and neg slope have different length).Clearly, this does not ho the case considered.As long as the frequency is in the c of an attractor interval, hyperbolic BVPs behave simila elliptic ones, that is small boundary deformations have small impact on the solutions.
It should be noted that the strategy sketched in the pr paper might also be useful to find wave attractor solu for more general models.For instance, Harlander and (2007) constructed solutions of the Tricomi equation b ing an expansion in terms of Airy functions that also fo complete system of orthogonal functions.
A final note is concerned with the sharp drop in the trum of the trapezium for large n, shown in Fig. 2b.In to improve the boundary collocation method used to o the spectrum, is is important to understand this featur view of the spectrum of example two (Fig. 4b) that not show this feature, it is likely that the drop is a num artefact that results from the truncation of the series (7 the subsequent numerical matrix inversion that gives th efficients.However, lacking an analytical expression o spectrum or its asymptotic behavior, this conjecture is yet.the oblique side wall responsible for the wave focusing.The shape of the boundary segment that leads to focusing is implicitly given by f but explicitly unknown.Figure 5 shows a magnification of the region of reflection at the wall y=π.We highlighted two zero levels that seem to be connected at (y, z)≈(2.6,0.53).The red and green curve correspond to flattened circles defined by The parameter α defines the flatness of the circle.For the red curve we use the plus sign and (α, s, r, y 0 , z 0 ) = (0.15, 3, π/4, π, 0), for the green curve the minus sign and (α, s, r, y 0 , z 0 ) = (0.4, 1, π/2, π/2, π/2).Closer to the attractor, the zero level (shown by blue lines) can accurately be described by hyperbolas For the upper part we use the positive sign and (a, b, y 0 , z 0 ) = (0.35, 1/2, −π/2, −π−0.35), for the lower part the negative sign and (a, b, y 0 , z 0 ) = (1/2, 0.35, 1/2, −π/2).Now we can approximately follow the sloping boundary, that is responsible for the focusing: starting on the red curve at z=0 and following the curve until it is connected with the green curve.Then switch to the green curve and follow it upward towards the blue curve.Finally, follow the blue curve until y=π is reached.Note first that the focusing rate is decreasing the closer we are to the attractor.This is different from the case with a linear slope (Harlander and Maas, 2007), for which the focusing rate is constant.The space dependent focusing is responsible for the rather broad not-resolved region surrounding the attractor (see Fig. 4).Focusing is already very strong at a distance of the attractor.Clearly, this makes it difficult to resolve its neighboring region.Note further that attractor solutions do not depend much on the precise shape of the focusing boundary as long as the frequency is not at the edge of an attractor interval.Attractor solutions at such critical frequencies have a large or small aspect ration (i.e. the parts of the attractor with positive and negative slope have different length).Clearly, this does not hold for the case considered.As long as the frequency is in the center of an attractor interval, hyperbolic BVPs behave similar like elliptic ones, that is small boundary deformations have just a small impact on the solutions.
It should be noted that the strategy sketched in the present paper might also be useful to find wave attractor solutions for more general models.For instance, Harlander and Maas (2007) constructed solutions of the Tricomi equation by using an expansion in terms of Airy functions that also form a complete system of orthogonal functions.
A final note is concerned with the sharp drop in the spectrum of the trapezium for large n, shown in Fig. 2b.In order to improve the boundary collocation method used to obtain the spectrum, is is important to understand this feature.In view of the spectrum of example two (Fig. 4b) that does not show this feature, it is likely that the drop is a numerical artefact that results from the truncation of the series (7) and the subsequent numerical matrix inversion that gives the coefficients.However, lacking an analytical expression of the spectrum or its asymptotic behavior, this conjecture is open yet.

Fig. 1 .
Fig.1.The domain M , given by a square with a sloping sidewall (dashed lines).The two fundamental intervals are shown by bold solid lines along the upper boundary.The set of points where data are prescribed are indicated by dotted lines along the fundamental intervals and the lower part of the slope.In total there are N + J 1 + J 2 data points.The wave attractor is shown by the dashed square.The figure is adopted fromHarlander and Maas (2007).

Fig. 1 .
Fig.1.The domain M, given by a square with a sloping sidewall (dashed lines).The two fundamental intervals are shown by bold solid lines along the upper boundary.The set of points where data are prescribed are indicated by dotted lines along the fundamental intervals and the lower part of the slope.In total there are N +J 1 +J 2 data points.The wave attractor is shown by the dashed square.The figure is adopted fromHarlander and Maas (2007).