1 Introduction
Recent advances in the field of network analysis have revealed the structures of internet, technological, social and biological networks (Albert & Barabási Reference Albert and Barabási2002; Newman Reference Newman2003, Reference Newman2010; Barrat et al. Reference Barrat, Barthélemy, Pastor-Satorras and Vespignani2004a ). Having characterized these networks, we are able to study dynamics such as disease outbreak and information propagation on networks, and analyse the resilience of network-based activities (Albert, Jeong & Barabási Reference Albert, Jeong and Barabási2000; Barrat, Barthélemy & Vespignani Reference Barrat, Barthélemy and Vespignani2008). These analysis techniques are founded on graph theory, dynamical systems, and operator theory, but place unique emphasis on interactions and connectivity amongst the elements that establish a network. Thus far, most of the applications of network analysis have been concerned with discrete settings in which nodes are individualized quantities, such as people, organisms, equipment, or stations (Caldarelli Reference Caldarelli2007; Newman Reference Newman2010). In this paper, we extend the network analysis to continuous representation of physical phenomena, in particular two-dimensional turbulence.
The chaotic motion of a large number of vortices in turbulent flows is caused by the induced velocities of the vortices themselves. What makes turbulence rich and complex are the vortical interactions in the flow field, which take place over a wide range of length scales (Tennekes & Lumley Reference Tennekes and Lumley1972; Hinze Reference Hinze1975; Frisch Reference Frisch1995; Pope Reference Pope2000; Davidson Reference Davidson2004; Lesieur Reference Lesieur2008). Thus, complete understanding of turbulence has remained a challenge to this day because of its high-dimensionality, multiscale interactions, nonlinearity and the resulting chaos. Network science provides an alternative view of complex fluid flows in terms of a network of vortex interactions (Nair & Taira Reference Nair and Taira2015), and this perspective illuminates the underlying structure and organization of turbulent flows. In this work, we show that two-dimensional isotropic turbulence (Kraichnan & Montgomery Reference Kraichnan and Montgomery1980; McWilliams Reference McWilliams1984; Benzi, Paladin & Vulpiani Reference Benzi, Paladin and Vulpiani1990; Benzi & Colella Reference Benzi and Colella1992; Davidson Reference Davidson2004; Boffetta & Ecke Reference Boffetta and Ecke2012) has a scale-free network structure reminiscent of other networks found in nature (Barabási & Albert Reference Barabási and Albert1999; Caldarelli Reference Caldarelli2007). While most of the attention has been placed on unweighted scale-free networks, we consider the use of weighted scale-free networks to describe the variations in the strength of interactions or connectivities (Barrat, Barthélemy & Vespignani Reference Barrat, Barthélemy and Vespignani2004b ). Upon identifying the network structure of turbulence, physical insights can be obtained as to which vortical interactions are important in capturing the overall physics and how it may be possible to control the dynamics of turbulent vortices (Farazmand, Kevlahan & Protas Reference Farazmand, Kevlahan and Protas2011; Liu, Slotine & Barabási Reference Liu, Slotine and Barabási2011; Brunton & Noack Reference Brunton and Noack2015).
2 Problem description and approach
To extract the network structure of the flow, we quantify the interactions between fluid elements based on the vortical interactions. The velocity
$\boldsymbol{u}$
at position
$\boldsymbol{x}$
induced by the vorticity distribution
${\bf\omega}$
of the flow is

In this study, we focus on unforced two-dimensional isotropic turbulence in a periodic box and assess the influence of the vorticity distribution over a Cartesian domain. Here, the two-dimensional vorticity field reduces to
${\bf\omega}(\boldsymbol{x},t)={\it\omega}(x,y,t)\hat{\boldsymbol{e}}_{z}$
, with
$\hat{\boldsymbol{e}}_{z}$
denoting the unit normal plane vector. Modelling the vortical component for each discrete Cartesian element as a line vortex, we can evaluate how fluid elements influence each other, as depicted in figure 1. Here, the magnitude of the induced velocity from fluid element
$i$
on another element
$j$
reduces from (2.1) to

