1 Motivation

Mixing closure in computational models of turbulent combustion is typically implemented by partially or fully intermixing pairs or groups of notional fluid parcels selected from a parcel population that discretely instantiates the joint probability distribution function (PDF) of the thermochemical variables that are time advanced by the model [10]. In the absence of specially formulated constraints, this approach allows the intermixing of parcels with highly dissimilar states, which is both unphysical in principle and detrimental to model performance in practice.

One such constraint that has proven effective is to intermix only parcel pairs that are close, by some criterion, in a metric space defined on the manifold of thermochemical states [26]. This constraint avoids the noted artifact, but the resulting formulation lacks other potentially important features of turbulent mixing. One indication of this limitation is that it allows intermixing of any pair of thermochemically close parcels, but in reality, parcels that are thermochemically close are not necessarily spatially adjacent. Consequently, individual fluid parcels do not necessarily exhibit physically realistic mixing time-histories and multi-parcel joint statistics are not necessarily well represented.

Guided by turbulent cascade phenomenology, the present study introduces a minimal representation of the time-evolving occurrences of adjacency of parcel pairs (or more generally, p-tuples) in order to capture the leading-order effects of this phenomenology on mixing time-histories, multi-parcel joint statistics, and other such observables. Operationally, the role of this representation is to partition the parcel population into pairs at each instant, such that the parcels in each pair are deemed to be adjacent and therefore subject to intermixing at a rate that is based on the molecular diffusivity and a specified local length scale.

Details of these and other features of the mixing representation can be formulated in various ways, resulting in choices to be made in order to instantiate the described mixing representation as a complete model suitable for analysis or numerical implementation. Such model instantiations and their application to cases of interest are the ultimate goal of the research, but here the emphasis is on the generic representation and its portrayal of turbulent mixing. This emphasis is intended to convey the breadth of possibilities within the mathematical framework that is described. Application-specific details and associated analysis and numerical results will be presented elsewhere.

The partitioning of parcels into adjacent pairs is done through the assignment of parcel indices, such that each even-numbered parcel and the odd-numbered parcel that follows form a pair. Time advancement reorders parcel indices, an algorithmically simple operation that adds negligible cost to a typical numerical simulation of turbulent non-reactive or reactive mixing yet, as will be shown, can substantially augment the representation of turbulence phenomenology in a manner that improves the fidelity of the mixing simulation.

Though designed to improve turbulent combustion closures, it is shown that the method introduced here, termed the hierarchical pair-swapping (HiPS) representation of turbulent mixing, can also be useful as a self-contained model, not only to simulate turbulent mixing, but also, with the incorporation of velocity (see Sect. 5.1), to simulate turbulence more generally. In both contexts, the model has features in common with previous hierarchical models of turbulent flow, and also some essential differences designed to address the noted objectives.

2 Hierarchy Geometry and Associated System Variables

The scale-space representation that is adopted involves an acyclic graph, or tree. Tree geometry has been used for reduced modeling of turbulence intermittency [2, 4, 5, 14, 22], in some instances formulated as generalizations of the zero-dimensional shell models [6, 15] in which the scale space of turbulence is parameterized by wavenumber modulus. In such formulations, the system state consists of real variables representing velocity values that reside at each node of the tree and are time advanced. In the present approach, the system state consists of thermochemical state variables (and optionally velocity; see Sect. 5.1) associated with fluid parcels that reside at tree endpoints (the nodes at the base of the tree, below which there are no further links or nodes). Variables residing at other nodes of the tree are auxiliary variables in the present context. They are used to specify operations performed on the tree structure as time advances. In this regard, the formulation is fundamentally different from previously formulated hierarchical models in which dynamical variables at all levels, representing either filtered quantities or terms in a notional mode decomposition, are time advanced.

It is convenient to specialize to a binary tree. (Generalization to branching factors p>2 is straightforward; see Sect. 5.2.) The apex of the tree corresponds to the largest length scale l 0 that is represented. l 0 can be the turbulence integral scale or the spatial resolution scale of a coarse-grained flow computation in which the tree is imbedded to provide subgrid closure. (Other application-specific choices of l 0 are possible.)

Below the apex, successive levels of the tree, with integer labels n>0 (n=0 denotes the apex) correspond to length scales l n =A n l 0, where A<1 is a parameter of the representation. Here, increasing n corresponds to decreasing length scales.

For now it is assumed (but not required in all instances) that all paths downward from the apex terminate at the same level n=N, which is thus the base of the tree. Counting the apex, there are thus N+1 tree levels. Additionally, all fluid parcels are assumed to have the same nominal volume V. Sub-trees at a given level n connect downward to half as many fluid parcels as sub-trees at the level above (level n−1), so nominally they subsume half as much fluid volume. This interpretation suggests using A=2−1/3 so that the cube of the length reduction per level step matches the volume reduction.

This requirement can be relaxed if the collection of parcels is interpreted as corresponding to a lower-dimensional cut through 3D space. Symmetry considerations require this interpretation in some cases (see Sect. 5.1). Moreover, smaller A reduces computational cost by reducing the number of parcels.

To proceed further, a labeling of the nodes and links of the tree is needed. Level n of a p-adic tree has p n nodes. For p=2 (dyadic or binary tree), two links extend downward from each node. These links connect the node to two nodes at the next level down. Thus, the coordination number of the tree network is 3, except for the apex node (which has no link extending downward to it, so its coordination number is 2) and the nodes at the base of the tree, whose coordination number is 1.

Notation introduced by Cleve et al. [7] is adopted. At each level, nodes i are numbered sequentially from zero going from left to right. For given n>0, if i is written in binary notation with enough leading zeros to form an n-digit binary number, then this numbering encodes the path from the apex to node i. Namely, reading this binary label from left to right, successive values 0 (1) are interpreted as choosing the leftward (rightward) downward fork, starting from the apex. (On this basis, the consistent notation for the apex is the empty-set symbol ∅.) Because there is a unique path leading to each node, it is valid and convenient to assign the node label to the link connecting it to the preceding (higher) level.

The variable controlling system advancement is a node time scale τ. In the inertial range, τ=(l n /l 0)2/3 τ 0 is assumed, where τ 0 is a representative flow time scale corresponding to length scale l 0 and the exponent enforces the inertial-range scaling relationship between eddy length and time scales (see Sect. 3.1). An extension to node-dependent, time-varying τ is described in Sect. 5.1.

In order to represent mixing of scalars with high Prandtl number Pr or Schmidt number Sc (defined in Sect. 4.2.2), the tree extends to length scales l as small as needed to resolve the Batchelor scales of all scalars. (Additionally, the possibility of using finer resolution than this is discussed in Sect. 4.2.1.) Thus the tree might need to resolve dissipative as well as inertial-range flow scales.

In Sect. 4.2.2, the determination of the Kolmogorov microscale η and the associated time scale τ η marking the transition from the inertial range to the dissipative range is discussed. Owing to viscous suppression of flow perturbations in the dissipative range, flow structure smaller than scale η is negligible for present purposes. Accordingly, τ η is the relevant flow time scale in the dissipative range, so dissipative-range tree levels are assigned the time scale τ η =(η/l 0)2/3 τ 0 (subject to round-off of η due to the discreteness of hierarchy levels).

3 Representation of Turbulent Advection

3.1 The Turbulent Cascade from a Mixing Perspective

