The 3D Elliptical Parcel-In-Cell (EPIC) method

We present the three-dimensional version of the Elliptical Parcel-In-Cell (EPIC) method for the simulation of ﬂuid ﬂows and analogous continuum systems. The method represents a ﬂow using a space-ﬁlling set of ellipsoidal parcels, which move, rotate and deform in the ﬂow ﬁeld. Additionally, parcels may carry any number of attributes, such as vorticity, density, temperature, etc, which generally evolve in time on the moving parcels. An underlying grid is used for eﬃciency in computing the velocity ﬁeld from the interpolated vorticity ﬁeld, and in obtaining parcel attribute tendencies. Mixing is enabled by permitting parcels to split when excessively deformed, and by merging very small parcels with the nearest other parcel. Several tests are provided which illustrate the behaviour of the method and demonstrate its effectiveness in modelling complex, buoyancy-driven turbulent ﬂuid ﬂows. The results are compared with a large eddy simulation (LES) and a direct numerical simulation (DNS) model. © 2023


Introduction
In a previous paper, we developed the Elliptical Parcel-In-Cell (EPIC) method for two-dimensional (2D) flows, and applied it to complex density-stratified turbulent flows [1].In particular, we showed that it can faithfully represent turbulent mixing without explicit diffusion or a subgrid mixing parametrisation.Unlike previous Particle/Parcel-in-Cell (PIC) methods, the EPIC method preserves the volume occupied by parcels in each grid cell over long integration times without ever needing to regrid the parcels.Here, we present the corresponding three-dimensional (3D) method and apply it to a variety of test cases, including a moist buoyant bubble studied in [2,3].
PIC methods are hybrid Lagrangian-Eulerian numerical methods that use moving elements to represent certain flow properties (Lagrangian), and use a grid (Eulerian) for operations like constructing the velocity field.The grid also plays a role in computing the tendencies of certain parcel properties, such as vorticity.This hybrid approach combines accurate and efficient aspects of each representation (Lagrangian and Eulerian), and typically benefits from strong numerical stability without explicit stabilisers.A more detailed discussion can be found in [1]; here we focus specifically on the 3D extension of EPIC.
The plan of the paper is as follows.In Section 2, we discuss the flow model and provide an overview of the numerical approach.We then discuss the various new aspects of the 3D EPIC method in detail, some of which differ non-trivially from the 2D method in [1].In Section 3, the method is then applied to several demanding test cases and compared with standard Eulerian (pseudo-spectral and LES) methods.Our conclusions are offered in Section 4.

The flow model
We consider an inviscid incompressible fluid having variable density, under the Boussinesq approximation, and include a simple model for the condensation and evaporation of water vapour.This model was used previously by [2,3] to examine the evolution of an idealised cloud in the Moist Parcel-In-Cell (MPIC) method.The MPIC method uses non-deformable spherical parcels rather than ellipsoids, and differs in a number of other respects as discussed below.The flow model considered here is sufficiently complex to illustrate all of the key features of the EPIC method; we anticipate however that the method may be much more widely applicable, e.g. to geophysical and astrophysical fluid dynamics and beyond (see Section 4).
In a Cartesian coordinate system x = (x, y, z) with z pointing upwards, the prognostic (evolution) equations for velocity u(x, t) = (u, v, w), liquid-water buoyancy b l (x, t) and total water (i.e.water vapour plus liquid water) specific humidity q(x, t) (ignoring precipitation) are Du Dt where t is time, DF /Dt = F t + u • ∇F is the material derivative of any quantity F (subscripts t, x, y and z denote partial differentiation), p is the pressure, ρ 0 is a constant reference density, is the (constant) background rotation vector, and êz is the vertical unit vector.These are supplemented by the relations b = b l + b c max 0, q/q 0 − e −λ c z (5) where λ c is the reciprocal of the condensation scale height, q 0 is the surface saturation humidity, and b c is a constant, in units of buoyancy (acceleration), proportional to the latent heat of condensation L v .Specifically, b c = gL v q 0 c p θ v0 (6) where g is the acceleration due to gravity, c p is the specific heat at constant pressure and θ v0 is a constant reference liquidwater potential temperature.Note, the liquid-water buoyancy is defined by b l = g(θ l − θ v0 )/θ v0 where θ l is the (variable) liquid-water potential temperature.The form of the total buoyancy b in (5) models the latent heating (increase in buoyancy) associated with condensation as a parcel rises above its condensation level z c (q) = λ −1 c log(q 0 /q) (assuming the ground is at z = 0), and likewise evaporative cooling when z falls below z c (q).Further discussion of the physical model can be found in [3].
We can eliminate the pressure p by instead evolving the vorticity ω = ∇ × u, which satisfies Dω Dt = (2 + ω) • ∇u + ∇b × êz .(7) Notably, buoyancy gradients alter only the horizontal vorticity components since ∇b × êz = (b y , −b x , 0).In variables ω, b l and q ≡ q/q 0 , the flow model depends only on the parameters , λ c and b c (besides domain-related parameters).In the simplest case of a dry, non-rotating stratified flow, we have q = 0; then b l = b is the total buoyancy and the variable q is not needed.
We consider a domain which is periodic in both x and y, and bounded by free-slip surfaces in z.In general, x min ≤ x ≤ x max , and similarly for y and z.Along the z boundaries, the normal flow must vanish: if u = (u, v, w) in Cartesian coordinates, then w = 0 at z = z min and z max .There is no restriction on the values of ω, b l and q on these boundariesin particular, the boundaries are not in general shear-free [4].Indeed, variations of buoyancy along the boundaries imply non-zero tendencies for the horizontal components of the vorticity, and hence non-zero values of u z and v z .However, since the boundaries are free-slip in the absence of viscosity, there can be no net generation of vorticity, implying that the domain integral of ω remains constant in time (see [5,6] and references therein).This is trivially satisfied by the vertical vorticity component ζ = v x − u y on account of the periodic boundary conditions in x and in y (in fact the horizontal average of ζ is zero at each height z).For the other two components, constancy of their domain integrals implies that the difference between the values of the horizontal momentum on the two boundaries remains constant.

Vorticity inversion
The use of vorticity as a prognostic variable entails recovering the velocity u by inverting the definition ω = ∇ × u, subject only to homogeneous boundary conditions on the vertical velocity component w.In particular, one cannot generally assume any symmetry properties for the other field components when considering impermeable free-slip boundaries, and this significantly complicates matters.The procedure developed to find u, detailed in [4] and summarised below, is based on a novel field decomposition that permits an arbitrary vertical dependence of fields near the boundaries.
Working in Cartesian coordinates, we take ω = (ξ, η, ζ ) and u = (u, v, w).We consider a horizontally-periodic domain and Fourier transform all field components in x and y: f (x, y, z, t) → f (k, l, z, t) where k and l are the x and y wavenum- bers, respectively.Derivatives in x and y then become multiplications by ik and il in spectral space (here i = √ −1).Since the relation ω = ∇ × u is linear, the component equations can be readily combined to form a Poisson equation for ŵ : subject to homogeneous boundary conditions, ŵ = 0 at z = z min and z max (here a prime denotes z differentiation).From ŵ and the vertical vorticity component ζ , the horizontal velocity components are recovered from Notably, using the field decomposition introduced in [4], the vertical derivative of ŵ is obtained from a wavenumber mul- tiplication of the sine series part followed by a cosine transform, and from the analytic differentiation of the harmonic part.
For (k, l) = (0, 0), i.e. for the part of the flow independent of x and y, instead ŵ = 0 and the horizontal velocity components are found by integrating û = η and v = − ξ .Again, this is done analytically.Full details and tests can be found in [4].

