printacmref=false ifaamas doi [AAMAS’20]Proc. of the 19th International
Conference on Autonomous Agents and Multiagent Systems (AAMAS 2020), B. An, N.
Yorke-Smith, A. El Fallah Seghrouchni, G. Sukthankar (eds.)May 2020Auckland,
New Zealand 2020 2020

0000-0001-6705-8694 USC Information Sciences Institute Marina Del Rey
California 90292 USC Information Sciences Institute Marina Del Rey California
90292 USC Information Sciences Institute Marina Del Rey California 90292
Pytho, LLC. Brooklyn New York USC Information Sciences Institute Marina Del
Rey California 90292

# Anchor Attention for Hybrid Crowd Forecasts Aggregation

Yuzhong Huang<EMAIL_ADDRESS>,  Andres Abeliuk<EMAIL_ADDRESS>,  Fred
Morstatter<EMAIL_ADDRESS>,  Pavel Atanasov<EMAIL_ADDRESS>and  Aram
Galstyan<EMAIL_ADDRESS>
###### Abstract.

Forecasting the future is a notoriously difficult task. To overcome this
challenge, state-of-the-art forecasting platforms are ”hybridized”, they
gather forecasts from a crowd of humans, as well as one or more machine
models. However, an open challenge remains in how to optimally combine
forecasts from these pools into a single forecast. We proposed anchor
attention for this type of sequence summary problem. Each forecast is
represented by a trainable embedding vector, and use computed anchor attention
score as the combined weight. We evaluate our approach using data from real-
world forecasting tournaments, and show that our method outperforms the
current state-of-the-art aggregation approaches.

###### Key words and phrases:

Aggregation; Crowd Sourcing; Embedding; Attention Model

##  1\. Introduction

The “wisdom of crowd effect” has been demonstrated as a successful approach in
diverse domains Surowiecki (2004), extending to complex problem solving tasks
such as reconstructing gene regulatory networks Marbach et al. (2012) or
geopolitical forecasts Tetlock (2017).

Forecasting the outcome of rare events is a notoriously difficult problem,
where even experts in the domain area relevant to the forecasting problem fail
to significantly outperform a random baseline Tetlock (2017). One way to
tackle these problems is to have a multitude of forecasters working on the
same problem, in effect creating a crowdsourcing framework around these
problems.

Recent attention has been devoted to producing geopolitical forecasting
systems using various crowdsourcing methods. These efforts were impelled by
the Intelligence Advanced Research Project Agency (IARPA), who in 2011
launched a research program called Aggregate Contingent Estimation (ACE) to
enhance the accuracy and precision of crowdsourced geopolitical forecasting
systems.111https://www.iarpa.gov/index.php/research-programs/ace The wisdom of
the crowd approach was able to generate fairly accurate forecasts across a
wide range of forecasting problems Tetlock and Gardner (2016)

Recent efforts have sought to hybridize these forecasting endeavors. In
hybridized forecasting, human forecasters are paired with machine models in
the effort of generating more accurate forecasts about these events. The
general idea is that humans are more capable of reasoning over qualitative
data (e.g., news articles), while machine models are capable of processing
numerical data at a scale beyond human capacity. Combining these sources
should optimally lead to an _aggregated forecast_ that is more accurate than
that produced by either source independently.

While assessing the quality of each source is important to understanding the
best way to combine them, the strategy by which these sources of input should
be combined is not clear. Previous work Ungar et al. (2012); Atanasov et al.
(2017); Satopää et al. (2014a) has identified linear combinations among
forecaster estimates. However, such an approach is not optimal, since it
assigns a single static weight to each forecaster towards a variety of
forecasting problems. In the model proposed in this work, the combination
weight is conditioned on the question, forecaster, and time, and is more
flexible than a single static weight per forecaster.

The key insight to this method is the anchor attention models learn a
representation that could recall the best human and machine forecasters for a
given question. In essence, the model assesses the quality of each source
based on its past performance. Specifically, our contributions are the
following:

  1. (1)

Propose Anchor Attention, a novel attention architecture that is not sensitive
to the last input element and more suitable for sequence summary tasks.

  2. (2)

Experiment results on real-world forecast tournaments which show Anchor
Attention outperforms linear baselinesAtanasov et al. (2017).

(a)

(b)

Figure 1. Mean probability of daily forecasts for two sample questions.

##  2\. Related work

The general framework of crowdsourced platforms for forecast generations
consists of two main components. First, elicitation techniques are developed
to acquired information from several human experts regarding the predictions
about future events. Eliciting crowd information techniques range from
prediction markets Wolfers and Zitzewitz (2004); Yeh (2006) to prediction
polls Mellers et al. (2015); Tetlock and Gardner (2016). Second, aggregation
algorithms are developed to generate a prediction about future events based on
all the information gathered from the crowd. Atanasov et al. (2017) presents
an experimental comparison between prediction markets and aggregated
prediction polls. The prediction market was more accurate than the simple mean
of forecasts from prediction polls. However, prediction polls outperformed
prediction markets when forecasts were aggregated using a linear weighted
combination with temporal subsetting and extremization. We will compare our
methods to those proposed by Atanasov et al. (2017), which are described in
detail in Section 5.1.

Even though linear aggregation methods are the most popular approach for
combining probability forecasts, Ranjan and Gneiting (2010) proved
theoretically that this approach is sub-optimal on its own. Thus, any
technique based on linear pooling requires a non-linear recalibration or
extremization to improve performance Baron et al. (2014). This has motivated
other methods for aggregating probabilities in a non-linear way. Satopää et
al. (2014a) introduces a model-based aggregator that maps the probabilities
into the log-odds space, allowing to model the probabilities using the Normal
distribution. This approach has been extended in Satopää et al. (2014b) to
exploit the temporal correlation between forecasts, similar to standard time-
series approaches like autoregressive integrated moving average (ARIMA), where
future forecasts are a linear combination of past forecasts.

Pertaining to neural network based aggregation methods, to the best of our
knowledge, there are two relevant papers. Gaunt et al. (2016) implemented a
network of three fully connected layers, tailored for tasks such as data
annotation or labeling, to aggregate crowd opinions. Given the nature of the
task in mind, their method is unable to learn the temporal relation and
historical performance of the forecasters. Besides, their method requires that
each participant responds to all questions. Nebbione et al. (2019) proposed a
method capable of learning from historical performance, but it is only able to
rank individual forecasts without aggregating them.

The attention model proposed by Vaswani et al. (2017) is widely used in NLP
communities as an efficient and versatile sequence to sequence modeling
architecture. We borrow some elements from it when designing our anchor
attention sequence summary architecture.

##  3\. Crowd Forecasting Platform

In this work, we study forecasts about geopolitical events. These forecasts
are created on a forecasting platform we developed, Synergistic Anticipation
of Geopolitical Events (SAGE)Morstatter et al. (2019)222https://sage-
platform.isi.edu/. Our hybridized forecasting platform contains not only human
input, but also that from machine models. The system uses recruited
participants from Amazon Mechanical Turk. The forecasters are paid $16 for two
hours of work and are invited back each week to create new forecasts and
update existing ones. Forecasters are paid for their forecasting effort: they
are paid for participation, not for accuracy. The dataset consists of the
forecasts, as well as the forecasts generated by machine models on the same
questions. Upon publication, we will release an anonymized version of these
forecasts as well as the questions themselves.

Table 1. High-level statistics of our data. Name | Value
---|---
Number of Questions | 375
Number of Human Forecasters | 2240
Number of Machine Models | 3
Table 2. Summary statistics of the questions and forecasts in our dataset. STD
is the sample standard deviation.

Statistic | Min | Median | Mean (STD) | Max
---|---|---|---|---
Forecasts per Question | 9 | 192 | 265.2 (204.6) | 1029
Forecasts per User | 1 | 32 | 44.4 (63.2) | 1339
Users per Question | 8 | 108 | 173.6 (148.9) | 669
Days Question is Open | 1 | 42 | 53.3 (42.2) | 184

###  3.1. Forecasting Tournament Overview

All of the forecasting problems take the following form. Each question
pertains to categorical or ordinal geopolitical events. For example, “What
will be the long-term interest rate for Hungary in July 2018?” (ordinal) or
“Will UK’s Prime Minister Theresa May vacate office between 6 December 2018
and 15 January 2019?” (categorical). Figure 1 depicts the answer options for
both questions. Questions remain open (available for forecasting) for a
predetermined amount of time, except for “contingent stop” questions. The
question about Theresa May is contingent stop as it will close either on 15
January 2019 or in the case that she vacates office anytime before that.
During the time that the question is open, forecasters can generate forecasts.
Both users and machines enter their forecasts by assigning a probability to a
set of predefined answer options.

Each question contains 2 to 5 answer options. The answer options are non-
overlapping and are generated based upon historical data values for the
quantity in question. Finally, the resolution criteria provide the user with
an explicit link to resolution criteria, including the source that will be
used to decide the correct answer to the question. In this example, the
question is resolved based on data reported by the Organization for Economic
Cooperation and Development (OECD). Figure 1 illustrates two examples showing
the daily average forecast for reach question. Notice that in both questions
as they became closer to their resolution date, forecasters have more and new
information available, thus the consensus shifts towards the correct answer.
For this reason, a single static weight is not optimal.

###  3.2. Human Forecasting

The participants of the study were recruited through the Amazon Mechanical
Turk platform as well as through advertisements conducted on blogs and
periodicals related to geopolitical forecasting. Recruitment was done
exclusively before the beginning of the study. Forecasters were free to
provide forecasts to any number of open questions that were available. To
maintain user engagement, new questions were added to the platform each week.
The duration that a question was active varied in length from 3 days to 184
days, and some contingent stop questions closed early when the event occurred
before the close date. Forecasters are allowed and encouraged to update their
predictions as they receive new information over time. If a user updates
his/her forecasts over time, both old and new forecasts by the user on the
question are used.

