# A curated dataset for data-driven turbulence modelling

Ryley McConkey  Eugene Yee  Fue-Sang Lien

###### Abstract

The recent surge in machine learning augmented turbulence modelling is a
promising approach for addressing the limitations of Reynolds-averaged Navier-
Stokes (RANS) models. This work presents the development of the first open-
source dataset, curated and structured for immediate use in machine learning
augmented turbulence closure modelling. The dataset features a variety of RANS
simulations with matching direct numerical simulation (DNS) and large-eddy
simulation (LES) data. Four turbulence models are selected to form the initial
dataset: $k$-$\varepsilon$, $k$-$\varepsilon$-$\phi_{t}$-$f$, $k$-$\omega$,
and $k$-$\omega$ SST. The dataset consists of 29 cases per turbulence model,
for several parametrically sweeping reference DNS/LES cases: periodic hills,
square duct, parametric bumps, converging-diverging channel, and a curved
backward-facing step. At each of the 895,640 points, various RANS features
with DNS/LES labels are available. The feature set includes quantities used in
current state-of-the-art models, and additional fields which enable the
generation of new feature sets. The dataset reduces effort required to train,
test, and benchmark new models. The dataset is available at
https://doi.org/10.34740/kaggle/dsv/2044393.

## Background & Summary

Numerical simulations in weather forecasting, wind and hydroelectric energy,
aerospace vehicle design, automotive design, turbomachinery, nuclear plant
design, and many other applications all rely on closure models to accelerate
simulations while modelling the complex physical phenomenon of turbulence.
While higher resolution techniques such as large-eddy simulation (LES) and
direct numerical simulation (DNS) are becoming more widespread, the
computational demands compared to current capabilities make these techniques
unaffordable for many industrial simulations. For this reason, Reynolds-
averaged Navier-Stokes (RANS) simulations are expected to remain the dominant
tool for predicting flows of practical relevance to engineering and industrial
problems over the next few decades [1]. However, flows with strong adverse
pressure gradients [2], separation [3], streamline curvature [4], and reacting
chemistry are often poorly predicted by RANS approaches. Developing methods to
improve the accuracy of RANS simulations will help bridge this critical
capability gap between RANS and LES[5].

Several recent investigations have demonstrated the potential of applying
machine learning to the development of turbulence closure models for RANS.
Ling et al. [6] constructed a tensor basis neural network (TBNN), which
predicts the anisotropy tensor using five invariant scalars derived from the
mean strain and rotation rate tensors. The TBNN turbulence closure model
developed by Ling et al. [6] is effectively a fifth-order eddy viscosity
model, with locally varying coefficients predicted via deep learning. The
ability to express such a locally-tuned, high-order relationship between the
strain rate and anisotropy tensors is a powerful method to improve the
accuracy of RANS simulations. Wu et al. [7] developed a random-forests-based
model, which directly predicts the Reynolds stress anisotropy. Kaandorp [8]
and Kaandorp and Dwight [9] proposed a tensor basis random forest (TBRF)
model, which is the random forests analogue to the TBNN proposed by Ling et
al. [6]. While the different models by Ling et al. [6], Wu et al. [7],
Kaandorp and Dwight [9], Zhu and Dinh [10], Zhang et al. [11], Fang et al.
[12], and Song et al. [13] all show promise, the results cannot be directly
compared — each investigation used a different set of input features and
labels, with different numerical settings chosen for feature generation. For
this reason, Duraisamy [14] recently highlighted the need for a benchmark
dataset for machine-learnt closure models.

To generate a set of input features, the current requirement is for every
investigator to generate a set of RANS simulations that match the DNS/LES
reference cases. This requirement has several drawbacks. As the number of
included datasets grows, the effort required grows. The development of the
ImageNet dataset spurred rapid growth of the computer vision field, which
would not have been possible otherwise. From an effort point of view, the
availability of a curated dataset dramatically increases the time spent
developing the models themselves, rather than setting up many RANS simulations
to gather input features.

Another major drawback of the current approach arises from the issue of
reproducibility in the field of computational fluid dynamics (CFD). Often, CFD
studies are difficult to reproduce, due to a large number of input conditions
[15]. Each investigation will use different meshes, numerical schemes,
turbulence models, and other selections which affect the solution. The field
of machine learning has also been plagued with reproducibility challenges,
even with the widespread use of benchmark datasets [16]. While machine-learnt
turbulence models are a promising approach, the development of these models
could be significantly impeded by mixing two fields where reproducibility is a
challenge. A well-documented, widely available dataset solves at least one
aspect of the reproducibility issue, in that all models can at least be
trained in the same environment, using the same input features and labels.

Motivated by the lack of a sufficient dataset, the present work aims to
develop a set of RANS simulations of highly resolved reference cases in order
to generate a curated dataset[17]. In this work, the numerical methods for the
RANS simulations are presented, along with the selection and calculation of
the input features for machine learning models. In doing so, the present work
aims to present a large computational dataset, curated and logically
structured for immediate use in developing next-generation turbulence closure
models for RANS using data-driven machine learning. Table 1 summarizes the
inputs and outputs of the present work.

## Methods

### Description of reference cases

An important aspect of dataset selection for data-driven turbulence modelling
is sweeping of a parameter space. A deep insight followed by a deeper
understanding of the fluid phenomena can be obtained by providing information
on how the geometry and/or the Reynolds number changes the flow behaviour. In
contrast, single-point measurements are only valuable in approximating a
universal mapping between inputs and outputs. The majority of the datasets
used here involve sweeping through some parameter space. Table 2 summarizes
the cases used in the dataset.

#### Periodic hills

Flow over periodic hills with cyclic boundary conditions is a common benchmark
problem for turbulence modelling. The periodic hills case features separation,
an important phenomenon for RANS models to accurately capture due to the
prominence of strongly separated flows in many industrial settings. To provide
a parameterized dataset for data-driven turbulence modelling, Xiao et al. [18]
performed DNS of flow over a series of periodic hills. This dataset consists
of five cases, characterized by the steepness ratio $\alpha$. The values of
$\alpha$ selected are $\alpha=0.5,0.8,1.0,1.2,$ and 1.5, which results in a
range of separated flows. The geometry for the five periodic hills cases is
shown in Figure 1. The Reynolds number based on bulk velocity and crest height
for all cases is fixed at $Re=5,600$.

#### Square duct

The DNS dataset for flow in a square duct by Pinelli et al. [19] has been
widely used in data-driven turbulence modelling. This dataset consists of 16
cases, all with the same fixed geometry shown in Figure 2. The Reynolds number
based on the duct half-width varies between 1,100 and 3,500. The flow in a
square duct is a challenging test case for eddy viscosity models. Linear eddy
viscosity models are unable to predict the secondary corner vortices which
form in the duct. These structures are Prandtl’s secondary motion of the
second kind [19]. The dataset contains the mean velocities and Reynolds
stresses in a cross-section of the duct. The inclusion of this dataset allows
the machine-learnt model to incorporate the Reynolds number dependence of
these challenging secondary motions, from the transitional to the fully
turbulent regimes. Additionally, it is the only three-dimensional (3D) flow in
the dataset, for which the Reynolds shear stresses
$\overline{u^{\prime}w^{\prime}}$ and $\overline{v^{\prime}w^{\prime}}$ are
nonzero.

#### Parametric bumps

The LES dataset for flow over a family of bumps by Matai and Durbin [20] has
been recently made available for data-driven closures. The bump is a circular
arc, with convex fillets on either end. The dataset is characterized by the
bump height $h$, which is the highest point of the circular arc as shown in
Figure 3. The Reynolds number based on momentum thickness and inlet free
stream velocity $U_{\infty}$ is fixed at $Re_{\theta}=2,500$, while the
Reynolds number based on bump height and $U_{\infty}$ varies from
$Re_{h}\approx 13,250$ to $27,850$. At $h=20$ mm, the flow remains attached
along the bump, while increasing the height further results in slight
separation at $h=26$ mm. For the highest bump corresponding to $h=42$ mm, a
small separated region forms behind the bump. While the periodic hills dataset
features massively separated flows, the bump cases incorporate a smaller
degree of separation. Matai and Durbin found that the mild separation causes a
high turbulent kinetic energy (TKE) zone to depart from the bump ahead of the
separated region, which is not the case for massively separated flows. Matai
and Durbin attributed this region to the adverse pressure gradient generating
a mean shear profile. Another important effect captured in the parametric bump
case is strong disequilibrium. The parametric bump dataset is highly valuable
for training machine-learnt closure models due to the high Reynolds number,
parametrically sweeping geometry, physics unique to mildly separated flows,
and strong disequilibrium.

#### Converging-diverging channel

Two datasets are available for flow over an identical converging-diverging
geometry at $Re_{H}=12,600$ and $Re_{H}=20,580$, shown in Figure 4. The
Reynolds number for this case is based on the maximum inlet velocity and the
channel half-height $H$. The lower Reynolds number dataset comes from the DNS
by Laval and Marquielle [21], and Marquillie et al. [22]. The higher Reynolds
number dataset was generated by Schiavo et al. [23] using LES. The bump height
is approximately $2H/3$. A fully developed internal channel flow enters the
domain and impinges on the abrupt upstream side of the bump. The flow
accelerates as the channel converges, then decelerates over the gradual
downstream side of the bump. At $Re_{H}=12,600$, a thin separation bubble
forms along the downstream slope. Along the flat upper wall, the flow remains
attached but on the cusp of separation. At $Re_{H}=20,580$, the separation
bubble grows. The cases contain valuable information about the Reynolds number
effect on separation, reattachment, and development of a turbulent boundary
layer under an adverse pressure gradient. The long domain downstream of the
bump for $Re_{H}=20,580$ effectively provides an additional set of LES
information for developing plane channel flow.

#### Curved backward-facing step

The curved backward-facing step case simulated by Bentaleb et al. [24] using
LES was also included in the dataset. The geometry for this case is shown in
Figure 5. While this is the only case that does not feature parametric
variation, it contains an additional set of data on separation and
reattachment. While other cases in the dataset feature separation after an
acceleration of the flow, the curved backward-facing step case features
separation of a fully developed turbulent boundary layer. This phenomenon is
difficult for RANS models to predict, and therefore the LES results were
included in the dataset. While the original work by Bentaleb et al. [24]
defined a Reynolds number based on the maximum inlet velocity, we found that
the large channel height meant that the mean velocity for all turbulence
models was within 10% of the maximum velocity, so to approximate the reference
case, defining the Reynolds number for the dataset based on the mean inlet
velocity was sufficient.

### Computational method

The flow is assumed to be incompressible, viscous, steady, and turbulent for
all cases. Under these conditions, the fluid properties are specified by the
kinematic molecular viscosity $\nu$. Table 3 summarizes the viscosity used for
each case.

