# Geometry and organization of coherent structures in stably stratified
atmospheric boundary layers

Abhishek Harikrishnan1, Cedrick Ansorge2, Rupert Klein1, Nikki Vercauteren3

(1Institute of Mathematics, Free University of Berlin
2Institute of Geophysics and Meteorology, University of Cologne
3Department of Geosciences, University of Oslo)

###### Abstract

Global intermittency is observed in the stably stratified Atmospheric Boundary
Layer (ABL) and corresponds to having large nonturbulent flow regions to
develop in an otherwise turbulent flow. In this paper, the differences between
continuous and intermittent turbulence are quantified with the help of
coherent structures. Eight classes of coherent structures are identified from
literature, most of which are indicated by scalar criteria derived from
velocity fields. A method is developed to geometrically classify structures
into three categories: blob-like, tube-like or sheet-like. An alternate
definition of the intermittency factor $\gamma$ based on coherent structures
is introduced to separate turbulent and nonturbulent parts of a flow. Applying
this conditioning technique and the geometrical characterization on direct
numerical simulations (DNS) of an Ekman flow, we find the following: (i)
structures with similar geometries (either tube-like or sheet-like) are found
regardless of the strength of stratification; (ii) global intermittency
affects all regions of the ABL - viscous sublayer, buffer layer, inner, and
outer layer; (iii) for the highly stratified case, sweep/ejection pairs form
well-separated clusters within the viscous sublayer which can possibly explain
the abundance of hairpin-like vortices with a particular orientation; (iv)
nonturbulent regions are occupied with streamwise velocity fluctuations and
there is a switch between high- and low-speed streaks at a particular height
for all stratified cases.

##  1 Introduction

The Atmospheric (or Planetary) Boundary Layer (ABL) is the lowest part of the
troposphere (up to 2 km), i.e., the region directly responding to changes in
surface conditions. It shares many similarities with the turbulent boundary
layer in a simulated channel flow but differs in two key aspects: First, the
structure of this layer is influenced by the presence of clouds and strong
surface heating/cooling. Second, the outer layer (also known as the Ekman
layer) is affected by the Coriolis effect due to the rotation of the Earth.
Energy exchange at the surface leads to density variations penetrating
throughout the ABL. Strong heating, encountered for example in daytime
conditions, leads to unstable stratification whereas strong cooling,
encountered for example in night-time conditions, yields a stably stratified
ABL (Garratt, 1994). The former is characterized by continuous turbulence
whereas, in the latter, turbulence may be intermittent or even absent (Van de
Wiel et al., 2012). The organization of turbulence in intermittent patches is
denoted by Mahrt (1989) as global intermittency. Our goal is to study the
differences between continuous and intermittent turbulence by comparing the
geometry of the coherent structures.

The demarcation of the flow into turbulent and nonturbulent flow regions was
first investigated by Townsend (1948) in the wake of a circular cylinder. Away
from the center of the wake, he noted that regions of fully developed
turbulent flow and almost laminar flow were separated by a sharp boundary.
This phenomenon was quantified with an intermittency factor, denoted by
$\gamma$, and it was measured with the kurtosis of velocity derivatives.
Similar probe-based measurements were employed by Townsend (1949); Corrsin and
Kistler (1955); Fiedler and Head (1966); Kovasznay et al. (1970); Rao et al.
(1971); Falco (1977). Since these techniques have limited spatial resolution,
the three-dimensional geometry of the turbulent/nonturbulent patches and the
corresponding dynamics are not adequately captured. Da Silva et al. (2014)
used a technique based on vorticity magnitude to segregate the turbulent
(rotational) and nonturbulent (irrotational) regions on the direct numerical
simulations (DNS) of boundary layers, jets and shear-free turbulence. However,
this approach inherits the well-known drawbacks of vorticity magnitude: its
reliance on a subjective threshold and producing false-positives in shear
dominated regions (Lugt, 1979). Ansorge and Mellado (2016) addressed the
latter problem by applying two horizontal high-pass filters (for streamwise
and spanwise directions) on velocity fields close to the wall. They found that
the high-pass filtered velocity fields were able to correctly detect the
absence of turbulence for most of the inner layer. The latter shortcoming of
vorticity magnitude can also be overcome by considering another vortex
criterion like $Q$ (Hunt et al., 1988), which identifies a purely spinning
region by subtracting the strain-rate tensor from the vorticity tensor.

The adoption of such vortex-based partition techniques forms a useful
connection to the notion of coherent structures. Numerous types of coherent
structures such as streaks, sweeps, ejections, vortical structures have been
identified and studied in wall-bounded flows by Robinson (1991), hereafter
referred to as Robinson structures. Previous investigations into coherent
structures such as sweeps and ejections attribute these entities to momentum
transfer in the boundary layer (Katul et al., 1997; Li and Bou-Zeid, 2011). In
their study, Li and Bou-Zeid (2011) report a dissimilarity in the transport
efficiency of scalars and momentum between the neutral and unstable ABL. They
speculate that this is most likely caused by a change in topology of coherent
structures and stress the necessity of studies considering their three-
dimensional extent. Plenty of field measurements and laboratory studies exist
which confirm the presence of hairpin vortices in the neutral ABL (Hommema and
Adrian, 2003; Carper and Porté-Agel, 2004; Huang et al., 2009; Inagaki and
Kanda, 2010). However, their role in stratified ABLs is not well known. All of
this taken together presents a need to study the role of Robinson structures
in the neutrally and stably stratified ABL.

The geometry of coherent structures has also attracted considerable interest.
Moisy and Jiménez (2004), hereafter MJ2004, applied the classical technique of
box counting to 3D scalar fields of vorticity and strain-rate magnitude. They
defined individual structures as a connected set of points satisfying a
thresholding condition and studied their geometry by identifying three
characteristic lengths for each structure. This enabled them to relate the
structures as being close to a sphere, tube, ribbon or a flat sheet. They
concluded that the strain-rate structures change from being flat sheets to
ribbons with increase in threshold. Vorticity strutcures, on the other hand,
shift from ribbons to long tubes with increase in threshold. Bermejo-Moreno
and Pullin (2008), hereafter B2008, developed a non-local methodology to study
the geometry of structures by applying a curvelet transform on a passive
scalar from the DNS of forced, isotropic turbulence. They utilized
differential geometry properties such as Shape Index and Curvedness to
characterize the structures and concluded that blob-like and tube-like
structures prevail in the inertial range and sheet-like structures in the
dissipation range. In an extension of the study, Bermejo-Moreno et al. (2009)
applied the same technique to study the geometry of structures at different
resolutions of the numerical simulation and concluded that the coarsest
resolution was unable to detect sheet-like structures at small scales in 3D
scalar fields of enstrophy and kinetic energy dissipation-rate. Similar
characizerization was also carried out for lagrangian structures (Yang et al.,
2010; Yang and Pullin, 2011) and with the adoption of Minkowski functionals
(Leung et al., 2012). In our current work, we aim to study the geometry of
Robinson structures in neutral and stably stratified ABLs by developing a
methodology derived from the work of MJ2004 and B2008. This will enable us to
establish a detailed comparison on the geometry, distribution and organization
of coherent structures. To this end, we pose the following questions,

(a) How can we extract and study coherent structures in the ABL?

(b) How different are the properties of coherent structures at increasing
levels of stratification?

(c) How does global intermittency affect the self-organization of structures
within the flow in the boundary layer?

The paper is organized as follows: We start by discussing the details of our
DNS database in section 2. Section 3 details our methodology to extract and
geometrically characterize coherent structures. We start with a brief review
of all coherent structures and the indicators used to identify them. Following
this, we present a modified neighbor scanning (NS) algorithm to extract
structures from 3D scalar fields. We discuss the technique of percolation
analysis to objectively choose an appropriate threshold for the entire
boundary layer. Finally, we discuss the methodology derived from MJ2004 and
B2008 for geometrical characterization of the structures. In section 4, we
detail the results for the different layers of the ABL, namely the viscous,
buffer, inner, and outer layers. A significant number of figures produced to
support our remarks are presented in a supplementary file. These figures are
denoted with C#, where # denotes a section number. A special study is also
carried out on the geometry of hairpin-like structures. In section 5, we
investigate the dataset from the viewpoint of conditional one-point statistics
and compare the results to section 4. All observations are summarized in
section 6.

##  2 Overview of the dataset

Case | Line specification | $Ri_{B}$ | $Fr$ | $Re$ | $L_{x}$ | $L_{y}$ | $L_{x}^{+}$ | $L_{y}^{+}$
---|---|---|---|---|---|---|---|---
N | ———– | 0 | $\infty$ | 1000 | 20.4$\delta$ | 20.4$\delta$ | 27 040 | 27 040
S_1 | ………… | 2.64 | 0.02
S_2 | \- - - - - - - | 0.76 | 0.07
S_3 | – — – — | 0.58 | 0.09
Table 1: Parameters of the numerical simulations used in the work. The cases
with the prefix $S$ are stratified and $N$ is the neutrally stratified case.
Reynolds number can be defined with $\delta$ and the laminar boundary layer
depth. The value of the latter case is shown here. If we define $Re$ with the
former, then the value is $26\,450$. Further simulation details can be found
in Ansorge and Mellado (2016).

Figure 1: The schematic shows the three types of subgrids - A, B and C used in
the paper. $x,y,z$ correspond to the streamwise, wall-normal and spanwise
directions respectively. $(\cdot)^{+}$ indicates viscous units as described in
equation 4.

We study the ABL with a simplified physical configuration, namely Ekman flow:
the flow of a stratified incompressible viscous fluid over a smooth flat plate
driven by a uniform pressure gradient which experiences steady rotation around
the wall-normal axis. It can be ideally compared to the ABL in the limit of
neutral stratification (Coleman et al., 1990). We inherit the simulation
datasets from the work of Ansorge and Mellado (2016). Their set-up uses the
incompressible (divergence-free) Navier-Stokes equations with the addition of
a Coriolis component. It was solved under the Boussinesq approximation, where
density variations were neglected except in determining the buoyancy. They
also impose the f-plane approximation where the Coriolis parameter acts in a
horizontal plane along the vertical axis and is assumed to be constant. A no-
slip boundary condition is applied at the wall and the domain is doubly
periodic in the horizontal direction. Different levels of density
stratification are considered and stratification is imposed via a Dirichlet
(fixed-value) boundary condition at the top and bottom boundaries. Further
details regarding the simulation can be found in Ansorge and Mellado (2014,
2016) and Ansorge (2016).

Figure 2: Horizontal vorticity magnitude slices normalized with its RMS over
that slice are shown here at $y^{+}\approx 100$. (a, b, c, d) correspond to
cases N, S_1, S_2, S_3 respectively.

Our database is comprised of three stratified cases (S_1, S_2, S_3) along with
a neutrally stratified case (N). All simulation parameters are summarized in
table 1. The computational domain has $3072\times 512\times 6144$ points
corresponding to the streamwise ($x$), wall-normal ($y$) and spanwise ($z$)
directions respectively. For convenience, we define three horizontal grid
sizes to communicate our results. Subgrid A represents the entire
computational domain whereas subgrids B and grid C show $1/3$ and $1/6$ of the
domain respectively as illustrated in figure 1.

A classification based on the time evolution of vertically integrated
turbulence kinetic energy (TKE) suggests that the three cases S_1, S_2 and S_3
are all very stable (Ansorge and Mellado, 2014). In such a regime, turbulence
dies out in some regions showing large regions of little to no activity. This
is evident from figure 2(b, c, d) where slices of vorticity magnitude
normalized with its root mean square (RMS) over the horizontal slice is shown.

For all cases, the Reynolds number ($Re$) is $26\,450$ in terms of the
geostrophic wind velocity $G$, the boundary layer height under neutral
conditions $\delta$, and viscosity $\nu$:

$Re\equiv G\delta/\nu$ (1)