Elliptical parcels: basic representation and dynamics
As in the 2D model [1] we propagate the centre of an elliptical Lagrangian parcel i, located at x i (t), and with velocity u i (t) ≡ u(x i , t) at time t, by the advection equation dx i dt = u i (t) , (10) and evolve its shape by (see [7]) Here the 3 × 3 symmetric positive-definite matrix B (upon dropping the subscript i) defines the shape of the ellipsoidal parcel through x (t)B −1 (t)x(t) = 1 , (12) with denoting a transpose, while S is the local velocity gradient matrix, (Note that x in Eq. ( 12) refers to points on the ellipsoid and is not to be confused with the parcel centre x i .)Thanks to incompressibility and symmetry, we need to store only five elements of B per parcel.This corresponds to three additional parcel attributes compared to the 2D model [1].
Notably, the evolution of materially-conserved quantities like b l and q in Eqs. ( 2) and ( 3) is done almost for free: the corresponding parcel attributes are simply constant -they never need to be updated.Hence, as the parcels move and deform, they exactly preserve all conserved attributes.Only merger (see below) mixes parcel attributes.Critically, the EPIC method is automatically monotone [1].This is in sharp distinction to Eulerian advection methods, where significant effort is often required to preserve monotonicity, besides numerical stability [8].Moreover, the advantages of Lagrangian advection grow as more attributes are added, as would be required in more realistic atmospheric models.
Lagrangian advection is also advantageous for attributes which are not materially conserved yet are advection dominated (e.g. chemical and biological species, etc).In the present model, the vector vorticity is a prime example.As in [2], we evolve the vorticity ω i (t) = (ξ i , η i , ζ i ) on each parcel i using the mathematically-equivalent flux form of Eq. ( 7) (see also [5,6]), (15) dζ i dt = ∇ • (ω a w) , (16) where ω a ≡ ω + 2 is the absolute vorticity, and the right-hand sides are interpolated from the grid to the parcels (see Section 2.5).This appears to be the most stable approach based on extensive numerical tests.We have tried using the usual expression for the vortex stretching term in Eq. ( 7), but this is numerically unstable.The flux form above was previously advocated by [9] in Lagrangian PIC methods; other equivalent formulations discussed in [9] were found to be numerically unstable.
Additionally, to reduce aliasing errors when computing products of gridded quantities, we filter the vorticity ω (before inverting to find the velocity u) using the filter recommended by [10].This filter takes the form exp(−36(k/k max ) 36 ) in each wavenumber direction (the total filter is a product of these one-dimensional filters).Here, it is applied to the decomposed vorticity field (see details in [4]), so is 2D for the surface-associated part and 3D for the remaining interior part.This filter is much less aggressive than the usual de-aliasing filter where the upper third of all modes are set to zero (the latter was used previously in MPIC [2]).Our results however are not sensitive to the choice of the filter.Both filters preserve numerical stability.
An important feature of the algorithm is the preservation of the domain mean vorticity, which is required for mathematical consistency [5,6].We have found that numerical errors result in a drift of ω , which may strongly affect energy conservation -as found in the Beltrami test case in Section 3.3.In EPIC, we remove the domainmean parcel vorticity at the beginning of each time step.Finally, the gridded vorticity field obtained by interpolation from the parcels (see Section 2.5 below) is not exactly solenoidal, as required by the definition ω = ∇ × u.Moreover, evolving all three components of vorticity is redundant since the constraint ∇ • ω = 0 is not guaranteed.A way around this is to evolve the horizontal curl of the horizontal vorticity, χ = η x − ξ y , alongside the vertical component ζ , then recover ξ and η from the relations which become simple algebraic relations (at each z independently) in semi-spectral space, and are easily solved (details omitted).Note that because ζ = v x − u y , both ζ and ζ z have zero horizontal mean values, as does χ .Thus Eq. ( 18) can always be inverted for horizontal wavenumbers k and l not both zero (we do not use the above relations for the horizontal mean parts of ξ and η; these are not modified).However, the dynamical equation satisfied by χ is much more complicated than those for ξ and η.Instead, we continue evolving all three components of vorticity, but correct the horizontal vorticity whenever needed by forming the horizontal curl χ = η x − ξ y then inverting the relations in Eq. ( 18) for ξ and η (excluding the horizontal mean part, k = l = 0).For the derivative ζ z , we use centred differences with linear extrapolation at the boundaries as in [4].

Characteristics of the shape matrix
For various procedures, the EPIC model needs to obtain the elliptical semi-axis lengths, as well as their orientations, from the shape matrix B. Conveniently, the eigenvalues of the shape matrix B are the squared values of the semi-axes a 2 , b 2 and c 2 , where we assume a ≥ b ≥ c without loss of generality [7].The axis orientations are the corresponding normalised eigenvectors â, b and ĉ.
While in principle the eigenvalues a 2 , b 2 and c 2 of the 3 × 3 symmetric matrix B may be found using the exact solution to the cubic first recorded by Cardano [11] (see also [12] for some of the fascinating history), it turns out not to be stable enough numerically due to occasionally severe round-off errors.Typically, we require a large number of eigenvalue/eigenvector solves each time step, and our tests show that some of the matrices having two nearby eigenvalues result in unacceptable errors, far exceeding numerical precision.This problem is well known, and has led to alternative iterative approaches e.g. using the Jacobi algorithm [13,14].Here we use the approach of Scherzinger and Dohrmann [15] that combines the efficiency of analytical methods with the stability and accuracy of numerical methods.

Interpolation
A core feature of PIC methods is the transfer of parcel properties to the underlying grid, and vice-versa.Gridded fields are obtained by summing weighted parcel properties at each grid point, while the velocity of a parcel for example is interpolated from nearby grid points.In EPIC, parcels are ellipsoids of finite volume, and thus the contribution of a given parcel depends on its specific shape.In principle, one must integrate over each parcel, since different parts of it contribute differently to nearby grid points (moreover a parcel may straddle grid cell boundaries).In 2D, this integration can be done efficiently using 2-point Gaussian quadrature, where the points lie along the major axis of the ellipse at half the focal distance from the ellipse centre [1].In 3D, the analogous quadrature uses 4 points positioned within the plane containing the major and middle axes of the ellipsoid [7].This is in fact the simplest quadrature (apart from using a single point); more accurate quadrature formulae use 7, 13 or more points.However, as a compromise between accuracy and efficiency, only 4 points are used here.
The 4 points inside the ellipsoid are called "support points".For parcel i centred at x i , they are located in space at the positions , and ϑ m ≡ mπ /2 − π/4 for m = 1, 2, 3 and 4. The support points are shown in the plane spanned by the eigenvectors âi and bi in the right panel of figure 2 in [7] (note the reference uses the reverse axis order c i ≥ b i ≥ a i ).An alternative 3D view is provided in Fig. 1.The support points always lie inside the ellipsoid.

