{'fraction_non_alphanumeric': 0.02524364914523087, 'fraction_numerical': 0.00556736760603193, 'mean_word_length': 5.009009674322428, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 1}, '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}}
Large Eddy Simulation of Flow Interactions Between a Turbulent Free-Stream and a Permeable Bed

Systems comprising a turbulent channel flow overlaying a permeable bed can be found in a variety of industrial and natural applications (e.g. urban planning, fracking, submerged vegetation). One important realization of this system is at the bottom of rivers, where surface waters of the river permeate the deposited sediment and exchange their contents with subsurface waters. Obtaining a complete picture of this process, known as hyporheic exchange, is important for the understanding and maintenance of water quality and river ecology. Given the wide range of length scales (mm to km) and corresponding time scales of interest, reduced-order numerical models (e.g. stochastic particle tracking, advection-dispersion equations) are often used to study the transport of mass and momentum in this system. However, the many scales associated with surface/subsurface mixing make the development of a high-fidelity model for turbulent scalar transport a formidable challenge. In an effort to capture transport behavior at the pore scale, this work is focused on development of a detailed, three-dimensional flow model for the Large Eddy Simulation (LES) of turbulent flow over a permeable bed. A rigorous assessment of the computational mesh, domain size and turbulence model is performed. Concurrently, the double-averaging methodology is used to study the large-scale, persistent flow behavior in the presence of spatio-temporal heterogeneity. It is found that accurate resolution of the largest turbulent scales is necessary to develop a complete picture of interfacial transport, and that these large structures are felt deeper within the bed. Additionally, the development of a particle tracking software module is detailed. This module provides a foundation for the extraction of Lagrangian dispersion information from the LES, which will eventually be used to physically inform reduced-order models.


Abstract
Systems comprising a turbulent channel flow overlaying a permeable bed can be found in a variety of industrial and natural applications (e.g. urban planning, fracking, submerged vegetation). One important realization of this system is at the bottom of rivers, where surface waters of the river permeate the deposited sediment and exchange their contents (e.g. nutrients) with subsurface waters. Obtaining a complete picture of this process, known as hyporheic exchange, is important for the understanding and maintenance of water quality and river ecology [1].
Given the wide range of length scales (mm to km) and corresponding time scales of interest, reduced-order numerical models (e.g. stochastic particle tracking, advection-dispersion equations) are often used to study the transport of mass and momentum in this system [2,3].
However, the many scales associated with surface/subsurface mixing make the development of a high-fidelity model for turbulent scalar transport a formidable challenge [4].
In an effort to capture transport behavior at the pore scale, this work is focused on development of a detailed, three-dimensional flow model for the Large Eddy Simulation (LES) of turbulent flow over a permeable bed. A rigorous assessment of the computational mesh, domain size and turbulence model is performed. Concurrently, the double-averaging methodology is used to study the large-scale, persistent flow behavior in the presence of spatio-temporal heterogeneity. It is found that accurate resolution of the largest turbulent scales is necessary to develop a complete picture of interfacial momentum transport, and that these large structures are felt deeper within the bed.
Additionally, the development of a particle tracking software module is detailed. This module provides a foundation for the extraction of Lagrangian dispersion information from the LES, which will eventually be used to physically inform reduced-order models.

Background and Motivation
The work in this thesis is focused on the numerical study of systems in which a rigid permeable bed is bounded on one side by an unobstructed surface flow. Here, a permeable bed is defined as some structure containing adjacent, connected pores throughout which fluid may flow, and the term rigid is used to specify that the solid constituents of this permeable bed do not move in time. Such a system is the subject of study in myriad environmental and engineering applications, including submerged vegetation [6], urban planning [7] and hydraulic fracturing [8].
Another application, which has proved to be a fundamental influence of the work presented in this thesis, appears in the study of fluvial systems. Throughout rivers, mass and momentum are exchanged between the surface waters of the river channel and the groundwater, or subsurface flow, when fluid permeates the bedforms and deposed sediment at the bottom of the channel. Quantifying the complex interactions between the surface and subsurface flows and their effects on scalar transport within this region, referred to as the hyporheic zone and shown in Figure 1.1, is critical to the understanding and study of river ecology and water quality [1]. Boano et al. recently authored a thorough review of physical mechanisms, numerical models and environmental implications associated with hyporheic exchange [4].
A dominant theme of this work is that the many temporal and spatial scales associated with phenomena of interest (e.g. nutrient delivery to sediments) make the development of a detailed numerical model which accounts for the broad range of physical activity seen in this region challenging.
Profiles of concentration break through curves, which generally exhibit heavy tailing, have been modeled with some success using stochastic particle tracking models and fractional advection-dispersion equations (fADE) [2,9,3]. Additionally, there has been work using the results of particle tracking simulations to parameterize lower-order models, known as upscaling [10]. However, such numerical models often must be parameterized ad hoc and lack physical information at the individual pore scale, where geometry-induced flow structure (e.g. vortices) are thought to be a primary cause of the heavy tailing seen in experiments [11].
Significant effort has been allocated towards both understanding and quantifying the interactions between surface and subsurface waters, and more broadly the interfacial momentum transport, at the pore scale. An illustration of the driving processes behind momentum transport at the surface/subsurface interface is shown in Figure 1.2. Interest is often directed at the length over which turbulent structures penetrate into the bed, δ e , and the depth-wise variation in an averaged streamwise velocity profile, u(z). This average is frequently taken to be a time-space average, rigorously defined in [12], which eases the study of large-scale, persistent structure in such a heterogeneous environment. While in a laminar flow regime δ e has been found to be on the order of a single grain diameter and insensitive to the bulk Reynolds number, derived from the mean surface velocity and surface flow height [13], the introduction of free-stream turbulence adds complexity to the structure and development of these features.
Experiments have shown that the time-space averaged velocity is found to decrease dramatically from the surface into the transition region, achieving a minimum in the first layer of pores, then increasing towards a constant velocity deep within the bed [14,15]. Additionally, the friction factor, or ratio of frictional to inertial forces, in flows over permeable beds has been found to increase with Re far beyond the plateau seen in flows over rough walls [16], and is suspected to be a result of increased turbulent penetration into the bed as turbulent structures grow in size [5]. This penetration significantly alters the near-wall, or more truly the near-interface, turbulent structure from what is seen in rough wall flows, prompting the construction of several modified laws to better describe velocity scaling in the logarithmic layer of the flow. It also raises concern that the viscous wall unit may not be appropriate for scaling inner variables in such a system [17,18,19].
In addition to these experiments, simulations have proved to be an invaluable tool for Figure 1.2: A conceptual schematic of momentum transport in a permeable bed bounded by a free surface flow. Eddies (curled arrows) of different scales are seen moving fluid in and out of the the transition region (a layer of length δ e experiencing high turbulent activity) with a penetration velocity v p . The average streamwise velocity changes from a constant surface magnitude V to a constant Darcy velocity U d within the bed. Taken from [5].
investigating the driving processes behind surface/subsurface flow coupling. Direct numerical simulation (DNS) of turbulent flow over a channel of regularly packed cubes confirmed that large vortical turbulent structures near the permeable wall are generated from Kelvin-Helmholtz (KH) type instabilities, which significantly contribute to the skin friction felt by the surface flow [19]. The same study, in concert with previous work [20], also promotes the use of the permeability Reynolds number, Re K , and the roughness Reynolds number, Re D , as means to quantify the influence of bed permeability and roughness elements on the flow field, respectively. DNS has been used to compare flow over permeable and impermeable beds of staggered cubes, also confirming the prevalence of interfacial KH instabilities.
Although, a significant contribution from sweep events (i.e. downward fluid motions with high streamwise momentum) to the Reynolds shear stress occurs just below the interface, the Reynolds shear stress decays almost entirely over the first grain diameter [21,22]. The influence of these hydrodynamic characteristics on heat transfer has recently been studied in a similar packed cube geometry, finding enhanced turbulent heat transfer above the porous interface, due to the presence of large eddies, and large temperature fluctuations deep within the bed, caused by pressure waves propagating from the transition region [23].
While tremendous progress has been been made in understanding these systems, challenges moving forward are evident. Experimental investigation has large, physical constraints on its ability to capture spatial and temporal information. Although simulations may compliment experimental work, having a much higher spatio-temporal resolution, the computational resources necessary to perform DNS limit the system of interest to low Reynolds numbers and relatively simple geometries. One alternative to DNS, which offers a high degree of detail in the resolved flow field without imposing the same intractability problems of DNS, is Large Eddy Simulation (LES) [24,25].
The driving principle behind LES is to apply a low-pass filter to the Naiver-Stokes equations and model the effect of the unresolved motions with some closure scheme. Given that the fine scale motions are not directly computed, the case for LES as a useful tool grows stronger when momentum and mass transport is driven by the larger, resolved flow structures [26]. This is certainly suggested by the previously mentioned literature for the case of flows over permeable beds. Thus, LES seems to be an reasonable candidate tool for investigating a system so rich in structure and host to various phenomena (e.g. flow separation and reattachment, surface/subsurface mixing).
There has been limited use of LES to examine flows over fixed, coarse-gravel beds. Notably, Stoesser et al. simulated turbulent flow over a three-layer deep bed of spheres, providing evidence that LES is able to predict several experimental and DNS results for this type of system (e.g. streamwise velocity profiles and penetration of pressure waves deep into the bed) [27]. However, this work does not address the impact of many numerical modeling decisions (e.g. turbulence model, domain size) on the resultant flow statistics. Although LES is generally tractable, yet expensive to run, the cost associated with rigorously testing the numerical model may be orders of magnitude more than a single simulation with respect to compute time [26]. Regardless of this expense, several studies have cemented the need to examine variations in the flow field due to changes in turbulence model, domain size, mesh size, etc., particularly when using periodic boundary condition in both the streamwise and spanwise directions [28,29,30].
The goal of the work presented in this thesis is to advance the understanding of the driving mechanisms behind interfacial momentum exchange in turbulent flows over permeable beds, and how these mechanisms may impact scalar transport. Primarily, progress towards this goal is achieved through the development of a suitable numerical model for a highly detailed LES of the surface/subsurface system and analysis of its predicted dynamics. Secondary to this project is the creation of a Lagrangian particle tracking model and its implementation within an open source fluid dynamics code base, with an eye on extracting and upscaling dispersion information from the detailed LES.