where $\delta\equiv u_{\tau}/f$, $u_{\tau}$ is the friction velocity and $f$
is the Coriolis parameter. We define the strength of stratification with two
dimensionless parameters namely the global bulk Richardson number ($Ri_{B}$)
and Froude number ($Fr$). The former is defined as,

$Ri_{B}\equiv B_{0}\delta/G^{2}$ (2)

where $B_{0}$ quantifies the buoyancy difference between the top and bottom
layer. The Froude number is defined as

$Fr=\frac{G^{2}}{B_{0}\Lambda}$ (3)

where $\Lambda$ is the Rossby deformation radius. $Ri_{B}$ is related to $Fr$
as $Ri_{B}=Fr^{-1}\frac{\delta}{\Lambda}$. Throughout the paper, we consider
the vertical direction in terms of viscous units such that

$y^{+}\equiv yu_{\tau}/\nu$ (4)

where $u_{\tau}$ is the friction velocity.

##  3 Characterization of individual structures

###  3.1 Identification

Type | Coherent structure | Indicator | Region | Color specification
---|---|---|---|---
Quantitative | High-speed streaks | $u^{\prime}>0$ | $y^{+}<40$ |
| Low-speed streaks | $u^{\prime}<0$ | $y^{+}<40$ |
| Sweeps | $u^{\prime}>0,v^{\prime}<0$ | $y^{+}<1550$ |
| Ejections | $u^{\prime}<0,v^{\prime}>0$ | $y^{+}<1550$ |
| Vortices | $\frac{1}{2}(||\Omega||^{2}-||S||^{2})>0$ | $y^{+}<1550$ |
| Shear layers | $|\omega|$ | $y^{+}<80$ |
| Backs | $|\omega|$ | $y^{+}<1550$ |
Qualitative | Pockets | streamlines | $y^{+}<5$ | -
| Bulges | $|\omega|$ | $y^{+}<1550$ | -
Table 2: Types of coherent structures considered in the study. The regions
indicate where the dataset has been cut in the $y$-direction. In the case of
sweeps, ejections, vortices, backs and bulges where the structures are found
throughout the boundary layer, the region is limited to a $y^{+}<1550$ as it
already includes a good portion of the outer layer. The maximum $y^{+}$ value
is $6849.54$.

Even though a precise mathematical formulation for a “coherent structure” has
not yet surfaced, we adopt the idea from Pope (2001) and understand that it is
a connected region in space which persists in time. We refer to these
connected regions as a “structure” or “entity” throughout the paper, and we
will use the taxonomy of boundary layer structures introduced by Robinson
(1991), to study them in an organized manner. We divide the Robinson
structures into two major categories: (a) quantitative, where structures can
be extracted and geometrically characterized and (b) qualitative, where they
can only be observed. In this section, we present a brief overview on all
Robinson structures, their purpose, which region of the boundary layer they
are prominently observed in and how they can be extracted/visualized.

The seminal work of Kline et al. (1967) has experimentally revealed streaks
close to the wall of turbulent boundary layers. These regions of relatively
slow-moving fluid can be clearly seen until $y^{+}\approx 40$ in flat plate
boundary layers (Kline, 1978). Streaks are extremely helpful in understanding
the nature of flow within the viscous sublayer ($y^{+}<5$). Within the context
of ABL, however, these structures are usually studied with Large Eddy
Simulations (LES) above the viscous sublayer as the first grid point in LES is
relatively far away from the wall (Khanna and Brasseur, 1998; Drobinski and
Foster, 2003; Jayaraman and Brasseur, 2021). Computing the streamwise velocity
fluctuation $u^{\prime}$, both low-speed ($u^{\prime}<0$) and high-speed
($u^{\prime}>0$) streaks can be readily observed and used to quantify the
effect of global intermittency within this region.

Kim et al. (1971) noted that streaks not only move downstream but also away
from the wall. This movement, which is initially slow, extends rapidly outward
after a critical height. This process is described as ‘streak lifting’ or
ejection. Conversely, Corino and Brodkey (1969) identified a stream of fluid
moving towards the wall called sweep. The quadrant technique (Wallace et al.,
1972) is well suited to detect these entities. It splits the product of
streamwise and wall-normal velocity fluctuations into four categories: Q1
$(+u^{\prime},+v^{\prime})$, Q2 $(-u^{\prime},+v^{\prime})$, Q3
$(-u^{\prime},-v^{\prime})$ and Q4 $(+u^{\prime},-v^{\prime})$. Out of the
four, Q2 and Q4 are used to identify ejection and sweep events respectively.
They can be used to quantify momentum transport throughout the ABL (Li and
Bou-Zeid, 2011; Katul et al., 1997, 2006; Narasimha et al., 2007).

Vortices were described by Küchemann (1965) as the ‘sinews and muscles’ of
fluid flows, and they have been studied extensively over the years. Despite
their immense importance, there exists no universally accepted definition and
a vast range of criteria are currently used. For a review on various criteria
see Cucitore et al. (1999); Günther and Theisel (2018). They fall under two
broad categories, one of which does point-wise characterization on the
velocity gradient tensor ($\nabla v$) for every instant in time (Eulerian),
and the other follows fluid particle trajectories (Lagrangian). While both
have their merits, we choose the former technique as the latter would be
computationally expensive for our dataset. In choosing an appropriate
criterion, we follow the work of Chakraborty et al. (2005), who found that
popular vortex criteria such as the $Q$ (Hunt et al., 1988), $\lambda_{2}$
(Jeong and Hussain, 1995) and $\Delta$ criterion (Chong et al., 1990) identify
very similar regions as vortices. They also define an equivalent threshold
which allows for the visualization of similar structures among the criteria at
any non-zero threshold. Therefore, a choice among these three will not affect
our qualitative results. We choose the $Q$-criterion which is defined as,

$Q=\frac{1}{2}(||\Omega||^{2}-||S||^{2})>0$ (5)

where $\Omega=\frac{1}{2}[\nabla v-(\nabla v)^{T}]$ is the vorticity tensor,
$S=\frac{1}{2}[\nabla v+(\nabla v)^{T}]$ is the strain rate tensor. The
$Q$-criterion identifies vortices as regions where vorticity dominates over
the strain. They can be visualized only above the viscous sublayer and are
thought to play a central role in turbulence production (Robinson, 1991).

Closely related to the lift-up of low-speed streaks, the wall-normal shear
layers are known to exist in the buffer layer ($5<y^{+}<30$) and beyond until
$y^{+}\approx 80$. They are characterized by high values of the instantaneous
velocity gradient ($\partial u^{\prime}/\partial y$) and are capable of
retaining coherence for up to $1000^{+}$ in the streamwise direction
(Johansson et al., 1987). These structures may roll up into vortices
(Robinson, 1991) and may be responsible for near-wall turbulence production
(Alfredsson et al., 1988). Since vorticity magnitude ($\omega$) is known to
highlight shear regions, it can be used to study the near-wall shear layers
and the shear layers on the order of $\delta$-scale termed as backs.

Pockets have been documented by Falco (1977, 1980); Kim et al. (1987), but
their contribution to turbulence production is not fully understood yet. These
structures are understood as ‘footprints’ of outer region motions in the inner
layer and are thought to contribute to the generation of hairpin vortices (Chu
and Falco, 1988). They are known to exist within the viscous sublayer and are
observable with diverging streamlines (Robinson, 1991).

Another coherent structure frequently referred to in literature are the bulges
which are very large scale motions (VLSMs) or superstructures existing between
the free stream and the outer edge of the boundary layer. They show numerous
narrow incursions, sometimes extending into the buffer layer (Corrsin and
Kistler, 1955; Kovasznay et al., 1970). They also attenuate small-scale
fluctuations close to the wall (Marusic et al., 2010) are also documented in
ABLs (Shah and Bou-Zeid, 2014; Katul, 2019) and can instantly be observed with
vorticity magnitude (Robinson, 1991).

A summary of the Robinson structures, the chosen indicators and the region of
identification is presented in table 2. Pockets and bulges are qualitatively
understood whereas the others can be extracted and geometrically characterized
with the technique presented in the following sections.

###  3.2 Extraction

Following MJ2004, we define a structure as a connected set of points for which
the indicator is greater than a given threshold. Therefore, a point $x$
belongs to a coherent structure if,

$\alpha(x)>\tau_{p}$ (6)

where $\alpha$ is the indicator for extracting coherent regions and $\tau_{p}$
is an appropriate threshold. The process of selecting $\tau_{p}$ objectively
is described in section 3.3. The algorithm used to extract the structures
employs the neighbor scanning (NS) approach of MJ2004 where all points in the
domain are scanned successively, once $\tau_{p}$ is identified. Every point is
considered to be the center of a $3\times 3$ cube and the neighbors are any of
the 26 points (faces, edges and vertices of the cube) that surround it.
Although this approach yields good results at high threshold values where the
structures are spaced apart, we notice that some structures tend to be grouped
together at smaller threshold values. An example is shown in figure 3 where
the NS algorithm incorporates all neighbors close to the center point even
though it can be clearly seen upon visualization that some neighbors belong to
two different structures. This issue arises due to the fact that the NS
algorithm doesn’t know how the surface mesh is being generated by the
visualization algorithm. Therefore, we introduce an extension to the NS
approach by correcting the identified neighbors with the marching cubes
algorithm (Lorensen and Cline, 1987) so that it is accurate for visualization
purposes. The details of the NS+MC approach is presented in Appendix A.

Figure 3: Drawback of the NS approach is illustrated here. An extracted
structure is shown in (a). If we look at a plane as shown in (b), we can
instantly observe where the issue occurs. Blue and red scatter points belong
to two different structures. When the NS algorithm considers the center point
indicated by the tail end of arrows as the center of a $3\times 3$ cube, it
assumes that points belonging to both structures are neighbors.

###  3.3 On the choice of a threshold

As noted by Green et al. (2007), the crux when dealing with any Eulerian
criterion for structure extraction is its reliance on a user-defined
threshold. If the chosen value is too low, a complex, interconnected, sponge-
like structure can be seen. On the other hand, a large value can result in too
few structures being identified. MJ2004 consider this problem to be analogous
to that of a percolation transition, where the complex, interconnected
structure only appears above a critical threshold. This threshold, henceforth
known as the percolation threshold $(\tau_{p})$, presents a natural, non-
subjective way of choosing a global threshold value. Del Álamo et al. (2006)
have pointed out that $\tau_{p}$ remains unchanged with increasing Reynolds
number in channel-flow turbulence.

Figure 4 (a-g) shows the computation of the ratio $V_{max}/V$ for increasing
values of $\tau$ for all indicators in table 2, where $V_{max}$ is the volume
of largest structure at that threshold and $V$ is the volume of all
structures. When this ratio is $1$, i.e., $V_{max}=V$, we see a large,
interconnected structure. As $\tau$ is varied, the magnitude of the slope of
the ratio $V_{max}/V$ increases and at $\tau_{p}$, it reaches a maximum. At
this point, we see more disconnected structures. Due to the high cost of this
computation, especially at low values of $\tau$, we had to ensure that this
transition was seen with as few values of $\tau$ as possible. The entire range
of $\tau$ is split with a chosen number (usually $1000$) of evenly spaced
values. According to our experience, a smooth plot can be generated with the
first 50 values of $\tau$. If the transition doesn’t lie within this region,
the order of magnitude of the chosen number can be changed.

For wall-bounded flows, inhomogeneity in the wall-normal direction must be
taken into account when choosing a global threshold (Del Álamo et al., 2006).
For instance, a threshold chosen to specifically highlight structures close to
the wall may show an incomprehensible amount of structures away from it and
vice versa. This was clearly elucidated by Nagaosa and Handler (2003) who
found that by nondimensionalizing the $Q$-criterion with its RMS for every
wall-normal plane, the structures appeared to be more uniformly spread
throughout the channel. Therefore, we can rewrite equation 6 as,

$\alpha(x)>\tau_{p}\;\overline{\alpha^{{}^{\prime}2}(y)}^{1/2}$ (7)