Interpolation from parcel values to the grid
As in the 2D model, we use linear interpolation in each coordinate direction (here tri-linear interpolation) to obtain a gridded field quantity, say q, from the parcel attributes q i and volumes V i .Here, there are 4 support points per parcel (rather than 2 in the 2D model) and each such point gets an equal weight of 1/4.First, we initialise the gridded volume V = 0 and an auxiliary field m = 0; after interpolation the ratio of these defines q.Second, we sum over all parcels and support points, interpolate, and accumulate V and m.For each support point x m i , we locate the grid box containing it -this is easily done on a regular grid.Let x denote the corner of this grid box having the smallest x, ȳ and z.The interpolation weight between x m i and x is The analogous weights at the other 7 corners of the grid cell are found by replacing x with x + x, or ȳ with ȳ + y, or z with z + z above.These weights are multiplied by the parcel volume V i and added to V at each corner of the grid cell.Likewise the weights are multiplied by q i V i and added to m.Once we have summed over all parcels and support points, the gridded quantity q is found simply by the ratio m/ V at all grid points x.Note that the same weights φ m i can be used to generate all fields at once for efficiency.
As in 2D, we must double the gridded values at the boundaries i z = 0 and n z to ensure the total gridded volume is conserved and equal to the domain volume.For other gridded quantities this factor of two cancels out, as required, when taking the ratio m/ V .A downside to this interpolation is that fields having a linear variation in z near the boundaries are in general poorly represented, unless the parcel attributes are chosen carefully.But we can only choose the parcel attributes at the initial time to match input gridded values (see Section 2.10 below), and have no control at later times.All of our attempts to use higher-order extrapolation led to numerical instability (divergence of kinetic energy) in the Beltrami test case discussed in Section 3.3.This boundary interpolation is over-ridden for the vertical vorticity component ζ when it is initially zero and when it cannot become non-zero theoretically (i.e. when the vertical component of the background rotation vector is zero).In this case, we set the gridded values of ζ as well as its gridded tendency Dζ /Dt at i z = 0 and n z to zero.

Interpolation from the grid to parcel values
The reverse operation, interpolating grid values q to each parcel support point, is done using the same weights in Eq. ( 20), now with contributions from all 8 corners of the grid box containing the support point.Then the values accumulated at the 4 support points per parcel are summed together to define the parcel attribute q i , usually a parcel tendency or velocity.

Parcel splitting
A parcel is split if either the parcel aspect ratio λ (the ratio of the largest to smallest semi-axis length) exceeds a user-specified parameter λ max , or the maximum parcel semi-axis length a (suppressing the index i) exceeds a max = 3 √ 3/4π min{ x, y, z}.In 3D, the aspect ratio λ = a/c, where the semi-axis lengths satisfy a ≥ b ≥ c.The value of a max is chosen to limit the maximum parcel volume V to the grid cell volume V cell = x y z.This is a change from the 2D model (see [1]) which simply limited the parcel volume directly.
After splitting, the shape matrices of the two new parcels are given by just as in the 2D model.Only the semi-major axis length a of the original parcel is divided in two; the other axis lengths remain the same.Each new parcel has half the volume of the original parcel.The only difference with the 2D model is that the new parcel centres are separated by a distance 2h = √ 3/5 a rather than √ 3/4 a in 2D.The new parcel centres lie at x ± h â (details may be found in Appendix A).Note, parcel splitting preserves the zeroth, first and second-order spatial moments.

Parcel merging
The criterion used for parcel merging is identical to that used in the 2D model.A parcel m is marked available for where Ṽ min is another user-controlled parameter (all numerical parameters may be found in Table 1 below or in [1]).Moreover, the procedure to find and merge groups of merging parcels remains the same as that described in appendix D of [1].
Merger preserves the total volume (zeroth moment) , approximately, the second-order spatial moments, where i ranges over all parcels in the group of merging parcels.The only difference in the entire merging procedure stems from the fact that the parcels are ellipsoids rather than ellipses.This means that the (tentative) shape matrix of a merged parcel must now be computed from where xi ≡ x i − x (these are the only unique components since B is symmetric).The final shape matrix is found from which ensures that the volume of the merged parcel equals the sum of the volumes of the merging parcels (det(B) = (abc) 2   is proportional to the squared parcel volume).

Enforcing boundary conditions on parcel operations
In the 3D model, the domain is periodic in x and in y, and periodicity is handled in precisely the same manner as in the 2D model [1] (where only x is periodic).At both free-slip z boundaries, the grid extends beyond the domain by one halo grid layer (i z = −1 and i z = n z + 1) to accommodate parcel centres and support points that may appear outside during an advection step.These boundaries are also handled the same way apart from one difference: the halo grid points are filled by extrapolation.In particular, all velocity and strain components are extrapolated, rather than mirrored assuming symmetry or anti-symmetry as in the 2D model.Tests using extrapolation in the 2D model (see [1]) in fact show generally negligible differences.This is because the halo grid points are both rarely needed and receive a small weight when they are.Nonetheless, in the 3D model, we have opted for extrapolation based on the exact test flow discussed in Section 3.3, for which neither the horizontal velocity components nor their z derivatives vanish at the boundaries.As the vertical velocity component w = 0 on the boundaries, extrapolation is equivalent to assuming anti-symmetry across the boundaries.

Incompressibility
As in the 2D model, at each time step we correct for errors in the gridded volume V arising from the parcel motion.Any finite number of parcels cannot in general ensure V = V cell at all times, and a novel feature of the EPIC method is that these errors are corrected every time step by a minor adjustment of the parcel centres (see further discussion in [1]).This adjustment is done in two complementary ways to accelerate the convergence of V to V cell .First, a displacement field ∇φ is found by solving Poisson's equation ∇ 2 φ = δ ≡ ( V − V cell )/V cell , subject to the boundary condition φ z = 0 at z = z min and z max to prevent parcels from leaving the domain.In semi-spectral space, after a horizontal Fourier transform, each Fourier This is solved numerically using centred differences resulting in a tridiagonal problem analogous to that described in the appendix of [1], only not using compact differencing.
Once the parcel centres are adjusted by ∇φ, the gridded volume V is recomputed for the next stage of the parcel cor- rection algorithm.In this stage, parcels are displaced in a way which depends on their location within the grid cell and the discrepancies between V and V cell at the eight corners of the grid cell.Each parcel coordinate is displaced independently; so for x i , the displacement is C x xi (1 − xi ) x where xi = (x i − x)/ x varies from 0 to 1, and C x is a dimensionless pre-factor dependent on V − V cell at the corners of the grid cell.For a cell edge which is oriented in the x-direction, between grid points xn and xn+1 , C x is given by The default value of β is set to 1.8, as in [1].For a point in the cell interior, linear interpolation is used to determine C x , again following [1].The magnitudes of the pre-factors C x , C y and C z are limited to C max , a user-specified numerical parameter with default value 0.5.These two stages are repeated once to keep volume errors small, typically below 10 −3 V cell .This is demonstrated below (see Fig. 8 in Section 3.1).

Initialisation
A simulation normally begins by reading in field data from a grid for each prognostic quantity in the model, here b l , q and ω.These fields are then used to define the corresponding parcel attributes by interpolation from the grid as explained in Section 2.5.2.By default, we place a regular 2 × 2 × 2 array of parcels in each grid cell, each having a volume of V cell /8.As the parcels are typically spherical, we simply interpolate the gridded values to the parcel centres to define the parcel attributes.Note that the initialisation does not guarantee that the reverse interpolation, from parcel attributes to gridded values in Section 2.5.1, matches the given field values.An algorithm designed to match field values was proposed in section 3.1 of [1], but in applications here we have found that it can lead to overshoots and undershoots of parcel attributes.Such overshoots and undershoots have been found to degrade energy conservation in the test flows examined below.Exceptionally, we may bypass this algorithm and directly assign parcel attributes from analytical functions as in [2,3].