Overall, we had $2,240$ participants producing $98,258$ forecasts over $375$
questions. Table 1 shows a summary of the data we collected. There was a wide
variance in the number of forecasts a given question received. This was due to
difficulty, which primarily stems from the amount of research required to
generate a forecast.

###  3.3. Machine Predictions

Given that questions covered a broad set of geopolitical topics and different
data sources, we used a general approach to forecast time series to produce
machine forecasts from the historical data relevant to each question.
Specifically, we used three types of machine models: 1) AutoRegressive
Integrated Moving Average (ARIMA) Hyndman and Khandakar (2008); 2) M4-Meta
Makridakis et al. (2019); and 3) Arithmetic Random Walk (RW).

##  4\. Methodology

###  4.1. Preliminary

We consider a set of $Q$ ($|Q|=q$) multiple-choice questions. For each
question $i\in Q$, it has $n_{i}$ answer options, and these answers may be
categorical or ordinal. There is exactly one correct answer $o_{i}$. Each
question has a start date and an end date. The system needs to make an
aggregated forecast each day between the start and end date. The questions are
posed to a set of $U$ ($|U|=u$) participants. Each participant $j\in U$
assigns probabilities for all answer options in a question of their choice,
denoted as vector $P_{ij}$. A forecaster may enter and update his/her estimate
at any time while the question is open. The same procedure applied to machine-
generated forecasts. In the rest of the paper, we use forecasters to refer to
both human forecasters and machine models.

The target is to recover the correct answer, which is represented as a one-hot
vector with weight on the correct answer option. This could be recovered from
a series of noisy observations coming from the human forecasters and machine
models. At time $t^{\prime}$, the aggregator needs to make a prediction
$\mathbb{O}_{t^{\prime}}=f(\mathbf{P_{t}}|t\leq t_{0})$. The number of
aggregated forecasts to be made is usually much smaller than the number of
individual forecasts.

As mentioned above, we consider several additional features when training the
model. These include static features and time-variant features. Examples of
static features include the question’s high-level category (e.g., politics,
economics), and the number of answer options. Examples of time-variant
features include linear and sinusoidal position embeddingVaswani et al.
(2017), the number of forecasters, and the variance of the average forecast.
Vaswani et al. (2017).

###  4.2. Attention Model

Attention modelVaswani et al. (2017) has shown to be an efficient and
versatile sequence modeling architecture. But to make it works in forecasts
aggregation task, some modification is needed, especially on its self-
attention layer. Otherwise, it’s prone to overfit and oscillating during
training.

####  4.2.1. Issue with Self Attention

Following the notation in Vaswani et al. (2017), we denote the input vector as
$X$, and three trainable weight matrices as $W^{Q},W^{K},W^{V}$. The common
practice in self-attention models is to use the last hidden state $H_{i}$ as a
representation of information in $X$ until time step $i$.

$Q=X\cdot W^{Q},K=X\cdot W^{K},V=X\cdot W^{V}.$ (1)
$\alpha_{ij}=align(Q_{i},K_{j})=\frac{exp(Q_{i}\cdot
K_{j}^{T}/\sqrt{d_{k}})}{\sum_{j^{\prime}}exp(Q_{i}\cdot
K_{j^{\prime}}^{T}/\sqrt{d_{k}})}.$ (2)
$H_{i}=(V_{1},V_{2},\cdots,V_{i})\cdot(\alpha_{i1},\alpha_{i2},\cdots,\alpha_{ii})^{T},$
(3)

We can see $H_{i}$ is most influenced by $X_{i}$ (the last input), more
specifically its projection $Q_{i}$, because $Q_{i}$ will be used as the query
vector to compute the alignment score with input from previous time steps. For
language tasks, this is desired, as the last token usually contains important
information (e.g., “?” implies a sentence is a question). In the task of
forecast aggregation, this is _not_ desired. Any forecast can be the last
forecast at the moment of making an aggregated forecast. The output of the
aggregation system should not be oversensitive to the last forecast.

####  4.2.2. Anchor Attention

Anchor attention use an anchor vector $A$, which is independent of input
sequence X, to replace the query vector $Q$. Here $A$ is the sentence
embedding of question text that captures the semantics of question. Now the
alignment score is calculated as:

$\alpha_{ij}=align(A,K_{j})=\frac{exp(A\cdot
K_{j}^{T}/\sqrt{d_{k}})}{\sum_{j^{\prime}}exp(A\cdot
K_{j^{\prime}}^{T}/\sqrt{d_{k}})}.$ (4)

With this modification, we didn’t give forecasts at a particular position any
advantages.

Figure 2. Overview of proposed model.

####  4.2.3. Walk through of the model

Denote $P_{i}\in R^{5}$ is i-th forecast for a question. $U_{i}$ is its
corresponding forecaster embedding. $F_{i}$ is its corresponding additional
feature vector. (eg: number of forecasters, and variance of the forecast)
$P_{i},U_{i},F_{i}$ are concatenated to formulate a rich representation of
that forecast. $W_{k},W_{v}$ are two trainable matrix, $K_{i},V_{i}$ are two
different projections of X, $K_{i}=X_{i}\cdot W^{K},V_{i}=X_{i}\cdot W^{V}$.
$A$ is the sentence embedding of that question. $I_{i}$ is the time embedding
of the moment making a aggregated forecast. $K_{i}$ and $[A,I_{i}]$ have the
same dimension. $\alpha_{i}\in R=K_{i}\cdot[A,I_{i}]^{T}$. Similiar to the
look backward attention Radford et al. ([n. d.]), if an aggregated forecast is
made at time $D_{j}$, only forecast occurred before $D_{j}$ will be selected.
Note that $i,j$ are index for individual forecasts and aggregated forecasts,
respectively. They are not aligned, for example. when the first aggregation is
made at $D_{1}$, it’s possible that there are 3 individual forecasts have
occurred. Then the first aggregation result
$O_{1}=\sigma(\sigma([\alpha_{1},\alpha_{2},\alpha_{3}])\cdot[V_{1};V_{2};V_{3}]\cdot
W_{o})$, where $W_{o}$ is an trainable matrix in fully connected layer. In
general, it could be described as
$O_{i}=\sigma(\sigma([\alpha_{1},\alpha_{2},\cdots\alpha_{k}])\cdot[V_{1};V_{2};\cdots
V_{k}]\cdot W_{o})$, where $P_{k}$ is the last forecasts occurred before
$D_{i}$.

###  4.3. Brier Score Loss

We adopt a loss function that is based on the Brier score Brier (1950), which
is naturally used to assess the quality of a probabilistic forecast. First, we
define the Brier score and continue to describe how we adapt it as a loss
function in our approach.

######  Definition  (Unordered Brier Score).

Given a forecast $P_{ij}=(p_{1},p_{2},\ldots,p_{n})$, and the actual outcome
$o_{i}\in\Delta^{n}$, we use the Brier score $B(p)$ as a measure of accuracy:

$B(p)=\sum_{i=1}^{n}\left(p_{i}-o_{i}\right)^{2}.$ (5)

While the above scoring is only for categorical questions, to address the fact
that some questions have ordered outcomes, i.e. $p_{k-1}$ is more favorable to
$p_{k-2}$ if the actual outcome index is $k$, we use the ordered Brier score,
which is a variant that uses the cumulative probabilities instead of the
densities Jose et al. (2009).

######  Definition  (Ordered Brier Score).

Applying the cumulative sum to the unordered Brier score yields the ordered
Bothrier score.

$B(p)=\frac{1}{n-1}\sum_{i=1}^{n-1}\left[\left(\sum_{j=1}^{i}p_{j}-\sum_{j=1}^{i}o_{j}\right)^{2}+\left(\sum_{j=i}^{n}p_{j}-\sum_{j=i}^{n}o_{j}\right)^{2}\right]$

(6)

Given a question $j$ of duration $T_{j}$, and daily aggregated forecasts
$\mathbb{P}_{t},t=1,2,\ldots,T_{j}$, the score of the question is the mean
daily Brier score (MDB), which is an unweighted average of the Brier score
across days within a question. Finally, the overall evaluation for the
forecast quality of a method is the mean of mean daily Brier scores (MMDB),
i.e., the unweighted average of MDB among all questions. We implemented MMDB
as our loss function. E The ordered scoring rule is used for ordered
categorical questions. Using the ordered Brier score enables the model to
learn the ordinal relation between answer options.

###  4.4. Embedding of question text and forecaster id

The nature of the question may influence the performance of forecasters. We
use word embeddings Mikolov et al. (2013) to represent the semantic meaning of
the question’s text. We use pre-trained embeddings from fastText Mikolov et
al. (2018). For each word $w_{i}$ in the question text, we find a vector
$v_{i}\in\mathbb{R}^{300}$ as its representation. Then we take the average of
the respective $v_{i}$ for each word in the question as the question
embedding, which will be used as anchor vector $A$ for that question.

To build an embedding that considers the forecasters’ identity and their
similarity, we add a trainable embedding as part of the network. The embedding
is randomly initialized with uniform noise between $(-\sqrt{3},\sqrt{3})$ with
zero mean, unit variance, and is trained jointly with the entire network. For
each forecaster id $u_{i}$, we learn a vector $U_{i}$ as its embedding.

###  4.5. Model Configuration Detail

This model has an adaptive learning rate decay and weight reset mechanism. If
the validation set loss continues to increase for 5 epochs, we multiply the
current learning rate by 0.95. If the validation set loss keeps increasing for
20 epochs, we reset the network weight to the last state having the lowest
validation set loss.

To avoid overfitting, we utilize dropout Srivastava et al. (2014), Leaky ReLU
Maas et al. (2013) in our implementation. Other configuration parameters are
listed below.

Table 3. Hyper parameters and training details. Name | Value
---|---
Embedding Dimension | 300
Key, Value Dimension | 64
Batch size | 256
Training epoch | 100 epochs, with weight reset
Learning rate |  $10^{-4}$, with adaptive decay
Optimizer | Adam

