∎

11institutetext: J.S. Cheng 22institutetext: Department of Mechanical
Engineering
University of Rochester
Rochester, NY 14627
22email<EMAIL_ADDRESS>33institutetext: D. H. Kelley
44institutetext: Department of Mechanical Engineering
University of Rochester
Rochester, NY 14627
44email<EMAIL_ADDRESS>
# Laboratory model of electrovortex flow with thermal gradients, for liquid
metal batteries

J.S. Cheng  B. Wang  I. Mohammad  J.M. Forer  D.H. Kelley

(Received: date / Accepted: date)

###### Abstract

We present a novel laboratory setup for studying the fluid dynamics in liquid
metal batteries (LMBs). LMBs are a promising technology suited for grid-scale
energy storage, but flows remain a confounding factor in determining their
viability. Two important drivers of flow are thermal gradients, caused by
internal heating during operation, and electrovortex flow (EVF), induced by
diverging current densities. Our setup explores thermal gradients and
electrovortex flow separately and in combination in a cylindrical layer of
liquid gallium, simulating the behavior in a single layer of an LMB. In this
work, we discuss the design principles underlying our choices of materials,
thermal control, and current control. We also detail our diagnostic tools —
thermocouple measurements for temperature and Ultrasonic Doppler Velocimetry
(UDV) probes for velocities — and the design principles which go into choosing
their placement on the setup. We also include a discussion of our post-
processing tools for quantifying and visualizing the flow. Finally, we
validate convection and EVF in our setup: we show that scaling relationships
between the nondimensional parameters produced by our data agree well with
theory and previous studies.

###### Keywords:

Energy storage Liquid metals Convective heat transfer Magnetohydrodynamics

††journal: myjournal

##  1 Introduction

The electricity generated by renewable energy sources almost doubled over the
past decade in the United States (USe, 2021) and is likely to continue to rise
in the near future. However, the intermittent nature of renewable sources such
as solar and wind presents a challenge for widespread integration into the
electrical grid: the demand for electrical power can vary by up to 50% in a
given region of the US over the course of a day (Hodge, 2020), but this
variation rarely lines up favorably with peak renewable production. To
facilitate the adoption of renewable sources, grid scale energy storage
solutions are needed. The liquid metal battery (LMB) is one promising
technology due to its low cost and long life (Wang et al., 2014; Perez and
Kelley, 2015).

LMBs are electrochemical cells that consist of two liquid metal electrodes
separated by a layer of molten salt electrolyte, and with all three layers
stably stratified. The more electropositive metal, an electron donor, is the
negative electrode (anode), while the more electronegative metal, an electron
acceptor, is the positive electrode (cathode). The interaction between the two
electrodes creates a thermodynamical cell voltage. During discharge, the
negative electrode metal is oxidized into cations and electrons. The molten
salt electrolyte only allows cations through to the positive electrode metal
forcing the electrons through an external circuit thus providing electrical
power (Kim et al., 2013a). The reverse happens during charging where
electrical power is absorbed. Much of the cost effectiveness of LMBs is owed
to the liquid nature of the battery: the layers naturally separate by density
which leads to simpler construction and improved ability to scale to larger or
smaller cells, the liquid electrodes are immune to the microscale damage
suffered by solid electrodes during charge and discharge, and the candidate
materials can be relatively inexpensive (Bradwell et al., 2012; Kim et al.,
2013a; Wang et al., 2014; Li et al., 2016; Dai et al., 2018). LMBs are a
developing technology, with various container geometries and chemical
compositions under exploration (Bradwell et al., 2012; Poizeau et al., 2012;
Kim et al., 2013b; Ning et al., 2015; Gong et al., 2020).

Figure 1: (a) Schematic demonstrating thermal gradients and electromagnetic
forces in a liquid metal battery and (b) our experimental setup for studying
these forces and their interactions. The diagram in panel (a) shows the
variation in density $\rho$ with vertical position $z$, demonstrating that the
cathode, electrolyte, and anode layers are separated by stable stratification.
Each layer also contains density gradients due to heating in the electrolyte
layer. The experimental setup can model either the cathode or anode of the
battery, depending on whether we impose a stable or unstable temperature
gradient. The velocity field shown in panel (b) is for illustration only and
was adapted from Personnettaz et al. (2019).

While most earlier analyses focused on the electrochemistry and economic
viability, recent years have brought attention to the fluid dynamics occurring
within the battery during operation (Kelley and Weier, 2018; Weber et al.,
2018). The fluid nature of LMBs lends them a strong propensity for flow: a
variety of thermal, chemical, and electromagnetic forces play against one
another with serious consequences, both helpful and harmful, for the battery.
Helpfully, flow can prevent buildup of intermetallics, a primary source of
failure in LMB test cells. During discharge, flow redistributes ion-saturated
parcels of cathode material away from the interface faster than diffusion and
can thus improve efficiency. Harmfully, sufficiently vigorous flows can
disturb the interfaces between different fluid layers and cause a short,
destroying the battery. The fluid dynamics in LMBs can be paramount to their
operation and efficiency. It is thus important to develop a broad
understanding of how various flow drivers interact with one another via
scaling laws between governing parameters that can be generalized.

Here, we present a novel laboratory device for exploring the interaction
between two prominent drivers of fluid motion in LMBs – thermal gradients and
electro-vortex flow (EVF) – over broad ranges of the governing parameters. By
covering a large swath of parameter space, we can construct scaling laws and
regime diagrams to extrapolate results to real battery systems. We apply these
forces to a layer of liquid gallium to simulate behaviors in a single LMB
layer, collecting temperature data using thermocouples and velocity data using
Ultrasonic Doppler Velocimetry (UDV) probes.

In this paper we will discuss the design principles going into the
construction of the experimental setup, laboratory techniques for optimal
acquisition of data, and novel post-processing techniques for interpreting
experimental results. In section 2, we describe the framework of flow physics
we intend to capture with our device and how these lead to our overarching
design choices. In section 3, we delve into detail about the design,
explaining the specifics of how we control the flow with precision and ensure
the setup is robust. In section 4, we describe the expected behavior under the
imposed flow drivers and how our diagnostics are arranged to perform
meaningful measurements. In section 5, we validate our design by showing early
data and comparison with predictions from the literature. In section 6, we
provide concluding remarks.

Table 1: Estimated material properties of liquid gallium at 43∘C (Okada and Ozoe, 1992; Brito et al., 2001; Aurnou et al., 2018) compared to ranges of properties for LMB candidate materials at their melting temperatures (Kelley and Weier, 2018). These properties are comparable, particularly with the cathode, though it should be noted that broad ranges of values are possible in LMBs. Anode materials considered are Li, Mg, and Na (Davidson, 1968; Iida and Guthrie, 2015a, b) while cathode materials considered are Bi, Pb, Sb, and Zn (Iida and Guthrie, 2015a, b; Fazio et al., 2015). Properties listed are mass density $\rho$, kinematic viscosity $\nu$, thermal diffusivity $\kappa$, volumetric expansion coefficient $\alpha$, Prandtl number $Pr$, electrical conductivity $\sigma$, and speed of sound $c$. Material |  $\rho$ [kg/m3] |  $\nu$ [m2/s] |  $\kappa$ [m2/s] |  $\alpha$ [1/K] | $Pr$ |  $\sigma$ [S/m]
---|---|---|---|---|---|---
Gallium | 5870 | $3.588\times 10^{-7}$ | $1.345\times 10^{-5}$ | $1.25\times 10^{-4}$ | 0.027 | $3.88\times 10^{6}$
Anode | 500–1600 | 7.5–12$\times 10^{-7}$ | 2–7$\times 10^{-5}$ | 1.6–2.5$\times 10^{-4}$ | 0.01–0.06 | 3.6–10$\times 10^{6}$
Cathode | 6500–10000 | 1.6–5.3$\times 10^{-7}$ | 1.0–1.6$\times 10^{-5}$ | 1.2–1.5$\times 10^{-4}$ | 0.01–0.03 | 0.8–2.7$\times 10^{6}$

##  2 Design Motivation

The LMB environment is host to many flow drivers and interactions. We choose
to focus on two prominent contributors to fluid motion: thermal gradients and
electrovortex flow (Bojarevičs et al., 1989; Davidson, 2001; Kelley and
Sadoway, 2014; Shen and Zikanov, 2016; Personnettaz et al., 2018; Weber et
al., 2020; Keogh et al., 2021). Figure 1a is a schematic of a liquid metal
battery during operation, demonstrating how each flow manifests. Thermal
gradients across the layer cause stable stratification or convection, while a
diverging current density across the layer causes electrovortex flow. We
describe the physics behind each flow below in sections 2.1 and 2.2.

Our experimental setup — shown in Figure 1b — simulates the behavior in the
cathode or the anode of the LMB by producing thermal gradients and
electrovortex flow together in an enclosed fluid layer. We focus on each layer
independently, precluding interfacial flows that make it difficult to pinpoint
the source of a given observed motion (Köllner et al., 2017). We chose liquid
gallium as the primary working fluid: it offers similar properties to LMB
candidate materials, but near room temperature, rather than at hundreds of
degrees Celcius (see Table 1). Of particular relevance to the fluid dynamics
is the Prandtl number $Pr=\nu/\kappa$, where $\nu$ is the kinematic
diffusivity and $\kappa$ is the thermal diffusivity, which tends to have
values between 0.01 and 0.03 for LMB layers and has value 0.027 in our
experimental setup. Operating around room temperature affords us a greater
degree of thermal control, since water cooling and heating become viable
(discussed further in Section 3.1).