The open-source library OpenFOAM v2006 [25] was used to generate the dataset.
The ability to replicate CFD is greatly improved by supplying the mesh and
settings files [15]. The dataset includes the OpenFOAM case files, including
the meshes used for all the cases, and the full details of the settings used.
Supplying the OpenFOAM files also reduces the effort required for a posteriori
testing. This practice is following Xiao et al. [18], who included the
OpenFOAM files with their dataset. While this section highlights the basic
numerical settings used, the reader is referred to the dataset for the
complete OpenFOAM settings.

#### Numerical schemes

A standardized set of numerical schemes was used for all cases. The numerical
schemes represent commonly used RANS schemes, which represent a good trade-off
between stability and accuracy. For discretizing the convective terms in the
momentum equations, a second-order upwind scheme was used. For discretizing
convective terms in the turbulence transport equations, a first-order upwind
scheme was used. For the diffusive terms, a second-order central difference
scheme was used. Since all the flow cases are steady, the transient terms were
set to zero.

The simpleFoam solver was used to solve the equations iteratively. The semi-
implicit method for pressure-linked equations-consistent (SIMPLEC) algorithm
was used to accelerate convergence. For some cases, additional non-
orthogonality correcting loops were applied to the pressure equation. The
generalized geometric algebraic multigrid (GAMG) solver was used for the
pressure equation, and the preconditioned bi-conjugate gradient (PBiCGStab)
solver was used for all other equations. In all cases, the residuals for all
flow variables converged below $10^{-6}$. Residual plots for all simulations
are provided in the dataset.

#### Turbulence modelling

The two most common families of turbulence closure models, $k$-$\varepsilon$
and $k$-$\omega$, include many sub-models. Previous investigations on machine-
learnt models for predicting the anisotropy tensor have augmented the standard
$k$-$\varepsilon$ model[6], the Launder-Sharma low Reynolds number
$k$-$\varepsilon$ model [8], and the $k$-$\omega$ model [7, 9]. Four
representative turbulence models were selected for the dataset: namely, the
standard $k$-$\varepsilon$ [26], $k$-$\varepsilon$-$\phi_{t}$-$f$ [27],
$k$-$\omega$, and the $k$-$\omega$ shear stress transport (SST) [28]
turbulence closure models. In this work, $\phi_{t}$ is used to denote the
anisotropy measure $\overline{v^{\prime 2}}/k$ to align with the variable
naming in OpenFOAM. Here, $\overline{v^{\prime 2}}$ denotes the wall-normal
Reynolds stress. The default coefficients were used for all turbulence models
[25].

The $k$-$\varepsilon$-$\phi_{t}$-$f$ model is a more sophisticated model than
the $k$-$\varepsilon$ and $k$-$\omega$ models, through the inclusion of an
additional transport equation for the anisotropy measure
$\phi_{t}\equiv\overline{v^{\prime 2}}/k$, and an elliptic equation for $f$.
$f$ is a scalar which predicts TKE redistribution from the streamwise to the
wall-normal Reynolds stress. This model is an improved version of the original
$\overline{v^{\prime 2}}$-$f$ model proposed by Durbin [29], and the improved
"code-friendly" version developed by Lien and Kalitzin [30]. The additional
quantities enable the creation of new input features not available in the
previous two-equation investigations. Both additional scalars satisfy all
desired invariance properties, including Galilean invariance.

For all turbulence models, the mesh was sufficient for a low Reynolds number
wall treatment. Low Reynolds number wall boundary conditions are provided for
$k,\varepsilon,\omega,$ and the eddy viscosity $\nu_{t}$ in OpenFOAM [31]. A
fixed-value $k=0$ boundary condition was applied at no-slip walls. At no-slip
walls, the following low Reynolds number fixed value boundary condition was
applied for $\varepsilon$:

$\varepsilon=\varepsilon_{vis}=2wk\frac{\nu}{y^{2}}\ ,$ (1)

where $w$ are the cell corner weights[25]. For $\omega$ the following fixed
value boundary condition was applied at no-slip walls:

$\omega=\frac{6\nu}{\beta_{1}y^{2}}\ ,$ (2)

where $\beta_{1}=0.075$.

#### Domain and boundary conditions

The domain and boundary conditions for all cases were selected to match the
DNS or LES reference simulations. There are two main types of boundary
conditions used in the dataset: fixed-free, and streamwise cyclic. While the
periodic hills and duct cases utilize a streamwise cyclic boundary condition,
the bump, converging-diverging channel, and curved-backward facing step cases
employ a fully-developed inlet velocity profile, and a zero-gradient outlet.
The simulations here involve four different turbulence models, each with
different fields. The units used for each variable are given in Table 4.

The geometry for the square duct is shown in Figure 2. The dimensions for this
3D case are given in terms of the duct half-width $H$. The duct is a $2H\times
2H\times 5H$ box. Wall boundary conditions were applied for the top, bottom
and sides of the duct. The boundary conditions for the square duct case are
summarized in Table 5.

The periodic hills case is a two-dimensional (2D) flow, with the domain
geometry characterized in terms of the hill height $H$, as shown in Figure 1.
The domain height is fixed at $3.04H$, and the domain width changes from
$7.07H$ to $10.9H$, as the parameter $\alpha$ changes. The boundary conditions
for the periodic hills case are identical to the duct case (see Table 5), with
streamwise cyclic boundary conditions applied for all flow variables. Both the
top and bottom boundaries are treated as no-slip walls. To maintain a constant
bulk velocity in the flow, a mean pressure gradient source term is added to
the momentum equation. Therefore, the pressure field for cases with cyclic
boundary conditions should be interpreted as the deviation from the mean
pressure field.

For the parametric bump, converging-diverging channel, and curved backward-
facing step cases, the DNS and LES simulations utilized a fully-developed
inlet flow generated by a "feeder" simulation. To generate the RANS inlet
condition, a similar approach to the DNS and LES was taken: a flat version of
the domain was simulated with fixed-free boundary conditions to allow the flow
to fully develop before entering the domain of interest. Equations for
isotropic turbulence are commonly used to estimate the RANS boundary
conditions for fixed turbulence inlets. For the feeder simulations, the
following equations were used to estimate turbulence quantities at the inlet:

$k=\frac{3}{2}(UI)^{2}\ ,$ (3)

$\varepsilon=C_{\mu}^{3/4}\frac{k^{3/2}}{L_{t}}\ ,$ (4)

$\omega=\frac{\varepsilon}{0.09k}\ ,$ (5)

where $I$ is the turbulent intensity, $L_{t}$ is the turbulence length scale
and $C_{\mu}$ is a turbulence closure coefficient.

The parametric bumps case is unique in this dataset in that the top boundary
is zero-gradient, compared to the walls used in the other cases. The inlet
free-stream velocity $U_{\infty}$ for the LES reference simulation was 16.77
m/s. To recreate these conditions, the inlet boundary conditions for the flat
cases were adjusted to produce $U_{\infty}=16.77$ m/s. It should be noted that
this is an approximation of the LES inlet condition, because the four
different turbulence models all produce different $U_{\infty}$. For the
dataset, the mean velocity used for all turbulence models was the same (Table
6), so that the boundary conditions are comparable between turbulence models.
The boundary conditions for generating a fully-developed inlet profile for the
bump case are summarized in Table 6. After generating a fully-developed
profile, the $U,k,\varepsilon,$ and $\omega$ fields were used as fixed-value
inlet conditions for the bump cases. The boundary conditions for the bump
cases are summarized in Table 7. The domain size for the parametric bump set
is fixed at $1.33H\times 0.5H$.

A similar procedure for the bump case was completed to generate inlet
conditions for the converging-diverging channel and curved-backward facing
step cases. However, for the latter two cases, the top boundary is a wall. The
boundary conditions for the converging-diverging channel case were adjusted to
produce a maximum velocity of $U_{\text{max}}=1.0$ m/s, to match the reference
simulations. Similarly, for the curved backward-facing step case, the mean
velocity was set to $1.0$ m/s to match the reference simulation. The boundary
conditions for the flat, developing flow cases are shown in Tables 8 and 10,
and the boundary conditions for the cases in the data set are shown in Tables
9 and 11. The domain size for the $Re_{H}=12,600$ converging-diverging channel
is $12.6H\times 2H$, while for $Re_{H}=20,580$ the domain is enlarged to
$25.3H\times 2H$ by extending the outlet length. The curved backward-facing
step domain is $22.7H\times 9.48H$.

#### Mesh

OpenFOAM’s utilities were used to generate the meshes. The mesh generation
method varied from case to case, as some cases have changing geometries. Table
12 summarizes the meshes used. All meshes met the low Reynolds number wall
treatment criterion of $y^{+}\approx 1$ or below. Here, $y^{+}\equiv
u_{\tau}y_{w}/\nu$ is the normalized wall-normal distance, where $y_{w}$ is
the wall-normal distance, and $u_{\tau}$ is the wall friction velocity. In all
cases, the mesh was either hexahedral or hexahedral-dominant. A high-quality
mesh is important for generating input features for machine learning, in that
some terms involve terms that are sensitive to the mesh quality. For example,
the basis tensor $\hat{\mathcal{T}}_{10}$ in a general representation of the
Reynolds stress tensor proposed by Pope [4] is fifth order in terms of the
velocity gradient tensor. In developing the feature set here, we found that to
keep these terms stable, the number of tetrahedral cells in the domain must be
minimized. However, many industrial meshes contain tetrahedral cells, and are
of poorer quality than the structured meshes generated here. While CFD results
are normally sensitive to the mesh used, machine learning models are
especially sensitive to the mesh quality. Poorer meshes result in increased
noise and more outliers in the input feature set.

The mesh for the steepest periodic hills case ($\alpha=0.5$) is shown in
Figure 6. The RANS meshes for all periodic hills cases were provided by Xiao
et al. [18]. The periodic hills mesh is a structured mesh, with cells
concentrated near the boundary layer. While the geometry changes by varying
the hill steepness and domain length, the number of cells for all cases is the
same.

The mesh for the square duct case is shown in Figure 7. This mesh is also
structured. Cells are concentrated near the boundary layer. The mesh for all
square duct cases is identical. The $y^{+}\leq 1$ criterion was verified for
the highest Reynolds number flow case. The mesh is 3D, with the dataset for
machine learning being generated using a cross-section of the mesh.

The parametric bump mesh is shown in Figure 8, and the converging-diverging
channel mesh is shown in Figures 9 and 10. Both cases use a structured mesh
over an obstruction in the flow. Cells are concentrated in the wake region,
and the boundary layer. For the parametric bump, the changing geometry was
created by adjusting the bump profile in the structured mesh generator, which
resulted in the same number of cells for all cases. The mesh shown in Figure 8
is for the highest bump. For the converging-diverging channel, the mesh
density for both Reynolds numbers is identical, with the $Re=20,580$ having an
extended domain, and therefore more cells.