##  5\. Experimental Results

We utilized human and machine forecasts on a wide range of forecasting
problems from a hybrid forecast competition in our evaluation. Some statics of
this dataset are listed in Table 1.

###  5.1. Baseline Methods

We compare our approach to Atanasov et al. (2017), who proposed an aggregation
method using temporal decay, differential weighting based on past performance,
and extremization. The approach outlined by Atanasov et al. (2017) provides a
simple yet strong performance baseline that outperforms prediction markets for
distilling the wisdom of crowds. The approach is the weighted aggregation of
individual predictions for each question into a single forecast:

$\overline{p}_{t,q}\propto\sum_{t,i}d_{t}\times w_{t,i}^{\gamma}\times
p_{t,i,q},$

where, $\overline{p}_{t,q}$ is the weighted aggregate forecast for question
$q$ across all individual forecasts. $p_{t,i,q}$ coming from user $i$ at time
$t$ for question $q$; $d_{t}$ is a temporal decay parameter which gives more
weight to recent forecasts; $w_{t,i}$ are individual weights for each
forecaster $i$ and time $t$ based on their past performance; and $\gamma$ is
the exponent parameter for the individual weights. Finally, the aggregate
forecast $\overline{p}_{t}$ is extremized based on Karmarkar’s equation
Karmarkar (1978) by a linear transformation in the log-odds scale to move the
distribution away from the uniform distribution $1/a$, as follows:

$\log\left(\frac{(a-1)\hat{p}_{t,k}}{1-\hat{p}_{t,k}^{\alpha}}\right)=\alpha\left(\frac{(a-1)\overline{p}_{t,k}}{1-\overline{p}_{t,k}}\right),$

where $\hat{p}_{t}$ is the final forecast probability distribution; $\alpha$
is the extremization parameter; $a$ is the number of answer options, such that
$k\in\\{1,2,3,\cdots,a\\}$. For ordinal questions, we use the cumulative
probabilities instead of the densities.

We consider three variants based on Atanasov et al. (2017):

M0::

Unweighted average with temporal decay.

M1::

Weighted average with temporal decay and extremization.

M2::

Weighted average with temporal decay, differential weighting based on past
performance and extremization.

All parameters were trained using the training data: accuracy weights for each
participant were set to the inverse of their average Brier across the
questions they answered; the rest of the parameters were estimated using grid
search.

###  5.2. Performance Comparison

####  5.2.1. Average Brier score comparison

We perform 5-fold cross-validation. Cross-validation is done at the question
level, meaning that all of the forecasts for a particular question are
constrained within their respective fold. We record the mean and standard
deviation of the Brier score across folds, and the quantile range of all Brier
scores in Table LABEL:tab:scores. Our proposed model has the lowest average
Brier score and quantile scores, along with the lowest variance, suggesting
our model is very stable.

We performed a one-sided test between M2 (strongest baseline) and the proposed
anchor attention method. The null hypothesis is that the population mean Brier
score of the Attention method is not lower than the M2 method. Both methods
have $N=375$ questions. According to Table LABEL:tab:scores, we could compute
final test statistic is z=2.445 which translates to p = 0.0071. This p-value
allows us to reject the null hypothesis at alpha=0.01. This test concludes the
proposed anchor attention-based method gets a lower Brier score than M2
methods, so our results are significant.

Table 4. Results of 5-fold cross-validation across different methods. We report the mean Brier score as well as the quantile range (e.g., for our proposed anchor attention based model, the forecast at the 25th percentile has a Brier score of 0.037). Method | Mean (Variance) | 25% | 50% | 75%
---|---|---|---|---
M0 | 0.321 (0.075) | 0.126 | 0.249 | 0.432
M1 | 0.319 (0.100) | 0.088 | 0.212 | 0.422
M2 | 0.304 (0.103) | 0.073 | 0.192 | 0.419
Attention |  0.251 (0.071) | 0.037 | 0.173 | 0.382

####  5.2.2. Analysis of the Distribution of Brier Scores

Figure 3. The distribution of the Brier scores for each method.

Figure 3 shows a violin plot depicting the distribution of Brier scores of
each method. We evaluate the Brier score distribution of aggregated forecasts.
The width of each violin is proportional to the ratio of forecasts having that
Brier score. The dotted line represents 25%, 50%, and 75% quantile. We
observed that our proposed achieving the lowest average Brier scores and
quantile scores. Also, the highest (worst) Brier score is also comparable to
the baseline model, which means our model is robust.

Appendix Figures 9 and 10 provide alternative ways, that are independent of
Brier, to compare the methods based on calibration and discrimination (AUC and
ROC Curve) which are common in classification tasks. The deep aggregation
methods have better discrimination with an AUC of 0.9, while the baseline
methods have 0.85.

####  5.2.3. Percentile of aggregated forecast among all individual
forecasters

Next, we use the quantile metric Han and Budescu (2019) to compare the
performance of each ensemble method by their relative position on the
cumulative distribution of the corresponding individual scores. We rank each
aggregated forecast among all the individual forecasts. The percentile of the
ranking is plotted.

In Figure 4 we plot the distribution of the ranking of aggregated forecasts
among all individual forecasts. The smaller the ranking, the lower the Brier
score. We observe that for our proposed aggregation methods, the aggregated
forecasts are better than most individuals. For example, half of the time, the
proposed model beats 70% of the individual forecasters, while baseline methods
beat approximately 50% of the individual forecasters.

Figure 4. The distribution of ranks. A method that has more forecasts in lower
positions is better than one that has more forecasts in higher ones.

####  5.2.4. Performance over time

We also compare the performance of our methods as a function of how much time
a question has been open. As expected, the closer we are to the resolution
date, the better _all_ aggregation methods perform. However, our methods
converge to the correct solution faster. Figure 5 shows the average Brier
score across all questions as a function of time.

Figure 5. Average Brier score across all questions along the forecast
progress, for each method. Average Brier score over progress of question,
normalized from 0 to 1.

We normalize the progress in the forecast period as 0 to 1 and calculate the
average Brier score across all questions along with the forecast progress. We
use 100 bins and the results are smoothed with Gaussian kernel, sigma = 5.
Figure 5 shows the average Brier score across all questions along with the
forecast progress, for each method. As expected, the longer a question is
being open (i.e., closer to the resolution date), all the aggregation methods
perform better. However, we observed that our models quickly reduce their
Brier score at the early stages of the question, and its overall Brier score
is also lower. The insight is our model could learn to find similarities
between questions and forecasters, so when a new question arrives, it could
quickly find the right forecasters to trust.

###  5.3. Visualization of Question Embedding

Figure 6. Question Embedding with t-SNE dimension reduction

From the first figure, we find questions are not uniformly distributed in the
embedding space, they form various clusters. So we made the second figure zoom
in to the red rectangle region in the first figure. We find those closely
located questions are actually similar, they are all asking about the closing
price of some subject on different dates. It demonstrates that our model could
find a good representation of questions, and is aware of the similarity
between them.

###  5.4. Visualization of User Embedding

Figure 7. User embedding with t-SNE dimension reduction

We find the forecasters are also not uniformly distributed, they roughly
formed two clusters. Interestingly, two machine models, Auto ARIMA and
Arithmetic RW are closely located in the embedding space, while they are
intrinsically similar (Arithmetic RW is a special case of Auto ARIMA model).
If we define good forecasters as its top 10% Brier score is smaller than 0.1,
we find these good forecaster are mostly located in the upper cluster. It
demonstrates that our model could learn a good representation of forecasters.

###  5.5. Analysis of Attention Scores

Figure 8 shows forecasts’ attention scores against their Brier scores. We
observe that attention score is negatively (coef=-0.089, p=1.2e-195)
correlated to Brier score, which means good forecasts (low Brier) have higher
weights in aggregation. Also, we observe there are no points in the upper
right corner, which means forecasts with a high Brier score do not receive a
high attention score. This suggests that the model is learning representations
that are useful to distinguish good forecasts from bad ones.

Figure 8. Relation between attention score to Brier score. The trend line
shows a downward trend, indicating that higher attention scores are given to
forecasts with lower Brier scores.

##  6\. Conclusion

In this work, we proposed a novel approach based upon deep neural networks to
aggregate a hybrid set of forecasts. Using data from real-world forecasting
competitions, we have demonstrated that this approach outperforms the current
state-of-the-art in forecast aggregation. Moreover, we showed preliminary
evidence that through learning an embedding of question and forecaster, the
proposed approach learns to pay attention to more accurate forecasts while
simultaneously combining human forecasters with machine models (See Appendix
Table 5).

While this work is undeniably focused on identifying quality forecasters in a
geopolitical event forecasting context, the results here can be applied in
other settings. For example, this approach could be used to identify which
workers are most appropriate to choose for a given task in a micro tasking
environment. Future work will be to better understand the implications of this
method in other environments and to adapt our approach to identify quality
workers with a minimal amount of data.

This research is based upon work supported in part by the Office of the
Director of National Intelligence (ODNI), Intelligence Advanced Research
Projects Activity (IARPA), via 2017-17071900005. The views and conclusions
contained herein are those of the authors and should not be interpreted as
necessarily representing the official policies, either expressed or implied,
of ODNI, IARPA, or the U.S. Government. The U.S. Government is authorized to
reproduce and distribute reprints for governmental purposes notwithstanding
any copyright annotation therein.

##  Appendix A Ablation Study

This section is devoted to answering the following questions

  1. (1)

To what extent the machine forecasts help the aggregation?

  2. (2)

Does knowing the past performance of forecasters help?

We restrict the analyses to the subset of questions that have both machine and
human forecasts available (even in configurations that only use humans or
machine models). This subset is different from the set of questions used in
previous experiments. We conducted a 2x3 factorial design experiment. The
first factor is whether or not we can identify (with an id) when different
forecasts are produced by the same forecaster; the second factor is the type
of forecasts use as the input of the aggregation: only human forecasts, only
machine forecasts and both.