We study cylindrical LMBs to avoid introducing another free parameter in the
aspect ratio between length versus width. Cylinders are also the most common
shape for models in the literature (Ning et al., 2015; Li et al., 2016;
Beltrán, 2017; Weber et al., 2018). Our cell diameter $D=10$ cm is chosen to
achieve relevant values of the nondimensional parameters (detailed below in
sections 2.1 and 2.2) and ensuring scalability to LMBs while also keeping
experiments manageable in time and length scale. The depths of the cathode and
anode layers change during battery operation, as material passes from one
layer to the other. Different cylindrical aspect ratios $\Gamma=D/H$ (where
$H$ is height) can contain significantly different flows (Wagner and
Shishkina, 2013; Vogt et al., 2018; Wang et al., 2020), and these flows may
vary in their effect on LMB mixing. We employ interchangeable side wall
segments ranging from height $H=3.3-7.1$ cm ($\Gamma=D/H=\sqrt{2}-3$) to
address this.

###  2.1 Thermal gradients

In liquid metal batteries, thermal gradients occur because of the high
electrical resistivity of the electrolyte layer — typically several orders of
magnitude more resistive than the other layers (Janz et al., 1968). The
electrolyte preferentially heats up as cations flow through it during both
charge and discharge of the battery. The anode experiences a destabilizing
thermal gradient that leads to convection, whose forcing strength is
characterized by the nondimensional Rayleigh number

$Ra=\frac{\alpha g\Delta TH^{3}}{\nu\kappa}\,.$ (1)

Here, $\alpha$ is the coefficient of thermal expansion, $g$ is gravitational
acceleration, and $\Delta T$ is the magnitude of the temperature difference.
Beyond a critical value $Ra\sim 2000$ (Chandrasekhar, 1961) and prior to
$Ra\sim 10^{9}$ convective flows remain within a regime of dynamical
similarity (Cioni et al., 1997; Ahlers et al., 2009). Liberal estimates for
typical LMB temperature gradients and dimensions place them squarely within
this regime: even for an upper bound LMB size of 50 cm, temperature gradients
upwards of $\sim 100~{}^{\circ}$C would be required to reach a higher regime
of convection. Conversely, even for a lower bound LMB size of 5 cm, a
temperature gradients of order 0.1 ∘C is enough to reach the onset of the
convection regime we explore. With this in mind, our setup is designed to
access $Ra$ values spanning two orders of magnitude, from $10^{4}$ to
$10^{6}$. Since $Ra$ scales linearly with $\Delta T$, we need detailed
temperature measurements to characterize $Ra$. We discuss these in Section 4.

One way to quantify the resultant flows is with the Reynolds number
$Re=UL/\nu$, where $L$ is a typical length scale and $U$ is a typical flow
speed. The Reynolds number serves as an overall descriptor of system dynamics,
with literature predictions for its scaling in both convection and EVF. In
convection, $Re$ is expected to follow a power law against $Ra$ (Ahlers et
al., 2009):

$Re\propto Ra^{2/5}Pr^{-3/5}\,,$ (2)

for the parameter ranges relevant to LMBs. Velocities therefore scale with the
temperature gradient, and $Ra$ can provide insight on whether the flow is in
the range of helpful mixing or strong enough to disrupt fluid interfaces.

We use this to validate our setup in section 5. However, the Reynolds number
alone does not tell us the morphology of the flow, which can vary depending on
flow forcing and geometry of the system. Velocity measurements that capture
the flow from various angles and across different spatial positions are
therefore a priority. More details about our velocity diagnostics are
discussed in section 4.

We note that recent studies estimate that compositional gradients will
dominate thermal gradients in the cathode during charge and discharge
(Personnettaz et al., 2019; Herreman et al., 2020). The nondimensional
parameters for solutal convection are structurally similar to thermal
convection, with a solutal Prandtl number $Pr_{s}=\nu/\mathfrak{d}$, analogous
to the thermal Prandtl number $Pr$ ($\mathfrak{d}$ is the mass diffusion
coefficient). This similarity means our experimental setup may also provide
insight into compositional gradients when water ($Pr\simeq 7$) is the working
fluid instead of gallium.

###  2.2 Electrovortex flow

Electrovortex flow occurs when a current density diverges radially (Shercliff,
1970), as is the case in almost any liquid metal battery configuration (shown
by the blue arrows in Figure 1a). The induced magnetic field from a given
current line interacts with nearby currents and generates Lorentz forces which
drive flow radially toward the smaller current collector. The effect
strengthens closer to the current collector, resulting in a downward-pointing
jet under the current collector, and induced flow throughout the rest of the
fluid layer. The strength of EVF is characterized by a so-called EVF parameter
(Lundquist, 1969):

$S=\frac{\mu_{0}I^{2}}{4\pi^{2}\rho\nu^{2}}$ (3)

where $\mu_{0}$ is the vacuum permeability and $I$ is the imposed current. An
especially vigorous jet has potential to disrupt the interfaces between LMB
layers. A well-documented transition toward an azimuthal swirling flow also
occurs (Bojarevičs et al., 1989; Davidson, 2001), which we describe in more
detail in section 4. The swirling flow may be more well-suited to mixing
within the cathode with less risk of interface disruption (Weber et al.,
2020). Ashour et al. (2018) observed this transition at a current of 20 A or
$S=1.7\times 10^{4}$ in an experiment using PbBi eutectic. Simulations by
Herreman et al. (2019) find different velocity scalings on either side of this
transition, with $Re\sim S$ for $S\lesssim 10^{3}$ and $Re\sim S^{1/2}$ for
$S\gtrsim 10^{5}$. To ensure we can explore the transition, our setup is
designed to access $S$ values from 0 to $3\times 10^{5}$.

The combination of EVF and convection has been studied in the past. Kelley and
Sadoway Kelley and Sadoway (2014) found that introducing EVF in a convecting
fluid causes the flow to become more ordered even at small input currents, an
observation corroborated by Ashour et al. (2018) and Keogh et al. (2021).
Predictions also exist for the transition between convection-dominated and
EVF-dominated flow: Davidson (2001) derived a criterion where convection
begins to dominate when $\left(RaPr\right)^{3/7}\gtrsim S^{1/2}Pr$, although
this is for a hemispherical geometry rather than cylindrical. Ashour et al.
(2018) argued that a transition occurs when the characteristic velocities for
each forcing become comparable. Our accessible $Ra$ and $S$ ranges ensure that
we can test these predictions over broad parameter ranges.

During LMB operation, a stabilizing thermal gradient manifests in the cathode
layer that may suppress other sources of fluid motion. To our knowledge,
stable thermal stratification combined with EVF has not been previously
explored although Herreman et al. (2020) examine EVF with solutal
stratification. We therefore make it possible to reverse the thermal gradient
in our setup to study this problem. Data produced by our apparatus will be a
source of novel information about the LMB cathode.

Figure 2: Current and thermal gradient controls for the experimental setup.
Heat exchangers are composed of a copper plate and copper coil. Water is
recirculated through the coils by a thermal bath and a chiller, which can be
interchanged to induce stable or unstable stratification. The current pathway
is drawn in blue. The bottom thermal block also acts a current collector. The
wires connected to the current collectors are routed vertically for a large
distance such that the magnetic field induced by the current does not
interfere with the flow.

##  3 Experimental design

###  3.1 Controls

Figure 2 is a diagram of the experimental setup and its associated external
components for controlling temperature and current. The gallium layer is
bounded above and below by copper plates, each plate soldered to a copper coil
which serves as a heat exchanger: two thermal baths recirculate water through
each of these coils such that either a stabilizing or destabilizing
temperature gradient can be imposed on the gallium layer. The coils are wound
in a non-inductive double spiral pattern to minimize horizontal temperature
gradients while accommodating instruments and the negative electrode.

Figure 3: Coating the thermal blocks where they contact gallium. (a), A (3 cm
$\times$ 5 cm) sheet of copper corroded by contact with liquid gallium over
the course of two days. (b), Before CuO coating. (c), After CuO coating. (d),
The boundary between coated and uncoated regions on a piece of Si, coated
using the same procedure and imaged via electron microscopy, which showed good
coverage. (e–f), Energy-dispersive x-ray spectroscopy indicated nearly pure Si
in the uncoated region and CuO in the coated region.

Copper is used for the heat exchangers due to its high thermal and electrical
conductivity. However, liquid gallium attacks and damages copper, as shown in
Figure 3a. Such corrosion can be counteracted by a protective coating;
following Zhang et al. (2007), we choose copper oxide for the coating layer. A
thermal evaporation method (Ozer and Tepehan, 1993) was performed using a Ladd
Evaporator (Ladd Research Industries, U.S.A). During the evaporation process,
copper oxide pellets (CuO, 99.9% purity) were placed in a tungsten basket and
evaporated in the vacuum chamber under a base pressure of $5\times 10^{-4}$
Pa, depositing copper oxide particles on the heat exchanger suspended above.
The resulting coating is shown in Figure 3c, and the X-ray spectrum in Figure
3d demonstrates that the process was successful. In situ experiments confirm
that this layer prevents gallium corrosion. Based on the evaporation rate, we
estimate the thickness of the coating layer as 0.6 µm. A 1-D Fourier’s law
approximation (assuming the copper plate and coating layer as resistors in
series) finds the layer’s effect on heat transfer to be negligible, with only
a 0.15% change in $\Delta T$ compared to a pure copper heat exchanger.