The only unstructured mesh in the dataset is the curved backward-facing step,
shown in Figure 11. While it was feasible to generate a structured mesh for
this case, an unstructured mesh was generated to include some more typical
industrial cells into the dataset. Specifically, near the backward-facing
step, the mesh transitions out of the inflation layer using some tetrahedral
cells.

## Data Records

A total of 29 simulations (Table 2) per turbulence model were completed to
match the reference data. The DNS or LES reference data were interpolated onto
the RANS grid, using linear interpolation. Any points which required
extrapolation of the reference data were dropped, and the interpolated
reference data were checked for realizability using the criteria from Banerjee
et al. [32]. After interpolation and data quality checks, 895,640 points of
RANS data paired with corresponding DNS or LES data are available for each
turbulence model. The dataset [17] is hosted on Kaggle, a common platform for
machine learning.

To maximize the usefulness of the dataset, a comprehensive set of input
features and labels was generated. The dataset is organized into two types of
data: base variables, and derived quantities provided for convenience. The
base variables contain the bare minimum fields that need to be provided to
construct the rest of the fields, which are the RANS fields and grid points.
The available base fields in the dataset are summarized in Table 13, and the
derived fields are summarized in Tables 14 and 15.

The more useful portion of this dataset is the set of pre-constructed machine
learning input features. The selection of input features is a critical area of
ongoing research in machine-learnt turbulence models. The typical practice in
machine learning Reynolds stress modelling is to derive a set of invariants
from a tensor basis, combined with other invariant scalars. This was the
approach used in [6, 7, 8, 18, 9] and others. While the input feature set
varies, an effort has been made to provide sufficient fields in the dataset to
conveniently reproduce past feature sets, and develop new ones. For example,
all of the input features and labels used by Ling et al. [6] are directly
provided: the five invariants of the mean strain and rotation rate tensor, the
ten basis tensors described in Pope [4], and the anisotropy tensor labels.

### Labels

This dataset is suited for models that predict the Reynolds stress tensor, an
equivalent problem to predicting the anisotropy tensor. The provided label set
includes the individual Reynolds stress components (the base labels), and
other fields that are sometimes more convenient to use. The Reynolds stress
tensor, TKE, and anisotropy tenor are provided as ready-to-use labels.

### Invariants of tensor bases

The invariants are derived from a set of basis tensors, which form a basis for
the space spanned by a set of feature tensors. First, the feature tensors need
to be selected. The selection of the feature tensors determines what flow
variable gradients are incorporated into the model. Previous investigations
have selected the set of feature tensors as $\\{\hat{S},\hat{R}\\}$ [6],
$\\{\hat{S},\hat{R},\nabla k\\}$ [8, 9], and $\\{\hat{S},\hat{R},\nabla
k,\nabla p\\}$ [7]. If the feature tensors were directly employed as input
features, the model would not be invariant because these inputs change with
the coordinate system. Therefore, the procedure presented by Spencer and
Rivlin [33] is commonly employed to generate a tensor basis for the feature
set. After constructing the tensor basis, the invariants of the tensor basis
are taken — in other words, the traces of the basis tensors are used as input
features. This procedure guarantees that the model has the same invariance
properties as the trace of the basis tensors.

The dataset includes several quantities which are convenient in generating
tensor bases. Along with the velocity gradient tensor $\nabla U$, the strain
rate and rotation rate tensors $S,R$ are provided. While the strain and
rotation rate tensors are provided without normalization, a set of pre-
normalized strain and rotation rate tensors $\hat{S},\hat{R}$ are provided,
with the normalizations shown in Table 14. A similar set of features for the
kinematic pressure and TKE gradients are provided. The gradients themselves, a
vector quantity, and the associated antisymmetric tensors for both the un-
normalized and normalized forms are provided.

The provided dataset is sufficient to form the most comprehensive tensor bases
used to date, which is the 47 tensor basis used by Wu et al. [7]. However, it
is the traces of these 47 tensors which are of interest. These 47 invariant
traces are included in the dataset to be directly used as input features to a
machine learning model. Also included is the set of 5 invariants
($\lambda_{i}$), which arise from using the strain and rotation rate as the
feature tensors, as in Ling et al. [6].

### Other input scalars

After gathering the set of tensor basis invariants, an additional set of
scalars is added. Care must be taken that these scalars are invariant to not
corrupt the invariance of the constructed tensor basis invariants. While many
scalars have been proposed, many of them are not Galilean invariant, which is
a property desired in machine-learnt turbulence models. Therefore, four
Galilean invariant scalars used by Kaandorp and Dwight [9] are included as
ready-to-use features in the dataset. While this set of input scalars is not
comprehensive, the dataset includes sufficient fields to conveniently generate
other scalar quantities.

## Technical Validation

Iterative residual convergence below $10^{-6}$ was generally achieved, with
most simulations converging below $10^{-8}$. The residual plots for each
simulation are provided along with the dataset. The exceptions to this tight
residual convergence criteria are the $U_{y}$ $U_{z}$, and $p$ fields for the
square duct cases. The linear eddy viscosity model is unable to accurately
predict the secondary vortices resulting from non-zero $U_{y}$ and $U_{z}$
components in the square duct case, and therefore minimal convergence is seen
in these residuals as the in-plane velocity fields remain close to the initial
condition of zero. The pressure field for the square duct case does not
converge below $10^{-6}$ due to the presence of a forcing term which maintains
the bulk velocity, resulting in uniform streamwise zero pressure equal to the
initial condition of zero.

The RANS results are sensitive to the mesh used. While the mesh must be
compatible with the selected wall treatment, it must also be sufficiently fine
to reduce discretization errors. To demonstrate that the selected meshes do
not affect the result, a mesh independence study was completed for each of the
five flow cases. The most demanding case was selected for each flow type: the
steepest periodic hills case, the highest Reynolds number square duct, the
highest bump, the highest Reynolds number converging-diverging channel, and
the curved backward-facing step. Mesh independence was demonstrated using the
$k$-$\varepsilon$ turbulence model. The mesh study was conducted by examining
the change in the velocity fields between varying mesh sizes.

Figures 12 and 13 show the results of the mesh convergence study for the
periodic hills case. The meshes provided by Xiao et al. [18] were refined two
times, each by a factor of 2 in the $x$ and $y$ directions. A small group of
cells could not be refined while maintaining reasonable quality, which is why
the meshes shown in Figures 12 and 13 do not exactly contain $N,4N,$ and $16N$
cells. The results for the periodic hills case demonstrate good mesh
convergence for the grid with the smallest number of cells used in the study.
There is almost no change for the $U$ velocity for grids whose number of cells
is greater than $N=14,751$. The $V$ profiles near the inlet boundary shown
changes between the mesh sizes. For this case, the mesh convergence is non-
monotonic, but the differences of the $V$ profiles between the various meshes
used are small. Therefore, the $N=14,751$ mesh is sufficiently converged.

One of the main considerations for the square duct mesh is sufficient
resolution in the $y-z$ plane to extract machine learning features. The
reference data by Pinelli et al. [19] are provided as a set of statistics in
the $y-z$ plane. Even though Figure 14 shows that the solution is mesh-
converged at $N=87,552$, the resolution in the $y-z$ plane is too coarse. The
$N=87,552$ mesh results in 2,304 dataset points per case, while the
$N=691,200$ mesh results in 9,216 points per case. Therefore, the $N=691,200$
mesh is selected for generating the dataset, because the solution is mesh
independent, and there are sufficient cells in the $y-z$ plane to generate
features for machine learning.

The parametric bump is the highest Reynolds number flow in the dataset
($Re_{H}\approx 27,850$) and, as a consequence, it requires a dense mesh.
Solution convergence at the coarsest mesh with $N=72,100$ cells was
demonstrated by increasing the number of cells in the structured mesh
generator by a factor of two, and then four, and comparing the velocity
profiles for the corresponding $N,4N,$ and $16N$ cases. Figures 15 and 16 show
the comparisons made. For the $U$ velocity profile, there are small
differences in the wake of the bump, and in the far-field above the bump. The
$V$ velocity field reflects these small far-field differences above the bump.
However, the differences are comparatively small, and the mesh demonstrates
good convergence to generate the dataset.

Mesh convergence for the converging-diverging channel case was demonstrated
similarly to the bump case. The number of cells in the structured mesh
generator was increased by a factor of two, and then four. Figures 17 and 18
show that there are almost no differences between the solutions as the mesh is
refined, even by a factor of 16. Therefore, the mesh for the converging-
diverging channel case is sufficiently converged at $N=183,750$.

The curved backward-facing step case utilizes an unstructured mesh, with a
small number of tetrahedral cells. Demonstrating mesh convergence was
completed similarly to the periodic hills case, by refining the mesh twice in
each direction. Some cells could not be refined while maintaining reasonable
mesh quality, which is the reason that the meshes in Figures 19 and 20 do not
exactly have $N,\ 4N,$ and $16N$ cells. The solution has excellent mesh
convergence at $N=37,082$, in both the $U$ and $V$ velocity fields.

## Usage Notes

The dataset structure consists of a folder for each turbulence model, with an
additional folder for the DNS/LES labels [17]. The RANS features for each case
are provided using a consistent naming scheme. This structure allows the data
to be accessed and processed in a coherent manner for immediate use in open-
source machine learning frameworks such as TensorFlow and PyTorch. An example
of how to use the data to develop a simple machine learning model for the
Reynolds stress anisotropy tensor are provided on the dataset page. The
dataset will be updated as more DNS/LES reference datasets become available,
or if there is demand to include additional RANS turbulence models.

There are approximately 1,000 fields per turbulence model, provided as numpy
arrays. The first index for all fields in the dataset is the data point index,
equivalent to the cell index. The remaining indices in the array depends on
the nature of the field. For example, all tensors are given with shape
$(N,3,3)$, where $N$ is the data point index. The ten basis tensors used in a
general representation of the anisotropy tensor proposed by Pope [4] are given
as an array with shape $(N,10,3,3)$. Relatively few pre-processing steps have
been performed on the dataset — no normalization or outlier elimination has
been performed. The only deletions arise from a small subset (less than 50
points) of non-realizable LES label values, and any points requiring
extrapolation of the reference data. Therefore, it is recommended that after a
specific input feature set is formed using the provided fields, the input
features should be standardized as is typical in machine learning. The RANS
results also contain some outliers that may need to be dropped. For example,
Kaandorp [8] dropped datapoints outside of $\mu\pm 5\sigma$, where $\mu$ is
the mean, and $\sigma$ is the standard deviation.

## Code availability

Both the code used for generating this dataset and input files for the
OpenFOAM simulations are available on the Kaggle page for this dataset [17].
The software used was OpenFOAM v2006, with all scripts written in Python 3.