Time stepping
Time stepping is carried out by a low-storage 4th-order Runge-Kutta method, exactly as described in [1].The time step t is adapted at the beginning of each time step to ensure that the maximum (generalised) buoyancy frequency N max and maximum rate of deformation γ max are both resolved in time.Here N ≡ ∇b while γ is the maximum eigenvalue of the symmetrised velocity gradient matrix (S + S )/2, as used in the 2D method [1].Note that N and γ are evaluated at each grid point, and the maximum is taken over all grid points.

Table 1
Numerical parameters and their default values.The default value of the initial number of parcels per cell and the minimum parcel volume fraction are different compared to the 2D method [1].Also, the 3D model has a slightly different parcel splitting criterion which does not require the maximum parcel volume fraction parameter anymore, see Section 2.6 for further details.

Numerical parameter Value
Initial number of parcels per cell 8 Maximum aspect ratio, λ max The time step is chosen as where α is a prescribed dimensionless numerical constant.The recommended default, α = 0.2, is the result of extensive tests carried out in [1].Tests with the 3D method continue to show that α = 0.2 is adequate -a compromise between accuracy and efficiency.

Results
We next illustrate the performance of the EPIC method in four distinct test cases.We assess accuracy by the rate of convergence of early-time solutions, gridded-volume errors in time, probability density functions of vorticity, etc.We also compare with the Moist Parcel-In-Cell (MPIC) method previously developed to study moist convection [2,3], as well as with both a standard large eddy simulation (LES) method and a pseudo-spectral method (PS3D) [4,16] to exhibit some of the distinct advantages of EPIC.For reference, we list all of the numerical parameters used in Table 1.In Appendix E, the cost and accuracy of the model is examined with respect to the minimum parcel volume fraction and the maximum parcel aspect ratio using the rotating Rayleigh-Taylor test case.For PS3D, we use α = 0.1 for time stepping unless otherwise indicated.
In several of the test cases, we assess convergence through the evolution of the kinetic and potential energy components, K and P, as well as the enstrophy ϒ.These are defined per unit volume by where V dom = L x L y L z is the volume of the domain.Since the potential energy P may be dominated by the background hydrostatic equilibrium state (see in particular the internal wave test case in Section 3.2), it is more useful here to consider the available potential energy A = P − P ref , where P ref is the constant potential energy of the buoyancy (i.e.density) field re-arranged into hydrostatic equilibrum (see [17][18][19] and §2.9 in [1]).For small departures from hydrostatic equilibrium, numerical inaccuracies in the discrete calculation of the reference state can swamp the actual available potential energy, a known issue [18,19].
Where possible, we make use of the known hydrostatic reference state, b ref (z), to accurately compute A from the form proposed in [20].There, the available potential energy density, a(b, z), is defined by
In this test case, the analytical form of the available potential energy density a(b, z) can be found from Eq. ( 30 Since the initial vertical buoyancy profile has density increasing with height z, the flow overturns in an attempt to restore hydrostatic equilibrium.During this process the fluid undergoes strong mixing which decreases the buoyancy magnitude and leads to a sharp decay in kinetic energy (see below).Fig. 3 illustrates y = 0 cross sections of the gridded buoyancy field b computed by EPIC at several representative times.For comparison, Fig. 4 shows results at approximately the same times, but now from a 128 3 simulation using the Eulerian pseudo-spectral method 'PS3D' described in [4].Notably, the detailed flow structure obtained with EPIC is closely similar to that obtained with PS3D at early times.During the turbulent phase when strong mixing occurs, the simulation results diverge as expected in part due to loss of predictability, but nonetheless remain qualitatively similar.As hydrostatic equilibrium is approached at later times, the results more closely resemble each other again.
We next examine the energy and enstrophy evolution, comparing EPIC and PS3D both for a grid resolution of 128 3 .
Enstrophy ϒ peaks at a substantially larger value in EPIC than in PS3D -see the lower right panel of Fig. 5.The enstrophy growth corresponds to the destabilisation of the unstable base flow.Otherwise, the kinetic energy K (shown in the upper left panel of the same figure) exhibits closely similar behaviour, with the steep decay occurring at nearly the same time in both methods (slightly later in EPIC).The initial decay in available potential energy A closely compensates the growth in K, as required for total energy conservation.In Table 2 we list the energy loss at time t = 4 relative to the initial total energy at t = 0 as a function of the grid resolution.For EPIC, the energy loss reduces by at least √ 8 when doubling the grid resolution in all dimensions.
Notably, the PS3D simulation exhibits buoyancy overshoots and undershoots as found previously in 2D simulations [1], while EPIC by construction exhibits monotonically increasing and decreasing buoyancy extreme values b min (t) and b max (t), respectively -see Fig. 6.The non-physical undershoots and overshoots in PS3D are caused by hyperdiffusion, and are most prominent near the boundaries, clearly visible in Fig. 4.This feature is especially undesirable in situations where it is important to preserve monotonicity, e.g. when a field like specific humidity q must remain positive definite (see Section 3.4 below).While the pseudo-spectral method normally has excellent energy conservation properties, these may be outweighed by the method's inability to preserve monotonicity.Alternatively, one can use an explicitly monotonic advection scheme as in the large eddy simulation (LES) model featured in Section 3.4, but this comes with a substantial increase in numerical  cost.The vertical profile of mean buoyancy (mean with respect to the horizontal directions x and y) in Fig. 7 shows good agreement between EPIC and PS3D at early times, but at later times (t > 5) PS3D exhibits larger -and unphysical -mean buoyancy values at the bottom and top boundaries compared to EPIC, see in particular times t = 6 and 10.
Finally, we verify that EPIC closely maintains incompressibility (local volume conservation) in this demanding test case.As explained in Section 2.9, we apply two parcel correction methods after every time step to maintain a homogeneous and space-filling parcel distribution inside the domain.To demonstrate the effectiveness of these parcel correction methods, we evaluate the (relative) r.m.s.volume error as a function of time.As shown in Fig. 8, EPIC keeps the Note that we used the parcels to evaluate these quantities in EPIC.

Table 2
Relative total energy loss for the Rayleigh-Taylor test case between the initial time t = 0 and time t = 4 as a function of the grid resolution.At this early time, energy is still well conserved.

Number of grid cells
Energy loss (%) EPIC PS3D 32 3  1.45 1.55 48 3  0.81 0.69 64 3  0.51 0.43 96 3  0.25 0.24 128 3  0.15 0.18 384 3  0.02 relative error below 1.5 • 10 −3 throughout the simulation, despite the strong mixing that takes place.This error is similar to what was found previously in [1] for the 2D method.Note that mixing in EPIC is fully represented by explicit splitting and merging of parcels, and no additional parameters (such as a Smagorinsky constant or turbulent Prandtl number) are needed.With the exception of the initial parcel number per cell, these parameters are the same as used in the 2D method [1], where a thorough parameter-sensitivity analysis was carried out.In Appendix C, a stochastic model is used to argue that the choice of a max is not important for a well-developed flow.

Internal wave test case
In this test case we consider a resting basic state with linear stratification b = N 2 z in a uniformly rotating flow with The linearised equations for vorticity and buoyancy are We now seek eigensolutions of the form Fig. 6.Evolution of the gridded buoyancy and parcel buoyancy extrema ( bmin , bmax ) and (b min , b max ), respectively.The theoretical bounds [−1, 1] are indicated by the dashed lines.PS3D exhibits large undershoots and overshoots, greatly exceeding these bounds.Fig. 7. Gridded mean buoyancy vertical profile at various times for the rotating Rayleigh-Taylor test case.Fig. 8. Gridded relative r.m.s.volume error for the rotating Rayleigh-Taylor test case with different grid resolutions.EPIC manages to preserve incompressibility very well -the error remains smaller than 1.5 • 10 −3 -despite the strongly turbulent overturning flow.Note that the simulation having 384 3 grid cells does not use the default value of the minimum parcel volume fraction (cf.Table 1), but Ṽ min = 1/40.where ϕ = kx + ly − σ t.These forms are consistent with w = 0 at z = ±L z /2, as well as the definition of vorticity and the incompressibility relations for velocity and vorticity, so long as m = (2 j − 1)π /L z with j a positive integer.
Inserting these forms into Eqs.( 35) to (38) and performing much algebra leads to the dispersion relation Moreover, if we start with w = ŵ cos mz cos(kx + ly) and take ŵ to be real, then the other fields for any t ≥ 0 are We next use these solutions to test the numerical method.We initialise the flow with the above vorticity and buoyancy fields at t = 0, taking a small value of ŵ, here 10 −3 , so that the linearised approximation is valid (nonlinear effects are expected to arise at O( ŵ2 )).These fields are then used to obtain the initial parcel attributes as described in Section 2.10.
In the test case illustrated, we have taken k = l = 1/2 and m = 1 in a domain of extent 4π × 4π × π centred on the origin, and have used an isotropic grid resolution of 48 × 48 × We have also taken N = 2 and f = 1, giving an eigenfrequency of σ = √ 2. The wave period is 2π / √ 2 ≈ 4.44.
In linear theory, the flow evolves by uniformly translating in x and y at constant speeds σ /k and σ /l, both equal to 2 √ 2 in the example illustrated.Simulations using both EPIC and a pseudo-spectral method, PS3D, originally developed in [4] and extended to include buoyancy effects (version 0.0.13 [16]), closely approximate the exact linear solution despite the low resolution employed.Note that for PS3D we evolve the buoyancy perturbation b instead of the total buoyancy b = b + b , which here is dominated by b = N 2 z.The perturbation buoyancy evolves according to b t + u • ∇b = −N 2 w.Since only the horizontal gradient of b is needed for the vorticity tendency, in EPIC we recover b on the grid by interpolating the parcel values of b i − N 2 z i .When the total buoyancy is required for output purposes, we simply add N 2 z to b .This is done to avoid boundary interpolation errors associated with the dominant mean part of b.
A convergence study of the kinetic energy K, available potential energy A, total energy T and enstrophy ϒ evolution is shown in Fig. 9 for both PS3D and EPIC.We make use of the known reference state, b ref (z) = N 2 z, to accurately compute A from the form proposed in [20], i.e. as an integral over the available potential energy density a, here given by For the linearised solution for b above, one can show that A = N 2 ŵ2 /(8σ 2 ), or 2.5 × 10 −7 for the parameters chosen.Moreover, the ratio of kinetic to available potential energy is ), or 2 for the parameters chosen.In EPIC, K, A and ϒ are calculated from sums over parcels.In PS3D, A uses a limited form of b, namely min(b max , max(b min , b)), restricted to lie between the theoretical minimum and maximum buoyancy values b min and b max , here −2π and +2π respectively.This is not necessary in this test case, but is essential in other cases where PS3D produces unphysical overshoots and undershoots, see in particular Section 3.1.EPIC, by design, does not generate such errors.
As indicated by the dashed lines in Fig. 9, the convergence rate is quadratic in K, A and ϒ for PS3D and EPIC.The figure shows the relative error with respect to the time-averaged quantities.For the theoretically-conserved quantity T , the total energy, cubic convergence is found, presumably because the other quantities are only conserved by the linear dynamics.In Table 3 we also list the relative total energy loss in ‰ between the initial time t = 0 and the final time t = 4π / √ 2 as a function of the grid resolution.Again, we observe an approximate cubic convergence for EPIC and PS3D.Fig. 10 shows the time evolution of the relative errors in K, A and ϒ, for various grid resolutions, now with respect to the exact theoretical values.For EPIC, these errors feature a wave period that closely matches the internal wave period, 2π / √ 2. For PS3D on the other hand, the frequency is twice as high.Notably, in this test case PS3D requires a time step 20 times smaller (i.e.α = 0.01) than EPIC to achieve the level of accuracy shown.Otherwise the relative errors increase with a linear trend and do not converge with spatial resolution.The errors reduce by approximately a factor 4, 4 and 3.6 for K, A and ϒ, respectively, when doubling the grid resolution in EPIC.The higher errors in EPIC are due to the parcel initialisation which cannot exactly represent the initial input fields.Regarding computational cost, the runtime for both EPIC and PS3D grows as the cube of the vertical grid resolution as shown in Fig. 11.Note that we can only compare the growth rate, but not individual execution times since EPIC was run with 8 MPI (Message Passing Interface) ranks (distributed memory parallelism) and PS3D with 8 OpenMP threads (shared memory parallelism).
Another way to quantify errors is provided in Fig. 12, which shows the vertical vorticity difference ζ − ζ between the numerical and exact (linear) solution at several times over two wave periods (using 256 2 × 64 grid cells).Here this  difference is scaled by ζ rms , the r.m.s.value of ζ .The larger error in EPIC (at maximum 8.5 times that of PS3D) is partly due to the parcel initialisation, specifically the tri-linear interpolation, which cannot perfectly match the initial fields.Note that the r.m.s.error (lower panel of Fig. 12) remains less than 1.1% of ω rms .The growing error in EPIC may be a phase error caused by the imperfect initialisation.This error decreases quadratically with spatial resolution, and is nearly absent in PS3D, which is why the error is so much smaller in those simulations.