To ensure temperature uniformity across each plate, we must account for the
impacting and detaching of convective plumes on the surface: these should not
create significant temperature gradients within the plate. We address this
concern by considering an estimate of the Biot number for convection

$Bi=\frac{q}{\Delta
T}\frac{h_{\mathrm{plate}}}{k_{\mathrm{plate}}}=\frac{k}{k_{\mathrm{plate}}}\frac{h_{\mathrm{plate}}}{H}Nu,$
(4)

where $h_{\mathrm{plate}}$ is the thickness of the plate, $k_{\mathrm{plate}}$
is the thermal conductivity of the plate, $k$ is the thermal conductivity of
gallium, and the Nusselt number $Nu$ is a nondimensional parameter
representing the convective heat transfer efficiency. Typically, $Bi<0.1$
means that the temperature inside the plate remains uniform. Assuming that
$Nu$ scales with $Ra$ as shown in the literature (Ahlers et al., 2009), we
estimate that a plate thickness $h_{\mathrm{plate}}<0.8$ cm would give uniform
temperature.

We must also ensure that the temperature at the surface in contact with the
gallium layer is insensitive to the structure of the cooling coil. We address
this by creating a simplified, two-dimensional model of heat transfer —
considering a cross-section of the copper plate with the coil as a series of
hot and cold patches — and finding the steady-state temperature distribution
at the plate-gallium interface. See the Appendix for a more thorough
description. We found that when the plate thickness is 0.5 cm, in all cases
the temperature deviation is $<2$%. Ultimately, 0.5 cm was chosen as the plate
thickness to minimize $Bi$, maximize temperature uniformity, and also optimize
machining efficiency.

The side wall of the vessel is a 1.1 cm thick hollow cylinder made of Delrin
(polyoxymethylene), a thermally resistive material ($k_{\mathrm{Delrin}}=0.37$
W/mK). This ensures that heat flows almost exclusively through the gallium
layer: using a 1-D Fourier’s law approximation (assuming the Delrin sidewall
and gallium layer as resistors in parallel) yields that less than 0.6% of the
applied heat flux travels through the sidewalls for all experiments (and is a
notable overestimate because convection has not been taken into account).

The entire experimental setup is surrounded by a layer of foam rubber
insulation that is at least 1 inch thick to minimize heat losses, ensuring
that nearly all the heat introduced by the thermal bath is extracted by the
chiller, and maximizing the range of achievable temperature differences
$\Delta T$.

Figure 4:  The copper coil and plate which form the top thermal block, seen
from the side (a) and from the top (b). Dimensions are given in inches. We
minimize spacing between coils and ensure non-inductive winding while
accommodating the UDV probes, filling ports, and top electrode. The dashed
line marks the outside of the vessel side wall while the gray region indicates
the position of the gallium layer.

We induce currents of 0 to 90 A through the system with a TDK-Lambda power
supply. Delrin’s low electrical conductivity (0.2 W/mK) ensures all currents
pass through the gallium layer instead of traveling along the side wall. The
negative electrode is a copper rod of multiple possible diameters (0.318 cm,
0.476 cm, 0.635 cm) attached to the top plate by a plastic fitting,
electrically isolating it from the top plate. The bottom plate is mated to a
copper rod connected to the power supply: with the low electrical resistivity
of copper — roughly a factor of 15 lower than liquid gallium (Ginter et al.,
1986) — the current spreads out uniformly and the entire bottom plate forms
the positive electrode. The large diameter ratio between electrodes ensures a
large horizontal divergence in the current paths to trigger EVF.

Electrovortex flow is sensitive to external magnetic fields (Davidson, 2001).
To minimize the induced field from horizontal currents at the gallium layer,
we ensure that the currents flowing in and out of the setup are oriented only
vertically in the vicinity of the gallium: the vertical rods attached to the
top and bottom plates for carrying current are at least 0.5 m long (see Fig.
2). Using the Biot-Savart law, we estimate that the induced lateral magnetic
field strength is $<20$ µT at the fluid layer (smaller than Earth’s magnetic
field).

##  4 Diagnostics

Temperature data are collected by K-type thermocouple probes placed on the top
and bottom heat exchangers. They are distributed around the plates at various
radial distances and separated by $90^{\circ}$ in angle (find distances).
Small blind holes are drilled into the thermal plates to assign reproducible
probe locations and to bring the tips of the probes into closer vicinity of
the fluid layer (within 3 mm). Signals are digitized and recorded by National
Instruments NI 9211 cDAQ thermocouple input modules, operated by a LabVIEW
acquisition program. Temperature measurements show that the spread in
temperatures between probes on the same plate is typically less than 10% of
the temperature difference across the gallium layer in all cases, and usually
close to 5%. Furthermore, the spread in temperatures on a given plate is only
double the standard deviation of an individual thermocouple probe. We believe
this demonstrates a successful degree of temperature uniformity over each heat
exchanger.

Liquid gallium is opaque, which greatly limits our velocity measurement
options: particle tracking, Schlieren, and dye based techniques are not
viable. UDV is therefore our method of choice. In this technique, pulsed
ultrasound waves are emitted by transducers, fluid provide an instantaneous
velocity profile along the sound propagation direction. This technique is non-
invasive and has been validated for use with various liquid metals (Takeda,
1987; Brito et al., 2001; Eckert et al., 2007).

For gallium experiments, we do not introduce any tracer particles. Naturally
existing gallium oxide inclusions are capable of scattering the ultrasound
waves. Due to the density difference between liquid gallium and gallium oxide,
though, scattering particles tend to settle over time, reducing scattering
strength. The wetting and acoustic coupling conditions between the UDV
transducer and gallium are also observed to deteriorate over time (Wang and
Kelley, 2021). We therefore limit the duration of each set of UDV measurements
to 1 hour for optimal signal quality.

The working principle of UDV has been introduced in detail in our previous
work (Perez and Kelley, 2015). Here, we use nine ultrasound transducers
(Signal Processing, Switzerland) to measure flow velocities, each of which
measure the component of velocity towards and away from the probe. Placements
and orientation are shown in Figure 5. Because of the sparse, ‘gappy’ nature
of UDV data, we must carefully choose probe locations to best capture the
expected flows. We also have to take into account the geometric limitations
imposed by the setup: each probe is held in place by a compression fitting of
finite size (SwageLok 8 mm), limiting the density with which probes can be
arranged. For probes in direct contact with the gallium, the fitting seals
around the probe and provides a flush surface with the inner walls of the
vessel, but also requires extra clearance. Our UDV measurements are validated
via comparison to particle tracking velocimetry measurements — see Wang and
Kelley (2021) for details.

Figure 5:  UDV probe locations on the 50 mm high setup (other tank heights use
analogous locations), seen from the side (a) and from above (b). Each probe
measures one component of velocity along the probe beam. Probes 5,6 and 7 are
placed at chord positions and angled 5∘ away from the y-direction, such that
their beams are oriented in the y-direction after passing through the sidewall
(ie. parallel to probes 1,2 and 3).

Convection and EVF are both expected to take the form of overturning rolls in
the accessible parameter ranges and geometries (Cioni et al., 1997; Davidson,
2001). In convection, a single convection roll is expected to fill the
geometry whereas in EVF, a vertical jet descends along the centerline inducing
rising flow along the side walls (Davidson, 2001; Ashour et al., 2018). We
therefore place two transducers with working frequency 8 MHz (probes 8 and 9)
at the top lid and in direct contact with gallium for measuring vertical
velocities at different radial positions. One is 4 mm from the side wall
($x=0.92D/2$) and the other is at the center or 12.5 mm from the center
($x=0.25D/2$), depending on the case (the center position is occupied by the
negative electrode during EVF cases, but is available during convection-only
cases).

Figure 6: UDV measurements of convection at $Ra=3.7\times 10^{5}$ in the
$\Gamma=1.7$ setup. (a), Velocity measured by probe 2 (see Fig. 5), varying
over space and time. (b), Time-averaged velocity varying with distance from
probe 2. The green highlighted region is used for $Re$ calculations to produce
Figure 8. (c), Spatially averaged mean velocity varying over time. (d), The
corresponding power spectrum has a prominent peak at 0.033 Hz.

Four transducers with working frequency 8 MHz are placed at the side wall of
the vessel, oriented in the radial direction and in direct contact with the
gallium. Convection rolls are likely to span the whole vessel in the vertical
direction in our $\Gamma>1$ geometries, and strong horizontal flows should be
situated near the top and bottom of the tank. In addition, horizontal
velocities are expected to show vertical asymmetry when stable stratification
is introduced (Herreman et al., 2020). We therefore place horizontal probes at
$z=0.75H$ and $z=0.25H$. Since flows are not necessarily axisymmetric (Vogt et
al., 2018; Keogh et al., 2021), we place probe 4 on the x-z plane instead of
the y-z plane.

Flows with strong azimuthal components, such as the swirl flow in EVF, cannot
be detected by radial probes. We therefore use three horizontal probes
oriented along chord positions $\left(x=-0.25D,0.25D\right)$ and at heights
$z=\left(0.25H,0.5H,0.75H\right)$. These transducers (probes 5 through 7) have
a working frequency of 4 MHz and measure flow velocities through the wall of
the vessel (no direct contact with gallium). An 8 MHz ultrasound beam has a
smaller wavelength and thus gives better spatial resolution. A 4 MHz
ultrasound beam, on the other hand, has a lower attenuation rate and thus
allows stronger signals to be transmitted through the vessel wall. To maximize
the transmission of ultrasound signal, the area of the side wall where
transducers make contact is thinned to 0.3 mm.