Outline
This thesis is the culmination of two currently distinct projects that are united under the purpose of interrogating the hydrodynamics governing hyporheic exchange. Thus, the first two chapters following this introduction may be read independently from the third, although future work will integrate these components.
Chapter 2 is concerned with developing and performing a Large Eddy Simulation (LES) of turbulent flow over a permeable bed. The characteristic geometry used in all computational experiments will be presented, followed by a discussion on the domain discretization and numerical treatment of the governing fluid equations. Finally, the double-averaging methodology for quantities of interest is introduced and initial simulation results are discussed.
Chapter 3 follows Chapter 2 chronologically, focusing on several challenges encountered while performing the LES. Particular focus is paid to the mesh refinement, issues associated with obtaining a representative volume element (RVE) and selection of turbulence closure.
Concurrent with these discussions, results from several different runs of the LES are examined to both gain insights into the processes governing momentum transport and highlight the lingering hurdles in modeling flow over permeable beds.
Disparate from the previous two chapters, Chapter 4 introduces the addition of a particle tracing module to the Nalu code base. A brief introduction to Nalu is given, followed by a thorough presentation of the particle module's organization, key algorithms and current performance.
Finally, Chapter 5 summarizes primary conclusions drawn from the numerical experiments and provides suggestions for future work on both the LES and particle module.  Each sphere has a 2 mm gap, d g , between itself and its nearest neighbor, which leads to a simple cubic packed structure for the porous medium. This gap is useful to avoid challenges associated with singularities in the geometry at contact points while remaining small enough to allow for comparing obtained simulation results against a vast body of experimental literature (e.g. [15,32,33]). To simplify discussions concerning bed geometry moving forward, a unit cell is defined to be a cube of length L = 4.0 cm which is concentric with a given sphere when in the bed.
The task of determining a suitable domain size to capture all of the relevant physics at play is not a trivial one. A primary consideration in selecting appropriate domain dimensions is the need to resolve the largest turbulent scales while maintaining a feasible computational work load. Although initial geometries presented in this chapter reflect a dominant focus on mimicking the experimental configuration mentioned above, proper resolution of turbulence will be discussed thoroughly in Section 3.3. Figure 2.1 shows several example configurations which will be examined throughout the following chapter in an effort to address this concern.
As shown in Figure 2.1, the right-handed coordinate system is adopted here, where the  Several studies (e.g. [33,27]) find that the influence of turbulence on the subsurface flow is minimal beyond a depth of two to three layers of spheres, so four layers are cautiously taken to be sufficient for resolving all dominant behavior while limiting computational cost.

Details Regarding the Computational Mesh
With an aim to elucidate the hydrodynamic processes driving the interfacial transport, an unstructured, conforming mesh consisting of both wedge and tetrahedral element types was applied to the computational domain. To mitigate the challenges associated with incorporating grains, or spherical voids, into the model (e.g. enforcing the fluid-solid interface and resolving a thin viscous sub-layer), a conforming mesh has been applied to the spherical voids, which allows one to minimize the complexity of the numerical model by directly applying no-slip boundary conditions to the spheres' surfaces, rather than applying these conditions through an immersed-boundary or interface capturing method. Additionally, several layers of wedge elements envelop the spheres in an effort to obtain high-resolution LES near the solid boundaries.
In order to determine how many layers of wedge elements are necessary to resolve the boundary layer activity while maintaining a feasible numerical work load, the viscous length scale, l * , which governs the viscous sub-layer may be computed by Equation 2.1, where τ w and u * are the wall shear-stress and shear velocity, respectively.
Literature recommends that elements nearest to the wall have a height no greater than [28,30,27]. Calculation of a meaningful wallunit relies on the use of a τ w which accurately represents the wall shear-stress experienced by fluid particles. Using a control-volume argument for a body-force driven surface flow While this formulation provides an initial guess at τ w , it is important to note that the presence of wall permeability and only tangential contact with the spheres significantly complicates the calculation of a wall-unit and obscures the result's meaning, discussed further in [17].
In the following simulations, a per volume body force of 1.1 kg/m 2 s 2 is used to drive the flow. All models have a surface-flow height of H = 0.161m, a molecular viscosity of µ = 1.002 · 10 −3 kg/(m · s) and a density of ρ = 1 · 10 3 kg/m 3 . Thus, a viscous wall-unit of l * = 7.58 × 10 −5 is obtained. Applying a four-element thick boundary layer to each sphere, where each layer grows in height by a factor of 1.2 and the initial element height is set to d g /20, the wall-nearest element height and boundary layer mesh height are z + = 1.32 and z + = 7.09, respectively. An alternative method for computing this quantity comes from taking the two wall-nearest nodes, computing a shear velocity gradient between them and obtaining a local wall shear stress using the molecular viscosity. Doing this several times across multiple datasets also recovers z + = O(1). As both methods respect the suggestions of previously mentioned literature, this boundary-layer mesh is taken to be satisfactory for resolving near-wall detail.