In the context of the goal of simulating parcel-pairing evolution (Sect. 1), advection is an auxiliary operation serving to modify fluid-parcel pairings as time advances. Nevertheless, by incorporating turbulence phenomenology needed to do this realistically, the chosen approach embodies, in a parameterized way (and in a more self-contained way in the extension described in Sect. 5.1), key features of the inertial-range turbulent cascade.

Two features of the turbulent cascade that are crucial from a mixing viewpoint, as well as more generally, are the gradual reduction of fluid-parcel length scales, which is a down-scale process, and the gradual increase of the rate of separation of a pair of initially close parcels, which is an up-scale process. These two features are manifestations of a fundamental cascade property, the so-called sweeping of the small scales by the large scales [8].

This sweeping is one of the consequences of the non-dissipative, quasi-steady down-scale transport of turbulence kinetic energy (TKE), such that the rate ϵ of TKE dissipation per unit mass is scale invariant in the inertial range. ϵ is termed a dissipation rate because it corresponds to the ultimate dissipation of TKE at viscous scales, but in the inertial range it represents the down-scale rate of transfer of TKE and at the integral scale it represents the production of TKE by external mechanisms.

The scale invariance of ϵ enables the association of a time scale τ with any given inertial-range length scale l. Based on the dimensional relation ϵ∼l 2/τ 3 and the lack of l dependence of ϵ, τ∼l 2/3 is obtained. Likewise, a velocity v∼l/τ∼l 1/3 associated with length scale l is obtained. Here, powers of ϵ needed for dimensionally consistent scalings are omitted.

A fluid parcel of given size l can thus be viewed as moving relative to its size-l neighborhood at a velocity v(l)∼l 1/3. Additionally the parcel is an element of larger parcels corresponding to length scales l +>l that move relative to their size-l + neighborhoods at speeds v(l +)>v(l). This inequality indicates that scale-l relative motions are slow relative to those at scale l +, so the larger-scale motions appear to sweep smaller-scale structures past an observation point with relatively little disturbance of those structures by small-scale relative motions.

Another implication of the scaling τ∼l 2/3 is that small-scale motions, though occurring at low velocity, respond to larger-scale input quickly relative to the longer time scales τ associated with that input. This self-consistently justifies the assumed quasi-steadiness of the inertial range, even when the turbulent cascade is subject to large-scale transient forcing. This rapid local equilibration allows the possibility of universal self-similar turbulent-cascade structure, a feature that has been extensively demonstrated and is the basis of the present approach.

The disparity of time scales at widely separated length scales suggests relatively weak dynamical coupling of widely separated (vs. less separated) length scales. This supports the notion of gradual breakdown of length scales through a sequence of moderate scale reductions.

The separation of parcel pairs can likewise be understood on the basis of this picture. Denoting their separation by x, a diffusion process d〈x 2〉/dt∼D is assumed, where angle brackets denote an ensemble average and D is a separation-dependent diffusion coefficient. Dimensionally, D∼l 2/τ∼l 4/3, where l denotes a representative parcel-pair separation (l 2∼〈x 2〉), yielding the Richardson dispersion law l 2∼t 3 [9]. Here, length scales larger than l are assumed to make no leading-order contribution to dispersion because they sweep the parcel pair without contributing significantly to their separation, and length scales smaller than l are deemed unimportant because their associated velocities are smaller than the scale-l velocity.

This heuristic description of the inertial-range cascade is supported by analysis and empirical data that has been discussed extensively, e.g., in [11]. The manner in which this phenomenology is represented within the present framework is explained next. In this regard, mention of fluid parcels in the discussion of the inertial range is intended to be illustrative rather than to imply that fluid parcels are associated with nodes throughout the inertial range of the tree. As noted in Sect. 2, fluid parcels reside only at the base of the tree.

3.2 The Advection Mechanism

Within the hierarchical geometry, fluid volumes associated with various length scales (tree levels) are not individual parcels, but rather, collections of fluid parcels at the bases of the corresponding sub-trees. In this context, a mechanism for implementing advection is suggested by the intended use of the formulation. Namely, the goal is to specify the time evolution of parcel indexing for the purpose of identifying adjacent pairs. In effect, a process involving permutations of parcel indices is needed.

These considerations suggest the following representation of advection. Every node that is not at the base of the tree is subject to the following process applied to its associated sub-tree, meaning the sub-tree whose apex is the node. A time sequence of permutations is applied to the sub-tree, where permutation occurrences are sampled from a Poisson process in time with mean time τ between occurrences. Here, τ is the time scale associated with the node, as described in Sect. 2.

The instantiation of this method in a numerical simulation of mixing may involve an empirical coefficient in the expression for τ. Henceforth, the possible need for empirical coefficients in various expressions is understood but not mentioned further.

Links downward from any node i (in the binary notation of Sect. 2) connect it to nodes i0 and i1, which in turn connect downward as follows: from i0 to nodes i00 and i01 and from i1 to nodes i10 and i11. (To interpret this notation, the binary representation of a particular i value should be substituted, e.g., for i=0010101110, i01 denotes i=001010111001. Recall that leading zeroes must be included as needed to specify uniquely the path to a node.) Owing to the unique association of a node with its sub-tree, the index i is used to refer to either.

In this notation, the permutation applied to sub-tree i consists of the random selection (by a Bernoulli trial) of either sub-tree i00 or sub-tree i01 with equal probability of either, another such random selection of either sub-tree i10 or sub-tree i11, and finally, swapping of the two selected sub-trees. Denoting the selected sub-trees as i0q and i1r (where q and r are each determined by the aforementioned independent random sampling of either 0 or 1), then the formal implementation of the permutation consists of the relabeling of every node i0qs, where s is any binary integer including any number of leading zeroes, as i1rs, and the concurrent relabeling of every node i1rs as i0qs. This operation is compactly denoted i0qs↔i1rs, where s in any such expression has the same meaning as in the previous sentence. Where a designated path extends below the base of the tree and hence to a non-existent node, the relabeled path again leads to a non-existent node, so no inconsistencies are introduced by ignoring termination of paths at the base of the tree. Figure 1 shows a five-level tree structure and representative swaps that are possible within this structure.

Fig. 1
figure 1

A five-level binary tree and representative permutations. Squares are parcels at the base of the tree (level 4, based on top-down level numbering starting from 0) and triangles are nodes at other levels. Line segments are node connections. Circles around individual nodes indicate nodes selected for illustrative permutations. One permutation is illustrated per level at which permutations are possible in a five-level tree. At level 0, the only available node is ∅ (the apex of the tree). The illustrated level-0 permutation is the swap of level-2 sub-trees 01 and 10, which are circumscribed by the large dashed ovals. A level-1 permutation is similarly shown, where the level-1 node index is 0 and the swapped (level-3) sub-trees are 000 and 011 (in mid-size dashed ovals). Finally, the illustrated level-2 permutation at node 11 consists of a swap of level-4 sub-trees 1101 and 1111 (in small dashed ovals), which are individual parcels. Only swaps of individual parcels can modify parcel pairings

Full size image

If the only system property of interest is the time evolution of the collection of fluid-parcel states irrespective of parcel labeling (where, like a node, a parcel is labeled by the path leading to it), then either q or r can be set to a fixed value, say 0, without loss of generality, so only one Bernoulli trial is needed. This is because the additional degree of freedom is equivalent to allowing the two sub-trees i0 and i1 to be swapped, which affects parcel indexing but not parcel pairing, where only the latter affects mixing-induced evolution of parcel states. However, configurations are possible in which it can be meaningful to associate parcel indexing with spatial ordering, so this simplification is not always applicable.