The speed of sound in Delrin (2430 m/s) is slower than in gallium (2860 m/s),
and refraction occurs when the transducer beam passes into the fluid layer. In
order to keep the chord probe beams parallel to probes 1–3, we used Snell’s
law to determine that the probes should be tilted 5∘ away from the y-direction
(see Fig. 5b).

All transducers are connected to a DOP 3010 Velocimeter (Signal Processing,
Switzerland) and operated in emit/receive mode for data acquisition. The Pulse
Repetition Frequency (PRF) and number of emissions per profile together
determine the sampling frequency and measurement quality of UDV. During
measurements, those two parameters are adjusted according to the flow
conditions – for example, it is necessary to ensure that the sampling
frequency is always faster than the variation of the flow velocities. To
compensate for ultrasound attenuation in gallium, the Time Gate Control (TGC)
function was used. Rather than a uniform TGC value, we apply multiple custom-
defined TGC values at different distances along the beam in order to measure
flow velocities as close as possible to the transducer surface. Even so, a
small region close to the transducer surface (usually 5-15 mm) is unmeasurable
due to ringing in the transducer tip. When the current is turned on for EVF
measurements, the signal becomes noticeably noisier due to interference from
induced magnetic fields. However, high quality data are still possible.

Finally, for both thermocouples and UDV probes, we hope to capture all the
relevant spectral information, requiring that time between samples be shorter
than the typical flow timescales. For convection, the buoyancy turnover
timescale can be estimated from balancing buoyant and inertial terms in the
Navier-Stokes equation: $t_{\mathit{ff}}=\left(H/\alpha g\Delta
T\right)^{1/2}$ (Glazier et al., 1999). For electrovortex flow, the velocity
scaling from Herreman et al. (2019) can be rearranged to give
$t_{\mathrm{EVF}}\sim H/S$ for $S\lesssim 10^{3}$ and $t_{\mathrm{EVF}}\sim
H/S^{1/2}$ for $S\gtrsim 10^{5}$. For our most vigorous achievable flows, we
estimate a minimum timescale of 1.2 seconds. Both thermal and velocity
measurements are acquired at multiple samples per second.

Figure 7:  a) Time-averaged UDV measurements from multiple probes in a
convection case at $Ra=3.7\times 10^{5}$ and $\Gamma=1.7$. Colors correspond
to the probe positions denoted in Fig. 5. b) a frame of a movie produced from
the same data (Replace with time-averaged version instead). Velocity vectors
are represented by cones pointing in the direction of flow. Cone size and
color represent velocity magnitude. Cropping, smoothing, and interpolation in
time and space eliminate spurious vectors and noise. See Supplementary
Materials for a time-dependent movie in the style of panel b.

##  5 Typical measurements

Experimental measurements are shown in Figure 6. Noisy regions corresponding
to reflections from the back wall of the vessel have been truncated.
Measurements from each UDV probe are represented as a Dopplergram in panel a),
where velocities along the beam are plotted versus time and space. Panels b)
and c) are time- and space-averaged data from the Dopplergram, respectively,
and panel d) is the FFT spectrum, where a peak likely associated with the
recently characterized jump-rope vortex appears (Vogt et al., 2018). In panels
a) and b), gaps in the velocity observed near the probe surface are an
artifact of the probe tip ringing effect.

Raw UDV datasets as in Figure 6 do not lend themselves to easy visualization,
and their relationship with the actual flow field may not be obvious. In
Figure 7, we apply a novel post-processing tool to address this issue: time-
dependent velocity data from each probe are mapped to their physical position
on an analogue to the apparatus and recorded in a movie. Figure 7 shows the
time-averaged velocity field; see Supplementary Materials for an animation. We
apply a spatial moving average of 30 bins (7 mm) to reduce noise and account
for velocity gaps.

Figure 8:  a) Reynolds number $Re$ plotted versus Rayleigh number $Ra$ for a
set of convection cases in the $\Gamma=2$ setup. Results from Takeshita et al.
(1996); Cioni et al. (1997); Scheel and Schumacher (2017); Zürner et al.
(2019) are also included for comparison. The best-fit scaling, $Re\propto
Ra^{0.40}$, matches the exponent 0.4 predicted by theory in Eq. 2 (Ahlers et
al., 2009). b) Reynolds number $Re$ plotted versus EVF parameter $S$ for a set
of EVF cases in the $\Gamma=2$ setup. The best-fit scaling $Re\propto
S^{0.52}$ matches closely with the exponent $0.5$ predicted by theory
(Bojarevičs et al., 1989).

We examined how $Re$ scales with $Ra$ in order to validate that the convective
flows we drive have their expected properties, as shown in Figure 8a. We
performed a series of convection experiments, with no current, by setting
adverse temperature gradients of $1-20~{}^{\circ}$C between the thermal bath
and chiller. To obtain a precise value of $\Delta T$ for calculating $Ra$ via
Eq. 1, we averaged the temperature difference between top and bottom
thermocouples, over both space and time. To derive a $Re$ value from
velocities, we first observed from Figure 7 that a single large scale
circulation (LSC) appears to fill the container. We therefore adopted the
method of Zürner et al. (2019) to extract a characteristic LSC velocity from
UDV probes: the magnitude of the velocity at
$\left(x,y,z\right)=\left(0,0,0.75H\right)$ is
$U_{\mathrm{LSC}}=\langle\left(u_{3}^{2}+u_{4}^{2}\right)^{1/2}\rangle$, where
$u_{3}$ and $u_{4}$ are the velocities measured by probes 3 and 4,
respectively, and $\langle\rangle$ indicates time averaging. Unlike Zürner et
al. (2019), we did not average this velocity near the top with a velocity near
the bottom, since there is only one probe beam intersecting
$\left(x,y,z\right)=\left(0,0,0.25H\right)$. Figure 8a shows that our data are
well-fit by a $Re$–$Ra$ power law with exponent 0.40, in close agreement with
the theoretically predicted 0.4, and also in agreement with scaling results
from previous convection studies.

We validated EVF results through a similar procedure, as shown in Figure 8b.
We set a small stabilizing temperature gradient $\Delta T=-0.5~{}^{\circ}$C
and varied the current from 5 to 40 A. We calculated $S$ following Eq. 3. To
estimate $Re$, we modified the methodology of Herreman et al. (2019): the
characteristic EVF velocity was chosen by finding the peak velocity out of all
UDV measurements at each time step, and then averaging this value over time,
i.e. $U_{\mathrm{EVF}}=\langle u_{\mathrm{max}}\rangle$. We find that results
bear close agreement to the predicted $Re$–$S$ power law exponent of 0.5.
These fits to theoretical predictions and literature values suggest that our
thermal controls, current controls, and diagnostics are robust.

##  6 Conclusion

Liquid metal batteries are an evolving technology. Since it is unclear which
geometries and materials will ultimately have industrial relevance, our
understanding of flows in LMBs must therefore be general enough to accommodate
these possibilities. Our new experimental apparatus addresses this problem by
isolating two primary flow drivers — thermal gradients and electrovortex flow
— in order to predict their effects and interactions over a range of governing
parameters scalable to LMBs. Compared to a full LMB prototype, we can more
clearly pinpoint the causes for each observed flow behavior. This is aided by
treating a simplified cylindrical geometry, looking at each layer
independently, and conducting experiments near room temperature.

To produce well-characterized flows, we ensure that the container materials
provide the correct boundary conditions and that the control systems provide
precise temperatures and currents. We validate these systems with simplified
theoretical models and experimental tests. To adequately quantify the flows,
we inform the placement of temperature and velocity probes with predictions
from the literature. Using temperature and velocity measurements, we can
calculate governing parameter values and partially reconstruct the flow field.
We validate these measurements through comparison to previous results in the
literature.

Future studies in our setup will involve surveys of convection, EVF, combined
convection and EVF, and combined stable stratification and EVF. In combined
surveys, we will check existing predictions for the transition between EVF-
dominated and thermal-dominated flows (Davidson, 2001; Ashour et al., 2018).
Our study of EVF combined with stable thermal stratification will be, we
believe, the first of its kind.

While our simplified model has many advantages, care must also be taken when
applying our results to LMBs. Thermal gradients could also manifest in other
patterns, such as with horizontal components (Shishkina et al., 2016). The
Lorentz forces driving EVF occur in close proximity to the narrower electrode,
meaning geometrical or material differences in that region could have a
significant effect on the resultant flow. Finally, convection and EVF must be
weighed against the multitude of other possible flows. As mentioned in section
2, solutal convection can dominate thermal effects. Fluid-fluid interfaces can
serve as a source for many flows: for example, tilted interface surfaces
between fluids of different conductivity lead to the metal pad roll
instability (Sele, 1977; Zikanov, 2015; Horstmann et al., 2018); variations in
surface tension along interfaces lead to Marangoni flow (Köllner et al.,
2017). The relative strengths of each effect are not fully characterized, but
Kelley and Weier (2018) contains an extensive discussion about this topic.

Advancing the study of LMBs requires distinguishing the contributions of
myriad flows in a controlled setting. Our novel setup provides an important
step toward understanding the fluid behaviors in LMBs and thereby moving them
ever closer to commercial viability.

##  7 Appendix: Theoretical model for heat transfer inside heat exchanger
plate

We consider a two-dimensional radial cross-section of the copper plate placed
in an $(x,z)$ Cartesian coordinate system with $z$ pointing upwards, $x$ along
the diameter of the plate, and the origin at the midpoint of the plate’s
bottom. Since we are interested in the steady-state temperature distribution
of the plate, the governing equation is the Laplace equation:

