# Detecting Anomalies using Overlapping Electrical Measurements in Smart Power
Grids

Sina Sontowski1, Nigel Lawrence2, Deepjyoti Deka3, and Maanak Gupta4
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>14Dept. of Computer Science, Tennessee Technological
University, Cookeville, Tennessee, USA
23Los Alamos National Laboratory, Los Alamos, New Mexico, USA

###### Abstract

As cyber-attacks against critical infrastructure become more frequent, it is
increasingly important to be able to rapidly identify and respond to these
threats. This work investigates two independent systems with overlapping
electrical measurements with the goal to more rapidly identify anomalies. The
independent systems include HIST, a SCADA historian, and ION, an automatic
meter reading system (AMR). While prior research has explored the benefits of
fusing measurements, the possibility of overlapping measurements from an
existing electrical system has not been investigated. To that end, we explore
the potential benefits of combining overlapping measurements both to improve
the speed/accuracy of anomaly detection and to provide additional validation
of the collected measurements. In this paper, we show that merging overlapping
measurements provide a more holistic picture of the observed systems. By
applying Dynamic Time Warping more anomalies were found – specifically, an
average of 349 times more anomalies, when considering anomalies from both
overlapping measurements. When merging the overlapping measurements, a percent
change of anomalies of up to 785% can be achieved compared to a non-merge of
the data as reflected by experimental results.

###### Index Terms:

Unsupervised Anomaly Detection, Electrical Measurements, Dynamic Time Warping

##  I Introduction

In recent years, cyber attacks in critical infrastructures have become more
frequent [1, 2, 3, 4] and increased in complexity and sophistication. Smart
power grids are vulnerable to these attacks in large part due to an increase
in cyber-physical connectivity in Supervisory Control and Data Acquisition
(SCADA) systems that are used to monitor and control parts of the power grid
[5]. While newer digital devices, even at the edge of the grid such as micro-
phasor measurement units ($\mu$-PMUs) [6, 7], can increase connectivity,
control, and estimation tasks used in power grid operation, they also increase
the risk from potential cyber threats [8]. An unobserved or hidden cyber
attack on the power grid can de-energize power system components and aggregate
operating conditions by causing overloading and instability [9]. The
2015111https://www.wired.com/2016/03/inside-cunning-unprecedented-hack-
ukraines-power-grid/ Ukraine power grid hack exemplifies the threat that cyber
attacks can pose to the power grid. Many consumers in Ukraine temporarily lost
power due to a cyber attack exploiting Windows vulnerability CVE-2014-4114. A
phishing email was opened that installed BlackEnergy malware on the system.
This was not a zero-day vulnerability and the operator should have been aware
of the vulnerability in the system and taken action beforehand [8]. While
anomaly detection solutions would have not prevented this attack, they would
have alerted the operator about unusual activity in the system, potentially
allowing them to respond more rapidly. Accordingly, anomaly detection
identifies potential issues in the system and can therefore enable operators
to respond to malicious activity as it happens [10]. Thus, it allows for an
early detection of cyber-intrusions, and can be understood as an early warning
mechanism [9]. In general, anomaly detection refers to the problem of finding
patterns in data that do not conform to expected behavior, and it is being
widely used in different domains [11, 12, 13, 14, 15, 16].

The objective of this paper is to develop improved techniques for anomaly
detection in smart power grids using electrical measurements by which
operators can be alerted early of any abnormal activity, and action can be
taken to ameliorate from that state. Specifically, we present an approach to
improve the speed and accuracy of anomaly detection by using overlapping
electrical measurements from two independent systems taken at different points
throughout the United States Department of Energy Los Alamos National
Laboratory’s power grid. Our approach consists of three parts as illustrated
in Figure 1. First, we determine the overlapping measurements by using Dynamic
Time Warping (DTW). Then unsupervised anomaly detection algorithms including
autoregression, level shift, and rolling average are applied to the
overlapping measurements. This approach differs from existing work (as
discussed in related work section) as the overlapping electrical measurements
used for anomaly detection are from a real power grid. To the best of our
knowledge, no existing work has applied DTW as a way to improve anomaly
detection. In the final stage, we analyze the anomalies that were discovered
through merging the overlapping measurements, and we compare them to their
individual anomalies.

Figure 1: Overview of Proposed Approach

The main contributions of this paper are as follows:

  * •

We elaborate approaches used for dealing with overlapping electrical
measurements, in addition to verifying them.

  * •

We compare results of different unsupervised anomaly detection algorithms
applied on real-life electrical measurements.

  * •

We conjecture that merging overlapping measurements for anomaly detection
gives a more holistic view and finds more anomalies.

The rest of the paper has been organized as follows. In Section II, the paper
explores related works and their significance. In Section III, background
information and preliminaries are provided for a better understanding of our
implementation described in Section IV and V. A conclusion is provided and
possible future work approaches are mentioned in Section VI and Section VII,
respectively.

##  II Related Work

A neuro-inspired architecture called Hierarchical Temporal Memory (HTM) was
developed to perform unsupervised anomaly detection [17]. HTM learns sparse
distributed temporal representation of sequential data which is shown to
achieve competitive scores for real-time anomaly detection compared to state-
of-the-art approaches. Anomaly detection is also performed on unlabeled
electrical measurements ($\mu$-PMU) from the smart grid. HTM was able to
capture spatial and temporal anomalies in the $\mu$-PMU data from the power
grid. While electrical measurements ($\mu$-PMU data) were used for this
approach, the anomaly detection is real-time, involves hierarchical temporal
memory, but does not include overlapping electrical measurements as proposed
in our work.

Researchers have used dynamic time warping (DTW) related to the power grid;
however, Elafoudi et al. [18], [19] apply it to smart meter readings for
lowering the complexity of power disaggregation. Ausmus et al. [19] apply DTW
to electric utility data to cluster the electric utility net data based on the
distance measure to improve operation for the next day. None of these papers
apply DTW with the goal to improve anomaly detection. Another work by Diab et
al. uses DTW for anomaly detection [20]. Instead of using electrical
measurements, it used network traffic data to detect network anomalies. The
network traffic is decomposed into control and data planes. Based on the DTW
distance between these two, the network activities are classified as either
benign or anomalous. The goal of Zheng’s et al.’s work [21] is to detect road
anomalies based on acceleration data. The data windows of various length are
compared with DTW. Anomalies are automatically identified by machine learning
algorithms, and the types are distinguished with DTW. Their method improved
the time consumption of a random forest filter. Although, this work uses DTW
and has a similar goal of improving the speed of anomaly detection, it uses
road acceleration data and not electrical measurements.

Autoregressive, rolling average, and level shift processes are common anomaly
detection approaches. In Zhou’s and Li’s work [22] multilevel autoregression
is applied on network traffic for anomaly detection. This method proved to be
successful by correctly detecting more than 95% of network anomalies. Often,
autoregressive processes and rolling average are combined as illustrated in
[23] and [24]. Schmidt et al. [23] apply the ARIMA algorithm, which is a
combination of autoregression and rolling average, on real-time cloud
monitoring data. Yaacob et al. [24] use ARIMA for network anomaly detection.
Solomentsev et al. [25] apply level shift detection during radar operation,
while Geng and Lai [26] apply it on sensor networks.

In summary, the related work does not address the challenge of improving
anomaly detection based on overlapping electrical measurements, the core idea
of our proposed work.

##  III Preliminaries and Background

###  III-A Dataset

The data set used for this paper originates from Los Alamos National
Laboratory’s power grid from two independent overlapping electrical
measurement systems: ION and HIST. HIST is a SCADA historian and records all
of the data measured by the SCADA system. ION is an automatic meter reading
(AMR) system that pulls electrical metering information from different meters
and records them. As a result, the SCADA system does not measure all of the
ION points, while the ION system does not measure all of the SCADA points.
Therefore, they are two separate systems that share measurements for some
meters, but not all of them. The data set was collected in December 2020 and
contains 744-756 ION data points and 534,686 HIST data points for each time
series. HIST data is collected at a higher frequency (approximately every 5
seconds) while ION is collected at a much lower frequency (approximately every
hour). As illustrated in Figure 2, a subset of measurement points are common
to both systems and are therefore said to be overlapping.

Figure 2: Overview of Dataset

###  III-B Dynamic Time Warping

The DTW algorithm is utilized to identify the overlapping measurements. It
finds the similarity between two time series by calculating the distance
between them. The lower the distance, the more similar the two time series are
to each other. DTW involves a non-linear optimal alignment which ensures that
similar time series match each other even if they are out of phase on the
x-axis or require compression or expansion. This is different from Euclidean
distance, where the i-th point in one time series is aligned with the i-th
point in the other (one-to-one mapping). Instead, the DTW algorithm can make
one-to-many mappings between points [27]. To find the best alignment, the DTW
applies a warping function which minimizes the distance between the two time
series. In our case, DTW is implemented to account for the difference in
frequency between the two systems as well as any potential minor discrepancies
in clock synchronization. DTW has applications in many domains such as
robotics, data mining, and manufacturing [28].