Figure 4: Percolation analysis on seven quantitative indicators as shown in
table 2. From left to right each row are the (a) high-speed streaks, (b) low-
speed streaks, (c) sweeps, (d) ejections, (e) Q-criterion, (f) vorticity
magnitude restricted to $y^{+}\approx 80$ and (g) vorticity magnitude until
$y^{+}\approx 1550$. $V_{max}$ refers to the volume of the largest structure
in the domain and $V$ is the volume of all structures at a particular
threshold $\tau$. $|\cdot|$ is the modulus of the threshold value and is used
for low-speed streaks, sweeps and ejections (b-d) where all thresholds are
negative. A legend corresponding to our different ABL cases in shown in (a).
The thick red line highlights the region of maximum slope.

where $\overline{\alpha^{{}^{\prime}2}(y)}^{1/2}$ is the RMS of the indicator
over wall-normal planes. It is important to note that this normalization is
necessary only when dealing with a large domain. If we focus on a specific
region as in the case of streaks where we confined ourselves to the region
$y^{+}<40$, then the normalization is no longer useful as first plane and last
wall-normal planes are not too far from each other and a single threshold
value will highlight structures uniformly.

This strategy works well at neutral and intermediate levels of stratification
i.e., cases N, S_2 and S_3. However, for the highly stratified case, we notice
that structures tend to aggregate into clusters at non-zero thresholds (see
figure 20 (a)). If $\tau_{p}$ is still used to extract structures, then an
entire cluster would be identified as a single structure. To circumvent this,
we developed a technique where percolation analysis can be applied iteratively
to break down a complex cluster into simpler components. This enables every
structure to be identified at a unique value of $\tau$. In other words, every
structure has its own threshold value. The technique is hence a multilevel
percolation threshold approach and is presented in Appendix B.

###  3.4 Geometrical characterization

Once $\tau_{p}$ is identified, the structures are extracted at this threshold
with the NS+MC algorithm described in section 3.2 and we now proceed to
geometrically characterize them with an approach similar to B2008. The main
steps along with our modifications are presented below:

Figure 5: Panel (a) shows the visualization space with a three-dimensional
projection and a set of two-dimensional orthogonal projections composed of
$\lambda$ and $\hat{C}$ and $\hat{S}$ and $\hat{C}$. Figure adapted from
Bermejo-Moreno and Pullin (2008). Panel (b) shows $\hat{S}$ for some commonly
observed structures (A - D) and Robinson structures (E - K) extracted from
case N. A and B are blob-like structures, C, D, E, G, H and I are tube-like
structures, F, J and K are sheet-like structures. Color coding for the
Robinson structures is according to table 2.

  * $(a)$

Fractal dimension: The first step is to filter structures which are sparser
than lines i.e., structures having a mean fractal dimension $\langle
D_{\alpha}\rangle<1$. This is important for two reasons. First, some
structures extracted with the NS+MC approach are noise-like consisting of one
or a few connected points and cannot be meaningfully interpreted. Second, the
computation of curvedness and stretching (described in point b) for
geometrical characterization is dependent on the area and volume of a
structure and therefore requires that we consider only closed surfaces.
Structures having $1\leq\langle D_{\alpha}\rangle\leq 2$ are excluded because
they tend to have complicated geometries and closed, stretched, tube-like
surfaces which include a large number of grid points in one dimension and very
few in the other two get classified under this category (see C1 for an example
where a significant amount of structures are shown to fall under this
category). The fractal dimension of a structure is computed with the box-
counting approach described in MJ2004 - the extracted structure is placed in
the domain which embeds it. Then, we use a fixed grid scan approach where a
box of size $r$ is moved across the domain without overlapping and the number
of boxes $N_{\alpha}(r;\tau_{p})$ containing a value is counted. This is done
for numerous values of $r$ which are chosen as $2^{-n}L$ where $n$ is varied
from 0 to 9 and $L$ corresponds to the length of the largest dimension. It
should be noted that values of $r$ smaller than the size of a single grid cell
are ignored. The fractal dimension is defined as follows,

$D_{\alpha}(r)=-\frac{\text{d ln }N_{\alpha}(r)}{\text{d ln }r}$ (8)

  * $(b)$

Shape index, Curvedness and Stretching: The iso-surfaces of the structures at
$\tau_{p}$ which have $\big{\langle}D_{\alpha}(r)\big{\rangle}>1$ are
extracted and we proceed to geometrically characterize them with three
parameters. The first two, shape index $(S)$ and curvedness $(C)$ are
differential-geometry properties which represent the local shape and the
intensity of curvature (Koenderink and Van Doorn, 1992). In terms of the
principal curvatures $\kappa_{1}$ and $\kappa_{2}$ (see Do Carmo (2016) for a
description of principal curvatures), they can be written as,

$S=\bigg{|}-\frac{2}{\pi}arctan\bigg{(}\frac{\kappa_{1}+\kappa_{2}}{\kappa_{1}-\kappa_{2}}\bigg{)}\bigg{|},\quad
C=\mu\sqrt{\frac{\kappa_{1}^{2}+\kappa_{2}^{2}}{2}}$ (9)

where $|.|$ indicates the absolute value. $\mu=3V/A$ with $V$ and $A$ the
volume and area of the structure respectively, is used to nondimensionalize
$C$. The area is computed by summing up the area of each triangle in the
surface mesh. For accurate volume calculation, the signed volume of
tetrahedron as shown in Zhang and Chen (2001) is used. By introducing a
feature center for the joint probability density function (jpdf) of $(S,C)$,
B2008 have represented these local properties in a non-local sense. For
asymmetrical jpdfs, the feature centers, henceforth denoted as $\hat{S}$ and
$\hat{C}$, will lie closer to the region of higher density. The formula to
compute the feature center is shown in Appendix C of B2008. A third parameter,
denoted by $\lambda$ characterizes the amount of stretching experienced by the
structure and is denoted by

$\lambda=\sqrt[3]{36\pi}\frac{V^{2/3}}{A}$ (10)

Algorithm 1 Geometrical characterization of structures extracted from scalar
fields
---
Input: Scalar field of the indicator, Percolation threshold $\tau_{p}$
| (i) | At $\tau_{p}$, extract structures with algorithm A.
---|---
(ii) | Embed every structure in a box with the maximum extent along every direction.
(iii) | Compute fractal dimension with equation 8.
(iv) | If $\langle D_{\alpha}\rangle>1$, compute Shape Index, Curvedness and Stretching.
Output: $\hat{S}$, $\hat{C}$ and $\lambda$.
  * (c)

Visualization space: Once the three parameters are computed for every
structure, they are plotted onto a visualization space which provides a
graphical means of educing clusters of similar structures. All structures can
be broadly classified into three categories as blob-like, sheet-like or tube-
like depending on where they lie in the visualization space. For instance, a
sphere which has $(\hat{S},\hat{C},\lambda)=(1,1,1)$ is a blob-like structure
and any other structures which are close to this region can be thought of as
approaching the geometry of a sphere. Similarly, tube-like structures can be
identified close to a cylinder which has
$(\hat{S},\hat{C},\lambda)=(0.5,1,\lambda)$, where $\lambda$ can indicate the
stretching experienced by the tube. The parameter $\hat{C}$ solely determines
if a structure is sheet-like ($\hat{C}=0$) or transitioning towards it
($\hat{C}\rightarrow 0$). Figure 5 shows examples of commonly encountered
structures and an example of each quantitative coherent structure extracted
from case N. Intuitively, one can understand from figure 5(a) that A, B are
blob-like, C, D, E, G, H and I are tube-like and the rest are sheet-like
structures. This is confirmed with a K-means clustering algorithm initialized
with 3 cluster centers. The results are shown in C2.

To summarize, we identify the Robinson structures with indicators as described
in section 3.1 and split them into two categories: quantitative and
qualitative. For structures belonging to the former, we compute the global
percolation threshold $\tau_{p}$ by applying percolation analysis as described
in section 3.3 and then extract the structures with the NS+MC approach
described in section 3.2. Finally, we filter all structures with
$\big{\langle}D_{\alpha}(r)\big{\rangle}<1$ and compute the stretching and
feature centers of shape index and curvedness (section 3.4). This summary is
also presented as a four-step algorithm (Algorithm 1) and results of this
analysis on the ABL dataset are shown in the subsequent section.

##  4 Geometry of structures in the ABL

The results are organized into five sub-sections, four of which correspond to
different regions in the boundary layer namely, the viscous sub-layer
$(y^{+}<5)$, buffer layer $(5<y^{+}<30)$, inner layer $(y^{+}<1000)$, and
outer layer $(y^{+}>50)$. This allows us to study the geometrical
characteristics and interactions among coherent structures within each region.
The final sub-section compares the organization of hairpin structures between
the neutral and stably stratified ABL. Although results for all cases are
discussed here, the figures for cases S_2 and S_3 are shown in the
supplementary.

###  4.1 Viscous sub-layer $(y^{+}<5)$

Figure 6: Isosurfaces of $u^{\prime}$ in the viscous sublayer for grid B. The
structure highlighted in light green is the longest low-speed streak within
the domain and the one in pale yellow is the longest high-speed streak. (a)
corresponds to case N and (b) to case S_1. The color specification is
according to table 2 where the regions shaded in red correspond to
$u^{\prime}>0$ and the ones in blue are $u^{\prime}<0$.

The coherent structures of interest in this region are the low-and high-speed
streaks, sweeps, ejections and pockets. Other Robinson structures are not
observable as the flow in this region is locally laminar.

First, we discuss the low-and high-speed streaks. Both are identified with the
indicator $u^{\prime}=u-\langle u\rangle$. Here, $\langle\cdot\rangle$ denotes
averaging over wall-normal planes for $5$ time steps. Figures 6 and C3.1 show
the isosurfaces of $u^{\prime}$ for all cases at the global percolation
threshold $\tau_{p}$. Even qualitatively, one can observe that the
distribution and geometry of structures in all cases is different from each
other, particularly between cases N and S_1.

Case N shown in figure 6(a) has thin, elongated streaks and their distribution
is most comparable to the ones described in flat plate boundary layers (FPBL)
(Kline et al., 1967; Kline, 1978; Robinson, 1991). Both low-and high speed
streaks tend to be coherent streamwise for several hundred viscous lengths
with the longest spanning $6883^{+}$ and $2650^{+}$ units respectively. These
numbers are much higher than those reported in R1991. The strength of
stratification appears to play an important role in the streamwise and
spanwise coherence of low-speed streaks. We note that their streamwise and
spanwise coherence tends to increase with an increase in the strength of
stratification (from S_3 to S_1). In case S_1, the longest low-speed streak
tends to occupy a significant portion of the domain as highlighted in figure
6(b). If we consider only the high-speed streaks for this case, we would
essentially see large empty patches reminiscent of the nonturbulent patches
seen in figure 2(b) at $y^{+}\approx 100$. This suggests that global
intermittency causes a deceleration of the flow in the viscous sublayer. These
large low-speed regions force the high speed streaks to form tight clusters.
The nature of the high speed streaks, compared to the ones observed in FPBLs,
doesn’t appear to change much with their mean streamwise and spanwise
coherence ranging between $255^{+}$ to $313^{+}$ and $44^{+}$ to $50^{+}$
units respectively. The volume of occupancy of these entities falls with
increase in stratification with the converse being true for low speed streaks.
For S_3, low-and high speed streaks occupy $3.1\%$ and $16.3\%$ of grid B
respectively, S_2 with $4.8\%$ and $14.3\%$ respectively and S_1 with
$18.77\%$ and $8.76\%$ respectively. It can be qualitatively observed through
figures 6 and C3.1.

Figure 7: Contour plots of (a) ejections on top of low-speed streaks and (b)
sweeps and low-speed streaks at $y^{+}\approx 3.58$ for case N and grid B. The
color specification is according to table 2 where low-speed streaks are
colored blue, sweeps colored green and ejections colored yellow.