Figure 5 compares the Brier scores for this experiment. First, the model’s
performances are improved with the addition of forecaster id, suggesting that
the model learns to accurately weigh forecasters based on their past
performance. Second, with humans only or machine models only, the score is not
as good as having both, however, this is only true if the model can
differentiate a human forecast from a machine based on the ids.

Table 5. Performance results for different inputs. Feature | Forecaster type | Brier score
---|---|---
None | Human | 0.303
None | Machine | 0.319
None | Human + Machine | 0.303
Forecaster id | Human | 0.297
Forecaster id | Machine | 0.318
Forecaster id | Human + Machine | 0.290

##  Appendix B Calibration and Discrimination

Calibration refers to the ability to make forecasts that coincide with the
observed empirical frequencies of events being predicted. Calibration error
would be zero if the events with predicted probabilities p actually occur
about p percentage of the time. Therefore, any deviation from the diagonal
line suggests poor calibration. Figure 9 shows the calibration curves for each
method. The x-axis is divided into 10 bins, each spanning 10% on the
probability scale. The y-axis denotes the observed proportion of events that
fall in each bin.

Figure 9. Calibration Plots

Discrimination refers to the capacity to distinguish between true positive
outcomes and true negative ones. We use the AUC (Area Under The Curve) ROC
(Receiver Operating Characteristics) curve to evaluate the discrimination of
the models. The ROC curve plots the probability of a true positive against
that of a false positive. Perfect resolution means all of the probability mass
is under the curve. Figure 10 shows the ROC curves for all methods.

Figure 10.  ROC Curves

## References

  * (1)
  * Atanasov et al. (2017) Pavel Atanasov, Phillip Rescober, Eric Stone, Samuel A Swift, Emile Servan-Schreiber, Philip Tetlock, Lyle Ungar, and Barbara Mellers. 2017\.  Distilling the wisdom of crowds: Prediction markets vs. prediction polls.  Management science 63, 3 (2017), 691–706.
  * Baron et al. (2014) Jonathan Baron, Barbara Mellers, Philip Tetlock, Eric Stone, and Lyle Ungar. 2014.  Distilling the wisdom of crowds: Prediction markets vs. prediction polls.  Decision analysis 11, 2 (2014), 133–145.
  * Brier (1950) Glenn W Brier. 1950\.  Verification of forecasts expressed in terms of probability.  Monthly weather review 78, 1 (1950), 1–3.
  * Gaunt et al. (2016) Alex Gaunt, Diana Borsa, and Yoram Bachrach. 2016.  Training deep neural nets to aggregate crowdsourced responses. In Proceedings of the Thirty-Second Conference on Uncertainty in Artificial Intelligence. AUAI Press, 242251\.
  * Han and Budescu (2019) Ying Han and David Budescu. 2019.  A universal method for evaluating the quality of aggregators.  Judgment and Decision Making 14, 4 (2019), 395\.
  * Hyndman and Khandakar (2008) Rob Hyndman and Yeasmin Khandakar. 2008.  Automatic Time Series Forecasting: The forecast Package for R.  Journal of Statistical Software, Articles 27, 3 (2008), 1–22.  https://doi.org/10.18637/jss.v027.i03
  * Jose et al. (2009) Victor Richmond R Jose, Robert F Nau, and Robert L Winkler. 2009\.  Sensitivity to distance and baseline distributions in forecast evaluation.  Management Science 55, 4 (2009), 582–590.
  * Karmarkar (1978) Uday S Karmarkar. 1978\.  Subjectively weighted utility: A descriptive extension of the expected utility model.  Organizational behavior and human performance 21, 1 (1978), 61–72.
  * Maas et al. (2013) Andrew L Maas, Awni Y Hannun, and Andrew Y Ng. 2013\.  Rectifier nonlinearities improve neural network acoustic models. In Proc. icml, Vol. 30. 3.
  * Makridakis et al. (2019) Spyros Makridakis, Evangelos Spiliotis, and Vassilios Assimakopoulos. 2019.  The M4 Competition: 100,000 time series and 61 forecasting methods.  International Journal of Forecasting (2019).  https://doi.org/10.1016/j.ijforecast.2019.04.014
  * Marbach et al. (2012) Daniel Marbach, James C Costello, Robert Küffner, Nicole M Vega, Robert J Prill, Diogo M Camacho, Kyle R Allison, Andrej Aderhold, Richard Bonneau, Yukun Chen, et al. 2012\.  Wisdom of crowds for robust gene network inference.  Nature methods 9, 8 (2012), 796.
  * Mellers et al. (2015) Barbara Mellers, Eric Stone, Terry Murray, Angela Minster, Nick Rohrbaugh, Michael Bishop, Eva Chen, Joshua Baker, Yuan Hou, Michael Horowitz, et al. 2015\.  Identifying and cultivating superforecasters as a method of improving probabilistic predictions.  Perspectives on Psychological Science 10, 3 (2015), 267–281.
  * Mikolov et al. (2018) Tomas Mikolov, Edouard Grave, Piotr Bojanowski, Christian Puhrsch, and Armand Joulin. 2018\.  Advances in Pre-Training Distributed Word Representations. In Proceedings of the International Conference on Language Resources and Evaluation (LREC 2018).
  * Mikolov et al. (2013) Tomas Mikolov, Ilya Sutskever, Kai Chen, Greg Corrado, and Jeffrey Dean. 2013.  Distributed Representations of Words and Phrases and Their Compositionality. In Proceedings of the 26th International Conference on Neural Information Processing Systems - Volume 2 (NIPS’13). Curran Associates Inc., USA, 3111–3119.  http://dl.acm.org/citation.cfm?id=2999792.2999959
  * Morstatter et al. (2019) Fred Morstatter, Aram Galstyan, Gleb Satyukov, Daniel Benjamin, Andres Abeliuk, Mehrnoosh Mirtaheri, KSM Tozammel Hossain, Pedro Szekely, Emilio Ferrara, Akira Matsui, Mark Steyvers, Stephen Bennet, David Budescu, Mark Himmelstein, Michael Ward, Andreas Beger, Michele Catasta, Rok Sosic, Jure Leskovec, Pavel Atanasov, Regina Joseph, Rajiv Sethi, and Ali Abbas. 2019.  SAGE: A Hybrid Geopolitical Event Forecasting System. In Proceedings of the Twenty-Eighth International Joint Conference on Artificial Intelligence, IJCAI-19. International Joint Conferences on Artificial Intelligence Organization, 6557–6559.  https://doi.org/10.24963/ijcai.2019/955
  * Nebbione et al. (2019) Giuseppe Nebbione, Derek Doran, Srikanth Nadella, and Brandon Minnery. 2019.  Deep Neural Ranking for Crowdsourced Geopolitical Event Forecasting. In Machine Learning for Networking, Éric Renault, Paul Mühlethaler, and Selma Boumerdassi (Eds.). Springer International Publishing, Cham, 257–269.
  * Radford et al. ([n. d.]) Alec Radford, Karthik Narasimhan, Tim Salimans, and Ilya Sutskever. [n. d.].  Improving language understanding by generative pre-training.  ([n. d.]).
  * Ranjan and Gneiting (2010) Roopesh Ranjan and Tilmann Gneiting. 2010.  Combining probability forecasts.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 72, 1 (2010), 71–91.
  * Satopää et al. (2014a) Ville A Satopää, Jonathan Baron, Dean P Foster, Barbara A Mellers, Philip E Tetlock, and Lyle H Ungar. 2014a.  Combining multiple probability predictions using a simple logit model.  International Journal of Forecasting 30, 2 (2014), 344–356.
  * Satopää et al. (2014b) Ville A Satopää, Shane T Jensen, Barbara A Mellers, Philip E Tetlock, Lyle H Ungar, et al. 2014b.  Probability aggregation in time-series: Dynamic hierarchical modeling of sparse expert beliefs.  The Annals of Applied Statistics 8, 2 (2014), 1256–1280.
  * Srivastava et al. (2014) Nitish Srivastava, Geoffrey Hinton, Alex Krizhevsky, Ilya Sutskever, and Ruslan Salakhutdinov. 2014.  Dropout: A Simple Way to Prevent Neural Networks from Overfitting.  Journal of Machine Learning Research 15 (2014), 1929–1958.  http://jmlr.org/papers/v15/srivastava14a.html
  * Surowiecki (2004) James Surowiecki. 2004\.  The wisdom of crowds: Why the many are smarter than the few and how collective wisdom shapes business.  Economies, Societies and Nations 296 (2004).
  * Tetlock (2017) Philip E Tetlock. 2017\.  Expert political judgment: How good is it? How can we know? Princeton University Press.
  * Tetlock and Gardner (2016) Philip E Tetlock and Dan Gardner. 2016.  Superforecasting: The art and science of prediction.  Random House.
  * Ungar et al. (2012) Lyle Ungar, Barbara Mellers, Ville Satopää, Philip Tetlock, and Jon Baron. 2012\.  The good judgment project: A large scale test of different methods of combining expert predictions. In 2012 AAAI Fall Symposium Series.
  * Vaswani et al. (2017) Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Lawrence Jones, Aidan N. Gomez, Lukasz Kaiser, and Illia Polosukhin. 2017\.  Attention Is All You Need. In NIPS.
  * Wolfers and Zitzewitz (2004) Justin Wolfers and Eric Zitzewitz. 2004.  Prediction markets.  Journal of economic perspectives 18, 2 (2004), 107–126.
  * Yeh (2006) Puong Fei Yeh. 2006\.  Using prediction markets to enhance US intelligence capabilities.  Studies in Intelligence 50, 4 (2006), 137–149.


# Towards decoupling the effects of permeability and roughness on turbulent
boundary layers