(a) First Ranked Match

(b) Second Ranked Match

(c) Third Ranked Match

(d) Fourth Ranked Match

Figure 3: Top Four Results Step-size.

Each ranked match (a), (b), (c), and (d) includes four graphs. The top-left
graph includes both measurements drawn together, while the top-right graph
draws their step-size sample. The bottom two graphs for each match show the
actual measurements drawn independent from each other.

###  III-C Unsupervised Anomaly Detection

Unsupervised anomaly detection is executed on the overlapping electrical
measurements. Anomaly detection, also called outlier detection, is applied in
various fields to find irregularities within data. Anomaly detection
approaches slightly differ depending on the data set, for example, a time
series will have modified algorithms compared to a multidimensional data facet
[29]. The anomaly detection algorithms implemented in this paper to compare
anomalies between the overlapping time series include autoregression (AR),
level shift (LS), and rolling average (RA). All of these algorithms are
statistical unsupervised anomaly detection approaches, which means they are
used to find outliers within unlabeled time series data sets.

####  III-C1 Autoregression (AR)

AR looks at autoregressive behavior changes to detect anomalies. AR models are
used for time series analysis. Basically, linear regression is applied on the
current data series and predictions for future values are created based on
past values. Similar to linear regression, the outcome variable (Y) at a
certain point in time is related to the predictor variable (X). However, for
AR models past values of (Y) are factored into predicting (Y). Therefore,
behavior is modeled based on past data [30, 31].

####  III-C2 Level Shift (LS)

Level shift detects level shift outliers, which are often represented by a
step function. These outliers occur when a sudden change in the mean level
move an outlier and its following data onto a new level. This can be seasonal
or not. Level shift is often used when the data has a lot of outliers, because
it is not as sensitive to spikes in data. It works by taking two sliding
windows and comparing their median values to detect a shift of values [30].

####  III-C3 Rolling Average (RA)

Rolling average or Moving Average is similar to the AR model. However, instead
of using past values of (Y) to predict (Y), past forecast errors are used to
predict (Y) [30].

##  IV Dynamic Time Warping

###  IV-A Implementation

To determine similarity among the two different systems - HIST and ION, the
DTW algorithm is used. Due to the overhead of the traditional DTW, we instead
chose to use the fastDTW [28] implementation from the DTAIDistance library.
Additionally, several sampling strategies were implemented to reduce run-time
as opposed to using each time series in its entirety.

We considered three approaches to obtain a viable sample of the data:

  1. 1.

Step size, which includes collecting data points at different steps, and
examples can be seen in Figure 3d, illustrated in the top-right graph for each
overall result in (a)-(d),

  2. 2.

Certain amount of points, for example, collecting the first 100 points, and

  3. 3.

Range of dates, for example, a day’s worth of data.

####  IV-A1 Step Size

Run Number | HIST step size | ION step size | Run-time (s)
---|---|---|---
1 | 100 | 2 | 2318.6
2 | 1000 | 1 | 616.7
3 | 1000 | 2 | 342.8
4 | 2000 | 2 | 222.6
5 | 3000 | 4 | 108.6
6 | 5000 | 7 | 65.7
TABLE I: Step Size Running Times.

For step-size, we collected a total of 6 runs, which are depicted in Table I.
DTW can handle sequences of different lengths to the same ability as equal
length sequences. In fact, there is no significant difference in accuracy
between equal-length sequences and variable-length sequences [32]. Therefore,
the step-size was reduced as much as possible to include the most points
without regard to the measured frequencies. However, the smaller the steps,
the more data points, and the longer the DTW algorithm takes to compare all of
the time series.

####  IV-A2 Amount of Points

For certain amount of points, we completed two runs total, where we took the
first 100 and 200 points as illustrated in Table II. Similar to the step runs,
as the amount of points increase, the run-time increases, too.

Run Number | Point Amount | Run-time (s)
---|---|---
7 | 100 | 24.3
8 | 200 | 38.6
TABLE II: Point Amount Running Times. Run Number |  Date Range |  HIST step size |  ION step size |  Run-time (s)
---|---|---|---|---
9 |  3 days |  50 |  1 |  191.8
TABLE III: Date Range Running Times.

(a) step size 100 HIST, 2 ION

(b) three day range and step size

(c) first 200 points

Figure 4: Top ten results of the different runs.

####  IV-A3 Date Range

For date range, one run was completed including data over a period of three
days, as can be seen in Table III, instead of including four weeks of
measurements. However, this method took a long time compared to the others,
because each individual date had to be filtered beforehand by checking if it
is within the three day range. The additional computation time is not included
in the run-time listed in the table. However, by reducing the date range to
three days, we were able to reduce the step size and include more data points.

###  IV-B Experimental Results

For each run, the DTW calculates the distance for each ION and HIST pair and
sorts them based on distance. The lower the distance, the more similar the
data. The DTW distance results for the first 200 points stayed more stable in
distance compared to the step size and the three day range results due to
deviations in the data typically appearing later in the time series. The first
200 points and step size results have several of the same matches, while the
top two results remain the same. The three day results are not as similar as
the other two sample approaches, but they include some of the same top ten
results. The results of the runs are shown in Figure 4. Finally, we decided to
focus on the step size results as an accurate similarity ranking due to its
inclusion of the whole range of data and manageable run-time.

To validate the results of the DTW for step size, we look at the statistics
for the top four matches and lowest-ranked four matches ranked by the DTW,
seen in Table IV and Table V.

The count for both the lowest and highest-ranked matches is similar with HIST
possessing about 640 times more points compared to the ION measurements which
is expected because HIST is collected at a higher frequency. While the
standard deviation for the ION measurements for the top four matches is zero
because they form a straight line, the standard deviation for the ION
measurements for the lowest-ranked matches is higher. This implies that the
ION measurements for the lowest-ranked matches are more spread out.

Merging the ION and HIST data points for all the matches show that the four
highest-ranked matches are more similar to each other compared to the lowest-
ranked four matches. The standard deviation for the merged data for the top
four matches stays constant, while the standard deviation for the lowest-
ranked merged matches increases substantially compared to the individual
standard deviations. This is because the HIST and ION measurements for the
lowest-ranked four matches are less similar to each other and therefore the
values will be more spread out when the data is merged which also leads to an
increased data range. Based on these differences, we were able to verify that
the DTW ranked the measurements correctly.

measurement name |  count |  mean |  std |  min |  max
---|---|---|---|---|---
ION-4-3472 |  744 |  0 |  0 |  0 |  0
HIST-40-S |  481302 |  0.0098 |  8.8 |  -2048 |  1792
Merged |  481892 |  0.0098 |  8.8 |  -2048 |  1792
ION-5-139 |  739 |  0.013 |  0 |  0.013 |  0.013
HIST-40-S |  481302 |  0.0098 |  8.8 |  -2048 |  1792
Merged |  481896 |  0.0098 |  8.8 |  -2048 |  1792
ION-4-3472 |  744 |  0 |  0 |  0 |  0
HIST-44-S |  481301 |  0.014 |  9.4 |  -2048 |  1748
Merged |  481891 |  0.014 |  9.4 |  -2048 |  1748
ION-5-139 |  739 |  0.013 |  0 |  0.013 |  0.013
HIST-44-S |  481301 |  0.014 |  9.4 |  -2048 |  1748
Merged |  481895 |  0.014 |  9.4 |  -2048 |  1748
TABLE IV: Statistics Top Four Ranked Matches. measurement name |  count |  mean |  std |  min |  max
---|---|---|---|---|---
ION-4-198 |  756 |  $6.5x10^{7}$ |  188 |  $6.5x10^{7}$ |  $6.5x10^{7}$
HIST-23-S |  891237 |  -91.9 |  6.6 |  -99 |  100
Merged |  891730 |  54631 |  1878654 |  -99 |  $6.5x10^{7}$
ION-4-198 |  756 |  $6.5x10^{7}$ |  188 |  $6.5x10^{7}$ |  $6.5x10^{7}$
HIST-17-S |  891237 |  -565 |  46 |  -782 |  -465
Merged |  891730 |  54159 |  1878668 |  -782 |  $6.5x10^{7}$
ION-4-198 |  756 |  $6.5x10^{7}$ |  188 |  $6.5x10^{7}$ |  $6.5x10^{7}$
HIST-16-S |  891237 |  -628 |  47 |  -785 |  -538
Merged |  891730 |  54096 |  1878670 |  -785 |  $6.5x10^{7}$
ION-4-198 |  756 |  $6.5x10^{7}$ |  188 |  $6.5x10^{7}$ |  $6.5x10^{7}$
HIST-24-S |  891237 |  -2048 |  0 |  -2048 |  -2048
Merged |  891730 |  52677 |  1878711 |  -2048 |  $6.5x10^{7}$
TABLE V: Statistics Four Lowest-Ranked Matches.

###  IV-C Discussion