$\nabla^{2}T(x,z)=0.$ (5)

The side walls of the plate are in contact with a heat insulator, and we
assume that no heat escapes through the sides, i.e., the heat flux through the
side walls is zero, which is asserted as the Neumann boundary condition:

$\frac{\partial T}{\partial x}(x=-\frac{D}{2};\frac{D}{2},z)=0.$ (6)

At the plate-gallium interface heating (or cooling) convection exists between
the plate boundary and the liquid gallium, and the heat transfer through the
boundary will be proportional to the temperature difference between the
gallium and the boundary. We use Newton’s law of cooling and find the Robin
boundary condition:

$\frac{\partial T}{\partial
z}(x,z=0)=\frac{-Bi}{h_{\mathrm{plate}}}(T(x,z=0)-T_{Ga}),$ (7)

where $T_{Ga}$ is the gallium’s mean temperature.

The coil is in contact with the last plate boundary and we represented its
thermal effect as an evenly spaced distribution of hot and cold spots added to
a set temperature, which can be stated as the Dirichlet boundary condition:

$T(x,z=h_{\mathrm{plate}})=T_{\mathrm{set}}+\delta T\cos\left(\frac{2\pi
x}{\lambda}\right),$ (8)

where $T_{\mathrm{set}}$ is the set temperature of the plate, $\delta T$ is
the temperature difference between a hot spot and a cold spot, and $\lambda$
is the distance between two hot or two cold spots.

We solve Eq. 5 with boundary conditions given by Eqs. 8, 6, and 7 numerically
using MATLAB PDEModel with a mesh of $\sim$4800 two-dimensional quadratic
triangular elements. We use the parameters $Bi=0.005$, $\lambda=2.5~{}cm$,
$D=10~{}cm$, $T_{Ga}=43~{}^{\circ}C$, and a range of plate thicknesses
$<0.6~{}cm$ with $T_{\mathrm{set}}$ values equivalent to the ones used in
experiment. We set $\delta T$ to 5% of the thermal gradient applied between
top and bottom $\Delta T$. We quantify the temperature variation at the
boundary in contact with gallium as the standard deviation of the temperature
there, reported as a percentage of $\Delta T$.

Figure 9:  a) Radial cross-section of the copper plate attached to heat
exchanger. b) Steady state temperature distribution inside the plate for
$\Delta T=12$. c) The temperature variation at the plate-gallium boundary for
different $\Delta T$ values where the plate thickness is $0.5$ cm. As shown,
the variation is $<2\%$ for all relevant $\Delta T$ cases.

###### Acknowledgements.

The authors wish to thank Mike Pomerantz and Jim Alkins for design
consultation and construction of the setup. The authors also thank Gerrit M.
Horstmann for productive discussions and the idea of writing this paper, as
well as Norbert Weber and Tobias Vogt for fruitful scientific discussions. We
thank Alex M. Grannan for inspiring the 2D model described in the Appendix,
which was built off of his PhD work. This work was supported by the National
Science Foundation under award number CBET-1552182.

## Conflict of interest

The authors declare that they have no conflict of interest.

## References

  * USe (2021) (2021) Electricity explained: Electricity in the United States. URL https://www.eia.gov/energyexplained/electricity/electricity-in-the-us.php
  * Ahlers et al. (2009) Ahlers G, Grossmann S, Lohse D (2009) Heat transfer and large scale dynamics in turbulent Rayleigh–Bénard convection. Rev Mod Phys 81:503–537
  * Ashour et al. (2018) Ashour RF, Kelley DH, Salas A, Starace M, Weber N, Weier T (2018) Competing forces in liquid metal electrodes and batteries. J Power Sources 378:301–310
  * Aurnou et al. (2018) Aurnou JM, Bertin V, Grannan AM, Horn S, Vogt T (2018) Rotating thermal convection in liquid gallium: multi-modal flow, absent steady columns. J Fluid Mech 846:846–876
  * Beltrán (2017) Beltrán A (2017) MHD natural convection flow in a liquid metal electrode. Appl Therm Eng 114:1203–1212
  * Bojarevičs et al. (1989) Bojarevičs V, Freibergs Y, Shilova EI, Shcherbinin EV (1989) Electrically Induced Vortical Flows. Kluwer Academic Publishers, Dordrecht, The Netherlands
  * Bradwell et al. (2012) Bradwell DJ, Kim H, Sirk AHC, Sadoway DR (2012) Magnesium–antimony liquid metal battery for stationary energy storage. J Am Chem Soc 134(4):1895–1897
  * Brito et al. (2001) Brito D, Nataf HC, Cardin P, Aubert J, Masson JP (2001) Ultrasonic Doppler velocimetry in liquid gallium. Expt Fluids 31(6):653–663
  * Chandrasekhar (1961) Chandrasekhar S (1961) Hydrodynamic and Hydromagnetic Stability, 1st edn. Oxford University Press
  * Cioni et al. (1997) Cioni S, Ciliberto S, Sommeria J (1997) Strongly turbulent Rayleigh–Bénard convection in mercury: comparison with results at moderate Prandtl number. J Fluid Mech 335:111–140
  * Dai et al. (2018) Dai T, Zhao Y, Ning XH, Narayan RL, Li J, Shan Zw (2018) Capacity extended bismuth–antimony cathode for high-performance liquid metal battery. J Power Sources 381:38–45
  * Davidson (1968) Davidson HW (1968) Compilation of thermophysical properties of liquid lithium. Report No. NASA TN-D-4650
  * Davidson (2001) Davidson PA (2001) An Introduction to Magnetohydrodynamics, 1st edn. Cambridge Texts in Applied Mathematics
  * Eckert et al. (2007) Eckert S, Cramer A, Gerbeth G (2007) Velocity measurement techniques for liquid metal flows. In: Magnetohydrodynamics, Springer, pp 275–294
  * Fazio et al. (2015) Fazio C, Sobolev VP, Aerts A, Gavrilov S, Lambrinou K, Schuurmans P, Gessi A, Agostini P, Ciampichetti A, Martinelli L, et al. (2015) Handbook on lead-bismuth eutectic alloy and lead properties, materials compatibility, thermal-hydraulics and technologies-2015 edition. Tech. rep., Organisation for Economic Co-Operation and Development
  * Ginter et al. (1986) Ginter G, Gasser JG, Kleim R (1986) The electrical resistivity of liquid bismuth, gallium and bismuth-gallium alloys. Philos Mag B 54(6):543–552
  * Glazier et al. (1999) Glazier JA, Segawa T, Naert A, Sano M (1999) Evidence against ‘ultrahard’ thermal turbulence at very high Rayleigh numbers. Nature 398:307–310
  * Gong et al. (2020) Gong Q, Ding W, Bonk A, Li H, Wang K, Jianu A, Weisenburger A, Bund A, Bauer T (2020) Molten iodide salt electrolyte for low-temperature low-cost sodium-based liquid metal battery. J Power Sources 475:228674
  * Herreman et al. (2019) Herreman W, Nore C, Ramos PZ, Cappanera L, Guermond JL, Weber N (2019) Numerical simulation of electrovortex flows in cylindrical fluid layers and liquid metal batteries. Phys Rev Fluids 4(11):113702
  * Herreman et al. (2020) Herreman W, Bénard S, Nore C, Personnettaz P, Cappanera L, Guermond JL (2020) Solutal buoyancy and electrovortex flow in liquid metal batteries. Phys Rev Fluids 5(7):074501
  * Hodge (2020) Hodge T (2020) Hourly electricity consumption varies throughout the day and across seasons - Today in Energy - U.S. Energy Information Administration (EIA). URL https://www.eia.gov/todayinenergy/detail.php?id=42915
  * Horstmann et al. (2018) Horstmann GM, Weber N, Weier T (2018) Coupling and stability of interfacial waves in liquid metal batteries. J Fluid Mech 845:1–35
  * Iida and Guthrie (2015a) Iida T, Guthrie RIL (2015a) The Thermophysical Properties of Metallic Liquids: Volume 1: Fundamentals. Oxford University Press, Oxford, UK
  * Iida and Guthrie (2015b) Iida T, Guthrie RIL (2015b) The Thermophysical Properties of Metallic Liquids: Volume 2: Predictive models. Oxford University Press, Oxford, UK
  * Janz et al. (1968) Janz GJ, Dampier FW, Lakshminarayanan GR, Lorenz PK, Tomkins RPT (1968) Molten salts: Volume 1. Electrical conductance, density, and viscosity data. National Bureau of Standards
  * Kelley and Sadoway (2014) Kelley DH, Sadoway DR (2014) Mixing in a liquid metal electrode. Phys Fluids 26(5):057102
  * Kelley and Weier (2018) Kelley DH, Weier T (2018) Fluid mechanics of liquid metal batteries. Appl Mech Rev 70(2):020801
  * Keogh et al. (2021) Keogh DF, Timchenko V, Reizes J, Menictas C (2021) Modelling Rayleigh-Bénard convection coupled with electro-vortex flow in liquid metal batteries. J Power Sources 501:229988
  * Kim et al. (2013a) Kim H, Boysen DA, Newhouse JM, Spatocco BL, Chung B, Burke PJ, Bradwell DJ, Jiang K, Tomaszowska AA, Wang K, Wei W, Ortiz LA, Barriga SA, Poizeau SM, Sadoway DR (2013a) Liquid metal batteries: past, present, and future. Chem Rev 113(3):2075–2099
  * Kim et al. (2013b) Kim H, Boysen DA, Ouchi T, Sadoway DR (2013b) Calcium–bismuth electrodes for large-scale energy storage (liquid metal batteries). J Power Sources 241:239–248
  * Köllner et al. (2017) Köllner T, Boeck T, Schumacher J (2017) Thermal Rayleigh-Marangoni convection in a three-layer liquid-metal-battery model. Phys Rev E 95(5):053114
  * Li et al. (2016) Li H, Yin H, Wang K, Cheng S, Jiang K, Sadoway DR (2016) Liquid metal electrodes for energy storage batteries. Adv Energy Mater 6(14):1600483
  * Lundquist (1969) Lundquist S (1969) On the hydromagnetic viscous flow generated by a diverging electric current. Arkiv fr Fysik 40(5):89––95
  * Ning et al. (2015) Ning X, Phadke S, Chung B, Yin H, Burke P, Sadoway DR (2015) Self-healing li-bi liquid metal battery for grid-scale energy storage. J Power Sources 275:370–376
  * Okada and Ozoe (1992) Okada K, Ozoe H (1992) Experimental heat transfer rates of natural convection of molten gallium suppressed under an external magnetic field in either the X, Y, or Z direction. J Heat Transfer 114(1):107–114
  * Ozer and Tepehan (1993) Ozer N, Tepehan F (1993) Structure and optical properties of electrochromic copper oxide films prepared by reactive and conventional evaporation techniques. Sol Energy Mater Sol Cells 30(1):13–26
  * Perez and Kelley (2015) Perez A, Kelley DH (2015) Ultrasound velocity measurement in a liquid metal electrode. J Vis Exp (102)
  * Personnettaz et al. (2018) Personnettaz P, Beckstein P, Landgraf S, Köllner T, Nimtz M, Weber N, Weier T (2018) Thermally driven convection in Li$\|$Bi liquid metal batteries. J Power Sources 401:362–374
  * Personnettaz et al. (2019) Personnettaz P, Landgraf S, Nimtz M, Weber N, Weier T (2019) Mass transport induced asymmetry in charge/discharge behavior of liquid metal batteries. Electrochem Commun 105:106496
  * Poizeau et al. (2012) Poizeau S, Kim H, Newhouse JM, Spatocco BL, Sadoway DR (2012) Determination and modeling of the thermodynamic properties of liquid calcium–antimony alloys. Electrochim Acta 76:8–15
  * Scheel and Schumacher (2017) Scheel JD, Schumacher J (2017) Predicting transition ranges to fully turbulent viscous boundary layers in low Prandtl number convection flows. Phys Rev Fluids 2(12):123501
  * Sele (1977) Sele T (1977) Instabilities of the metal surface in electrolytic alumina reduction cells. Metall Trans B 8(3):613–618
  * Shen and Zikanov (2016) Shen Y, Zikanov O (2016) Thermal convection in a liquid metal battery. Theoret Comput Fluid Dynamics 30(4):275–294
  * Shercliff (1970) Shercliff JA (1970) Fluid motions due to an electric current source. J Fluid Mech 40(2):241–250
  * Shishkina et al. (2016) Shishkina O, Grossmann S, Lohse D (2016) Heat and momentum transport scalings in horizontal convection. Geophys Res Lett 43(3):1219–1225
  * Takeda (1987) Takeda Y (1987) Measurement of velocity profile of mercury flow by ultrasound Doppler shift method. Nucl Technol 79(1):120–124
  * Takeshita et al. (1996) Takeshita T, Segawa T, Glazier JA, Sano M (1996) Thermal turbulence in mercury. Phys Rev Lett 76(9):1465
  * Vogt et al. (2018) Vogt T, Horn S, Grannan AM, Aurnou JM (2018) Jump rope vortex in liquid metal convection. Proc Natl Acad Sci USA 115(50):12674–12679
  * Wagner and Shishkina (2013) Wagner S, Shishkina O (2013) Aspect-ratio dependency of Rayleigh-Bénard convection in box-shaped containers. Phys Fluids 25(8):085110
  * Wang and Kelley (2021) Wang B, Kelley DH (2021) Microscale mechanisms of ultrasound velocity measurement in metal melts. Flow Meas Instrum p 102010
  * Wang et al. (2014) Wang K, Jiang K, Chung B, Ouchi T, Burke PJ, Boysen DA, Bradwell DJ, Kim H, Muecke U, Sadoway DR (2014) Lithium–antimony–lead liquid metal battery for grid-level energy storage. Nature 514(7522):348–350
  * Wang et al. (2020) Wang Q, Verzicco R, Lohse D, Shishkina O (2020) Multiple states in turbulent large-aspect-ratio thermal convection: What determines the number of convection rolls? Phys Rev Lett 125(7):074501
  * Weber et al. (2018) Weber N, Nimtz M, Personnettaz P, Salas A, Weier T (2018) Electromagnetically driven convection suitable for mass transfer enhancement in liquid metal batteries. Appl Therm Eng 143:293–301
  * Weber et al. (2020) Weber N, Nimtz M, Personnettaz P, Weier T, Sadoway D (2020) Numerical simulation of mass transfer enhancement in liquid metal batteries by means of electro-vortex flow. J Power Sources Adv 1:100004
  * Zhang et al. (2007) Zhang H, Charmchi M, Veilleux D, Faghri M (2007) Numerical and experimental investigation of melting in the presence of a magnetic field: Simulation of low-gravity environment. J Heat Transfer 129(4):568–576
  * Zikanov (2015) Zikanov O (2015) Metal pad instabilities in liquid metal batteries. Phys Rev E 92(6):063021
  * Zürner et al. (2019) Zürner T, Schindler F, Vogt T, Eckert S, Schumacher J (2019) Combined measurement of velocity and temperature in liquid metal convection. J Fluid Mech 876:1108–1128, DOI 10.1017/jfm.2019.556


