{'fraction_non_alphanumeric': 0.035612711759614196, 'fraction_numerical': 0.021392358105601583, 'mean_word_length': 4.475964793500339, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Direct simulation of turbulent entrainment due to a plume impinging on a density interface

The law for turbulent entrainment due to plumes and jets impinging on a density interface is subject to significant uncertainty, with reported differences in entrainment rates up to a factor of 10. We report preliminary results obtained by Direct Numerical Simulation which are part of a PRACE project on turbulent entrainment carried out on JUGENE at Jülich, Germany. Various interface tracking methods are discussed and the entrainment coefficient is determined.


Introduction
Entrainment laws, which couple an entrainment velocity to a turbulence intensity and density contrast between the layers, are subject to very significant uncertainties. Indeed, laboratory measurements of entrainment rates in identical apparatus vary by orders of magnitude (Fernando, 1991;Coffey & Hunt, 2010). This throws into question many models and modeling approaches that rely on a specification of entrainment rate. Due to the continuous increase in computational capacity, supercomputers have now become sufficiently powerful to perform Direct Numerical Simulation of turbulent entrainment at realistic Reynolds numbers. The objective of this work is to identify a reliable and efficient method to determine the entrainment rate using DNS for a plume impinging on a density interface. In principle, the availability of the entire three-dimensional velocity and temperature field allows for an accurate estimation of the rate at which fluid is entrained. We explore several different methods to determine the entrainment rate and show preliminary results for the interface position as a function of time.

Simulation details
A coordinate system is adopted in which the x and y coordinates represent the lateral directions, and the z coordinate is vertical and points upwards. The domain is of size L × L × H, where L = 72b 0 , H = 36b 0 , and b 0 is the radius of the buoyancy source situated in the centre of the bottom boundary at (0, 0, 0). The stably stratified environment prescribed at time t = 0 consists of two layers with buoyancies ∆ a and ∆ b > ∆ a , respectively. The sharp density interface is initially positioned at h 0 = 24b 0 and ∆ a = 0 without loss of generality. The code for DNS solves the Navier-Stokes equations in the Boussinesq approximation and the enthalpy equation is formulated in terms of the relative temperature θ which is related to buoyancy by ∆ = βgθ, where β and g are the expansion coefficient and gravitational acceleration respectively. Details of the code can be found in van Reeuwijk (2007); van Reeuwijk et al. (2008). The boundary conditions are no-slip and insulating. The heat source on the bottom boundary maintains a constant temperature flux φ 0 , the associated integral source buoyancy flux of which is B 0 = πβgφ 0 b 2 0 . The plume which forms as a result of the heat source is tripped to turbulence by adding small random fluctuations to the velocity field in the cells immediately above the source (Plourde et al., 2008). The source Reynolds number Re 0 based on B 0 and source diameter 2b 0 is where ν is the kinematic viscosity. We further define the Prandtl number Pr = ν/κ where κ is the thermal diffusivity, and the interfacial Richardson number Here, ∆ 1 , b 1 and w 1 are the plume buoyancy, radius and velocity at the position of the density interface z = h, respectively.
Results are presented of a simulation with Re 0 = 1000, Pr = 1 and Ri = 4.7. The simulation is part of a project on turbulent entrainment supported by the Partnership for Advanced Computing in Europe (PRACE), had a resolution of 1536 × 1536 × 1024 cells and was performed on 32,768 processes on JUGENE at Forschungszentrum Jülich, Germany. The simulation cost was 6 × 10 5 CPU-h. The grid is equidistant with ∆x = ∆y = ∆z, where ∆x i is the increment in the x i direction. To avoid overturning of the plume outflow in the region below the interface (cf. Kaye & Hunt, 2007), the ratio of layer width L to depth h should be no less than 2. As an immediate consequence, the cross-sectional area of the plume πb 2 1 ≪ L 2 (for the current simulation, πb 2 1 /L 2 ≈ 0.01). The plume rising from the heat source and impinging on the density interface was extensively compared with plume theory (Morton et al., 1956;Turner, 1979). Both the entrainment coefficient and the powerlaw dependencies of volume, momentum and buoyancy fluxes agreed well with plume theory.

Determination of the entrainment rate
For a plume impinging on a density interface, the entrainment law is assumed to be of the form (Turner, 1968;Baines, 1975;Kumagai, 1984;Lin & Linden, 2005) where Q * is the entrained volume flow rate divided by π and E the entrainment coefficient, which is assumed to depend on Ri only. In laboratory experiments, Q * is usually inferred from the evolution of the layer depth h Q * = L 2 π dh dt .
(2) An advantage of DNS over experiments is the exact control of boundary conditions, and the availability of the entire three-dimensional velocity and temperature fields. Consequently, with DNS Q * can be determined in several ways (Sullivan et al., 1998), e.g. by assuming that: (i) the interface h is located at the inflection point of the horizontally averaged temperature profile; (ii) the interface is given by the average of the individual z-locations (for each x, y) at which ∂θ/∂z is maximal; and (iii) the interface is given by the average of the individual z-locations which correspond to the first maximum in ∂θ/∂z (searching downwards from z = H) exceeding a threshold value 1 . Method (iii) was developed because method (ii) sometimes misidentifies the location of the interface in the plume impingement area due to large temperature gradients near the plume source, and (iii) is therefore considered to be more reliable.
1 The threshold is necessary to avoid misidentifications due to very small spurious local maxima from contaminating the position of the interface. Tests show that method (iii) is insensitive to the exact value of the threshold. ]. The grey band shows the spatial standard deviation of the interface obtained with method (iii).
The evolution of h as a function of time is shown in Fig. 1. Method (i) and (iii) are in good agreement, although the signal from method (i) is noisier than from method (iii). Method (ii) underestimates the location of the interface due to misidentifications in the plume impingement area. Four stages can be identified: 1) the interface is stagnant until the plume has risen through the lower layer and impinges on the interface at t/t * = 5, where t * = b 1 /w 1 ; 2) the outflow of the plume and interfacial gravity waves spread radially until the edges of the domain are reached around t/t * = 30; 3) the interfacial waves reflect off the side walls, which results in an approximately static position of h for about 10 t * ; and 4) a secondary layer forms just below the interface for t * > 40. It is the stage 4 development that corresponds to the true turbulent entrainment across a density interface. The grey band around the average interface position for method (iii) in Fig. 1 is the spatial standard deviation of the interface position at each time. This clearly shows the difficulties in pinpointing the location of the interface, as the spread at a single time instant is comparable to the change in the interface location over 70 t * . The entrainment coefficient E is the slope in Fig. 1; for 50 < t/t * < 70 we obtain E = 0.12. As a comparison, for Ri = 4.7, the entrainment relation from Baines (1975), Kumagai (1984) and Lin & Linden (2005) predict E = 0.07, 0.05 and 2.0, respectively.
In conclusion, we have performed direct simulation of a plume impinging on a density interface and tested three interface detection methods. Each method centres around tracking an interface position informed by the temperature field. We note that even in this highly controlled environment, the entrainment rate, a quantity which can only be inferred, is sensitive to the way in which the interface is defined.