With figures 7, 8 and C3.2, we study sweeps, ejections and their interactions
with low speed streaks at $y^{+}\approx 3.58$. Sweeps and ejections were
computed with the quadrant technique and correspond to
Q2($-u^{\prime},+v^{\prime}$) and Q4($+u^{\prime},-v^{\prime}$) events.
$v^{\prime}$ is computed the same way as that of $u^{\prime}$. Akin to
streaks, sweeps and ejections in case N (see figure 7) behave similarly to
FPBLs. Confirming previous studies (Corino and Brodkey, 1969; Bogard and
Tiederman, 1986), we observe that ejection events show less streamwise
coherence than low speed streaks and several ejection events may arise from a
single low speed streak (see figure 7(a)). This is possible since both are
regions of negative $u^{\prime}$. This character does not appear to change for
cases S_2 and S_3. For case S_1 as seen in figure 8(a), some ejections arise
out of the edges of the long low speed streak seen in figure 6(b). However,
most of the structures in the inner region show no interactions. This
reinforces the idea that nonturbulent patches may simply be low-speed streaks.

Sweep events in all cases do not interact with the low speed streaks and
appear to be aligned to ejections and low-speed streak. This suggests that
sweeps and ejections often exist as a pair and each pair is generally
associated with the existence of a vortex above the viscous sublayer
(Robinson, 1991). With our extraction algorithm, we counted the number of
sweep/ejection structures for each case. For case N, the ratio of
sweep/ejection structures is $3033/6676\,(0.45)$, case S_3 with
$3104/4842\,(0.64)$, case S_2 with $2798/4980\,(0.56)$ and case S_1 with
$2194/3054\,(0.72)$. These numbers illustrate the following: (i) not all
sweeps have an ejection associated with it, (ii) although case S_2 and S_3 do
not show much of a difference, an increase strength of stratification clearly
reduces the number of sweeps and ejections. These results together with
figures 7 and 8 suggest the following: (i) case N will have relatively the
least amount of vortices whereas case S_1 will have the highest close to the
wall and (ii) case S_1 will show that vortices are organized in tight
clusters. These views will be revisited in the next sub-section and sub-
section 4.5.

Pockets are visualized with streamlines by the line integral convolution
technique (Cabral and Leedom, 1993) in figures 9 and C3.3. Since pockets are
thought to be the ‘footprints’ of a structure aloft the viscous sublayer that
induces wallward fluid motion, the outer structure is often associated with
vortices. Wallward motion is also known to cause streamlines to diverge, hence
we have marked 3 regions of diverging streamlines in each case. Apart from
case S_1, we note that most regions of diverging streamlines can be associated
with a sweep/ejection pair thereby suggesting the presence of a vortex. Along
with sweeps and ejections, pockets are reviewed again in sub-section 4.5.

Figure 8: Similar to figure 7, contour plots of (a) ejections on top of low-
speed streaks and (b) sweeps and low-speed streaks at $y^{+}\approx 3.58$ are
shown for case S_1 and grid B.

###  4.2 Buffer layer $(5<y^{+}<30)$

Figure 9: Pockets are shown with diverging streamlines for (a) case N and (b)
case S_1. Three regions are highlighted in each case which show examples of
pocket-like regions. These are overlayed with sweeps and ejections. The
combination of pocket-like region and a sweep/ejection can possibly highlight
the presence of a vortex above the viscous sublayer. Streamlines are
visualized with Line integral convolution for the grid C. The color
specification is according to table 2 where sweeps are colored green and
ejections colored yellow.

The reason why we didn’t characterize the geometry of structures in the
viscous sublayer is because very few of them, which have a fractal dimension
greater than 1, are fully formed within this region. Most of these structures
extend into the buffer layer and sometimes beyond. In this section, we first
establish a comparison between the structures formed within the buffer layer
itself i.e., structures which start and end within the buffer layer and
structures which extend from the viscous sublayer to this region. Next, we
compare the geometry of these structures subject to different levels of
stratification.

Figure 10: The visualization space for all Robinsion structures except
pockets, backs and bulges are shown here with joint pdfs. (a - f) correspond
to high-speed streaks, low-speed streaks, sweeps, ejections, vortices and
shear layers. Filled contours are used for case S_1 (with dark shade of red
showing the region of high density) whereas unfilled contours with dashed
lines are for case N. The number of structures between contours are also
indicated - dark blue for case S_1 and black for case N. Additionally, numbers
in parenthesis indicate structures which start within the viscous sublayer and
end in the buffer layer. Numbers outside parenthesis indicate structures
within the buffer layer itself.

We summarize briefly the procedure to obtain the figures 10 and C3.6. All
indicators described in table 2 are subjected to percolation analysis and an
objective threshold $\tau_{p}$ is obtained for each of them. This step is done
constrained by grid A for $y^{+}<1550$. At $\tau_{p}$, all structures are
extracted with the algorithm described in sub-section 3.2. Extraction is done
on grid B to keep the cost of computation low. Next, two filters are applied
to ensure that the surfaces are compatible with our geometrical
characterization. The first applies the box counting technique (cf. 3.4(a)) to
filter all structures with $\langle D_{\alpha}(r)\rangle<1$. The second one
removes all structures close to the side walls as they are incomplete, open
surfaces. Finally, the parameters characterizing the geometry of a structure
are computed. These are shape index, curvedness and stretching (see sub-
section 3.4(b,c)). Since we have characterized a large number of structures,
we represent the visualization space with joint pdfs rather than scatter plots
as shown in figure 5(a).

We make a couple of general observations from figures 10 and C3.6: (i) within
each category, no obvious distinction can be made in the distribution of
geometry of structures across different levels of stratification, i.e.,
strength of stratification does not appear to impact the geometry, (ii) few
structures in the buffer layer are sheet-like ($\hat{C}\rightarrow 0$) and the
remaining are in transition ($0<\hat{C}<1$).

Results from figures 10(a,b) and C3.6(a,b) suggest that high-and low-speed
streaks are mostly tube-like structures for all cases. Few high-speed streaks
exist in the buffer layer itself whereas low-speed streaks show more activity.
This view is supported by Khanna and Brasseur (1998) who found that the
intensity of high-speed regions decreases much more rapidly with height than
that of the low-speed regions for the near-neutral ABL. Our results suggest
that this scenario is true for all cases, regardless of the strength of
stratification. Most structures also seem to extend outward from the viscous
sublayer (indicated by the number of structures within parenthesis). This is
contrary to the view suggested by R1991, who pointed out that high-speed
streaks are localized structures extending a limited distance from the wall.

Willmarth and Lu (1972) state that streaks and ejections are the dominant
contributor to Reynolds shear stress $\rho\overline{u^{\prime}v}$ and show a
high activity in the buffer layer of FPBL. They note that sweeps are dominant
when $y^{+}<15$ and ejections dominate when $y^{+}>15$. This view is supported
by Kim et al. (1987) where DNS of channel flow at multiple Reynolds numbers
are investigated. We can interpret a similar picture from figures 10(c,d) and
C3.6(c,d). The number of sweep events within the buffer layer itself
(indicated by the number of structures outside parenthesis) is much higher
than the ejection events for all cases. This suggests a higher sweep activity
within the buffer layer. Since sweeps are responsible for wallward flow, these
results suggest a high wallward flow in the sweep regions within the buffer
layer.

Vortices and shear layers begin their journey within the buffer layer. We
surmised in the previous sub-section that case N will exhibit the least vortex
activity and case S_1 the highest based on the sweep/ejection pair. Results
from figure 10(e) indicate the complete opposite with substantial vortex
activity for case N. However, less vortex activity in case S_1 may be due to
the clustering effect we described earlier where a very large region is
extracted as a single structure (see figure 20(a)). Such structures which span
the domain are excluded from the geometrical characterization due to
computational time. Even if included, they will be characterized as a sheet-
like structure since the area is quite large i.e., the ratio $\mu$ becomes
very small and $\hat{C}\rightarrow 0$. Apart from this, one can note that
vortices identified by the $Q-$criterion are mostly tube-like structures in
the buffer layer for all cases. With similar geometry, shear layers are found
in less abundence within this region. This is interesting as both
$Q-$criterion and $|\omega|$ are vortex indicators suggesting different levels
of activity within the buffer layer.

Figure 11: Similar to figure 10, the visualization space for high-speed
streaks, low-speed streaks, sweeps, ejections, vortices and shear layers (a -
f) are shown here. In this case, numbers in parenthesis indicate structures
which start within the viscous sublayer and end in the inner layer. Numbers
outside parenthesis indicate structures beyond the buffer layer.

###  4.3 Inner layer $(y^{+}<1000)$

The formation of high-and low-speed streaks within the inner layer itself,
i.e., beyond $30<y^{+}<1000$ is insignificant and at most times nonexistent
(see figures 11(a,b) and C3.7(a,b)). Although we say inner layer, streaks are
only studied until a $y^{+}\approx 40$ (as shown in table 2). Therefore, we
can only conclude that there is very little activity of streaks in the region
$30<y^{+}<40$. Most of the streaks originate from the viscous sublayer. This
further supports the notion that high-speed streaks extend futher outward from
the wall than suggested by R1991. While these structures are mostly tube-like,
low-speed streaks show a downward shift in both $\hat{S}\hat{C}$ and
$\lambda\hat{C}$ planes showing a clear transition towards sheet-like
structures.

Sweeps and ejections in figures 11(c,d) and C3.7(c,d) exhibit a distinguishing
feature from all other structures we have examined so far. When visualized
with scatter points, one can observe two distinct clusters (tube-like and
sheet-like) for the structures formed beyond the buffer layer. By K-means
clustering algorithm, we segregate the points into two clusters (see C3.4).
All structrues conforming to each cluster are then extracted and added to an
empty 3D scalar field so that they can be reconstructed and visualized (see
C3.5). We see that tube-like and sheet-like structures are distributed quite
randomly throughout the domain, even along the wall-normal direction. Although
higher activities of sweep and ejection events are evident from the figures,
much more sweep events exist than ejection events for all cases except S_3
where the converse is observed. The ratio of sweep/ejection events is also
much more consistent here than in the viscous sublayer. For case N, the ratio
of sweep/ejection events is $892/820(1.09)$, case S_1 with $296/201(1.47)$,
case S_2 with $614/520(1.18)$ and case S_3 with $440/635(0.69)$. If we ignore
case S_3, we can still observe the same trend spotted in the viscous sublayer
that an increase in the strength of stratification reduces sweep and ejection
events. With all these information, we can now comment on the overall behavior
of these events: sweeps dominate the inner and buffer layer suggesting high
wallward flow within the sweep regions whereas ejections dominate the viscous
sublayer suggesting high outwards flow within the ejection regions from the
wall.

A significant number of structures can be seen for figures 11(e,f) and
C3.7(e,f) showing high vortex and shear layer activity in the inner layer. For
shear layers, we have limited the region of study to $y^{+}\approx 80$ (see
table 2). Any structures which extend beyond this wall-normal height are
classified as backs and are discussed in section 4.4. Interestingly, both
vortices and shear layers exhibit tube-like and sheet-like geometry with the
latter having more sheet-like structures. The reconstructed scalar field for
shear layer structures in case S_1 (see C3.11) traces out a vague outline of
the turbulent/non-turbulent interface suggesting that the outer edges of
turbulent regions exhibit high amounts of shear.

Figure 12: Visualization space for backs for case S_1 (filled contours) and
case N (unfilled contours). The number of structures between contours are also
indicated - dark blue for case S_1 and black for case N.

###  4.4 Outer layer $(y^{+}>50)$

Although this region presents a lot of overlap with the inner layer, we focus
only two Robinson structures namely, extended shear layers also known as backs
and bulges which are often referred in literature as very large scale motions
(VLSM).