3.3 Parcel-Swapping Phenomenology

The manner in which the representation of turbulent advection effects on mixing described in Sect. 3.2 reflects various aspects of turbulent cascade phenomenology is now explained. By construction, HiPS implements fluid-parcel rearrangements over the full range of relevant length scales l, where rearrangement occurrences are governed by a commensurate range of time scales τ with the appropriate l dependence. Details of the mixing representation are deferred until Sect. 4. Here it suffices to say that the contribution of turbulent advection to mixing is to break down large-scale structures and to bring the resulting small-scale structures into close enough proximity to each other to allow prompt molecular mixing.

To illustrate the HiPS analog of this mechanism, assume that fluid parcels within one of the two sub-trees (labeled 0 and 1) of the apex all have the same given initial composition, while the parcels within the other sub-tree all have some other given composition, as in Fig. 2. (Motivated by laboratory analogs of this initial-value problem, it is termed the scalar-mixing-layer configuration.) Then the only permutation that changes the system state is the permutation associated with node ∅ (the level-0 apex node). The associated time scale τ is the slowest time scale of the hierarchy. After the passage of an order-τ time interval, such a permutation occurs with order-one probability. Then level-1 permutations governed by correspondingly smaller time scales can change the system state. By induction, a permutation one level lower (meaning n larger by unity) than any prior permutation is the only mechanism that can induce further scale breakdown. As in the inertial-range cascade, the largest-scale breakdown steps are rate-limiting due to the decrease of τ with increasing n.

Fig. 2
figure 2

Top: Mixing-layer initial configuration, involving two parcel states denoted by striped and open rectangles, respectively. Also shown is the level-0 permutation illustrated in Fig. 1. Any level-0 permutation changes the system state, but permutations at lower levels (higher level indices) leave the system state unaffected. Bottom: Configuration after swapping the sub-trees involved in the level-0 permutation, showing a level-1 permutation that can now change the system state. After this permutation occurs, level-2 permutations can change the system state in a way that creates pairs whose parcels have different states, thereby allowing mixing to begin changing parcel states

Full size image

After the first level-n permutation, a second permutation at the same level-n node before a level-(n+1) permutation occurs downward from that node can undo the effect of the first permutation at that node, which would then require a third such level-n permutation to enable further scale breakdown in the sub-tree. Such a reversal is qualitatively analogous to backscatter, but the fidelity of this backscatter representation needs further investigation.

In the absence of molecular mixing, this breakdown process can proceed ad infinitum, corresponding to a non-terminating tree, which is a mathematically well-behaved limiting case of the tree geometry, though not necessarily a physically realistic representation of the purely advective regime (see Sect. 4.2.1). Only in the presence of molecular mixing is it valid to impose a finite termination of the tree based on sufficient resolution of all physically significant scales of parcel-property variation.

Next, consider the dispersion of parcels that are initially close in the sense that they are contained within a large-n (small-scale) sub-tree that is termed the minimal spanning sub-tree if there is no larger-n sub-tree containing the parcels. The level n of the minimal spanning sub-tree is termed the parcel proximity, where larger n corresponds to closer proximity. (Proximity defined in this way is equivalent to the ‘ultrametric distance’ between tree nodes [7].) As illustrated in Fig. 3, the separation of parcels with proximity n can be increased only by level-(n−1) permutations because larger-n permutations leave them within the same minimal spanning sub-tree though possibly differently positioned within the sub-tree, and any smaller-n permutation repositions the entire minimal spanning sub-tree without modifying its internal structure. Separation increase by a level-(n−1) permutation results in parcel proximity n−1, allowing further separation only by a level-(n−2) permutation, and so forth (but a subsequent level-(n−1) permutation can reduce the separation, resulting in particle proximity n). This stepwise process, in conjunction with the inertial-range dependence of τ on l, implies consistency of HiPS with the heuristic picture of Richardson scaling.

Fig. 3
figure 3

Top: The proximity of the striped and checkerboard parcels can be changed from 3 to 2 by any swap associated with node 01 (circled). One such swap is indicated by the dashed ovals. A pair swap associated with any other node would not change their proximity. Bottom: After the indicated swap, the proximity of these parcels can be changed from 2 to 1 by any swap associated with node 0 (circled), or can be changed from 2 to 3 by some but not all of the swaps associated with node 01. These are the only swaps that can change their current proximity

Full size image

The diffusion picture of dispersion outlined in Sect. 3.1 corresponds to a pair-separation random-walk process in which parcels can approach as well as separate from each other, though with increasing separation on average. In HiPS, a permutation can change parcel proximity only by ±1. Thus, both approach and separation are stepwise in n. For any dependence of τ on n, parcel proximity then evolves as a random walk with n-dependent rates of upward and downward stepping in n.

With regard to pair dispersion, HiPS thus corresponds to a discrete-state stochastic process whose continuum limit can be formulated as a separation-dependent diffusion process describing the time evolution of the PDF of the pair separation. In this respect, it is consistent with the conventional theoretical approach to pair-dispersion analysis [9]. Owing to the dependence of HiPS dispersion only on the assumed hierarchical geometry and the n dependence of τ, comparison to dispersion data is a suitable basis for determining the coefficient that is implicit in the dependence of τ on l.

Multi-parcel multi-time statistics of turbulent dispersion are a more subtle matter. Non-trivial extremal properties of multi-parcel turbulent dispersion have been identified and analyzed [24]. The HiPS representation of these phenomena will be investigated in future work. The purely advective phenomenon of dispersion illustrates the potential relevance of HiPS to phenomenology extending beyond turbulent mixing.

In comparison to continuous-in-time fluid flow, the instantaneous nature of swaps is unphysical, but it does not necessarily impair the formulation relative to its intrinsic idealizations such as the use of tree rather than physical-space geometry. Indeed, the stepwise progression of both breakdown and dispersion in HiPS mitigates possible artifacts, with the notable feature that this stepwise behavior mimics salient features of turbulent cascade phenomenology.

Though swaps are instantaneous, those at scales well above the base of the tree have no immediate effect on mixing. Their effects are felt gradually as the modified large-scale tree connectivity affects the outcomes of successively smaller-scale swaps until parcel pairings at the base of the tree, which influence mixing (see Sect. 4), are affected. Swaps at the base of the tree occur at a time scale of the same order as the molecular-mixing time scale, so algorithmically their effect is not qualitatively different from conventional time discretization of the evolution equation governing physical-space scalar advection and diffusion. Thus, no major artifacts are anticipated to result from the time discontinuities per se, but these are not the only significant idealizations within the overall framework.

4 Representation of Molecular Mixing

4.1 The Mixing Mechanism

The implementation of molecular mixing in HiPS is basically the same as in conventional mixing models used in engineering applications. This has the practical benefit that HiPS can readily be incorporated into such models with little modification or added computational cost, while offering potential performance improvement stemming from its additional physical content. The additional physical content is the introduction of parcel proximity, albeit not in the usual spatial coordinates, and time evolution of proximity that reflects mixing-relevant turbulence phenomenology, as explained in Sect. 3.3.