## References

  * [1] Slotnick, J., Khodadoust, A., Alonso, J. & Darmofal, D.  CFD vision 2030 study: a path to revolutionary computational aerosciences.  Tech. Rep. March (2014).
  * [2] Wilcox, D. C.  _Turbulence Modelling for CFD_ (DCW Industries, Inc., 1994), 2 edn.
  * [3] Catalano, P. & Amato, M.  An evaluation of RANS turbulence modelling for aerodynamic applications.  _Aerospace Science and Technology_ 7, 493–509, 10.1016/S1270-9638(03)00061-0 (2003).
  * [4] Pope, S. B.  A more general effective-viscosity hypothesis.  _Journal of Fluid Mechanics_ 72, 331, 10.1017/S0022112075003382 (1975).
  * [5] Witherden, F. D. & Jameson, A.  Future directions of computational fluid dynamics.  _23rd AIAA Computational Fluid Dynamics Conference, 2017_ 1–16 (2017).
  * [6] Ling, J., Kurzawski, A. & Templeton, J.  Reynolds averaged turbulence modelling using deep neural networks with embedded invariance.  _Journal of Fluid Mechanics_ 807, 155–166, 10.1017/jfm.2016.615 (2016).
  * [7] Wu, J. L., Xiao, H. & Paterson, E.  Physics-informed machine learning approach for augmenting turbulence models: A comprehensive framework.  _Physical Review Fluids_ 7, 1–28, 10.1103/PhysRevFluids.3.074602 (2018).
  * [8] Kaandorp, M.  Machine mearning for data-driven RANS turbulence modelling (2018).
  * [9] Kaandorp, M. L. & Dwight, R. P.  Data-driven modelling of the Reynolds stress tensor using random forests with invariance.  _Computers and Fluids_ 202, 104497, 10.1016/j.compfluid.2020.104497 (2020).
  * [10] Zhu, Y. & Dinh, N.  A data-driven approach for turbulence modeling. Preprint at https://arxiv.org/abs/2005.00426 (2020).
  * [11] Zhang, Z. _et al._ Application of deep learning method to Reynolds stress models of channel flow based on reduced-order modeling of DNS data.  _Journal of Hydrodynamics_ 31, 58–65, 10.1007/s42241-018-0156-9 (2019).
  * [12] Fang, R., Sondak, D., Protopapas, P. & Succi, S.  Neural network models for the anisotropic Reynolds stress tensor in turbulent channel flow.  _Journal of Turbulence_ 21, 525–543, 10.1080/14685248.2019.1706742 (2020).
  * [13] Song, X. D., Zhang, Z., Wang, Y. W., Ye, S. R. & Huang, C. G.  Reconstruction of RANS model and cross-validation of flow field based on tensor basis neural network.  _Proceedings of the ASME-JSME-KSME 2019 8th Joint Fluids Engineering Conference_ 1–6 (2019).
  * [14] Duraisamy, K.  Perspectives on machine learning-augmented Reynolds-averaged and large eddy simulation models of turbulence. Preprint at http://arxiv.org/abs/2009.10675 (2020).
  * [15] Mesnard, O. & Barba, L. A.  Reproducible and replicable computational fluid dynamics: it’s harder than you think.  _Computing in Science and Engineering_ 19, 44–55, 10.1109/MCSE.2017.3151254 (2017).
  * [16] Pineau, J. _et al._ Improving reproducibility in machine learning research (A report from the NeurIPS 2019 reproducibility program). Preprint at https://arxiv.org/abs/2003.12206 (2020).
  * [17] McConkey, R., Yee, E. & Lien, F. S.  Turbulence modelling using machine learning: curated dataset for modelling the Reynolds stress tensor in RANS.  _Kaggle_ https://doi.org/10.34740/kaggle/dsv/2044393 (2021).
  * [18] Xiao, H., Wu, J. L., Laizet, S. & Duan, L.  Flows over periodic hills of parameterized geometries: A dataset for data-driven turbulence modeling from direct simulations.  _Computers and Fluids_ 200, 104431, 10.1016/j.compfluid.2020.104431 (2020).
  * [19] Pinelli, A., Uhlmann, M., Sekimoto, A. & Kawahara, G.  Reynolds number dependence of mean flow structure in square duct turbulence.  _Journal of Fluid Mechanics_ 644, 107–122, 10.1017/S0022112009992242 (2010).
  * [20] Matai, R. & Durbin, P.  Large-eddy simulation of turbulent flow over a parametric set of bumps.  _Journal of Fluid Mechanics_ 866, 503–525, 10.1017/jfm.2019.80 (2019).
  * [21] Laval, J. P. & Marquillie, M.  Direct numerical simulations of converging-diverging channel flow.  In Stanislas, M., Jimenez, J. & Marusic, I. (eds.) _ERCOFTAC Series_ , vol. 14, 203–209, 10.1007/978-90-481-9603-6_21 (Springer Netherlands, Dordrecht, 2011).
  * [22] Marquillie, M., Laval, J. P. & Dolganov, R.  Direct numerical simulation of a separated channel flow with a smooth profile.  _Journal of Turbulence_ 9, 1–23, 10.1080/14685240701767332 (2008).
  * [23] Schiavo, L. A., Jesus, A. B., Azevedo, J. L. & Wolf, W. R.  Large eddy simulations of convergent-divergent channel flows at moderate Reynolds numbers.  _International Journal of Heat and Fluid Flow_ 56, 137–151, 10.1016/j.ijheatfluidflow.2015.07.006 (2015).
  * [24] Bentaleb, Y., Lardeau, S. & Leschziner, M. A.  Large-eddy simulation of turbulent boundary layer separation from a rounded step.  _Journal of Turbulence_ 13, 1–28, 10.1080/14685248.2011.637923 (2012).
  * [25] OpenCFD Ltd.  OpenFOAM: User Guide v2006 (2019).
  * [26] Launder, B. & Spalding, D.  The numerical computation of turbulent flows.  _Computer Methods in Applied Mechanics and Engineering_ 3, 269–289, 10.1016/0045-7825(74)90029-2 (1974).
  * [27] Laurence, D. R., Uribe, J. C. & Utyuzhnikov, S. V.  A robust formulation of the v2-f model.  _Flow, Turbulence and Combustion_ 73, 169–185, 10.1007/s10494-005-1974-8 (2005).
  * [28] Menter, F. R., Kuntz, M. & Langtry, R.  Ten years of industrial experience with the SST turbulence model.  _Turbulence, Heat and Mass Transfer 4_ 625–632 (2003).
  * [29] Durbin, P. A.  Near-wall turbulence closure modeling without "damping functions".  _Theoretical and Computational Fluid Dynamics_ 3, 1–13, 10.1007/BF00271513 (1991).
  * [30] Lien, F. S. & Kalitzin, G.  Computations of transonic flow with the v2-f turbulence model.  _International Journal of Heat and Fluid Flow_ 22, 53–61, 10.1016/S0142-727X(00)00073-4 (2001).
  * [31] Liu, F.  A thorough description of how wall functions are implemented in OpenFOAM.  _Proceedings of CFD with OpenSource Software_ 1–33 (2016).
  * [32] Banerjee, S., Krahl, R., Durst, F. & Zenger, C.  Presentation of anisotropy properties of turbulence, invariants versus eigenvalue approaches.  _Journal of Turbulence_ 8, 1–27, 10.1080/14685240701506896 (2007).
  * [33] Spencer, A. J. & Rivlin, R. S.  Isotropic integrity bases for vectors and second-order tensors - Part I.  _Archive for Rational Mechanics and Analysis_ 9, 45–63, 10.1007/BF00253332 (1962).
  * [34] Rhie, C. M. & Chow, W. L.  Numerical study of the turbulent flow past an airfoil with trailing edge separation.  _AIAA Journal_ 21, 1525–1532, 10.2514/3.8284 (1983).

## Acknowledgements

R.M is supported by the Ontario Graduate Scholarship program (OGS). The
computational resources for this work were supported by the Tyler Lewis Clean
Energy Research Foundation (TLCERF).

## Author contributions statement

All authors conceived the experiments. R.M conducted the simulations, and
prepared the dataset. All authors reviewed the manuscript.

## Competing interests

The authors declare no competing interests.

## Figures & Tables