Beltrami flow test case
Next we consider the same Beltrami flow studied in [4].In the domain [−π /2, π/2] 3 , the unperturbed velocity field is while the corresponding unperturbed vorticity is ω = 3u.This theoretically steady flow is unstable to small perturbations and breaks down into strongly anisotropic turbulence, with maximum intensity near the boundaries [4].To force a similar initial instability development for EPIC and PS3D, we explicitly perturb the initial horizontal vorticity components ξ and η by ξ = A ξ cos 2 y cos z and where we have taken A ξ = 1/5 and A η = 1/10 in the results below.This perturbation ensures ∇ • ω = 0. Fig. 13.R.m.s.error in the initial fields compared to the (exact) input fields.For PS3D (rightmost panel), the error in the horizontal vorticity components is due to the solenoidal correction.The vertical vorticity component is not changed by the solenoidal correction and the error remains at machine precision.For EPIC (left and centre panel), this error is mainly due to the parcel-to-grid interpolation.The centre panel excludes the lower and upper surfaces in the r.m.s.error calculation.The error in the horizontal vorticity components is then lower than the error in the vertical component.

Table 4
Comparison of the enstrophy peak value and the slope of the kinetic energy decay between EPIC and PS3D.EPIC is split into parcel (p) and gridded field (f) evaluation.The slope is obtained by a least-squares fit from the data which is 25% around the middle value of kinetic energy.