where
${\it\gamma}_{i}={\it\omega}(\boldsymbol{x}_{\boldsymbol{i}}){\rm\Delta}x{\rm\Delta}y$
is the circulation of fluid element
$i$
with side lengths of
${\rm\Delta}x$
and
${\rm\Delta}y$
. The superposition of the induced velocity from all other fluid elements provides the advective velocity of the fluid element. Detailed discussions on using point vortices to develop the network-theoretic framework for describing unsteady vortical flows can be found in Nair & Taira (Reference Nair and Taira2015). Note that adjacency matrices are commonly defined with positive weights as considered here, but they can be relaxed to accommodate positive and negative weights within the context of vortical interactions. This point will be revisited later.

Figure 1. Interaction of fluid elements in two-dimensional turbulence. The strength of the vortical interaction between elements
$i$
and
$j$
having vorticity
${\it\omega}_{i}$
and
${\it\omega}_{j}$
is quantified through the induced velocities
$u_{i\rightarrow j}$
and
$u_{j\rightarrow i}$
, respectively. For discretizing the Cartesian domain, we take
$n_{x}$
and
$n_{y}$
points in the horizontal and vertical directions, respectively, providing the adjacency matrix
$\unicode[STIX]{x1D63C}$
of size
$n\times n$
with
$n=n_{x}n_{y}$
. Shown in the background with a contour plot is the corresponding vorticity field with initial
$Re(t_{0})=814$
at
$t=18$
.
To assess and describe the vortical interactions in the flow field, we utilize a weighted network (graph). The definition of a network (graph)
$\mathscr{G}$
requires sets of vertices (nodes)
$\mathscr{V}$
, edges
$\mathscr{E}$
, and weights
$\mathscr{W}$
(Newman Reference Newman2010). With these three components defined, a graph can be uniquely determined, i.e.
$\mathscr{G}=\mathscr{G}(\mathscr{V},\mathscr{E},\mathscr{W})$
. The nodes
$\mathscr{V}$
in this study are taken to be the vortical elements residing within the Cartesian cells, while the edges
$\mathscr{E}$
represent the vortical interactions between those vortical elements. The edge weights
$\mathscr{W}$
quantify the strengths of the vortical interactions. Given
$n$
nodes, a collection of the weights
$w_{ij}$
in the form of a matrix
$\unicode[STIX]{x1D63C}\in \mathbb{R}^{n\times n}$
with

is called the adjacency matrix and is used to describe the network connectivity. In the above definition,
$\unicode[STIX]{x1D608}_{ij}$
is set to the edge weight
$w_{ij}$
if there exists an edge (interaction) between nodes
$i$
and
$j$
. Details on the fundamental concepts involved in network theory can be found in Dorogovtsev (Reference Dorogovtsev2010) and Newman (Reference Newman2010), with descriptions of vortical-interaction networks in Nair & Taira (Reference Nair and Taira2015).
Based on (2.2), we define the network adjacency matrix as the average induced velocity

to quantify the magnitude of interaction between fluid elements
$i$
and
$j$
(Nair & Taira Reference Nair and Taira2015). Note that an element cannot impose velocity upon itself, which is captured by the null entry along the diagonal of the adjacency matrix. In the present study, the influence from the neighbouring periodic vortex images are also accounted for in the analysis. This formulation yields a full matrix, except for its diagonal entries, which are identically zero. In assessing the strength of the vortical interaction between two fluid elements, we utilize (2.4) to perform network analysis to extract the spatial connectivity structure. This approach has been successful in capturing the nonlinear vortex dynamics and modelling the trajectories of vortex clusters (Nair & Taira Reference Nair and Taira2015). The adjacency matrix considered here is symmetric to quantify the average interaction strength. Note that the geometric mean can be alternatively chosen and yields similar results. In general, the adjacency matrix can be formulated in an asymmetric manner:

Here the parameter
${\it\phi}$
takes a value between 0 and 1. For the aforementioned symmetric formulation in (2.4),
${\it\phi}$
is selected as
$1/2$
. When
${\it\phi}=0$
and 1, the adjacency matrix
$\unicode[STIX]{x1D608}_{ij}$
are defined by the velocity imposed to the other elements (
$\unicode[STIX]{x1D608}_{ij}=u_{j\rightarrow i}$
) and upon themselves (
$\unicode[STIX]{x1D608}_{ij}=u_{i\rightarrow j}$
), respectively, for
$i\neq j$
. We mainly focus on the use of the symmetric adjacency matrix in this work, but will consider the asymmetric formulation briefly to highlight the difference from a physical point of view in the next section. We note in passing that the theoretical tools for symmetric adjacency matrices are more widely available compared to the asymmetric matrices.
The flow field analysed in this study is obtained from direct numerical simulation on a square biperiodic computational domain
$(x,y)\in [0,L]\times [0,L]$
with a grid size of
$m_{x}\times m_{y}=1024\times 1024$
. The unforced two-dimensional incompressible isotropic turbulent flow is simulated by numerically solving the two-dimensional vorticity transport equation