D. D. Wangsawijaya1<EMAIL_ADDRESS>P. Jaiswal1,2  B.
Ganapathisubramani 1  1Aeronautical and Astronautical Engineering, University
of Southampton SO17 1BJ, United Kingdom 2Department of Mechanical Engineering,
University of Sherbrooke, Sherbrooke, QC, J1K 2R1, Canada

###### Abstract

Boundary layer flow over a realistic porous wall might contain both the
effects of wall-permeability and wall-roughness. These two effects are
typically examined in the context of a rough-wall flow, i.e., by defining a
“roughness” length or equivalent to capture the effect of the surface on
momentum deficit/drag. In this work, we examine the hypothesis of Esteban et
al. (2022), that a turbulent boundary layer over a porous wall could be
modelled as a superposition of the roughness effects on the permeability
effects by using independently obtained information on permeability and
roughness. We carry out wind tunnel experiments at high Reynolds number
($14400\leq Re_{\tau}\leq 33100$) on various combinations of porous walls
where different roughnesses are overlaid over a given permeable wall.
Measurements are also conducted on the permeable wall as well as the rough
walls independently to obtain the corresponding lengthscales. Analysis of mean
flow data across all these measurements suggests that an empirical formulation
can be obtained where the momentum deficit ($\Delta U^{+}$) is modelled as a
combination of independently obtained roughness and permeability lengthscales.
This formulation assumes the presence of outer-layer similarity across these
different surfaces, which is shown to be valid at high Reynolds numbers.
Finally, this decoupling approach is equivalent to the area-weighted power-
mean of the respective permeability and roughness lengthscales, consistent
with the approach recently suggested by Hutchins et al. (2023) to capture the
effects of heterogeneous rough surfaces.

###### keywords:

##  1 Introduction

Flow over porous walls comprises an extensive number of natural phenomena,
ranging from blood vessels to the atmospheric boundary layer (ABL) developing
over a forest canopy. The latter can be considered as a turbulent boundary
layer (TBL) developing over a porous wall. An example of porous walls,
constructed from packed spheres, is shown in figure 1. Here, both walls
approximately have the same permeability (resistance of a substrate to fluid
flows), but the wall shown in figure 1(a) has a roughness interface above the
porous wall, while in figure 1(b) this interface comprises of a flat surface.
Thus, a realistic representation of a porous wall is _both_ permeable _and_
rough to some extent; for the wall in figure 1(a), the effects of both
permeability and roughness have to be considered for characterisation of a TBL
developing over such wall.

(a)(b) Figure 1: Illustrations of realistic permeable walls with: (a) a rough
wall on the interface and (b) flat interface. Adapted from Rosti et al.
(2015).

###  1.1 Rough walls

The presence of rough walls increases skin friction from that of smooth walls
and thus the time-averaged streamwise velocity $U$ of TBLs developing over
rough walls can be written as a downward shift of the logarithmic region
$\Updelta U^{+}$ from the smooth wall velocity profile

$U^{+}=\frac{1}{\kappa}\ln\left[\frac{(y+d)U_{\tau}}{\nu}\right]+B-\Updelta
U^{+}=\frac{1}{\kappa}\ln(y+y_{0})^{+}+B-\Updelta U^{+}$ (1)

The viscous-scaled velocity is defined as $U^{+}\equiv U/U_{\tau}$, where
$U_{\tau}\equiv\sqrt{\tau_{w}/\rho}$ is the friction velocity, $\tau_{w}$ and
$\rho$ are the wall shear stress (WSS) and the density of fluid, respectively.
The logarithmic profile is defined as follows: $\kappa$ is the von Kármán
constant, $d$ is the zero-plane displacement, $\nu$ is the kinematic viscosity
of fluid, and $B$ is the log-law intercept. The logarithmic shift $\Updelta
U^{+}$, also known as the Hama roughness function (Hama, 1954), is defined as

$\Updelta
U^{+}=\frac{1}{\kappa}\ln\left(\frac{k_{s}U_{\tau}}{\nu}\right)+B-B_{FR}=\frac{1}{\kappa}\ln
k_{s}^{+}+B-B_{FR}$ (2)

where $k_{s}$ is the equivalent sand grain roughness and $B_{FR}=8.5$ is the
‘fully’ rough intercept of the velocity profile of sand grain roughness. It
should be noted that $k_{s}$ is a measure of roughness effect on the flow
relative to that of a uniform sand grain roughness (Nikuradse, 1933); it can
only be determined by flow measurements. The mean streamwise velocity profile
of a rough wall can therefore be written as a function of its $k_{s}$

$U^{+}=\frac{1}{\kappa}\ln\left(\frac{y+d}{k_{s}}\right)+B_{FR}$ (3)

###  1.2 Permeable walls

Earlier studies of TBLs developing over permeable walls involved various types
of such walls, namely: packed spheres (Zagni & Smith, 1976), perforated sheets
(Kong & Schetz, 1982), bed of grains (Zippe & Graf, 1983), and multi-layered
walls (Manes et al., 2009). Permeable walls have been found to increase drag
(i.e. skin friction coefficient $C_{f}\equiv 2(U_{\tau}/U_{\infty})^{2}$) from
that of solid, impermeable walls, which is attributed to the increase of
dissipation, momentum flux, and Reynolds shear stress on the interface between
the fluid and the substrates (Zagni & Smith, 1976; Shimizu et al., 1990).
Similar to the rough wall TBLs, the increase in $C_{f}$ of permeable walls is
also characterised by a downward shift in the logarithmic region from that of
a solid, smooth wall (Hahn et al., 2002; Efstathiou & Luhar, 2018). Thus, it
is probable that there is a “universal" parameter that characterise permeable
walls – possibly equivalent to $k_{s}$ for rough wall TBLs. It was suggested
by Manes et al. (2011) that the logarithmic region of a TBL developing over a
permeable wall scales on permeability $K$, and this can be further rearranged
to show a downward logarithmic shift,

$U^{+}=\frac{1}{\kappa}\ln\left(\frac{y+d}{\sqrt{K}}\right)+c_{1}=\frac{1}{\kappa}\ln(y+d)^{+}-\frac{1}{\kappa}\ln
Re_{K}+c_{1}$ (4)

where $Re_{K}\equiv\sqrt{K}U_{\tau}/\nu$. Previous studies observed a wide
range of magnitude of $\kappa$ (Breugem et al., 2006; Manes et al., 2011; Suga
et al., 2010), which is largely attributed to the low $Re$ at which these
studies were conducted and thus there was not enough separation between the
inner and outer layers of the wall-bounded flows (Manes et al., 2011). A more
recent experimental work by Esteban et al. (2022), conducted at a higher order
of magnitude of $Re$ ($2000\leq Re_{\tau}\leq 18000$), observed that
$\kappa=0.39$, similar to that of smooth and rough wall TBLs. It was further
hypothesised by Esteban et al. (2022) that $c_{1}$ in equation (4) is an
additive constant related to the blockage effect of a porous substrate, as a
realistic porous wall comprises of both permeable matrix and solid substrates
(see, for example, figure 1), whose size does not permit the full isolation of
permeability effect (Breugem et al., 2006). This blockage effect might be
represented by a roughness function $\Updelta U_{b}^{+}\equiv 1/\kappa\ln
k_{s_{b}}^{+}+B-B_{FR}$, similar to that of equation (2)

$\displaystyle U^{+}$
$\displaystyle=\frac{1}{\kappa}\ln(y+d)^{+}-\frac{1}{\kappa}\ln
Re_{K}-\frac{1}{\kappa}\ln
k_{s_{b}}^{+}+B_{FR}=\frac{1}{\kappa}\ln\left(\frac{y+d}{Re_{K}k_{s_{b}}}\right)+B_{FR}$
(5)

where subscript ‘$b$’ denote the blockage effect of the porous substrate. By
comparing equation (3) and (5), the ‘equivalent sand grain roughness’ $k_{s}$
of a porous wall (subscript ‘$p$’) can be defined as

$k_{s_{p}}=Re_{K}k_{s_{b}}$ (6)

In this study, we take the first steps towards exploring the possibility of
decoupling permeability and roughness effects on turbulent boundary layers. We
have constructed three different porous test surfaces where we maintain the
permeability ($Re_{K}$) but the blockage is systematically altered by adding
roughness on to a given permeable substrate. Detailed hot-wire and drag-
balance measurements are taken for the permeable, rough and the combination of
permeable-rough surfaces over a wide range of Reynolds numbers (ensuring the
flow is in fully-rough regime for all cases). We use the experimental data to
explore two different approaches for decoupling the two effects. First, we
extend the above-mentioned framework to include the effects of additional
roughness. Second, we use a power-mean averaging approach proposed for
heterogeneous roughness (Hutchins et al. 2023) to find an equivalent roughness
lengthscale that can capture the combined effect of permeability and
roughness.

##  2 Experimental Setup

(a)(b)(c)(d)

| $k$ | $s$ | $\varepsilon$ | $K$ | $A_{o}$ | $l$ | $w$ | $t$
---|---|---|---|---|---|---|---|---
Surface | [mm] | [mm] | [–] | [$\times 10^{-8}$ m2] | [–] | [mm] | [mm] | [mm]
P | 15 | 0.89 | 0.859 | 3.62 | – | – | – | –
R1 | 3.0 | – | – | – | 0.73 | 62 | 30 | 4
R2 | 3.5 | – | – | – | 0.81 | 70 | 53 | 5

Figure 2: (a) Illustration of a test surface laid inside the test section of
the boundary layer wind tunnel. Combination of porous-rough test surfaces: (b)
porous P–rough wall R1, and (c) porous P–rough wall R2. (d) Geometric
parameters of P, R1, and R2, where $k$ is the thickness of the surfaces. For
P: $s$ is the average pore size of the foam, $\varepsilon$ is the porosity,
and $K$ is the permeability, obtained from Esteban et al. (2022). For R1 and
R2: $A_{o}$ is the open area ratio of the mesh, $l$ is the length of the
longway of the mesh, $w$ is the length of the shortway of the mesh, and $t$ is
the width of the mesh strand, as illustrated in (c).