The HiPS (and conventional) implementation of molecular mixing involves the exchange of some portion of the contents of a pair (or more generally, a p-tuple; see Sect. 5.2) of fluid parcels. In itself, this is not qualitatively different from typical model-free numerical implementation of flux-based diffusion, and in fact has been used as a numerical method for solving the diffusion equation [27]. The modeling arises from the treatment of parcels as Lagrangian entities whose encounters are governed by a stochastic process rather than the exact equations of fluid motion.

Like the occurrence of swaps, the occurrence of these exchanges between a given parcel pair is governed by a Poisson process in time, with a mean time between occurrences that is based on a diffusional time scale. In the simple case that diffusive transport is governed by a single diffusion coefficient κ, the time scale τ d =l 2/κ is suitable, where l is the length scale of the level at which the parcels reside, specifically, l N for the tree structure that has been assumed.

As in existing mixing models [10], details of the exchange mechanism can be tailored to the complexity of the case-specific molecular-transport properties and to the degree to which representation of this complexity is needed in order to achieve the desired fidelity. Because suitable methods are available, the possibilities are not discussed here. For now, the simplest option, which is to fully mix the two parcels so that their final compositions are identical, is adopted. To obtain a simple idealization of mixing effects, this option can be further specialized by assuming τ d to be smaller than advection time scales, such that it can be assumed to a good approximation that parcels fully mix at the instant that they first become a pair and are therefore eligible to mix. Then mixing becomes a direct adjunct to the advection process.

4.2 Turbulent Mixing Phenomenology

4.2.1 Determination of the Number of Tree Levels

Fluid parcels reside at the base of the tree, and as noted in Sect. 1, an adjacent pair consists of an even-numbered parcel and the parcel that follows it in the indexing scheme. In a tree with its base at level n=N, the adjacent pairs are those with proximity N−1 (as defined in Sect. 3.3). It is noted in Sect. 2 that the choice of N depends on the spatial refinement needed to resolve the scale at which the assumption of parcel homogeneity is valid to a sufficient approximation, i.e., there are no significant unresolved small-scale property fluctuations. The determination of N on this basis is considered, for now with the restriction that the tree base is at a single level, i.e., that all branches terminate at the same level N. This ostensibly straightforward question touches on subtle aspects of turbulent mixing phenomenology and therefore is examined in detail, with emphasis on the HiPS representation of this phenomenology.

In the numerical advancement of the exact equations of fluid motion, increase of the spatial refinement of the computational mesh causes no loss of fidelity, so the only practical restriction is computational cost. However, the analogous increase of N in HiPS can potentially reduce fidelity due to its hierarchical structure. For example, consider the transient downscale development of the scalar mixing layer discussed in Sect. 3.3. Suppose that the relationship between the n-dependent τ and τ d values is such that the crossover from τ d >τ (advection dominance) to τ d <τ (diffusion dominance) occurs in the vicinity of some n=n B (B for the level corresponding to the Batchelor length scale, which in HiPS notation is \(l_{n_{B}}\), but for brevity will be denoted l B ). At scales below l B , diffusional dissipation of advection-induced state differences between adjacent parcels is increasingly effective at suppressing state fluctuations, such that these fluctuations rapidly become negligible with decreasing scale. Thus, any value of N that is substantially larger than n B will generally provide sufficient resolution, subject to intermittency considerations.

This applies also if N≫n B , but then an artifact can arise. In the HiPS scalar mixing layer, mixing of dissimilar states occurs only when advective length-scale breakdown brings dissimilar parcels into adjacency at the base of the tree. Thus, breakdown to level n B , which should induce mixing physically, does not automatically do so in HiPS. HiPS requires breakdown to level N in order for dissimilar states to become adjacent according to the HiPS definition and therefore subject to mixing. Therefore the choice N≫n B adversely impacts the HiPS representation of scalar spectral scalings (see Sect. 4.2.2). Also, though mixing occurs to the physically required extent, its occurrence is not at the time of advective breakdown to level n B but at a later time due to the additional required breakdown to level N. In general, this additional time is small relative to overall rate-controlling large-scale advective times. However, there can be exceptions, albeit of little practical relevance. Namely, in the dissipative (viscosity-dominated) scale range, τ is independent of n, which can result in an arbitrarily large time for breakdown from level n B to level N for sufficiently large N. This is noted mainly as a point of principle to show that arbitrary increase of N is not an unequivocally safe numerical procedure, and additionally, that the formal limit N→∞, even if mathematically well behaved, is generally not a physically valid regime of HiPS.

This raises the question of the degree to which HiPS is uniquely defined for this illustrative case and more generally. For the scalar mixing layer, if there is a range of N values greater than n B in which features of interest are insensitive to N, then this range is effectively an ‘intermediate asymptote’ in which HiPS is unique within some tolerance. Which features are of interest is somewhat subjective as well as case dependent. In this regard the situation is analogous to subgrid closures of under-resolved 3D turbulence simulations, in which physical modeling often cannot be decoupled from numerical implementation details.

Extending this analogy, it is possible that mesh adaptivity used in 3D simulations would likewise be advantageous in numerical implementations of HiPS. This might involve adding levels to individual termination points of the tree by adding two links below a given parcel, then removing the parcel and replacing it by two parcels, each half the volume, at the newly created termination points.

Non-diffusive parcel dispersion, considered in Sect. 3.3, is revisited with reference to the determination of N. In physical-space simulations, it can be useful for the purpose of identification of dispersion patterns to track a large number of initially closely spaced parcels in each simulated realization [24]. For HiPS, this consideration indicates a sense in which taking the limit N→∞ can be physically valid and possibly advantageous. Suppose some fluid volume V corresponding to a length scale l n and associated tree level n is deemed to be marked initially by being assigned a different composition than the remaining fluid. In HiPS, this is represented by marking all the fluid parcels in one level-n sub-tree. Suppose that N is not much greater than n, so there are only a few marked parcels. The trajectories of these parcels can characterize statistically the rate of parcel separation, but they cannot adequately reflect the degree to which parcel dispersion fills space at any instant, as measured, say, by the probability that any chosen fluid volume is free of marked fluid at that instant (where lower probability corresponds to a greater tendency to fill space). For this purpose, the number of marked parcels must be at least large enough so that every chosen fluid element of given volume can contain at least one marked parcel (so that a completely uniform spatial distribution of parcels is detectable). By this criterion, the needed number diverges as the chosen volume vanishes, so an infinite number of parcels is ultimately needed to detect arbitrarily close approach to a space-filling distribution of parcels. In this context, use of the largest feasible value of N is desirable for dispersion studies.

The need for an unbounded number of parcels to detect the approach to a space-filling distribution reflects an important feature of turbulent mixing that has not yet been discussed, namely, the stretching, or extensive strain, that is a necessary adjunct to the compressive strain that reduces fluid-parcel length scales in the inertial range (and at smaller scales). When a fluid parcel of given volume is compressed in one or more directions, it must be stretched in at least one direction due to volume conservation. This stretching is the mechanism by which a fixed-volume fluid parcel can undergo an unbounded increase of its surface area and thus ultimately approach a space-filling morphology. HiPS has no explicit representation of stretching. As illustrated by the dispersion scenario, subdivision of a fluid volume into multiple parcels by choosing a large value of N enables a surrogate representation of stretching effects, at least for some purposes.