where
$\boldsymbol{u}$
and
${\it\omega}$
are the velocity and vorticity variables, respectively. The simulation is performed with the Fourier spectral method and the fourth-order Runge–Kutta time integration scheme (Canuto et al.
Reference Canuto, Hussaini, Quarteroni and Zang1988). The vorticity field is initialized with a smooth distribution comprised of a large number (
${\approx}100$
) of superposed vortices (Taylor Reference Taylor1918) with random strengths, core sizes, and locations chosen such that the kinetic energy spectra satisfies
$E(k)\propto k\exp (-k^{2}/k_{0}^{2})$
, where
$k_{0}=26.5$
, following the set-up by Kida (Reference Kida1985) and Brachet, Meneguzzi & Sulem (Reference Brachet, Meneguzzi and Sulem1986). The initial core sizes are selected to be sufficiently small compared to the size of the computational domain (McWilliams Reference McWilliams1984) arranged in random positions. The velocity variable is normalized by the square root of the spatial average of the initial kinetic energy
$u^{\ast }(t_{0})\equiv [\overline{u^{2}}(t_{0})]^{1/2}$
, where the overline denotes the spatial average. The spatial length and time scales are non-dimensionalized by the initial integral length scale
$l^{\ast }(t_{0})\equiv [2\overline{u^{2}}(t_{0})/\overline{{\it\omega}^{2}}(t_{0})]^{1/2}$
and the initial eddy turnover time
$t_{0}^{\ast }\equiv l^{\ast }(t_{0})/u^{\ast }(t_{0})$
, respectively. The Reynolds number is defined accordingly as
$Re\equiv u^{\ast }l^{\ast }/{\it\nu}$
, where
${\it\nu}$
is the kinematic viscosity. In this study, turbulent flows with initial Reynolds numbers of
$Re(t_{0})=75$
, 439, 814, 1607 and 2485 are selected.
3 Results
3.1 Network-based characterization
We identify the underlying network structure and characteristics of two-dimensional turbulence based on the aforementioned symmetric adjacency weights. The time-evolving vorticity field is obtained from a two-dimensional incompressible biperiodic direct numerical simulation (Canuto et al.
Reference Canuto, Hussaini, Quarteroni and Zang1988) for unforced isotropic turbulence. Given the vorticity field over a Cartesian grid, each fluid element is considered to be connected to all other elements through vortical network edges. The resulting fluid flow network can in fact be described by a complete graph with a range of weights. Next, we visualize the network edges with a transparent grey scale corresponding to the adjacency weight, as shown in figure 2(a). The captured structure reveals the turbulent network. Some regions in the flow have a large number of strong connections, corresponding to larger stronger vortices seen in red, serving as primary network hubs. Note that these strong vortices induce velocities over long distances. Moderate size vortices that act as secondary hubs also possess dominant connections to primary hubs and other secondary hubs. In contrast, fluid elements corresponding to smaller, weaker eddies, shown in blue, generally have influence only in their vicinity. The node strength distribution (
$s_{i}=\sum _{j}\unicode[STIX]{x1D608}_{ij}$
) over space shows that the vortices with large circulation have larger strength, as illustrated in figure 2(b). The node strength distribution over space enables us to distinguish secondary and primary hubs, which may not be easily differentiated from simply visualizing the vorticity field or the
$Q$
criterion in a traditional manner. For instance, see the green vortices in figure 2(b), which can appear similar to primary ones in vorticity level.