## (1) Overview

## Title

An ObsPy library for event detection and seismic attribute calculation:
preparing waveforms for automated analysis

## Paper Authors

1\. Turner, Ross J;
2\. Latto, Rebecca B;
3\. Reading, Anya M

## Paper Author Roles and Affiliations

1\. School of Natural Sciences (Physics), University of Tasmania, Private Bag
37, Hobart, 7001, Australia;
2\. School of Natural Sciences (Physics), University of Tasmania, Private Bag
37, Hobart, 7001, Australia;
3\. School of Natural Sciences (Physics), University of Tasmania, Private Bag
37, Hobart, 7001, Australia; and Institute for Marine and Antarctic Studies,
University of Tasmania, Private Bag 129, Hobart, TAS 7001, Australia.

## Abstract

We have implemented an extension for the observational seismology _obspy_
software package to provide a streamlined tool tailored to the processing of
seismic signals from non-earthquake sources, in particular those from
deforming systems such as glaciers and landslides. This _seismic attributes_
library provides functionality to: (1) download and/or pre-process seismic
waveform data; (2) detect and catalogue seismic events using multi-component
signals from one or more seismometers; and (3) calculate characteristics
(‘attributes’/‘features’) of the identified events. The workflow is controlled
by three main functions that have been tested for the breadth of data types
expected from permanent and campaign-deployed seismic instrumentation. A
selected STA/LTA-type (short-term average/long-term average), or other, event
detection algorithm can be applied to the waveforms and user-defined functions
implemented to calculate any required characteristics of the detected events.
The code is written in Python 2/3 and is available on GitHub together with
detailed documentation and worked examples.

## Keywords

geophysics, seismology, data processing, reproducibility, Python

## Introduction

Seismology provides an attractive tool to investigate physical processes in
deforming systems. The seismic signals from active glaciers, for example,
could enable monitoring of mechanisms including basal sliding [6, 13, 24],
fracturing [9], melt water drainage [4], and iceberg calving [14, 15]. The
detection of seismic events from the recorded continuous seismic waveform data
is a vital first step in any analysis. Event catalogues thus constructed are
needed for local seismicity studies, comparisons between locations, or
detection of change over extended time periods. However, the automated
detection of seismic events is complicated in environmental and geotechnical
seismology by the diverse populations of signal generation mechanisms. Those
generated by active glaciers, for example, can be expected to span several
orders of magnitude in duration and amplitude [16]. Volcanoes [10], landslides
[19] and mining activity [25] similarly produce a broad range of seismic
signals. As a result, the majority of cryoseismology studies to date use
manual identification of events [6, 17, 18]. Manual techniques are not readily
scalable nor exactly reproducible: in particular, they are not a first choice
for monitoring applications nor for the data-driven detection of change.

Established event detection algorithms have largely been developed for
earthquake seismology [1, 2, 3, 7]. These algorithms, including STA/LTA
(short-term average/long-term average) [1] and template matching [3], are
applied in real-time to seismic data to detect earthquakes and produce event
catalogues [2, 7, 8, 23]. The core classes in the _obspy_ software package are
therefore designed assuming event metadata is available online alongside the
waveform data [22]; additional functions are included outside the core classes
to select events directly from waveform data. These standard algorithms are
generally only applicable to non-earthquake signals by employing an
experimental approach to parameter selection [5, 16]. The _seismic attributes_
library provides software tools to download seismic waveform data from online
repositories and detect environmental and geotechnical seismic events using a
choice of algorithms. Algorithm options include the classic, recursive and
delayed STA/LTA algorithms [22], and the newly developed multi-STA/LTA
algorithm [11]. The multi-STA/LTA can simultaneously extract both short and
long duration events of very different signal-to-noise levels and potentially
enables real-time monitoring of cryoseismic events.