###  2.1 Facility

Measurements are conducted inside the closed return boundary layer wind tunnel
(BLWT) at the University of Southampton. The flow passes through a contraction
section of 6:1 ratio before entering the 12 m $\times$ 1.2 m $\times$ 1 m
(length $\times$ width $\times$ height) test section. The boundary layer
develops over the floor (bottom surface) of the BLWT. A 220 mm-long smooth
wall ramp is installed at the end of the contraction section (figure 2(a)
WAH1ga) to match the thickness of the test surfaces. This ramp marks the inlet
towards the test section of the BLWT and the streamwise datum ($x=0$). The
boundary layer is tripped by a 8.5 mm-wide, 0.4 mm-thick turbulator tape
attached on the ramp (figure 2(a) WAH1gb). The tunnel has the freestream
turbulence level of $\sigma_{u^{\prime}}/U_{\infty}\approx 0.1$%, where
$\sigma_{u^{\prime}}$ is the standard deviation of streamwise turbulent
fluctuation at the freestream and $U_{\infty}$ is the freestream velocity. The
BLWT is equipped with a water-cooled heat exchanger, maintaining the flow
temperature variation of 1% for the longest measurements ($\sim$8 hours).

###  2.2 Test surfaces

Case | $\delta$ | $\theta$ | $U_{\infty}$ | $U_{\tau}$ | $Re_{x}$ | $Re_{\theta}$ | $Re_{\tau}$ | $Re_{K}$ | $C_{f}$ | $\Pi$ | $\Updelta U^{+}$ | $k_{s}$ | $k_{s}^{+}$ | Sym.
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
| [mm] | [mm] | [ms-1] | [ms-1] | [$\times 10^{7}$] | [$\times 10^{4}$] | [$\times 10^{4}$] | [–] | [$\times 10^{-3}$] | [–] | [–] | [mm] | [–] |
Porous walls
| 201.11 | 24.55 | 19.84 | 1.06 | 1.18 | 3.29 | 1.44 | 13.62 | 5.72 | 0.27 | 12.01 |  | 560 |
| 199.96 | 24.45 | 25.45 | 1.36 | 1.52 | 4.21 | 1.85 | 17.57 | 5.74 | 0.27 | 12.66 |  | 723 |
P | 208.26 | 25.37 | 29.06 | 1.54 | 1.76 | 5.07 | 2.21 | 20.22 | 5.65 | 0.31 | 13.02 | 7.83 | 832 |
| 211.27 | 25.88 | 35.46 | 1.87 | 2.14 | 6.30 | 2.72 | 24.50 | 5.59 | 0.34 | 13.52 |  | 1008 |
| 221.25 | 26.43 | 40.05 | 2.13 | 2.39 | 7.18 | 3.20 | 27.52 | 5.67 | 0.24 | 13.81 |  | 1132 |
Rough walls
| 232.40 | 28.51 | 14.90 | 0.86 | 0.87 | 2.82 | 1.33 | – | 6.69 | 0.19 | 12.88 |  | 773 |
| 236.84 | 28.78 | 19.93 | 1.15 | 1.17 | 3.83 | 1.82 | – | 6.70 | 0.18 | 13.64 |  | 1041 |
R1 | 239.38 | 28.70 | 25.20 | 1.46 | 1.48 | 4.82 | 2.32 | – | 6.67 | 0.19 | 14.23 | 13.52 | 1310 |
| 243.78 | 28.84 | 29.58 | 1.71 | 1.73 | 5.68 | 2.77 | – | 6.69 | 0.18 | 14.65 |  | 1539 |
| 254.33 | 29.52 | 35.21 | 2.04 | 2.05 | 6.87 | 3.42 | – | 6.69 | 0.15 | 15.08 |  | 1821 |
| 241.11 | 29.41 | 15.10 | 0.89 | 0.89 | 2.96 | 1.44 | – | 7.01 | 0.19 | 13.44 |  | 963 |
R2 | 229.33 | 28.56 | 19.74 | 1.18 | 1.15 | 3.75 | 1.79 | – | 7.09 | 0.18 | 14.14 | 16.17 | 1263 |
| 229.34 | 28.25 | 24.90 | 1.49 | 1.46 | 4.68 | 2.27 | – | 7.12 | 0.18 | 14.74 |  | 1599 |
| 230.03 | 28.64 | 29.83 | 1.78 | 1.75 | 5.69 | 2.73 | – | 7.13 | 0.16 | 15.21 |  | 1920 |
Porous-rough walls
| 225.41 | 26.90 | 17.28 | 1.06 | 1.01 | 3.08 | 1.58 | 13.32 | 7.48 | 0.18 | 14.08 |  | 1235 |
| 218.12 | 27.17 | 21.47 | 1.31 | 1.26 | 3.90 | 1.92 | 16.74 | 7.49 | 0.15 | 14.66 |  | 1552 |
PR1 | 227.93 | 27.28 | 26.04 | 1.60 | 1.54 | 4.78 | 2.46 | 20.55 | 7.60 | 0.14 | 15.19 | 17.65 | 1906 |
| 226.64 | 27.33 | 31.49 | 1.94 | 1.85 | 5.74 | 2.94 | 24.65 | 7.60 | 0.09 | 15.65 |  | 2286 |
| 229.11 | 27.28 | 35.01 | 2.17 | 2.05 | 6.36 | 3.31 | 27.46 | 7.67 | 0.10 | 15.93 |  | 2547 |
| 219.42 | 27.29 | 15.06 | 0.95 | 0.90 | 2.79 | 1.41 | 12.27 | 7.94 | 0.13 | 14.23 |  | 1254 |
PR2 | 221.73 | 27.41 | 19.97 | 1.26 | 1.20 | 3.75 | 1.91 | 16.40 | 7.93 | 0.13 | 14.98 | 19.45 | 1677 |
| 224.41 | 27.61 | 25.20 | 1.59 | 1.52 | 4.77 | 2.45 | 20.76 | 7.98 | 0.13 | 15.58 |  | 2122 |
| 222.96 | 28.06 | 28.26 | 1.79 | 1.69 | 5.40 | 2.71 | 23.16 | 8.00 | 0.14 | 15.86 |  | 2368 |
Table 1: Summary of all test surfaces: baseline smooth walls (S), porous (P),
rough (R1 and R2), and porous-rough (PR1 and PR2) walls. The statistics are
obtained from HWA measurements: $\delta$ is the 99% boundary-layer thickness
and $\theta$ is the momentum thickness. Reynolds number definitions are:
$Re_{x}\equiv xU_{\infty}/\nu$, $Re_{\theta}\equiv\theta U_{\infty}/\nu$,
$Re_{\tau}\equiv\delta U_{\tau}/\nu$, and $Re_{K}\equiv\sqrt{K}U_{\tau}/\nu$.
Coefficient of friction is defined as $C_{f}\equiv
2(U_{\tau}/U_{\infty})^{2}$. The last column shows the symbol associated with
each test surface. Colours denote test cases with approximately matched
$Re_{\tau}$ and $Re_{K}\approx$ 14400 and 13 (), 18600 and 17 (), 23400 and
20.5 (), 27700 and 24 (), 33100 and 27.5 ().

Five test surfaces are constructed for this study: porous wall (denoted by
‘P’, figure 2(a) WAH1gc), two types of rough walls (‘R1’ and ‘R2’, figure
2b,c), and two combinations of rough walls on top of the porous wall (‘PR1’
and ‘PR2’, figure 2b,c). All test surfaces are assembled on the bottom surface
of the BLWT, downstream of the ramp and trip (figure 2a). The ramp and test
surfaces only cover the first 10.8 m-long part of the test section, while the
section from $x=10.8$ m to $x=12$ m comprises of a smooth wall (figure 2a).
All measurements are conducted at a constant $x$ location at a range
freestream velocity $U_{\infty}$, corresponding to a set of matched friction
Reynolds number $Re_{\tau}\equiv\delta U_{\tau}/\nu\pm 1600$ (here, $\delta$
is the 99% boundary-layer thickness) between P, R1, and PR1, as well as
another set of matched $Re_{\tau}$ for P, R2, and PR2. The first set of test
surfaces has 5 matched $Re_{\tau}$ cases ($Re_{\tau}\approx 14400$, 18600,
23400, 22900, and, 33100), while the second set has four matched cases
($Re_{\tau}\approx 14400$, 18600, 23400, and 27700). Details of each test
surfaces, including the statistics obtained from HWA and WSS measurements, are
given in table 1. Throughout this study, each test surface is denoted with a
symbol (see the last column in table 1). The same colour denotes surfaces with
matched $Re_{\tau}$ and $Re_{K}$ (within $Re_{K}\pm 10\%$), which runs from
lighter towards darker colour as $Re_{\tau}$ and $Re_{K}$ increases.

The porous walls are constructed from sheets of 15 mm-thick ($1060\leq
k^{+}\equiv kU_{\tau}/\nu\leq 2130$), 45 ppi (pores per inch) polyurethane
reticulated foam. The substrate has the average pore size of $s=0.89$ mm
($63\leq s^{+}\equiv sU_{\tau}/\nu\leq 126$) and permeability of $K=3.62\times
10^{-8}$ m2, corresponding to $13\leq Re_{K}\leq 27.5$. The rough walls are
constructed from two types of diamond-shaped, expanded metal mesh sheets
denoted by R1 and R2 (figure 2b,c). The mesh sheets have the thickness of
$k=3$ and 3.5 mm ($178\leq k^{+}\leq 476$) for R1 and R2, respectively, and
the open area $A_{o}$ (ratio of empty to total area in $xz$–plane) of 0.73 and
0.81. Both mesh sheets only cover a 1 m-wide portion (in $z$) of the working
section (figure 2(a) WAH1gd). Details of the relevant parameters regarding the
porous substrate and mesh geometries are given in figure 2(d).