The inability of vorticity magnitude to segregate vortex and shear structures
is used to find backs, specifically looking at structures which originate
within the buffer layer and extend beyond $y^{+}>80$. Although plenty of
$|\omega|$ structures can be seen in the inner layer (indicated by the number
of structures in parenthesis in figures 11(f) and C3.7(f)), only a handful of
structures fulfill the criteria for backs. The highest amount of backs out of
the total number of structures $(165/105119)$ is detected for case N, followed
by $(35/47391)$, $(11/29342)$, $(8/59906)$ for cases S_2, S_1 and S_3
respectively. Figure 12 and C3.8 show the geometrical characterization of
backs for all cases. The results suggest no explicit dependence on the
strength of stratification and similar geometries exist for all cases. Also,
their extent suggests that these structures are not on the order of $\delta-$
scale as indicated in previous works (Chen and Blackwelder, 1978; Robinson,
1991). The longest streamwise extent for all cases varies from $26^{+}$ to
$959^{+}$ units and the longest wall-normal extent varies from $15^{+}$ to
$110^{+}$ units. Previous research with these structures also show that these
structures are inclined in the streamwise direction. We computed the angle of
inclination with $\theta=arctan(\Delta y^{+}/\Delta x^{+})$ for all cases and
found that some structures are inclined in the streamwise direction with a
mean inclination angle for each case ranging between $50^{\circ}$ to
$67^{\circ}$. Although structures in case S_2, S_3 and N show a wide
distribution of inclination angles ranging between $16^{\circ}$ and
$80^{\circ}$, case S_1 has no structure with an inclination angle smaller than
$45^{\circ}$. This can possibly indicate a weak dependence on the strength of
stratification.

The presence of $\delta$-scale structures is most apparent in figure 13(a),
where contours of vorticity magnitude can be seen interacting with the outer
layer. These large, $\delta$-scale structures are referred in literature as
bulges. We identify these structures by visualizing contours of $|\omega|$.
One should note that the $|\omega|$ has been normalized with its RMS as given
in equation 7. This evenly highlights the entire boundary layer including the
large scale structures in the outer layer. We also note that bulges exist in
all cases regardless of the strength of stratification and an example has been
shown for each case in figures 13 and C3.8. The point of interaction of the
outer layer and the boundary layer shows a spike in $|\omega|$ suggesting the
presence of large scale vortices.

###  4.5 Special study: Hairpin-like structures

In this section, we pursue a specific category of vortical structures which
are shaped like a hairpin (or horseshoe). This is done for several reasons:
(i) upon visualization of case S_1 (see figure 14(a)), we see a remarkable
amount of vortices shaped like hairpins organized in clusters in the region
$30<y^{+}<200$ (Harikrishnan et al., 2020), (ii) these hairpin-like structures
seem to be oriented in a particular direction and (iii) as pointed out
previously, a large region spanning the whole domain in case S_1 is identified
as a single structure and not used for geometrical characterization, thereby
restricting our analysis in this region.

To extract the hairpin-like structures, we apply multi level percolation
analysis (MLP) as discussed in Appendix B on a small domain of size $300\times
90\times 600$ (streamwise $\times$ wall-normal $\times$ spanwise) as shown in
figures 20(a,b). A small domain along with $V_{max}/V=0.5$ as the stopping
criterion is chosen in order to keep the computational costs bearable. As seen
in figure 20(c), this stopping criterion is sufficient to educe individual
structures. All structures within this region were extracted and from these
100 typical structures were chosen manually. We also chose $100$ hairpin-like
structures from other cases as well.

Figure 13: Contours of vorticity magnitude along the XY plane highlighting the
presence of $\delta$-scale bulges. Here, the full domain in the streamwise
direction is visualized until $y^{+}<1550$. (a) shows case S_1 and (b) shows
case N. The wall-normal direction, which is $1.75\delta$, has been exaggerated
3 times to show the structures clearly.

Figures 15 and C3.9 show the geometrical characterization of hairpin-like
structures for all cases. Interestingly, no sheet-like structures can be found
for case S_1 and the curvedness ranges between $0.5$ and $1$. This can be
attributed to the application of MLP which produces ‘simple structures’
existing at their own threshold. For other cases, the structures can be
complex, sometimes with several structures interacting at once which tends to
have an impact on the curvedness value. Example visualizations on the top row
of figure 15 highlight an intriguing aspect of hairpin-like structures - the
head of the hairpins seem to be oriented in a similar direction. This view can
be further confirmed with figure 14(a), where a majority of hairpin-like
structures are turned at some angle to the streamwise direction. This behavior
is not observed in other cases.

Finally, the importance of pockets and sweep/ejection pairs in vortex
generation is reviewed. Figure 14(b) shows hairpin-like structures, pockets,
sweeps and ejections superimposed over one another. Upon close inspection, we
find some structures, close to the wall, which can be directly associated with
diverging streamlines and a sweep/ejection pair. This reaffirms the idea that
pockets are caused by the presence of an outer rotating structure and the
presence of a sweep/ejection pair can be associated with the leg of a hairpin-
like structure. These results illustrate the geometry and distribution of
Robinson structures quite well, despite the analysis being restricted to one-
third of the entire domain.

Figure 14: (a) shows a volume rendering of case S_1 with $Q-$criterion at
$\tau_{p}$. The domain is restricted to a wall-normal height of $y^{+}\approx
370$ to remove large scale structures. The white arrow indicates the
approximate orientation of a majority of hairpin-like structures and the blue
arrow indicates the streamwise velocity direction. (b) shows the hairpins
chosen for analysis after MLP. It is projected over streamlines showing
pockets, sweeps and ejections at a $y^{+}\approx 3.58$ in the same region
where MLP was performed.

Figure 15: Visualization space with for 100 hairpin-like structures are shown
for case S_1 (filled contours) and case N (unfilled contours). Two examples
for each case are visualized on the right hand side. The subscripts indicate
the case to which the structure belongs.

##  5 Discussion

This section outlines some practical implications of our work: First, in sub-
section 5.1, we introduce an alternative quantification of global
intermittency that overcomes the known shortcomings of vorticity magnitude.
Next, in sub-section 5.2, we compare the results from section 4 to the ones
obtained with a conditional analysis of one-point statistics that does not
make use of feature extraction.

###  5.1 Global intermittency

Figure 16: Top panel shows vertical profiles of $\gamma(y^{+})$ for case N:
(a) vorticity magnitude without normalization with RMS over every wall-normal
plane, (b) vorticity magnitude with normalization with RMS over every wall-
normal plane and (c) $Q-$criterion after normalization with RMS for every
wall-normal plane. Bottom panel shows vertical profiles of $\gamma(y^{+})$ for
$Q-$criterion for case (d) S_1, (e) S_2 and (f) S_3. All figures show profiles
for increasing values (left to right in each figure) of threshold. The first
dashed line corresponds to $0$ threshold whereas the solid line indicates
percolation threshold.

The intermittency factor according to Pope (2001) is given by

$\gamma(x,t)=\langle H(|\omega(x,t)|-\omega_{\tau})\rangle$ (11)

where $H$ is a Heaviside function, $|\omega|$ is the vorticity magnitude,
$\omega_{\tau}$ is a small threshold, $x$ is a streamwise position in the
flow, $t$ is time and $\langle\cdot\rangle$ represents averaging over
horizontal planes. The equation 11 is generalized to include other indicators
as follows,

$\gamma(y)\equiv\langle H(\alpha(x,y,z)-\alpha_{\tau_{p}})\rangle$ (12)

where $\alpha$ is an indicator and $\alpha_{\tau_{p}}$ is the global
percolation threshold for the indicator. Time $t$ has been excluded since we
are dealing with a single timestep in our current work. The choice of $\alpha$
and $\alpha_{\tau_{p}}$ are critical as they differentiate the turbulent and
nonturbulent parts. This definition of $\gamma$ is general as $\alpha$ can be
any vortex indicator.

In previous works, Da Silva et al. (2014); Ansorge and Mellado (2016) used
vorticity magnitude $|\omega|$ as the indicator. Figure 16(a, b) show
$|\omega|$ applied to the case N. The normalization with RMS over wall-normal
planes is not applied for figure 16(a), whereas it is done for figure 16(b).
The profiles of $\gamma(y^{+})$ in figure 16(a) is similar to Ansorge and
Mellado (2016); Kovasznay et al. (1970). However, this clearly highlights the
problems associated with $|\omega|$. First, the profiles indicate that for low
values of $\tau$, the region close to the wall, where shear is dominant, is
filled with structures. Second, the choice of $\tau_{p}$ is not applicable for
the entire boundary layer. The latter issue can be fixed with a normalization
of wall-normal planes with its RMS as shown in Figure 16(b). Here, the
profiles are more or less vertical for high values of $\tau$. It should be
noted that very close to the wall ($y^{+}<30$), minor deviations occur. The
former issue can be fixed by using an alternate indicator. Figure 16(c) shows
case N with $Q-$criterion. The profile correctly shows that regions close to
the wall are relatively ‘less filled’ with structures when compared with
$|\omega|$.

Figures 16(d, e, f) show the intermittency factor profiles for the three
stratified cases with $Q-$criterion. Unlike case N, the profiles show
significant variations beyond $y^{+}=30$. We justify the choice of a single
threshold by noting that the standard deviation of the profiles is small. The
implication of a modified method to quantify global intermittency on
conditional statistics is shown in the next subsection.

###  5.2 A conditional statistics perspective

Sub-section 5.1 suggests that the $Q-$criterion can be used to partition the
turbulent flow in analogy to the turbulent/nonturbulent partitioning used for
the conditional study of large-scale intermittency in turbulent flows (cf.
Shah and Bou-Zeid, 2014; Ansorge and Mellado, 2016). We present the
conditional analysis for high-and low speed streaks, sweeps and ejections
here.

The flow is split into turbulent and nonturbulent regions with equation 12
using the global percolation threshold $Q_{\tau_{p}}$ as discriminator. The
wall-normal planes are averaged such that every height is represented by a
single value (cf. figure 17). The vertical profiles show that the effect of
low-speed streaks (cf. figure 17(b)) can extend well over $y^{+}>100$ for all
cases indicating that an analysis up to $y^{+}\approx 40$ (used in the
previous section) may not be sufficient to understand the role of these
structures.

Figure 17: Vertical profiles of streamwise velocity fluctuation for (a) high-
speed streaks, (b) low-speed streaks, (c) sweeps and (d) ejections. Line
specification is according to table 1 for the different cases. Black color
indicates turbulent regions and teal shows the nonturbulent regions.

When conditioning to high-speed ($u^{\prime}>0$; figure 17(a)) and low-speed
streaks ($u^{\prime}<0$; figure 17(b)), we find a distinct difference between
the lowest part of the boundary layer $(y^{+}<~{}30)$ and the region aloft. In
vicinity of the surface, the positive fluctuations in the turbulent part are
larger than in the laminar one; on the contrary, when looking at the low-
speed streaks, the negative perturbations are larger in the laminar parts. In
other words, intense high-speed streaks tend to belong to the turbulent
partition, intense low-speed to the laminar one. This is in accordance with
the dominant origin of turbulence in the buffer layer, namely shear
production: given the boundary condition $u(y^{+}=0)=0$, a larger velocity
immediately implies more shear, thus more intense turbulence production, and
larger likelihood of the corresponding region to be detected as turbulent by
the conditioning method introduced in section 5.1 (and vice versa in low-speed
streaks).

This picture is inverted aloft $y^{+}\approx 30$: There, the nonturbulent
partition contains faster moving fluid than the turbulent one in high-speed
streaks and in low-speed streaks (the negative perturbation is smaller, so the
fluid is actually faster in an absolute sense). In this context, we have to be
aware that the simple characterization of ‘streaks’ by the sign of
$u^{\prime}$ bears no geometric characterization whatsoever, i.e., above
$y^{+}=40$, we characterize regions with positive or negative fluctuation of
streamwise velocity which is not necessarily a ‘wall streak’ in the classical
sense. In this context, the non-locality of turbulence above the buffer layer
becomes important. Turbulent fluid mostly originates from below and is thus
slower, whereas as non-turbulent fluid originates from above and is therefore
faster Ansorge and Mellado (2016). This is manifested in larger positive
velocity fluctuations in the laminar partition of high-speed streaks and less
negative, i.e., also faster in streamwise direction, velocity perturbations in
low-speed streaks.