Comparison to the treatment of stretching in fully resolved direct numerical simulation (DNS) of turbulent scalar mixing provides a perspective on stretching in HiPS. In DNS, it is possible to follow the trajectory of advected marker particles (flow-following points) in a Lagrangian manner by interpolating the resolved flow to obtain velocities at the marker locations. Owing to the high-order accuracy of DNS numerics, this can provide Lagrangian resolution much finer than the mesh scale provided that the turbulent motions are well resolved on the mesh. In contrast, HiPS is not a physical-space flow solution and thus offers no means of resolving length scales finer than the parcel size. From this viewpoint, DNS is a more effective tool than HiPS for studies of stretching effects. Nevertheless, it will be interesting as a point of principle to determine the extent to which the HiPS stochastic representation of these effects emulates stretching in physical-space turbulence.

The comparison to DNS highlights an additional point. Provided that DNS resolves the Batchelor scale(s) of the advected scalar(s), DNS accurately time advances the system without the need to consider any smaller-scale effects of stretching. This is because, by its formulation, the Batchelor scale is the smallest scale at which such features can develop without being dissipated by the smoothing effects of molecular diffusion. Provided that HiPS likewise resolves the Batchelor scale, then in the context of mixing studies of diffusive advected scalars, theoretical concerns about purely advective stretching effects are of no more relevance than in DNS.

4.2.2 Flow and Scalar Sub-ranges

Much of the phenomenology of scalar mixing in turbulence is reflected in the scalar spectral array, i.e., the family of power spectra of scalar fluctuations parameterized by Reynolds number Re and Prandtl number Pr or Schmidt number Sc. For present purposes, Pr and Sc are mathematically equivalent, the distinction between them being that Pr and Sc refer to the diffusion of heat and species respectively. Following convention, Pr is used here generically to designate either.