Figure 2. The scale-free network of vortical interactions in two-dimensional turbulence with initial
$Re(t_{0})=814$
. (a) Turbulent network structure overlaid on the vorticity field, with the darkness of the network edges corresponding to the values of the adjacency weights (
$t=18$
). (b) Contour plot of the node strength
$s$
distribution. Vortex cores having high degree of connectivity act as hubs in the turbulent vortical network. (c) The corresponding node strength probability distribution, exhibiting the scale-free characteristics with
$P\sim s^{-2.7}$
. The same contour level is shared by (b,c). Also shown in the background of (c) in grey are the out and in-degree distributions (
${\it\phi}=0$
and 1, respectively). The network visualized in (a) does not show interactions from periodic images, and uses
$32\times 32$
nodes for graphical clarity.
Plotting the probability of the strength distribution
$P(s)$
over the strength
$s$
of fluid elements in figure 2(c), we find that two-dimensional isotropic turbulence network has a power-law distribution
$P(s)\sim s^{-{\it\gamma}}$
, with
${\it\gamma}=2.7$
, at the time shown. This tells us that the vortex interactions in turbulence can be characterized by a weighted scale-free network. This realization enables the interaction-based analysis of turbulent flows from a new perspective through network theory (Cohen & Havlin Reference Cohen and Havlin2010; Newman Reference Newman2010). In particular, this type of network is known to have certain resilience properties, as we will explore later in this section. Also shown in figure 2(c) in grey are the degree distributions for asymmetric adjacency formulations. The out- and in-degree distributions can be found by setting
${\it\phi}=0$
and 1, respectively, in (2.5). It can be observed that the scale-free symmetric distribution is mostly comprised of the out-degree components, which describe how each vortical element influences all other elements (i.e.
$u_{j\rightarrow i}$
). In contrast, we find that the in-degree distribution has a single peak, which conveys that all fluid elements receive a similar amount of collective influence from vortices in the flow field. We have found that the scale-free property of two-dimensional isotropic turbulence is most well captured by the symmetric weights compared to the other asymmetric formulations. It is also possible to examine the strength distribution taking positive and negative values of circulations, as we have briefly discussed in § 2. Utilizing positive and negative weights, their strength distribution can also exhibit a scale-free behaviour, but with network strength having both negative and positive values. This leads to a symmetric strength distribution over the strength, with a resemblance to the probability density function of scaled displacements (Weiss, Provenzale & McWilliams Reference Weiss, Provenzale and McWilliams1998). In what follows, results based on the symmetric adjacency matrix (using the magnitude of induced velocity) are presented.
Let us further examine the time-varying properties of the turbulent network. In unforced turbulence, the kinetic energy of the flow decreases over time due to viscous dissipation, as shown in figure 3(a). The strength distribution
$P(s)$
of the turbulence network and the corresponding flow field snapshots are presented in figure 3(b). Turbulent flow is comprised of vortical structures over a wide range of spatial scales initially. The distribution
$P(s)$
exhibits scale-free characteristics with
$P(s)\sim s^{-{\it\gamma}}$
, where
${\it\gamma}\approx 2.7$
, when the kinetic energy spectra exhibits the
$k^{-3}$
profile for
$t\lesssim 30$
. For the flow under consideration, a bend in the strength distribution appears for
$t\gtrsim 30$
as the system starts to exhibit scale separation and loses the
$k^{-3}$
energy spectra. This is caused by the diffusion of smaller-scale structures and their merging with other structures. Over time, viscous dissipation removes kinetic energy through the smaller eddies and leaves only the larger vortices. This behaviour can be described by two power laws,
$P(s)\sim s^{-{\it\gamma}_{1}}$
and
$P(s)\sim s^{-{\it\gamma}_{2}}$
, which capture the weaker fluid elements and the larger stronger vortices, respectively. The bifurcation of these power laws is shown in figure 3(c), indicated by the vertical dashed line. We note that, regardless of the initial condition used, the turbulent interaction network maintains the scale-free behaviour in the present investigation as long as the energy spectra relaxes to the
$k^{-3}$
profile (Kida Reference Kida1985; Brachet et al.
Reference Brachet, Meneguzzi and Sulem1986; Benzi et al.
Reference Benzi, Paladin and Vulpiani1990). This scale-free behaviour may be observed during the initial transient, but is not a guaranteed common feature without the
$k^{-3}$
energy spectra being present.