The figures 17(c,d) partially confirm the results from section 4 where sweep
events appear to dominate over ejection events in the buffer layer with the
highest activity in the turbulent regions for case S_1. Sweeps also show
significant activity beyond $y^{+}\approx 200$ suggesting strong wallward flow
in these regions. Correspondingly, the region beyond the buffer layer and
within $y^{+}\approx 200$ shows significant ejection activity for all cases in
the turbulent regions. This was not seen in section 4. The results indicate a
switching behavior for case S_1 with sweeps and ejections in turbulent
regions. This suggests alternating regions of wallward and outward flow at
different heights for case S_1.

We define the nonturbulent fluctuating regions as $\langle
u^{\prime}_{nonturbulent}\rangle\equiv\langle u^{\prime}\lessgtr
u^{\prime}_{\tau_{p}}\rangle$ for $Q<Q_{\tau_{p}}$. Interestingly, for case
S_1, $\langle u^{\prime}_{nonturbulent}<0\rangle$ shows significant activity
in the buffer layer. We examine this by counting the number of points which
satisfy the aforementioned condition over the counts of all nonturbulent
points at every wall-normal plane. This result is shown in figure 18 and C4.1.
For case S_1, we observe that low-speed streaks dominate over the high-speed
streaks until a $y^{+}\approx 31.26$, after which there is a switch and the
high-speed streaks dominate. This behavior is observable for other stratified
cases as well but effect is not as significant. However, the switch takes
place at a much lower height at $y^{+}\approx 6.18$ and $y^{+}\approx 7.55$
for cases S_2 and S_3 respectively.

Figure 18: For case S_1, vertical profiles indicating the number of points $C$
satisfying the condition $\langle u^{\prime}_{nonturbulent}>0\rangle$ and
$\langle u^{\prime}_{nonturbulent}<0\rangle$ over the total nonturbulent
points $C_{total,nonturbulent}$ are shown on the left plot. Line specification
is according to table 1. On the right, horizontal slices of $\langle
u^{\prime}_{nonturbulent}>0\rangle$ (red), $\langle
u^{\prime}_{nonturbulent}<0\rangle$ (blue) and $Q-$criterion (black) are shown
at (a) $y^{+}\approx 3.58$, (b) $y^{+}\approx 13.45$, (c) $y^{+}\approx 31.26$
(d) $y^{+}\approx 104$.

##  6 Conclusions

We have established a detailed comparison on the geometry of Robinson
structures among three stably stratified cases and a neutrally stratified case
of an Ekman flow. For this purpose, we have developed a methodology derived
from the works of MJ2004 and B2008. This employs a modified neighbor scanning
algorithm to extract structures from 3D scalar fields. The identified
neighbors are corrected with the marching cubes visualization algorithm.
Performance of this extraction can be further improved by considering an
efficient implementation of marching cubes such as the Lewiner algorithm
(Lewiner et al., 2003). Structures having a fractal dimension less than 1 are
not considered. We also exclude structures which are attached to the domain
walls. For all remaining structures, three parameters are computed, namely
shape index, curvedness and stretching. This allows us to geometrically
classify a structure as blob-like, tube-like or sheet-like. Due to the
computationally expensive nature of this methodology, the analysis was limited
to one-third of the entire domain i.e., grid B. To understand the behavior of
structures on the entire domain i.e., grid A, we apply conditional one-point
statistics on four Robinson structures namely high-and low-speed streaks,
sweeps and ejections. The conditioning is based on a new definition of
intermittency factor introduced in equation 12. The following conclusions are
reached:

  1. (a)

For all Robinson structures, similar geometries (mostly tube-like and sheet-
like) are observed across all levels of stratification. Blob-like structures
are nonexistent.

  2. (b)

Even if structures are clustered, as shown for sweeps and ejections in the
inner layer, no relation can be made between a type of geometry such as tube-
like or sheet-like and its distribution across the spatial domain.

  3. (c)

Although the viscous sublayer is locally laminar and characterization with
turbulent and nonturbulent subvolumes is not applicable in this region, the
effect of global intermittency can be clearly seen here (see figure 6(b)). The
defining characteristic is the association of an ejection with a low speed
streak. Large regions where no ejections arise out of a low speed streak (as
seen in figure 8(a)) are reminiscent of the nonturbulent region observed in
the buffer layer and beyond. Therefore, global intermittency can be
generalized to indicate active/inactive regions rather than
turbulent/nonturbulent regions to include its effect in the viscous sublayer.

  4. (d)

Sweep/ejection pairs are constrained to the active regions and form tight
clusters in the viscous sublayer. We hypothesize that this can possibly lead
to the formation of vortical structures shaped like hairpins above the buffer
layer.

In future work, we will focus on the dynamics of the Robinson structures. With
a suitable tracking algorithm, the Robinson structures can be followed in time
and changes in the geometry can be studied. One can take this further by
considering the visualization space as a phase space which hosts all possible
geometrical states of the system. This allows for the estimation of two
instantaneous metrics namely the local dimension and inverse persistence
(Faranda et al., 2017) which may offer insight into the preferred geometrical
shape of a structure.

The methodology described in section 3 requires significant computational
power. In total, $10^{6}$ structures were geometrically characterized with
this technique. To ensure reproducibility of the results, all codes used in
the course of the work are made available here:
https://github.com/Phoenixfire1081/GeometryCharacterization.git.

Acknowledgements

The authors are grateful to Elie Bou-Zeid for fruitful discussions during the
course of the work. Due to the size of the database, computing was performed
on the Jülich Supercomputing Centre under the project hku24 and stadit. The
data was stored under the project dns2share. Additionally, data was also
stored at the Zuse Institut Berlin (ZIB). This work has been supported by the
Deutsche Forschungsgemeinschaft (DFG) through grant CRC 1114 “Scaling Cascades
in Complex Systems”, Project Number 235221301, Project B07: Self-similar
structures in turbulent flows and the construction of LES closures. Cedrick
Ansorge acknowledges funding through the European Commission’s ERC Starting
Grant no. 851347.

##  Appendix A Neighbor Scanning approach with Marching Cubes neighbor
correction

For any indicator $\alpha$, we apply the threshold $\tau_{p}$ and scan all
resulting points. With the neighbor scanning (NS) approach, all points
surrounding the center of a $3\times 3$ cube are added as neighbors. This
process is repeated until no new neighbors can be found. To ensure that the
extracted structures are accurate to visualization, we correct the identified
neighbors with the marching cubes (MC) algorithm. Although numerous
visualization algorithms are available, we choose MC as it is employed by our
visualization software Amira (Stalling et al., 2005). This algorithm uses one
grid cell (eight vertices of a cube) to determine how the triangles should be
constructed. Depending on which of the eight vertices satisfy the threshold,
there are $2^{8}$ ways in which the surface mesh can be represented. For
further details, the reader is referred to the original paper by Lorensen and
Cline (1987). Out of the 256 combinations, there are 92 cases where the
surface mesh is not continuous and in which only 2-4 nonadjacent points of the
grid cell satisfy the threshold. When these points belong to different
structures, the NS approach is unable to distinguish between them as it
doesn’t know how the surface mesh is going to be constructed. We, therefore,
present the following modification to the NS approach:

Figure 19: Examples of surface mesh construction with the marching cubes
algorithm. (a) shows a continuous mesh which will be a part of only one
structure. (b, c, d) show examples where multiple structures can be a part of
a single cube.

  1. $(a)$

Construction of a lookup table: We visualize every single case out of the 256
combinations and group indices of the cube. For cases with a continuous mesh
as in figure 19(a), the lookup table is assigned an empty array. This is also
done for the case where all or none of vertices satisfy the threshold as it
corresponds to all the points being inside the structure or outside
respectively. When more than one mesh element is present as shown in Figure 19
(b, c, d), the indices belonging to each mesh element are grouped together.
The lookup table is shown in table 3.

  2. $(b)$

Finding the case in the lookup table: Once the neighbors are identified with
the NS approach, all $2^{3}$ subcubes need to be matched with a case in the
lookup table. If $V=\\{v_{0},\ldots,v_{7}\\}$ is the set of vertices of the
subcube, each vertex is assigned a value $1$ when the neighbor is added and
$0$ otherwise. The corresponding case can be found by:

$case=\sum_{i=1}^{8}2^{i}[V(i)=1]$ (13)

  3. $(c)$

Correcting the neighbors: Once the appropriate case is found, the grouped
indices can be retrieved from the lookup table. No changes are made to the
neighbors if the array is empty. For a non-empty array, the ‘right group’ of
indices have to be kept as the neighbors and the remaining are discarded. This
is chosen as the one which contains the center point of the $3\times 3$ cube.
Since this is our current point of interest and is common to all the subcubes,
we can be certain that the mesh element corresponding to this is a part of the
structure. For implementation purposes, these steps are shown as an algorithm
(see Algorithm A).