Grid
Enstrophy  [4,16], uses a mixed spectral space representation of the prognostic fields in order to represent shear (non-zero horizontal vorticity) on the free-slip boundaries.The cascade of enstrophy to small scales is diffused in 2D only (by hyperdiffusion).In EPIC, there is no need for hyperdiffusion as the mixing associated with parcel splitting and merging has a weak diffusive effect.Moreover, advection in EPIC is monotone and all fields remain within their theoretical bounds even during mixing.
In both methods, a solenoidal correction is required each time step to maintain ∇ • ω = 0 to machine precision.This correction only modifies (weakly) the horizontal vorticity components [see 4, for details], but it means that even the initial fields (after correction) differ from the input fields.In EPIC, a further error in initialisation occurs due to assigning parcel attributes directly from interpolated input fields.After a parcel-to-grid interpolation, there is generally a slight discrepancy between the initial and input fields, which is most pronounced at the z boundaries (the interpolation is only first-order accurate there, see Section 2.10).This error in initialisation is quantified as a function of resolution for both numerical methods in Fig. 13.PS3D makes essentially no error in ζ (only round-off); otherwise both methods exhibit errors scaling like the square of the grid spacing, as expected.PS3D is not 'spectrally accurate' due to the need to do some grid operations in z by finite differences in order to handle shear at the z boundaries [4].The left and middle panels show that the boundary interpolation in EPIC is significantly less accurate (first-order) than the interior interpolation (second-order), degrading the overall accuracy in capturing the input fields.In summary, the two methods start from slightly different initial fields, and one may anticipate subsequent differences in the flow evolution.
Fig. 14 compares the evolution of kinetic energy and enstrophy between the methods and across a range of resolutions.There is broad agreement in the timing of the instability (i.e. the peak in enstrophy) and in the overall decay in kinetic energy.The PS3D simulations all start from nearly the same value of kinetic energy for all grid resolutions, while the EPIC simulations exhibit lower kinetic energy at coarse resolution -a consequence of the parcel initialisation.Nevertheless, qualitatively the results compare well, with the most significant difference being the magnitude of enstrophy, which is much higher in EPIC than in PS3D (see Table 4 for a quantitative comparison).The middle panel shows the effect of interpolating parcel values to compute these diagnostics from the gridded fields.This interpolation averages the locally intense subgrid vorticity values, resulting in much weaker enstrophy, comparable to that found in PS3D at the same grid resolution.The parcel-based diagnostics in the left panel, however, are appropriate, since they measure the additional degrees of freedom afforded by the parcel representation of subgrid scales.
Enstrophy growth occurs as the flow destabilises, characterised by vorticity stretching and intensification.Cross sections of the vorticity magnitude |ω| in the horizontal mid-plane z = 0 are shown in Fig. 15 at three characteristic times around the peak in enstrophy, and for both the EPIC and PS3D methods.Both methods use 128 3 grid cells.The computed fields are closely similar despite the differences in the initial fields.In EPIC, however, the vorticity field is both smaller scale and significantly more intense, consistent with the much higher enstrophy growth seen previously in Fig. 14.A second view is provided in Fig. 16, showing |ω| in the xz-plane at y = 0 at the same times as Fig. 15.EPIC and PS3D again display a similar instability development, breaking down into disorganised turbulence around the time of peak enstrophy.Notably, PS3D exhibits significant vertical Gibbs fringing due to the spectral filter.A final view is provided in Fig. 17 at the bottom surface of the domain (z = −π /2).The structure of the vorticity field is broader scale but becomes organised into intense fronts (shear line) where the horizontal velocity rapidly switches direction.The PS3D simulation here has more intense boundary vorticity (only at these early to intermediate times), yet the EPIC and PS3D fields are broadly similar in structure.At later times (not shown), the bottom and top surfaces contain the most intense vorticity, mainly because the boundaries restrict vortex line bending and hence dissipation as discussed in [4].

