11institutetext: NORCE Norwegian Research Centre, Bergen, Norway
22institutetext: School of Energy, Geoscience, Infrastructure and Society,
Heriot-Watt University, Edinburgh, United Kingdom

# Deep learning for prediction of complex geology ahead of drilling

Kristian Fossum 11 Sergey Alyaev 11 Jan Tveranger 11 Ahmed Elsheikh 22

###### Abstract

During a geosteering operation the well path is intentionally adjusted in
response to the new data acquired while drilling. To achieve consistent high-
quality decisions, especially when drilling in complex environments, decision
support systems can help cope with high volumes of data and interpretation
complexities. They can assimilate the real-time measurements into a
probabilistic earth model and use the updated model for decision
recommendations.

Recently, machine learning (ML) techniques have enabled a wide range of
methods that redistribute computational cost from on-line to off-line
calculations. In this paper, we introduce two ML techniques into the
geosteering decision support framework. Firstly, a complex earth model
representation is generated using a Generative Adversarial Network (GAN).
Secondly, a commercial extra-deep electromagnetic simulator is represented
using a Forward Deep Neural Network (FDNN).

The numerical experiments demonstrate that the combination of the GAN and the
FDNN in an ensemble randomized maximum likelihood data assimilation scheme
provides real-time estimates of complex geological uncertainty. This yields
reduction of geological uncertainty ahead of the drill-bit from the
measurements gathered behind and around the well bore.

###### Keywords:

Geosteering Machine Learning Deep Neural Network Generative Adversarial
Network Ensemble randomized maximum likelihood

##  1 Introduction

The process of drilling wells for hydrocarbon production represents a major
cost in petroleum reservoir development. However, drilling of new wells is
necessary to increase the total oil recovery. To maximize the value for each
drilled well it is necessary to optimize the placement of the well within the
reservoir structure. An optimally placed well will mobilize more of the
petroleum resources, and reduce the need for injected water  – reducing the
environmental impact of oil production.

To place a well in its optimal position, operators apply geosteering. Here,
the well trajectory is adjusted while drilling in response to real-time
measurement of the geology surrounding the drill bit. The value of geosteering
has been well documented in the literature [1, 17, 19].

The main objective with geosteering is to utilize the information in the
measurements to make optimal decisions. Hence, geosteering can be seen as a
sequential decision process under uncertainty and should be treated in a
probabilistic framework [20]. Recently, a workflow based on the Ensemble
Kalman Filter (EnKF) [12] has been employed to condition the geological model
on measurements acquired while drilling [9, 22]. In the EnKF the uncertainty
is represented by an ensemble of equiprobable realizations. This workflow has
then been combined with a global optimization method and applied as a Decision
Support System (DSS)[3].

The DSS framework provides high quality decisions on synthetic cases, but
practical challenges should be addressed for it to be applicable to real
operations [3]. Firstly, to our knowledge, there is no studies which combine
the ensemble update with a commercial tool for simulation of measurements.
Secondly, the earth model utilized in the published studies does not represent
a realistic geological complexity. Conceptually, it is easy to insert any
numerical model for simulating the measurements into the DSS workflow.
Similarly, there is nothing that prohibits the use of complex earth models.
However, as the complexity increase, the numerical run-time also increases
hindering real-time performance. Moreover, complex earth models can not
typically be represented using a Gaussian distribution, and consequently EnKF
updates will not retain the geological complexity [28].

Figure 1: The proposed DSS workflow. Green boxes highlight the new elements
introduced in this paper.

In this paper, we introduce important elements to make the DSS better suited
to real operations, see Figure 1. The main novelty in our approach is to
introduce a machine-learning method to represent both the earth model and the
forward model of extra-deep borehole electromagnetic (EM) measurements. Within
this setting we demonstrate that a real-time ensemble-based inversion can
indeed predict the distribution of non-trivial geology ahead of bit.

To construct a reference earth model we generate realizations of a fluvial
geological environment using a commercial software. These realizations are
then sub-sampled to form a training dataset for the offline training of a
Generative Adversarial Network (GAN). The GAN is then used, online, to
generate plausible geological realizations from a low-dimensional Gaussian
input vector. The earth modeling is described in Section 2. For modeling the
extra-deep EM measurements we use a forward deep neural network (FDNN) trained
on a dataset generated using a commercial simulator (Section 3). In Section 4
we discuss the exact and the fast data assimilation (DA) methods. The
numerical results, showing the applicability of our proposed method are given
in Section 5. Finally, we summarize and conclude the paper in Section 6.

##  2 Earth modeling using GAN

GANs are a class of unsupervised machine learning methods which can learn to
generate new formatted data with the same statistics as the training set.
Motivated by successful applications of GANs for modeling channelized
structures for reservoir simulation [8, 7], we use a GAN for efficient earth
modeling.

The GAN consists of two deep neural networks (DNNs): a generator and a
discriminator. The generator takes a random Gaussian low dimensional vector as
input and generates a realization of formatted data: geological realization.
The discriminator takes the formatted data and gives a probability of it being
’real’, i.e., belonging to the training set. During training the DNNs contest
each other in a min-max game. They are trained simultaneously. On each
training step the generator creates (fake) geological realization from random
vectors. Fake geological realizations are combined with random samples of real
earth model and are fed to the discriminator. The loss function for the
generator is proportional to number of ’fakes’ correctly identified by the
discriminator. The loss function for the discriminator is proportional to the
total misjudged data samples. In our study we use an adapted Wasserstein GAN
with hierarchical convolutional networks for the generator and the
discriminator, see [4] for implementation details.

Figure 2: The original earth model generated by the commercial tool.

For geosteering we want to reproduce likely geological realizations of facies
and porosity distributions on a 2D vertical geological section along the well
to identify the oil-bearing sands ahead of bit. For training of GAN we use a
large (compare to the area of prediction) reference earth model, which should
provide a realistic test case for the present study in terms of scale and
actual geological features and properties. The reference earth model is
constructed using a commercial software that models a synthetic structural
framework, a facies model set-up derived from outcrop analogue data, and
synthetic petrophysical properties of individual facies derived from published
literature. The resulting model measures 4000m x 1000m x 200m (xyz) with cell
dimensions set to 10x10x0.5 m, yielding a regular corner-point grid of size
400x100x400, see Figure 2.

The constructed facied model represents a low net/gross fluvial depositional
system. It was chosen since it provides complex 3D architectures comprising a
limited number of facies, which form contrasting geometries, see Figure 2.
Input numbers for statistical generation of facies and geometries are derived
from a well-documented outcrop of the Cretaceous lower Williams Fork Formation
(Mesa Verde Group) at Coal Canyon, Colorado, USA [26, 27, 30].

Key parameters of the facies model set-up are listed in Table 1. The model is
not intended as a rendering of the outcrop itself and is consequently
simplified compared to descriptions of the original outcrop [11, 24, 25, 27].
The model contains three facies: Background/shale, Channels and Crevasse
splays. The probability distribution of channel width in the model is adapted
to include “narrow channel bodies”, and stacking of channels accounts for
multi-story channels which comprise more than 80% of the observed channel
bodies. The flow direction of the channel system is set towards $45\pm 10$
degrees. No trends were used to condition the spatial distribution of
channels.

Table 1: Parameter settings for facies models.

Volumetric fraction | Value | Tolerance | Comments
---|---|---|---
Channel system volume fraction | 0.3 | 0.05 |
Channel positioning | 1 |  | No trends
Crevasse volume fraction | 0.1 | 0.03 | Of channel system vol. frac.

Channel geometry | Value | SD | Min. | Max.
---|---|---|---|---
Thickness | 4.2 | 1.5 |  |
Width | 155 | 50 | 20 | 500
Correlation W/T | 36 |  |  |
Amplitude | 400 | 50 |  |
Sinuosity | 1.3 |  |  |
Azimuth | 45 | 10 |  |

Form/repulsion | Setting
---|---
Cross-section geometry | Parabolic, basic variability
Channel form | Rigid
Repulsion | None

The geological realization is parameterized by a vector of 60 independent
parameters. For each 60-dimensional vector, the generator outputs a 64x64 grid
with three values in each grid block. For a grid block (with dimensions 10.0m
along-well and 0.5m thickness) the three values, ’channels’, represent the
probability of the grid-block belonging to the respective facies class:
Background/Channel/Crevasse. Our generator is also predicting
porosity/resistivity distribution within the geo-bodies, but in this initial
study only the facies classes are used.

For training, the original 3D earth model is sampled as 64x64 2D images with
three channels. The facies index from the training set is converted into one-
hot three-dimensional vector. That is, the vector represents the probability
of facies: the value of the true index is set to one and other channels to
zero. During evaluation, the resistivity of the facies with highest
probability is applied.

##  3 LWD Neural Network

To maintain real-time performance of a data assimilation workflow the forward
model should be fast and support batch, preferably parallel execution.
Proprietary forward models provided by measurement instrument vendors provide
the most accurate results, but they are often not sufficiently fast, and not
always optimized for batch execution. In [2], the authors developed a DNN
approximation of such a forward model [29], which we abbreviate FDNN.

The model approximates the output of the ultra-deep electromagnetic well-bore
logging instrument. The instrument is configured to transmit four shallow and
nine pairs of deep directional measurements, and has sensitivity to boundaries
up to 30 meters to the side from the well bore. We emphasize that the tool
provides information around, but not ahead of the drilling position. An
illustration of the deep measurements is provided in Figure 3.

The input to the FDNN model is a layered geological media with up to three
boundaries above and below the measurement instrument as well as the
resistivity values of all seven layers. In this study we assume that the layer
resistivity is isotropic and that the well is aligned with the horizontal
axis.

We produce one synthetic set of measurements for every horizontal position of
the gridded model which we ’drill’ through. First, we choose most probable
facies for each ’pixel’ and substitute it with the corresponding resistivity
value:

  1. 1.

Background, $R=220.0$ Ohm m;

  2. 2.

Channel, $R=3.6$ Ohm m;

  3. 3.

Crevasse, $R=4.1$ Ohm m.

Second, we find boundaries between layers composed of pixels with equal
resistivities and use them as the input to the forward model.