The four most similar results from the DTW 100 and 2 step size are depicted in
Figure 3d. Most of the points for these graphs are concentrated on y = 0, with
some HIST outlier points deviating from the axis. The ION measurements on the
bottom left of each graph are a straight line while the HIST points on the
bottom right have more spikes. The DTW algorithm matches them, because the
HIST measurement points are concentrated on the same line as the ION
counterpart, but with a few points that fall outside that line. The difference
in pattern between the time series when looking at the bottom two graphs
within each figure could mean one of two things, they could be different
measurements; or they are the same but the ION measurements are not measured
at a high enough frequency to record the spikes that appear in their HIST
counterpart. The representative sample for the HIST measurements also does not
include all the spikes that the actual measurements have because of the
increased step size. This was the case with all the different methods of
choosing a representative sample. Therefore, the HIST sample could have
negatively affected the results of the DTW algorithm, because it does not
include the entirety of spikes from the original data. The sample data is
drawn in the top right corner for each graph in Figure 3d and clearly shows
less spikes compared to the actual HIST measurement. After we obtained these
results, we used the four most similar matches from the DTW 100 HIST, 2 ION
steps for anomaly detection as illustrated in Figure 3d and Figure 4(a) (up to
and including ION-5-139 and HIST-44-S).

measurement name |  Rolling Average |  Autoregression |  Level Shift
---|---|---|---
ION-4-3472 |  0 |  0 |  0
HIST-40-S |  94 |  276 |  0
Merged |  94 |  276 |  0
ION-5-139 |  0 |  0 |  0
HIST-40-S |  94 |  276 |  0
Merged |  832 |  274 |  4
ION-4-3472 |  0 |  0 |  0
HIST-44-S |  221 |  804 |  6
Merged |  221 |  804 |  4
ION-5-139 |  0 |  0 |  0
HIST-44-S |  221 |  804 |  6
Merged |  959 |  2269 |  16
TABLE VI: Number of Anomalies Top Four Matches.

##  V Anomaly Detection

###  V-A Implementation

We took the top four matches of the DTW results and applied autoregression
(AR), level shift (LS), and rolling average (RA) on them. Each anomaly
detection algorithm is used separately on each matched pair. The data of the
matched pairs are then combined into one time series to provide the same
context and to possibly reduce the run-time of anomaly detection by running
the algorithms on one time series compared to two. After combining the data,
anomaly detection is applied on the time series as seen in Table VI.

Figure 5: Autoregression HIST 40 S Figure 6: Autoregression ION 4 3472

###  V-B Results

For the top four matches, up to 94 to 804 times more anomalies were found when
considering anomalies detected by both overlapping measurements compared to
looking at a single measurement. These numbers depend on the anomaly detection
algorithm applied on the data. Most anomalies are missed when only the ION
measurements are considered instead of both HIST and ION measurements. When
merging the data of the top four matches, a percent change of up to 182 to
785% in anomalies can be detected compared to not merging the data. Merging
the data therefore provides a more holistic view, because more anomalies are
detected. However, identifying critical anomalies and false positives from the
detected anomalies is a challenge because the data set is unlabeled. When
looking at the number of anomalies, all ION measurements have no anomalies
because they are a straight line or are close to forming a straight line. For
the top four anomaly matches, autoregression discovers the most anomalies
after rolling average, and then level shift. Level shift is less susceptible
to noisy data and therefore detects the least number of anomalies. When the
ION and HIST measurements were combined for the top four matches, the amount
of anomalies either stayed the same or increased except for ION-4-3472 and
HIST-44-S during level shift anomaly detection. This supports the point that
merging the measurements leads to a more holistic view. The differences in
anomalies among DTW matches are displayed in Figure 5 and 6. Figure 5 has
anomalies mostly along its spikes. In comparison, Figure 6 is a straight line
with no deviations and therefore has no anomalies.

###  V-C Discussion

Since autoregression works by comparing its previous value to the current one,
the spikes should appear one at a time which also highlights that the step
size sample approach might have easily missed the spikes because they are
relatively rare. However, the mismatch in anomalies highlights the need to
look at both HIST and ION measurements instead of only looking at one of them.
The spikes, which were identified as anomalies, would have been missed by just
looking at the ION measurements. Additionally, using both HIST and ION for
anomaly detection increases the amount of data points and decreases the time
between measurements. Therefore, using both measurements can increase the
accuracy of anomaly detection, because they provide a more holistic view of
the same measurement. The challenge with evaluating the accuracy of the
different anomaly detection algorithms is that the data is unlabeled and it is
unknown if the anomalies that are identified by the anomaly detection
algorithms are cause for concern (i.e., ”critical”). Critical anomalies can
help with detecting cyber attacks or unusual activity in the system, while
naturally occurring anomalies are common to the system. Further research needs
to be done on how to evaluate the anomaly detection models and how to identify
critical anomalies and false positives.

##  VI Conclusion

Overall, this paper examines approaches used with overlapping electrical
measurements from a power grid with the goal of improving anomaly detection.
While there were challenges with the evaluation of the unsupervised anomaly
detection, it was demonstrated that merging the overlapping measurements
extends its time frame and provides a more holistic view of the data. In
general, this research identified ways to incorporate overlapping measurements
for unsupervised anomaly detection. The DTW algorithm matched the most similar
measurements, and we were able to validate the matching results by looking at
the graphs and statistics of the highest-ranked and lowest-ranked matches. By
using DTW, we were able to identify more anomalies. The number of anomalies
identified by the anomaly detection algorithms differ among the overlapping
measurement matches which provide additional challenges that need to be
addressed. The difference in frequency between the two independent system
introduces challenges that are most likely impacting the DTW results and the
difference in anomalies. Overall, this work is important, because anomaly
detection mechanism can serve as an early warning system when trying to detect
cyber attacks.

##  VII Future Work

Future work includes analyzing other statistical methods that assess when the
differences in anomalies are relevant enough that the measurements are not the
same. Especially, investigating the mismatch between anomalies for the same
measurements, which could indicate certain types of data tampering. In
addition, exploring different ways to filter the data before running DTW to
reduce run-times and improve the scalability of the approach. Using the
pearson correlation coefficient would be a possible option to find similar
time series within each system and then filter based on the results.

Another possible area for future work includes determining which anomaly
detection algorithm is the most accurate one. A possible approach would be to
simulate an artificial anomaly and see which algorithm detects it. If the zero
measurements are excluded from the overall data, a Denial-of-service (DoS)
attack can be inserted. A DoS attack is one of the most common threats to
synchrophasor systems [33] and overwhelms a PMU with bogus frames so that
legitimate frames are lost, delayed, denied, or dropped [34]. The consequences
include that the real-time measurements would be delayed or dropped. For this
case, we considered inserting zero measurements over a time period of about
six to eight seconds to indicate a DoS attack. The zero measurements indicate
dropped packages during transmission and would show up when looking at the
database of sensor measurements. Ideally, both overlapping measurements would
more accurately identify the anomalies and allow us to compare the accuracy of
the different anomaly detection algorithms.

The same could be done for a data integrity attack, which includes tampering
with the signal measurements units of devices through interference, or
changing calibration, forging data, or even GPS spoofing [34]. We considered
adding artificial anomalies, for example, adding noisy data points from the
Gaussian distribution. Because the HIST measurements are already noisy due to
the many spikes, adding the artificial points would make it difficult to
distinguish between the actual attack and the occurring spikes. The cause for
the spikes in the data needs to be defined first to be able to include a data
integrity attack. If the spikes are naturally occurring noisy data points due
to poor calibration, these spikes could be removed and then a data integrity
attack can be inserted.

Evaluating our approach using overlapping measurements on a labeled data set
to confirm the effectiveness is another direction of future work. One possible
approach is generating a synthetic data set and labeling it to change it to a
supervised anomaly detection problem which has more known methods. However,
the shape of the data needs to be explained, in addition to defining the
anomalies in more detail to make sense with the shape.

## Acknowledgment

The authors would like to thank the Information Science and Technology
Institute (ISTI) at Los Alamos National Laboratory. This work is supported in
part by the National Science Foundation under Grant 2043324 and 2025682.