Event catalogues constructed using the _seismic attributes_ library would be
expected to comprise diverse signals generated by a range of mechanisms as
noted above. The various signals typically need to be separated into related
clusters based on the characteristics of their waveforms in order to study the
events further [10, 19]. For example, Provost et al. [19] consider 71
attributes based on waveform data in their study of landslide seismicity.
These are broadly split into four categories: (1) waveform attributes (e.g.
duration, energy, kurtosis); (2) spectral and spectrogram attributes (e.g.
discrete Fourier transform); (3) network attributes (e.g. station with maximum
amplitude); and (4) polarity attributes (e.g. azimuth, inclination). The
_seismic attributes_ library includes functions to calculate a number of
standard signal properties: duration, ratio between ascending and descending
time, energy in the autocorrelation function, energy in the frequency filtered
spectrum, and the direction of wave propagation. These are provided in three
bundles of attribute functions describing the waveform, spectrum and polarity
of the signal; we do not include network attributes as these are unordered,
discrete variables (and largely application dependent). User-defined functions
can be added to derive customised characteristics within our software
architecture. The correlation between attributes can be investigated using a
plotting function to select the best subset for the subsequent application of
(e.g.) clustering to a given set of events, informing the removal of redundant
variables if appropriate. Data-driven techniques more generally, such as
machine learning algorithms, may be readily applied to the calculated
attributes to inform the current and future state of the glacier (i.e.
identify signals in the lead-up to a large event [20]).

## Implementation and architecture

The analysis of waveform data from environmental and geotechnical seismic
deployments requires multiple distinct steps. Following the principle that one
should write programs that do one thing well, and write programs that work
together [21], the _seismic waveforms_ library is correspondingly split into
three primary functions to streamline the workflow. The name and purpose of
these functions are as follows:

  * •

get_waveforms(); this function downloads waveform data from an online
repository or alternatively reads data stored locally.

  * •

get_events(); this function uses an event triggering algorithm to produce an
event catalogue based on the seismic signals in one or more components of one
or more seismic stations.

  * •

get_attributes(); this function produces a pandas DataFrame of the attributes
for each event using functions included in the library, or user-defined
attribute functions.

The library includes the numerous attribute functions that are called by
get_attributes(). Additional functions are also included to analyse the output
of the primary functions (e.g. plotting), together with several private
functions that support the primary workflow (documented in the code itself).
The implementation of the three primary functions is outlined below.

get_waveforms
Seismic datasets can build to a high data volume (several terrabytes) if high
sampling rates are required over extended periods of time, or if using data
from seismic arrays with multiple seismometers. This can result in memory
problems when processing the waveforms and also disk space limitations. The
first function in our workflow therefore acts as a gateway, to manage the
volume of data downloaded from online repositories using _obspy_ clients.
Waveform data from each seismometer is written into separate files comprising
only a single day of data for a single component (e.g. vertical, north, east).
The user can specify the write location of the files to split different time
periods across several external drives. This same function checks if the
requested data have previously been downloaded (by default) and can optionally
not check for, or not download missing data (e.g. in the case of known data
gaps). The waveform data are written as .mseed files using standard _obspy_
functions.

The get_waveforms function recombines downloaded or locally stored files into
a single _obspy_ Trace object for the requested time period, thus preparing
the downloaded files for the next steps in the workflow. The user should
examine small subsets of the time period of interest (e.g. a week at a time)
based on factors including data volume and computer hardware performance. The
seismic waveforms for the requested time period and components of a given
seismometer are output as an _obspy_ Stream object. The streams from different
seismometers are combined using the ‘+’ operator into a single object.

get_events
As noted in the introduction, several algorithms exist to trigger events from
single component waveform data. The get_events code first separates the
signals from multiple seismometers and combines their (usually) three
components into a single waveform. The user can optionally specify whether the
component waveforms are combined as the sum-of-squares of their amplitude
(i.e. Euclidean norm) to give the wave energy, or the absolute value of the
wave amplitude. Further, to ensure that taking the absolute value of the
amplitude does not affect the results (e.g. doubling frequency) we tested a
computationally intensive option to fit the (time-varying) principle component
of the direction of wave oscillation to obtain a signed amplitude; this gave
identical event detections to the absolute amplitude and so is not included in
our published version of the code due to the significantly longer computation
time.

Figure 1: Top: Diagrammatic representation of algorithm used to identify
events from a single seismometer comprising one or more channels (traces show
their Euclidean norm). The characteristic function for an STA/LTA algorithm is
used to ‘trigger’ events (shown in red) and small gaps between these events
(shorter than thr_event_join) are ignored (orange). No events are present at
other times (shown in green). Bottom: Representation of algorithm used to
identify events in the ‘reference event’ and ‘trace’ catalogues for seismic
arrays (with multiple seismometers). The example shows an indicative array
with $N=3$ seismometers and thr_coincidence_sum $\geqslant 2$ simultaneous
detections. Events identified for the single seismometers traces are shown in
red (as given in the top panel). The duration of these events is extended at
each end by half the delay in arrival time of a wave between the most distant
seismometers in the array (shown in deep red). The reference event is
identified as the times when thr_coincidence_sum $\geqslant 2$ seismometers
have a detection (shown in red or deep red); small gaps with fewer
seismometers are joined over (shown in orange). The events at each seismometer
are simply the events identified from the single seismometer traces (red but
not deep red), but with any times between events that occur during the
reference event also included (shown in grey).

We use an adapted version of the _obspy_ coincidence_trigger function to
create a reference catalogue of events based on the STA/LTA characteristic
function at each seismometer for their combined component waveforms. Our
version of the coincidence_trigger function includes adjustments in the
algorithm to better align with the outputs needed in our workflow and
improvements to computational performance that are possible for our narrower
use of this function (see Figure 1 for schematic of algorithm). This function
can use the standard STA/LTA algorithms available within obspy [22], and
optionally the recently developed multi-STA/LTA algorithm [11]. In all cases,
location-specific algorithm parameters are likely to be more successful for
environmental seismology applications. The start and stop times of the
detected event from a single seismometer can be extended to consider small
gaps in between triggers (up to length thr_event_join), creating a single,
longer-duration record. The reference event catalogue is created by finding
times when $n$ (i.e. thr_coincidence_sum) of the $N$ seismometers in the array
have temporally coincident records. The physical dimensions of the seismic
array are calculated at this step to estimate an upper bound on the delay in
arrival time of the signal at different seismometers. This delay is added to
the duration of the single seismometer records to ensure the reference events
are not artificially shortened. The reference event catalogue, which includes
a reference start time and duration for each event, is output as a _pandas_
DataFrame to provide a convenient format to write to file and interrogate.

Our adapted version of the coincidence_trigger function provides additional
utility for other researchers by outputting a secondary catalogue of trace
metadata, for each seismometer, for all identified reference events. The
triggered records at each seismometer that occur during (at least) part of the
time period covered by a given reference event are joined together (if not
already a single record). This record may extend beyond the time period of the
reference event, be contained entirely within it, or have no detection at all
in the case of weaker events. This trace (metadata) catalogue, which includes
the trigger start time and duration for all seismometers, is similarly output
as a _pandas_ DataFrame. The trace metadata can be used to extract waveform
data for identified events in the catalogue in a format that will be familiar
to seismologists, for example, as shown in Figure 2 for a low-frequency event
detected using four seismometers on the Whillans Ice Stream in Antarctica from
13:35:37 on 16 December 2010.

Figure 2: Waveforms of an event detected using four seismometers (BB01, BB03,
BB04 and BB06) comprising part of a seismic array on the Whillans Ice Stream
in Antarctica from 13:35:37 on 16 December 2010. The top, middle and bottom
panels show the vertical, north and east component of the signal respectively.
The waveforms for each seismometer start 30 seconds prior to the start time in
the trace catalogue and terminate 60 seconds after the stop time. The code to
reproduce this plot is provided as a worked example.

get_attributes
The get_attributes function is written following further principles [21], that
one should write flexible and open programs. This function extracts the
waveform data for a given seismometer for the duration of a given event listed
in the reference event, or trace (metadata) catalogue. The waveforms for each
component are stored as separate Trace objects in a single stream. This Stream
object is passed to one or multiple attribute functions to derive the
characteristics of the given event as measured by a given seismometer. The
values of the requested attributes for the array of seismometers are output as
a _pandas_ DataFrame. The correlation between the spectral attributes included
in the library are shown in Figure 3 for events detected at Ilulissat,
Greenland (DK.ILULI) on 1 January 2018 using the recursive STA/LTA algorithm.

Optionally, custom attribute functions may be added to calculate any chosen
characteristic of a waveform. Custom functions must take a Stream object
containing one or more components as an input, and output the attribute name
and value. Public functions to combine the component waveforms into the wave
energy or absolute value of the wave amplitude are included in the _seismic
attributes_ library to aid the user. The attribute functions (as many as are
required) are passed to the get_attributes function as optional parameters,
with return values stored as above.

Figure 3: Corner (or pair) plot illustrating the correlation between spectral
attributes for the events detected at Ilulissat, Greenland on 1 January 2018
using the recursive STA/LTA algorithm. The attribute names are defined to
match those used by Provost et al. [19]. In this example, attribute 26
considers a frequency range close to the sampling rate of the recorded signal,
therefore, it has null value for all events. The code to reproduce this plot
is provided as a worked example in the detailed documentation provided with
the software.