Algorithm A Extracting structures from scalar fields with Neighbor Scanning
and Marching
---
Cubes
Input: Scalar field of the indicator, threshold $\tau$
| (i) | Convert scalar field with floating-point numbers as boolean data type where a point
---|---
|  $x==True$ when $x>\tau$.
(ii) | For every point in the domain where $x==True$,
|  | (a) | Assign $x$ with a structure label.
---|---
(b) | Check neighbors of the point by considering the point to be at the center
| of a $3\times 3$ cube.
(c) | If neighbors have a $True$ value, label them with the same one as assigned before.
(d) | Split the $3\times 3$ cube as $8$ subcubes.
(e) | For each cube, determine which Marching Cubes case it belongs to by using
| equation 13.
(f) | With a table look-up (table 3), the grouping can be seen. If table look-up is
| empty, the labels are not altered. If not, the group with the center point
of the
| cube is kept and the remaining labels are altered.
(g) | For all remaining neighbors, repeat from (b).
Output: Scalar field with labels for structures, $V_{max}$, $V$.
Case | Indices | Case | Indices | Case | Indices | Case | Indices | Case | Indices
---|---|---|---|---|---|---|---|---|---
1 | - | 53 | [2], [4, 5] | 105 | [3], [5, 6] | 157 | - | 209 | -
2 | - | 54 | [2], [0, 4, 5] | 106 | [0, 3], [5, 6] | 158 | - | 210 | -
3 | - | 55 | - | 107 | [3], [1, 5, 6] | 159 | - | 211 | [1], [4, 6, 7]
4 | - | 56 | - | 108 | - | 160 | - | 212 | -
5 | - | 57 | [3], [4, 5] | 109 | - | 161 | [5], [7] | 213 | -
6 | [0], [2] | 58 | - | 110 | - | 162 | [0], [5], [7] | 214 | -
7 | - | 59 | [3], [1, 4, 5] | 111 | - | 163 | [1, 5], [7] | 215 | -
8 | - | 60 | - | 112 | - | 164 | [0, 1, 5], [7] | 216 | [3], [5]
9 | - | 61 | [2, 3], [4, 5] | 113 | - | 165 | [2], [5], [7] | 217 | -
10 | - | 62 | - | 114 | - | 166 | [0], [2], [5], [7] | 218 | -
11 | [1], [3] | 63 | - | 115 | - | 167 | [1, 2, 5], [7] | 219 | [1], [3, 4, 6, 7]
12 | - | 64 | - | 116 | - | 168 | [0, 1, 2, 5], [7] | 220 | -
13 | - | 65 | - | 117 | - | 169 | [3, 7], [5] | 221 | -
14 | - | 66 | [0], [6] | 118 | - | 170 | [0, 3, 7], [5] | 222 | -
15 | - | 67 | [1], [6] | 119 | - | 171 | [1, 5], [3, 7] | 223 | -
16 | - | 68 | [0, 1], [6] | 120 | - | 172 | - | 224 | -
17 | - | 69 | - | 121 | [3], [4, 5, 6] | 173 | [2, 3, 7], [5] | 225 | -
18 | - | 70 | [0], [2, 6] | 122 | - | 174 | [0, 2, 3, 7], [5] | 226 | [0], [5, 6, 7]
19 | [1], [4] | 71 | - | 123 | [3], [1, 4, 5, 6] | 175 | - | 227 | -
20 | - | 72 | - | 124 | - | 176 | - | 228 | -
21 | [2], [4] | 73 | [3], [6] | 125 | - | 177 | - | 229 | -
22 | [0, 4], [2] | 74 | [0, 3], [6] | 126 | [1], [7] | 178 | - | 230 | [0], [2, 5, 6, 7]
23 | [1, 2], [4] | 75 | [1], [3], [6] | 127 | - | 179 | - | 231 | -
24 | - | 76 | [0, 1, 3], [6] | 128 | - | 180 | - | 232 | -
25 | [3], [4] | 77 | - | 129 | - | 181 | [2], [4, 5, 7] | 233 | -
26 | - | 78 | - | 130 | [0], [7] | 182 | [2], [0, 4, 5, 7] | 234 | -
27 | [1], [3], [4] | 79 | - | 131 | [1], [7] | 183 | - | 235 | -
28 | - | 80 | - | 132 | [0, 1], [7] | 184 | - | 236 | [2], [4]
29 | [2, 3], [4] | 81 | [4], [6] | 133 | [2], [7] | 185 | - | 237 | -
30 | - | 82 | [0, 4], [6] | 134 | [0], [2], [7] | 186 | - | 238 | -
31 | [1, 2, 3], [4] | 83 | [1], [4], [6] | 135 | [1, 2], [7] | 187 | - | 239 | -
32 | - | 84 | [0, 1, 4], [6] | 136 | [0, 1, 2], [7] | 188 | - | 240 | -
33 | - | 85 | [4], [2, 6] | 137 | - | 189 | - | 241 | -
34 | [0], [5] | 86 | [0, 4], [2, 6] | 138 | - | 190 | - | 242 | -
35 | - | 87 | [1, 2, 6], [4] | 139 | [1], [3, 7] | 191 | [0], [6] | 243 | -
36 | - | 88 | - | 140 | - | 192 | - | 244 | -
37 | [2], [5] | 89 | [3], [4], [6] | 141 | - | 193 | - | 245 | -
38 | [0], [2], [5] | 90 | [0, 3, 4], [6] | 142 | - | 194 | [0], [6, 7] | 246 | -
39 | - | 91 | [1], [3], [4], [6] | 143 | - | 195 | [1], [6, 7] | 247 | -
40 | - | 92 | [0, 1, 3, 4], [6] | 144 | - | 196 | [0, 1], [6, 7] | 248 | -
41 | [3], [5] | 93 | [2, 3, 6], [4] | 145 | - | 197 | - | 249 | -
42 | [0, 3], [5] | 94 | - | 146 | - | 198 | [0], [2, 6, 7] | 250 | -
43 | [1, 5], [3] | 95 | [1, 2, 3, 6], [4] | 147 | [1], [4, 7] | 199 | - | 251 | -
44 | - | 96 | - | 148 | - | 200 | - | 252 | -
45 | [2, 3], [5] | 97 | - | 149 | [2], [4, 7] | 201 | - | 253 | -
46 | [0, 2, 3], [5] | 98 | [0], [5, 6] | 150 | [2], [0, 4, 7] | 202 | - | 254 | -
47 | - | 99 | - | 151 | [1, 2], [4, 7] | 203 | [1], [3, 6, 7] | 255 | -
48 | - | 100 | - | 152 | - | 204 | - | 256 | -
49 | - | 101 | - | 153 | - | 205 | - |  |
50 | - | 102 | [0], [2, 5, 6] | 154 | - | 206 | - |  |
51 | - | 103 | - | 155 | [1], [3, 4, 7] | 207 | - |  |
52 | - | 104 | - | 156 | - | 208 | - |  |
Table 3: Grouped indices for every case in MC.

##  Appendix B Multilevel percolation analysis

Figure 20: $Q-$criterion applied to the case S_1 is shown here. (a) The
highlighted blue region shows a single structure. A small domain from this
structure as seen in (b) is subjected to MLP and the result can be seen in
(c).

Figure 21: A schematic for MLP is shown here. $\tau_{p}$ is the percolation
threshold, $\langle D_{\alpha}\rangle$ is the mean fractal dimension,
$V_{max}$ is the volume of the largest structure and $V$ is the sum of the
volume of all structures in the domain.

In this section, we introduce an extension of the percolation analysis
technique described in section 3.3 to specifically deal with highly stratified
flows. Since global intermittency causes large nonturbulent patches to develop
within the flow, a global percolation threshold $\tau_{p}$ is not sufficient
to educe individual structures. For instance, consider the case S_1 as shown
in Figure 20(a). Here we are visualizing with the $Q-$criterion indicator and
thresholded at $\tau_{p}$. If one proceeds with $\tau_{p}$, the entire region
highlighted in blue will be extracted as a single structure and geometrically
characterized. We can observe in Figure 20(b) that this large blue structure
is actually composed of numerous smaller structures. In order to appropriately
educe these structures, we propose a novel technique called multilevel
percolation (MLP) analysis where percolation analysis is applied in a
repetitive manner on all structures at $\tau_{p}$ until all large clusters are
broken down into individual structures. The procedure for MLP is presented in
the schematic 21. It involves the following steps:

Figure 22: (a, b) show the behavior of the structure when the minimum value of
the ratio $V_{max}/V$ is chosen as $1$ and $0.5$ respectively. (c) shows a
proof-of-concept for MLP applied to a subset of the data from case S_1. Three
points are highlighted here corresponding to the following thresholds
$\tau=0.0625,0.33455,3.01098$.

  1. (a)

The global percolation threshold $\tau_{p}$ is computed with the procedure
described in section 3.3.

  2. (b)

All structures are extracted. For each structure percolation analysis is
repeated and a new value of $\tau_{p}$ is obtained for every structure.

  3. (c)

At the new $\tau_{p}$, all structures are extracted and fractal dimension is
computed with the method described in section 3.4 (a).

  4. (d)

Structures with $\langle D_{\alpha}\rangle<1$ are removed. This step is
crucial to obtain a noise-free result.

  5. (e)

Finally, the filtered structures are checked against a stopping criterion.
This criterion applies percolation analyis on the structure again and checks
for the minimum value of the ratio $V_{max}/V$ for the entire range of
threshold values. If the minimum value of the ratio $V_{max}/V$ is always
above $0.5$, we classify it as a simple structure and it is subjected to no
further processing. Essentially the new $\tau_{p}$ calculated in step (c) is
associated with the structure i.e. this structure exists at the new
$\tau_{p}$. On the other hand, if the minimum value of the ratio $V_{max}/V$
falls below $0.5$ at some threshold, we denote it as a complex structure. A
new $\tau_{p}$ is computed and the entire procedure is repeated from step (c).

  6. (f)

This process is repeated on all structures until all complex clusters are
broken into individual strutcures.

Rationale for the stopping criterion:

Adopting a minimum value of the ratio $V_{max}/V$ over the entire threshold
range is chosen as a stopping criterion due to its simplicity. One can extract
completely individual structures by setting the minimum value of the ratio
$V_{max}/V=1$. This condition ensures that at any increasing threshold the
volume of the largest structure is always equal to the total volume of all
structures i.e., there is exactly one structure (see Figure 22(a)).

We note that setting $V_{max}/V=1$ is not always practical due to the large
computation time requirement. Therefore, we identify an alternate minimum
value of $0.5$ for our stopping criterion. When $V_{max}/V>0.5$, at most two
structures can exist where one has a larger volume than the other (see Figure
22(b)). Two structures with exactly the same volume exist when
$V_{max}/V=0.5$. The choice of a minimum value other than $1$ is subjective
and must be considered only for larger or complex datasets.

A proof-of-concept for a simple and complex structure is illustrated in Figure
22(c). MLP is applied to a small region of the blue structure shown in Figure
20(a). The size of the domain is $300\times 90\times 600$. The minimum value
for the ratio was chosen as $0.5$. In its current implementation, the
computation took approximately 5 days and the result is shown in Figure 22(c).