\begin{overpic}[width=433.62pt]{fig-True_model_res_c_aspect} \put(0.75,1.5){
\leavevmode\hbox to345.26pt{\vbox
to53.12pt{\pgfpicture\makeatletter\hbox{\hskip 2.76054pt\lower-3.68073pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{}{{}}{}
{}{}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{2.13396pt}\pgfsys@invoke{ }{}{{
{\pgfsys@beginscope\pgfsys@setlinewidth{1.70717pt}\pgfsys@setdash{}{0.0pt}\pgfsys@roundcap\pgfsys@roundjoin{}
{}{}{} {}{}{}
\pgfsys@moveto{-2.76054pt}{3.68073pt}\pgfsys@curveto{-2.53049pt}{2.30045pt}{0.0pt}{0.23004pt}{0.69012pt}{0.0pt}\pgfsys@curveto{0.0pt}{-0.23004pt}{-2.53049pt}{-2.30045pt}{-2.76054pt}{-3.68073pt}\pgfsys@stroke\pgfsys@endscope}}
}{}{}{{}}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@lineto{339.88939pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{
}{{}{{}}{}{}{{}}{{{}}{{{}}{\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{339.88939pt}{0.0pt}\pgfsys@invoke{
}\pgfsys@invoke{ \lxSVG@closescope }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}{{}}}} \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{{}}{} {}{}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{2.13396pt}\pgfsys@invoke{
}{}{}{}{}{{}}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@lineto{0.0pt}{46.82588pt}\pgfsys@stroke\pgfsys@invoke{
}{{}{{}}{}{}{{}}{{{}}{{{}}{\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{0.0}{1.0}{-1.0}{0.0}{0.0pt}{46.82588pt}\pgfsys@invoke{
}\pgfsys@invoke{ \lxSVG@closescope }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}{{}}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}} } \put(95.0,-2.0){640
m.} \put(-5.0,11.0){32 m.} \put(-0.5,-0.5){0 m.} \put(2.0,3.0){
\leavevmode\hbox to1.42pt{\vbox to32.72pt{\pgfpicture\makeatletter\hbox{\hskip
0.71132pt\lower-0.71132pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{}{{}}{}
{}{}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{1.42264pt}\pgfsys@invoke{
}\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\pgfsys@color@rgb@stroke{0}{1}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{1}{0}\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{0,1,0}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@lineto{0.0pt}{31.2982pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}} \put(3.56,3.0){
\leavevmode\hbox to1.42pt{\vbox to32.72pt{\pgfpicture\makeatletter\hbox{\hskip
0.71132pt\lower-0.71132pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{}{{}}{}
{}{}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{1.42264pt}\pgfsys@invoke{
}\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\pgfsys@color@rgb@stroke{0}{1}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{1}{0}\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{0,1,0}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@lineto{0.0pt}{31.2982pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}} \put(5.12,3.0){
\leavevmode\hbox to1.42pt{\vbox to32.72pt{\pgfpicture\makeatletter\hbox{\hskip
0.71132pt\lower-0.71132pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{}{{}}{}
{}{}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{1.42264pt}\pgfsys@invoke{
}\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\pgfsys@color@rgb@stroke{0}{1}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{1}{0}\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{0,1,0}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@lineto{0.0pt}{31.2982pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}} \put(6.68,3.0){
\leavevmode\hbox to1.42pt{\vbox to32.72pt{\pgfpicture\makeatletter\hbox{\hskip
0.71132pt\lower-0.71132pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{}{{}}{}
{}{}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{1.42264pt}\pgfsys@invoke{
}\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\pgfsys@color@rgb@stroke{0}{1}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{1}{0}\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{0,1,0}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@lineto{0.0pt}{31.2982pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}} \put(8.25,3.0){
\leavevmode\hbox to1.42pt{\vbox to32.72pt{\pgfpicture\makeatletter\hbox{\hskip
0.71132pt\lower-0.71132pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{}{{}}{}
{}{}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{1.42264pt}\pgfsys@invoke{
}\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\pgfsys@color@rgb@stroke{0}{1}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{1}{0}\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{0,1,0}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@lineto{0.0pt}{31.2982pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}} \put(9.81,3.0){
\leavevmode\hbox to1.42pt{\vbox to32.72pt{\pgfpicture\makeatletter\hbox{\hskip
0.71132pt\lower-0.71132pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{}{{}}{}
{}{}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{1.42264pt}\pgfsys@invoke{
}\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\pgfsys@color@rgb@stroke{0}{1}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{1}{0}\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{0,1,0}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@lineto{0.0pt}{31.2982pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}} \put(11.375,3.0){
\leavevmode\hbox to1.42pt{\vbox to32.72pt{\pgfpicture\makeatletter\hbox{\hskip
0.71132pt\lower-0.71132pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{}{{}}{}
{}{}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{1.42264pt}\pgfsys@invoke{
}\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\pgfsys@color@rgb@stroke{0}{1}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{1}{0}\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{0,1,0}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@lineto{0.0pt}{31.2982pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}} \put(12.9375,3.0){
\leavevmode\hbox to1.42pt{\vbox to32.72pt{\pgfpicture\makeatletter\hbox{\hskip
0.71132pt\lower-0.71132pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{}{{}}{}
{}{}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{1.42264pt}\pgfsys@invoke{
}\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\pgfsys@color@rgb@stroke{0}{1}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{1}{0}\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{0,1,0}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@lineto{0.0pt}{31.2982pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}} \put(14.5,3.0){
\leavevmode\hbox to1.42pt{\vbox to32.72pt{\pgfpicture\makeatletter\hbox{\hskip
0.71132pt\lower-0.71132pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{}{{}}{}
{}{}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{1.42264pt}\pgfsys@invoke{
}\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\pgfsys@color@rgb@stroke{0}{1}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{1}{0}\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{0,1,0}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@lineto{0.0pt}{31.2982pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}} \put(1.7,7.25){
\leavevmode\hbox to47.8pt{\vbox to2.28pt{\pgfpicture\makeatletter\hbox{\hskip
1.1381pt\lower-1.1381pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{}{{}}{}
{}{}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{2.27621pt}\pgfsys@invoke{
}\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\pgfsys@color@rgb@stroke{1}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{1}{0}{0}\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{1,0,0}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@lineto{45.52458pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}} \put(15.0,7.5){
\leavevmode\hbox to34.54pt{\vbox to15.6pt{\pgfpicture\makeatletter\hbox{\hskip
0.2pt\lower-0.2pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{}{{}}{}{{}}{}{{}}{}{}{}{{}}{}{{}}{}{}{}\pgfsys@beginscope\pgfsys@invoke{
}\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\pgfsys@color@rgb@stroke{1}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{1}{0}{0}\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{1,0,0}\pgfsys@setdash{3.0pt,3.0pt}{0.0pt}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@curveto{5.69046pt}{0.39812pt}{11.38092pt}{1.59941pt}{17.07182pt}{3.62866pt}\pgfsys@curveto{22.76228pt}{6.52707pt}{28.45276pt}{10.35588pt}{34.14322pt}{15.20148pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}} } \put(15.0,3.1){
\leavevmode\hbox to34.54pt{\vbox to15.6pt{\pgfpicture\makeatletter\hbox{\hskip
0.2pt\lower-15.40147pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}{}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}{}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}{}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}{}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}{}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}{}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}{}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{{}}{}{{}}{}{{}}{{}}{{}}{{}}{}{{}}{}{}{}{{}}{{}}{{}}{{}}{}{{}}{}{}{}\pgfsys@beginscope\pgfsys@invoke{
}\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\pgfsys@color@rgb@stroke{1}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{1}{0}{0}\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{1,0,0}\pgfsys@setdash{3.0pt,3.0pt}{0.0pt}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@curveto{5.69046pt}{-0.39812pt}{11.38092pt}{-1.59941pt}{17.07182pt}{-3.62866pt}\pgfsys@curveto{22.76228pt}{-6.52707pt}{28.45276pt}{-10.35588pt}{34.14322pt}{-15.20148pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}} } \end{overpic}
Figure 3: Resistivity of earth model plotted in 1:2 aspect ratio. The green
lines shows the measurements and their extent illustrate the maximum
sensitivity depth. The full red line is the drilled well, and the dashed red
lines indicate the potential for geosteering.

##  4 Data assimilation

In the DSS for geosteering [3], one uses data assimilation to condition the
earth model to measurements made while drilling. The fundamental idea is that
if a poorly known earth model can be made consistent with measurements it will
provide more accurate forecasts, and, hence, provide a better basis for
decisions.

In this paper, the emphasis is placed on the data assimilation part of the
DSS. Especially data assimilation utilizing an efficient neural network model
for the synthetic logs, and an efficient GAN-generator for representing
complex earth models. To check whether this setup provides useful conditioned
models we consider two data assimilation algorithms. Firstly, the Markov Chain
Monte Carlo (MCMC), which is considered as a gold standard method for
sampling. Secondly, the ensemble randomized maximum likelihood (EnRML), an
approximate method suitable for DSS.

###  4.1 MCMC

A reliable method for sampling from a complex posterior distribution is the
MCMC technique. MCMC relates to the general framework of methods introduced in
[23] and [18] for Monte Carlo (MC) integration. One designs a Markov chain
that produce samples from the desired posterior distribution, and subsequently
utilize these samples for MC estimation. In this section, the adaptive
Metropolis-Hastings method, utilized in the numerical study, is introduced.
For more information on MCMC we refer the reader to [6], and references
therein.

Suppose we want samples from the un-normalized posterior distribution $F$,
which is the general case with the Bayesian method where the normalizing
factor often is very difficult to calculate. Assume that the current element
is $\vec{m}$, and that the chain proposes a move to $\vec{m}^{*}$, with
conditional probability density $q\left(\vec{m}^{*}|\vec{m}\right)$. The move
is performed with probability

$b\left(\vec{m},\vec{m}^{*}\right)=min\left(1,r\left(\vec{m},\vec{m}^{*}\right)\right)$
(1)

where the Hastings ratio is defined as

$r\left(\vec{m},\vec{m}^{*}\right)=\frac{F\left(\vec{m}^{*}\right)q\left(\vec{m}|\vec{m}^{*}\right)}{F\left(\vec{m}\right)q\left(\vec{m}^{*}|\vec{m}\right)}.$
(2)

If the move is not made $\vec{m}^{*}=\vec{m}$. This is the basis for the
Metropolis-Hastings method, and it can be shown that the method will generate
samples from the posterior distribution $F$.

The Metropolis-Hastings algorithm requires a choice of proposal distribution,
and some distributions will work better than others. Intuitively, one would
like to draw proposal samples from $F$. However, this is not possible since we
cannot sample from this distribution. However, one idea is to consider the
previous samples from the algorithm as approximate samples from $F$. With this
approach proposal samples are drawn from

$\vec{m}^{*}\sim\left(1-\beta\right)\mathcal{N}\left(\vec{m},\left(\frac{2.38^{2}}{N_{m}}\right)\tilde{C}_{\vec{m}}\right)+\beta\mathcal{N}\left(\vec{m},Q_{\vec{m}}\right),$
(3)

where $\tilde{C}_{\vec{m}}$ is the empirical covariance matrix calculated
utilizing all the preceding iterations of the Markov Chain, $Q_{\vec{m}}$ is
some fixed non-singular matrix and $0<\beta<1$. Note that $\beta=1$ until
$C_{\vec{m}}$ is well defined. This sampling method was applied in [13, 14].
It is well known that the MCMC requires a certain burn-in period, since the
initial samples are not from the posterior distribution. Hence, it is
necessary to monitor the convergence of the method. In this work, convergence
is monitored by assessing the maximum root statistic of the multivariate
potential scale reduction factor [5].

###  4.2 EnRML

The EnRML [16] has recently become one of the most successful methods for
automatic history matching of petroleum reservoirs. The EnRML is based on
minimization of an objective function using the ensemble approximation of the
sensitivity matrix. Hence, the EnRML can be formulated in many different ways.
In this study we utilize the approximate form of the Levenberg-Marquardt
method, introduced in [10].

Iteration number $i$ of the Levenberg-Marquardt method is given as

$\displaystyle\delta\vec{m}_{i}=$
$\displaystyle-\left[\left(1+\lambda_{i}\right)C^{-1}_{\vec{m}}+G^{T}_{i}C^{-1}_{\vec{d}}G_{i}\right]^{-1}$
(4)
$\displaystyle\times\left[C^{-1}_{\vec{m}}\left(\vec{m}_{i}-\vec{m}_{prior}\right)+G^{T}_{i}C_{\vec{d}}^{-1}\left(g\left(\vec{m}\right)-\left(\vec{d}_{obs}+\vec{\epsilon}\right)\right)\right]$
(5)

where $\lambda_{i}$ is the Levenberg-Marquardt multiplier, $G$ is the
sensitivity of data to the parameters, and
$\vec{\epsilon}\sim\mathcal{N}\left(0,C_{\vec{d}}\right)$ is a realization of
the measurement observation noise.

In the ensemble framework, we approximate $C_{\vec{m}}$ and $G$ using the
ensemble. To this end we define

$\tilde{G}=C^{1/2}_{sc}\Delta\vec{d}\left(\Delta\vec{m}\right)^{-1}$ (6)
$\tilde{C}_{\vec{m}}=\Delta\vec{m}\Delta\vec{m}^{T}$ (7)

where

$\Delta\vec{m}=\left[\vec{m}_{1},\dots,\vec{m}_{j},\dots,\vec{m}_{N}\right]\left(I_{N}-\frac{1}{N}11^{T}\right)/\sqrt{N-1},$
(8)
$\Delta\vec{d}=C_{sc}^{-1/2}\left[g\left(\vec{m}_{1}\right),\dots,g\left(\vec{m}_{j}\right),\dots,g\left(\vec{m}_{N}\right)\right]\left(I_{N}-\frac{1}{N}11^{T}\right)/\sqrt{N-1},$
(9)

$N$ denotes the ensemble size, and $C_{sc}$ is a diagonal matrix for scaling
the data, typically containing the measurement variance on the diagonal. We
get the approximate version of the Levenberg-Marquardt update equation by
inserting ensemble approximations of $G$ and $C_{\vec{m}}$, neglecting the
updates from the model mismatch term, substituting the prior precision matrix
$C^{-1}_{m}$ with $\tilde{C}^{-1}_{\vec{m}_{i}}$, and rewriting the equation
using the Sherman-Woodbury-Morrison matrix inversion formula [15]

$\delta\vec{m}_{i}=-\tilde{C}_{\vec{m}_{i}}\tilde{G}_{i}^{T}\left[\left(1+\lambda_{i}\right)C_{\vec{d}}+\tilde{G_{i}}\tilde{C}_{\vec{m}_{i}}\tilde{G}^{T}_{i}\right]^{-1}\left(g\left(\vec{m}\right)-\left(\vec{d}_{obs}+\vec{\epsilon}\right)\right).$
(10)

The update equation is simplified by calculating the truncated singular value
decomposition of $\Delta d$

$\Delta\vec{d}=U_{p}S_{p}V_{p}^{T},$ (11)

where the subscript $p$ indicates the number of singular values that are kept.
In this work, we define $p$ such that the cumulative sum of the $p$ singular
values equals 99% of the cumulative sum of all the singular values. Further,
to allow for correlated measurement errors, we substitute $C_{D}$ with the
ensemble approximation $\tilde{C}_{D}$

$\tilde{C}_{\vec{d}}=\Delta\vec{\epsilon}\Delta\vec{\epsilon}^{T},$ (12)

where

$\Delta\vec{\epsilon}=\left[\vec{\epsilon}_{1},\dots,\vec{\epsilon}_{j},\dots,\vec{\epsilon}_{N}\right]\left(I_{N}-\frac{1}{N}11^{T}\right)/\sqrt{N-1}.$
(13)

Inserted into (10) gives

$\displaystyle\delta\vec{m}_{i}=-$
$\displaystyle\Delta\vec{m}_{i}V_{p}\left[\left(1+\lambda_{i}\right)S_{p}^{-1}U_{p}^{T}C^{-1/2}_{scl}\Delta\vec{\epsilon}\Delta\vec{\epsilon}^{T}C_{scl}^{-1/2}U_{p}S_{P}^{-T}+I\right]$
$\displaystyle\left(U_{p}S^{-1}_{p}\right)^{T}C^{-1/2}_{sc}\left(g\left(\vec{m}\right)-\left(\vec{d}_{obs}+\vec{\epsilon}\right)\right)$
(14) $\displaystyle=-$
$\displaystyle\Delta\vec{m}_{i}V_{p}Z\left[\left(1+\lambda_{i}\right)\zeta+I\right]^{-1}\left(U_{p}S^{-1}_{p}Z\right)^{T}C^{-1/2}_{sc}\left(g\left(\vec{m}\right)-\left(\vec{d}_{obs}+\vec{\epsilon}\right)\right),$

where $Z$ and $\zeta$ are the eigenvectors and eigenvalues of
$S_{p}^{-1}U_{p}^{T}C^{-1/2}_{scl}\Delta\vec{\epsilon}\Delta\vec{\epsilon}^{T}C_{scl}^{-1/2}U_{p}S_{P}^{-T}$.

The iterative scheme is run until it is converged. Here, we consider the
method to be converged when the relative difference in the data misfit is
below a given threshold.

##  5 Numerical Results

We, throughout, utilize the generative neural network, introduced in section
2, to represent the poorly known earth model. Hence, realizations of the earth
model is generated by applying the generative network to parameters sampled
from the multivariate standard distribution,
$\vec{m}\sim\mathcal{N}\left(0,C_{\vec{m}}\right)$. The numerical
investigation considers a well drilled horizontally, approximately, in the
center of the model and through the 9 first grid-cells. For each drilled grid-
cell we simulate measurements using the model introduced in section 3.

Throughout the investigation, we consider a diagonal $C_{\vec{m}}$ with
elements equal to $1\times 10^{-6}$. Figure 4 shows two random earth model
realizations from the prior model. From the figure, we observe that this setup
provides significant variation in the earth model.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

Figure 4: Rows 1-3 show three realizations from the prior ensemble. Each row
shows the three output channels (e.g. fig. a-c) of GAN corresponding to
probability of facies (Background/Channel/Crevasse), and the derived
resistivity image (e.g. fig. d). Red stars indicate measurement position.

The synthetic true earth model is also drawn from the prior model. All 3
channels and the derived resistivity of the synthetic truth is illustrated in
figure 5. The true observations are simulated using the true earth model.

(a)

(b)

(c)

(d)

Figure 5: (a)-(c): Facies probabilities from GAN for the synthetic truth. (d):
Derived resistivity for the synthetic truth. Red stars indicate measurement
positions.

For each of the 9 measurements positions the measurement standard deviation is
given as 5% of the measurement value. In addition, we let the measurements at
each position be correlated. Assuming equidistance between measurements, the
correlation length equals 10 times the inter-measurement distance. We further
assume that the measurements at two different well positions are uncorrelated.

We conduct two numerical experiments. Firstly, we conduct a MCMC run to
properly characterize the posterior distribution. Here, 8 chains are run in
parallel for $10^{6}$ iteration. At that point the chains were converged, and
we extract samples from the posterior by, for each of the 8 chains, removing
the first half of the chain and retaining every 100 iteration from the second
half of the chain. Hence, leaving $4\times 10^{4}$ samples from the posterior
distribution. Secondly, we estimate the posterior distribution using the EnRML
method introduced in section 4.2. Due to the fast simulation time we utilized
an ensemble size of $N=500$, and in addition we applied the correlation based
localization technique introduced in [21]. The method is allowed to iterate
until the relative improvement of the updates is less than $1\%$. When showing
the numerical results we will only plot the values of the derived resistivity.

###  5.1 MCMC

Based on the posterior realizations obtained by the MCMC, we calculate the
posterior mean and posterior standard deviation, shown in figure 6 (a) and
(b). In figure 6 (c) and (d) we plot two random realizations from the
posterior.

(a)

(b)

(c)

(d)

Figure 6: (a): MCMC mean of the resistivity. (b): MCMC standard deviation of
the resistivity. (c): Resistivity for MCMC posterior realization 1. (d):
Resistivity for MCMC posterior realization 2. Red stars indicate measurement
positions.

The MCMC result illustrates that the generative neural network can be utilized
for data assimilation, and be successfully conditioned to measurements
utilizing a neural network proxy model. The model standard deviation is
significantly reduced close to the drill bit, and also ahead of the drill bit
position. Moreover, the mean value shows that the correct resitivity is
identified in these areas. Note that the posterior still has significant
variance in most parts of the field. The areas directly around the drill bit
have not obtained sufficient reduction of the standard deviation. We do not
properly understand why this is so, however, it indicates that there
measurements are less sensitive to the region near the drill bit.

###  5.2 EnRML

After the EnRML has converged we calculate the mean and standard deviation
from the ensemble. These are approximations to the true posterior mean and
standard deviations, and are shown in figure 7 (a) and (b). In addition, the
results from two random realizations are illustrated in figure 7 (c) and (d).

(a)

(b)

(c)

(d)

Figure 7: (a): Mean resistivity from the EnRML. (b): Standard deviation of the
resistivity from the EnRML. (c): Resistivity for EnRML realization 1. (d):
Resitivity for EnRML realization 2. Red stars indicate measurement positions.

The results for the EnRML show a significant reduction in the standard
deviation around the drill bit, and from the mean model we observe that the
correct facies is identified for these regions. There is some reduction of the
standard deviation ahead of the drill bit, but we observe that the variance is
retained in most of the field.

##  6 Summary and conclusions

In this paper, we have demonstrated that two essential parts, the earth model
and the simulated log, of an ensemble based DSS system can be substituted with
neural networks. For the earth model we utilize the GAN trained with images
from a realistic geological setting, for the simulated log we use a deep
neural network trained using a large set of simulations from a commercial
tool. The setup redistributes the computational cost from on-line to off-line
calculations, enabling complex earth models utilizing simulated logs with high
accuracy to be used in the real time DSS. The numerical results illustrate
that DSS, equipped with the GAN, provides good predictions ahead of drilling
when conditioning to only measurements with sideways sensitivity.

The proposed approach has many beneficial factors. Firstly, a GAN provides
large flexibility for defining the geological setting. Here, we consider three
different facies, but one can easily imagine the inclusion of features like
faults and pinch-outs as well as smoothly-varying properties. Secondly, we
only need to condition a few parameters with Gaussian distribution to the
measurements, which is very beneficial for the ensemble based DA approach.
Thirdly, since we are utilizing a neural network model to generate the
simulated log the computational cost of simulating a single ensemble member is
very low. Hence, the proposed approach can utilize a large ensemble for the DA
part.

The numerical experiments illustrated that the setup provides a reasonable
posterior distribution, and we can estimate this using the EnRML
approximation. In this study we have only considered a single part of the DSS,
namely the conditioning of the earth model to measured data. However, due to
the promising results, the developments shown in this paper will integrate
with the framework developed in [3]. Hence, allowing DSS under much more
complex geological setting. Demonstrating the complete DSS with the proposed
setup is left for future work.

## Acknowledgments

The authors are supported by the research project ’Geosteering for IOR’ (NFR-
Petromaks2 project no. 268122) which is funded by the Research Council of
Norway, Aker BP, Equinor, Vår Energi and Baker Hughes Norway.

## References

  * [1] Al-Fawwaz, A., Al-Yosef, O., Al-Qudaihy, D., Al-Shobaili, Y., Al-Faraj, H., Maeso, C., Roberts, I.: Increased Net to Gross Ratio as the Result of an Advanced Well Placement Process Utilizing Real-Time Density Images. In: IADC/SPE Asia Pacific Drilling Technology Conference and Exhibition. pp. 151–160. Society of Petroleum Engineers (apr 2004). https://doi.org/10.2118/87979-MS, http://www.onepetro.org/doi/10.2118/87979-MS
  * [2] Alyaev, S., Shahriari, M., Pardo, D., Omella, A.J., Larsen, D., Jahani, N., Suter, E.: Modeling extra-deep em logs using a deep neural network. arXiv preprint arXiv:2005.08919, accepted in SEG Geophysics (2020)
  * [3] Alyaev, S., Suter, E., Bratvold, R.B., Hong, A., Luo, X., Fossum, K.: A decision support system for multi-target geosteering. Journal of Petroleum Science and Engineering 183(August), 106381 (dec 2019). https://doi.org/10.1016/j.petrol.2019.106381, https://doi.org/10.1016/j.petrol.2019.106381%**** ̵ms.bbl ̵Line ̵25 ̵****https://linkinghub.elsevier.com/retrieve/pii/S0920410519308022
  * [4] Arjovsky, M., Chintala, S., Bottou, L.: Wasserstein GAN. arXiv (jan 2017), http://arxiv.org/abs/1701.07875
  * [5] Brooks, S.P., Gelman, A.: General Methods for Monitoring Convergence of Iterative Simulations. J. Comput. Graph. Stat. 7(4), 434 (dec 1998). https://doi.org/10.2307/1390675, http://www.tandfonline.com/doi/abs/10.1080/10618600.1998.10474787http://www.jstor.org/stable/1390675?origin=crossref
  * [6] Brooks, S.P., Gelman, A., Jones, G.L., Meng, X.L. (eds.): Handbook of Markov Chain Monte Carlo. Chapman and Hall / CRC (2011)
  * [7] Chan, S., Elsheikh, A.H.: Parametric generation of conditional geological realizations using generative neural networks. Computational Geosciences 23(5), 925–952 (oct 2019). https://doi.org/10.1007/s10596-019-09850-7, http://link.springer.com/10.1007/s10596-019-09850-7
  * [8] Chan, S., Elsheikh, A.H.: Parametrization of Stochastic Inputs Using Generative Adversarial Networks With Application in Geology. Frontiers in Water 2(March), 1–21 (mar 2020). https://doi.org/10.3389/frwa.2020.00005, https://www.frontiersin.org/article/10.3389/frwa.2020.00005/full
  * [9] Chen, Y., Lorentzen, R.J., Vefring, E.H.: Optimization of Well Trajectory Under Uncertainty for Proactive Geosteering. SPE Journal 20(02), 368–383 (apr 2015). https://doi.org/10.2118/172497-PA, https://onepetro.org/SJ/article/20/02/368/206467/Optimization-of-Well-Trajectory-Under-Uncertainty
  * [10] Chen, Y., Oliver, D.S.: Levenberg–Marquardt forms of the iterative ensemble smoother for efficient history matching and uncertainty quantification. Comput. Geosci. (may 2013). https://doi.org/10.1007/s10596-013-9351-5, http://link.springer.com/10.1007/s10596-013-9351-5
  * [11] Cole, R.D., Cumella, S.: Sand-body architecture in the lower Williams Fork Formation (Upper Cretaceous), Coal Canyon, Colorado, with comparison to the Piceance Basin subsurface. Mountain Geologist 42, 85–107 (2005)
  * [12] Evensen, G.: Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics. J. Geophys. Res. 99(C5), 10143 (1994). https://doi.org/10.1029/94JC00572, http://doi.wiley.com/10.1029/94JC00572
  * [13] Fossum, K., Mannseth, T.: Parameter sampling capabilities of sequential and simultaneous data assimilation: II. Statistical analysis of numerical results. Inverse Problems 30(11), 114003 (nov 2014). https://doi.org/10.1088/0266-5611/30/11/114003, http://stacks.iop.org/0266-5611/30/i=11/a=114003?key=crossref.2e2654e480ea9cb9eab6d8c245d1dc46
  * [14] Fossum, K., Mannseth, T.: Assessment of ordered sequential data assimilation. Computational Geosciences 19(4) (2015). https://doi.org/10.1007/s10596-015-9492-9, http://link.springer.com/10.1007/s10596-015-9492-9
  * [15] Golub, G.H., Van Loan, C.F.: Matrix Computations. Johns Hopkins series in the mathematical sciences, The Johns Hopkins University Press, Baltimore (1983)
  * [16] Gu, Y., Oliver, D.S.: An Iterative Ensemble Kalman Filter for Multiphase Fluid Flow Data Assimilation. SPE J. 12(4), 438–46 (nov 2007). https://doi.org/10.2118/108438-PA, http://www.spe.org/ejournals/jsp/journalapp.jsp?pageType=Preview–&˝jid=ESJ–&˝mid=SPE-108438-PA–&˝pdfChronicleId=090147628014cce3
  * [17] Guevara, A.I., Sandoval, J., Guerrero, M., Manrique, C.A.: Milestone in Production Using Proactive Azimuthal Deep-Resistivity Sensor Combined With Advanced Geosteering Techniques: Tarapoa Block, Ecuador. In: SPE Latin America and Caribbean Petroleum Engineering Conference. vol. 2, pp. 1508–1520. Society of Petroleum Engineers (apr 2012). https://doi.org/10.2118/153580-MS, http://www.onepetro.org/doi/10.2118/153580-MS
  * [18] Hastings, W.K.: Monte Carlo Sampling Methods Using Markov Chains and Their Applications. Biometrika 57(1), 97 (apr 1970). https://doi.org/10.2307/2334940, http://www.jstor.org/stable/2334940?origin=crossref
  * [19] Janwadkar, S., Thomas, M., Privott, S., Tehan, R., Carlson, L., Spear, W., Setiadarma, A.: Reservoir-Navigation System and Drilling Technology Maximize Productivity and Drilling Performance in the Granite Wash, US Midcontinent. SPE Drilling & Completion 27(01), 22–31 (mar 2012). https://doi.org/10.2118/140073-PA, https://onepetro.org/DC/article/27/01/22/198159/Reservoir-Navigation-System-and-Drilling
  * [20] Kullawan, K., Bratvold, R., Bickel, J.E.: A decision analytic approach to geosteering operations. SPE Drilling & Completion 29(01), 36–46 (2014)
  * [21] Luo, X., Bhakta, T., Nævdal, G.: Correlation-Based Adaptive Localization With Applications to Ensemble-Based 4D-Seismic History Matching. SPE JOURNAL 23(2), 396–427 (apr 2018). https://doi.org/10.2118/185936-PA
  * [22] Luo, X., Eliasson, P., Alyaev, S., Romdhane, A., Suter, E., Querendez, E., Vefring, E.: An ensemble-based framework for proactive geosteering. SPWLA 56th Annual Logging Symposium 2015 (2015)
  * [23] Metropolis, N., Rosenbluth, A.W., Rosenbluth, M.N., Teller, A.H., Teller, E.: Equation of State Calculations by Fast Computing Machines. The Journal of Chemical Physics 21(6), 1087 (1953). https://doi.org/10.1063/1.1699114, http://link.aip.org/link/JCPSA6/v21/i6/p1087/s1–&˝Agg=doi
  * [24] Panjaitan, H.: Sand-body dimensions in outcrop and subsurface, Lower Williams Fork Formation, Piceance Basin, Colorado. Master’s thesis, Colorado School of Mines (2006)
  * [25] Pranter, M.J., Cole, R.D., Panjaitan, H., Sommer, N.K.: Sandstone-body dimensions in a lower coastal-plain depositional setting: Lower Williams Fork formation, Coal Canyon, Piceance Basin, Colorado. American Association of Petroleum Geologists Bulletin 93(10), 1379–1401 (2009). https://doi.org/10.1306/06240908173
  * [26] Pranter, M.J., Hewlett, A.C., Cole, R.D., Wang, H., Gilman, J.: Fluvial architecture and connectivity of the Williams Fork Formation: use of outcrop analogues for stratigraphic characterization and reservoir modelling. Geological Society, London, Special Publications 387(1), 57–83 (2014). https://doi.org/10.1144/sp387.1
  * [27] Pranter, M.J., Sommer, N.K.: Static connectivity of fluvial sandstones in a lower coastal-plain setting: An example from the Upper Cretaceous lower Williams Fork Formation, Piceance Basin, Colorado. AAPG Bulletin 95(6), 899–923 (2011). https://doi.org/10.1306/12091010008
  * [28] Sebacher, B., Stordal, A.S., Hanea, R.: Bridging multipoint statistics and truncated Gaussian fields for improved estimation of channelized reservoirs with ensemble methods. Computational Geosciences pp. 341–369 (2015). https://doi.org/10.1007/s10596-014-9466-3, http://link.springer.com/10.1007/s10596-014-9466-3
  * [29] Sviridov, M., Mosin, A., Antonov, Y., Nikitenko, M., Martakov, S., Rabinovich, M.B., Baker Hughes, BP: New Software for Processing of LWD Extradeep Resistivity and Azimuthal Resistivity Data. SPE Reservoir Evaluation & Engineering 17(May) (2014). https://doi.org/10.2118/160257-PA
  * [30] Trampush, S.M., Hajek, E.A., Straub, K.M., Chamberlin, E.P.: Identifying autogenic sedimentation in fluvial-deltaic stratigraphy: Evaluating the effect of outcrop-quality data on the compensation statistic. Journal of Geophysical Research: Earth Surface 122(1), 91–113 (2017). https://doi.org/10.1002/2016JF004067


#  Investigation of Görtler vortices in high-speed boundary layers via an
efficient numerical solution to the non-linear boundary region equations

Omar Es-Sahli email<EMAIL_ADDRESS>Department of Aerospace Engineering,
Mississippi State University, US  Adrian Sescu email<EMAIL_ADDRESS>Department of Aerospace Engineering, Mississippi State University, US
Mohammed Afsar  Department of Mechanical & Aerospace Engineering, Strathclyde
University, UK  Yuji Hattori  Institute of Fluid Science, Tohoku University,
Japan

###### Abstract

Streamwise vortices and the associated streaks evolve in boundary layers over
flat or concave surfaces due to disturbances initiated upstream or triggered
by the wall surface. Following the transient growth phase, the fully-developed
vortex structures become susceptible to inviscid secondary instabilities
resulting in early transition to turbulence via ‘bursting’ processes. In high-
speed boundary layers, more complications arise due to compressibility and
thermal effects, which become more significant for higher Mach numbers. In
this paper, we study Görtler vortices developing in high-speed boundary layers
using the boundary region equations (BRE) formalism, which we solve using an
efficient numerical algorithm. Streaks are excited using a small transpiration
velocity at the wall. Our BRE-based algorithm is found to be superior to
direct numerical simulation (DNS) and ad-hoc nonlinear parabolized stability
equation (PSE) models. BRE solutions are less computationally costly than a
full DNS and have a more rigorous theoretical foundation than PSE-based
models. For example, the full development of a Görtler vortex system in high-
speed boundary layers can be predicted in a matter of minutes using a single
processor via the BRE approach. This substantial reduction in calculation time
is one of the major achievements of this work. We show, among other things,
that it allows investigation into feedback control in reasonable total
computational times. We investigate the development of the Görtler vortex
system via the BRE solution with feedback control parametrically at various
freestream Mach numbers $M_{\infty}$ and spanwise separations $\lambda$ of the
inflow disturbances.

##  1 Introduction

Streaks formation in pre-transitional boundary layer flows over flat or curved
surfaces occur when the amplitude of the freestream disturbances is
sufficiently large, or the height of the roughness elements of the surface
exceeds a certain critical value. The streamwise velocity component exhibits
elongated, so-called ‘streaky’, features characterized by adjacent regions of
acceleration (high-speed) and deceleration (low-speed) of fluid particles
(Kendall [15], Matsubara & Alfredsson [25], or Landahl [17]). Elongated
streaks in the form of streamwise (Görtler) vortices also appear inside a
boundary layer flow along a concave surface due to the imbalance between
radial pressure gradients posed by the wall, and centrifugal forces (e.g.,
Görtler [8], Hall [9, 10, 11], Swearingen & Blackwelder [36], Malik & Hussaini
[22], Li & Malik [19], Wu et al. [40], Sescu et al. [31, 32], Marensi & Ricco
[23], Xu et al. [41]). For highly curved walls, for example, vortex formation
occurs more rapidly and can significantly alter the mean flow causing the
laminar flow to breakdown into turbulence. The literature treating and
discussing Görtler vortices evolving in incompressible boundary layer flows is
rich; we encourage the reader to consult the aforementioned studies and many
others that are included in their respective references lists.

Görtler vortices in compressible boundary layers have been studied for quite a
while motivated by numerous engineering applications such as high-speed flow
in engine intakes, flows over the concave surface a turbomachinery blades, or
flows over the walls of supersonic nozzles, for example. They were studied by
using parallel flow theory by Kobayashi and Kohama [16]. Non-parallel effects
were studied by El-Hady and Verma [6], Hall and Fu [13], and Hall and Malik
[14]. Spall and Malik [35] further improved the parallel eigenvalue framework
by adding initial conditions to the partial differential equations, assuming
zero amplitude perturbations in the external boundary layer (see slight
modifications of this approach in Wadey [38] or Dando and Seddougui [4]). The
number of experiments involving Görtler vortices developing in compressible
boundary layers is not as large as the number of experiments performed in the
incompressible regime. Worth to mention are the experiments of De Luca et al.
[5], Ciolkosz and Spina [3], or Wang et al. [39].

While the aforementioned theoretical and numerical studies are relatively old,
there has been a resurgence of interest in Görtler vortices evolving in high-
speed boundary layers in recent years. For instance, Li et al. [20] studied
the linear and nonlinear growth of Görtler vortices in hypersonic boundary
layers using the parabolized stability equations (PSE), linear stability
analysis, and direct numerical simulations (DNS). They identified multiple
sets of unstable secondary instability modes and investigated their linear and
nonlinear spatial development. DNS was used to explore the onset of transition
in order to determine the most important physical mechanisms associated with
high-speed boundary layers. Ren and Fu [26] conducted a series of numerical
computations to investigate the fundamental, subharmonic and detuned secondary
instabilities of Görtler vortices in high-speed boundary layer flows with a
focus on the Mach number effect. They found that the growth rate associated
with Görtler vortices decreases with the Mach number and contributes to the
appearance of the trapped-layer mode in the primary instability. Chen et al.
[2] employed DNS and linear stability analysis to explore the transition of
stationary Görtler vortices in high-speed boundary layer flows by exciting the
instability using steady blowing and suction on the wall, similar to the
disturbance employed in this work. It was shown that, depending on its
frequency and wavelength, the first Mack mode can turn into either a varicose
or sinuous mode streak instability, while the second Mack mode turns into a
varicose mode.

In Li et al. [21], DNS and linear secondary instability theory were utilized
to study Görtler vortices and their associated secondary instabilities in the
hypersonic boundary layer flow over a cone featuring a concave aft body. It
was revealed that the secondary instability is predominantly characterized by
sinuous modes that concentrates in the wall-normal internal shear layer in the
lower portion of the mushroom shapes. A form of compressible boundary region
equations has been used by Viaro and Ricco [37] to study the evolution of
Görtler vortices in compressible boundary layer flows. They used a previous
approach that was derived by Ricco and Wu [28] as an extension of the
incompressible approach of Leib et al. [18] to the compressible regime. Within
this formalism, both the effect of the initial conditions and the boundary
condition at the top of the boundary layer are correctly accounted for. They
also explained the formation and growth of thermal streaks, which are thought
to play a significant role in the secondary instability process (see also
Ricco [27]). Ricco, Tran & Ye [30] and Ricco, Shah & Hicks [29] further
studied the influence of wall heat transfer and wall suction, respectively, on
the thermal streaks. Finally, we mention the work of Song et al. [34], where
the evolution of first and second Mack mode inside a compressible boundary
layer involving Görtler vortices is studied in detail via DNS and PSE.
Nevertheless, DNS falls short on computation efficiency for large Reynolds
numbers that are typically of interest, whereas PSE-based approaches remain
largely ad-hoc and suffer from convergence problems [1]. This justifies
investigation into a more efficient numerical solution based on a more robust
theoretical foundation.

The present paper focuses on investigating the streamwise vortices and the
associated Görtler vortices that develop in high-speed boundary layers over
concave surfaces using the full nonlinear BRE. The BRE represent the high
Reynolds number asymptotic limit of the Navier-Stokes (N-S) equations under
the assumption that the streamwise wavenumbers of the disturbances are much
smaller than those associated with the crossflow disturbances. This set of
equations is parabolic in the streamwise direction allowing for the
application of a straightforward marching procedure along the streamwise
direction. To this end, upstream conditions have to be imposed to start the
calculation, and here we accomplish this by a small disturbance at the wall in
the form of wall transpiration. At the upstream boundary, we impose mean flow
profiles obtained from a similarity solution applied to compressible boundary
layer equations (equivalent to the Blasius solution for an incompressible
flow). Since a concave wall is considered in the analysis, the disturbances
take the form of Görtler instabilities featuring counter rotating pairs of
vortices and associated streaks, with streamwise velocity contours plotted in
crossflow sections resembling mushroom shapes. We analyze and quantify the
evolution of these streaks via contour plots of velocity and temperature in
crossflow planes as well as vortex energy, wall shear stress, and wall heat
flow distributions versus the streamwise coordinate. We emphasize that the
numerical algorithm that is applied to BRE is very efficient, making it
suitable for timely parametric studies.

The rest of the paper is organized as follows: In section II, we introduce and
describe the mathematical model as we discuss the scaling of various
independent/dependent variables, the appropriate initial and boundary
conditions, as well as the numerical algorithm. In section III, we report and
discuss results for various freestream Mach numbers and spanwise separations,
in supersonic and hypersonic regimes. Section IV includes concluding remarks.

##  2 Problem formulation and numerical algorithm

We consider a compressible flow of uniform velocity $V_{\infty}^{*}$ and
temperature $T_{\infty}^{*}$ over a curved surface. The air is treated as a
perfect gas such that the speed of sound $c_{\infty}^{*}=\sqrt{\gamma
RT_{\infty}^{*}}$, where $\gamma$ = 1.4 is the ratio of the specific heats,
and $R=287.05Nm/(kgK)$ is the universal gas constant. The Mach number,
$M_{\infty}=V_{\infty}^{*}/c_{\infty}^{*}$, is assumed to be of order one.
Note that the superscript * and the subscript $\infty$ symbols indicate
dimensional freestream quantities. The flow is divided into four regions as in
Leib et al. [18], Ricco & Wu [28], and Marensi et al. [24] (see figure 1).
Region I is in proximity to the the leading edge, outside of the boundary
layer; the flow is assumed inviscid and the disturbances are treated as small
perturbations of the base flow. Region II is the boundary layer in the
vicinity of the leading edge with thickness much smaller than the spanwise
separation associated with the freestream disturbances; the disturbances are
governed by the linearized boundary region equations and the diffusion in the
spanwise direction is of the same order of magnitude as that in the wall-
normal direction. Region III is the viscous region that follows downstream
region II; the boundary layer thickness is of the same order of magnitude as
the spanwise separation and the flow is governed by the BRE derived from the
full N-S equations by neglecting the streamwise pressure-gradient and the
streamwise viscous diffusion (since they evolve on a slow streamwise scale
consistent with the vortex structure being elongated in that direction).
Region IV, outside the boundary layer above region III, is assumed inviscid
since the viscous effects are negligible; the flow is influenced at leading
order by the displacement effect due to the increased thickness of the viscous
layer underneath it.

Figure 1:  The flow domains illustrating the asymptotic structure and the four
regions of the flow domain. The box distinguishes regions I and II with I
inside the box and outside the boundary layer and II inside the box and inside
the boundary layer.

The focus in this paper is on region III, where the streamwise wavenumber of
disturbances are expected to be small, as previous incompressible Görtler flow
studies suggest (e.g. Wu et al [40], Marensi & Ricco [23]), and the flow is
governed by the compressible boundary region equations.

###  2.1 Scalings

All dimensional spatial coordinates $(x^{*},y^{*},z^{*})$ are normalized by
the spanwise separation $\lambda^{*}$ of the disturbances, while the dependent
variables by their respective freestream values, and the pressure field is
normalized by the dynamic pressure:

$\displaystyle\bar{t}=\frac{t^{*}}{\lambda^{*}/V_{\infty}^{*}};\hskip
11.38109pt\bar{x}=\frac{x^{*}}{\lambda^{*}};\hskip
11.38109pt\bar{y}=\frac{y^{*}}{\lambda^{*}};\hskip
11.38109pt\bar{z}=\frac{z^{*}}{\lambda^{*}};\hskip
11.38109pt\bar{u}=\frac{u^{*}}{V_{\infty}^{*}};\hskip
11.38109pt\bar{v}=\frac{v^{*}}{V_{\infty}^{*}};\hskip
11.38109pt\bar{w}=\frac{w^{*}}{V_{\infty}^{*}};$
$\displaystyle\bar{\rho}=\frac{\rho^{*}}{\rho_{\infty}^{*}};\hskip
11.38109pt\bar{p}=\frac{p^{*}-p_{\infty}^{*}}{\rho_{\infty}^{*}V_{\infty}^{*2}};\hskip
11.38109pt\bar{T}=\frac{T^{*}}{T_{\infty}^{*}};\hskip
11.38109pt\bar{\mu}=\frac{\mu^{*}}{\mu_{\infty}^{*}};\hskip
11.38109pt\bar{k}=\frac{k^{*}}{k_{\infty}^{*}},$ (1)

where $(u^{*},v^{*},w^{*})$ are the velocity components, $\rho^{*}$ the
density, $p^{*}$ is pressure, $T^{*}$ temperature, $\mu^{*}$ dynamic
viscosity, and $k^{*}$ thermal conductivity. All quantities with $\infty$ at
the subscript represent conditions at infinity.

Reynolds number based on the spanwise separation, Mach number and Prandtl
number are defined as

$\displaystyle
R_{\lambda}=\frac{\rho_{\infty}^{*}V_{\infty}^{*}\lambda^{*}}{\mu_{\infty}^{*}},\hskip
14.22636ptM_{\infty}=\frac{V_{\infty}^{*}}{c_{\infty}^{*}},\hskip
14.22636ptPr=\frac{\mu_{\infty}^{*}C_{p}}{k_{\infty}^{*}}$ (2)

where $\mu_{\infty}^{*}$, $c_{\infty}^{*}$ and $k_{\infty}^{*}$ are freestream
dynamic viscosity, speed of sound and thermal conductivity, respectively, and
$C_{p}$ is the specific heat at constant pressure. For boundary layer flows
over curved surfaces, we define the global Görtler number as
$G_{\lambda}=R_{\lambda}^{2}\lambda^{*}/r^{*}$, where $r^{*}$ is the radius of
the curvature.

###  2.2 Compressible boundary-region equations

For a full compressible flow, the primitive form of the N-S equations in non-
dimensional variables, are considered here in the form

$\displaystyle\frac{D\bar{\rho}}{Dt}+\rho\left(\frac{\partial\bar{u}}{\partial\bar{x}}+\frac{\partial\bar{v}}{\partial\bar{y}}+\frac{\partial\bar{w}}{\partial\bar{z}}\right)=0$
(3)
$\bar{\rho}\frac{D\bar{u}}{D\bar{t}}=-\frac{\partial\bar{p}}{\partial\bar{x}}+\frac{1}{Re_{\lambda}}\frac{\partial}{\partial\bar{x}}\left[\frac{2}{3}\mu\left(2\frac{\partial\bar{u}}{\partial\bar{x}}-\frac{\partial\bar{v}}{\partial\bar{y}}-\frac{\partial\bar{w}}{\partial\bar{z}}\right)\right]+\frac{\partial}{\partial\bar{y}}\left[\mu\left(\frac{\partial\bar{u}}{\partial\bar{y}}+\frac{\partial\bar{v}}{\partial\bar{x}}\right)\right]+\frac{\partial}{\partial\bar{z}}\left[\mu\left(\frac{\partial\bar{w}}{\partial\bar{x}}+\frac{\partial\bar{u}}{\partial\bar{z}}\right)\right]$
(4)
$\bar{\rho}\frac{D\bar{v}}{D\bar{t}}=-\frac{\partial\bar{p}}{\partial\bar{y}}+\frac{1}{Re_{\lambda}}\frac{\partial}{\partial\bar{y}}\left[\frac{2}{3}\mu\left(2\frac{\partial\bar{v}}{\partial\bar{y}}-\frac{\partial\bar{u}}{\partial\bar{x}}-\frac{\partial\bar{w}}{\partial\bar{z}}\right)\right]+\frac{\partial}{\partial\bar{x}}\left[\mu\left(\frac{\partial\bar{v}}{\partial\bar{x}}+\frac{\partial\bar{u}}{\partial\bar{y}}\right)\right]+\frac{\partial}{\partial\bar{z}}\left[\mu\left(\frac{\partial\bar{v}}{\partial\bar{z}}+\frac{\partial\bar{w}}{\partial\bar{y}}\right)\right]$
(5)
$\bar{\rho}\frac{D\bar{w}}{D\bar{t}}=-\frac{\partial\bar{p}}{\partial\bar{z}}+\frac{1}{Re_{\lambda}}\frac{\partial}{\partial\bar{z}}\left[\frac{2}{3}\mu\left(2\frac{\partial\bar{w}}{\partial\bar{z}}-\frac{\partial\bar{u}}{\partial\bar{x}}-\frac{\partial\bar{v}}{\partial\bar{y}}\right)\right]+\frac{\partial}{\partial\bar{x}}\left[\mu\left(\frac{\partial\bar{w}}{\partial\bar{x}}+\frac{\partial\bar{u}}{\partial\bar{z}}\right)\right]+\frac{\partial}{\partial\bar{y}}\left[\mu\left(\frac{\partial\bar{v}}{\partial\bar{z}}+\frac{\partial\bar{w}}{\partial\bar{y}}\right)\right]$
(6)
$\displaystyle\bar{\rho}\frac{D\bar{T}}{D\bar{t}}=\frac{1}{PrRe_{\lambda}}\left[\frac{\partial}{\partial\bar{x}}\left(k\frac{\partial\bar{T}}{\partial\bar{x}}\right)+\frac{\partial}{\partial\bar{y}}\left(k\frac{\partial\bar{T}}{\partial\bar{y}}\right)+\frac{\partial}{\partial\bar{z}}\left(k\frac{\partial\bar{T}}{\partial\bar{z}}\right)\right]$
$\displaystyle\quad-(\gamma-1)M_{\infty}^{2}\left[p\left(\frac{\partial\bar{u}}{\partial\bar{x}}+\frac{\partial\bar{v}}{\partial\bar{y}}+\frac{\partial\bar{w}}{\partial\bar{z}}\right)-\frac{2}{3}\mu\left(\frac{\partial\bar{u}}{\partial\bar{x}}+\frac{\partial\bar{v}}{\partial\bar{y}}+\frac{\partial\bar{w}}{\partial\bar{z}}\right)^{2}\right]+(\gamma-1)M_{\infty}^{2}\frac{\mu}{Re_{\lambda}}\Psi$
(7)

where

$\displaystyle\Psi=2\left(\frac{\partial\bar{u}}{\partial\bar{x}}\right)^{2}+2\left(\frac{\partial\bar{v}}{\partial\bar{y}}\right)^{2}+2\left(\frac{\partial\bar{w}}{\partial\bar{z}}\right)^{2}+\left(\frac{\partial\bar{u}}{\partial\bar{y}}+\frac{\partial\bar{v}}{\partial\bar{x}}\right)^{2}+\left(\frac{\partial\bar{w}}{\partial\bar{x}}+\frac{\partial\bar{u}}{\partial\bar{z}}\right)^{2}+\left(\frac{\partial\bar{v}}{\partial\bar{z}}+\frac{\partial\bar{w}}{\partial\bar{y}}\right)^{2}$
(8)

is the dissipation function, and
$D/D\bar{t}=\partial/\partial\bar{t}+\bar{u}\partial/\partial\bar{x}+\bar{v}\partial/\partial\bar{y}+\bar{w}\partial/\partial\bar{z}$
is the substantial derivative (in what follows, we consider the steady-state
case for the N-S equations where all time derivative terms drop out, i.e.
$\partial/\partial\bar{t}=0$). The pressure, $\bar{p}$, the temperature,
$\bar{T}$ and the density, $\bar{\rho}$, of the fluid are combined in the
equation of state in non-dimensional form, $\bar{p}=\bar{\rho}\bar{T}/\gamma
M_{\infty}^{2}$, assuming non-chemically-reacting flows. The dynamic viscosity
and thermal conductivity $k$ is linked to the temperature using the power law
in dimensionless form: $\bar{\mu}=\bar{T}^{b};\hskip
2.84526pt\bar{k}=C_{p}\bar{\mu}/Pr$, where $b=0.76$ (Ricco & Wu [28]),
$C_{p}=\gamma R/(\gamma-1)$, $\gamma=1.4$, and $Pr=0.72$ for air.

In region III, $x/R_{\lambda}=O(1)$ implying that the streamwise ellipticity
is weak and can be negligible (see Ricco & Wu [28], Ricco [27], or Marensi et
al. [24]). Following this assumption, the streamwise distance can be re-scaled
as $x=\bar{x}/R_{\lambda}$ ($x=O(1)$), whereas the cross-stream coordinates
$(y,z)$ are fixed at $O(1)$, i.e $y=\bar{y}=O(1)$ and $z=\bar{z}=O(1)$. The
temporal variation is also ‘re-scaled’ as $t=\bar{t}/R_{\lambda}$. Finally, we
take advantage of the crossflow components of velocity being asymptotically
small compared to the streamwise component, and of the leading order pressure
field being $O(R_{\lambda}^{2})$, which suggests the asymptotic expansions
with $1/R_{\lambda}$ as the small parameters:

$\displaystyle\bar{u}(x,y,z)=u(x,y,z)+...;\hskip
11.38109pt\bar{v}(x,y,z)=R_{\lambda}^{-1}v(x,y,z)+...;\hskip
11.38109pt\bar{w}(x,y,z)=R_{\lambda}^{-1}w(x,y,z)+...;$
$\displaystyle\bar{\rho}(x,y,z)=\rho(x,y,z)+...;\hskip
11.38109pt\bar{p}(x,y,z)=R_{\lambda}^{-2}p(x,y,z)+...;\bar{T}(x,y,z)=T(x,y,z)+...$
(9)

Inserting 2.2 into the full Navier-Stokes equations in curvilinear
coordinates, and retaining the first order terms in the expansions, we obtain
the nonlinear compressible boundary region equations (NCBRE) in the form

$\textbf{V}\cdot\nabla\rho+\rho\nabla\cdot\textbf{V}=0$ (10)
$\rho\textbf{V}\cdot\nabla u=\nabla_{c}\cdot\left(\mu\nabla_{c}u\right)$ (11)
$\rho\textbf{V}\cdot\nabla v+G_{\lambda}u^{2}=-\frac{\partial p}{\partial
y}+\frac{\partial}{\partial y}\left[\frac{2}{3}\mu\left(3\frac{\partial
v}{\partial y}-\nabla\cdot\textbf{V}\right)\right]+\frac{\partial}{\partial
x}\left(\mu\frac{\partial u}{\partial y}\right)+\frac{\partial}{\partial
z}\left[\mu\left(\frac{\partial v}{\partial z}+\frac{\partial w}{\partial
y}\right)\right]$ (12) $\rho\textbf{V}\cdot\nabla w=-\frac{\partial
p}{\partial z}+\frac{\partial}{\partial
z}\left[\frac{2}{3}\mu\left(3\frac{\partial w}{\partial
z}-\nabla\cdot\textbf{V}\right)\right]+\frac{\partial}{\partial
x}\left(\mu\frac{\partial u}{\partial z}\right)+\frac{\partial}{\partial
y}\left[\mu\left(\frac{\partial v}{\partial z}+\frac{\partial w}{\partial
y}\right)\right]$ (13) $\rho\textbf{V}\cdot\nabla
T=\frac{1}{Pr}\nabla_{c}\cdot\left(k\nabla_{c}T\right)+(\gamma-1)M_{\infty}^{2}\mu\left[\left(\frac{\partial
u}{\partial y}\right)^{2}+\left(\frac{\partial u}{\partial
z}\right)^{2}\right]$ (14)

where V is the velocity vector, $\nabla_{c}$ is the crossflow gradient
operator, and $G_{\lambda}u^{2}$ is the term accounting for surface curvature.

$\textbf{V}=u\vec{i}+v\vec{j}+w\vec{k};\hskip
17.07164pt\nabla_{c}=\frac{\partial}{\partial
y}\vec{j}+\frac{\partial}{\partial z}\vec{k}$ (15)

This set of equations is parabolic in the streamwise direction and elliptic in
the spanwise direction. Appropriate initial/upstream and boundary conditions
are necessary to close the problem. The initial/upstream conditions could, for
example, be the same as those used by Ricco & Wu [28] (see also Ricco [27]).
However, for simplicity, in this work we excite the boundary layer using a
small amplitude transpiration velocity ($v_{w}$) at the wall, in the form:

$\displaystyle
v_{w}=A\sin\left[\pi\frac{(x-x_{s})}{(x_{e}-x_{s})}\right]^{2}\cos\left(\pi\frac{z}{\lambda}\right)$
(16)

where $A$ is the amplitude, $x_{s}$ and $x_{e}$ are the start and end
locations of the blowing and suction, respectively (in this study,
$x_{s}=1.5\lambda^{*}$ and $x_{e}=4.5\lambda^{*}$), and $\lambda$ is the
$O(1)$ spanwise separation. The NCBRE were numerically solved using an
algorithm similar to the one employed by Sescu & Thompson [32] in the
incompressible regime. We used second and fourth-order finite-difference
schemes to, respectively, discretize the spatial derivatives in the wall-
normal and spanwise directions. We employed periodic conditions in the
spanwise direction to avoid compromising the stability and a staggered
arrangement in the wall-normal direction to avoid decoupling between the
velocity and pressure. We applied a first-order finite-difference marching
scheme in the streamwise direction and converged the set of equations using a
nonlinear time relaxation method. We found that the same numerical algorithm
applied to the compressible regime is much faster since the continuity
equation is no longer representing the divergence-free condition (making the
set of equations stiff), but it is an equation for density. We placed the
upper boundary far away from the wall ($10$ times the spanwise separation) and
imposed vanishing gradients for all dependent variables.

The mean inflow condition is generated from a similarity solution obtained
using the Dorodnitsyn-Howarth coordinate transformation
$\bar{Y}(x,y)=\int_{0}^{y}\rho(x,\tilde{y})d\tilde{y}$. A similarity variable
is defined as $\eta=\bar{Y}\left(Re_{x}/2x\right)^{1/2},$ where $Re_{x}$ is
the Reynolds number calculated based on the freestream velocity and the
distance from the leading edge. The base velocity and temperature can be
expressed as

$U=F^{\prime}(\eta),\hskip
11.38109ptV=(2xRe_{x})^{-1/2}(\eta_{c}TF^{\prime}-TF),\hskip
11.38109ptT=T(\eta)$ (17)

where prime represents differentiation with respect to $\eta$, and
$\eta_{c}=1/T\int_{0}^{\eta}T(\tilde{\eta})d\tilde{\eta}$. $F$ and $T$ satisfy
the following coupled equations

$\displaystyle\left(\frac{\mu}{T}F^{\prime\prime}\right)^{\prime}+FF^{\prime\prime}=0,$
$\displaystyle\frac{1}{Pr}\left(\frac{\mu}{T}T^{\prime}\right)^{\prime}+FT^{\prime}+(\gamma-1)M^{2}\frac{\mu}{T}F^{\prime\prime
2}=0,$ (18)

subjected to the boundary conditions $F(0)=F^{\prime}(0)=0,\hskip
11.38109ptT^{\prime}(0)=0,\hskip 11.38109ptF^{\prime}\rightarrow 1,\hskip
11.38109ptT\rightarrow 1\rightarrow\hskip 5.69054ptas\hskip
5.69054pt\eta\rightarrow\infty$.Equations (2.2) were solved numerically to
determine $F$ and $T$ and then used in equations (17) to obtained the mean
inflow condition.

##  3 Results and Discussion

We consider Görtler vortices developing in high-speed boundary layers, with
the Mach number ranging from supersonic, $M_{\infty}=2$, to hypersonic
conditions, $M_{\infty}=6$. We neglect chemical reactions inside the boundary
layer at hypersonic speeds. We vary the spanwise separation of the vortices
between $0.3$ cm and $0.7$ cm. We calculate the similarity solution imposed at
the upstream boundary using the Reynolds number based on the distance from the
leading edge, which is maintained constant for a given Mach number (note that
the Reynolds number based on the spanwise separation differs for each
$\lambda$). The radius of curvature of the concave surface is $1.6$ m for all
cases involved in the parametric study. The boundary layer flow is excited by
a small disturbance applied to the vertical velocity imposed at the wall, with
an amplitude of $0.2$% of the freestream velocity (the form of this
disturbance is given in equation (16)). Figure 2 shows a sketch of the flow
domain, with region III highlighted in blue. We impose the similarity type
velocity and temperature profiles for a compressible boundary layer at the
upstream boundary, located on the left hand side of the sketch. At the wall,
we impose the no-slip boundary condition for the velocity and either an
isothermal or adiabatic boundary condition for the temperature field. We
impose vanishing gradients at the top boundary, and a symmetry condition along
the spanwise direction since we only simulate one streamwise vortex
(corresponding to half of the mushroom shape), belonging to the pair of
counter-rotating vortices.

Figure 2:  Flow domain.

We validate the accuracy of our method by comparing our results to the DNS
results of Song et al. [34]. We use a disturbance similar to that employed in
their study with an amplitude of $0.5$% for the transpiration velocity; note,
however, that they imposed the disturbance at the inlet boundary, whereas in
our simulation, we impose it at the wall. In Song et al., the freestream Mach
number is $6.5$ and the radius of curvature is $1.6$ m, the same as in our
parametric study. In figure 3, we plot the scaled amplitude,
$A=\max_{y,z}(T^{\prime})$, based on the temperature disturbance calculated
from each $(y,z)$-plane at every fixed streamwise location ($A_{0}$ in figure
3 is the same amplitude calculated at the $x$-coordinate of the imposed
disturbance). Our results in figure 3 compare very well with the DNS in Song
et al. (note that the streamwise coordinate has been scaled properly to match
the DNS range).

Figure 3: Comparison to DNS results of Song et al. [34].

We focus on the parametric study by varying the Mach number and the spanwise
separation. In figures 4 and 5, we show contour plots of temperature in
consecutive cross-stream planes, for a Mach number $M_{\infty}=3$ and a
spanwise separation $\lambda=0.4$ cm. The contour plots in figure 4 correspond
to the adiabatic wall-condition, while those in figure 5 to the isothermal
condition. They both illustrate the streamwise development of the velocity
magnitude and the temperature field of the Görtler vortices, which display the
mushroom-shape characteristic structures. For both wall conditions, blue and
green colors in the velocity magnitude contours correspond to the low-speed
streaks, while the red color regions are associated with high-speed streaks.
The temperature contour plots of the adiabatic wall-condition indicate that
the low-speed streaks are associated with an increase in temperature (in red
or yellow). The high-speed streaks are characterized by a lower temperature
level (shown here in green). For the isothermal wall-condition case, however,
it appears as if the high temperature region of the structures is sandwiched
between two cooler regions belonging to both low and high speed streaks. This
‘sandwich’ phenomenon leads to a major difference between the two wall
conditions in terms of the height that these mushroom-shape structures take.
In the adiabatic wall-condition case, for example, a significant increase in
temperature at the wall generates more convection, which contributes to the
growth of ‘thermals’ inside the structures. In contrast, the cooler regions at
the wall result in less convection and therefore smaller vortex (i.e mushroom)
structures in the isothermal case.

(a)

(b)

Figure 4: Contour plots of velocity magnitude (a) and temperature (b) in
crossflow planes for $M_{\infty}=3$, $\lambda=0.4$, and adiabatic wall
condition.

(a)

(b)

Figure 5: Contour plots of velocity magnitude (a) and temperature (b) in
crossflow planes for $M_{\infty}=3$, $\lambda=0.4$, and isothermal wall
condition.

Figures 6 and 7 show selected velocity and temperature profiles along the
wall-normal direction at the spanwise location corresponding to the center of
the mushroom-shape structures and the streamwise location corresponding to the
highest vortex energy (see figure 8). In figure 6, we plot the profiles at
different Mach numbers (keeping the spanwise separation constant at
$\lambda=0.4$ cm), whereas in figure 7, we repeat this plot at different
spanwise separations (keeping the Mach number constant at $M_{\infty}=2$). As
we increase the Mach number, the ‘height’ of the streamwise vortices (i.e the
size of the mushroom-shape structures) increases (figs. 6a and 6c) as a result
of the increase in the temperature close to the wall (figs. 6b and 6d). The
increase in size is more evident in the adiabatic case, confirming the
aforementioned observations. For the highest Mach number ($M_{\infty}=6$) the
temperature near the wall increases by a factor of $6$ for the adiabatic case
and $2.5$ in isothermal conditions, compared to the lowest considered Mach
number of $M_{\infty}=2$. It is evident from figure 7 that as we increase the
spanwise separation, the vertical size of the mushroom-shape structures
decreases (almost linearly in fact) for both streamwise velocity and
temperature fields.

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

Figure 6:  Velocity and temperature boundary layer profiles in $z=0$, for
different Mach numbers: (a) velocity, adiabatic wall; (b) temperature,
adiabatic wall; (c) velocity, isothermal wall; (d) temperature, isothermal
wall.

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

Figure 7:  Velocity and temperature boundary layer profiles in $z=0$, for
different spanwise separations: a) velocity, adiabatic wall; b) temperature,
adiabatic wall; c) velocity, isothermal wall; d) temperature, isothermal wall.

We quantify the vortex energy as

$E(x)=\intop_{z_{1}}^{z_{2}}\intop_{0}^{\infty}\left[\left|u(x,y,z)-u_{m}(x,y)\right|^{2}+\left|v(x,y,z)-v_{m}(x,y)\right|^{2}+\left|w(x,y,z)-w_{m}(x,y)\right|^{2}\right]dzdy,$
(19)

where $u_{m}(x,y)$, $v_{m}(x,y)$, and $w_{m}(x,y)$ are the spanwise mean
components of velocity, and $z_{1}$ and $z_{2}$ are the coordinates of the
spanwise domain boundaries.

In figure 8, we plot the vortex energy from equation $(13)$ against the
streamwise coordinate, $x$. For each sub-figure, we fix the Mach number at a
constrant value and vary the spanwise separation, $\lambda$. At all Mach
numbers we find that, as we increase the spanwise separation, there is a
noticeable reduction in the scaled energy for both isothermal (in black) and
adiabatic (in red) wall conditions; the energy saturation location moves
downstream as the spanwise separation increases. The drop in the vortex energy
associated with the increase in the spanwise separation becomes less evident
as the Mach number increases, particularly in the adiabatic wall-condition as
the curves of the different $\lambda$ values coalesce until they almost fall
on top of each other at the highest Mach number (figure 8e). One possible
explanation of this observation is that the reduction in the vortex energy
caused by increasing $\lambda$ diminishes when the convection levels are
increased at the wall, particularly for high Mach numbers (note that as the
Mach number increases, the temperature in proximity to the wall increases
considerably; see figure 6b). Another interesting aspect that we can extract
from figure 8 is the effect that an adiabatic wall condition has on the vortex
energy development: i.e., it appears that the energy growth is delayed in the
adiabatic case, while the peak energy appears to be fixed at the level for the
isothermal case in relatively low supersonic Mach number conditions (figures
8a and 8b). On the other hand, this level is slightly higher at high
supersonic or hypersonic flow conditions (figures 8c, 8d and 8e).

(a) (b) (c)

(d) (e)

Figure 8:  Vortex energy distribution along the streamwise direction for both
isothermal (in black) and adiabatic (in red) wall conditions: a)
$M_{\infty}=2.0$; b) $M_{\infty}=3.0$; c) $M_{\infty}=4.0$; d)
$M_{\infty}=5.0$; e) $M_{\infty}=6.0$

The spanwise averaged wall shear stress evaluated using the integral

$\tau_{w}(x)=\frac{1}{(z_{2}-z_{1})}\intop_{z_{1}}^{z_{2}}\left.\frac{\partial
u}{\partial y}\right|_{y=0}(x,0,z)dz$ (20)

is plotted in figure 9 against the streamwise coordinate for both isothermal
and adiabatic conditions. For each case, we fix the Mach number at a constant
value and vary the spanwise separation. The results show that, in all
considered cases, the wall shear stress increases as the spanwise separation
increases. Also, as expected, the wall shear stress of the adiabatic wall-
condition is lower than that of the isothermal case as a result of the high
level of heating in proximity to the wall (this was observed in previous
studies, such as Spall and Malik [35], Elliot [7], and Sescu et al. [33]
etc.). The jump in the wall shear stress coincides (approximately) with the
same location associated with the the energy saturation initiation point. The
effect of increasing $\lambda$ on the shear stress becomes less apparent in
the case of the highest Mach number (i.e $M_{\infty}=6$) in adiabatic wall-
conditions.

(a) (b) (c)

(d) (e)

Figure 9:  Spanwise averaged wall shear stress distribution along the
streamwise direction: a) $M_{\infty}=2.0$; b) $M_{\infty}=3.0$; c)
$M_{\infty}=4.0$; d) $M_{\infty}=5.0$; e) $M_{\infty}=6.0$; isothermal (in
black) and adiabatic (in red) wall conditions

The spanwise averaged wall heat flux is calculated according to

$q_{w}(x)=-\frac{1}{(z_{2}-z_{1})}\intop_{z_{1}}^{z_{2}}\left.\frac{\partial
T}{\partial y}\right|_{y=0}(x,0,z)dz,$ (21)

which is plotted in figure 10 for the isothermal wall-condition (the wall heat
flux is zero for the adiabatic counterpart) as a function of the streamwise
coordinate. It is evident from this figure that the wall heat flux decreases
as the spanwise separation is increased. Moreover, we notice a decay of the
wall heat flux in the streamwise coordinate range where the energy saturation
takes place.

(a) (b) (c)

(d) (e)

Figure 10:  Spanwise averaged wall heat flux distribution along the streamwise
direction: a) $M_{\infty}=2.0$; b) $M_{\infty}=3.0$; c) $M_{\infty}=4.0$; d)
$M_{\infty}=5.0$; e) $M_{\infty}=6.0$; isothermal (in black) and adiabatic (in
red) wall conditions

The parabolic character of the NCBRE framework allows the solution to be
determined efficiently by a marching procedure in the streamwise direction.
This makes the numerical algorithm very fast, proving it suitable for
parametric studies that can be conducted in a timely manner. In table 1, we
show the CPU time for different cases for a grid resolution of $600$ points in
the streamwise direction, $201$ points along $y$, and $81$ points along $z$;
the spanwise length is limited to the spanwise separation of the vortices,
while the in the wall normal direction the domain size is $5$ times the
spanwise separation, with the grid stretched toward the upper boundary. In
addition, the method does not require large computational resources (in this
paper the simulations were all run on a laptop computer with Intel Core i7
processor).

Table 1: CPU time for different cases Mach number | CPU time (isothermal) | CPU time (adiabatic)
---|---|---
2.0 | 1.9 min | 2.1 min
3.0 | 2.2 min | 2.3 min
4.0 | 2.6 min | 2.7 min
5.0 | 2.7 min | 2.9 min
6.0 | 2.9 min | 3.1 min

##  4 Conclusions

In this paper, Görtler vortices that develop in high-speed boundary layer
flows over concave surfaces are investigated using a numerical solution to the
nonlinear compressible boundary region equations (NCBRE). We targeted the
nonlinear development of the centrifugal instabilities that develop on the
surface, by varying the spanwise separation that of the upstream disturbance
(this dictates the spanwise separation of the downstream Görtler vortices) and
upstream inflow the Mach number. The boundary layer was excited using a small
disturbance at the wall in the form of steady blowing and suction.

We considered a wide range of spanwise separations and the Mach numbers
covering both supersonic and hypersonic regimes. Contours of velocity and
temperature at various crossflow planes showed the vortex development in the
form of mushroom shape structures evolving in the streamwise direction. The
kinetic energy of the primary instability was calculated and plotted against
the streamwise coordinate. It was observed that as the upstream spanwise
separation is increased, the scaled kinetic energy maximum increases, and that
the streamwise location where the energy saturation sets in moves further
downstream. We also calculated the wall shear stress and the wall heat flux
and observed that the wall shear stress increases as the spanwise separation
increases, and - as expected - the wall shear stress is smaller for the
adiabatic condition, as a result of the significant increase in temperature in
proximity to the wall. There was also a jump of the wall shear stress
approximately at the streamwise location corresponding to the point of energy
saturation initiation. For the isothermal wall condition, the wall heat flux
showed a characteristic decrease as the spanwise vortex separation is
increased; there was also decay of the wall heat flux in this case at the
streamwise location corresponding to the point where the energy saturation is
initiated.

The framework of the NCBRE mathematical model is robust and the associated
numerical algorithm provides results in a very short time compared to other
mathematical models such as DNS or PSE, for example. The simulations were run
on a single processor, which makes our BRE approach suitable for multiple
parametric studies.

## References

  * [1] Bagheri, S. and Hanifi A. (2007) The stabilizing effect of streaks on Tollmien-Schlichting and oblique waves: a parametric study, Phys. Fluids, 19, 078103.
  * [2] Chen, X., Huang, G.L., and Lee, C.B. (2019) Hypersonic boundary layer transition on a concave wall: stationary Görtler vortices, J. Fluid Mech., Vol. 865, pp. 1-40.
  * [3] Ciolkosz, L. D. & Spina, E. F. (2006) An experimental study of Görtler vortices in compressible flow. AIAA Paper 4512.
  * [4] Dando, A. H. & Seddougui, S. O. (1993) The compressible Görtler problem in two-dimensional boundary layers. IMA J. Appl. Maths, Vol. 51 (1), 27-67.
  * [5] De Luca, L., Cardone, G., Aymer De La Chevalerie, D. & Fonteneau, A. (1993) Görtler instability of a hypersonic boundary layer. Exp. Fluids, Vol. 16, 10-16.
  * [6] El-Hady N. M. & Verma, A. K. (1983) Growth of Görtler vortices in compressible boundary layers along curved surfaces. J. Eng Appl. Sci., Vol. 2 (3), 213-238.
  * [7] Elliot, J.W. and Bassom, A.P. (2000) The effect of wall cooling on compressible Görtler vortices, Eur. J. Mech. B - Fluids, Vol. 19, pp. 37-68.
  * [8] Görtler, H . (1941) Instabilita-umt laminarer Grenzchichten an Konkaven Wanden gegenber gewissen dreidimensionalen Storungen, ZAMM, Vol. 21, pp. 250–52; english version: NACA Report 1375 (1954)
  * [9] Hall, P. (1982) Taylor-Görtler vortices in fully developed or boundary-layer flows: linear theory, J. Fluid Mech., Vol. 124, pp. 475-494.
  * [10] Hall, P. (1983) The linear development of Görtler vortices in growing boundary layers, J. Fluid Mech., Vol. 130, pp. 41-58.
  * [11] Hall, P. (1989) The nonlinear development of Görtler vortices in growing boundary layers. J. Fluid Mech.. Vol. 193, pp. 243-266.
  * [12] Hall, P. and Horseman. N. (1991) The linear inviscid secondary instability of longitudinal vortex structures in boundary layers. J. Fluid Mech.. Vol. 232, pp. 357-375.
  * [13] Hall, P. & Fu, Y. 1989 On the Görtler vortex instability mechanism at hypersonic speeds. Theor. Comput. Fluid Dyn., Vol. 1 (3), 125-134.
  * [14] Hall, P. & Malik, M. 1989 The growth of Görtler vortices in compressible boundary layers. J. Engng Maths, Vol. 23 (3), 239-251.
  * [15] Kendall, J.M. (1998) Experiments on boundary-layer receptivity to freestream turbulence. AIAA Paper 2004-2335.
  * [16] Kobayashi, R. & Kohama, Y. (1977) Taylor-Görtler instability of compressible boundary layers. AIAA J., Vol. 15 (12), 1723-1727.
  * [17] Landahl, M.T. (1980) A note on an algebraic instability of inviscid parallel shear flows. J. Fluid Mech., Vol. 98, pp. 243-251.
  * [18] Leib, S.J., Wundrow, W. & Goldstein, M. (1999) Effect of freestream turbulence and other vortical disturbances on a laminar boundary layer. J. Fluid Mech., Vol. 380, pp. 169-203.
  * [19] Li, F. and Malik, M. (1995) Fundamental and subharmonic secondary instabilities of Görtler vortices, J. Fluid Mech., Vol. 297, pp. 77-100.
  * [20] Li, F., Choudhari, M., Chang, C.-L., Greene, P., and Wu, M. (2010) Development and Breakdown of Gortler Vortices in High Speed Boundary Layers", 48th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition, Aerospace Sciences Meetings.
  * [21] Li, F., Choudhari, M., and Paredes, P. (2019) Nonlinear Görtler Vortices and Their Secondary Instability in a Hypersonic Boundary Layer, AIAA Paper 2019-3216.
  * [22] Malik, M.R. & Hussaini, M.Y. (1990) Numerical simulation of interactions between Görtler vortices and Tollmien-Schlichting waves. J. Fluid Mech., Vol. 210, pp. 183-199.
  * [23] Marensi, E. & Ricco, P. (2017) Growth and wall-transpiration control of nonlinear unsteady Görtler vortices forced by free-stream vortical disturbances. Physics of Fluids. Vol. 29. pp. 114106.
  * [24] Marensi, E., Ricco, P. & Wu, X. (2017) Nonlinear unsteady streaks engendered by the interaction of freestream vorticity with a compressible boundary layer. J. Fluid Mech.. Vol. 817, pp. 80-121.
  * [25] Matsubara, M. and Alfredsson, P.H. (2001) Disturbance growth in boundary layers subjected to free stream turbulence, J. Fluid Mech., Vol. 430, pp. 149.
  * [26] Ren, J. & Fu, S. (2017) Secondary instabilities of Görtler vortices in high-speed boundary layer flows, J. Fluid Mech., Vol. 781, pp. 388-421.
  * [27] Ricco, P. (2006) Response of a compressible laminar boundary layer to freestream turbulence. PhD thesis, University of London.
  * [28] Ricco, P. & Wu, X. 2007 Response of a compressible laminar boundary layer to freestream vortical disturbances, J. Fluid Mech., Vol. 587, pp. 97-138.
  * [29] Ricco, P., Shah, D. & Hicks, P. D. (2013) Compressible laminar streaks with wall suction. Phys. Fluids 25, 054110.
  * [30] Ricco, P., Tran, D.-L. & Ye, G. (2009) Wall heat transfer effects on Klebanoff modes and Tollmien- Schlichting waves in a compressible boundary layer. Phys. Fluids 21, 024106.
  * [31] Sescu, A., Pendyala, R. and Thompson, D. (2014) On the Growth of Görtler Vortices Excited by Distributed Roughness Elements, AIAA Paper 2014-2885.
  * [32] Sescu, A. and Thompson, D. (2015) On the Excitation of Görtler Vortices by Distributed Roughness Elements, Theoretical and Computational Fluids Dynamics, Vol. 57, pp. 1159-1171.
  * [33] Sescu, A., Alaziz, R., and Afsar, M. (2019) Effect of Wall Transpiration and Heat Transfer on Nonlinear Görtler Vortices in High-speed Boundary Layers, AIAA J. 9, 797-803.
  * [34] Song, R., Zhao, L., and Huang, Z. (2020) Secondary instability of stationary Görtler vortices originating from first/second Mack mode. Phys. Fluids, 32, 034109.
  * [35] Spall, R.E. and Malik, M.R. (1989) Görtler vortices in supersonic and hypersonic boundary layers, Physics of Fluids A: Fluid Dynamics, Vol. 1, pp. 1822.
  * [36] Swearingen, J.D. and Blackwelder, R.F. (1987) The growth and breakdown of streamwise vortices in the presence of a wall. J. Fluid Mech., Vol. 182, pp. 255-290.
  * [37] Viasro, S. and Ricco, P. (2019) Compressible unsteady Görtler vortices subject to freestream vortical disturbances, J. Fluid Mech., Vol. 867, pp. 250-299.
  * [38] Wadey, P. D. (1992) On the linear development of Görtler vortices in compressible boundary layers. Eur. J. Mech. (B/Fluids), Vol. 11, 705-717.
  * [39] Wang, Q.-C., Wang, Z.-G. & Zhao, Y.-X. (2018) Visualization of Görtler vortices in supersonic concave boundary layer. J. Vis., Vol. 21 (1), 57-62.
  * [40] Wu, X., Zhao, D. & Luo, J. (2011) Excitation of steady and unsteady Görtler vortices by freestream vortical disturbances. J. Fluid Mech., Vol. 682, 66-100.
  * [41] Xu, D., Zhang, Y. & Wu, X. (2017) Nonlinear evolution and secondary instability of steady and unsteady Görtler vortices induced by freestream vortical disturbances. J. Fluid Mech., Vol. 829, 681-730.