## References

  * [1] S. Sontowski, M. Gupta, S. S. L. Chukkapalli, M. Abdelsalam, S. Mittal, A. Joshi, and R. Sandhu, “Cyber attacks on smart farming infrastructure,” in _2020 IEEE 6th International Conference on Collaboration and Internet Computing (CIC)_. IEEE, 2020, pp. 135–143.
  * [2] A. Hahn and M. Govindarasu, “Cyber attack exposure evaluation framework for the smart grid,” _IEEE Transactions on Smart Grid_ , vol. 2, no. 4, pp. 835–843, 2011.
  * [3] M. Gupta and R. Sandhu, “Authorization framework for secure cloud assisted connected cars and vehicular internet of things,” in _Proceedings of the 23nd ACM on symposium on access control models and technologies_ , 2018, pp. 193–204.
  * [4] T. D. Le, A. Anwar, S. W. Loke, R. Beuran, and Y. Tan, “Gridattacksim: A cyber attack simulation framework for smart grids,” _Electronics_ , vol. 9, no. 8, p. 1218, 2020.
  * [5] B. Zhu, A. Joseph, and S. Sastry, “A taxonomy of cyber attacks on scada systems,” in _2011 International conference on internet of things and 4th international conference on cyber, physical and social computing_. IEEE, 2011, pp. 380–388.
  * [6] Y. Zhou, R. Arghandeh, I. Konstantakopoulos, S. Abdullah, A. von Meier, and C. J. Spanos, “Abnormal event detection with high resolution micro-pmu data,” in _2016 Power Systems Computation Conference (PSCC)_. IEEE, 2016, pp. 1–7.
  * [7] D. Deka, R. Baldick, and S. Vishwanath, “Optimal data attacks on power grids: Leveraging detection & measurement jamming,” in _2015 IEEE International Conference on Smart Grid Communications (SmartGridComm)_. IEEE, 2015, pp. 392–397.
  * [8] V. Venkataramanan, A. K. Srivastava, A. Hahn, and S. Zonouz, “Measuring and enhancing microgrid resiliency against cyber threats,” _IEEE Transactions on Industry Applications_ , vol. 55, no. 6, pp. 6303–6312, 2019.
  * [9] C.-W. Ten, J. Hong, and C.-C. Liu, “Anomaly detection for cybersecurity of the substations,” _IEEE Transactions on Smart Grid_ , vol. 2, no. 4, pp. 865–873, 2011.
  * [10] V. Chandola, A. Banerjee, and V. Kumar, “Anomaly detection: A survey,” _ACM computing surveys (CSUR)_ , vol. 41, no. 3, pp. 1–58, 2009.
  * [11] R. Moghaddass and J. Wang, “A hierarchical framework for smart grid anomaly detection using large-scale smart meter data,” _IEEE Transactions on Smart Grid_ , vol. 9, no. 6, pp. 5820–5830, 2017.
  * [12] D. Gupta, M. Gupta, S. Bhatt, and A. S. Tosun, “Detecting anomalous user behavior in remote patient monitoring,” in _IEEE IRI_ , 2021.
  * [13] G. Fenza, M. Gallo, and V. Loia, “Drift-aware methodology for anomaly detection in smart grid,” _IEEE Access_ , vol. 7, pp. 9645–9657, 2019.
  * [14] A. A. Cook, G. Mısırlı, and Z. Fan, “Anomaly detection for iot time-series data: A survey,” _IEEE Internet of Things Journal_ , vol. 7, no. 7, pp. 6481–6494, 2019.
  * [15] M. Adkisson, J. C. Kimmel, M. Gupta, and M. Abdelsalam, “Autoencoder-based anomaly detection in smart farming ecosystem,” _arXiv preprint arXiv:2111.00099_ , 2021.
  * [16] D. Gupta, O. Kayode, S. Bhatt, M. Gupta, and A. S. Tosun, “Hierarchical federated learning based anomaly detection using digital twins for smart healthcare,” _arXiv preprint arXiv:2111.12241_ , 2021.
  * [17] A. Barua, D. Muthirayan, P. P. Khargonekar, and M. A. Al Faruque, “Hierarchical temporal memory based machine learning for real-time, unsupervised anomaly detection in smart grid: Wip abstract,” in _2020 ACM/IEEE 11th International Conference on Cyber-Physical Systems (ICCPS)_ , 2020, pp. 188–189.
  * [18] G. Elafoudi, L. Stankovic, and V. Stankovic, “Power disaggregation of domestic smart meter readings using dynamic time warping,” in _2014 6th International Symposium on Communications, Control and Signal Processing (ISCCSP)_ , 2014, pp. 36–39.
  * [19] J. R. Ausmus, P. K. P. Sen, T. Wu, U. Adhikari, Y. Zhang, and V. Krishnan, “Improving the accuracy of clustering electric utility net load data using dynamic time warping,” in _2020 IEEE/PES Transmission and Distribution Conference and Exposition (T D)_ , 2020, pp. 1–5.
  * [20] D. M. Diab, B. AsSadhan, H. Binsalleeh, S. Lambotharan, K. G. Kyriakopoulos, and I. Ghafir, “Anomaly detection using dynamic time warping,” in _2019 IEEE International Conference on Computational Science and Engineering (CSE) and IEEE International Conference on Embedded and Ubiquitous Computing (EUC)_ , 2019, pp. 193–198.
  * [21] Z. Zheng, M. Zhou, Y. Chen, M. Huo, L. Sun, S. Zhao, and D. Chen, “A fused method of machine learning and dynamic time warping for road anomalies detection,” _IEEE Transactions on Intelligent Transportation Systems_ , pp. 1–13, 2020.
  * [22] Y. Zhou and J. Li, “Research of network traffic anomaly detection model based on multilevel autoregression,” in _2019 IEEE 7th International Conference on Computer Science and Network Technology (ICCSNT)_ , 2019, pp. 380–384.
  * [23] F. Schmidt, F. Suri-Payer, A. Gulenko, M. Wallschläger, A. Acker, and O. Kao, “Unsupervised anomaly event detection for cloud monitoring using online arima,” in _2018 IEEE/ACM International Conference on Utility and Cloud Computing Companion (UCC Companion)_ , 2018, pp. 71–76.
  * [24] A. H. Yaacob, I. K. Tan, S. F. Chien, and H. K. Tan, “Arima based network anomaly detection,” in _2010 Second International Conference on Communication Software and Networks_ , 2010, pp. 205–209.
  * [25] O. Solomentsev, M. Zaliskyi, and T. Gerasymenko, “Change-point detection during radar operation,” in _2016 IEEE First International Conference on Data Stream Mining Processing (DSMP)_ , 2016, pp. 295–298.
  * [26] J. Geng and L. Lai, “Bayesian quickest change point detection and localization in sensor networks,” in _2013 IEEE Global Conference on Signal and Information Processing_ , 2013, pp. 871–874.
  * [27] A. Mueen and E. Keogh, “Extracting optimal performance from dynamic time warping,” in _Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , ser. KDD ’16. New York, NY, USA: Association for Computing Machinery, 2016, p. 2129–2130. [Online]. Available: https://doi.org/10.1145/2939672.2945383
  * [28] S. Salvador and P. Chan, “Fastdtw: Toward accurate dynamic time warping in linear time and space,” in _KDD workshop on mining temporal and sequential data_. Citeseer, 2004.
  * [29] E. H. Budiarto, A. Erna Permanasari, and S. Fauziati, “Unsupervised anomaly detection using k-means, local outlier factor and one class svm,” in _2019 5th International Conference on Science and Technology (ICST)_ , vol. 1, 2019, pp. 1–5.
  * [30] R. Medico, “rob-med/awesome-ts-anomaly-detection,” Aug 2020.
  * [31] C. Hannon, D. Deka, D. Jin, M. Vuffray, and A. Y. Lokhov, “Real-time anomaly detection and classification in streaming pmu data,” in _2021 IEEE Madrid PowerTech_. IEEE, 2021, pp. 1–6.
  * [32] C. Ratanamahatana and E. Keogh, “Three myths about dynamic time warping data mining,” 04 2005.
  * [33] S. Grewal, M. Soni, and D. Jain, “Cyber security threats in synchrophasor system in wide area monitoring system,” _TELKOMNIKA Indonesian Journal of Electrical Engineering_ , vol. 15, pp. 436–444, 09 2015.
  * [34] A. Sundararajan, T. Khan, A. Moghadasi, and A. I. Sarwat, “Survey on synchrophasor data quality and cybersecurity challenges, and evaluation of their interdependencies,” _Journal of Modern Power Systems and Clean Energy_ , vol. 7, no. 3, pp. 449–467, 2019.


Reynolds number scaling of influence of boundary layers on the
global behavior of laboratory quasi-Keplerian flows

Eric M. Edlund* and Hantao Ji
Princeton Plasma Physics Laboratory

September 1, 2015

We present measurements of quasi-Keplerian flows in a Taylor-Couette device
that identify the boundary conditions required to generate near-ideal flows
that exhibit self-similarity under scaling of the Reynolds number. These
experiments are contrasted with alternate boundary configurations that result
in flows that progressively deviate from ideal Couette rotation as the
Reynolds number is increased. These behaviors are quantitatively explained in
terms of the tendency to generate global Ekman circulation and the balance of
angular momentum fluxes through the axial and radial boundary layers.

*current email address<EMAIL_ADDRESS>
##  1 Introduction

That the global properties of an extended system may be mapped to the
boundaries is an idea that has found success in holographic theories of
general relativistic systems [1, 2], and in magnetically confined plasmas [3,
4]. We report on a similar behavior observed in incompressible hydrodynamic
flows in a Taylor-Couette apparatus where it is observed that certain
characteristics of the global flow are largely dictated by the boundaries.
This finding is particularly relevant for experiments that examine quasi-
Keplerian (QK) flows, that is, rotation satisfying $0<q<2$ where
$q=-d\ln\Omega/d\ln r$, $\Omega$ is the fluid angular velocity and $r$ is the
radial coordinate, as models of astrophysical systems, namely accretion disks.
Numerous recent studies have commented on the hydrodynamic stability of such
systems [5, 6, 7, 8], with extensions to magnetohydrodynamics in electrically
conducting fluids [9, 10, 11, 12, 13].