## References

  * Garratt (1994) John Roy Garratt.  The atmospheric boundary layer.  _Earth-Science Reviews_ , 37(1-2):89–134, 1994\.
  * Van de Wiel et al. (2012) BJH Van de Wiel, AF Moene, HJJ Jonker, P Baas, S Basu, JMM Donda, Jielun Sun, and AAM Holtslag.  The minimum wind speed for sustainable turbulence in the nocturnal boundary layer.  _Journal of the Atmospheric Sciences_ , 69(11):3116–3127, 2012.
  * Mahrt (1989) Larry Mahrt.  Intermittency of atmospheric turbulence.  _Journal of the Atmospheric Sciences_ , 46(1):79–95, 1989.
  * Townsend (1948) AA Townsend.  Local isotropy in the turbulent wake of a cylinder.  _Australian Journal of Chemistry_ , 1(2):161–174, 1948.
  * Townsend (1949) AA Townsend.  The fully developed wake of a circular cylinder.  _Australian Journal of Chemistry_ , 2(4):451–468, 1949.
  * Corrsin and Kistler (1955) Stanley Corrsin and Alan L Kistler.  Free-stream boundaries of turbulent flows.  1955\.
  * Fiedler and Head (1966) H Fiedler and MR Head.  Intermittency measurements in the turbulent boundary layer.  _Journal of Fluid Mechanics_ , 25(4):719–735, 1966.
  * Kovasznay et al. (1970) Leslie SG Kovasznay, Valdis Kibens, and Ron F Blackwelder.  Large-scale motion in the intermittent region of a turbulent boundary layer.  _Journal of Fluid Mechanics_ , 41(2):283–325, 1970.
  * Rao et al. (1971) K Narahari Rao, R Narasimha, and MA Badri Narayanan.  The ‘bursting’phenomenon in a turbulent boundary layer.  _Journal of Fluid Mechanics_ , 48(2):339–352, 1971.
  * Falco (1977) RE Falco.  Coherent motions in the outer region of turbulent boundary layers.  _The Physics of Fluids_ , 20(10):S124–S132, 1977\.
  * Da Silva et al. (2014) Carlos B Da Silva, Rodrigo R Taveira, and Guillem Borrell.  Characteristics of the turbulent/nonturbulent interface in boundary layers, jets and shear-free turbulence.  In _Journal of Physics: Conference Series_ , volume 506, page 012015\. IOP Publishing, 2014.
  * Lugt (1979) Hans J Lugt.  The dilemma of defining a vortex.  In _Recent developments in theoretical and experimental fluid mechanics_ , pages 309–321. Springer, 1979.
  * Ansorge and Mellado (2016) Cedrick Ansorge and Juan Pedro Mellado.  Analyses of external and global intermittency in the logarithmic layer of ekman flow.  _Journal of Fluid Mechanics_ , 805:611–635, 2016.
  * Hunt et al. (1988) Julian CR Hunt, Alan A Wray, and Parviz Moin.  Eddies, streams, and convergence zones in turbulent flows.  1988\.
  * Robinson (1991) Stephen Kern Robinson.  The kinematics of turbulent boundary layer structure.  _NASA STI/Recon Technical Report N_ , 91, 1991.
  * Katul et al. (1997) Gabriel Katul, Greg Kuhn, John Schieldge, and Cheng-I Hsieh.  The ejection-sweep character of scalar fluxes in the unstable surface layer.  _Boundary-Layer Meteorology_ , 83(1):1–26, 1997\.
  * Li and Bou-Zeid (2011) Dan Li and Elie Bou-Zeid.  Coherent structures and the dissimilarity of turbulent transport of momentum and scalars in the unstable atmospheric surface layer.  _Boundary-layer meteorology_ , 140(2):243–262, 2011.
  * Hommema and Adrian (2003) Scott E Hommema and Ronald J Adrian.  Packet structure of surface eddies in the atmospheric boundary layer.  _Boundary-Layer Meteorology_ , 106(1):147–170, 2003.
  * Carper and Porté-Agel (2004) Matthew A Carper and Fernando Porté-Agel.  The role of coherent structures in subfilter-scale dissipation of turbulence measured in the atmospheric surface layer.  _Journal of Turbulence_ , 5(1):040, 2004.
  * Huang et al. (2009) J Huang, M Cassiani, and JD Albertson.  Analysis of coherent structures within the atmospheric boundary layer.  _Boundary-layer meteorology_ , 131(2):147–171, 2009.
  * Inagaki and Kanda (2010) Atsushi Inagaki and Manabu Kanda.  Organized structure of active turbulence over an array of cubes within the logarithmic layer of atmospheric flow.  _Boundary-layer meteorology_ , 135(2):209–228, 2010.
  * Moisy and Jiménez (2004) Frédéric Moisy and Javier Jiménez.  Geometry and clustering of intense structures in isotropic turbulence.  _Journal of fluid mechanics_ , 513:111–133, 2004.
  * Bermejo-Moreno and Pullin (2008) Ivan Bermejo-Moreno and DI Pullin.  On the non-local geometry of turbulence.  _Journal of Fluid Mechanics_ , 603:101–135, 2008.
  * Bermejo-Moreno et al. (2009) Ivan Bermejo-Moreno, DI Pullin, and Kiyosi Horiuti.  Geometry of enstrophy and dissipation, grid resolution effects and proximity issues in turbulence.  _Journal of fluid mechanics_ , 620:121–166, 2009.
  * Yang et al. (2010) Yue Yang, DI Pullin, and Ivan Bermejo-Moreno.  Multi-scale geometric analysis of lagrangian structures in isotropic turbulence.  _Journal of fluid mechanics_ , 654:233–270, 2010.
  * Yang and Pullin (2011) Yue Yang and DI Pullin.  Geometric study of lagrangian and eulerian structures in turbulent channel flow.  _Journal of Fluid Mechanics_ , 674:67–92, 2011.
  * Leung et al. (2012) T Leung, N Swaminathan, and PA Davidson.  Geometry and interaction of structures in homogeneous isotropic turbulence.  _Journal of Fluid Mechanics_ , 710:453–481, 2012.
  * Coleman et al. (1990) Gary Neil Coleman, JH Ferziger, and PR Spalart.  A numerical study of the turbulent ekman layer.  _Journal of Fluid Mechanics_ , 213:313–348, 1990.
  * Ansorge and Mellado (2014) Cedrick Ansorge and Juan Pedro Mellado.  Global intermittency and collapsing turbulence in the stratified planetary boundary layer.  _Boundary-layer meteorology_ , 153(1):89–116, 2014.
  * Ansorge (2016) Cedrick Ansorge.  _Analyses of turbulence in the neutrally and stably stratified planetary boundary layer_.  Springer, 2016.
  * Pope (2001) Stephen B Pope.  Turbulent flows, 2001.
  * Kline et al. (1967) Stephen J Kline, William C Reynolds, FA Schraub, and PW Runstadler.  The structure of turbulent boundary layers.  _Journal of Fluid Mechanics_ , 30(4):741–773, 1967.
  * Kline (1978) SJ Kline.  The role of visualization in the study of the structure of the turbulent boundary layer.  _Coherent Structure of Turbulent Boundary Layers_ , pages 1–26, 1978\.
  * Khanna and Brasseur (1998) Samir Khanna and James G Brasseur.  Three-dimensional buoyancy-and shear-induced local structure of the atmospheric boundary layer.  _Journal of Atmospheric Sciences_ , 55(5):710–743, 1998.
  * Drobinski and Foster (2003) Philippe Drobinski and Ralph C Foster.  On the origin of near-surface streaks in the neutrally-stratified planetary boundary layer.  _Boundary-layer meteorology_ , 108(2):247–256, 2003.
  * Jayaraman and Brasseur (2021) Balaji Jayaraman and James G Brasseur.  Transition in atmospheric boundary layer turbulence structure from neutral to convective, and large-scale rolls.  _Journal of Fluid Mechanics_ , 913, 2021.
  * Kim et al. (1971) HoT Kim, SJ Kline, and WC Reynolds.  The production of turbulence near a smooth wall in a turbulent boundary layer.  _Journal of Fluid Mechanics_ , 50(1):133–160, 1971.
  * Corino and Brodkey (1969) Edward Robert Corino and Robert S Brodkey.  A visual investigation of the wall region in turbulent flow.  _Journal of Fluid Mechanics_ , 37(1):1–30, 1969\.
  * Wallace et al. (1972) James M Wallace, Helmut Eckelmann, and Robert S Brodkey.  The wall region in turbulent shear flow.  _Journal of Fluid Mechanics_ , 54(1):39–48, 1972\.
  * Katul et al. (2006) Gabriel Katul, Davide Poggi, Daniela Cava, and John Finnigan.  The relative importance of ejections and sweeps to momentum transfer in the atmospheric boundary layer.  _Boundary-layer meteorology_ , 120(3):367–375, 2006.
  * Narasimha et al. (2007) Roddam Narasimha, S Rudra Kumar, A Prabhu, and SV Kailas.  Turbulent flux events in a nearly neutral atmospheric boundary layer.  _Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 365(1852):841–858, 2007.
  * Küchemann (1965) D Küchemann.  Report on the iutam symposium on concentrated vortex motions in fluids.  _Journal of Fluid Mechanics_ , 21(1):1–20, 1965\.
  * Cucitore et al. (1999) R Cucitore, Maurizio Quadrio, and Arturo Baron.  On the effectiveness and limitations of local criteria for the identification of a vortex.  _European Journal of Mechanics-B/Fluids_ , 18(2):261–282, 1999.
  * Günther and Theisel (2018) Tobias Günther and Holger Theisel.  The state of the art in vortex extraction.  In _Computer Graphics Forum_ , volume 37, pages 149–173. Wiley Online Library, 2018.
  * Chakraborty et al. (2005) Pinaki Chakraborty, S Balachandar, and Ronald J Adrian.  On the relationships between local vortex identification schemes.  _Journal of fluid mechanics_ , 535:189, 2005.
  * Jeong and Hussain (1995) Jinhee Jeong and Fazle Hussain.  On the identification of a vortex.  _Journal of fluid mechanics_ , 285:69–94, 1995.
  * Chong et al. (1990) Min S Chong, Anthony E Perry, and Brian J Cantwell.  A general classification of three-dimensional flow fields.  _Physics of Fluids A: Fluid Dynamics_ , 2(5):765–777, 1990.
  * Johansson et al. (1987) AV Johansson, PH Alfredsson, and J Kim.  Shear-layer structures in near-wall turbulence.  1987\.
  * Alfredsson et al. (1988) P Henrik Alfredsson, Arne V Johansson, and John Kim.  Turbulence production near walls: the role of flow structures with spanwise asymmetry.  _Studying Turbulence Using Numerical Simulation Databases-I1_ , 89:131, 1988.
  * Falco (1980) R Falco.  The production of turbulence near a wall.  In _13th Fluid and PlasmaDynamics Conference_ , page 1356, 1980.
  * Kim et al. (1987) John Kim, Parviz Moin, and Robert Moser.  Turbulence statistics in fully developed channel flow at low reynolds number.  _Journal of fluid mechanics_ , 177:133–166, 1987.
  * Chu and Falco (1988) CC Chu and RE Falco.  Vortex ring/viscous wall layer interaction model of the turbulence production process near walls.  _Experiments in Fluids_ , 6(5):305–315, 1988\.
  * Marusic et al. (2010) I Marusic, R Mathis, and N Hutchins.  Predictive model for wall-bounded turbulent flow.  _Science_ , 329(5988):193–196, 2010.
  * Shah and Bou-Zeid (2014) Stimit Shah and Elie Bou-Zeid.  Very-large-scale motions in the atmospheric boundary layer educed by snapshot proper orthogonal decomposition.  _Boundary-layer meteorology_ , 153(3):355–387, 2014.
  * Katul (2019) GG Katul.  The anatomy of large-scale motion in atmospheric boundary layers.  _Journal of Fluid Mechanics_ , 858:1–4, 2019.
  * Lorensen and Cline (1987) William E Lorensen and Harvey E Cline.  Marching cubes: A high resolution 3d surface construction algorithm.  _ACM siggraph computer graphics_ , 21(4):163–169, 1987.
  * Green et al. (2007) Melissa A Green, Clarence W Rowley, and George Haller.  Detection of lagrangian coherent structures in three-dimensional turbulence.  _Journal of Fluid Mechanics_ , 572:111, 2007.
  * Del Álamo et al. (2006) Juan C Del Álamo, Javier Jimenez, Paulo Zandonade, and Robert D Moser.  Self-similar vortex clusters in the turbulent logarithmic region.  _Journal of Fluid Mechanics_ , 561:329–358, 2006.
  * Nagaosa and Handler (2003) Ryuichi Nagaosa and Robert A Handler.  Statistical analysis of coherent vortices near a free surface in a fully developed turbulence.  _Physics of Fluids_ , 15(2):375–394, 2003.
  * Koenderink and Van Doorn (1992) Jan J Koenderink and Andrea J Van Doorn.  Surface shape and curvature scales.  _Image and vision computing_ , 10(8):557–564, 1992.
  * Do Carmo (2016) Manfredo P Do Carmo.  _Differential geometry of curves and surfaces: revised and updated second edition_.  Courier Dover Publications, 2016.
  * Zhang and Chen (2001) Cha Zhang and Tsuhan Chen.  Efficient feature extraction for 2d/3d objects in mesh representation.  In _Proceedings 2001 International Conference on Image Processing (Cat. No. 01CH37205)_ , volume 3, pages 935–938. IEEE, 2001.
  * Bogard and Tiederman (1986) DG Bogard and WG Tiederman.  Burst detection with single-point velocity measurements.  _Journal of Fluid Mechanics_ , 162:389–413, 1986.
  * Cabral and Leedom (1993) Brian Cabral and Leith Casey Leedom.  Imaging vector fields using line integral convolution.  In _Proceedings of the 20th annual conference on Computer graphics and interactive techniques_ , pages 263–270, 1993.
  * Willmarth and Lu (1972) WW Willmarth and SS Lu.  Structure of the reynolds stress near the wall.  _Journal of Fluid Mechanics_ , 55(1):65–92, 1972\.
  * Chen and Blackwelder (1978) Chyan-Hai P Chen and Ron F Blackwelder.  Large-scale motion in a turbulent boundary layer: a study using temperature contamination.  _Journal of Fluid Mechanics_ , 89(1):1–31, 1978\.
  * Harikrishnan et al. (2020) A Harikrishnan, C Ansorge, R Klein, and N Vercauteren.  The curious nature of hairpin vortices. gallery of fluid motion.  2020\.
  * Lewiner et al. (2003) Thomas Lewiner, Hélio Lopes, Antônio Wilson Vieira, and Geovan Tavares.  Efficient implementation of marching cubes’ cases with topological guarantees.  _Journal of graphics tools_ , 8(2):1–15, 2003\.
  * Faranda et al. (2017) Davide Faranda, Gabriele Messori, and Pascal Yiou.  Dynamical proxies of north atlantic predictability and extremes.  _Scientific reports_ , 7(1):1–10, 2017.
  * Stalling et al. (2005) Detlev Stalling, Malte Westerhoff, and Hans-Christian Hege.  Amira: a highly interactive system for visual data analysis.  In _The Visualization Handbook_ , pages 749–767. Elsevier, 2005.