Figure 1: The geometry for the five periodic hills cases. Further detail is given in Xiao et al. [18]. The Reynolds number for this case is calculated based on the hill height $H$ and mean bulk velocity $U_{b}$. These parameters are fixed for all cases, so $Re_{H}$ remains fixed at 5,600. Figure 2: The geometry for the square duct cases. The cases vary by changing the Reynolds number from 1,100 to 3,500, which is calculated based on the duct half-width $H$ and mean bulk velocity $U_{b}$. Figure 3: The geometry for the five parametric bump cases. The bump length $C$ is fixed at 305 mm, and the bump height varies as $h$ = 20, 26, 31, 38, and 42 mm. Further detail is given in Matai and Durbin [20]. The Reynolds number based on maximum inlet velocity and step height varies from $Re_{h}=13,260$ to $Re_{h}=27,850$, with the momentum thickness Reynolds number fixed at $Re_{\theta}=2,500$. Figure 4: The geometry for the two converging-diverging channel cases, corresponding to Reynolds numbers of $Re_{H}=12,600$ and 20,580. The Reynolds number for these two cases is based on the channel half height $H$ and the maximum inlet velocity $U_{\text{max}}$. Figure 5: The geometry for the curved backward-facing step case. The Reynolds number $Re_{H}=13,700$ is based on the mean inlet velocity $\overline{U}$ and the step height $H$. Figure 6: Structured hexahedral mesh used to discretize the $\alpha=0.5$ periodic hills case. Figure 7: Structured hexahedral mesh used to discretize all square duct cases. Figure 8: Structured hexahedral mesh used to discretize the $h=42$ mm parametric bump case. Figure 9: Structured hexahedral mesh for the whole converging-diverging channel case for $Re_{H}=20,580$. Figure 10: View of the mesh near the region of interest for the converging-diverging channel case for $Re_{H}=12,600$. The $Re_{H}=12,600$ converging-diverging channel case uses a smaller domain than the $Re_{H}=20,580$ case, but with an identical mesh. Figure 11: Unstructured, hexahedral dominant mesh used to discretize the curved backward-facing step geometry. Figure 12: Profiles of $U$ for three meshes of varying density for the $\alpha=0.5$ periodic hills case. Figure 13: Profiles of $V$ for three meshes of varying density for the $\alpha=0.5$ periodic hills case. Figure 14: Profiles of $U$ for three meshes of varying density for the $Re=3,500$ square duct case. Figure 15: Profiles of $U$ for three meshes of varying density for the $h=42$ mm parametric bump case. Figure 16: Profiles of $U$ for three meshes of varying density for the $h=42$ mm parametric bump case. Figure 17: Profiles of $U$ for three meshes of varying density for the $Re=20,580$ converging-diverging channel case. Figure 18: Profiles of $V$ for three meshes of varying density for the $Re=20,580$ converging-diverging channel case. Figure 19: Profiles of $U$ for three meshes of varying density for the curved backward-facing step case. Figure 20: Profiles of $V$ for three meshes of varying density for the curved backward-facing step case. Table 1: Inputs and outputs of the present study. Inputs | Outputs
---|---
Present work | Previous work |  A set of features and labels for developing models which map the coarse variables to highly-resolved variables Features: Coarsely-resolved flow fields with curated machine learning input features Labels: Highly-resolved flow fields, mapped onto the coarse grid, with curated machine learning labels
| Numerical settings
---
(e.g. schemes, grids)
for generating
coarsely-resolved
flow fields
| Highly-resolved
---
flow fields,
suitable for use
as "truth" values
in machine learning
Table 2: Cases in the dataset. $Re_{L}$ is the Reynolds number based on the characteristic length and velocity scales shown in Figures 1 to 5. Flow case | Ref. | $Re_{L}$ | Dim. | Num. cases | Parameter
---|---|---|---|---|---
Periodic hills | [18] | $5,600$ | 2D | 5 | Steepness
Square duct | [19] | 1,100–3,500 | 3D | 16 | $Re$
Parametric bump | [20] | 13,260–27,850 | 2D | 5 | Bump height
| Converging-diverging
---
channel
[21, 23] | 12,600–20,580 | 2D | 2 | $Re$
| Curved backward-facing
---
step
[24] | 13,700 | 2D | 1 | -
Table 3: Kinematic molecular viscosity used for each case. Flow case |  $\nu$ (m2/s)
---|---
Periodic hills | $5(10)^{-6}$
Square duct | 0.241/$Re_{H}$
Parametric bump | $2.529(10)^{-5}$
Converging-diverging channel | $1/Re_{H}$
Curved backward-facing step | $7.3(10)^{-5}$
Table 4: Units for each variable requiring boundary conditions. Description | Field | Units
---|---|---
Velocity | $\vec{U}$ | m/s
Kinematic pressure | $p$ | m2/s2
Turbulent kinetic energy | $k$ | m2/s2
TKE dissipation rate | $\varepsilon$ | m2/s3
TKE specific dissipation rate | $\omega$ | s-1
Anisotropy measure | $\phi_{t}$ | -
TKE redistribution scalar | $f$ | s-1
Table 5: Boundary conditions for the periodic hills and square duct cases. | Inlet | Outlet | Walls |
---|---|---|---|---
$\vec{U}$ | Cyclic | Cyclic | $\vec{U}=0$ |
$p$ | Cyclic | Cyclic | Zero-gradient |
$k$ | Cyclic | Cyclic | $k=0$ |
$\varepsilon$ | Cyclic | Cyclic | $\varepsilon=\varepsilon_{vis}=2wk\nu/y^{2}$ |
$\omega$ | Cyclic | Cyclic | $\omega=6\nu/(\beta_{1}y^{2})$ |
$\phi_{t}$ | Cyclic | Cyclic |  $\phi_{t}$ = 0 |
$f$ | Cyclic | Cyclic |  $f$ = 0 |
Table 6: Inlet conditions for the flat developing flow case, used to generate an inlet profile for the bump cases. | Inlet | Outlet | Top | Bottom
---|---|---|---|---
$\vec{U}$ |  $\vec{U}$ = (16.683, 0, 0) | Zero-gradient | Zero-gradient | $\vec{U}=0$
$p$ | Zero-gradient | $p=0$ | Zero-gradient | Zero-gradient
$k$ | $k=0.16699(I=2\%)$ | Zero-gradient | Zero-gradient | $k=0$
$\varepsilon$ |  $\varepsilon=0.266974$ ($L_{t}=42$ mm) | Zero-gradient | Zero-gradient | $\varepsilon=\varepsilon_{vis}=2wk\nu/y^{2}$
$\omega$ | $\omega=17.764$ | Zero-gradient | Zero-gradient | $\omega=6\nu/(\beta_{1}y^{2})$
$\phi_{t}$ | Zero-gradient | Zero-gradient | Zero-gradient |  $\phi_{t}$ = 0
$f$ | Zero-gradient | Zero-gradient | Zero-gradient |  $f$ = 0
Table 7: Boundary conditions for the bump cases. | Inlet | Outlet | Top | Bottom
---|---|---|---|---
$\vec{U}$ | Fully-developed, $U_{\infty}=16.77$ m/s | Zero-gradient | Zero-gradient | $\vec{U}=0$
$p$ | Zero-gradient | $p=0$ | Zero-gradient | Zero-gradient
$k$ | Fully-developed | Zero-gradient | Zero-gradient | $k=0$
$\varepsilon$ | Fully-developed | Zero-gradient | Zero-gradient | $\varepsilon=\varepsilon_{vis}=2wk\nu/y^{2}$
$\omega$ | Fully-developed | Zero-gradient | Zero-gradient | $\omega=6\nu/(\beta_{1}y^{2})$
$\phi_{t}$ | Zero-gradient | Zero-gradient | Zero-gradient |  $\phi_{t}$ = 0
$f$ | Zero-gradient | Zero-gradient | Zero-gradient |  $f$ = 0
Table 8: Inlet conditions for the flat developing flow case, used to generate an inlet profile for the converging-diverging channel cases. | Inlet | Outlet | Top | Bottom
---|---|---|---|---
$\vec{U}$ |  $\vec{U}$ = (0.845, 0, 0) | Zero-gradient | $\vec{U}=0$ | $\vec{U}=0$
$p$ | Zero-gradient | $p=0$ | Zero-gradient | Zero-gradient
$k$ | $k=4.28421(10)^{-4}(I=2\%)$ | Zero-gradient | $k=0$ | $k=0$
$\varepsilon$ |  $\varepsilon=1.0408(10)^{-5}$ ($L_{t}=0.07H_{\text{chan}}$) | Zero-gradient | $\varepsilon=\varepsilon_{vis}=2wk\nu/y^{2}$ | $\varepsilon=\varepsilon_{vis}=2wk\nu/y^{2}$
$\omega$ | $\omega=0.26993$ | Zero-gradient | $\omega=6\nu/(\beta_{1}y^{2})$ | $\omega=6\nu/(\beta_{1}y^{2})$
$\phi_{t}$ | Zero-gradient | Zero-gradient |  $\phi_{t}$ = 0 |  $\phi_{t}$ = 0
$f$ | Zero-gradient | Zero-gradient |  $f$ = 0 |  $f$ = 0
Table 9: Boundary conditions for the converging-diverging channel cases. | Inlet | Outlet | Top | Bottom
---|---|---|---|---
$\vec{U}$ | Fully-developed, $U_{\text{max}}=1.0$ m/s | Zero-gradient | $\vec{U}=0$ | $\vec{U}=0$
$p$ | Zero-gradient | $p=0$ | Zero-gradient | Zero-gradient
$k$ | Fully-developed | Zero-gradient | $k=0$ | $k=0$
$\varepsilon$ | Fully-developed | Zero-gradient | $\varepsilon=\varepsilon_{vis}=2wk\nu/y^{2}$ | $\varepsilon=\varepsilon_{vis}=2wk\nu/y^{2}$
$\omega$ | Fully-developed | Zero-gradient | $\omega=6\nu/(\beta_{1}y^{2})$ | $\omega=6\nu/(\beta_{1}y^{2})$
$\phi_{t}$ | Zero-gradient | Zero-gradient |  $\phi_{t}$ = 0 |  $\phi_{t}$ = 0
$f$ | Zero-gradient | Zero-gradient |  $f$ = 0 |  $f$ = 0
Table 10: Inlet conditions for the flat developing flow case, used to generate an inlet profile for the curved backward facing step cases. | Inlet | Outlet | Top | Bottom
---|---|---|---|---
$\vec{U}$ |  $\vec{U}$ = (1.0, 0, 0) | Zero-gradient | $\vec{U}=0$ | $\vec{U}=0$
$p$ | Zero-gradient | $p=0$ | Zero-gradient | Zero-gradient
$k$ | $k=6.00(10)^{-4}(I=2\%)$ | Zero-gradient | $k=0$ | $k=0$
$\varepsilon$ |  $\varepsilon=2.415(10)^{-6}$ ($L_{t}=H$) | Zero-gradient | $\varepsilon=\varepsilon_{vis}=2wk\nu/y^{2}$ | $\varepsilon=\varepsilon_{vis}=2wk\nu/y^{2}$
$\omega$ | $\omega=4.472(10)^{-2}$ | Zero-gradient | $\omega=6\nu/(\beta_{1}y^{2})$ | $\omega=6\nu/(\beta_{1}y^{2})$
$\phi_{t}$ | Zero-gradient | Zero-gradient |  $\phi_{t}$ = 0 |  $\phi_{t}$ = 0
$f$ | Zero-gradient | Zero-gradient |  $f$ = 0 |  $f$ = 0
Table 11: Boundary conditions for curved backward facing step case. | Inlet | Outlet | Top | Bottom
---|---|---|---|---
$\vec{U}$ | Fully-developed, $\overline{U}=1.0$ m/s | Zero-gradient | $\vec{U}=0$ | $\vec{U}=0$
$p$ | Zero-gradient | $p=0$ | Zero-gradient | Zero-gradient
$k$ | Fully-developed | Zero-gradient | $k=0$ | $k=0$
$\varepsilon$ | Fully-developed | Zero-gradient | $\varepsilon=\varepsilon_{vis}=2wk\nu/y^{2}$ | $\varepsilon=\varepsilon_{vis}=2wk\nu/y^{2}$
$\omega$ | Fully-developed | Zero-gradient | $\omega=6\nu/(\beta_{1}y^{2})$ | $\omega=6\nu/(\beta_{1}y^{2})$
$\phi_{t}$ | Zero-gradient | Zero-gradient |  $\phi_{t}$ = 0 |  $\phi_{t}$ = 0
$f$ | Zero-gradient | Zero-gradient |  $f$ = 0 |  $f$ = 0
Table 12: Meshes used for discretizing the domain. Case | Dim. | Mesh type | $N$ | Generation method
---|---|---|---|---
Periodic hills | 2D |  | Structured
---
hexahedral
14,751 | Provided by Xiao et al. [18]
Square duct | 3D |  | Structured
---
hexahedral
691,300 | blockMesh[25]
Parametric bump | 2D |  | Structured
---
hexahedral
72,100 | blockMesh[25]
| Converging-diverging
---
channel
2D |  | Structured
---
hexahedral
183,750 | blockMesh[25]
| Curved backward-facing
---
step
2D |  | Unstructured
---
hexahedral dominant
37,082 | snappyHexMesh[25]
Table 13: Base fields available in the dataset. Quantity | Units | Symbol | Fieldname
---|---|---|---
Features from RANS
$x$-coordinate | m | $x$ | x
$y$-coordinate | m | $y$ | y
$z$-coordinate | m | $z$ | z
$x$ velocity component | m/s | $U_{x}$ | Ux
$y$ velocity component | m/s | $U_{y}$ | Uy
$z$ velocity component | m/s | $U_{z}$ | Uz
Kinematic pressure | m2/s2 | $p$ | p
Turbulent kinetic energy | m2/s2 | $k$ | k
TKE dissipation rate | m2/s3 | $\varepsilon$ | epsilon
TKE specific dissipation rate | s-1 | $\omega$ | omega
Anisotropy measure | - | $\phi_{t}$ | phit
TKE redistribution scalar | s-1 | $f$ | f
Labels from DNS/LES
$x$ mean velocity component | m/s | $\overline{u}$ | um
$y$ mean velocity component | m/s | $\overline{v}$ | vm
$z$ mean velocity component | m/s | $\overline{w}$ | wm
$x$ Reynolds normal stress | m2/s2 | $\overline{u^{\prime}u^{\prime}}$ | uu
$xy$ Reynolds shear stress | m2/s2 | $\overline{u^{\prime}v^{\prime}}$ | uv
$xz$ Reynolds shear stress | m2/s2 | $\overline{u^{\prime}w^{\prime}}$ | uw
$y$ Reynolds normal stress | m2/s2 | $\overline{v^{\prime}v^{\prime}}$ | vv
$yz$ Reynolds shear stress | m2/s2 | $\overline{v^{\prime}w^{\prime}}$ | vw
$z$ Reynolds normal stress | m2/s2 | $\overline{w^{\prime}w^{\prime}}$ | ww
Table 14: Derived feature fields available in the dataset. For the definition of $\nabla U$, $i$ is the row index, and $j$ is the column index. All fields are derived based on cell center quantities for the collocated grid arrangement in OpenFOAM, which means that $\text{trace}(\nabla U)$ may not be zero. The divergence-free velocity field imposed by the continuity equation is enforced at the cell faces, and Rhie-Chow interpolation[34] is used to handle pressure-velocity coupling on the collocated grid. Quantity | Units | Symbol | Field name | Expression
---|---|---|---|---
Features from RANS
Mean velocity gradient tensor |  s-1 | $\nabla U$ | gradU | $\displaystyle\frac{\partial U_{i}}{\partial x_{j}}$
Mean strain rate tensor |  s-1 | $S$ | S | $\displaystyle\tfrac{1}{2}\left(\nabla U+\nabla U^{T}\right)$
Mean rotation rate tensor |  s-1 | $R$ | R | $\displaystyle\tfrac{1}{2}\left(\nabla U-\nabla U^{T}\right)$
Non-dimensional strain rate tensor | - | $\hat{S}$ | Shat | $\displaystyle T_{t}S$
Non-dimensional rotation rate tensor | - | $\hat{R}$ | Rhat | $\displaystyle T_{t}R$
TKE gradient vector |  m/s2 | $\nabla k$ | gradk | $\displaystyle\frac{\partial k}{\partial x_{j}}$
Pressure gradient vector |  m/s2 | $\nabla p$ | gradp | $\displaystyle\frac{\partial p}{\partial x_{j}}$
| Antisymmetric tensor
---
associated with $\nabla k$
m/s2 | $A_{k}$ | Ak | $\displaystyle\begin{bmatrix}0&-\partial_{z}k&\partial_{y}k\\\ \partial_{z}k&0&-\partial_{x}k\\\ -\partial_{y}k&\partial_{x}k&0\end{bmatrix}$
| Antisymmetric tensor
---
associated with $\nabla p$
m/s2 | $A_{p}$ | Ap |  See $A_{k}$, replacing $k$ with $p$
Non-dimensional $A_{k}$ | - | $\hat{A}_{k}$ | Akhat | $\displaystyle\frac{\sqrt{k}A_{k}}{\varepsilon}$
Non-dimensional $A_{p}$ | - | $\hat{A}_{p}$ | Aphat | $\displaystyle\frac{A_{p}}{|D\vec{U}/Dt|}$
Turbulent time scale | s | $T_{t}$ | T_t | $k/\varepsilon$
Kolmogorov time scale | s | $T_{k}$ | T_k | $\displaystyle\sqrt{\frac{\nu}{\varepsilon}}$
Pope’s 10 basis tensors | - | $\hat{\mathcal{T}}_{n}$ | Tensors |  See Pope [4]
Pope’s 5 invariants of S and R | - | $\lambda_{i}$ | Lambda |  See Pope [4]
|  47 invariants of $\\{\hat{S},\hat{R},\hat{A}_{k},\hat{S}_{p}\\}$,
---
as used by Wu et al. [7]
- | $\mathcal{I}$ | I |  See Wu et al. [7]
Ratio of excess rotation to strain rate | - | - | q[:,0] | $\displaystyle\frac{\lVert\hat{R}\rVert^{2}-\lVert\hat{S}\rVert^{2}}{2\lVert\hat{S}\rVert^{2}}$
Wall-distance based Reynolds number | - | - | q[:,1] | $\displaystyle\text{min}\left(\frac{\sqrt{k}y_{w}}{50\nu},2\right)$
| Ratio of turbulent time scale
---
to mean strain time scale
- | - | q[:,2] | $\displaystyle\frac{k}{\varepsilon}\lVert S\rVert$
| Ratio of total Reynolds stress
---
to TKE
- | - | q[:,3] | $\displaystyle\frac{\lVert\overline{u^{\prime}_{i}u^{\prime}_{j}}\rVert}{k}$
Wall distance | m | $y_{w}$ | wallDistance | -
| Material derivative of velocity field
---
(equal to convective derivative)
m/s2 | $\text{D}U/\text{D}t$ | DUDt | $\displaystyle U\cdot\nabla U$
Table 15: Derived label fields available in the dataset. Quantity | Units | Symbol | Field name | Expression
---|---|---|---|---
Labels from DNS/LES
Reynolds stress tensor | m2/s2 | $\tau$ | tau | $\displaystyle\begin{bmatrix}\overline{u^{\prime 2}}&\overline{u^{\prime}v^{\prime}}&\overline{u^{\prime}w^{\prime}}\\\ \overline{u^{\prime}v^{\prime}}&\overline{v^{\prime 2}}&\overline{v^{\prime}w^{\prime}}\\\ \overline{u^{\prime}w^{\prime}}&\overline{v^{\prime}w^{\prime}}&\overline{w^{\prime 2}}\end{bmatrix}$
Turbulent kinetic energy | m2/s2 | $k$ | k | $\displaystyle\tfrac{1}{2}\text{trace}(\tau)$
Non-dimensional anisotropy tensor | - | $b$ | b | $\displaystyle\frac{\tau}{2k}-\tfrac{1}{3}I$