While there is some disagreement between studies as to whether hydrodynamic
turbulence can be induced in QK flows, the balance seems to lean toward the
negative, at least insofar as incompressible turbulence is considered, and
points to the important role of magnetohydrodynamic effects in astrophysical
systems. However, while it is known that QK flows are linearly stable it
remains unknown whether there exists a nonlinear transition to turbulence,
even for incompressible hydrodynamic systems. Some experiments [5, 6, 7] and
simulations [14] indicate that such a transition is not likely, while others
present evidence that suggests that a subcritical transition may exist [8] and
some simulations find significant transient growth of perturbations that may
allow for nonlinear effects to enter [15]. Fluid experiments in other regimes
of operation that are not astrophysically- elevant have observed bi-stability
[16, 17], suggesting that should a similar mechanism exist for QK systems then
a subcritical pathway to turbulence may explain angular momentum transport in
accretion disks [18]. We show in this work that the influence of the
boundaries is intimately connected to the global structure of flows in Taylor-
ouette experiments and, by extension, is also related to the tendency of these
systems to generate and sustain turbulence.

One of the long-standing challenges of Taylor-Couette experiments in the quest
to understand angular momentum transport in astrophysically-relevant flows has
been the parasitic presence of Ekman circulation (secondary circulation)
induced by the mismatch between the fluid velocity and the solid body rotation
of the axial boundaries. A significant reduction in Ekman circulation has been
realized in experiment by using axial boundaries that are split into multiple
rings capable of differential rotation. Under particular boundary conditions,
azimuthal velocity profiles of the fluid can be generated that very nearly
match that of ideal Couette rotation [5, 7, 19, 20], the rotation profile that
is expected in the absence of axial influences for a constant radial flux of
angular momentum, and has been observed to hold over a wide range of Reynolds
numbers [20]. In contrast, studies in the “classical” configuration where the
axial boundaries rotate with the outer cylinder have shown performance that
further deviates from ideal Couette as the Reynolds number is increased [21].
Such trends are revealing of whether these systems are dominated by boundary
interactions or internal dynamics, a distinction with important consequences
for the applicability of such experiments to interpretation of astrophysical
systems, especially at large Reynolds numbers. First, through the experiments
reported here we identify two necessary criteria that define constraints on
the boundary configurations that allow near-ideal flows to develop. We then
discuss the competing roles of radial (Stewartson) boundary layers and axial
(Ekman) boundary layers, from which we develop a model that describes the
quantitative departure of the rotation profiles from ideal Couette flow as a
function of the angular momentum fluxes through the boundaries.

##  2 Experimental Apparatus

A Taylor-Couette (TC) device is a system of coaxial cylinders that rotate
independently of each other with the experimental fluid region between. The TC
apparatus used in these studies, called the Hydrodynamic Turbulence Experiment
(HTX), is a modified version of the classical device in that the axial
boundaries in HTX are segmented to allow differential rotation across the
boundaries [7]. The inner cylinder radius is $r_{1}=6.9$ cm and the outer
cylinder radius is $r_{2}=20.3$ cm. The inner radius and outer radius of the
independent rings are defined by the parameters $r_{3}$ and $r_{4}$,
respectively. The axial length of the experimental volume is $L=39.8$ cm,
giving an aspect ratio of $\Gamma=L/(r_{2}-r_{1})=2.97$ (see Fig. 1).
Corresponding components on the top and bottom are driven by the same motor so
that the system is up-down symmetric. The angular velocities of the inner
cylinder, outer cylinder and rings are identified by $\Omega_{1}$,
$\Omega_{2}$ and $\Omega_{3}$, respectively. Rotary encoders report the speed
of the motors to the control system.

A laser Doppler velocimeter (LDV) diagnostic system is used to measure the
local, azimuthal velocity ($v_{\theta}$), which in the experiments reported
here were measured at the midplane of the device. The LDV system is calibrated
by measuring fluid flow in solid body for which, after spin-up, the only
velocity component is $v_{\theta}$, which is a unique function of the motor
speeds. For these studies we define the global shear Reynolds number as
$\text{Re}_{s}=r_{g}^{2}\Delta\Omega/\nu$, where $r_{g}=(r_{1}r_{2})^{1/2}$ is
the geometric- ean radius, $\Delta\Omega=\Omega_{1}-\Omega_{2}$ and $\nu$ is
the kinematic viscosity, approximately $1\times 10^{-6}$ m2/s for water.

Figure 1:  (color online) Illustration of the HTX device at PPPL showing the
segmented axial boundaries (yellow, blue and green segments), the inner
cylinder in black, the outer cylinder in green and the working fluid in dark
blue.

##  3 Observations

The core observation of this work is summarized in Fig. 2, presenting flows
under two different values of normalized shear ($q=1.8$ and $q=1.5$) and three
different boundary conditions: Split, Optimized, and Ekman. For the $q=1.8$
cases these configuration speeds, reported as
$\Omega_{1}$-$\Omega_{3}$-$\Omega_{2}$, are multiples of 350-350-50 RPM
(Split), 350-185-50 RPM (Optimized), and 350-50-50 RPM (Ekman). For the
$q=1.5$ these are 250-250-50 RPM (Split), 250-140-50 RPM (Optimized), and
250-50-50 RPM (Ekman). The reference profile for these studies is described by
ideal Couette rotation, defined in terms of angular velocity as
$\Omega_{\text{C}}(r)=\Omega_{A}+\Omega_{B}(r_{g}/r)^{2}$, where
$\Omega_{A}=(r_{2}^{2}\Omega_{2}-r_{1}^{2}\Omega_{1})/(r_{2}^{2}-r_{1}^{2})$,
$\Omega_{B}=r_{g}^{2}(\Omega_{1}-\Omega_{2})/(r_{2}^{2}-r_{1}^{2})$. As a
function of azimuthal velocity the Couette solution is
$v_{\text{C}}(r)=r\Omega_{\text{C}}$. It is interesting that while the Ekman
and Split configurations exhibit progressive departure from ideal Couette as
the Reynolds number is increased, the shape of the Optimized cases is nearly
invariant with respect to scaling of the Reynolds number.

Figure 2:  (color online) Scaled measurements of $v_{\theta}$ at various
$\text{Re}_{s}$ for (a) $q=1.8$ and (b) $q=1.5$.

###  3.1 Necessary conditions for near-ideal flows

Recognizing that the axial boundaries in the Split cases will tend to increase
the angular momentum flux to the bulk, and conversely for the Ekman cases, we
begin by considering the ansatz that the balance of angular momentum fluxes
across the axial boundaries determines the deviation from ideal Couette. A
continuous variation of $\Omega_{3}$ from $\Omega_{2}$ to $\Omega_{1}$ would
then suggest that there is some intermediate state for which the profile must
pass close to ideal Couette rotation. We now address whether there is a method
for predicting this optimal value of $\Omega_{3}$.

The transition between the bulk flow and the walls of the TC device, which may
move at very different speeds, occurs over thin boundary layers. The structure
of these boundary layers depends on whether we are considering the balance of
forces in the axial or radial directions. The bulk flow transitions to the
axial boundary speeds over Ekman boundary layers whose thickness scales like
$\delta_{\text{E}}\sim r\text{Re}_{b}^{-1/2}$, and in the ideal model of
radial boundaries, over Stewartson boundary layers that scale like
$\delta_{\text{S}}\sim r\text{Re}_{b}^{-1/4}$, where
$\text{Re}_{b}=r^{2}\Omega_{b}/\nu$ is a local Reynolds number particular to
the boundary location $r$ with boundary speed $\Omega_{b}$ [22, 23, 24]. The
fluxes of angular momentum across these boundary layers, being inversely
proportional to the boundary layer thickness, do not scale proportionally, and
hence we anticipate that self-similarity of the global properties need not be
preserved as the Reynolds number is scaled. Thus, while the observed variation
in the shape of the Ekman and Split profiles is expected, it comes as some
surprise that the Optimized profiles are effectively independent of the
Reynolds number.

The experiments reported here showed no significant temporal variation in the
mean values, hence, it can be stated that in steady-state the net flux of
angular momentum into these flows must sum to zero, that is,
$\Phi_{1}+\Phi_{2}+2\Phi_{z}=0$ where $\Phi_{1}$ and $\Phi_{2}$ are the fluxes
integrated over the inner and outer cylinders, respectively, and $\Phi_{z}$ is
the axial flux integrated over each axial boundary (with inward towards the
fluid defined as a positive flux). The ideal Couette profile is the response
to a constant radial flux of angular momentum, implying that
$\Phi_{2}=-\Phi_{1}$ and that the axial fluxes of angular momentum are
everywhere zero. Such a state can be imagined in a TC system with free-slip
conditions on the axial boundaries, or with a continuously variable boundary
that can perfectly match the ideal Couette profile, conditions that will
result in vanishing stress at the axial boundaries. The radial flux of angular
momentum under these conditions (which we call $\Phi_{\text{C}}$, the Couette
flux) has a magnitude equal to