###  2.3 Hot-wire anemometry

Hot-wire anemometry (HWA) measurements are conducted for all test surfaces
listed in table 1 by traversing the wire in wall-normal direction $y$ across
35–40 logarithmically-spaced points from the wall towards the freestream. All
measurements are conducted at the centreline of the tunnel at $x=8.8$ m from
the datum of the tunnel test section (figure 2(a) WAH1gf). A modified Dantec
55P05 single-sensor with a boundary-layer type probe (and an appropriate probe
support) are used with a Streamline Pro Constant Temperature Anemometer (CTA)
with an overheat ratio of 0.8. The sensor is a 5 $\upmu$m diameter, 1-mm long
tungsten wire soldered to the tip of the hot-wire prong, satisfying the
recommended wire length-to-diameter ratio of 200 (Ligrani & Bradshaw, 1987)
and corresponding to the viscous-scaled sensor lengths of $57\leq
l_{w}^{+}\equiv l_{w}U_{\tau}/\nu\leq 145$. The output signal is sampled from
the CTA at $f_{s}=30$ kHz, yielding viscous-scaled sampling interval of
$1.6\leq t^{+}\equiv U_{\tau}^{2}/(f_{s}\nu)\leq 10$. The sampling time
$T_{s}$ at each measurement point differs between test surfaces such that the
boundary-layer turnover time is maintained at $T_{s}U_{\infty}/\delta\geq
20000$ to allow convergence of turbulence energy spectra (Hutchins et al.,
2009). The sensor is traversed to the freestream and calibrated prior to and
after each measurement. The relation between freestream velocity $U_{\infty}$
and output voltage of the sensor $E_{\infty}$ is defined by King’s law
$E_{\infty}^{2}=C_{1}+C_{2}U_{\infty}^{C_{3}}$. Temperature compensation is
applied to the output signal to account for a slight sensor drift.

###  2.4 Wall shear stress (WSS) measurements

Parameter | PR1 | PR2
---|---|---
$k_{s_{pr}}/k_{s_{p}}$ | 2.25 | 2.48
$C$ | 0.17 | 0.15
$n$ | 1.42 | 1.54
(a)(b)(c)(d)(e) Figure 3: (a) $C_{f}$ as a function of $Re_{x}$ for all test
surfaces. Validation of the WSS measurements with smooth walls (): fitted to
the composite profile of Rodríguez-López et al. (2015) () and the analytical
solution of Monty et al. (2016) ( ). The same porous surface measured by
Esteban et al. (2022) (). Porous wall ($k_{s}=7.83$ mm) at constant
$Re_{k}\equiv kU_{\infty}/\nu$ ( ) and constant $k_{s}/x$ ( ), obtained from
Monty et al. (2016). (b) $U^{+}$ as a function of $y^{+}$ for matched
$Re_{\tau}\approx 27700$,  : $1/\kappa\ln y^{+}+B$. (c) $\Updelta U^{+}$ as
functions of $k_{s}^{+}$ (white-filled symbols, bottom axis) and $Re_{K}$
(filled symbols, top axis).  : $1/\kappa\ln k_{s}^{+}+B-B_{FR}$. Legends are
shown in table 1. (d) Velocity defect $U_{\infty}^{+}-U^{+}$ as a function of
$y/\delta$ for matched $Re_{\tau}\approx 27700$. (e) Decoupling of
permeability and roughness effect using methods described in §3.

The WSSs of all test surfaces shown in table 1 are measured using an in-house
floating element (FE) drag balance. The balance has a 200 mm $\times$ 200 mm
FE located at $x=8.6$ m (figure 2(a) WAH1ge), slightly upstream of the HWA
measurement station. A section of each test surfaces are cut according to the
size of the FE and mounted on top of the FE, leaving a 0.5 mm clearance
between the FE and its housing, allowing it to move freely for WSS
measurements. The housing is sealed to prevent airflow into the balance.
Measurements are acquired at $f_{s}=256$ Hz with 120 s sampling time for each
test surface. Calibration of the balance is performed before measurements by
loading a set of calibration weights to the balance via a pulley system. The
relationship between known calibration weights and time-averaged output
voltage (within $\pm 0.96\%$ uncertainties) from the balance is obtained by
fitting the calibration data into a first order polynomial.

Figure 3(a) shows the $C_{f}$ of all tests surfaces as a function of fetch
Reynolds number $Re_{x}\equiv xU_{\infty}/\nu$. For validation, $C_{f}$ of
smooth walls ($2930\leq Re_{\tau}\leq 12500$) are determined first using three
different methods: (i) direct WSS measurements using the FE drag balance (),
(ii) fitting velocity profiles obtained from HWA measurements to the composite
profile of Rodríguez-López et al. (2015) (), (iii) the analytical solution of
Monty et al. (2016) ( ). All three methods show a reasonable agreement with
each other to within 5%.

For the rest of the test surfaces P, R, and PR, $C_{f}$ are determined from
direct WSS measurements. Both rough walls (R1 and R2) have higher $C_{f}$ than
that of porous walls, which increase further as the walls are combined with
the porous walls (PR1 and PR2), as shown in figure 3(a). The logarithmic shift
$\Updelta U^{+}$ and $k_{s}$ for these surfaces are obtained by fitting the
mean profile to a modified method of Rodríguez-López et al. (2015), in which
the log-law shift term $\Updelta U^{+}$ is added, forming equation (1) while
maintaining the $U_{\tau}$ obtained from direct WSS measurements. The fit
yields $\kappa=0.39$ and $B=4.34\pm 0.08$ (consistent with those obtained for
TBLs developing over various porous walls measured by Esteban et al., 2022).
Porous wall $C_{f}$ (figure 3a) is compared to the analytical solution
proposed by Monty et al. (2016) over constant $k_{s}/x=8.9\times 10^{-4}$
($k_{s}=7.83$ mm). Additionally, $C_{f}$ obtained by Esteban et al. (2022) for
the same porous wall, $k_{s}/x=2.4\times 10^{-3}$ ($k_{s}=7.97$ mm) is also
compared with the analytical solution, both showing good agreement with these
solutions.

Figure 3(b) shows the vertical downward shift in the logarithmic region for
all test surfaces at a matched $Re_{\tau}\approx 27700$ (the highest matched
$Re_{\tau}$ for PR2 set and second highest for PR1) from a smooth wall TBL at
this matched $Re_{\tau}$. This shift, denoted by $\Updelta U^{+}$, is shown as
a function of $k_{s}^{+}$ in figure 3(c) for all test surfaces (white-filled
contours). All collapse to the logarithmic function $1/\kappa\ln
k_{s}^{+}+B-B_{FR}$ (equation 2), although this is given since $k_{s}$ is not
obtained independently from $\Delta U^{+}$. The magnitude of $k_{s}^{+}$
(table 1) exceeds the threshold defined by Flack & Schultz (2010),
$k_{s}^{+}\gtrsim 70$, ensuring that all test surfaces are ‘fully’ rough and
thus $\Updelta U^{+}$ depends solely on $k_{s}$. It should be noted that
$\Delta U^{+}$ here accounts for the total momentum deficit (from that of
smooth wall TBLs) in both porous and rough walls, and it shows that it is
possible to characterise porous wall with the same framework used for rough
wall characterisation (i.e. logarithmic shift that scales on $k_{s}$). The
effect of permeability on $\Updelta U^{+}$ is shown in figure 3(c) for test
surfaces comprises of porous walls (P, PR1, and PR2, filled contours).
Although it is clear that $\Updelta U^{+}$ is logarithmically scaled by
$Re_{K}$, universality, as shown in figure 3(c) with $k_{s}$, is not apparent.
This highlights the importance of accounting for the ‘blockage’ effect
(equation 5). At high range of $Re_{\tau}$ tested in this study, the outer
layer similarity is preserved for all test surfaces (i.e. velocity defect
profiles collapse beyond $y/\delta\gtrsim 0.2$, figure 3d), ensuring that a
given value of $k_{s}$ (for a surface) can be used to predict the drag at
higher Reynolds numbers.

##  3 Decoupling permeability and roughness effect

###  3.1 Additional blockage effect

We rewrite equation (6) for porous wall (case P) in a more general term

$k_{s_{p}}=Re_{K}f(k_{s_{b}})$ (7)

where $f(k_{s_{b}})$ is a function representing the blockage effect by the
substrate. As a roughness interface is added (case PR), a roughness term,
possibly a function of the characteristic lengthscale of roughness (i.e. rough
wall $k_{s}$) increases blockage

$k_{s_{pr}}=Re_{K}f(k_{s_{b}},k_{s_{r}})$ (8)

where subscripts ‘$r$’ denotes the rough wall and ‘$pr$’ the porous-rough
wall. When two substrates have an approximately matched permeability, then the
additional blockage effect can be approximated by
$f(k_{s_{b}},k_{s_{r}})/f(k_{s_{b}})\approx k_{s_{pr}}/k_{s_{p}}$. Case PR1
has 5 matched $Re_{K}$ with case P (within 5% difference, table 1). As the
$k_{s}$ for these surfaces are known from fitting $\Updelta U^{+}$ (§2.4), the
additional blockage is approximated by
$f(k_{s_{b}},k_{s_{r}})/f(k_{s_{b}})=2.25$ (figure 3e) within 10% error due to
$Re_{K}$ difference. This is interpreted as the rough wall R1 increases the
total blockage effect by approximately 2.25 times from that of porous wall.

At this point, the empirical formulation of $f(k_{s_{b}},k_{s_{r}})$ is
unknown. Let the roughness effect be represented solely by rough wall $k_{s}$

$f(k_{s_{b}},k_{s_{r}})=f(k_{s_{b}})Ck_{s_{r}}$ (9)