Moist bubble test case
In the final test case, we compare EPIC with MPIC [2], a PIC method using non-deformable parcels, and with MONC [24], a standard LES model.We examine the evolution of a rising moist buoyant bubble, whose setup is detailed in section 2.4 of [3].In MONC, we use a Smagorinsky subgrid scheme, and the model configuration is the same as described in [3].In EPIC, in contrast to those earlier studies, we use a dimensional (rather than dimensionless) setup.However, for ease of comparison with [2,3], we present non-dimensional results when we compare the models below.
We consider a non-rotating flow ( = 0) initially at rest within a cubic domain, with x, y and z ranging from 0 to 6280 m.The background state consists of a lower neutrally-stratified layer extending from z = 0 to z = z b , where the liquidwater buoyancy b l = 0, and a linearly-stratified region above, where b l = N 2 (z − z b ) (for constant buoyancy frequency N).In the neutral lower layer, the specific humidity is uniform, q = q n , while in the stratified region it decays exponentially, q = q n exp(−λ c (z − z b )).Recall, λ c is the inverse condensation scale height, see Eq. (5).
We next centre a moist, buoyant spherical bubble of radius R = 800 m at (x • , y • , z • ) = (3140, 3140, 800) m.Within the bubble, we consider a weakly-varying liquid-water buoyancy and specific humidity distribution which smoothly transitions to the environmental values at the edge of the bubble (the smooth edge differs from earlier work using the MPIC model, but makes convergence easier to assess).Specifically, for x < R, with x ≡ x − x • , the liquid-water buoyancy and specific humidity distributions are given by where the edge-smoothing function S(h) takes the form Here, f s = 0.8 gives the fraction of the radius of the bubble that is not affected by the edge smoothing.Note that S(h) is continuous at h = 0 and 1, and moreover the first and second derivatives of S(h) vanish there (also S (h) = −30h 2 (h − 1) 2 , implying S is a monotonically decreasing function of h).
All of the required parameters for this test case and for our model equations, namely Eqs. ( 2), (3), ( 5) and (7), are listed in Table 5.These directly correspond to the values used in [3] if one takes a characteristic time scale of 100 s and a characteristic length scale of 1000 m.We use these scales, and a surface humidity q 0 = 0.015 kg/kg, to diagnose the flow in terms of non-dimensional quantities, so that figures can be directly compared to those in [2,3].For MPIC, the parcel attributes are initialised directly from the analytical functions above, following [3].In EPIC simulations, we found the differences between both methods of initialisation to be minor.Fig. 18 shows a cross section of the specific humidity field, q, at non-dimensional time t = 6, where the left-most panels visualise the individual parcels (only a small part of the domain is shown).The total parcel volume in both simulations is initialised to equal the grid volume, yet notably there are both spaces between parcels and regions of overlap in both simulations.The MPIC simulation has larger empty regions at the edge of the updraft, whereas in EPIC the parcels remain fairly homogeneously distributed.Our earlier work on two-dimensional flows [1] found that both the ellipsoidal shape of the parcels and the use of corrections to the parcel position to enforce incompressibility contribute to the more homogeneous distribution in EPIC.
In general, EPIC resolves the structure of the updraft in significantly more detail.The other panels in the figure show a visualisation using a super-gaussian kernel (as in [1]).The kernel is applied inside an ellipsoid which is scaled by a factor three or five compared to the corresponding ellipsoidal parcel.For EPIC, there are no significant contributions to the grid values beyond three times the parcel edge position.For MPIC, on the other hand, there is a small part of the domain which receives no contributions when using only three times the parcel edge position, and the edge of the updraft looks much more blurred when using a kernel up to five times the parcel radius.We will nevertheless use the more compact kernel in the remainder of this section.
Fig. 19 shows the same cross section of specific humidity q at different grid resolutions.We also show a MONC simulation with the same number of grid points.Note that for large simulations, EPIC simulations are about three times more expensive than MPIC ones (cf.Table 6), whereas MONC simulations have similar cost to MPIC.However, the time step used here in MONC is about 8 times smaller in the 512 3 EPIC and MPIC simulations, and is a relatively conservative estimate (EPIC and MPIC use a fourth-order Runge-Kutta method, whereas MONC uses a leapfrog approach).The relative cost of EPIC and MPIC may increase in longer simulations, where the parcel size distribution will approach equilibrium in the entire domain.As in [3], MPIC and EPIC can resolve the circulation using less than half the number of grid points in each direction.Compared to Fig. 18.Cross sections of the specific humidity field q at t = 6 in the xz-plane (zoom) in EPIC and MPIC simulations with 32 3 grid points, rendered as parcels with sharp edges, or with the algorithm used in [1] using two different cutoff ellipsoids (see main text).MPIC, EPIC shows an even sharper contrast between regions of high and low humidity at small scales.Diffusive effects are clearly evident at resolutions less than 512 3 in MONC.
The (normalised) histogram of liquid water specific humidity, namely q l = max(0, q − q 0 exp(−λ c z)), at t = 6 in the three models is shown in Fig. 20.We used bilinear interpolation to obtain better sampling of the MONC simulations with 32 3  and 64 3 points, increasing the number of sampling points in each direction by a factor of 5 and 3 respectively, so all MONC fields are sampled using at least 128 3 points (this is a relatively simple approach; more advanced methods for obtaining histograms are discussed in [25]).For comparison, the results from all 512 3 simulations are shown in the last panel.Both EPIC and MPIC show a peak corresponding to parcels with high liquid water content q l at low resolution, although the height of this peak is lower in EPIC as compared to MPIC.This indicates there is still scope for improving the representation of small-scale mixing in the lower-resolution EPIC simulations.MONC has many grid points with small amounts of liquid water at low resolution, which are absent in EPIC and in high-resolution MONC reference simulations.This is an advantage of the Lagrangian representation of mixing in EPIC.At the highest resolution, EPIC has more parcels with high q l compared to MONC, but the distributions are beginning to look similar, despite the variations in EPIC being at much smaller scales, as seen in the cross-sections.MONC under-represents high q l , as can be seen by the distribution extending outwards with increasing resolution.The significantly greater diffusion in MONC compared to EPIC is visible in the lower right panel, where the MONC distribution is shifted to lower q l .
A quantitative assessment of the convergence of the liquid water specific humidity histogram at t = 6 is presented in Fig. 21.Here, the root-mean-square error over all bins in the liquid water specific humidity distribution is plotted (except for the bin which contains the 0 value, and using linearly spaced bins between q l = 0 and 0.1).In each of the panels of Fig. 21, a different high-resolution model is used as the reference model when calculating the root-mean-square error, so that both convergence of a single model and convergence between models can be assessed.Each model gradually converges to its own liquid water specific humidity histogram at high resolution.Surprisingly, at the highest resolution, EPIC and MONC have smaller differences using this metric than EPIC and MPIC, even though EPIC and MPIC are both PIC models.This indicates that the mixing occurring in EPIC is substantially different from that occurring in MPIC.The difference between the models at the highest resolution is generally smaller than the difference between a high-resolution and a low-resolution simulation of the same model, although there is still a difference between MPIC and MONC at the highest resolution.
The histogram of vorticity magnitude, |ω|, at t = 6 is plotted in Fig. 22.The vorticity in MONC is computed from gradients in the velocity field and no subsampling is applied to the MONC data.The maximum values of vorticity magnitude in MONC are significantly smaller than the maximum parcel-based values in MPIC and EPIC (the histograms are plotted on a log-log scale).MONC and MPIC have more very small values (< 10 −5 ) of |ω| than EPIC.For |ω| between 10 −3 and 1 (these values largely correspond to regions outside the cloud), the histograms are in good agreement, both between models and across resolutions.For high |ω|, the significantly greater diffusion in MONC is visible by the shift in the distribution to lower values, creating a bulge at moderate values of |ω|.Sampling on a grid also mixes weak vorticity more effectively, giving the impression that there are many weak vorticity values.In EPIC, we sample the PDF (probability density function) using parcels, and this avoids the averaging associated with grid sampling.The spike in the vorticity histogram of EPIC at low |ω| is due to the mean vorticity correction, cf.Eq. ( 17).This correction ensures that there is no global vorticity production, i.e. the mean vorticity remains constant.MPIC and MONC do not apply this correction, and therefore their histograms do not exhibit this spike.Nor does EPIC without this correction, and the histogram differs from that plotted in Fig. 22 only for extremely weak vorticity, |ω| < 10 −5 (not shown).These values have a negligible influence on the dynamics of the flow.
To emphasise the wide range of scales captured by EPIC, in Fig. 23 we show a cross section of the specific humidity field q in the xz-plane at t = 6 including the intersected ellipses of the corresponding ellipsoids.The figure also contains four cascaded zooms having widths of 64, 32, 16 and 8 grid cells.The individual ellipsoids (or their elliptical cross sections) are  visible in the finest zooms.Many ellipsoids fill each grid cell, providing an explicit subgrid model.This ability to resolve at, and below, the grid scale enables EPIC to employ modest grid resolution to capture highly-turbulent flow structures with minimal diffusion -and this process is entirely controlled by parcel splitting and merging, not by any explicit numerical diffusion.