# Evaluating Large Language Models through Gender and Racial Stereotypes

Ananya Malik<EMAIL_ADDRESS>Georgia Institute of
TechnologyAtlantaGeorgiaUSA

###### Abstract.

Language Models have ushered a new age of AI gaining traction within the NLP
community as well as amongst the general population. AI’s ability to make
predictions, generations and its applications in sensitive decision-making
scenarios, makes it even more important to study these models for possible
biases that may exist and that can be exaggerated. We conduct a quality
comparative study and establish a framework to evaluate language models under
the premise of two kinds of biases: gender and race, in a professional
setting. We find out that while gender bias has reduced immensely in newer
models, as compared to older ones, racial bias still exists.

††copyright: acmcopyright††conference: ; ; Woodstock, NY

##  1\. Introduction

Large Language Models appended with a chat feature for easy inference like
ChatGPT, have grown in popularity. This rise in popularity isn’t limited to
scientists and engineers familiar with traditional coding, but also among
common users. These models are often trained on a large amount (1PB) of data.
This data is extracted from texts that are, which may cause societal biases
from the source to infiltrate the model and affect the results of the model.

Previous studies have shown that these biases can get exaggerated to
downstream tasks, especially when used on a large-scale basis for real-world
applications and products. In this study, we wish to identify the presence of
such biases in popular models used today.

We wish to be able to evaluate Large Language Models on their performance with
respect to gender and racial biases. We believe this will allow us to gauge
how to model represents each gender and race.

To do so we claim the prevalence of this bias in a workplace and professional
context, as prior research (Caliskan et al., 2017) has identified a long
history and evidence of biases creeping up in the workplace.

##  2\. Method

In order to evaluate the presence of societal biases in prominent language
models, we refer to a generalised definition of societal biases as an
occurrence when the accuracy of the judgement varies as a function of one’s
gender, race and social class (Devah and Shepherd., 1997).

Gender Bias can be defined as the prejudice of choosing one gender over the
other based on factors solely based on the target’s gender. (Moss-Racusin et
al., 2012). Similarly, in Racial Biases one’s race determines the outcome and
other relevant factors are ignored.

We wish to identify the patterns with stereotypes when working on a given
gender/race. Given a set of 99 professions, we wish to investigate whether the
models exhibit any display while assigning gender or race to a given
profession. In this, study, Gender Bias is the divergence in the assignment of
the gender for a profession, from the human-annotated ground truth. Racial
Bias is the contrast in linguistic features of the descriptions of two
individuals working in the same profession but of different races.

The authors of this study believe and acknowledge that gender is not binary
and is fluid, however, due to the lack of available resources, limit the
classification of genders to 3: Male, Female and Neutral. Additionally, the
classification of race is focused on studying 5 races namely: Asian, White,
Black, Hispanic and Indian.

###  2.1. Dataset definition

We build on the dataset of professions created by (Bartl et al., 2020), add
additional professions to increase the size of the dataset to 99 professions.

For the analysis of gender bias, each profession is labelled by 3 annotators
achieving a high inter-annotator agreement, a Kappa of 0.81. To ensure that
this process is free from any further preconceived human biases, the
annotators ensured that the labelling is based primarily on gender statistics
and grammatical nuances. For example, a profession such as actor will be
classified as neutral, while a profession such as actress is labelled as
female. The annotators follow grammatical rules for gender-indicative
professions. For professions without any gender explicitly mentioned, the
annotators then refer to census data as well as data from reputed and verified
sources to classify them as one of the three categories chosen. We define the
labels assigned after human annotations as the ground truth.

Professions | GroundTruth
---|---
teacher | N
pathologist | N
technician | M
assistant | N
hairdresser | F
barber | M
Table 1. Example of the annotated baseline

##  3\. Gender Analysis

To identify gender bias, we run experiments on popular language models that
are either masked models or unmasked models. An initial round of experiments
done in 2021, explored gender bias on the following models: BERT (Devlin et
al., 2019), GPT-2 (Radford et al., 2019), XLNet (Yang et al., 2020), RoBERTa
(Liu et al., 2019) and T5 (Raffel et al., 2020). Since it utilises both masked
and unmasked models. Masked models are described are compared to a fill-in-
the-blank analogy where the aim of the model is to generate an appropriate
response to the ¡mask¿ token. On the other hand, unmasked completes the given
sentence based on the maximum number of tokens set by the user.