where $C$ is a constant. For an approximately matched $Re_{K}$, $C\approx
k_{s_{pr}}/(k_{sp}k_{s_{r}})$. We obtain $C=0.17$ for case PR1 (figure 3e),
and further test this formulation for the second roughness PR2. Here, there
are 4 matched $Re_{K}$ cases (within 10% difference, table 1). The total
blockage effect increases by approximately 2.5 times from the porous wall
(figure 3e), consistent with the increasing $C_{f}$ of case R2 from that of R1
(figure 3a). The constant $C$ appears to be similar to that obtained from R1,
$C=0.15$. This suggest that the formulation has the potential to be applicable
across different cases and is not dependent on the type of roughness.

###  3.2 Equivalent homogeneous roughness

A recent study by Hutchins et al. (2023) suggested that for a heterogeneous
rough wall (i.e. a rough wall constructed from various patches of homogeneous
roughness of roughness length scale $k_{s_{i}}$ covering an area $A_{i}$), the
equivalent homogeneous roughness length $k_{ehr}$ can be defined as:

$k_{ehr}=\left[\frac{1}{A}\sum_{i=1}^{N}k^{n}_{s_{i}}A_{i}\right]^{\frac{1}{n}}\quad\mathrm{where}\quad
A=\sum_{i=1}^{N}A_{i}$ (10)

In the present study, the drag penalty from all test surfaces can be
represented by the logarithmic shift $\Updelta U^{+}$ (figure 3b), suggesting
the possibility of characterising permeable walls with the same framework used
for rough wall TBLs. Thus, we might consider porous-rough walls as simply a
combination of two ‘homogeneous roughnesses’ overlaying each other. Equation
(10) is therefore reduced to

$k_{s_{pr}}=\left[k_{s_{p}}^{n}+k_{s_{r}}^{n}\right]^{\frac{1}{n}}$ (11)

For case PR1, $n$ is obtained by solving equation (11) with the known $k_{s}$
for each PR1, P, and R1, resulting in $n=1.42$. The same approach is applied
to case PR2, resulting in $n=1.54$, which is relatively consistent with that
obtained from case PR1.

It is still unknown, at this point, what the empirical formulation for the
roughness effect is. Present results suggest that it is possibly in terms of
an additional logarithmic shift similar to equation (5). Further, it is still
unclear what the empirical relation with the characteristic lengthscale of
roughness (rough wall $k_{s}$) is. For the two methods tested in §3, it is
neither understood what the physical interpretations of constants $C$ and $n$
are, nor whether these constants are universal for various types of rough
walls. We are unable to answer these questions in this study. However, present
results suggest a promising start towards decoupling permeability and
roughness effects.

##  4 Conclusions and future work

We conduct velocity and drag measurements of TBLs developing over various
porous-rough wall combinations, with the roughness effect systematically
varied while maintaining the permeability effect. Measurements are conducted
at relatively high Reynolds numbers($14400\leq Re_{\tau}\leq 33100$), with a
set of matched $Re_{\tau}$ (within 10%) and $Re_{K}$ ($13\leq Re_{K}\leq
27.5$) of each porous-rough combinations. Present results suggest that the
increase in drag over porous, rough, and porous-rough walls is characterised
by a downward shift in the logarithmic region $\Updelta U^{+}$ from that of
smooth wall TBLs and that the mean flow follows outer-layer similarity for
cases examined at these high Reynolds numbers (with substantial scale
separation). Further analysis shows that $\Updelta U^{+}$ from these surfaces
can be characterised by the roughness lengthscale $k_{s}$, suggesting the
possibility of characterising porous walls with the same framework used for
rough wall TBLs. Following the hypothesis of Esteban et al. (2022), it is
viable that for a porous-rough wall, permeability and roughness effects are
represented by $Re_{K}$ and a total blockage function
$f(k_{s_{b}},k_{s_{r}})$, respectively, with the total blockage consists of
the blockage effect from the porous substrate and additional roughness
interface above the porous wall. Present results suggest that the rough walls
increase the blockage of porous-rough walls by approximately 2.3–2.5 times
from that of the porous wall for both roughnesses tested in this study. We
observed that the additional roughness effect follows $Ck_{s_{r}}$, where
constant $C\approx 0.16$ for both rough wall test surfaces. Analysis of the
same data with the method suggested by Hutchins et al. (2023) shows that the
effect of porous and rough walls may be decoupled using an area-weighted
power-mean (with $n\approx$ 1.5 for both surfaces) to obtain an equivalent
roughness length for both rough wall test surfaces. The empirical formulation
as well as the physical interpretation of constants observed in this study
requires further data and should be the focus of future work.

[Funding]We gratefully acknowledge the financial support from EPSRC (Grant Ref
no: EP/S013296/1) and European Office for Airforce Research and Development
(Grant No: FA9550-19-1-7022, Programme Manager: Dr. Doug Smith).

[Declaration of interest]The authors declare no conflict of interest.

[Data availability statement]The data that support the findings of this study
will be made available upon publication.

[Author ORCID]D. D. Wangsawijaya: 0000-0002-7072-4245, P. Jaiswal:
0000-0002-5240-9911, B. Ganapathisubramani: 0000-0001-9817-0486

[Author contributions]PJ and DDW designed, carried out measurements, and post-
processed the data. DDW analysed the data and wrote the manuscript, BG was
responsible for conceptualisation, funding acquisition, editing of drafts and
project management.

## References

  * Breugem et al. (2006) Breugem, W. P., Boersma, B. J. & Uittenbogaard, R. E. 2006 The influence of wall permeability on turbulent channel flow. J. Fluid Mech. 562, 35–72.
  * Efstathiou & Luhar (2018) Efstathiou, Christoph & Luhar, Mitul 2018 Mean turbulence statistics in boundary layers over high-porosity foams. J. Fluid Mech. 841, 351–379.
  * Esteban et al. (2022) Esteban, L. B., Rodríguez-López, E., Ferreira, M. A. & Ganapathisubramani, B. 2022 Mean flow of turbulent boundary layers over porous substrates. Phys. Rev. Fluids 7 (9), 094603.
  * Flack & Schultz (2010) Flack, K. A. & Schultz, M. P. 2010 Review of hydraulic roughness scales in the fully rough regime. J. Fluids Eng. 132 (4), 041203.
  * Hahn et al. (2002) Hahn, Seonghyeon, Je, Jongdoo & Choi, Haecheon 2002 Direct numerical simulation of turbulent channel flow with permeable walls. J. Fluid Mech. 450, 259–285.
  * Hama (1954) Hama, F. R. 1954 Boundary-layer characteristics for rough and smooth surfaces. Trans. SNAME 62, 333–351.
  * Hutchins et al. (2023) Hutchins, N., Ganapathisubramani, B., Schultz, M. P. & Pullin, D. I. 2023 Defining an equivalent homogeneous roughness length for turbulent boundary layers developing over patchy or heterogeneous surface. Ocean Eng. 271, 113454.
  * Hutchins et al. (2009) Hutchins, N., Nickels, T. B., Marusic, I. & Chong, M. S. 2009 Hot-wire spatial resolution issues in wall-bounded turbulence. J. Fluid Mech. 635, 103–136.
  * Kong & Schetz (1982) Kong, F. & Schetz, J. 1982 Turbulent boundary layer over porous surfaces with different surface geometries. In the 20th AIAA Aerospace Sciences Meeting.
  * Ligrani & Bradshaw (1987) Ligrani, P. M. & Bradshaw, P. 1987 Subminiature hot-wire sensors: development and use. J. Phys. E: Sci. Instrum. 20, 323–332.
  * Manes et al. (2011) Manes, C., Poggi, D. & Ridolfi, L. 2011 Turbulent boundary layers over permeable walls: scaling and near-wall structure. J. Fluid Mech. 687, 141–170.
  * Manes et al. (2009) Manes, Costantino, Pokrajac, Dubravka, McEwan, Ian & Nikora, Vladimir 2009 Turbulence structure of open channel flows over permeable and impermeable beds: A comparative study. Phys. Fluids 21, 125109.
  * Monty et al. (2016) Monty, J. P., Dogan, E., Hanson, R., Scardino, A. J., Ganapathisubramani, B. & Hutchins, N. 2016 An assessment of the ship drag penalty arising from light calcareous tubeworm fouling. Biofouling 32 (4), 451–464.
  * Nikuradse (1933) Nikuradse, J. 1933 Strömungsgesetze in rauhen rohren. VDI Forschungsheft 361, English translation: 1950 Laws of flow in rough pipes. NACA Tech. Mem. 1292.
  * Rodríguez-López et al. (2015) Rodríguez-López, E., Bruce, P. J. K. & Buxton, O. R. H. 2015 A robust post-processing method to determine skin friction in turbulent boundary layers from the velocity profile. Exp. Fluids 56, 68.
  * Rosti et al. (2015) Rosti, Marco E., Cortelezzi, Luca & Quadrio, Maurizio 2015 Direct numerical simulation of turbulent channel flow over porous walls. J. Fluid Mech. 784, 396–442.
  * Shimizu et al. (1990) Shimizu, Y., Tsujimoto, T. & Nakagawa, H. 1990 Experiment and macroscopic modelling of flow in highly permeable porous medium under free-surface flow. J. Hydrosci. Hydraul. Eng. 8 (1), 69–78.
  * Suga et al. (2010) Suga, K., Matsumura, Y., Ashitaka, Y., Tominaga, S. & Kaneda, M. 2010 Effects of wall permeability on turbulence. Int. J. Heat Fluid Flow 31, 974–984.
  * Zagni & Smith (1976) Zagni, Anthony F. E. & Smith, Kenneth V. H. 1976 Channel flow over permeable beds of graded spheres. J. Hydraul. Div. 102 (2), 207–222.
  * Zippe & Graf (1983) Zippe, Hans J. & Graf, Walter H. 1983 Turbulent boundary-layer flow over permeable and non-permeable rough surfaces. J. Hydraul. Res. 21 (1), 51–65.