A conformal mapping approach to modelling two-dimensional stratiﬁed ﬂow

Herein we describe a new approach to modelling inviscid two-dimensional stratiﬁed ﬂows in a general domain. The approach makes use of a conformal map of the domain to a rectangle. In this transformed domain, the equations of motion are largely unaltered, and in particular Laplace’s equation remains unchanged. This enables one to construct exact solutions to Laplace’s equation and thereby enforce all boundary conditions. An example is provided for two-dimensional ﬂow under the Boussinesq approximation, though the approach is much more general (albeit restricted to two-dimensions). This example is motivated by ﬂow under a weir in a tidal estuary. Here, we discuss how to use a dynamically-evolving conformal map to model changes in the water height on either side of the weir, though the example presented keeps these heights ﬁxed due to limitations in the computational speed to generate the conformal map. The numerical approach makes use of contour advection, wherein material buoyancy contours are advected conservatively by the local ﬂuid velocity, while a dual contour-grid representation is used for the vorticity in order to account for vorticity generation from horizontal buoyancy gradients. This generation is accurately estimated by using the buoyancy contours directly, rather than a gridded version of the buoyancy ﬁeld. The result is a highly-accurate, eﬃcient numerical method with extremely low levels of numerical damping. © 2023 The Author(s). Published by Elsevier Inc. This is an open access article under the CC BY license (http://creativecommons .org /licenses /by /4 .0/).


Introduction
Density stratified flows occur widely in the environment, particularly in the Earth's atmosphere and oceans but also in other planetary atmospheres [20,16].Stratification provides a natural restoring force, tending to relax flows back to a stable hydrostatic equilibrium.Nonetheless, there are many factors which may disrupt this equilibrium, including surface forcing in the oceans, tides over topography, etc.Such factors commonly result in long-lived propagating internal waves (rather than decay to hydrostatic equilibrium), and these waves can reach large amplitudes, up to hundreds of metres in the oceans (see e.g.Ostrovsky and Stepanyants [15], Helfrich and Melville [12]).They impact natural and man-made submarine structures and contribute significantly to the mixing of heat, salt, chemical and biological constituents.
Over the past few decades, there has been much research conducted to investigate the characteristics, stability, evolution and mixing of internal waves (see Carr et al. [3] and references therein).Experimental, theoretical and computational approaches all have contributed to our present knowledge of the dynamics of these waves.The present work focuses on a new computational approach enabling one to simulate such waves, and density stratified flows more generally, in arbitrary two-dimensional domains.Such domains may, for example, include bottom topography to examine how internal waves are generated, transformed or scattered from topographic features, a process thought to be generic in the world's oceans (see e.g.Trulsen et al. [19]).
The new approach described in this paper generalises a previous one used to study internal waves in rectangular domains (periodic or aperiodic) [14].Both approaches use the same underlying method, 'contour advection' [5,6], in which the prognostic fields of density (or buoyancy) and vorticity are represented by a set of material contours, with additional auxiliary gridded fields used for vorticity to allow for vorticity generation from density gradients.The advecting velocity field is obtained by inversion of the vorticity field on a grid (here by spectral methods) after the vorticity contours are converted to gridded values by a fast-fill algorithm developed in Dritschel and Ambaum [5].Contour advection allows one to resolve a very wide range of scales efficiently and accurately, as has been shown repeatedly (see e.g.Dritschel and Tobias [7] and references therein), and in the present context has led to new insights into the complex nonlinear processes involved in internal wave breaking [3].
The generalisation to arbitrary two-dimensional domains makes use of a conformal map to a rectangle, wherein the equations are solved in transformed coordinates.Remarkably few changes are required to the basic algorithm, for a static domain, and the efficiency is comparable to that enjoyed by the original algorithm.For a time-dependent domain, additional considerations arise, the most important of which is the need to re-generate the conformal mapping at each time step (or when the domain has changed sufficiently).The efficiency of the algorithm in this case is limited by the cost of re-generating the conformal map, which at present is the dominant cost (we use the SCtoolbox of Driscoll and Trefethen [4]).Nonetheless, we provide the algorithm since it may become viable if improvements in the conformal mapping can be found.
The structure of the paper is as follows.In section 2, we review the governing flow equations, the changes arising from conformal mapping, and the so-called 'inversion problem' in which one recovers the velocity field from the vorticity field and boundary conditions.We also briefly review the numerical method, which is closely similar to that used in [14].In section 3, we illustrate an example of an exchange flow, a classical and extensively studied problem in density-stratified flows (see Shin et al. [18] and references therein).Here, we look at the much less well studied effect of the domain geometry, in particular obstacles.To this end, we consider a dam-break separating two regions of fluid of different density in a complex domain containing a weir, a sloping bottom boundary and a flat top boundary having different heights either side of the weir.This domain is relevant to investigating tidal flow in an estuary, though one would have to go further and consider time-dependent boundaries with inflow/outflow conditions at the sea and river 'edges'.There are many, potential applications of this method that relate specifically to stratified flows within complex domains.For example, oceanic flow over shelves and slopes [2], topographic effects on exchange flows and gravity currents [2], vertical shear processes in river plumes [1], vortex motion within a circular bay [17], urban atmospheres in complex terrain [8], estuarine circulation [11] and subglacial plumes [13], to name just a few.Further discussion may be found in the conclusions in section 4.

Methodology
We consider a two-dimensional stratified flow under the Boussinesq approximation.Taking x to be the horizontal coordinate and y to be the vertical coordinate with x = (x, y), the flow is governed by the following system of equations in the physical domain (see e.g.Fig. 1):

Dζ
where b(x, t) is the buoyancy, ζ(x, t) is the vorticity component pointing out of the plane of view, u = (u(x, t), v(x, t)) is the two-dimensional velocity field, and D/Dt = ∂/∂t + u •∇ is the material derivative (see e.g.King et al. [14] and references therein).The vorticity equation comes from taking the curl of the momentum equations (thereby eliminating pressure), while the buoyancy equation expresses conservation of mass in an incompressible fluid for which ∇ • u = 0.The velocity field may be found in terms of a scalar streamfunction ψ via which then automatically ensures ∇ •u = 0. Then using the definition of vorticity ζ = ∂ v/∂ x − ∂u/∂ y, we arrive at Poisson's equation for ψ : This may be solved once we specify the boundary values ψ = ψ , which are typically zero for a static, closed domain.Below, we describe how to solve these equations after a conformal mapping.We consider first the simpler static domain case, then discuss the new features required for a moving domain.

Static domain
To simplify notation in the remainder of this paper, we henceforth use (X, Y ) to denote a point in the physical domain, and (x, y) to denote a point in the conformal domain.The conformal domain is taken to be a rectangle as in Fig. 1, and without loss of generality we let 0 ≤ x ≤ w x and 0 ≤ y ≤ w y .Other conformal domain shapes can be chosen (e.g. a circle), but a rectangle is convenient for the present purposes.

Preliminaries
A conformal map consists of the relations giving a point (X, Y ) in the physical domain in terms of the conformal coordinates (x, y).The map must satisfy the Cauchy relations and moreover (for consistency) must satisfy Laplace's equation in the conformal domain: To solve Laplace's equation, normally one specifies boundary values, but here both equations must be satisfied simultaneously.The method of solution depends on whether the physical domain boundary is taken to be a polygon or a piecewise smooth curve.Here, we consider the first case in order to exploit the SCtoolbox software package in MATLAB [4].This package enables one to construct general conformal maps for a wide range of polygonal boundaries.Moreover, it provides the derivatives ∂ X/∂ x etc. required to recast the fluid dynamical equations ( 1)-( 3) in conformal coordinates.
The advantage of a conformal coordinate transformation is that the material derivative D/Dt remains unchanged in form, as long as the velocity field is defined by u = dx/dt as opposed to U = dX /dt in the physical domain.Poisson's equation (3) simply picks up a conformal factor λ: ∇ 2 ψ = λζ (7) where and ∇ is the gradient operator now in the conformal domain.Note: λ integrated over the area of the conformal domain gives the area of the physical domain.The X and Y derivatives of any field f transform as where we have used the chain rule together with the Cauchy relations (5) to solve for ∂ f /∂ X and ∂ f /∂ Y in (9).Finally, the physical velocity field U also transforms by the chain rule: Using ( 9) and (11) to transform the incompressibility relations (2), we obtain simply Note: ∇ • u = 0; the conformal transformation does not preserve incompressibility (or local area).

Velocity inversion
The core part of the numerical method consists of solving Poisson's equation (7), needed to obtain the velocity field u from (12).The rectangular shape of the domain allows one to exploit accurate and efficient spectral methods (Fast Fourier Transforms or FFTs) for this purpose.Considering a grid with n x and n y divisions in x and y, fields are expanded either as a finite cosine or a sine series in x and y, with discrete wavenumbers k m = mπ /w x and l n = nπ /w y for integers m ∈ [0, n x ] and n ∈ [0, n y ].
While in a static domain it is sufficient to impose ψ = 0 on the boundaries, we here consider the more general situation in which ψ = ψ on the boundaries, which is relevant to the moving domain case with possible inflow and outflow conditions discussed in the following section.
As ζ and λ in (7) can in general have arbitrary boundary values, the source term S = λζ is expanded in a cosine series in x and in y.We therefore seek a streamfunction ψ of the form where the components ψ p , ψ c and ψ h satisfy the following conditions: • ∇ 2 ψ p = S and ψ p has the same cosine series as S but generally does not satisfy ψ p = ψ on the boundaries; • ∇ 2 ψ c = 0 and ψ c = ψ ≡ ψ − ψ p on all corners, i.e. (13) where respectively ψ 00 , ψ 10 , ψ 01 and ψ 11 are the values of ψ at the lower-left, lower-right, upper-left and upper-right corners of the conformal domain, and where ξ = x/w x and η = y/w y are scaled coordinates in the range [0, 1]; This construction ensures that ψ = ψ on all boundaries and is compatible with a Fourier series solution, as explained next.
Since the auxiliary field ψ = ψ − ψ c vanishes at each corner by construction, the boundary values of ψ can be represented as a sine series (in either x or y) on each edge.Specifically, where the coefficients A ± m and B ± n are in practice determined by FFTs.The Nyquist frequency m = n x or n = n y is absent because the sine function is zero at all grid points.Consider then the following proposed solution for ψ h : where with α m = k m w y and β n = l n w x .Each term in both sums in ( 15) is a solution to Laplace's equation, and hence ∇ 2 ψ h = 0 as required.Moreover, the functions f ± m and g ± n are either 0 or 1 on the boundaries, and it is simple to check that ψ h = ψ on all edges, as required.Thus the combined solution ψ = ψ p + ψ c + ψ h satisfies ∇ 2 ψ = S together with the required boundary conditions.
The velocity field u is found by differentiating ψ in (12): where the partial derivatives of ψ p are evaluated in spectral space by wavenumber multiplication followed by inverse FFTs (taking into account the appropriate changes from cosine to sine series), while the partial derivatives of ψ c are computed analytically as and finally those of ψ h are computed semi-analytically followed by an appropriate inverse FFT: where 1 (21)

Moving domain
We next extend the method described in the previous subsection to allow the shape of the physical domain to vary in time and to allow general inflow/outflow boundary conditions.This permits one to study tides in a river estuary, including changing upstream flow conditions.
Here we restrict attention to the case where the river and sea surfaces, at Y = Y R (t) and Y = Y S (t) respectively (see Fig. 1, top panel), are the only parts of the domain which are varying in time t.
We generalise the conformal map (4) to where x is a scaled x coordinate defined by x = κ(t)x where κ(t) ≡ w x0 /w x (t) is a dilation factor and w x0 is a reference conformal domain length in x.The reason for this construction is that, in the conformal mapping solution, w x varies with the shape of the physical domain, whereas w y does not.Moreover, it is convenient for imposing boundary conditions to use a static, modified, conformal domain in x, y coordinates (see below).
The objective is to obtain the velocity û = dx/dt, v = d y/dt in the scaled conformal domain in terms of the streamfunction ψ .To this end, we start with the definition of the velocity in the physical domain and use the chain rule: where new terms appear due to the time-dependent map.Next using together with the Cauchy relations (5), we obtain Inverting these for û = dx/dt and v = d y/dt, we find where λ = ∇Y 2 is the conformal factor as before.
Next we consider derivatives of the streamfunction ψ .By the chain rule and similarly where we have used the Cauchy relations and the definitions of U and V .However, note that the groups of terms representing ∂ψ/∂x and ∂ψ/∂ y appear also in (25).Thus we can simplify the latter to obtain the sought-after solution for û and v: Comparing with the static domain case (12), besides the scaling factor κ(t), new terms arise due to the time-dependent conformal map.These terms ensure that there is no normal flow across any impermeable boundary, i.e. that the normal velocity in the static, modified, conformal domain is identically zero there.To see this, consider the situation depicted in Fig. 1 (top panel) where the sea surface at Y = Y S (t) and river surface at Y = Y R (t) may generally vary in time.These boundaries correspond to x = 0 and x = w x0 respectively in the (scaled) conformal domain.We next show that û vanishes, as required, on these boundaries.First consider the sea surface.There, X does not vary in time while the upward velocity of the sea surface.Hence, according to (28), using the Cauchy relations.However, along this boundary, ψ must be chosen to be consistent with the velocity boundary condition there, namely V = ∂ψ/∂ X = V S (t).This implies ψ = V S (t)X + C (t), where C (t) is an ignorable constant (as it does not alter the velocity field).Hence, and this is seen to exactly compensate the same term in (29), proving that û = 0 there as required.The same argument can be used to show that û = 0 along the river surface at Y = Y R (t).

Results
In this section, we perform a series of tests to demonstrate the validity and accuracy of the methods developed in the previous section.The core numerical problem is obtaining the velocity field from a given interior vorticity field and prescribed boundary velocities.We demonstrate that the method is second-order accurate in grid spacing, as expected.We go on to illustrate an example of the flow evolution resulting from a 'dam-break', in a static domain.The moving domain case is deferred to a future study due to the significant added complexity associated with the inflow/outflow boundary conditions.

Numerical tests
We consider the physical domain illustrated in the top panel in Fig. 1.Taking equally-spaced grid-lines in the conformal domain (the rectangle in the bottom panel of Fig. 1), the image in the physical domain is shown in Fig. 2 at coarse resolution.

Inversion in a static domain
A core part of the numerical algorithm, indeed one which requires the greatest computational resources, is determining the velocity field U from the vorticity field ζ and prescribed boundary conditions, a process called 'inversion'.In the domain shown in Fig. 2, we take the vorticity field to be the Gaussian function centred on a prescribed point X 0 and having a characteristic radius R. We take X 0 = (1.5, 1.0) and R = 0.5, though these specific choices do not affect the rates of convergence of our results.This vorticity field is shown in Fig. 3(a).
Consider first a fixed, static domain with streamfunction ψ = constant (= 0 without loss of generality) along the domain edges.Fig. 3 shows the streamfunction and velocity components at the default resolution, n x = 400 and n y = 200.The flow field is consistent with the prescribed boundary conditions, namely that U is everywhere tangent to the boundary.Notably, the flow field is well behaved near the lower rounded portion of the weir, where the curvature of the boundary is high.To assess accuracy, we examine next the convergence of the solution with increasing resolution.As an exact solution is not known, we use the solution on a 1600 × 800 grid as a reference, and compare it to solutions generated on a series of coarser grids having 2 m fewer grid points in each direction for m = 1, 2 and 3. Two measures of error are considered, each for the three fields ψ , u and v.The first measure is the r.m.s.error computed using all common grid points (i.e.every eighth grid point when comparing 200 × 100 with 1600 × 800).The second measure is the maximum absolute error over all common grid points.
The results are shown in Fig. 4 in logarithmic scaling.The dashed lines show linear least-squares fits of the errors plotted in solid lines.For ψ , the r.m.s.error decays like n −3.07±0.17 x while the maximum error decays more slowly, like n −1.90±0.09 x .For u and v, the errors decay like n −1.75±0.04 x while the maximum error again decays more slowly, like n −1.41±0.17 x .The error in ψ is mainly concentrated in the lower right corner of the domain, where the grid cell size is largest in the physical domain (see Fig. 2), but error also occurs very close to the upper right corners of the domain to the right of the weir, see Fig. 5.The error in u is mainly located at the upper right edge of the domain, while that in v is located at the two right vertical portions of the domain (not shown).
The main conclusion is that the solution converges with increasing resolution.In an r.m.s.measure, the convergence rate is nearly quadratic for the velocity components, and cubic for the streamfunction.Largest errors are found at edges or in regions occupied by the largest grid cells in the physical domain.

Inversion in a moving domain
Consider next inversion in a moving domain where the streamfunction ψ may vary along the domain edges.We examine a moment of time when the domain shape is the same as in the previous case.We suppose that there is a uniform flow U R at the river edge (independent of Y and less than zero if flowing from right to left in Fig. 1).Furthermore the upper surfaces at Y = Y R and Y S are moving at speeds V R and V S respectively.The flow through the sea edge U S (also assumed uniform for simplicity) is determined by mass conservation (or volume conservation for the incompressible fluid considered).Let Y B R be the Y coordinate at the bottom of the river edge on the right, and Y B S be the Y coordinate at the bottom of the sea edge on the left.Also, let the river edge be at X = X R and the sea edge be at X = X S .Finally, let the right (sea) edge of the weir be at X = X W S and the left (river) edge of the weir be at X = X W R .Then conservation of volume determines U S from These speeds on the various edges are needed to determine ψ as a continuous function around the entire domain boundary.Along the (impermeable) bottom of the domain, we take ψ = 0 without loss of generality.Then along the sea edge, by continuity ψ = U S (Y B S − Y ).Going clockwise in Fig. 1, along the sea surface ψ = V S (X − X W S ) + C W where is the constant value of ψ on the weir (another impermeable boundary).Continuing, on the river surface ψ Note, we are tacitly assuming that the free surfaces remain flat, free of waves or turbulence (e.g. as induced by a hydraulic jump).This requires sufficiently low flow speeds (small Froude numbers) in practical applications.Otherwise, a two-dimensional model like the one developed would not be suitable, or at least not accurate.
To test the inversion procedure, we choose (arbitrarily) U R = −0.09,V R = 0.12 and V S = −0.06,then determine U S ≈ −0.1053076 from (30) and the domain shape parameters A similar spatial pattern is found when comparing other grid resolutions.coordinate along the left, bottom and right edges of the physical domain monotonically increases from X R to X S .Along the left edge, X remains very close to X S until abruptly increasing (by a factor of about 10 6 ), and the last point for which X − X S < 8 × 10 −8 determines Y B S .Likewise, along the right edge (moving downwards), X remains very close to X R until abruptly decreasing, and the last point for which X R − X < 8 × 10 −8 similarly determines Y B R .Only at the highest resolution 1600 × 800 does the abrupt decrease in X near the river edge X R not occur at the same Y location -it occurs one grid point lower.However, to compare results between different resolutions, we need to ensure Y B R and Y B S are identical.This is why we chose the specific tolerance 8 × 10 −8 in the above procedure.
Fig. 6 shows the streamfunction and velocity components at the default resolution, n x = 400 and n y = 200, for the same vorticity field shown in Fig. 3(a).Now, where fluid is moving across a boundary, the computed solution matches the prescribed boundary velocities to high accuracy (with errors of O(10 −7 ) in U R , U S , V R and V S ).The small errors come from using a finite Fourier series (finite resolution) to represent ψ .Elsewhere, the fluid velocity is tangent to the boundaries, as required.Note that relatively high flow speeds occur near the lower part of the weir where the flow rapidly changes direction to remain tangent to the boundary.If the weir had higher curvature at its lower end, flow speeds would be correspondingly larger, possibly leading to flow separation in reality (i.e. with weak viscosity).
To assess accuracy, we examine next the convergence of the solution with increasing resolution, comparing as before with the solution on the highest resolution (1600 × 800) grid.The results are shown in Fig. 7 in logarithmic scaling.The dashed lines show linear least-squares fits of the errors plotted in solid lines.For ψ , the r.m.s.error decays like n −2.00±0.17 x while the maximum error decays more slowly, like n −1.17±0.20 x .For u and v, the errors decay like n −1.25±0.07 x and n −1.60±0.10 x while the maximum error again decays more slowly, like n −0.264±0.002 x and n −0.72±0.22 x .The convergence is not as rapid as in the static domain case, cf.Fig. 4.This is likely due to the discontinuity in velocity along the lower edge of the conformal domain shown in Fig. 1, as well as in the lower left and right hand corners.This is confirmed in Fig. 8, which shows the differences in ψ , U and V between the solutions on the 400 × 200 and 1600 × 800 grids.The errors are mainly located in the lower left and right corners of the domain.The error in ψ vanishes at the boundaries, as enforced by the solution procedure.The error in U has a stippled pattern, likely Gibbs fringes associated  with the discontinuity in U at the corners of the domain.The error in V is more concentrated near the corners and of larger magnitude (four times greater than U ).
Nevertheless, the solution converges with increasing resolution, albeit more slowly due to the inflow/outflow boundary conditions.In an r.m.s.measure, the convergence rate is between linear and quadratic for the velocity components, and quadratic for the streamfunction.Largest errors occur near the corners of the domain where the velocity components are discontinuous in the conformal domain.

Contour advection
In the static domain case, the flow evolution is computed in the closed rectangular conformal domain exactly as described in Carr et al. [3], section 3.3.Both vorticity ζ and buoyancy b are represented a set of contours which are advected by the velocity field u, ζ is additionally represented by two residual gridded fields designed to uptake the vorticity source ∂b/∂x in (1) and improve total energy conservation -full details of the method may be found in Dritschel and Fontane [6].Contour advection allows one to reach extraordinarily high resolution at modest computational expense, or to achieve accuracy much more efficiently than standard numerical methods such as pseudo-spectral.This is demonstrated in Dritschel and Tobias [7] for 2D magneto-hydrodynamics, a system which shares a similar property that the vorticity source is typically fine-scale.

Illustration: the dam break problem
We next illustrate the complete method for an initial condition consisting of two regions of different density (buoyancy) either side of a vertical line in a flow at rest.Subsequently, the heavier fluid will slump and try to occupy the lower part of the domain, lowering the potential energy.The kinetic energy will grow to compensate, eventually leading to a complex turbulent flow and strong mixing.
We employ the same irregular domain used throughout this paper.In this domain, which extends from X = 0 to X = 2, we position the initial density transition line at X = X 0 = 1.5.To the left, we take the buoyancy b = 1, while to the right we take b = 0 (the heavier fluid).A small transition is introduced near X = 1.5 of width w = 0.02 so that the initial buoyancy has the form b = −0.5 erfc((X − X 0 )/w) where erfc(z) is the complementary error function.The chosen location of the dam break leads to particularly strong turbulent mixing, and also illustrates how the weir largely confines the turbulence to the right-hand portion of the domain.where x and y are the grid lengths in the rectangular conformal domain.Above, the number 0.7 is the CFL limit.This ensures both numerical stability for the spectral evolution of vorticity ζ and accuracy in the advection of the b and ζ contours.
Buoyancy contours are equally spaced in b, with 100 contour intervals across the range b min to b max (here 0 to 1).This is more than twice the recommended minimum number of intervals [6], but is used here to better represent a continuous where f denotes the domain average of f .In this way, ζ increases as the flow develops localised high vorticity values, but the increase is slower than when taking the extreme values of ζ .This choice of ζ allows the number of vorticity contour levels to grow to adapt to increasing flow complexity, thus maintaining an adequate resolution of ζ at all times.Following Carr et al. [3], the same form and magnitude of hyperviscosity is applied to the vorticity equation, to control the inevitable cascade of vorticity to small scales.Without some sort of damping, vorticity can accumulate at small scales, whereas ideally vorticity should cascade to ever smaller scales, below those resolved and be removed.Hyperviscosity is not benign, but it is necessary to prevent this accumulation.(Notably, no damping of any kind is applied to buoyancy.)Hyperviscosity is included by adding a term ν h (−∇ 2 ) p ζ to the right hand side of the vorticity equation in (1).Here, as in Carr et al. [3], we take p = 3 and take the hyperviscosity coefficient ν h proportional to the instantaneous maximum value of |ζ |, specifically ν h = C h |ζ | max /k 2p max where k max is the maximum wavenumber magnitude in the conformal domain, and C h = 80.In spectral space, hyperviscous damping is simply incorporated in the time integration by dividing the undamped spectral tendency by the factor 1 + ν h k 2p t/2.
We turn next to the results.Fig. 10 shows the fields of buoyancy b and vorticity ζ at a few selected times, for a simulation conducted at the default resolution 400 ×200.Note that the effective resolution is 16 times finer in each direction [6].At early times, the heavier fluid slumps below the lighter fluid, creating a strong shear layer (high vorticity) within the region of strong buoyancy gradients.Recall that vorticity is generated by horizontal buoyancy gradients, see (1).This shear layer then destabilises despite being stretched, rolling up into a series of Kelvin-Helmholtz billows.These interact with each other and with the domain boundaries, and rapidly grow in complexity.This causes the two regions of different buoyancy (density) to partially mix.However, this mixing is largely confined to the region to the right of the weir.Note that the vorticity spreads over an increasing area of the domain, in response to the buoyancy mixing.Maximum vorticity values grow significantly, reaching over 200 in magnitude by the last time shown.(On a coarser grid of dimensions 200 × 100, the maximum vorticity magnitude reaches approximately 100.) Fig. 11 shows the evolution of the energy components, kinetic (red) and potential (blue), together with their sum (black).
The dashed curves show results for a coarser 200 × 100 simulation.Initially, the potential energy is taken to be zero without loss of generality.At first, this decreases sharply as the heavier fluid rushes under the lighter fluid: the centre of mass lowers.Then there is a rebound and oscillation around the hydrostatic equilibrium configuration, in which the region heavier fluid lies below the lighter one.The kinetic energy exhibits the opposite behaviour, consistent with conservation of total energy.Some energy is dissipated here due to the intense small-scale mixing, a generic feature of density-stratified turbulent flows.Notably, the results depend only weakly on resolution, especially at early times (t < 3).
Fig. 12 shows the evolution of two buoyancy norms, the linear norm b (red) and the quadratic b 2 (black) -again for the default resolution and a coarser resolution (solid and dashed lines, respectively).Here, the norms are defined as domain integrals over the physical domain.The method exactly conserves b -effectively the total mass, an important global quantity.The (monotonic) decay of b 2 is the result of mixing: in the method this is the result of removing very thin filaments of buoyancy by 'contour surgery' at 1/16 of the grid cell size in the conformal domain.Notably, the mixing is not strongly dependent on resolution.This is because mixing is especially vigorous in density-stratified flows, due to the   production of intense vorticity (see below).Note that mixing is slightly delayed at higher resolution (the strong downturn in the black solid occurs later in time), consistent with the greater time it takes for scales to cascade to the surgical scale.Fig. 13 shows the evolution of the vorticity norms ζ (red) and ζ 2 (black). (The square root of ζ 2 is used to be able to plot all results together.)Again results for the default resolution are compared with those obtained at a coarser resolution (solid and dashed lines, respectively).Note that ζ is the total circulation within the domain.This is captured well at both resolutions, and there is excellent agreement until t = 4.Even ζ 2 is not especially sensitive to resolution.
As expected, this grows to larger values at higher resolution.Yet, although the default resolution simulation has half the grid spacing of the coarser resolution simulation -and thus the capacity to develop twice greater buoyancy gradients in (1) -the maximum r.m.s.vorticity is only approximately 1.27 times larger.It also slightly delayed in time, consistent with the longer time it takes for scales to reach the grid scale.

Conclusions and outlook
We have formulated a method for simulating nearly-inviscid, density-stratified (Boussinesq) flows in complex twodimensional domains.This is done via conformal mapping, in which the physical domain is mapped into a rectangle.This

Fig. 1 .
Fig. 1.Schematic of the general domain under consideration (top) and the corresponding conformal domain (bottom).

Fig. 2 .
Fig. 2.The image of the conformal map grid-lines in the physical domain depicted in Fig.1.Note, tests described in this section typically use ten times as many grid-lines in each direction.The colours along the edges of the domain remain constant between specified vertices.Straight lines connect adjacent vertices so that the domain is actually a polygon.

Fig. 3 .
Fig. 3. Various fields as indicated at the default resolution, and found by inverting the vorticity field ζ in panel (a).

Fig. 4 .
Fig. 4. Errors in streamfunction ψ (left), and the velocity components u (middle) and v (right), together with linear least-square fits in dashed lines.See text for definitions.

Fig. 6 .
Fig. 6.Various fields as indicated at the default resolution, and found by inverting the vorticity field ζ in Fig. 3(a), here in a moving domain.2.0, X S = 0, X = X W R = 1.0 and X = X W S = 0.95.We also take the bottom of the weir at Y = Y W = 0.2 though this has no effect on U S .The values of Y B R and Y B S are determined from the discrete conformal mapping after assigning preliminary values Y B R = −0.5 and Y B S = 0.The mapping rounds off the lower left and lower right corners of the domain.The X

Fig. 7 .
Fig. 7. Errors in streamfunction ψ (left), and the velocity components u (middle) and v (right), together with linear least-square fits in dashed lines.See text for definitions.

Fig. 9 .
Fig. 9. Initial buoyancy field b(X , 0) used in the dam break simulation.The initial vorticity is identically zero, corresponding to a state of rest.
variation of b.Vorticity contours are also equally spaced in ζ , but their spacing changes in response to the changing range in ζ (by contrast the range in b never changes as b is materially conserved).We choose the contour interval in ζ from ζ = ζ 2 50 |ζ |

Fig. 10 .
Fig. 10.Buoyancy b (left column) and vorticity ζ (right column) at times t = 2, 5 and 10 (top, middle and lower rows) in the dam break simulation conducted on a 400 × 200 inversion grid.The initial conditions for b are shown in Fig. 9. Note, the fields here are shown on the ultra-fine grid (6400 × 3200) mapped to the physical domain.

Fig. 11 .
Fig. 11.Kinetic K (red), potential P (blue) and total energy E = K + P (black) as a function of time for the dam break simulation (see text for details).Solid lines correspond to the default resolution (400 × 200), while dashed lines correspond to a coarser resolution (200 × 100).

Fig. 12 .
Fig. 12. Linear and quadratic norms of buoyancy: the linear norm b is shown in red, while the quadratic norm b 2 is shown in black.Solid lines correspond to the default resolution (400 × 200), while dashed lines correspond to a coarser resolution (200 × 100).