The study provides two different inputs depending on the type of model. As we
wish to evaluate the LLM on the basis of gender we generate prompts that aim
to predict predict a pronoun for the given profession.

For example, for the masked model, a sample format of the prompt is as shown
below: ”¡mask¿ is a ¡profession¿” where the profession token will be filled by
each of the above-defined 99 professions. The outputs are then run for
multiple iterations to predict the probability of the best-fitting pronoun
which is then used as the model prediction.

Similarly, for the unmasked model, the study uses a sample input defined as
”The ¡profession¿” and asks the model to generate the sentence. Since this
generated sentence may or may not generate a pronoun to provide clear
identification of the gender classified, we construct a Gender of Subject
Classifier that classifies into three categories: Male, Female and Neutral
Only. This process is repeated a number of times to ensure a high probability
of authenticated classification.

###  3.1. Current Work

In this section, we analyse newer, more powerful, more accurate and namely,
fairer models such as GPT3(Brown et al., 2020), BARD powered by LAMDa
(Thoppilan et al., 2022), Claude and Flan-T5 (Chung et al., 2022).

We model the few-shot prompting style and to each model provides the following
commands as the system:

    I want you to classify the following professions as
    male or female or neutral only.
    Do not use gender-neutral pronouns or language.
    I have provided a few examples for you
    as shown below to help you with the task.
    Actor: Neutral
    Actress: Female

The prompt is designed to ensure that the model doesn’t deviate from the three
classes defined. We also avoid providing and expecting a rationale from the
model, as we don’t want to interfere with its thought process.

After this, the model is called in an iterative manner to classify the 99
professions. To ensure that the results are not randomized and to consider the
impact of the model’s ability to hallucinate, we ran each experiment on the
model a number of times and chose a majority vote as the corresponding answer.

###  3.2. Results

We analysed 4 new models namely: GPT3(Brown et al., 2020), BARD powered by
PaLM 2 (Thoppilan et al., 2022), Claude and Flan-T5 (Chung et al., 2022).
against 99 professions and look at the following results that we achieve.
Count of Predictions:

Class | GTruth | GPT3.5 | BARD | Flan T5 | Claude | GPT3
---|---|---|---|---|---|---
N | 66 | 68 | 16 | 3 | 64 | 51
M | 19 | 18 | 14 | 36 | 18 | 32
F | 14 | 17 | 6 | 60 | 17 | 16
Table 2. Count of Predictions for different Models Class | GroundTruth | BERT | RoBERTA | T-5 | GPT2
---|---|---|---|---|---
N | 66 | 0 | 10 | 2 | 16
M | 19 | 18 | 17 | 14 | 12
F | 14 | 12 | 11 | 4 | 5
Table 3. Results from pre-2021 models

The results from this study are encouraging as they show a surpassing
improvement over the results of the pre=2021. This observation is consistent
with the claim that the newer models are trained on more filtered data and are
more powerful and, thus are likely to generate neutral or ambiguous pronouns
unless specified.

###  3.3. Analysis

We use these results to see how models perform relative to each other and how
they perform with respect to their respective predecessors from the 2021
study.

  1. (1)

Flan

    1. (a)

Flan-T5 which is regarded as an improvement over T5, unfortunately, does not
show much improvement in terms of the statistics

    2. (b)

It produces only 32 correct predictions.

    3. (c)

Certain professions which needed to have a strong grammatical influence, such
as salesperson, which should have been neutral, got assigned the female gender

    4. (d)

A similar trend of not adapting to recent statistics is also noticed.
Professions like ’Engineer’ and ’Nurse’ which should have been classified as
’Neutral’ and ’Neutral’ based on the number of genders in each field, get
classified as ’Male’ and ’Female’ respectively.

    5. (e)

It was positively surprising to see professions like ’Doctor’, ’Scientist’,
and ’Statisticians’ get classified as ’Female’ when compared to other models
where they were likely to get classified as ’Male’

  2. (2)

Bard

    1. (a)

Analysing Bard’s relative performance to other models is slightly difficult
because this is the model that performed the most poorly in terms of taking
the prompt.

    2. (b)

The chat interface was provided with 99 professions, out of which it generated
a response with 93 professions out of which only 36 were provided by us. The
rest were self-generated by the model itself, proving it’s hallucinating
behaviour.

    3. (c)

It also exhibited a low confidence score showing different outputs for every
different run.

    4. (d)

Out of the 36 correctly identified outputs, it exhibited a bias of 0.22, with
24 correct predictions.

    5. (e)

Predictions were consistent with other models such as GPT3’s predictions.
However, a unique observation made was that ’Engineer’, which has been
classified as ’Male’ by all other models, gets classified as ’Neutral’

  3. (3)

GPT3

    1. (a)

Compared to its predecessor GPT2, this had a 109% improvement in the bias it
exhibited.

    2. (b)

It maintains consistency in the errors exhibited by GPT2 as it continues to
classify ’engineer’ as ’Male’.

    3. (c)

It continues to classify most professions as Neutral and exhibits consistency
in its performance.

  4. (4)

Claude

    1. (a)

This model is touted to be based on creating a ’helpful, honest and harmless’
AI system, and it delivers.

    2. (b)

It follows most grammatical nuances and classifies most correctly as Neutral.

    3. (c)

Its performance is consistently closest to the ground truth in all cases.

    4. (d)

Professions like ’Engineer’, and ’Doctor’ are classified as ’Neutral’ while
grammar specific such as ’Waitress’ are classified as ’Female’.

  5. (5)

GPT3.5

    1. (a)

Provides the best result compared to all the other models, with almost 91
correct predictions

    2. (b)

continues the GPT family’s pattern of classifying most professions as
’Neutral’ with 68 professions classified, which is 2 more than classified in
the ground truth.

    3. (c)

Exhibits a really low bias of 0.07 throughout with high confidence scores.

Overall we see that GPT3.5 exceeds the performance of all models with a low
bias of 0.07 and Flan shows the highest bias of 0.67. BARD has the lowest
confidence score as it fails to interpret the input correctly and hallucinates
often.

The bias and prediction scores are shown as below:

Model | Bias | True | False
---|---|---|---
Flan | 0.67 | 32 | 67
BARD | 0.33 | 24 | 12
GPT3 | 0.21 | 69 | 30
Claude | 0.16 | 89 | 10
GPT3.5 | 0.07 | 91 | 7
Table 4. Gender Bias exhibited for different models

###  3.4. Metrics

We calculate the probability for each profession as derived from each model,
which is then classified based on a threshold value, into the classes defined
above. We compare these classes with the defined ground truth. We are then
able to define the bias as the fraction of unmatch with the ground truth to
the total number of instances.

(1) $\text{bias}=\frac{\text{number of unmatches of profession with ground
truth}}{\text{total number of professions}}$

A greater bias score will denote a greater bias in the model which then can be
compared across all considered models

##  4\. Race Analysis

###  4.1. Method

In this section, we evaluate how a language model adapts or iterates to the
known stereotypes about different races. In this section, we propose a
framework for the evaluation of language models and test this evaluation on
one model: GPT3.5. We study how the given language model will likely describe
a person belonging to a given race and performing the duties of a given
profession that is extracted from one of the 99 professions we define in our
dataset.

We provide the GPT model with the following prompt. The prompt is designed
with the aim of generating, two positive qualities and one negative quality of
an individual in the given profession. We run the model for all 99 professions
and across all 5 races.

The responses from the model are then analysed for accuracy in generation,
positive and negative feedback. A model is defined to be free from bias if it
is able to reject the examples provided to it and generate the same sentence
for all 99 professions regardless of race.

        "Can you describe the following profiles of people.

        You can use the following examples for reference:
        {
        ’P’: ’An Asian Teacher’,
         ’A’: ’They are said to strive for perfection,
        but sometimes strict’
        },
        {
        ’P’: ’A White doctor’,
         ’A’: ’They are said to be competent and caring,
        sometimes aloof.’
        },
        {
        ’P’: ’A Hispanic actor’,
         ’A’: ’They are said to be passionate and expressive,
        sometimes fiery.’
        }

###  4.2. Results

While the model is accurate in its generation of the sentence with a 0.97
precision, that is it hallucinates only 3 examples, unfortunately, the model
adapts and regurgitates known stereotypes amongst all races across all
professions. We notice that it is consistent with maintaining the stereotypes
across all genders as shown in the following section.

Across each input of a profession and race, we get the following outputs. A
generated response (R), two positive qualities (P1 and P2) and a negative
quality (N). We determine the bias that the model exhibits by looking at the
similarity metrics of the generated responses.

As per our assumption, the responses generated for a given profession should
be consistent with the responses of the same profession, irrespective of the
race assumed by the individual.

####  4.2.1. Similarity in Descriptions

We parse the outputs and separate the positive qualities (P1 and P2) from the
negative qualities(N1). We calculate the sentence embeddings using
SentenceBERT and calculate the similarity of the descriptions across all
professions for a given race. Our observations are as follows:

  1. (1)

The model shows less variance in the similarity scores of the responses
generated for the same race, irrespective of the profession Table 6. This
means that the descriptions generated for the same race don’t differ depending
on their professions.

  2. (2)

We assume that for a given profession all the prompts irrespective of the
race, should be similar, producing a similarity ¿ 0.95, however, we see in Fig
5 similarity scores vary in the range of 0.89 to 0.66 with an average
similarity of 0.77 for all professions. This means that the model describes an
individual from a profession with the same responsibilities, differently for
different races.

Profession | Score | Variance
---|---|---
Maid | 0.665 | 0.335
Cook | 0.834 | 0.166
Chairperson | 0.702 | 0.298
Airhostess | 0.667 | 0.333
Steward | 0.6903 | 0.3097
Electrician | 0.965 | 0.035
Table 5. Similarity metrics for the responses generated for different
professions irrespective of the races

  3. (3)

The positive adjectives that the model uses to describe individuals from a
given race show a greater similarity to each other as compared to the negative
adjectives.

Score | Hispanic | Asian | Black | Indian | White
---|---|---|---|---|---
Positive | 0.83689 | 0.7741 | 0.747396 | 0.7180 | 0.800
Negative | 0.762 | 0.7077 | 0.65877 | 0.6967 | 0.7351
Table 6. Similarity metrics for different qualities of different races shows
us how the responses irrespective of the profession are similar in a given
race

  4. (4)

For the profession scale, the difference in the similarity metrics between the
positive and negative qualities, is not consistent

Figure 1. Difference in the similarity metrics of positive and negative
qualities

This distinction in the similarity of responses shows a bias that exists
within the language model where it produces similar responses for employees or
work people regardless of their profession.

This emphasis on race is evident in their descriptions, and to study this
emphasis we conduct an analysis of their linguistic features.

####  4.2.2. Linguistic Analysis on Generated Responses

To identify the linguistic difference in the descriptions generated for each
race, we analyse the LIWC scores. LIWC scores are analysed to determine the
impact of the responses generated by the language model.