$\Phi_{\text{C}}=\Phi_{0}\frac{r_{g}^{2}}{r_{2}^{2}-r_{1}^{2}}\text{Re}_{s},$
(1)

where $\Phi_{0}=4\pi\rho\nu^{2}L$. As an aside, it is interesting to note that
$\Phi_{C}=\Phi_{0}\;\text{Re}_{s}$ when the radius ratio ($r_{2}/r_{1}$) is
equal to the golden ratio.

Even though there is, in general, a large mismatch between fluid and axial
boundary speeds, we find that the flows under the Optimized boundary
conditions very closely approximate the ideal Couette profile, and therefore
have a nearly constant radial flux of angular momentum. We now show that
rather than satisfying the condition that the axial flux of angular momentum
everywhere vanishes, these flows satisfy a much weaker constraint: it is the
surface integral of the axial angular momentum flux that vanishes, that is,
$\Phi_{z}\approx 0$.

To calculate the axial flux of angular momentum, being proportional to
$d\Omega/dz$, we make the assumption that
$d\Omega/dz\approx\Delta\Omega_{E}/\delta_{E}$, where
$\Delta\Omega_{E}=\Omega_{b}-\Omega$ is the difference of the boundary angular
velocity ($\Omega_{b}$) and the fluid angular velocity just inside the Ekman
boundary layer ($\Omega$), and that
$\delta_{E}=\alpha_{E}\left(\nu/\Omega_{b}\right)^{1/2}$ represents the
thickness of an Ekman boundary layer with an unknown numerical constant
$\alpha_{E}$. Defining the normalized quantities $\omega=\Omega/\Delta\Omega$
(recalling, $\Delta\Omega=\Omega_{1}-\Omega_{2}$),
$\omega_{b}=\Omega_{b}/\Delta\Omega$, $s_{1}=r_{1}/r_{g}$ and
$s_{2}=r_{2}/r_{g}$, we have

$\Phi_{z}=\Phi_{0}\frac{r_{g}}{2\alpha_{E}L}\text{Re}_{s}^{3/2}\int_{s_{1}}^{s_{2}}\left(\omega_{b}-\omega\right)\omega_{b}^{1/2}s^{3}ds.$
(2)

That the radial flux of angular momentum under ideal Couette rotation [Eq. 1]
scales like $\text{Re}_{s}$, whereas $\Phi_{z}$ scales like
$\text{Re}_{s}^{3/2}$ means that unless steps are taken to force the integral
in Eq. 2 to zero, there will always exist a Reynolds number beyond which the
axial flux will overwhelm the Couette flux and cause the flow to depart from
ideal rotation, regardless of aspect ratio. From Eq. 2 we can immediately
conclude that TC devices of the “classical” geometry with the axial boundaries
co-rotating with either the inner cylinder or outer cylinder will always have
non-zero $\Phi_{z}$ and can therefore never be a good model for astrophysical
systems at sufficiently high Reynolds numbers to be of interest. In looking
for the conditions under which ideal Couette flow may be generated, we set
$\omega=\omega_{\text{C}}$ and search for the boundary conditions (the
$\omega_{b}$) that cause the integral of Eq. 2 to vanish.

Figure 3:  (color online) Experiments conducted in $q=1.8$ flows showing
$v_{\theta}$ as a function of $\Omega_{3}$ for the (a) HTX and (b) wide-ring
geometries, progressing from low ring speeds (bottom curves) to higher ring
speeds (upper curves). In panels (c) and (d) the experimental $\chi^{2}$
(dashed green), for the measured profiles relative to the ideal Coutte profile
for these boundary conditions, is plotted against the calculated axial flux of
angular momentum from Eq. 2 (red, increasing function) and the pressure
differential from Eq. 3 (blue, decreasing function) for the HTX and wide-ring
geometries, respectively. $\Phi_{z}$ is normalized by $\Phi_{C}$ and $\Delta
p$ is normalized by the kinetic energy density for ideal Couette flow. The
vertical dotted lines indicate the zeros of $\Delta p$ and $\Phi_{z}$.

Figure 3 summarizes fluid velocity measurements over a scan of $\Omega_{3}$
for $q=1.8$ flows for two axial boundary configurations: the HTX configuration
with $r_{3}=7.8$ cm and $r_{4}=13.5$ cm, and a wide-ring configuration with
rings that span the entire radial gap. These experiments show that the HTX
configuration produces an exceptional match to the ideal Couette profile for a
narrow range of ring speeds centered about $185$ RPM, with very low
fluctuation levels spanning the entire gap [7]. Importantly, the boundary
conditions that generate near-ideal flows for the HTX configuration include
the zero crossing of $\Phi_{z}$ [Fig. 3(c)]. However, the wide-ring case also
has a zero in $\Phi_{z}$, near $85$ RPM [Fig. 3(d)], yet its flows never
resemble ideal Couette, indicating that the vanishing of $\Phi_{z}$ is a
necessary but not sufficient condition for achieving near-ideal flows.

Figure 4:  (color online) Same as Fig. 3 except for $q=1.5$.

The role of the axial boundaries in dictating global performance has also been
interpreted through the competition of pressures from the bulk flow and from
the boundary flow that is viscously coupled to the axial boundaries [25, 26].
The tendency to drive secondary flows was shown to be consistent with whether
the bulk pressure is larger than (Ekman circulation) or smaller than (anti-
Ekman circulation) the boundary pressure. Following the intuition motivated by
these simulations, we define a function $\Delta p$ that characterizes the
average pressure difference between ideal Couette rotation and boundary
rotation,

$\Delta p=\frac{\rho r_{g}^{3}\Delta\Omega^{2}}{\Delta
r}\int_{s_{1}}^{s_{2}}\int_{s_{1}}^{s}\left(\omega_{\text{C}}^{2}-\omega_{b}^{2}\right)\>s^{\prime}\>ds^{\prime}\>ds.$
(3)

For the wide-ring configuration the zeros of $\Delta p$ and $\Phi_{z}$ are
widely separated, meaning that no circumstance exists where pressure balance
and zero net axial flux can be simultaneously satisfied. In contrast, the HTX
configuration has zeros of these functions that are nearly coincident and fall
within the operating range in which small fluctuations were observed in Ref.
[7]. Other boundary configurations with $q=1.5$ also have nearly coincident
zeros of $\Phi_{z}$ and $\Delta p$, in remarkably good agreement with
experiments; see Fig. 4. Note that the requirement of coincident or nearly-
coincident zeros represents an effective third constraint. These studies
suggest that the requirements of having nearly coincident zeros of $\Phi_{z}$
and $\Delta p$ for the generation of ideal flows may represent necessary and
sufficient conditions. Additional research conducted over a wider range of
geometries and shear conditions will provide a stronger test of this
hypothesis.

###  3.2 Departure of flows from ideal Couette

While the coincident vanishing of $\Phi_{z}$ and $\Delta p$ define what may be
necessary and sufficient conditions for ideal flows to develop, they do not
reveal how QK systems should behave under non-optimized conditions. A general
model of the fluid response to forcing by the boundaries must account for the
fluxes across both axial and radial boundary layers. As there does not exist a
theory of Stewartson boundary layers under the conditions of QK rotation at
large Reynolds numbers, we cannot rely on results derived from a linear
analysis of perturbative differences in rotational speeds [24], especially for
experimental conditions where the Stewartson boundary layers are turbulent
[7]. Instead, we assume a generalized scaling of the Stewartson boundary layer
thickness of the form
$\delta_{\text{S}}=\alpha_{\text{S}}r(r/L)\text{Re}_{b}^{\beta}$, taking the
exponent of the Reynolds number ($\beta$) as a free-parameter to be determined
from the observed scalings in Fig. 2. While the numerical factors
$\alpha_{\text{E}}$ and $\alpha_{\text{S}}$ are introduced through the
definitions of the boundary layer thicknesses we cannot measure the boundary
layers directly, and therefore we must interpret their meaning as a measure of
the effectiveness of angular momentum flux.

We begin this analysis by considering the Ekman configurations from Fig. 2
where we note that $v_{\theta}$ transitions to solid body rotation at
$\Omega_{2}$ near the outer cylinder, implying that $\Phi_{2}\approx 0$. With
a jump in azimuthal velocity (a negative $\Delta v$) occurring over a
Stewartson layer at the inner cylinder, the appropriate form for $\Phi_{1}$ is

$\Phi_{1}=-\Phi_{0}\frac{c_{1}}{2\alpha_{S}}\>\text{Re}_{s}^{\beta}\>\frac{L\Delta
v}{\nu},$ (4)