Discussion
In this paper, we have extended the Elliptical Parcel-In-Cell (EPIC) method, first developed for two-dimensional flows in [1], to three dimensions.This offers an accurate new, predominantly Lagrangian, approach to simulating very high Reynolds number turbulent, rotating stratified flows -including simple thermodynamics and moisture.Like in two dimensions, the ellipsoidal parcels used to represent key flow attributes (e.g.buoyancy and vorticity) generally deform in the local velocity  strain, and furthermore may split and merge at sub-grid scales, thereby providing a natural mixing mechanism.The threedimensional extension of the associated parcel operations is straightforward and involves only minor code changes.Also, the ellipsoidal parcel shape requires just three additional attributes per parcel.Interpolation from parcel attributes to gridded values (and vice-versa) uses four support points instead of the two needed in the two-dimensional model.Determining the location of these support points requires the solution of a 3 × 3 eigensystem per parcel, in effect a cubic equation for the squared semi-axis lengths.The famous analytical solution of Cardano (1545) [11] unfortunately does not provide sufficient accuracy for a typical distribution of parcel shapes.Instead, we use the algorithm by Scherzinger and Dohrmann [15] which is comparable in speed to the Cardano method and yet meets our accuracy and robustness criteria.
Here we have benchmarked the method against four diverse test cases that exhibit relevant characteristics of atmospheric and oceanic flows.We have also compared the results with a LES (Large Eddy Simulation), a DNS (Direct Numerical Simulation) and another PIC (Parcel-In-Cell) model.For comparable grid resolutions, EPIC exhibits closely similar conserva-tion properties, i.e. comparable decay of energy, but a significant improvement in modelling vorticity.As measured by the mean-square vorticity (enstrophy), EPIC permits much higher values (roughly two times larger), owing to its less aggressive treatment of mixing (and lack of explicit numerical diffusion).Another benefit of EPIC is that it preserves monotonicity of tracers whose extrema are bounded theoretically.By construction, the parcel splitting and merging responsible for mixing in EPIC either preserve attributes or create new values between existing ones.This is a particularly important property of a numerical method when tackling problems involving fields that must remain positive-definite (e.g.moisture here, but more generally chemical and biological tracers).
A major challenge in the development of the three-dimensional version of the EPIC model was the accurate assignment of z boundary grid point values during the parcel-to-grid interpolation.As in the two-dimensional model, the boundary values must be doubled after interpolating the parcel properties onto the grid to ensure that the total parcel volume, and the volume-weighted sum of each parcel attribute, are equal to their gridded counterparts.This however results in only first-order accuracy at the boundaries for fields whose derivatives are non-zero there (and not prescribed).Any attempt to improve the accuracy near the boundaries using higher order extrapolation -rather than linear -caused a non-physical energy growth in the Beltrami test case due to greatly increased near-boundary vorticity values that are subsequently swept into the domain.Moreover, such extrapolation appears to be inconsistent keeping sums over parcels equal to sums over grid points.Two other approaches we are considering to improve the boundary interpolation are (1) the use of 3D 'ghost' parcels in vertical halo cells, and (2) two-dimensional, compressible surface parcels (ellipses) that live in the boundary grid layer.We are currently developing the surface parcel approach, and first tests with the Beltrami flow are returning promising results.However, at present there is no mechanism to mix surface parcels with those in the interior.This may be problematic for parcel attributes such as buoyancy, where for example a descending cold region should be able to cool the surface.We are currently formulating a way to mix the attributes of parcels in the grid layer adjacent to each boundary with those of surface parcels.This requires specifying a new mixing rate parameter.
We also intend to improve the parcel initialisation developed in [1].This currently suffers from higher errors at the vertical boundaries, leading to a discrepancy in the input field values after interpolating parcel attributes to field values.This makes it difficult to compare with analytical solutions, like the linear internal wave test case presented in this paper.The surface parcel approach described above would appear to overcome this issue, leading to a consistent interpolation scheme, an integral part of the EPIC method.
Future developments include a more realistic moist physics model, the addition of precipitation, and a boundary layer to model surface fluxes.We also aim to move from the Boussinesq approximation to the anelastic formulation to be able to study atmospheric convection more realistically (including the near exponential decay of density with height).

Table E.7
The cost per time step and the average number of parcels per cell when either the maximum parcel aspect ratio λ max or the minimum parcel volume fraction Ṽ min is varied, keeping the other at its default value in Table 1.These results were generated for the rotating Rayleigh-Taylor flow test case using 64 3    rapidly intensifies at small scales, a known feature of buoyancy-driven flows.These small scales are systematically removed with increasing coarsening level.The coarsest grid resolution of 32 3 exhibits therefore the lowest r.m.s.error at late times.
Since a higher value of Ṽ min implies larger small parcels, small scales are removed, thus generally resulting in lower r.m.s.error as well.The differences between the various Ṽ min thresholds are practically insignificant, yet a value of Ṽ min = 1/10 is discouraged as volume conservation (cf.

Fig. 1 .
Fig. 1.Location of the four support points (black dots) in the major-middle plane (blue shaded ellipse) of a general ellipsoid.

Fig. 3 .
Fig. 3. Buoyancy cross sections in the xz-plane at y = 0 for various time steps of the Rayleigh-Taylor instability using EPIC at 128 3 grid resolution.

Fig. 4 .
Fig. 4. As in Fig. 3 but using the PS3D method.Here, we have used the same colourmap and range, except that values outside the initial range of buoyancy [−1, 1] are plotted with a uniform, different colour.These unphysical values can be far outside the initial range of buoyancy (see text for details).

Fig. 5 .
Fig. 5. Evolution of kinetic energy K, available potential energy A, total energy T = K + A and enstrophy ϒ for the rotating Rayleigh-Taylor test case.

Fig. 9 .
Fig. 9. Convergence analysis of the time-averaged relative error to the mean values of kinetic energy K, available potential energy A, total energy T = K + A and enstrophy ϒ. u = û sin mz e iϕ

Fig. 10 .
Fig. 10.Convergence analysis of the relative error between the numerical and the theoretical values of kinetic energy K, available potential energy A and enstrophy ϒ for the internal wave test case.The analytical values are K = 5 • 10 −7 , A = 2.5 • 10 −7 and ϒ = 7.5 • 10 −7 .

Fig. 11 .
Fig. 11.Cost in CPU seconds versus grid resolution for EPIC and PS3D for the internal wave test case.Note that EPIC is run with 8 MPI ranks and PS3D with 8 OpenMP threads.

Fig. 12 .
Fig. 12. Cross sections in the xz-plane at y = 0 (upper panel) showing the vertical vorticity difference from the analytical result ζ (t) normalised by ζ rms (t), and (lower panel) the r.m.s.error normalised by ζ rms (t).Here results are shown for the simulation with 256 2 × 64 cells.Note that the fields are interpolated to exactly match the labelled times.

Fig. 14 .
Fig. 14.Evolution of kinetic energy K and enstrophy ϒ for the Beltrami test case.

Fig. 15 .
Fig. 15.Cross sections of |ω| in the xy-plane at z = 0 at times 9, 10 and 11.The top row is for EPIC and the bottom row is for PS3D.

Fig. 19 .
Fig. 19.Cross sections of the specific humidity field q in the xz-plane (zoom) in EPIC, MPIC and MONC simulations with different numbers of grid points.

Fig. 20 .
Fig. 20.Liquid water specific humidity histogram in EPIC, MPIC and MONC simulations with different numbers of grid points at t = 6.The bin width is 0.002.

Fig. 21 .
Fig. 21.Root-mean-square error of the EPIC, MPIC and MONC liquid water specific humidity histograms at t = 6, as compared to the histogram in the highest-resolution reference EPIC, MPIC and MONC simulations.

22 .
Histogram of the vorticity magnitude, |ω|, in the EPIC, MPIC and MONC simulations with different numbers of grid points at t = 6.The bins widths correspond to 5% increases in vorticity magnitude (i.e.log spacing).

Fig. E. 25 .
Fig. E.25.Gridded relative r.m.s.volume error for different maximum parcel aspect-ratios λ max , at a fixed grid resolution of 64 3 and minimum parcel volume fraction Ṽ min = 1/20.

Fig. E. 26 . 2 z.
Fig. E.26.Gridded relative r.m.s.volume error for different minimum parcel volume fractions Ṽ min , at a fixed grid resolution of 64 3 and maximum parcel aspect-ratio λ max = 4.
Fig. E.26) is less well preserved compared to the others.The recommended choice that keeps the computational cost low but accuracy high is Ṽ min = 1/20.A maximum parcel aspect ratio of λ max = 4 also ensures a low cost-accuracy ratio.

Fig. E. 27 .
Fig. E.27.Cross section of the buoyancy field for the simulation with 384 3 grid cells at y = 0 in the xz-plane with intersected ellipses at t = 6.The zoomed windows correspond to widths of 192, 96, 48 and 24 grid cells in both directions.

Table 3
Relative total energy loss (in per mille) for the internal wave test case between the initial time t = 0.0 and the final time t =

Table 5
Approximate physical constants and parameters used in the moist bubble setup.The exact values are calculated in the programme setup routine.

Table 6
Runtime in seconds of the moist bubble simulations used in this study.The numbers annotated with indicate I/O-server cores.
grid cells.Each simulation ran with 64 cores.