## Quality control

We have tested the seismic attributes library to find any software bugs and to
ensure outputs are reproducible by other researchers running the code or using
different platforms [12]. We first created test cases based on glacier seismic
signals from a field campaign seismic array on the Whillans Ice Stream in West
Antarctica [24]. Plausible uses of the code were tested by selecting some or
all of the seismometers in the array, some or all of the three components of
the signal recorded at each seismometer (i.e. vertical, north and east), and
considering different durations of requested time (with and without small data
gaps). The optional function inputs were also tested in this manner. In
application, e.g. [11], the event catalogue produced using the _seismic
attributes_ workflow compared favourably to events that were visually
classified using the same set of seismometers [18]. Finally, the attribute
functions included in the library were tested on mock waveform data to ensure
they produce the theoretically expected results.

Testing was primarily carried out using a Python 3.8 installation on a macOS
10/11 system (Unix). The key features of the _seismic attributes_ library were
tested on other platforms for compatability, including both Python 2.7 and
3.7/8 for the most up to date version of _obspy_ , _numpy_ , _pandas_ and
_seaborn_ available for that installation. The code was further tested on both
macOS 10/11 and Windows operating systems, especially to verify functionality
of the different file systems. The reference event and trace (metadata)
catalogues were compared between these platforms to assess the consistency of
code functionality.

## (2) Availability

## Operating system

macOS 10/11 (including Unix and Linux) and Windows.

## Programming language

Python 2 or 3.

## Additional system requirements

Memory and disk space will limit the volume of data that can be processed in a
contiguous chunk; the code will work on any system that can support _obspy_.
Internet access is required to download new waveform data from online
repositories.

## Dependencies

The minimal dependency for use is Python 2/3 with _obspy_ , _numpy_ , _pandas_
and _seaborn_ packages installed.

## Software location:

Archive

Name:

GitHub

Persistent identifier:

https://github.com/rossjturner/seismic_attributes

Licence:

GPL version 3

Publisher:

Ross James Turner

Date published:

14/10/2020

Code repository

Name:

GitHub

Persistent identifier:

https://github.com/rossjturner/seismic_attributes

Licence:

GPL version 3

Date published:

14/10/2020

## Language

Python 2 or 3

## (3) Reuse potential

The main purpose behind the _seismic attributes_ library is to provide a
streamlined workflow for researchers to detect and categorise seismic events
from environmental and geotechnical sources in a rigorous and reproducible
manner. The code will find extensive use in areas of seismology where a
diverse population of seismic signals are present, including glaciers,
volcanoes, landslides and mine sites. The get_events function is expected to
be especially useful, in particular in glacier seismology, to provide a
consistent method for creating event catalogues for ongoing machine learning
and conventional seismological analysis. We anticipate that our get_attributes
function may also be useful in applications outside seismology to aid in the
construction of a catalogue of waveform attributes for use in machine
learning. The code would need only minor modifications to handle time series
data, outside of seismology, stored in different formats. It could conceivably
find applications other areas of science (e.g. astrophysics), economics and
many other potential applications. Limited support may be provided by
contacting the corresponding author.

## Funding statement

This research was supported under Australian Research Council’s (ARC) Special
Research Initiative for Antarctic Gateway Partnership (project ID
SR140300001), and ARC Discovery Projects DP190100418 and DP2101000834.

## Competing interests

The authors declare that they have no competing interests.

## References

  * [1] R. Allen.  Automatic phase pickers: Their present use and future prospects.  Bulletin of the Seismological Society of America, 72(6B):S225–S242, 1982.
  * [2] R. V. Allen.  Automatic earthquake recognition and timing from single traces.  Bulletin of the Seismological Society of America, 68(5):1521–1532, 10 1978.
  * [3] N. A. Anstey.  The sectional auto-correlogram and the sectional retro-correlogram* part i: The sectional auto-correlogram.  Geophysical Prospecting, 14(4):389–426, 1966.
  * [4] N. Aso, V. C. Tsai, C. Schoof, G. E. Flowers, A. Whiteford, and C. Rada.  Seismologically Observed Spatiotemporal Drainage Activity at Moulins.  Journal of Geophysical Research: Solid Earth, 122(11):9095–9108, 2017.
  * [5] R. Aster and J. Winberry.  Glacial seismology.  Reports on Progress in Physics, 80(12):126801, 2017.
  * [6] C. G. Barcheck, S. Tulaczyk, S. Y. Schwartz, J. I. Walter, and J. P. Winberry.  Implications of basal micro-earthquakes and tremor for ice stream mechanics: Stick-slip basal sliding and till erosion.  Earth and Planetary Science Letters, 486:54 – 60, 2018.
  * [7] P. S. Earle and P. M. Shearer.  Characterization of global seismograms using an automatic-picking algorithm.  Bulletin of the Seismological Society of America, 84(2):366–376, 04 1994.
  * [8] D. Houliston, G. Waugh, and J. Laughlin.  Automatic real-time event detection for seismic networks.  Computers & Geosciences, 10(4):431 – 436, 1984.
  * [9] J. Kavanaugh, R. Schultz, L. Andriashek, M. Baan, H. Ghofrani, G. Atkinson, and D. Utting.  A New Year’s Day Icebreaker: Icequakes on Lakes in Alberta, Canada.  Canadian Journal of Earth Sciences, 56:183–200, 01 2019.
  * [10] A. Köhler, M. Ohrnberger, and F. Scherbaum.  Unsupervised pattern recognition in continuous seismic wavefield records using Self-Organizing Maps.  Geophysical Journal International, 182(3):1619–1630, 2010.
  * [11] R. B. Latto, R. J. Turner, A. M. Reading, S. Cook, and P. Winberry.  Event Detection in Cryoseismology.  In presented at 2020 Fall Meeting, AGU, 1-17 Dec, pages [C032–07], 2020.
  * [12] R. J. LeVeque, I. M. Mitchell, and V. Stodden.  Reproducible research for scientific computing: Tools and strategies for changing the culture.  Computing in Science Engineering, 14(4):13–17, 2012.
  * [13] B. P. Lipovsky, C. R. Meyer, L. K. Zoet, C. McCarthy, D. D. Hansen, A. W. Rempel, and F. Gimbert.  Glacier sliding, seismicity and sediment entrainment.  Annals of Glaciology, 60(79):182–192, 2019.
  * [14] M. Nettles and G. Ekström.  Glacial Earthquakes in Greenland and Antarctica.  Annual Review of Earth and Planetary Sciences, 38(1):467–491, 2010\.
  * [15] K. G. Olsen and M. Nettles.  Constraints on Terminus Dynamics at Greenland Glaciers From Small Glacial Earthquakes.  Journal of Geophysical Research: Earth Surface, 124(7):1899–1918, 2019.
  * [16] E. A. Podolskiy and F. Walter.  Cryoseismology.  Reviews of Geophysics, 54(4):708–758, 2016.
  * [17] J. Pomeroy, A. Brisbourne, J. Evans, and D. Graham.  The search for seismic signatures of movement at the glacier bed in a polythermal valley glacier.  Annals of Glaciology, 54(64):149–156, 2013.
  * [18] M. J. Pratt, J. P. Winberry, D. A. Wiens, S. Anandakrishnan, and R. B. Alley.  Seismic and geodetic evidence for grounding-line control of Whillans Ice Stream stick-slip events.  Journal of Geophysical Research: Earth Surface, 119(2):333–348, 2014.
  * [19] F. Provost, C. Hibert, and J.-P. Malet.  Automatic classification of endogenous landslide seismicity using the Random Forest supervised classifier.  Geophysical Research Letters, 44(1):113–120, 2017.
  * [20] B. Rouet-Leduc, C. Hulbert, N. Lubbers, K. Barros, C. J. Humphreys, and P. A. Johnson.  Machine learning predicts laboratory earthquakes.  Geophysical Research Letters, 44(18):9276–9282, 2017.
  * [21] P. H. Salus.  A quarter century of UNIX / Peter H. Salus .  Addison-Wesley Pub. Co Reading, Mass, 1994.
  * [22] O. D. Team.  obspy.core - core classes of obspy¶.
  * [23] Y. Vaezi and M. Van der Baan.  Comparison of the STA/LTA and power spectral density methods for microseismic event detection.  Geophysical Supplements to the Monthly Notices of the Royal Astronomical Society, 203(3):1896–1908, 2015.
  * [24] J. P. Winberry, S. Anandakrishnan, D. A. Wiens, and R. B. Alley.  Nucleation and seismic tremor associated with the glacial earthquakes of Whillans Ice Stream, Antarctica.  Geophysical Research Letters, 40(2):312–315, 2013.
  * [25] Z. Zhou, R. Cheng, X. Cai, D. Ma, and C. Jiang.  Discrimination of Rock Fracture and Blast Events Based on Signal Complexity and Machine Learning.  Shock and Vibration, 2018:1–10, 02 2018.

Copyright Notice
Authors who publish with this journal agree to the following terms:

Authors retain copyright and grant the journal right of first publication with
the work simultaneously licensed under a Creative Commons Attribution License
that allows others to share the work with an acknowledgement of the work’s
authorship and initial publication in this journal.

Authors are able to enter into separate, additional contractual arrangements
for the non-exclusive distribution of the journal’s published version of the
work (e.g., post it to an institutional repository or publish it in a book),
with an acknowledgement of its initial publication in this journal.

By submitting this paper you agree to the terms of this Copyright Notice,
which will apply to this submission if and when it is published by this
journal.