Figure 3. The dynamics of a turbulent network with
$Re(t_{0})=814$
. (a) Kinetic energy and (b) strength distribution of two-dimensional isotropic turbulence for
$t=15$
, 30, 75, 150 and 300 (line colours represent time). The inset plots in (b) show the corresponding vorticity fields. The kinetic energy
$E(k)$
is shown over the wavenumber
$k$
, exhibiting the asymptotic profile of
$E(k)\sim k^{-3}$
. The strength distribution
$P(s)$
displays the scale-free property of
$P(s)\sim s^{-{\it\gamma}}$
over node strength
$s$
. (c) The corresponding exponents
${\it\gamma}$
,
${\it\gamma}_{1}$
and
${\it\gamma}_{2}$
are shown. Later in time the strength distribution exhibits the emergence of two distributions,
$P(s)\sim s^{-{\it\gamma}_{1}}$
and
$s^{-{\it\gamma}_{2}}$
.
We have considered a range of Reynolds numbers and observed that
${\it\gamma}$
takes values of
${\it\gamma}=2.7\pm 0.5$
. The variations observed in
${\it\gamma}$
,
${\it\gamma}_{1}$
and
${\it\gamma}_{2}$
shown in figure 3(c) are influenced by the chaotic nature of turbulence. These parameters, however, appear to exhibit a coalescing behaviour when they are plotted over the product of the characteristic velocity and length,
$u^{\ast }(t)l^{\ast }(t)$
. Here, we interpret
$u^{\ast }(t)l^{\ast }(t)$
as the circulation of vortices that have the characteristic velocity and length scales. As shown in figure 4, we observe that the turbulence network shows coalescence of the scale-free parameter
${\it\gamma}$
to
${\it\gamma}_{\mathit{cr}}\approx 2.7$
over time for different cases of turbulent flows. Once the flows reach a state where the characteristic strength of vortices is
$(u^{\ast }l^{\ast })_{\mathit{cr}}\approx 0.063$
, the network distribution bifurcates to display two different slopes with
${\it\gamma}_{1}$
and
${\it\gamma}_{2}$
, as previously illustrated in figure 3. This observation reveals that a scale-free turbulent network is present until the unforced turbulent flow field loses the smaller-scale vortices, and mostly contains vortices with strengths larger than the critical value of
$(u^{\ast }l^{\ast })_{\mathit{cr}}$
.

Figure 4. Exponent
${\it\gamma}$
for the network strength distribution
$P(s)\sim s^{-{\it\gamma}}$
plotted over
$u^{\ast }(t)l^{\ast }(t)$
with different initial Reynolds numbers (green:
$Re(t_{0})=75$
, purple:
$Re(t_{0})=439$
, yellow:
$Re(t_{0})=814$
, red:
$Re(t_{0})=1607$
, and blue:
$Re(t_{0})=2485$
). Scale-free distributions are observed with
${\it\gamma}$
coalescing to
${\it\gamma}_{\mathit{cr}}\approx 2.7$
up until a bifurcation at
$(u^{\ast }l^{\ast })_{\mathit{cr}}\approx 0.063$
. The grey box shows
${\it\gamma}=2.7\pm 0.5$
as reference.
3.2 Resilience of turbulence networks
Characterizing turbulent flow with a scale-free network enables us to view turbulent interactions in a systematic manner, and provides insights into how vortical structures influence each other. It is known from network analysis that scale-free networks are resilient to random perturbations, but attacks towards network hubs can affect network dynamics in a detrimental manner (Albert et al. Reference Albert, Jeong and Barabási2000). Network resilience for fluid flow translates to the difficulty of modifying the vortical-interaction network and, consequently, the collective behaviour of the vortices over time. To measure the change in vortical interaction caused by network disturbance, we can consider how the removal of turbulence network nodes (percolation) modifies the characteristic network length

which is the average shortest network distance
$d(i,j)$
between any two nodes on a network. Here, we perform node percolation by setting the vorticity values at the chosen nodes to be zero. The above metric quantifies how well vortical elements are connected within a turbulent network. Note that the distance here refers to network distance based on the adjacency matrix, and not the spatial distance. In particular, we take the inverse of each adjacency weight
$1/a_{ij}$
and evaluate the minimal sum

over a network path that connects nodes
$i$
and
$j$
for this metric (Rubinov & Sporns Reference Rubinov and Sporns2010). This metric
$l_{\mathit{network}}$
can be thought of as the average of the minimal characteristic advective (commute) time per unit length between every pair of fluid elements in the domain. This minimal network distance is determined using the Floyd–Warshall algorithm (Floyd Reference Floyd1962).