A Brief Discussion on Reynolds Number
When modeling flow over a permeable wall, several length scales must be taken into consid- Frequently, a bulk Reynolds number, defined by the mean flow, U b , surface height, H, and kinematic viscosity, ν, is used(e.g. [32,27,33] [19].
where K is the permeability and is the constant bed porosity. The permeability Reynolds number can be considered a comparison of the relative importance of an effective pore diameter and the length scale associated with the viscous sublayer along individual grains.
The roughness Reynolds number holds a similar meaning, substituting the height of the roughness elements for the pore diameter [19].
Given the previously defined geometric parameters, including the approximate calculation for u * , values of Re K ≈ 34 and Re D ≈ 504 are obtained. One can expect to isolate the effects of permeability if Re K > 1 and Re D << 70 [17]. As both Reynolds numbers are an order of magnitude larger than these limits, the effects of permeability and roughness should both be present and hard to distinguish from one another. Simulations with different geometric parameters will be necessary to further study the independent influences of roughness and permeability on momentum transport, though this work is beyond the scope of this thesis.

Discretization and Governing Equations
Throughout the entire domain, the fluid is assumed Newtonian and incompressible, governed by the LES model consisting of the continuity and momentum equations, respectively defined as: where u i stands for the space-filtered velocity, ρ is the fluid density and b i is the body force used to drive the flow. The deviatoric components of the resolved stress and unresolved subgrid stress are defined by τ ij and τ SGS ij , respectively. While these subgrid stresses are computed directly, the subgrid turbulent kinetic energy is taken to be an additional normal stress in the pressure term, as the majority of turbulent kinetic energy exists within the resolved scales.
The resolved shear stress is defined as τ ij = 2µS ij , where µ denotes the molecular viscosity of the fluid and S ij is the resolved rate of deformation, defined as

Turbulence Models
In this work, both the standard Smagorinsky turbulence model [24] and the Wall-Adapting Local Eddy-Viscosity (WALE) model [36] are applied to the flow model in an effort to clarify each closure's effect on the observed hydrodynamics. The constant-coefficient Smagorinsky model computes the turbulent viscosity as: where C s is set to 0.17 and |S| = 2S ij S ij . While relatively simple to implement, this sub-grid stress model is known to produce turbulent viscosity even in the viscous sub-layer of wall elements where flow is nearly laminar [28].
One alternative to the Smagorinsky model is the WALE model, specifically formulated to capture the proper scaling of turbulent viscosity in the near-wall limit, ν t ∼ z + 3 [36]. The WALE model computes the turbulent viscosity as: where the constant C w is set to 0.325. In both cases, the grid filter is computed as ∆ = V 1/3 where V denotes the volume of the local control volume.

The Double-Averaged Navier Stokes (DANS) Equations
Two decomposition methodologies will be applied to instantaneous variables obtained from the LES, both of which will aid in developing a meaningful analysis and interpretation of flow data. First, Reynolds decomposition is used to separate fields into mean and fluctuating components, which results in the Reynolds-Averaged Naiver-Stokes (RANS) equations.
Although the RANS equations provide a method for navigating temporal variations in the flow, the flow both near and within the bed exhibits a high degree of spatial heterogeneity, which makes the RANS analysis an imperfect tool. Therefore, the double decomposition methodology is also utilized. Using this technique, the time-averaged variables from the Reynolds decomposition are further decomposed into mean and fluctuating fields in space.
This analysis leads to Double-Averaged Naiver-Stokes (DANS) equations [12]. One particular benefit of the DANS analysis is that it allows for calculation of both form and viscous drags, which is critical to analysis of the subsurface flow. Each method of decomposition is introduced formally in the following subsections.

Reynolds Decomposition Methodology
For a general flow variable, θ(x, t) the Reynolds decomposition is defined as θ =θ + θ , where an overbar denotes time averaging, and the prime represents the deviation from the mean field. This temporal decomposition satisfies several well known rules (e.g. [37]), given in Equation 2.13.
where φ = φ(x, t) and s is an independent variable (i.e. x i or t). One should note that the final rule, θ · φ = θ · φ is only satisfied when the flow has sufficient scale separation [38]. Effectively, this means that the time scale over which the mean flow changes must be significantly larger than the time scale at which turbulent fluctuations persist. Given the tendency in the field to take the assumption that this criterion is met, along with the length of time over which time-averaged statistics will be developed in the following simulations, the assumption of sufficient temporal scale separation is taken here.
Applying the Reynolds decomposition to the instantaneous variables in time domain governed by Equations 2.7 and 2.8 leads to LES-based Reynolds-averaged Naiver-Stokes equations, which are referred to as RA-LES equations, as follows. and is the Reynolds stress andb i = b i when b i holds a constant value. Cross-terms produced from the decomposition (e.g. the Leonard stress) have been disregarded in accordance with the assumption that a satisfactory separation of scales exists.
Additionally, the divergence of resolved and modeled shear stresses, respectively, are computed as:

Computing a Running Time Average
Formulations for a running temporal mean must be developed for each term of interest.
The LES computes only the instantaneous fields, and maintaining the full time history for each term of interest over the course of the simulation would be prohibitive due to memory requirements. Here, an example of this derivation is given both for the mean velocity and Reynolds stress. While higher-order terms may be more complicated to compute, the underlying method for accumulating the time-averaged values is similar. In the following derivations, the superscript N , T N and ∆t N denote the time level, time at time level N and the time step size at time level N , respectively. A unit density is taken here for simplicity.
The mean fluid velocity may be expressed discretely as: Extracting the current velocity from the summation and performing some manipulation allows one to obtain the current mean velocity with compact temporal support, needing information from only the current and previous time step.
Calculation of the Reynolds stress requires a bit more work. Again, the term of interest is first presented in discrete form, assuming the time-averaged velocity has already been computed at the current time level. For the purpose of this derivation only, let τ ij = u i u j .
Now, the current velocities may be extracted from the summation, as done previously, to obtain the running mean.

Double Decomposition Methodology
When performing double decomposition, averaging fields in a different order (i.e. time-space or space-time) does not necessarily produce equivalent derivations of the DANS equations.
As noted in [12], the averaging methods commute in the specific case of a fixed bed, although time-space averaging is more suitable in general for describing rough-bed hydrodynamics, due to its consistency with traditions in turbulence research. Therefore, focus will be placed on the DANS equations derived with time-space averaging order in this work.
For our purposes, the spatial average is computed as the volume average over a thin slab, although formally the volume of interest may more general. For an arbitrary flow variable, θ, the spatial average is taken as: where the presence and lack of the subscript s denotes the superficial and intrinsic averages, respectively. These averages differ in that the superficial volume average normalizes the integration over the fluid volume by the entire averaging volume, while the intrinsic average normalizes by the fluid volume. All formulations produced here will focus on use of the intrinsic volume average, to be consistent with existing literature (e.g. [12,19]). By defining a geometry indicator function as the fluid volume fraction in the averaging domain, these averages may be related as: Similar to the traditional Reynolds decomposition, the spatial decomposition separates a variable into its spatial mean and deviation from this mean: where the θ andθ denote the mean and fluctuating components, respectively. With the spatial decomposition defined, it may be applied to a time-averaged flow variable to obtain the double decomposition:θ = θ +θ (2.26) whereθ indicates the spatial perturbation from the time-space averaged variable represented by θ [12]. One should note that the spatial average, just like the time average, is a linear operation, and thus satisfies the rule θ + φ = θ + φ .
Extending this framework to include differential operators, both in time and space, Whitaker's transport and spatial averaging theorems [39], respectively, enable one to study transport from a volume-averaged perspective. The transport theorem allows for the expansion of the temporal derivative as: where S int stands for the contact area between the fluid and the solid phase (spheres) inside the averaging volume V o , and n is the unit normal at S int that points from the solid phase into the fluid. On account of the no-slip condition applied at the surface of the spheres, v, the velocity at the surface of solid phase, is identically zero. Accordingly, the surface integral disappears and the transport equation may be simplified as: which highlights the commutability of spatial-averaging and temporal differentiation for a stationary bedform. The spatial-averaging theorem provides a similar expansion for the spatial derivative: Using Equation 2.24, relationships for the differentiation of double-averaged quantities which have been defined with respect to the intrinsic volume average are readily obtained from the above equations as: and

Dimension of the Averaging Domain
Within this work, the primary goal in taking a volume average is to capture the large-scale hydrodynamic features which arise from the presence of the bed topography, while smoothing out much of the small-scale flow structure coming from the spatial heterogeneity within the bed.
Considering the configurations shown in Figure 2.1 and the conceptual model provided in Figure 1.2, the variation in the flow field as a function of depth is of critical interest. An averaging volume enabling one to observe such behavior may be obtained by first taking a plane which covers the streamwise and spanwise dimensions in full. This surface may then be uniformly extruded a length L in the wall-normal direction to create the desired volume.
By setting L = L and taking volume averages throughout depth of the entire domain, three physical regions are recovered, distinct in their porosity.
The subscripts S, B and T are used to distinguish these surface, bed and transition regions, respectively. As shown in Figure 2.4, where the volume-averaged porosity is plotted against the depth of the averaging volume's centroid, the surface flow primarily holds a volume-averaged porosity of S = 1 and may be described as a homogeneous fluid region.
Deep within the bed, the subsurface domain is characterized by a porosity B = 0.5511, and can be referred to as the homogeneous porous region. When the centroid of an averaging volume reaches a depth z * ∈ (−0.5, 0.5), the volume enters a region of variable averaged porosity, termed the transition region. This region is described by averaged porosity T = (z), and will be of particular interest when analyzing simulation results. One should note that by taking an averaging volume of these dimensions, the three regions discussed in Useful properties are associated with data that has been averaged in this way. First and foremost, any flow variable subject to double-averaging becomes a function of only depth (i.e. θ(x, y, z, t) = θ (z). Additionally, it should be clear that as the geometry indicator takes a unit value in the homogeneous fluid region, the intrinsic and superficial average become equivalent (i.e. θ = θ s ).
Moving forward, double-averaged quantities will commonly be plotted against the depth of the centroid of the averaging volume. To simplify the description of this location, a depth z * may be defined in the following way: where

Double-Averaged LES Equations
Application of the volume-averaging operator and the theorems defined in Section 2.3 to each term in RA-LES equations (2.14 and 2.15), while enforcing the boundary conditions specified in Section 2.2, results in the double-averaged continuity and momentum equations, obtained, respectively, as: where f p i represents the form drag per unit fluid volume: and f v i denotes the viscous drag per unit fluid volume: Using 2.31, these drag terms may also be calculated in the form: and A quick note must be made concerning the volume average of the convection term in Equation 2.34. Using the double decomposition defined in Equation 2.26, the volume-averaged velocity product may be rewritten as: In line with the assumptions made for the calculation of a time average, a separation of scales is assumed to exist for the spatial average, such that θ φ = θ φ . Such assumptions allow for the simplification of Equation 2.39 as: where ũ iũj is the so-called form-induced stress [12]. It is important to note that because the difference in the length scale of the averaging volume and the extent over which geometric heterogeneity creates variation in the flow field is likely on the order of the number of streamwise or spanwise unit cells, this assumption is less valid for the spatial decomposition than the temporal decomposition. Currently, this assumption only affects the calculation of the form-induced stresses. Consequently, errors accrued by this simplification will need to be rigorously assessed if these stresses prove to be a dominant influence in the momentum transport equation, or if this assumption is later used to make different observations.

Simplification of the Double-Averaged Equations
where i denotes a lack of summation over the repeated index. Furthermore, to analyze functions of the streamwise velocity, the i = 1 case is examined: A simplified version of the resolved viscous shear stress may be obtained using the incompressiblity assumption (u i,i = 0), which will prove useful in the following analysis.
Additionally, using the periodic and no-slip boundary conditions in concert with the dimensions of the averaging volume, a simple argument can be constructed to show that w = 0.
Thus, for the i = 1 DANS momentum equation, Equation 2.40 may be simplified as: which provides a convenient method for computing the form-induced stresses, given the time-averaged velocity data. To perform this analysis, Equation 2.42 is first multiplied by φ s and then integrated from some datum, z 0 , to some point of interest, z:

Integral form of Double-Averaged Momentum Equation
By applying the fundamental theorem of calculus, several terms may be rewritten as a function evaluation at the limits of integration: where the constant C(z 0 ) is defined as: The resolved shear stress may be simplified significantly using the Equations 2.44, 2.31 and the no-slip boundary conditions to provide: which links the contribution of each term on RHS of above equation to the variation of φ s ū along the vertical direction. It is worthy to note that anytime the averaged porosity is constant (i.e. φ s = S or φ s = B ), the geometry indicator may be removed from the derivative and the change in u may be examined explicitly.

Initial Studies and Results
Early on in this work, two test cases with different turbulence models were examined. Simulations were conducted on a 15x5x4 bed using the Smagorinsky model and on a 10x5x4 bed using the WALE model. Figure 2.5 shows the double-averaged streamwise velocity u , for each set of results. One need not examine a thorough set of turbulent statistics to identify several differences between the two curves. Looking at the surface flow profile, one might expect to see u approach a constant value, indicating that the domain's surface flow height, H, is larger than that of a fully developed turbulent boundary layer, δ, as discussed in [40]. are seen in [14], suggesting that the Smagorinsky model may be more reliably representing the physics at play. However, given that the Smagorinsky model is known to have trouble accurately predicting near-wall dynamics, and that the WALE model does a particularly good job capturing this behavior [28], this conclusion must be met with some skepticism.
Before moving forward with further analysis, the decision was made to rigorously assess the quality of the mesh (note the inconsistent detail along the boundaries and surface in

Mesh Refinement and Time Averaging Study
Mesh Refinement Study on a Sample Geometry In order to determine a suitable mesh for the detailed modeling of the geometries in Figure   2   Each mesh was used to simulate the model described in Section 2.   in Figure 3.2a and 3.2b, respectively. As seen in Figure 3.

The Influence of Domain Size on the Flow Field
When performing LES, periodic boundary conditions enable one to obtain a fully developed turbulent inflow condition by recycling the fluctuating velocity field exiting the domain at the inlet. One must be careful using these conditions, however, as the constraint of perfect correlations at the periodic domain boundaries may produce nonphysical influences on the   Table 3.1.
flow field if the simulation box is too small [41]. In particular, the domain size directly limits the size of the largest turbulent structures which may be resolved [30,29]. As the larger structures are expected to play an important role in interfacial transport, the fidelity of the simulation depends on using a sufficiently large domain.
To clarify the impact of domain size on the LES results, all four geometric models presented in Figure    The 2x2x4 curve uses the same domain and meshing scheme as the Medium scheme seen in the mesh refinement study.
The correlation coefficient, R ij , may be computed as a function of some time lag, ∆t and some space lag ∆x, as described in Equation 3.1.
Because the correlation coefficient is the covariance of two random variables normalized by the root mean square of each variable, the correlation coefficient may hold values between −1 and 1, where the former denotes perfect anticorrelation and the latter denotes perfect correlation [40].   [19], and therefore corroborates the argument that adequate resolution of large turbulent scales, and consequently the use of a sufficiently large domain, is critical in performing high-fidelity LES on this type of system.

Comparison of The Smagorinsky and WALE Turbulence Models
Within this section, data sets resulting from simulations using the WALE and Smagorinsky turbulence models are compared. Both simulations were run on the 10x5x4 geometry discussed in Section 3.3 using the Medium mesh, detailed in Section 3.2. Analysis is split into two parts. First, double-averaged data is presented to look at large-scale persistent hydrodynamic information. These results are followed by instantaneous flow data, providing a more detailed look at the presence of turbulence near and within the bed.

Double-Averaged Flow Measurements
Double-averaged streamwise velocity profiles for simulations using the WALE and Smagorinsky closure models are shown in Figure 3.6. Given the well documented differences in the treatment of sub-grid stresses between the models (e.g. [28,36]), it is no surprise that there exists discrepancy between the velocity profiles, even when the simulations share identical domains and meshes. However, it is clear that the concerns discussed in Section 2.4 were not due to the mesh or geometric model, as velocity minima for the WALE and Smagorinsky models are found at nearly the same locations (z * = −2.19 and z * = −1.21, respectively).
The corresponding double-averaged Reynolds shear stresses (1,3) for the two closures are presented in Figure 3.7. Both curves assume nearly the same profile, growing linearly in magnitude when approaching the permeable wall from the surface flow and then decreasing rapidly with penetration into the bed. Although similar, the WALE model predicts a more  mild decay with depth while the Smagorinsky model predicts nearly laminar flow within the first pore, indicating that vortical motions produced by the WALE model are able to penetrate further beyond the permeable wall and have a more pronounced role in transport.
One notable feature belonging to these profiles is the rapid decrease in magnitude near the interface, which makes clear a challenge in linking surface and subsurface models via a boundary condition when studying this system. Knowing that the double-average effectively smears the presence of highly localized events across an entire averaging volume, one can expect both the peak and rate of decay of the stress to increase when examining only the time-averaged quantity (note the decay of the TKE below the permeable wall in Figure   3.13). Thus, this double-averaged profile may prove useful in the parameterization of such reduced-order models.
Recalling the formulation of the DANS equations presented in Section 2.3, the differential form of the DANS momentum balance for the i = 1 case is restated in Equation 3.2 for  reference.
The terms of this equation are shown in Figure 3.8. Simulations were performed on the 10x5x4 geometry with both the WALE (3.8a) and Smagorinsky (3.8b) turbulence models.
Simulations with the different models show agreement in regions of constant porosity.
Within the homogeneous fluid region, the body force, which drives the flow, is balanced by the derivative of the Reynolds shear stress. Additionally, viscous forces are seen to have little relevance away from the wall, as noted in [19] for walls with appreciable permeability.
Throughout the homogeneous porous region, the form and viscous drag forces act to balance the body force. According to arguments in [14], one should expect to see the driving force of the flow primarily balanced by form drag in the upper pores and viscous drag in the lower pores, but both turbulence models predict that form drag does a better job than viscous  velocity, which has been found to appear in the first layer of pores [14,15], is related to the penetration depth of turbulence and not solely described by the form and viscous drag terms.
Once again recalling the formulations laid out in Section 2.3, the integrated DANS momentum equation is given in Equation 3.3 for reference, and its constituent terms are plotted in Figure 3.9 for the same simulations just discussed.
As this formulation allows for the isolation of the derivative of u within the homogeneous porous region, where φ s = B , these plots are of interest primarily in the region z * ∈ [−3.5, −0.5]. In this region, results disagree with the conceptual model put forth in [14] and show that form drag extracts more momentum than viscous drag, even deep within the bed. However, due to a small but nonzero value in the momentum residuals shown in Figure   3 In an effort to better understand the contribution of modeled viscosity to the total viscous momentum transport, and thus measure the sensitivity of the flow to the turbulence model, the viscous drag force for each closure scheme has been separated into its resolved and modeled components, shown in suggests that the mesh is fine enough for high quality LES. To more carefully investigate this difference, another mesh refinement study may be useful, executed with a focus on the full double-averaged momentum balance, or at least a more thorough study of the resolved and modeled viscous forces.

Verification of the Drag Force Calculations
In an effort to verify the double-averaging methodology's ability to accurately recover the form and viscous drag terms given the proposed averaging domain and geometric model, both viscous and pressure based contributions to the mean drag force per particle have been computed in two ways. Each drag force was computed indirectly via the DANS equations, as described in Section 2.3, and directly via a surface integration within Nalu according to where Γ is the surface of each sphere and the viscous term, τ ij , uses the effective viscosity (i.e. µ e = µ + µ t ). This force may be decomposed into contributions from pressure and viscous terms, and then time-averaged to obtain the respective mean drag forces.

Instantaneous Flow Measurements
In addition to the use of the double-averaging formulation, analysis of the instantaneous flow data can help to elucidate the role of turbulence in interfacial transport.   Given the identical domain dimensions used between these two simulations, however, the weaker correlation may be explained by looking at Figure 3.8. Within the transition region, the Smagorinsky model predicts an increased and decreased loss of momentum due to viscous stresses and form drag, respectively, when compared to the WALE model. The Smagorinsky model is known to overcompensate for shear in the near-wall limit by producing an unphysical amount of turbulent viscosity [28], and consequently overdamping the fluctuating velocity   suggested that the bed geometry acts to transform Q4 events into Q3 events [33]. Given the increased strength of sweep events within the first pore layer predicted by the WALE model, a more thorough study of quadrant events within the pore spaces may help to elucidate the role of geometry in this proposed transition and clarify the link between the penetration depth of turbulence and the profiles seen in Figure 3.6.

Concluding Remarks
Within these last two chapters, the formulation of a detailed LES for studying turbulent flow over a permeable bed has been examined. The double-averaging methodology has been implemented to aid in the study of large-scale, persistent flow structure in the presence of temporal and spatial heterogeneity. Additionally, the influence of mesh size, domain size and turbulence model on the predicted hydrodynamics has been discussed.
Notably, results of the mesh refinement study show that achieving a fine resolution in

Chapter 4
Implementation of Particle Tracking within Nalu

Introduction
As discussed in Section 1.1, a primary interest of those studying hyporheic exchange, and more broadly, flows of permeable bed forms, is to understand how scalar quantities are exchanged between the surface and subsurface flows. This process is often modeled by lower-order stochastic particle tracking models or advection-dispersion equations [2,9,3], which could likely be made more accurate with the adoption of parameterizations derived from pore scale information [10]. Given the wealth of hydrodynamic detail available at this length scale when using Large Eddy Simulation (LES), provided the necessary information can be extracted, the detailed flow model proposed and evaluated in the previous chapters is a viable candidate for producing such parameters.
The primary challenge in implementing a Lagrangian particle tracking model within an Eulerian system is the continual need to locate the particles' host elements, which are neces-sary for accurate interpolation of the fluid fields to the particles' locations. Several methods have been proposed for handling these problems of search and interpolation. Notably, algorithms have been developed to efficiently handle particle tracking on unstructured and mixed element meshes [42,43,44]. Efforts have also been focused on reducing interpolation errors for particles crossing element boundaries [45].
However, these methods generally don't address the influence of domain decomposition and distributed memory on the efficiency/applicability of the proposed algorithms. As the implementation of this particle tracking module will be within an existing open source software, designed for efficient computation of fluid fields across unstructured meshes on massively parallel machines, host-element determination for each particle becomes a critical matter of computational efficiency. Additionally, the treatment of particles as classes within an object-oriented programming framework adds complexity to the inter-process communication of particle objects within a distributed memory environment.
The work in this chapter presents solutions to the posed challenges, providing the foundation for a particle tracking module. This module can then be incorporated within a detailed LES to extract Lagrangian dispersion statistics from the flow field. First, an introduction to the open source Nalu and STK code bases are presented. Following this, an overview of the particle tracking module and its associated search methods are given. Finally, particle evolution, parallel communication and boundary conditions are explained in detail. It is the author's intent that this chapter will serve as documentation for the current particle tracking implementation and present a thorough introduction for future developers of the module.

Brief STK Vocabulary Discussion
Part of what makes Nalu such a powerful tool is that it relies on several highly optimized software packages and third-party libraries, allowing the code base to continually improve in simulation fidelity, memory management, portability, etc. A particularly important package for development of Nalu's particle module is the Sierra ToolKit, or STK [46]. The STK was created by Sandia National Laboratories to aid in the construction of high-performance software for engineering analysis, and offers modules to support the use of unstructured meshes and geometric searching, among other utilities, in a distributed memory environment.
While most of the particle module development may be explained in the absence of a working knowledge of Nalu's usage of STK, defining a few terms up front may provide significant clarity in the following sections: • Mesh -The mesh is defined as a set of entities, parts, fields and field data, which is divided between MetaData and BulkData.
• Part -A part defines a subset of the mesh and its associated entities. Various kinds of parts exist within a mesh; examples include the locally-owned part, the aura part and the sideset part (if the mesh is read from an Exodus file).
• Field -A field represents a set of data associated with mesh entities. Fields may be created across the entire mesh or on subsets (i.e. parts) of the mesh. Examples include a nodal field which stores the volume of a control volume or an area vector field stored on element edges.
• Entity -An entity is a general term for an object in the discretization including the following types: node, edge, face, element and constraint.
• MetaData -MetaData is the component of the mesh which holds definitions of parts, fields and relationships among and between the two. This data is duplicated on every process and must be created before the BulkData is constructed.
• BulkData -BulkData is the component of the mesh which holds definitions of entities and their ownership (including ghosting information), connectivity data and data associated with the fields defined in the MetaData. Unlike the MetaData, BulkData is distributed amongst processes rather than duplicated, and must be constructed after the MetaData has been finalized.
• Parallel Consistency -STK Mesh is used for several engineering disciplines (e.g. thermal/fluid mechanics and structural dynamics). Because the application using the mesh is not known a priori by the developers, the mesh must be consistent, i.e. it must always observe specific rules, regardless of the governing application. Enforcing consistency becomes significantly harder in parallel, and doing this requires a strict set of rules regarding the distribution of mesh data. The primary rule the Nalu particle module is concerned with is that of parallel ownership, which ensures that each mesh entity is only owned by one process, and that every process with a copy that entity (whether from sharing or ghosting) knows the entity's owning process.
• Aura -The STK Mesh module offers an automatic single-element thick ghosting layer around each process, denoted as the aura. Though more intricate ghosting may be created through custom ghosting, the aura is turned on easily through the Nalu input file and is sufficient for the particle module.
Note that these definitions and more formal explanations for the above terms are available in [47], which provides documentation including conceptual overviews and code examples covering the STK modules.

An Informal Tour of Nalu
Before the particle module may be explained in detail, an elementary understanding of the Nalu code base must be developed. In this section, important structures within Nalu are discussed and the flow of data over the course of a simulation is presented at a highlevel. Developers and those interested in a more complete picture are referred to [35], which provides a more formal treatment of the topic.
Nalu is capable of modeling various flow phenomena, either in isolated or multiphysics simulation. Thus, the code base must support several equation systems (e.g. momentum transport, thermal heat conduction, filtered mixture fraction), which may need to be handled on different meshes with disparate material properties. To manage the organizational challenges associated with this task, Nalu implements realms and transfers to separate each physics' data and allow for coupling.
Nalu creates a realm for each type of physics to be modeled. Realms then develop their own understanding of the discretized domain and associated fields, equation systems and other information relevant to their simulation environment. This level of encapsulation not only provides ease of use for the user/developer, but allows for distinct boundary conditions, initial conditions and output frequency of data to be specified.
In a multiphysics simulation, the realms will have some degree of coupling. Nalu takes advantage of STK's transfer module to support sharing of data between realms. When coupling is required, the fields to be communicated and the sending/receiving realm pair are specified by the user. Transfers will then take place before, during and after time integration.
This ensures consistency of state in a realm's field data throughout system evolution. A

Overview of the Particle Classes
Within Nalu, the particles module exists as a collection of classes handling particle physics, evolution and I/O. A note on the role of each class is given in Table 4.1. This section is meant to be a brief overview of the module's structure, addressing where various functions are carried out and how data is handled over the course of a simulation at a high level.
As mentioned in Section 4.2, Nalu separates physical models by realm. The particle module implements a specific particle-based realm, ParticleRealm, derived from the standard Nalu realm, to handle the numerical environment for particles. While the particle realm follows a load procedure similar to the other realms, it is unique in that it doesn't have to manage assembly of an equation system. This allows the particle realm to have its own set of rules for initialization, directed towards constructing the mesh objects and loading particle data into the domain.
Managers are those objects which are only created once within the particle realm and are responsible for overseeing the most complicated tasks during the simulation. Of the managers in the particle module, ParticleManager has primary oversight. Owning the other particle managers (e.g. ParticleCommManager and ParticlePeriodicManager) and the methods responsible for both particle integration and host element determination, the particle manager is the first place to explore when studying the module in-depth.
Inter-process particle communication is governed by ParticleCommManager. At each timestep, this manager handles the packing, communication and unpacking of particles using tools from the STK Utils module. In the presence of periodic boundary conditions, the ParticleCommManager relies on the ParticlePeriodicManager to create a map for boundaryto-boundary communication. This mapping feature is still under development, so the user is recommended to approach the application of periodic boundary conditions with caution.
Finally, VtuManager handles the output of particle data. Currently, this manager is only able to output coordinate and timestep data in a format suitable for visualization, although extension to output additional data (e.g. particle temperature, unique identifier) is under development. Write operations and the structure of the output are further explained in Section 4.5. Governs simulation environment and multiphysics interactions. ParticleManager Governs particle evolution and host-cell determination. ParticleCommManager Handles particle multi-process particle communication. ParticlePeriodicManager Handles application of periodic boundary conditions. ParticleEquation Defines particle equations (TBI). ParticlePhysics Defines coupling and physics governing particle (TBI). ParticleVariables Handles variables defining particle state. ParticleType Defines material properties and physics (TBI). ParticleIntegrator Governs particle integration (TBI).

TracerParticlePhysics
Defines physics for tracer particles (TBI). ParticleInsertion Stores and loads particle data from input file. VtuManager Outputs particle data.

Host Element Determination
At all times during a simulation, a particle must be aware of its host element to ensure accurate tracing of its immediate environment. Thus, host-cell determination, as discussed in Section 4.1, is a critical component of initialization and time integration. The following subsections briefly discuss the methods used for particle search within Nalu, as well as discussing another application which may find such methods useful.

Face-crossing Search
While several methods exist to query an entity's location within a domain, the challenge of search is greatly simplified when the initial host element is known. Assuming a small time step such that a particle does not traverse several elements in a single update but has left its previous host-element, the particle may be tracked via a face-crossing search. This type of search method is the primary tool used within the particles module for host-element determination, due to its compact search domain.
In the face-crossing search, a particle's coordinates are known at both the previous and updated time levels. With this information, the face of the particle's initial host through which the particle exited may be determined. If this face only corresponds to one element, its clear that the particle has left the domain, and will be sent to the coarse search tool for deletion. If the exit face is shared by two elements, however, then the particle is moved to the element opposite the original host. If this element is the proper host, the particle will continue to integrate. In the instance where the particle has traveled through multiple elements, then the particle will continue through this process, moving from element to element without integrating, until it finds its proper host.

Nalu Coarse Search
During particle initialization, the host element of the particle is not known. In this case, or in the instance the face-crossing search fails to find a suitable host, a more powerful search tool is necessary. STK provides a suite of search tools for Nalu which is already optimized for use on decomposed, unstructured meshes. From this suite, the particle module implements an extension of the Boost R-Tree spatial index to execute a range search, herein deemed a domain search.
An R-tree is a tree-based data structure designed to organize a set of n-dimensional geometric entities into a set of n-dimensional minimum bounding rectangles (MBRs) [48].   host elements. A unique identifier for the object (each particle or element) is first paired with the object's owning process rank, creating an identifier, deemed an ident. Each ident is then paired with an STK geometric object, Boxes for elements and Points for particles in this case, which are the spatial constructs on which the domain search can operate. These pairs are then loaded into vectors, distinguished by the type of geometric object, and can be fed into the coarse search tool. Upon completion, a list of ident pairs is returned, providing a candidate host-element (or multiple candidates) for each particle, along with each object's owning process. With a more rigorous check, the proper host-element may then be found from the candidates put forth by the search.

Application to Fluid-Structure Interaction
Integration of Boost's RTree search method into STK and subsequently Nalu's utilities provides a powerful tool for several applications, including host-cell determination for particles and entity mapping for periodic boundary conditions. One application of the RTree search worth mentioning briefly here, which has a scope that goes beyond this thesis, is in the development of fluid-structure algorithms within Nalu.
During a fluid-structure coupling algorithm, force will need to be transferred from an

The Particle Evolution Algorithm
The primary utility gained from implementing a particle tracking framework in Nalu is to capture a Lagrangian perspective of the flow field. Thus, an update algorithm is needed to advance the particles in time, given the conditions of their local environment. While the algorithm detailed in this section may be adapted to handle multiple-way coupling and particles experiencing varied physics, the current implementation is designed for particles which are strictly one-way coupled to the fluid field, obtaining only the nodal velocity vectors from its environment.

Initialization and I/O
In order to integrate particles through time, the particles must first be created within the domain. Particle initialization currently occurs through the standard Nalu input file. Data corresponding to each particle is read in and stored off into memory during the ParticleManager's load step, then used to create the particles during initialization. Each process begins with the same list of particles to create. After a global domain search and some accounting, each process creates only those particles which reside within its piece of the distributed mesh.
If a particle is to be created on the boundary of multiple processes, it will be created only on the highest rank process.
Output of particle data is handled by the VtuManager class, which enables data to be visualized according to the Visualization Toolkit (VTK) standard [49]. VTK is a widely supported, xml-based format for data visualization. The pvtu format, the only format currently supported by the particles module, is meant for unstructured data created in a distributed memory computation environment, and thus provides an excellent solution for handling output from the particle module.
The process for particle output is straight forward. At each timestep, including initialization, a directory is created in a user-specified path. Within each directory, every process writes a .vtu file containing information regarding its locally-owned particles. Rank zero then creates a .pvtu file which links all of the .vtu files together. Finally, the path of this .pvtu file and the corresponding timestep are written to a .pvd (ParaView data) file, which may be read by ParaView alone to load in all particle data and corresponding timesteps for the simulation. An example directory structure for a two-timestep two-process simulation are shown in Figure 4.3 for reference.

Local Time Integration
Particle evolution can be thought of as a two-tier process. In the bottom tier, particle integration is a procedure local to each process, ending when all particles on the process are either finished integrating and have found their respective host cells or waiting to find a new host cell via global communication. The top tier looks at particle integration as a global procedure, governing local advancement, particle communication and searching, and completion of the integration step. In this section, focus will be given towards explaining the lower tier, i.e. the process of advancing particles locally. Following this discussion, explanations will detail how local evolution fits in to the bigger picture when running a multi-process simulation.
Time integration is carried out using a second order explicit Runge-Kutta scheme, defined calling method. In the more interesting case where active, found particles exist, particles are iterated over and integrated. For a particle going through this process, its containing element's information (e.g. nodal information related to interpolation) is gathered. The particle is then checked to see if it is contained by the element it believes it is owned by.
If not, a new host element must be determined. This is done through a method called handle particle exits element (HPEE), and is described below.
If a particle is indeed contained by the element it believes to be contained by, a check is done to see if the particle is finished integrating. Assuming this is not the case (the next particle is drawn on the contrary), the particle must be either active and predicted or active and yet to be predicted. If the former, the particle is updated according to  left the domain are then deleted, and any particles found locally by the HPEE method are moved into their new hosts. This procedure continues until the process has no remaining active, found particles, at which point it will return to the calling method and get involved with the top tier, global aspects of particle time integration.

Handling A Particle That Exits an Element
Whenever a particle exits its host element, it must enter the handle particle exits element method, which will attempt to find the particle's proper host. HPEE does this by executing a face-crossing search with the element's current and past coordinates, as described in Section 4.4. Based on the results of this search, the particle may take one of several paths.
If the face-crossing search cannot find any exit side or new containing element, the particle is added to a set of particles that will undergo either a local or global domain search in a final attempt to find the proper host. This setup allows the domain search to act only as a support to the face-crossing search, which minimizes the computational resources spent locating particles. Note that if a particle has no proper host (i.e. it has left the domain), it will be marked for deletion during the global domain search and deleted immediately after.
When the face-crossing search is successful in finding a new element, the lost particle is assigned the target element. Once all active, found particles have been iterated over, the particle is then moved into its new host and integration may resume. Note that it is possible that the new element is not the proper host, but an intermediate element between the previous containing element and the new host. In this case, the particle will continue in the integration loop, continually getting sent into HPEE until the proper host is determined.
According to this logic, every time a particle moves beyond its owning process' boundaries, it would be subjected to a global domain search in order to find its new host. Although this is possible with a small number of particles, this frequency of global communication is not feasible with particle numbers of engineering interest. To mitigate this expense, we can take advantage of the aura feature provided by STK. Since the aura provides a single-element layer of ghosting around process boundaries, a particle will always enter a ghosted element before leaving the owning process' domain, assuming a sufficiently small time step and that the boundary is shared. As mentioned in Section 4.1, STK ensures unique ownership of each entity (e.g. an element) within the discretization. Thus, every time the face-crossing search returns an element, the element's ownership is assessed and determined to be either locally owned or ghosted, providing a third path for particles passed to HPEE.
In the instance a particle enters a ghosted element, the element's owning process is queried from the STK entity database and both this process and the ghosted host's identifier are stored on the particle. By storing this data, the particle may be communicated to the correct process and assume integration in the locally owned copy of the ghost element entered previously. This particle is then added to a set of particles to be communicated once the set of active, found particles is empty, and integration at the current timestep continues. Host identification using the face-crossing search and aura layer is shown in Figure 4.6, and the communication procedure for a traveling particle is described in greater detail below.  Then, particles marked for communication are communicated via the ParticleCommManager. The manager does this by packaging up the critical particle information, including target host element and coordinates, and communicating the data with send/receive methods provided by the STK Utils module. Communicated particle contents are then unpacked on the receiving process, where integration of the traveling particle is resumed if necessary.

Global Time Integration
The original copy of the particle is marked for deletion, and any particles marked as such are deleted following the communication step.
If at any time a particle is found locally via the face-crossing search or is communicated to another process, the advancing status of the process owning the new host is set to 0 (i.e. not ready to advance). Parallel reduction is used following the communication step to determine if every process in the parallel group is finished integrating. If so, specified data is written to a file and control is returned to the fluid realm for the next timestep. If the group is not ready to advance, all processes continue through the integration loop, regardless of the number of particles owned, until all particles reach the proper time level.

Boundary Conditions
With the implementation detailed for particles existing in the interior of the domain, a quick discussion on boundary conditions is necessary. Currently, particles interact with domain boundaries in two ways. As a general rule, particles experience boundaries as one-way outflow conditions. Due to the requirement that a particle must always know its host element at the end of every timestep, a search is always performed at the conclusion of each time integration step. When this happens, any particle that has moved beyond the limits of the domain is left hostless, and thus gets deleted. This approach is suitable for simple models where flow is either one-dimensional in nature or in simulations where particles approach impenetrable walls. However, in more intricate geometries, integration error may allow for particles to move beyond walls in an unphysical manner, requiring a more involved treatment of the boundaries.
In addition to the general outflow condition, support for periodic boundary conditions within the particle module is currently under development. Periodic boundary conditions are necessary to track particle trajectories over multiple flow-through times, and thus critical in computing bed residence times or gathering accurate dispersion statistics. While a proof of concept implementation has been demonstrated (detailed below), further work must be done before it is production-ready.
Nalu enforces periodic boundary conditions through a master-slave mapping of boundary nodes, which is applied within the equation systems. Largely due to the particles being objects distinct from the STK Mesh, the particles cannot implicitly respect the periodic boundary condition applied to the fluid realm. Thus, these conditions must be enforced explicitly within the particle realm. Currently, this is achieved by constructing a face-element map, depicted in

Verification of the Particle Evolution Algorithm
In order to verify the evolution algorithm detailed in Section 4.5, particles have been simulated in flow through a laminar cylindrical channel such that their computed paths may be compared against exact theoretical trajectories. Although this test lacks an examination of numerical convergence and removes many of the physical and numerical considerations associated with problems of industrial or academic interest (e.g. boundary layers, turbulence), it provides an adequate arena to assess the essential components of the implementation discussed previously. Eventually, more rigorous testing (e.g. tracing of a rotating flow) will be useful in verifying extensions to the modules and its ability to perform in more dynamic environments.
The channel mesh is unstructured and consists of ≈ 30, 000 hexahedral elements, seen in    Results of this test are shown in Figure 4.9b for simulation taking one thousand time steps with a domain decomposed among ten processors. The relative error between the theoretical and computed positions is on the order of machine error, as expected from the constant velocity field. This simple test demonstrates the successful function of the particle module, showing that a particle is properly communicated between processors and interpolating its  velocity appropriately from an unstructured, distributed mesh.
In addition to the test above, images are shown in Figure 4.10 of particles tracking laminar flow past a sphere. Particles can be seen tracking the flow around the sphere, as well as getting stuck in the low-pressure zone at the upstream stagnation point. These results are shown to indicate qualitatively that the first order behavior of the fluid is captured by the particles, and that this module is on its way to being used in more applicable simulation environments.

Introduction
In this final chapter, the primary conclusions and their implications will be reaffirmed, recommendations for future work on both the detailed LES and the particle tracking module are discussed, and plans for their unity and a suggested path forward towards studying scalar transport are detailed.

Regarding the Detailed Large Eddy Simulation
Throughout Chapters 2 and 3 of this thesis, aspects of a detailed LES model for turbulent flow over a permeable bed have been discussed, examined and tested. Analysis of a mesh refinement study clearly demonstrates the feasibility of using a graded mesh to reduce computational cost, provided the mesh in the transition region is sufficiently resolved to avoid numerical instability. By comparing averaged velocity profiles and autocorrelations predicted by simulations with varying domain/bed size, periodic boundary conditions have been shown to artificially limit the development of large turbulent structures. The doubleaveraging methodology has been applied to study interfacial momentum transport in the presence of significant spatial heterogeneity, notably showing that form drag always dominates viscous drag within the bed and that the Reynolds shear stress rapidly decays within and below the transition layer. Lastly, two different turbulence closures have been studied, highlighting the role of turbulent viscosity in creating smaller, and possibly weaker, turbulent structures in the Smagorinsky model compared to the WALE model.
While these findings are largely focused on details regarding the development of a high fidelity numerical model, several aspects of the physical problem are still ripe for exploration.
First and foremost, a tool with which one can more closely investigate the presence of turbulence is needed. Whether this tool be turbulent kinetic energy budget calculations, sophisticated conditional sampling of velocity fluctuations or measures of vorticity and swirling strength, elucidating the size, strength and motion of turbulent structures within the transition layer will prove useful in quantifying interfacial momentum transport.
One way to more thoroughly probe the pore space is with quadrant analysis at locations that vary in both depth and streamwise position. Several studies have put effort towards studying the structure of turbulence within pore throats across different layers of the bed (e.g. [32,33]). By leveraging the hydrodynamic detail present in LES, one might be able to more clearly explain the origin, evolution and influence of bed turbulence at various depths.
Finally, parametric studies of several model parameters may prove useful in understanding interfacial transport. Such parameters include the body force, the sphere diameter, the interparticle gap spacing and the surface flow height. By examining the changes in the flow field as a function of these parameters, one may gain insight into which Reynolds number (e.g. Re K , Re D ) best characterizes the flow, or how the depth of the surface flow height and the particle diameter affects inner and outer scaling of the turbulent flow.
Although expensive, LES has been shown to be a suitable alternative to DNS, enabling one to resolve the dominant physics at play while reducing computational requirements.
While many questions remain open in the study of this type of system, the work in this thesis further demonstrates that LES can be a useful tool in helping to find answers.

Regarding Nalu's Particle Tracking Module
Within Chapter 4, the implementation of a particle tracking module within the Nalu code base has been proposed and detailed. Critical aspects of I/O, time integration and distributed memory computation have been documented. The module has been shown to successfully integrate a particle across a distributed, unstructured mesh with machine precision and capture the first order behavior of flow past a spherical obstruction with 5,000 particles.
Additionally, a proof-of-concept implementation of periodic boundary conditions has been demonstrated, laying the groundwork for simulations of academic and industrial interest. Several components of the particle tracking module must be improved before particles are ready to be used within an LES. Primarily, a robust algorithm for establishing periodic boundary conditions must be implemented (currently underway), which will enable particle tracking over time scales much larger than a single flow-through time. Once this is accomplished, a method to enable a particle to keep track of its previous positions will be needed.
Giving the particles this awareness will allow one to statistically evaluate particle trajectories in post-processing, providing for the extraction of dispersion information for upscaling to lower-order models. An example of this upscaling procedure is seen in [10], where detailed information is gathered from two-dimensional simulations. To the best of the author's knowledge, this procedure has not been done with detailed, three-dimensional turbulence modeling at the pore scale.
Eventually, the particle module should be expanded to allow for particles which may experience different kinds of physics (e.g. electrodynamics, thermodynamics). Further degrees of coupling may be explored, enabling particle-particle contact or computing heat transfer between the particles and fluid. Advanced extensions may be developed, such as linking against molecular dynamics libraries to study molecular gas dynamics within Nalu. Many possibilities for expansion exist, and it seems that virtually all of them lead to exciting methods of inquiry.

Regarding the Combination of the LES and Particles
One of the primary goals of developing the detailed LES and the particle tracking module is to incorporate them together. By doing this, an extensive data set, rich in both Eulerian and Lagrangian information, may be gathered. This data set can then be used both to interrogate the hydrodynamics throughout the domain and more easily upscale dispersion information to reduced-order models.
While this plan cannot be executed until the particle module has been given the necessary feature set, one alternative to this plan has yet to be explored. Knowing that the particles are purely passive, and their trajectories are entirely governed by the surrounding fluid motion, it might be possible, without the particles, to extract the same dispersion information that is intended to be captured by the particles. The path to do this is not immediately clear.
A reasonable place to start may be looking at the Reynolds shear stresses. As lower-order models often look only to treat the wall-normal and streamwise flow directions, perhaps the turbulent momentum transport described by τ R 13 can provide some sort of parameterization for the probability and magnitude of vertical displacement in a stochastic particle tracking model.
With either this method or the use of particle tracking, this suggested path forward offers both an opportunity to gain direct insight into the mechanisms driving interfacial transport in flows over permeable beds, and inform reduced-order, less expensive models which can be used to predict transport across a host of environmentally applicable length and time scales.