where $c_{1}=(r_{1}\Omega_{1}/r_{2}\Delta\Omega)^{\beta}$ is a result of
converting from a representation of the boundary layer thickness that depends
on the local $\text{Re}_{b}$ to a global $\text{Re}_{s}$. We approximate these
flows with a piecewise function of the form
$v_{\theta}(r_{a}<r<r_{t})=v_{\text{C}}+\Delta v$ and
$v_{\theta}(r_{t}<r<r_{2})=r\Omega_{2}$, where $r_{t}$ is the transition
radius, similar to the identification of separate flow regions in the recent
work of Nordsiek et al. [21]. Substituting this rotation profile into Eq. 2
(using $\omega=v_{\theta}/r\Delta\Omega$) and using Eq. 4 for the radial flux
at $r_{1}$, we solve the global balance of angular momentum,
$\Phi_{1}+2\Phi_{z}=0$, for the dimensionless $\Delta
v^{\prime}_{\text{neg}}=\Delta v/r_{g}\Delta\Omega$ as a function of
$\text{Re}_{s}$. Noting that $\Phi_{z}$ has two parts given this
representation, our solution has three terms and is of the form

$\Delta
v^{\prime}_{\text{neg}}=\frac{f_{1}}{\alpha\>c_{1}\>\frac{L^{2}}{2r_{g}^{2}}\>\text{Re}_{s}^{\beta-1/2}+f_{2}},$
(5)

where $\alpha=\alpha_{\text{E}}/\alpha_{\text{S}}$, and the integral
expressions $f_{1}$ and $f_{2}$ are defined as

$f_{1}=\int_{s_{1}}^{s_{t}}\left(\omega_{b}-\omega_{\text{C}}\right)\omega_{b}^{1/2}s^{3}ds$
(6)

$f_{2}=\int_{s_{1}}^{s_{t}}\omega_{b}^{1/2}s^{2}ds.$ (7)

Figure 5:  (color online) Scaling of the velocity differential $\Delta v$ for
the three groups of boundary conditions from Fig. 2 for (a) $q=1.8$ and (b)
$q=1.5$. The dashed lines in the figure are the solutions to Eqs. 5 and 9 with
$\beta=0.08$ and $\alpha_{\text{E}}/\alpha_{\text{S}}=15$.

Equation 5 is a transcendental expression in $\Delta v$ since the limits of
the integrals in $f_{1}$ and $f_{2}$ are functions of $r_{t}$ which itself
depends on $\Delta v$. Equation 5 has only two free parameters: $\beta$ and
$\alpha=\alpha_{E}/\alpha_{S}$. Comparison of Eq. 5 with experimental
measurements of $\Delta v$ is presented in Fig. 5 for the $q=1.8$ and $q=1.5$
cases, where values of $\beta\approx 0.08$ and $\alpha_{E}/\alpha_{S}\approx
15$ provide the best fit to the experimental measurements.

This analysis can be extended to the case of the Split configuration with
positive $\Delta v$ by accounting for the different boundary conditions and a
slightly more complex flow structure. The most important feature of the
positive $\Delta v$ cases is that there is a nearly linear decrease in $\Delta
v$ across the gap, where $\Delta v$ at $r_{2}$ is defined to be equal to
$\epsilon\Delta v$. In the $q=1.8$ cases $\Delta v$ decreases by about $60\%$
across the gap ($\epsilon=0.4$), and by $80\%$ for $q=1.5$ ($\epsilon=0.2$),
nearly independent of Reynolds number. The radial flux at $r_{2}$ in terms of
the $\epsilon$ reduced $\Delta v$ is

$\Phi_{2}=\Phi_{0}\>\frac{c_{2}}{2\alpha_{S}}\text{Re}_{s}^{\beta}\>\frac{L\>\epsilon\>\Delta
v}{\nu},$ (8)

where $c_{2}=(r_{2}\Omega_{2}/r_{1}\Delta\Omega)^{\beta}$. The solution for
positive, normalized $\Delta v$ is

$\Delta
v^{\prime}_{pos}=\frac{f_{3}}{\epsilon\>\alpha\>c_{2}\>\frac{L^{2}}{2r_{g}^{2}}\>\text{Re}_{s}^{\beta-1/2}+f_{4}},$
(9)

where, similar to Eq. 5, we have

$f_{3}=\int_{s_{t}}^{s_{2}}\left(\omega_{b}-\omega_{\text{C}}\right)\omega_{b}^{1/2}s^{3}ds$
(10) $f_{4}=\int_{s_{t}}^{s_{2}}\left[\frac{s_{2}-\epsilon
s_{t}-(1-\epsilon)s}{s_{2}-_{t}}\right]\omega_{b}^{1/2}s^{2}ds.$ (11)

For the Split cases, the term $s_{t}$ defines the transition point from
increasing to decreasing $\Delta v$, and the terms in square brackets in Eq.11
derive from the variation of $\Delta v$ with radius (note that if $\Delta v$
were constant across the gap then $\epsilon=1$ and the term in brackets
reduces to unity, recovering Eq. 7). Equation 9 is compared against the
measured $\Delta v^{\prime}$ in Fig. 5 and is also found also to be in good
agreement with experiment when using the same values of
$\alpha_{E}/\alpha_{S}$ and $\beta$ derived from the analysis of the Ekman
cases.

The success of this model in reproducing the general features of the
departures from ideal Couette is remarkable given that it represents a simple
accounting of the angular momentum fluxes across the boundaries, ignoring
completely the complex internal dynamics of these flows, especially at large
Reynolds number where a substantial fraction of the volume exhibits turbulent
fluctuations [7]. What disagreement exists between experiment and theory
should not overshadow the success of this model in reproducing the general
trend of the velocity deviations, the relative amplitude of the positive and
negative $\Delta v$ cases, and the change in the amplitude of the positive
$\Delta v$ cases between $q=1.8$ and $q=1.5$. While the model presented here
uses a single set of parameters for both values of $q$, better fits to the
data can be found if we let these parameters depend on $q$. It is also
possible the these parameters may depend on geometry and Reynolds number in a
way that is not captured by the power-law scaling used here.

Studies of the measured fluctuation levels arising from the Stewartson
boundary layers in Ref. [7] found that the transition to turbulent Stewartson
boundary layers was consistent with the Taylor model of low Reynolds number
centrifugal instability. So while cannot offer a precise explanation for why
$\beta$ should take a value close to $0.1$, we reiterate that the inferred
scaling applies to the angular momentum flux which is a combination of
variations in the boundary layer thickness and an effective viscosity within
these layers that is modified by turbulent fluctuations. It is beyond the
scope and ability of these experiments to identify the separate scalings of
the turbulent viscosity and the thickness of the Stewartson boundary layers,
though perhaps future experiments will be able to provide greater insight into
this problem.

##  4 Conclusions

We conclude with an outlook to future physical and numerical experiments. The
ratio of $\Phi_{z}$ and $\Phi_{\text{C}}$ explored earlier can be recast as a
constraint on the aspect ratio $L/\Delta r$ as a function of the Reynolds
number. The relative contribution of $\Phi_{z}$ can be made arbitrarily small
by increasing the axial size of the system, so that to make
$\Phi_{z}/\Phi_{\text{C}}\approx 10^{-2}$ for an Ekman configuration, for
example, one would need an aspect ratio of order $100$ at a Reynolds number of
$10^{6}$. Outside of using very large aspect ratios, one can employ mechanical
advantages as we have done in HTX, such as extensions of the inner and outer
cylinder or independent rings. Despite the desire to simplify the mechanical
design of such modified TC devices, further analysis shows that the only way
to achieve simultaneous zeros in $\Phi_{z}$ and $\Delta p$ is through a design
with at least one independent ring. Greater insight into the balance of forces
in quasi- eplerian flows could be explored in future experiments in modified
Taylor Couette devices like HTX. In particular, further exploration of the
Reynolds number dependence on these deviations, the dependence of the model
parameters on $q$, and the structure of the Ekman cells would tell us much
more about the influence of the boundaries on global behavior.

We have shown that nearly ideal flows exhibit profile shape invariance under
scaling of the Reynolds number, an effect we interpret through the dual
conditions of vanishing axial angular momentum flux and vanishing pressure
differential that are nearly simultaneously satisfied, offering predictive
capability for selecting optimized boundary conditions and in experimental
design. The strongest piece of evidence in support of this model is the
prediction of self-similarity of the profiles with respect to scaling of the
Reynolds number only for cases in which the axial flux of angular momentum and
the pressure differential vanish nearly coincidentally, a prediction in
excellent agreement with the observations presented in Fig. 2. It should be
reiterated that while the experiments for Optimized flows show very small
departures from ideal Couette, and whose interpretation is congruent with the
coincident zeros in $\Delta p$ and $\Phi_{z}$, this does not strictly prove
that these criteria represent true sufficient conditions. An interesting
problem for future studies would be to measure how the departures from ideal
Couette depend on the “distance” between these zeros.