Figure 5. The resilience of turbulence network against fraction of node removals,
$f$
, for
$t=15$
, 30 and 75, with
$Re(t_{0})=814$
. Shown are relative changes in the characteristic network length
${\rm\Delta}\tilde{l}_{\mathit{network}}$
of turbulent flow for random node and hub node removals. The colours of the curves represent the time when node removal is considered, and follows figure 3.
The changes in the turbulence network characteristic length
$l_{\mathit{network}}$
when network nodes are removed in a random fashion and a coordinated manner targeting hub nodes are summarized in figure 5. Here, the changes in the normalized characteristic network length

for varied fraction of node removal
$f$
are shown. While it would be difficult to completely remove nodes, as we have performed in this investigation, the present analysis sheds light on how external forcing or perturbations can alter the turbulent flow from an interaction-based analysis. We observe that turbulent flow is resilient against random forcing, as evident from the characteristic network length being unaffected even for a large fraction
$f$
of nodes being removed. This behaviour is consistently observed over time. On the other hand, we find that the global vortical-interaction network can be greatly modified by targeting large vortex cores (hubs), as exhibited by the substantial change in the characteristic length. It may be more energetically expensive to remove well-connected hub nodes, which often correspond to regions of concentrated vorticity. However, it is clear from figure 5 that even the smallest fraction of hub node removal can influence the overall interaction, which suggests that hub removal still provides a more effective and efficient way to modify the flow than random node removal.
When the vortical-interaction network is grossly altered, the dynamics of the collection of vortices would be significantly modified (Nair & Taira Reference Nair and Taira2015). These observations also agree with past studies in flow control that identified effective actuation frequencies to be associated with the length scale of the large coherent structures in turbulent flows (Gad-el-Hak Reference Gad-el-Hak2000; Joslin & Miller Reference Joslin and Miller2009). With increasing time, we can further notice that network connectivity decreases with hub removal, due to viscous dissipation of smaller vortical structures, and the influence of removing the core structures becomes more evident. The present network-based understanding reveals which types of flow structure should be targeted with flow control if we aim to alter the behaviour of the turbulent flow field in a global manner.
4 Concluding remarks
The approach presented in this paper is the initial effort in performing network-based analysis of complex turbulent flows. Using the mathematical toolsets from network theory, we have identified that the vortical interactions in two-dimensional decaying isotropic turbulence have a scale-free network structure. We have been able to reveal the structure by taking a continuous representation of the flow field and quantifying the network using a Cartesian discretization. For two-dimensional isotropic turbulence, the node strength distribution was uncovered to be
$P(s)\sim s^{-{\it\gamma}}$
, where
${\it\gamma}=2.7\pm 0.5$
. Furthermore, we have found that the unforced turbulent flow field possesses an underlying scale-free network structure until the circulation of vortices with characteristic velocity and length scales reach
$(u^{\ast }l^{\ast })_{\mathit{cr}}\approx 0.063$
. By noticing that the turbulence network has scale-free characteristics, we were able to systematically show that the turbulence network is resilient against random perturbations, but vulnerable against coordinated forcing on the hub vortices. It should be noted that estimating and controlling each and every vortical structure in a turbulent flow is most likely improbable and impractical. Instead, network analysis may provide a refreshing view point on how one can predict and modify the collective dynamics of vortices in the turbulent flow fields. We believe that network-based analysis and control (Mesbahi & Egerstedt Reference Mesbahi and Egerstedt2010; Liu et al.
Reference Liu, Slotine and Barabási2011; Cornelius, Kath & Motter Reference Cornelius, Kath and Motter2013; Kaiser et al.
Reference Kaiser, Noack, Cordier, Spohn, Segond, Abel, Daviller, Osth, Krajnović and Niven2014; Yan et al.
Reference Yan, Tsekenis, Barzel, Slotine, Liu and Barabási2015) will provide a novel mathematical fabric for paving the path towards network-based modelling and control of turbulent flows, which can potentially impact a wide spectrum of problems.
Acknowledgements
K.T. and A.G.N. acknowledge the support from the US Army Research Office (grant: W911NF-14-1-0386, Program manager: Dr S. Stanton) and the US Air Force Office of Scientific Research (grant: FA9550-13-1-0183, Program manager Dr D. Smith). S.L.B. acknowledges support by the Department of Mechanical Engineering at the University of Washington and as a Data Science Fellow in the eScience Institute.