A suitable starting point for further investigation of the HiPS representation of turbulent mixing is to determine the degree to which it reproduces the features of the scalar spectral array, which raises the question of how a mathematical analog of the scalar power spectrum can be defined in the hierarchical geometry. A property of the scalar power spectrum E(k) that suggests a suitable analogy is that the k dependence of \(\int_{k}^{\infty} E(k') \, dk'\) scales as the variance of the scalar fluctuations associated with wavenumbers exceeding k. Operationally, this variance can be defined as 〈(θ−θ k )2〉, where θ is the scalar field and θ k denotes the filtered scalar field, such that fluctuations at wavenumbers greater than k have been removed. A physical-space moving average of width w=2π/k is an example of such a filter (albeit one that is not spectrally sharp).

In HiPS, w can be associated with the level n at which w=l n , recognizing that the correspondence is rough due to the step-wise discretization of the length scales l n . The moving average is continuous in physical space, but in HiPS only discretized averaging is meaningful, so in HiPS, length scales smaller than w are smoothed by averaging the θ values of the fluid parcels within each level-n sub-tree, which is analogous to physical-space box filtering. (Now, the fluid-parcel contents are scalar values θ i referenced by the node index i defined in Sect. 2.) Thus, the HiPS analog of \(\int_{k}^{\infty} E(k') \, dk'\) is \(\operatorname{var}_{n} \theta \), where \(\operatorname{var}_{n} \theta \) denotes the average over level-n sub-trees of the variance of θ within each sub-tree.

This spectrum analysis is applicable only if sub-trees are statistically indistinguishable for any given n, which is the HiPS analog of spatial homogeneity. Assuming that homogeneous turbulence is represented operationally using finite trees, the formal definition of \(\operatorname{var}_{n} \theta \) in HiPS additionally invokes ensemble averaging. Alternatively, an approach involving a chain of finite trees and a way of relating ultrametric distance to physical distance [7] can be used.

Recalling that A<1, the HiPS analog of \(\int_{Ak}^{k} E(k') \, dk'\) is then \(\operatorname{var}_{n-1} \theta - \operatorname{var}_{n} \theta \). Approximating the integrand by E(k), this gives \(E(k) \sim k^{-1} (\operatorname{var}_{n-1} \theta - \operatorname{var}_{n} \theta )\), which is the desired HiPS analog of the scalar power spectrum (omitting numerical coefficients). This also shows that the HiPS analog of θ 2(l), which is the scalar variance contribution from an order-one geometric wavenumber interval in the vicinity of k=2π/l, is \(\operatorname{var}_{n-1} \theta - \operatorname{var}_{n} \theta\).

The physical-space analog of this procedure would use the same procedure to compute E(k), but using physical-space box filtering with box size l instead of averaging of θ within scale-l sub-trees. In physical space, box filtering is unsuitable because it removes all fluctuations within the box, not solely those of wavenumber k>2π/l. The problem is that a systematic θ variation across the box can be a portion of a trend that extends over much larger scales. In contrast, moving-average filtering, which preserves such trends, is a suitable low-pass filter.

In HiPS, such large-scale trends are reflected in statistical differences between sub-trees rather than in trends within sub-trees. In other words, such trends are themselves hierarchically partitioned. For example, consider the HiPS scalar mixing layer. There is always a systematic difference between two level-1 sub-trees as the system evolves, but the only manifestation of this difference is in the difference between the statistics of the two level-1 sub-trees. For the given initial condition, there is no statistical difference between higher-level sub-trees within a level-1 sub-tree. Because the higher-level sub-tree states have no low-wavenumber (k≪2π/l) content available for removal by filtering, the sub-tree averaging removes only high-wavenumber content, and the average then reflects the order 2π/l wavenumber content. On this basis, the HiPS analog of box filtering is functionally a low-pass filter and therefore is suitable for its intended use.

Evaluation of the scalar power spectrum (now not distinguishing between the true power spectrum and its HiPS analog where the meaning is clear) requires the determination of the n dependence of \(\operatorname{var}_{n} \theta\). To proceed further, the parameters of the scalar spectral array and some elementary relations among them are explained.

In Sects. 2 and 3.1 there are oblique references to viscous dissipation. Now the role of viscosity is considered explicitly. It governs the diffusive transport of momentum, and if HiPS included momentum as a fluid-parcel property, it would be time advanced, including diffusive exchanges, in a manner roughly like the HiPS scalar treatment. (Such an extension is described in Sect. 5.1.) Instead, flow effects on scalar advancement are parameterized through the assigned n dependence of τ.

As noted in Sect. 2, there is a length scale η and an associated time scale τ η (the Kolmogorov length and time scales) at which the turbulent cascade transitions from the inertial range to the dissipative range. Dimensionally, they are based on the product of unique powers of ϵ and the kinematic viscosity ν that yield a length and a time respectively. On this basis, η=(ν 3/ϵ)1/4 and τ η =(ν/ϵ)1/2. The same results are obtained heuristically by identifying the length scale l=η at which the advective time scale τ∼l 2/3/ϵ 1/3 equals the time l 2/ν for viscosity to smooth size-l flow structures. At l<η, viscosity suppresses the turbulent breakdown process, so the flow is smooth with associated strain 1/τ η .

Turbulence intensity is quantified using the Reynolds number \(\operatorname{Re} = \nu_{t} / \nu\), where the turbulent viscosity is expressed in terms of the length and time scales of turbulence energy input as \(\nu_{t} \sim l_{0}^{2}/t_{0}\). It follows that \(\eta = \operatorname{Re}^{-3/4} l_{0}\).

The family of scalar power spectra for given Re is parameterized by \(\operatorname{Pr} = \nu / \kappa\), where κ is the thermal diffusivity, or equivalently for present purposes, the species diffusivity. Analogous to the Kolmogorov scale, the Batchelor scale l B is the length scale at which the advective time scale equals the time \(l^{2} / \kappa = \operatorname{Pr} l^{2} / \nu\) for diffusivity to smooth size-l scalar variations. l B marks a transition between scalar spectral ranges and thus is a key signature of turbulent mixing phenomenology.

Different considerations apply for \(\operatorname{Pr} > 1\) and \(\operatorname{Pr} < 1\), so they are considered separately. For \(\operatorname{Pr} > 1\), l B <η. In the scale range l B <l<η, the diffusion time \(\operatorname{Pr} l^{2} / \nu\) is balanced against the advective time τ η that governs l<η. This gives \(l_{B} \sim \operatorname{Pr}^{-1/2} \eta\). For \(\operatorname{Pr} < 1\), l B is in the inertial range, so \(\operatorname{Pr} l^{2} / \nu\) is balanced against the inertial-range l dependence of τ to obtain \(l_{B} \sim \operatorname{Pr}^{-3/4} \eta\).

For non-unity Pr, η and l B divide scale space into three ranges. The HiPS representation of mixing phenomenology and scalar spectral scalings in these ranges is now considered.

For l exceeding both l B and η, inertial-range scale breakdown occurs and molecular diffusion has negligible influence. For this ‘inertial-convective’ range, the conventional (and extensively validated) picture of scalar length-scale breakdown is closely analogous to the TKE breakdown picture described in Sect. 3.1. Namely, the rate χ of dissipation of scalar variance governs the rate of down-scale transfer of scalar fluctuations, and the quasi-stationary, non-dissipative nature of this transfer requires χ to be independent of l. From χ∼θ 2(l)/τ and the inertial-range scaling τ∼l 2/3, it follows that the scalar fluctuation amplitude θ(l) associated with length scale l scales as θ(l)∼l 1/3. This is the same as the inertial-range scaling of v(l) because both results follow from the same physical picture.

In HiPS, that physical picture justifies the assignment of the scaling τ∼l 2/3, but scalar cascading in HiPS is an outcome rather than a prescribed behavior, so the nature of that cascading in HiPS must be ascertained. The relevant features of HiPS have already been noted. First, the swapping operation induces a length-scale breakdown that is local in scale space (Sect. 3.3). Second, the process is quasi-stationary because the governing time scale τ is assigned to be a decreasing function of hierarchy level n, and the small scales therefore respond quickly to the relatively slow time variation of input from large scales. Third, the scale range under consideration is non-dissipative with respect to the scalar field. Fourth, \(\operatorname{var}_{n-1} \theta - \operatorname{var}_{n} \theta \), which is the adopted measure of level-n scalar fluctuations, has been shown to be a valid analog of θ 2(l). Based on these considerations, the heuristics leading to θ(l)∼l 1/3 are applicable to HiPS, in particular implying the inertial-convective power-spectrum scaling E(k)∼k −5/3.

Intermittency can cause scaling exponents to deviate from values based on heuristic estimates. The nature of any such deviations may depend on details of the HiPS formulation that affect its intermittency properties, a topic that will be addressed in future work.

As noted in Sect. 4.2.1, HiPS behaves unphysically if it resolves scales much smaller than l B . Therefore the only additional scale range of interest is the intermediate range l B <l<η for \(\operatorname{Pr} > 1\). In this ‘viscous-convective’ range, scalar fluctuations are not dissipated, and all the features of inertial-range cascading apply except that τ is independent of l rather than scaling as l 2/3. Therefore the l dependence of the rate θ 2(l)/τ of scalar down-scaling is given by θ 2(l). Owing to the lack of molecular-dissipative influence at scales larger than l B , this must be independent of l. The earlier discussion of the scalar power spectrum showed that E(k)∼k −1 θ 2(l), so the Batchelor scaling [3] E(k)∼k −1 is obtained.

Batchelor’s analysis is more elaborate than this because it addresses additional features such as the functional form of the spectrum transition to the dissipative range l<l B . Though the present derivation is with reference to HiPS, the noted analogies to turbulent mixing phenomenology indicate that it is more generally applicable.

The artifact of HiPS when scales l<l B are resolved reflects the restriction of diffusive exchange to the smallest resolved scale l N , at which the fluid parcels reside. This excludes direct molecular-diffusive exchange at parcel separations in the range l N <l<l B over which it is physically the dominant transport mechanism. A possible mitigation of this artifact would be to allow molecular-diffusive exchange among all the parcels in any sub-tree whose advective time τ is larger than its molecular-diffusive time scale. Allowing such interaction among more than two parcels introduces additional modeling options, as discussed in a different context in Sect. 5.2. If the suggested mitigation is successful, it would be interesting to use HiPS to study the low-Pr ‘inertial-diffusive’ scalar spectral regime η<l<l B due to subtle considerations that raise some uncertainty about its scaling behavior [12].

4.2.3 Fractal Dimension of Iso-surfaces

The definition of proximity in Sect. 3.3 is shown there to be a valid analog of physical-space proximity in the context of parcel dispersion. The HiPS analog of the scalar power spectrum (Sect. 4.2.2) is likewise based on the association of sub-tree levels with physical-space proximity, thus further extending the analogy. This analogy is also suitable for defining and analyzing the HiPS analog of the fractal scaling of iso-surfaces of the scalar field θ.

Specifically, the fractal dimension D f of these iso-surfaces is defined operationally using the box-counting method. The level-n θ ∗ iso-surface area a n is defined as \(m_{n} l_{n}^{2}\), where m n is the number of level-n sub-trees that contain at least one parcel with θ>θ ∗ and one with θ<θ ∗. According to the definition of D f , \(a_{n} \sim l_{n}^{2-D_{f}}\), which gives \(a_{n} \sim A^{(2-D_{f})n}\).

It was noted in Sect. 2 that the requirement that the number of parcels in a level-n sub-tree scales as \(l_{n}^{3}\), so that sub-trees obey 3D volumetric scaling based on the number of parcels they contain, implies A=2−1/3. This gives \(a_{n} \sim 2^{(D_{f}-2)n/3}\). As an illustration, suppose that the iso-surface is space filling, corresponding to D f =3. Then all sub-trees intersect the iso-surface so m n ∼2n, giving \(a_{n} \sim m_{n} l_{n}^{2} \sim 2^{n} A^{2n} \sim 2^{n/3}\), which agrees with the specialization of the scaling expression for a n to D f =3.

The adopted definition of sub-tree occupancy by the scalar iso-surface introduces the following anomaly. In terms of parcels contained by a sub-tree, each sub-tree is itself the union of two sub-trees. If one of these sub-trees has only parcels with θ>θ ∗ and the other has only parcels with θ<θ ∗, then neither is occupied by the scalar iso-surface but the sub-tree containing these two sub-trees is occupied by the scalar iso-surface, a situation that does not occur when the box-counting method is applied to a spatially continuous scalar field. This anomaly reflects the restricted definition of parcel adjacency in the tree geometry relative to physical space.

In this regard, the fractal properties of HiPS should be interpreted as an abstract representation of the geometrical structure of the iso-surface rather than a complete mathematical analog of physical-space iso-surface geometry. Analogous considerations apply generally to hierarchical models of the turbulent cascade [22].

Thus far it has been assumed that the fluid parcels contained in the tree are intended to represent the fluid contents of a volume of physical space, giving A=2−1/3. As noted in Sect. 2, it is no less valid to view the collection of parcels as a representation of fluid states along a low-dimensional cut through 3D space. For a d-dimensional cut, with d=3 corresponding to the volumetric interpretation adopted thus far, the determination of A in Sect. 2 straightforwardly generalizes to give A=2−1/d. For general Euclidean dimension d, the scaling of a fractal subset generalizes to \(a_{n} \sim l_{n}^{d-1-D_{f}}\) in present notation. Substitution of l n =A n l 0 now gives \(a_{n} \sim 2^{(D_{f}-d+1)n/d}\). This expression allows D f to be determined from the quantities a n for any Euclidean dimension d implied by the choice of A. The d=3 definition of the HiPS scalar iso-surface generalizes directly to the scalar isopleth in any d-dimensional cut through 3D space, so a n is operationally defined in HiPS for any such isopleth. Moreover, the additive law D f (3)=D f (d)−d+3 expresses the scalar isopleth fractal dimension D f (3) in the entire 3D space in terms of the fractal dimension D f (d) of the isopleth in any d-dimensional cut [21, 22]. Thus, for any choice of A and associated d value, there is a meaningful definition of the fractal dimension of scalar isopleths that can be related quantitatively to observed and predicted 3D physical values.

5 Model Extensions

5.1 Flow Simulation Using HiPS

Parcel swapping is the HiPS representation of the advective term in the scalar advection-diffusion equation. The momentum equation governing flow-field evolution has roughly the same form, suggesting that the flow itself might in some manner be simulated using HiPS. Indeed, a minimal flow representation in which velocity is a scalar v treated in the same manner as θ, but with κ replaced by the kinematic viscosity ν, might be useful, as well as possible extensions incorporating vector velocity and effects of the pressure-gradient term in the momentum equation.

Here, the formulation of a minimal flow treatment to make the HiPS mixing representation more self-contained and physically complete is described. The adjunct to HiPS that enables this is the use of the flow state to evaluate the τ value at each node at each instant. In other words, τ becomes both node-dependent within a level and time-varying at each node.

The HiPS analog of the scalar variance at scale l is defined in Sect. 4.2.2, thus identifying a scale-dependent fluctuation θ(l) whose close connection to the notional scale-dependent velocity fluctuation v(l) has been noted. One way to arrive at a single representative time scale applicable to all level-l nodes is to substitute an expression for v(l) into the relation τ(l)∼l/v(l), in close correspondence to the reasoning in Sect. 3.1.

However, a richer physics representation is obtained by applying this approach locally, on an individual sub-tree basis, rather than averaging over sub-trees at a given level. This is done by constructing a node velocity v i for each node i from the set of parcel velocities within sub-tree i.

Reflecting the notion that v i , like v(l), represents a fluctuation at the level of node i with higher-level fluctuations filtered out, the first step is to compute the sub-tree i mean velocity \(\bar{v}_{i}\) at each node i. At level N, where parcels reside, \(\bar{v}_{i}\) is simply the parcel velocity at node i.

Then a definition of v i that is suitable for many cases is \(v_{i} = |\bar{v}_{i0} - \bar{v}_{i1}|\). However, for reasons that will be explained elsewhere, this definition is inadequate in some instances, so instead the more robust formulation \(v_{i} = |\max ( \bar{v}_{i0}, \bar{v}_{i1}, \bar{v}_{i00}, \bar{v}_{i01}, \bar{v}_{i10}, \bar{v}_{i11}) - \min ( \bar{v}_{i0}, \bar{v}_{i1}, \bar{v}_{i00}, \bar{v}_{i01}, \bar{v}_{i10}, \bar{v}_{i11}) |\) is used. The intent is to equate v i to the magnitude of the largest mean-velocity difference two levels below i within the sub-tree. \(\bar{v}_{i0}\) and \(\bar{v}_{i1}\) are superfluous at levels n<N−1 because their omission would not affect the result; they are never outside the range of the other arguments of the max and min functions. They are included only because evaluation of v i at level N−1 would otherwise require quantities at level N+1, which does not exist. In this case, i0 and i1 are level-N nodes at which parcels reside, so as specified in the preceding paragraph, \(\bar{v}_{i0}\) and \(\bar{v}_{i1}\) are assigned to be the parcel values at those nodes. Thus, a local analog of v(l) is obtained, giving finally the local analog τ i =l/v i of τ(l).

This reasoning is valid in the inertial range, in which advection at a given scale is governed by flow structures of comparable size. As noted in Sect. 2, the relevant time scale in the dissipative range is the one associated with the Kolmogorov microscale η. A practical way to implement this on a local basis is to assign τ i to be the minimum of the value computed locally, as above, and smallest value assigned to any node on the path from the apex of the tree to node i. (To avoid recursive evaluations, this requires top-to-bottom updates of τ i values.) This operationally defines η locally within the flow simulation as the level n along any given path at which l n /v i is smallest. This approach has the added benefit that inertial-range fluctuations producing unusually small τ i values at large scales directly modify time scales over some sub-range of smaller inertial-range scales, reflecting scale non-locality associated with large-scale extreme events.

These observations implicitly assume that basic inertial-range properties such as −5/3 power-law scaling of the energy spectrum (defined analogously to the scalar power spectrum definition in Sect. 4.2.2) will be reproduced by the flow simulation. Such behavior would be consistent with the key features of the HiPS cascade, namely the non-dissipative nature of swapping (because it re-indexes parcels without affecting their contents) and the step-wise, hence local in scale space, mechanism of length-scale breakdown.

This indicates the plausibility of inertial-range behavior in HiPS, but it does not constitute a demonstration that the flow simulation is uniquely constrained to behave in this manner. This can be definitively established only by numerical implementation, but some assurance is provided by comparison of the flow formulation to one-dimensional turbulence (ODT). ODT is formulated in one spatial dimension and represents advection by means of instantaneous maps that, in effect, swap groups of mesh cells in a manner that multiplicatively reduces the length scales of property fluctuations. In its original formulation [17], the time scale governing map occurrences is defined locally in space and time in a manner much like the approach described here. In this and other respects, it is roughly accurate to describe HiPS flow simulation as an implementation of ODT in a hierarchical geometry. ODT exhibits spectral and other properties of inertial-range turbulence, supporting the expectation that HiPS flow simulation will do likewise.

The relationship to ODT suggests that algorithmic features of ODT will also be pertinent. For example, the thinning method [19] that is used to sample ODT map occurrences will likewise be advantageous for sampling swap occurrences.

ODT was preceded the Linear Eddy Model (LEM) [16], in which mapping frequency is a specified function of map size, much as τ has a specified dependence on sub-tree level when the flow simulation is not used to determine τ locally. The HiPS formulation and analysis in Sects. 2 through 4 in many ways reflects analogous considerations pertinent to LEM.

The flow simulation will induce local fluctuations of τ that imply local transient variations of the length scales η and l B at which there is a balance of advective and, respectively, viscous and molecular-diffusive time scales. This further complicates the considerations in Sect. 4.2.1 determining the preferred number of tree levels N. Local adaptivity of the tree termination level might prove to be especially advantageous for HiPS flow simulation.

As noted in Sect. 4.2.3, the freedom to adjust a particular manifestation of dimensionality through the choice of A allows the possibility of formulating the HiPS flow simulation on a 1D or 2D sub-space through the choice A=2−1/d for a d-dimensional sub-space. In such sub-spaces, initial-value problems representing several canonical free-shear flows can be formulated in a manner that captures the similarity scalings of real flows obeying the corresponding symmetries.

The use of A to set the sub-space dimension and thus the flow symmetry is illustrated by the HiPS initial-value problems representing temporal free-shear flows. In physical space, an ideal temporal free-shear flow is a flow whose statistical properties do not vary in the streamwise direction. It can be initialized by setting fluid in mean streamwise motion, and possibly in zero-mean motion in other directions, such that the flow statistics vary only in the plane normal to the streamwise direction. For the round jet in quiescent surroundings, this initial motion is confined to a circular cross-section and for the planar jet it is confined to a slab or slot.

Each of these flows can be described statistically as a function of a single coordinate, radial and lateral respectively. However, their effective dimensions are different in the following sense. Lateral spreading of the planar jet causes property dilution corresponding to an average change of scalar values that is inversely proportional to the degree of spreading, and thus corresponding to an effective dimension d=1. However, the corresponding dilution effect of radial spreading of the round jet causes an average change of scalar values that is inversely proportional to the square of the radial spreading, corresponding to an effective dimension d=2.

This distinction is reflected in the similarity scalings of temporal turbulent planar and round jets, illustrated by constant-density momentum scaling. For the planar jet, assume characteristic velocity V(t) and characteristic width W(t) at time t. Then global momentum conservation implies constant VW, hence V∼1/W. The jet Reynolds number VW/ν is thus constant, and the eddy-diffusive scaling dW 2/dt∼VW gives W∼t 1/2 because VW is constant. For the round jet, W is replaced by characteristic jet radius R(t) and momentum conservation implies constant VR 2, hence V∼1/R 2. Then eddy-diffusive spread is governed by dR 2/dt∼VR∼1/R, so R∼t 1/3. Now, \(\operatorname{Re} \sim 1/R\) and V∼t −2/3.

In HiPS flow simulation, a simple initial condition that develops into a turbulent jet is obtained by setting one parcel velocity equal to some value V and the others to zero, where V is large enough so Re corresponding to developed turbulence will persist through some desired extent of the simulation. The dilution scaling corresponding to planar (cylindrical) symmetry is enforced by choosing d=1 (d=2) in A=2−1/d. In this manner, the similarity scalings of the respective cases are enforced.

In HiPS jet simulations, larger swaps are inherently less frequent because their time scales increase at large scales due to the factor l in the expression τ i ∼l/v i . This does not assure that their effects are small because turbulent diffusivity scales as l 2/τ i ∼lv i , which is negligible at large l only if v i decreases rapidly enough with increasing l.

For the temporal jets, the l dependence of v i for l≫L, where L generically denotes W or R, is estimated as follows. Consider a sub-tree corresponding to some scale l≫L that contains a jet of lateral span L. Based on the definition of v i , it scales as the difference of mean velocities of two sub-trees of order-l scale. Typically, only one of these sub-trees intersects the jet and therefore fully contains it, so its mean velocity scales as (L/l)d V, and the other sub-tree has mean velocity zero, giving v i ∼l −d.

For the planar jet, d=1 implies that swaps at all scales contribute comparably to the turbulent diffusivity lv i , so the turbulent diffusivity diverges when it is summed over scales. An analogous anomaly arises in ODT and is mitigated by any of several possible ‘large-eddy-suppression’ methods [1, 20] that are likewise applicable to HiPS. For the round jet (d=2), there is no such divergence anomaly, but even the occasional occurrence of swaps at scales l≫L might introduce unphysical artifacts that can also be mitigated using these methods.

Based on the noted considerations, the role of dimensionality in HiPS is fundamentally different than in other hierarchical approaches such as [2], in which details of the formulation affect whether the model exhibits 2D or 3D cascading. In all instances, HiPS represents 3D turbulence, but with the capability to choose the dimension of the sub-space of the flow domain that is used to simulate flow advancement.

5.2 p-Adic HiPS

Specialization to a binary (dyadic) tree is convenient for introducing HiPS but is not necessarily advantageous relative to p-adic trees with branching factors p>2. The HiPS formalism straightforwardly generalizes to arbitrary p. Now the node index i is a base-p number that uniquely specifies the path to node i and A=p −1/3 satisfies the requirement that the number of parcels at the base of the tree scales as \(l_{N}^{-3}\) for any p. For the interpretation of HiPS as a d-dimensional cut through 3D space, this generalizes to A=p −1/d.

For p>2, a new consideration is that parcel pairings for swapping and mixing occurrences are no longer uniquely specified by the tree geometry. One way to resolve this ambiguity is to introduce some form of random sampling to determine the pairings. At the base of the tree, this enables various methods of representing local structures such as flames within HiPS. For example, mixing can be restricted to pairs whose index values are adjacent (i.e., differ by unity). Also, swapping and mixing are no longer required to be pairwise. Swapping could become a permutation of any number m≤2p−2 of sub-trees based on random or index-conditioned selection of the affected sub-trees.

Likewise, mixing could involve an all-to-all coupling such as the interaction by exchange with the mean (IEM) model [10]. Because mixing occurs at the base of the tree and involves distinct modeling considerations relative to swapping, it might be advantageous to choose a different branching factor at the base than at other levels. The preferred choice among these and other possibilities is likely to be case-specific.

6 Discussion

Various approaches to modeling turbulence and turbulent mixing by coupling state variables defined on a set of spatial (or equivalently, wavenumber) levels within a geometrical hierarchy have been developed [2, 4, 5, 7, 10, 13, 14, 22]. A distinguishing feature of the present hierarchical approach is that the physical state is specified solely at the smallest scale and is thus fully resolved in space and time. Quantities defined at coarser scales are auxiliary variables used to specify the time advancement of the resolved parcels, whose contents constitute the physical state.

Effects of inertial-range phenomenology, especially small-scale intermittency, on turbulent reacting flows have been addressed in numerous modeling studies, and opportunities for further progress have been noted [18, 25]. The approach introduced here enables the incorporation of these effects into an important class of turbulent combustion models in greater detail than previously. Intermittency per se has not been discussed in detail here, but the features of HiPS that are analogous to other models (e.g., hierarchical and LEM) that provide useful intermittency representations suggest that HiPS will perform similarly in this regard.

In Sect. 5.1 it is shown that symmetry considerations constrain the spatial dimension d implied by the relation A=2−1/d to be either 1 or 2 for certain self-contained HiPS applications. In subgrid and other implementations that are not intended to capture large-scale flow inhomogeneity, no such considerations apply, so in principle any desired d can be used, including d<1, allowing arbitrary reduction of the number of parcels required to simulate a given configuration. Some commensurate loss of fidelity can be anticipated, but for d corresponding to the number of parcels used in the mixing models described in Sect. 1, HiPS introduces relevant turbulence phenomenology with little added cost, so cost/performance benefits are anticipated.

Though the proximate motivation of this formulation is the improvement of such mixing models, it has been shown that the approach is of potential interest more broadly. In some respects it is a bridge between fundamental statistical physics concepts and methodology and the needs of engineering applications. To emphasize this point, an overview of the approach has been presented that encompasses both realms and highlights the connections. The formulation has been outlined and features that can be inferred from analysis or heuristics have been described. Sequels to this overview will be narrower in scope, with greater focus on numerical implementation, including details of its use as a subgrid combustion closure. This requires fleshing out the details of HiPS to obtain a fully specified algorithm, in fact, several due to application-specific considerations that have been noted. Here, HiPS has been termed a representation rather than a model of turbulent mixing to emphasize that it encompasses various possible model instantiations.