It is interesting to note that in all cases presented here the local axial
fluxes of angular momentum may be quite large given the substantial difference
in speed between the bulk fluid and the axial boundaries, and that it is only
on summation over the axial boundaries that zero net axial flux is realized
for the Optimized cases. Recalling that multiple experiments [7, 20, 19] and
simulations [25, 27] have observed a nearly uniform axial structure through
the bulk of the fluid volume, the existence of the large axial fluxes
naturally raises the question: What allows the bulk flows to depart from the
solid body rotation forced by the boundaries? Intuition based on the Taylor-
Proudman theorem for Rayleigh-stable flows, that is $dv_{\theta}/dz\approx 0$,
would suggest that the bulk should tend to follow the boundary. However, it
should be recalled that the Ekman boundary layers do not need to satisfy the
Taylor- roudman theorem because the axial gradient of $v_{\theta}$ is balanced
with a viscous diffusion of vorticity over the scale of the boundary layer
thickness. Thus, the nearly uniform axial structure suggests that the large
fluxes of angular momentum that must be present due to the existence of the
the Ekman boundary layers are redistributed locally, perhaps by small Ekman
cells that do not extend far into the bulk of the fluid or even in the
boundary layers themselves. As has been shown in prior studies, the free-shear
layer that develops at the interface between the differentially rotating rings
on the axial boundaries is expected to be centrifugally unstable and may help
to enhance the redistribution of angular momentum locally [13].

Another problem for future experiments, both physical and numerical, will be
to explore in greater detail how the very large, local fluxes of angular
momentum are redistributed near the axial boundaries with only minor effect on
the bulk flow in the Optimized configurations. Non-optimized boundaries, like
the Split and Ekman configurations, show progressive departure from the ideal
Couette flow as the Reynolds number is increased, in agreement with the
expectations of a dominant axial flux of angular momentum. A good test for
numerical experiments will be to accurately model the global behavior of the
mean flows in TC experiments by properly accounting for the angular momentum
fluxes from the boundaries. We believe that the boundary layer scalings
presented here may enable simulations to bootstrap to larger effective
Reynolds numbers by using specified boundary flux models, thus bypassing the
need for very fine grids to resolve the boundary layer structure directly.
With the combination of recent simulations of QK flows that can attain
Reynolds numbers of the order of $10^{5}$ [14, 28], though with axially-
periodic boundary conditions, a specified-flux boundary model may allow
simulations to reach something resembling experimental conditions of large
Reynolds number TC flows.

## Acknowledgement

We thank J. Goodman for sharing his thoughts on angular momentum transport and
boundary layers in our studies, and E. Schartman, E. Gilson, and P. Sloboda in
helping to keep the experiments running smoothly. This work was supported by
the U.S. Department of Energy, Fusion Energy Sciences under contract DE-
AC02-09CH11466 through the Center for Momentum Transport & Flow Organization
in Plasmas and Magentofluids (CMTFO).

## References

  * [1] A. Adams, P. M. Chesler, and H Liu.  Holographic Turbulence.  Phys. Rev. Lett., 112:151602, 2014.
  * [2] F. Carrasco, L. Lehner, R. C. Myers, O. Reula, and A. Singh.  Turbulent flows for relativistic conformal fluids in $2\mathbf{+}1$ dimensions.  Phys. Rev. D, 86:126006, 2012.
  * [3] J. E. Rice, J. W. Hughes, P. H. Diamond, Y. Kosuga, Y. A. Podpaly, M. L. Reinke, M. J. Greenwald, Ö. D. Gürcan, T. S. Hahm, A. E. Hubbard, E. S. Marmar, C. J. McDevitt, and D. G. Whyte.  Edge temperature gradient as intrinsic rotation drive in alcator $c$-mod tokamak plasmas.  Phys. Rev. Lett., 106:215001, 2011.
  * [4] O. D. G urcan, P. H. Diamond, C. J. McDevitt, and T. S. Hahm.  A simple model of intrinsic rotation in high confinement regime tokamak plasmas.  Phys. Plasmas, 17:032509, 2010.
  * [5] H. Ji, M. Burin, E. Schartman, and J. Goodman.  Hydrodynamic turbulence cannot transport angular momentum effectively in astrophysical disks.  Nature, 444:343, 2006.
  * [6] E. Schartman, H. Ji, M. Burin, and J. Goodman.  Stability of quasi-Keplerian shear flow in a laboratory experiment.  Astron. Astrophys., 543:A94, 2012.
  * [7] E. M. Edlund and H. Ji.  Nonlinear stability of laboratory quasi-keplerian flows.  Phys. Rev. E, 89:021004, 2014.
  * [8] M. S. Paoletti and D. P. Lathrop.  Angular momentum transport in turbulent flow between independetly rotating cylinders.  Phys. Rev. Lett., 106:024501, 2011.
  * [9] M. Seilmeyer, V. Galindo, G. Gerbeth, T. Gundrum, F. Stefani, M. Gellert, G. R udiger, M. Schultz, and R. Hollerbach.  Experimental evidence for magnetorotational instability in a taylor-couette flow under the influence of a helical magnetic field.  Phys. Rev. Lett., 113:024505, 2014S.
  * [10] F. Stefani, T. Gundrum, G. Gerbeth, G. Rüdiger, M. Schultz, J. Szklarski, and R. Hollerbach.  Experimental evidence for magnetorotational instability in a taylor-couette flow under the influence of a helical magnetic field.  Phys. Rev. Lett., 97:184502, 2006.
  * [11] D. R. Sisan, N. Mujica, W. A. Tillotson, Y. Huang, W. Dorland, A. B. Hassam, T. M. Antonsen, and D. P. Lathrop.  Experimental observation and characterization of the magnetorotational instability.  Phys. Rev. Lett., 93:114502, 2004.
  * [12] A. H. Roach, E. J. Spence, C. Gissinger, E. M. Edlund, P. Sloboda, J. Goodman, and H. Ji.  Observation of a Free-Shercliff-Layer Instability in Cylindrical Geometry.  Phys. Rev. Lett., 108:154502, 2012.
  * [13] E. J. Spence, A. H. Roach, E. M. Edlund, P. Sloboda, and H. Ji.  Free magnetohydrodynamic shear layers in the presence of rotation and magnetic field.  Phys. Plasmas, 19:056502, 2012.
  * [14] R. Ostilla-Mónico, R. Verzicco, S. Grossman, and D. Lohse.  Free magnetohydrodynamic shear layers in the presence of rotation and magnetic field.  J. Fluid Mech., 748:R3, 2014.
  * [15] S. Maretzke, B. Hof, and M. Avila.  Free magnetohydrodynamic shear layers in the presence of rotation and magnetic field.  J. Fluid Mech., 742:254, 2014.
  * [16] D. S. Zimmerman, S. Andrés Triana, and D. P. Lathrop.  Free magnetohydrodynamic shear layers in the presence of rotation and magnetic field.  Phys. Fluids, 23:065104, 2011.
  * [17] S. G. Huisman, R. C. A. van der Veen, C. Sun, ad D. Lohse.  Free magnetohydrodynamic shear layers in the presence of rotation and magnetic field.  Nature Commun., 5:3820, 2014.
  * [18] G. Lesur, and P -Y. Longaretti.  Free magnetohydrodynamic shear layers in the presence of rotation and magnetic field.  Astron. Astrophys., 444:25, 2005.
  * [19] M. Burin, H. Ji, E. Schartman, R. Cutler, P. Heitzenroeder, W. Liu, L. Morris, and S. Raftopolous.  Reduction of Ekman circulation within Taylor-Couette flow.  Experiments in Fluids, 40:962, 2006.
  * [20] E. Schartman, H. Ji, and M. Burin.  Development of a Couette-Taylor flow device with active minimization of secondary circulation.  Rev. Sci. Instrum., 80:024501, 2009.
  * [21] F. Nordsiek, S. G. Huisman, R. C. A. van der Veen, C. Sun, D. Lohse, and D. P. Lathrop.  Azimuthal velocity profiles in Rayleigh-stable Taylor-Couette flow and implied axial angular momentum transport.  J. Fluid Mech., 774:342, 2015.
  * [22] D. A. Bennetts and L. M. Hocking.  On nonlinear Ekman and Stewartson layers in a rotating fluid.  Proc. R. Soc. A., 333:469, 1973.
  * [23] D. J. Baker.  Shear layers in a rotating fluid.  J. Fluid Mech., 29:165, 1967.
  * [24] K. Stewartson.  On almost rigid rotations.  J. Fluid Mech., 3:17, 1957.
  * [25] A. V. Obabko, F. Cattaneo, and P. F. Fischer.  The influence of horizontal boundaries on Ekman circulation and angular momentum transport in a cylindrical annulus.  Phys. Scr., T132:014029, 2008.
  * [26] O. Czarny, E. Serre, P. Bontoux, and R. Leuptow.  Interaction between Ekman pumping and the centrifugal instability in Taylor-Couette flow.  Phys. Fluids, 15:467, 2003.
  * [27] M. Avila.  Free magnetohydrodynamic shear layers in the presence of rotation and magnetic field.  Phys. Rev. Lett., 108:124501, 2012.
  * [28] L. Shi, M. Rampp, B. Hof, and M. Avila.  Free magnetohydrodynamic shear layers in the presence of rotation and magnetic field.  Comput. Fluids, 106:1, 2015.