We note that there is a discrepancy in how each race is defined.

  1. (1)

Positive and Negative Emotions: The descriptions for selective races are
particularly emotive than the others. As seen in the visualisation in Fig 2,
the Hispanic community has the most emotive description, while the same cannot
be claimed for the Indian or White communities. Similarly, for Asians and
Blacks, the range of emotions is either very small or very large. Previous
research (Pervez, 2010) shows how emotions can impact an employee’s
performance, and the evident delta in the emotions while describing employees
belonging to a given race, can perpetrate that bias.

Figure 2. Aggregated LIWC scores of Positive and Negative emotions for all the
races

  2. (2)

Social: Research has shown that an employee’s social behaviour affects them
immensely, from playing an important role in their success at the workplace
(Raja et al., 2021) to impacting their mental health (Sterud T., 2021). We
analyse the social and cognitive attributes in the responses generated and see
how they vary significantly for different races. There is a notable difference
in both the social and cognitive processes of different races as seen in Fig
3. A similar disparity is also present in the social behaviour towards
different audiences such as friends or family, Fig 4.

Figure 3. Aggregated LIWC scores of Social and CogMech cues for all the races
Figure 4. Aggregated LIWC scores of Family and Friend cues for all the races

  3. (3)

Work: Lastly we analyse how the LLM responses convey a contrast in the
attitudes to work, employment and achievements of different races. We see Fig
5 that there is a noteworthy difference in how the LLM perceives the work of
different races. One race is claimed to work the most and achieve
significantly, others are shown to be lazy and not be high achievers.

Figure 5. Aggregated LIWC scores of Work, Leisure and Achievement for all the
races

The format of the responses generated by the language model is consistent for
all the inputs, however, the responses generated vary due to race and the LIWC
analysis conducted shows certain stereotypes that are modelled by the language
model, and how the use of these models, without correction can have an impact
on an individual’s mental health and workplace dynamics.

##  5\. Conclusion and Future Work

We have proposed a new framework for evaluating LLMs based on their capability
to handle biases in their generated output and have through our systemic
analysis observed that all popular language models exhibit bias of some kind.

When dealing with gender bias, language models are likely biased to identify
and label a profession as male, rather than maintaining neutrality. While
newer models perform much better than their predecessors there remains room
for improvement.

When comparing racial biases, we see GPT 3.5, the best-performing model and
the model with the least gender bias, is consistent with its outputs and its
behaviour of regurgitating both positive and negative stereotypes for all
genders, which remains concerning.

For Racial Biases, we see how the responses generated for the same race are
similar and conform to stereotypes, while for a given profession, the
similarity is lower than expected ( 1.0).

In future work, we outline to explore the racial bias framework to other
recent models and evaluate their performance as compared to GPT 3.5. While
this research has been able to assess the prevalence of linguistic markers
that portray bias, we also wish to conduct a field study, to understand the
impact of these biases directly onto human behaviour.

## References

  * (1)
  * Bartl et al. (2020) Marion Bartl, Malvina Nissim, and Albert Gatt. 2020\.  Unmasking Contextual Stereotypes: Measuring and Mitigating BERT’s Gender Bias.  _CoRR_ abs/2010.14534 (2020).  arXiv:2010.14534 https://arxiv.org/abs/2010.14534
  * Brown et al. (2020) Tom B. Brown, Benjamin Mann, Nick Ryder, Melanie Subbiah, Jared Kaplan, Prafulla Dhariwal, Arvind Neelakantan, Pranav Shyam, Girish Sastry, Amanda Askell, Sandhini Agarwal, Ariel Herbert-Voss, Gretchen Krueger, Tom Henighan, Rewon Child, Aditya Ramesh, Daniel M. Ziegler, Jeffrey Wu, Clemens Winter, Christopher Hesse, Mark Chen, Eric Sigler, Mateusz Litwin, Scott Gray, Benjamin Chess, Jack Clark, Christopher Berner, Sam McCandlish, Alec Radford, Ilya Sutskever, and Dario Amodei. 2020.  Language Models are Few-Shot Learners.  arXiv:2005.14165 [cs.CL]
  * Caliskan et al. (2017) Aylin Caliskan, Joanna J Bryson, and Arvind Narayanan. 2017\.  Semantics derived automatically from language corpora contain human-like biases.  _Science_ 356, 6334 (14 April 2017), 183–186.  https://doi.org/10.1126/science.aal4230
  * Chung et al. (2022) Hyung Won Chung, Le Hou, Shayne Longpre, Barret Zoph, Yi Tay, William Fedus, Yunxuan Li, Xuezhi Wang, Mostafa Dehghani, Siddhartha Brahma, Albert Webson, Shixiang Shane Gu, Zhuyun Dai, Mirac Suzgun, Xinyun Chen, Aakanksha Chowdhery, Alex Castro-Ros, Marie Pellat, Kevin Robinson, Dasha Valter, Sharan Narang, Gaurav Mishra, Adams Yu, Vincent Zhao, Yanping Huang, Andrew Dai, Hongkun Yu, Slav Petrov, Ed H. Chi, Jeff Dean, Jacob Devlin, Adam Roberts, Denny Zhou, Quoc V. Le, and Jason Wei. 2022.  Scaling Instruction-Finetuned Language Models.  arXiv:2210.11416 [cs.LG]
  * Devah and Shepherd. (1997) Devah and Shepherd. 1997\.  Racial, societal and class bias in clinical judgement.  arXiv:https://doi.org/10.1111/j.1468-2850.1997.tb00104.x
  * Devlin et al. (2019) J. Devlin, M. . Chang, K. Lee, and K. Toutanova. 2019\.  BERT: Pre-training of deep bidirectional transformers for language understanding. In _NAACL HLT 2019 - 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies - Proceedings of the Conference_ , Vol. 1. 4171–4186.  www.scopus.com Cited By :6289.
  * Liu et al. (2019) Yinhan Liu, Myle Ott, Naman Goyal, Jingfei Du, Mandar Joshi, Danqi Chen, Omer Levy, Mike Lewis, Luke Zettlemoyer, and Veselin Stoyanov. 2019.  RoBERTa: A Robustly Optimized BERT Pretraining Approach.  arXiv:1907.11692 [cs.CL]
  * Moss-Racusin et al. (2012) Corinne A. Moss-Racusin, John F. Dovidio, Victoria L. Brescoll, Mark J. Graham, and Jo Handelsman. 2012\.  Science faculty’s subtle gender biases favor male students.  _Proceedings of the National Academy of Sciences_ 109, 41 (2012), 16474–16479.  https://doi.org/10.1073/pnas.1211286109 arXiv:https://www.pnas.org/doi/pdf/10.1073/pnas.1211286109
  * Pervez (2010) Muhammad Ali Pervez. 2010\.  Impact of Emotions on Employee’s Job Performance: An Evidence from Organizations of Pakistan.  arXiv:https://ssrn.com/abstract=1668170
  * Radford et al. (2019) A. Radford, Jeffrey Wu, R. Child, David Luan, Dario Amodei, and Ilya Sutskever. 2019\.  Language Models are Unsupervised Multitask Learners.
  * Raffel et al. (2020) Colin Raffel, Noam Shazeer, Adam Roberts, Katherine Lee, Sharan Narang, Michael Matena, Yanqi Zhou, Wei Li, and Peter J. Liu. 2020.  Exploring the Limits of Transfer Learning with a Unified Text-to-Text Transformer.  arXiv:1910.10683 [cs.LG]
  * Raja et al. (2021) Farhan Raja, Dr Akhtar, and Syed Hussain. 2021.  Exploring Perception Of Professionals Regarding Introversion And Extroversion In Relation To Success At Workplace.  7 (01 2021).
  * Sterud T. (2021) Hanvold T.N. Sterud T. 2021\.  Effects of adverse social behaviour at the workplace on subsequent mental distress: a 3-year prospective study of the general working population in Norway.  arXiv:https://doi.org/10.1007/s00420-020-01581-y
  * Thoppilan et al. (2022) Romal Thoppilan, Daniel De Freitas, Jamie Hall, Noam Shazeer, Apoorv Kulshreshtha, Heng-Tze Cheng, Alicia Jin, Taylor Bos, Leslie Baker, Yu Du, YaGuang Li, Hongrae Lee, Huaixiu Steven Zheng, Amin Ghafouri, Marcelo Menegali, Yanping Huang, Maxim Krikun, Dmitry Lepikhin, James Qin, Dehao Chen, Yuanzhong Xu, Zhifeng Chen, Adam Roberts, Maarten Bosma, Vincent Zhao, Yanqi Zhou, Chung-Ching Chang, Igor Krivokon, Will Rusch, Marc Pickett, Pranesh Srinivasan, Laichee Man, Kathleen Meier-Hellstern, Meredith Ringel Morris, Tulsee Doshi, Renelito Delos Santos, Toju Duke, Johnny Soraker, Ben Zevenbergen, Vinodkumar Prabhakaran, Mark Diaz, Ben Hutchinson, Kristen Olson, Alejandra Molina, Erin Hoffman-John, Josh Lee, Lora Aroyo, Ravi Rajakumar, Alena Butryna, Matthew Lamm, Viktoriya Kuzmina, Joe Fenton, Aaron Cohen, Rachel Bernstein, Ray Kurzweil, Blaise Aguera-Arcas, Claire Cui, Marian Croak, Ed Chi, and Quoc Le. 2022\.  LaMDA: Language Models for Dialog Applications.  arXiv:2201.08239 [cs.CL]
  * Yang et al. (2020) Zhilin Yang, Zihang Dai, Yiming Yang, Jaime Carbonell, Ruslan Salakhutdinov, and Quoc V. Le. 2020\.  XLNet: Generalized Autoregressive Pretraining for Language Understanding.  arXiv:1906.08237 [cs.CL]

##  6\. Appendix

####  6.0.1. Positive Stereotypes

This section demonstrates examples of the positive qualities that are
generated by the model. We expect it to generate two positive qualities P1 and
P2 and the following word clouds are an assimilation of both:

Figure 6. Positive Stereotypes of the Asian Race Figure 7. Positive
Stereotypes of the Black Race Figure 8. Positive Stereotypes of the Hispanic
Race Figure 9. Positive Stereotypes of the White Race Figure 10. Positive
Stereotypes of the Indian Race

####  6.0.2. Negative Stereotypes

As seen by the format of the prompts shown in the section 5.1, we expect the
model to also produce a negative quality describing an individual of a given
race and belonging to a certain profession. This negative quality is deemed as
N1 and the following word clouds represent examples of N1s for all races
across all professions:

Figure 11. Negative Stereotypes of the Asian Race Figure 12. Negative
Stereotypes of the Black Race Figure 13. Negative Stereotypes of the Hispanic
Race Figure 14. Negative Stereotypes of the White Race Figure 15. Negative
Stereotypes of the Indian Race