# Status of the detector setup for the FAMU experiment at RIKEN-RAL for a
precision measurement of the Zemach radius of the proton in muonic hydrogen

R. Rossini,**footnotetext: Corresponding author. A. Adamczak  D. Bakalov  G.
Baldazzi  S. Banfi  M. Baruzzo  R. Benocci  R. Bertoni  M. Bonesini  V.
Bonvicini  H. Cabrera  S. Carsi  D. Cirrincione  M. Clemenza  L. Colace  M. B.
Danailov  P. Danev  A. de Bari  C. de Vecchi  E. Fasci  K. S. Gadedjisso-
Tossou  R. Gaigher  L. Gianfrani  A. D. Hillier  K. Ishida  P. J. C. King  V.
Maggi  A. Menegolli  E. Mocchiutti  L. Moretti  G. Morgante  J. Niemela  C.
Petroselli  A. Pirri  C. Pizzolotto  M.C. Prata  A. Pullia  M. Pullia  R.
Ramponi  M. Rossella  R. Sarkar  A. Sbrizzi  M. Stoilov  J. J. Suarez-Vargas
G. Toci  L. Tortora  E. S. Vallazza  K. Yokoyama  A. Vacchi

###### Abstract

The FAMU experiment at RIKEN-RAL is a muonic atom experiment with the aim to
determine the Zemach radius of the proton by measuring the 1s hyperfine
splitting in muonic hydrogen. The activity of the FAMU Collaboration in the
years 2015-2023 enabled the final optimisation of the detector-target setup as
well as the gas working condition in terms of temperature, pressure and gas
mixture composition. The experiment has started its data taking in July 2023.
The status of the detector setup for the 2023 experimental runs, for the beam
characterisation and muonic X-ray detection in the 100-200 keV energy range,
is presented and discussed.

##  1 Introduction

The first measurements of the 2s-2p Lamb shift in muonic hydrogen ($\mu p$)
returned a value of proton charge radius equal to $0.84087(39)$ fm[1],
7$\sigma$ away from the accepted value from CODATA 2010 $0.8775(51)$ fm[2],
which gave rise to the so-called _proton radius puzzle_. Many efforts are
being made in order to shed light on this puzzle, also in the field of _ep_
scattering[3].

The FAMU (_Fisica degli Atomi MUonici_ \- Muonic Atom Physics) experiment aims
to measure the proton Zemach radius ($r_{Z}$) with an experimental uncertainty
in the order of 1%[4]. The experiment consists in measuring the hyperfine
splitting _hfs_ in 1s muonic hydrogen, which is strongly dependent on
$r_{Z}$[2, 4, 5, 6, 7, 8], with an accuracy better than $10^{-5}$. As of
today, this value has been extracted from _hfs_ in ordinary hyrogen $ep$[7, 8,
9] as well as from the 2s _hfs_ in $\mu p$[1], but not from $\mu p$ _hfs_.

The FAMU experimental methods consists in forming $\mu p$ atoms using a pulsed
low-momentum $\mu^{-}$ beam at the RIKEN-RAL Port1 beamline[10, 11], and
looking for a specific X-ray signature of the excitation of the _hfs_ as a
function of the wavelength from a tunable Mid-InfraRed (MIR) laser. This
signature consists of an excess of delayed muonic hydrogen X-rays in the
100-200 keV region, as explained in Section 2. The FAMU experimental method
yields some strict requirements for the experiment detectors.

The presence of a muon beam hodoscope in the FAMU setup is important for many
reasons. First, the beam shape had to be optimised in order to maximise the
amount of muons entering the hydrogen target. Moreover, the amount of delayed
oxygen X-rays depends on the number of $\mu p$ atoms in the chamber, and
therefore on the muon flux in each muon spill. As this value depends on the
ISIS synchrotron current, which can face small variations in time, muon beam
monitoring is crucial for data normalisation. The FAMU beam hodoscope was
calibrated at a low-rate proton beam at the CNAO synchrotron in Pavia, Italy.
Further details on the beam hodoscope are presented in Section 3.

Now focusing on the X-ray detectors, the FAMU experimental technique requires
a good compromise among time resolution performance, efficiency and energy
resolution at 100-200 keV. As a consequence, it was decided to focus on
inorganic scintillating crystals. The optimisation of timing performances
first led to read out the scintillation light using PhotoMultiplier Tubes
(PMTs). To increase the solid angle coverage with detectors with small-size
readout system, solutions with Silicon PhotoMultipliers (SiPM) were also
studied and implemented, with much effort being made to make the timing
performance comparable with the PMT-read detectors. More details on the FAMU
X-ray detection system and the related R&D can be found in Section 4.

In addition, the presence of a detector with good energy resolution, even with
worse timing performance, is important to always get a check on the X-ray
lines present in the spectra and to cross-check the spectra of the
scintillators. On this purpose, a coaxial High-Purity Germanium (HPGe)
detector has been installed. This detector is described in Section 5.

A sketch of the 2023 FAMU experimental setup is shown in Figure 1, where the
laser, target and detectors position are marked.

Figure 1: CAD drawing of the 2023 FAMU setup. The target is positioned over
the laser optical table, from which the laser is injected in the target
through a periscope. The muon beam enters the target after passing through the
beam hodoscope (muon beam monitor) and a lead collimator. All 35 X-ray
detectors (34 LaBr3 scintillators and 1 HPGe detector) are positioned around
the gas target.  Figure 2: Picture of the 2023 FAMU experiment setup. Left to
right: muon beamline collimator, beam hodoscope, FAMU lead collimator, three
annular detector holders around the target. The big metal boxes are the
detectors read by PMTs, whereas the ones in black plastic-printed cases are
those read by SiPMs. The periscope, the metal stucture entering the target
from the right side, is also visible.

##  2 Muon beamline, laser and target

A $\upmu^{-}$ pulsed beam is delivered to the FAMU experiment through the
RIKEN-RAL facility[10]. All three species of pions are produced by impinging
the 800 MeV proton beam produced by the ISIS synchrotron against a 1 cm thick
graphite target. The decay of negative pions gives rise to a negative muon
beam, which undergoes momentum selection and beam focusing in the RIKEN-RAL
beamline[11]. Eventually, the beam is directed against Port1 and reaches the
FAMU target. Muon beam momentum and focusing optimisation is crucial to
maximise the amount of muonic atoms formed in the target, and therefore to
increase the signal-to-noise ratio. The FAMU experiment working momentum is in
the range between 50 and 60 MeV/c. The beam is delivered with an average rate
of 40 Hz, and each injection consists in two $\sim 50$ ns muon spills
separated by $\sim 320$ ns. The average flux with a momentum of 50-60 MeV/c is
$\sim 7\cdot 10^{4}$ muons per second.

A pulsed MIR laser system with tunable wavelength around 6.78 $\upmu$m has
been developed for the experiment[12]. The laser beam is produced by combining
two laser beams of wavelength 1064 nm and 1262 nm and injecting them into a
Difference-Frequency Generation (DFG) crystal. The laser is capable of
injecting an energy of $\sim 1$ mJ per pulse, with a rate of 25 Hz. This rate
has been set in order to alternate muon spills with/without laser injection,
in order to minimise the contribution of systematics due to drifts in the
target and detector conditions.

The target[4] consists in a $\sim 10\times 10\times 4$ cm3 chamber filled with
7.5 bar of gas mixture (in weight: 98.5% hydrogen, 1.5% oxygen) cooled with
liquid nitrogen (LN2). The cooling system requires both the target chamber and
the LN2 tank to be stored in a cryostat holding $10^{-5}$ bar cryogenic
vacuum. A schematic of the target cryostat is shown in Figure 3. Furthermore,
the target chamber contains a $2.2\times 2.7\times 10$ cm3 optical cavity
aimed at maximising the laser interactions with the gas. A pressure-proof
optical window connects the target to the periscope. This window is the only
thermal-weak point of the cooling system, and a constant flow of dry air is
injected in the periscope to keep the humidity low and avoid condensation on
the window and absorption of the MIR light.

Figure 3: Side section view of the target cryostat (CAD drawing), where the
target chamber, the liquid nitrogen tank and the muon beam direction are
clearly marked. The three detector-holding rings are attached to the target
cryostat around the target chamber (green volume).

##  3 Beam hodoscope

A beam monitor plays a crucial role in allowing data normalisation based on
the muon flux impinging the gas target in each beam spill. In fact, the muon
flux is proportional to the number of muonic atoms and hence the number of
target atoms for the transition.

The FAMU beam monitor consists of a grid of $32\times 32$ polystyrene
scintillating fibres, coated with BaSO4[13, 14]. Each fibre has squared
section, 1 mm pitch, and it is interspaced with the adjacent fibres by 1 mm.
The readout is performed by one end of each fibre with a Hamamatsu S13360-1350
SiPM (active area $1.3\times 1.3$ mm2, pixel pitch 50 $\upmu$m).

As the high instantaneous rate of muons (see Section 2) does not allow single
particle detection, a charge deposit calibration using single particles has
proved to be crucial to determine the instant muon flux during each beam
spill. After an initial test using cosmic rays, the calibration was carried
out at the CNAO synchrotron in Pavia (Italy), with a 50 Hz proton beam with
kinetic energy 150 MeV. This beam energy was chosen from Monte Carlo
simulation in which the energy deposit of negative muons is in the momentum
range of interest for FAMU. Further detail on this procedure is presented in
[15, 16, 17].

##  4 Scintillating detectors

The FAMU experiment is featured with three different kinds of scintillating
detector, placed around the target chamber by three circular holders, in order
to exploit the different features of these crystals. All these detectors are
based on lanthanum bromide with cerium (LaBr3:Ce) scintillating crystals.

The first kind of detector featured in the experiment consists in a set of
cylindrical crystals, having 1" diameter and 1" thickness, read out by
Hamamatsu R11260-200 PMTs, which have a very short rise time ($\sim 12$ ns).
Details on these detectors are to be found in [18]. Six detectors of this kind
are mounted on the central detector holder, as shown in Figure 4 (long gray
boxes).

1" and 1/2" detectors with SiPM readout have also been developed. Crystals for
these detectors have two different geometries: cylinders with 1" diameter and
0.5" thickness (called 1" crystals for simplicity), and cubes with 0.5" side
(1/2" crystals, for the time being). They use Hamamatsu S14161-3050-AS and
S14161-6050-AS SiPM arrays.While 1/2" detectors have a standard parallel
ganging, 1" detectors use a custom 4-1 PCB developed with Nuclear Instruments
srl, see reference [19] for more details. For construction details both of
1/2" and 1" detectors see references [20] and [21]. Twelve 1/2" detectors were
mounted on the downstream detector holder, whereas sixteen 1" detectors were
positioned in the central (6 detectors) and upstream (10 detectors) rings, as
shown in Figure 4. The performances, measured in laboratory, of the 1/2" and
1" detectors are summarised in Table 1[19, 20].

Figure 4: X-ray detector positions in the three detector holders (from the FAMU CAD project). In particular, the front (upstream) ring holds 10 LaBr3 1" crystals read by SiPMs, the central ring holds 6 LaBr3 1" crystals read by SiPMs and 6 LaBr3 1" crystals read by PMTs, whereas the back (downstream) ring holds 12 LaBr3 smaller (1/2") crystals read by SiPMs and it also hosts the HPGe detector. The space left with no detectors on the right side of the central and front rings is occupied by the periscope, as one can see in Figure 1.  Table 1: Time and energy performance of SiPM-based detectors determined during laboratory tests with radioactive sources, focusing on the 662 keV peak from 137Cs and the 122 keV peak from 57Co. | rise time (ns) | fall time (ns) | Res. (%) @ 137Cs | Res. (%) @ 57Co
---|---|---|---|---
1" | $29.3\pm 1.5$ | $147\pm 13$ | $3.01\pm 0.16$ | $7.9\pm 0.4$
1/2" | $43\pm 5$ | $372\pm 17$ | $3.27\pm 0.11$ | $8.4\pm 0.6$

Energy resolutions on the beamline, obtained from the reconstructed muonic Ag
peak at 141 keV for detectors in a typical run on the full experiment setup
with 55 MeV/c impinging muons, is reported in Figure 5.

Figure 5: Distribution of FWHM energy resolutions at the muonic silver peak
corresponding to the 5f-5g transition (141 keV, theoretically calculated).
Preliminary plot: three detectors are still under investigation and were not
included; sample size: . The standard boxplot rules hold: a set of box,
errorbar and orange line mark the four quartiles in each dataset (i.e. the
orange line corrsponds to the median value).

##  5 Germanium detector

A commercial ORTEC GEM-S5020P4 High-Purity Germanium (HPGe) Detector has been
installed in the FAMU setup along with the scintillating detectors. The
detector core is is a p-type coaxial germanium crystal (length 20.8 mm,
diameter 30.4 mm) with a 0.9 mm carbon entrance window. The detector is
supplied with 1.5 kV positive bias and the signals are treated with an ORTEC
A257P preamplifier and an ORTEC 672 gaussian amplifier, with 3 $\upmu$s
shaping time.

From the same measurements as the ones presented in the previous section, the
FWHM energy resolution for the HPGe detector at the 141 keV muonic silver peak
is $\mathbf{(1.26\pm 0.17)\%}$. Although the shaping time and the typical
signal length of HPGe detectors do not allow to extract timing information in
the scale of hundreds of nanoseconds, the good energy resolution is a key
factor in the target monitoring and in the beam optimisation for the
experiment.

##  6 Summary

The FAMU experiment has started its data taking at the RIKEN-RAL Port1
$\upmu^{-}$ beam facility in mid 2023. The X-ray detector setup, presented in
detail in this paper, meets the detector requirements using a set of 35
detectors, in particular:

  * •

1 High-Purity Germanium detector, mainly for spectra and target check;

  * •

34 LaBr3 scintillators, of which:

    * –

6 cylindrical 1" diameter $\times$ 1" thickness crystals read by PMTs, all in
the central detector holder, which have the best time performance;

    * –

12 cylindrical 0.5" diameter $\times$ 0.5" thickness crystals read by SiPMs,
all in the downstream detector holder, with the best energy resolution;

    * –

16 cylindrical 1" diameter $\times$ 0.5" thickness crystals read by SiPMs, all
in the downstream detector holder, with the best energy resolution together
with the previous ones.

## References

  * [1] A. Antognini, F. Nez, K. Schuhmann et al., _Proton structure from the measurement of 2s-2p transition frequencies of muonic hydrogen_ , _Science_ 339 (2013) 417.
  * [2] C. Carlson, _The proton radius puzzle_ , _Progress in Particle and Nuclear Physics_ 82 (2015) 59.
  * [3] W. Xiong, A. Gasparian, H. Gao et al., _A small proton charge radius from an electron–proton scattering experiment_ , _Nature_ 575 (2019) 147.
  * [4] C. Pizzolotto, A. Adamczak, F. Bakalov, D. et al., _The famu experiment: muonic hydrogen high precision spectroscopy studies_ , _Eur. Phys. J. A_ 56 (2020) 185.
  * [5] C. Carlson, V. Nazaryan and K. Griffioen, _Proton-structure corrections to hyperfine splitting in muonic hydrogen_ , _Phys. Rev. A_ 83 (2011) 042509.
  * [6] R. Faustov, A. Martynenko, M. G.A. and V. Sorokin, _Radiative nonrecoil nuclear finite size corrections of order $\alpha$(z$\alpha$)5 to the hyperfine splitting of s-states in muonic hydrogen_, _Physics Letters B_ 733 (2014) 354.
  * [7] A. Dupays, A. Beswick, B. Lepetit, C. Rizzo and D. Bakalov, _Proton zemach radius from measurements of the hyperfine splitting of hydrogen and muonic hydrogen_ , _Phys. Rev. A_ 68 (2003) 052503.
  * [8] A. Volotka, V. Shabaev, G. Plunien et al., _Zemach and magnetic radius of the proton from the hyperfine splitting in hydrogen_ , _Eur. Phys. J. D_ 33 (2005) 23.
  * [9] S. Brodsky, C. Carlson, J. Hiller and D. Hwang, _Constraints on proton structure from precision atomic-physics measurements_ , _Phys. Rev. Lett._ 94 (2005) 022001.
  * [10] T. Matsuzaki, K. Ishida, K. Nagamine et al., _The riken-ral pulsed muon facility_ , _NIM A_ 465 (2001) 365.
  * [11] A.D. Hillier, J.S. Lord, K. Ishida and C. Rogers, _Muons at isis_ , _Philosophical Transactions of the Royal Society A_ 377 (2019) 20180064.
  * [12] M. Baruzzo, J. Suáres-Vargas, H. Cabrera et al., _A mid-ir laser source for muonic hydrogen spectroscopy: The famu laser system_ , _Optics & Laser Technology, SSRN_ (2023) .
  * [13] M. Bonesini, R. Bertoni, F. Chignoli et al., _The construction of the fiber-sipm beam monitor system of the r484 and r582 experiments at the riken-ral muon facility_ , _Journal of Instrumentation (JINST)_ 12 (2017) C03035.
  * [14] M. Bonesini, R. Benocci, R. Bertoni et al., _The upgraded beam monitor system of the famu experiment at riken–ral_ , _NIM A_ 936 (2019) 592.
  * [15] R. Rossini, R. Benocci, R. Bertoni et al., _Characterisation of muon and proton beam monitors based on scintillating fibres with a sipm read-out_ , _NIM A_ 1046 (2023) 167684.
  * [16] R. Rossini, R. Benocci, R. Bertoni et al., _Characterisation of a scintillating fibre-based hodoscope exposed to the cnao low-energy proton beam_ , _NIM A_ 1046 (2023) 167746.
  * [17] R. Rossini, R. Benocci, R. Bertoni et al., _Characterisation of a low-momentum high-rate muon beam monitor for the famu experiment at the cnao-xpr beam facility_ , _submitted to JINST as another proceeding for the IPRD23 Conference_ (2023) .
  * [18] G. Baldazzi, A. Vacchi, C. Labanti et al., _A small proton charge radius from an electron–proton scattering experiment_ , _Journal of Instrumentation (JINST)_ 12 (2017) C03067.
  * [19] M. Bonesini, R. Bertoni, A. Abba et al., _Improving the time resolution of large area labr3:ce detectors with sipm array readout_ , _submitted to Condensed Matter, MDPI_ (2023) .
  * [20] M. Bonesini, R. Benocci, R. Bertoni et al., _Large area labr3:ce crystals read by sipm arrays with improved timing and temperature gain drift control_ , _NIM A_ 1046 (2023) 167677.
  * [21] M. Bonesini, R. Bertoni, R. Benocci et al., _One inch labr3:ce detectors, with temperature control and improved time resolution for low energy x-rays spectroscopy_ , _Proceedings of Science_ 398 (2022) 770.


# Interpretable Data-driven Methods for Subgrid-scale Closure in
LES for Transcritical LOX/GCH4 Combustion

Wai Tong Chung<EMAIL_ADDRESS>Aashwin Ananda Mishra  Matthias Ihme
Department of Mechanical Engineering, Stanford University, Stanford, CA 94305,
USA  SLAC National Accelerator Laboratory, Menlo Park, CA 94025, USA

###### Abstract

Many practical combustion systems such as those in rockets, gas turbines, and
internal combustion engines operate under high pressures that surpass the
thermodynamic critical limit of fuel-oxidizer mixtures. These conditions
require the consideration of complex fluid behaviors that pose challenges for
numerical simulations, casting doubts on the validity of existing subgrid-
scale (SGS) models in large-eddy simulations of these systems. While data-
driven methods have shown high accuracy as closure models in simulations of
turbulent flames, these models are often criticized for lack of physical
interpretability, wherein they provide answers but no insight into their
underlying rationale. The objective of this study is to assess SGS stress
models from conventional physics-driven approaches and an interpretable
machine learning algorithm, i.e., the random forest regressor, in a turbulent
transcritical non-premixed flame. To this end, direct numerical simulations
(DNS) of transcritical liquid-oxygen/gaseous-methane (LOX/GCH4) inert and
reacting flows are performed. Using this data, a priori analysis is performed
on the Favre-filtered DNS data to examine the accuracy of physics-based and
random forest SGS-models under these conditions. SGS stresses calculated with
the gradient model show good agreement with the exact terms extracted from
filtered DNS. The accuracy of the random-forest regressor decreased when
physics-based constraints are applied to the feature set. Results demonstrate
that random forests can perform as effectively as algebraic models when
modeling subgrid stresses, only when trained on a sufficiently representative
database. The employment of random forest feature importance score is shown to
provide insight into discovering subgrid-scale stresses through sparse
regression.

###### keywords:

Turbulence modeling , Direct numerical simulation , Transcritical combustion ,
Random forests , Machine learning , Explainable Artificial Intelligence

††journal:

##  1 Introduction

The development of accurate computational tools has been essential for
studying non-premixed flames within high pressure combustors. Large-eddy
simulations (LES) provide a feasible computational approach in capturing the
behavior of flows within practical combustors. However, high pressure
combustors in rockets and diesel engines operate under conditions that exceed
the thermodynamic critical limits of both fuel and oxidizer. Consequently,
these conditions generate trans- and supercritical flows – with complex
behaviors that pose challenges for numerical modeling and simulations [1]. In
many simulations, Lagrangian droplet methods are typically employed for
simulating, which assumes the presence of liquid and gas phases. Simulations
employing the Lagrangian droplet method have been shown to have good agreement
with experimental measurements [2, 3]. However, these models often involve
careful selection of the breakup and evaporation models along with parameter
tuning.

Another approach for investigating high-pressure flows involves the use of the
diffuse-interface method. In contrast to sharp interface techniques, where
interfaces are explicitly tracked or resolved in the computational domain,
this method artificially diffuses the interface and treats the entire flow
with a single real-fluid state equation. While LES incorporating real-fluids
effects have successfully been employed to simulate transcritical combustion
[4, 5, 6], many of these simulations employ existing subgrid-scale (SGS)
models that have been typically developed for applications in subcritical
pressure conditions [7, 8]. As a consequence, the application of these models
to non-ideal flow regimes introduces uncertainties. One method for evaluating
SGS models involves _a priori_ analysis, where modeled SGS terms are compared
with exact unclosed terms extracted from filtered DNS. Selle et al. [9]
performed a priori analysis on a three-dimensional DNS database of
supercritical binary mixtures in turbulent mixing layers to demonstrate that
the Smagorinsky model [10] performed poorly when predicting SGS stresses,
while the gradient [8] and scale-similar [11] models performed well. In the
same work, previously neglected unclosed terms for pressure and heat flux were
shown to be essential under supercritical conditions. Unnikrishnan et al. [12]
performed a priori analysis on two-dimensional DNS of a transcritical reacting
liquid-oxygen/gaseous-methane (LOX/GCH4) mixture to demonstrate that the mixed
SGS model incorporating the dynamic Smagorinsky [13] was three times more
accurate than the sole use of the dynamic Smagorinsky.

One approach for developing closure models in turbulent reacting flows involve
the use of data-driven methods. _A priori_ studies have been performed to
demonstrate that convolutional neural-networks can provide accurate closure
for turbulent combustion models [14, 15]. Henry de Frahan et al. [16]
demonstrated that deep learning models can generate as accurate results as
random forests, another regression model, with 25 fold improvement in
computational costs when predicting the sub-filter probability density
function in another _a priori_ study. Ranade and Echekki [17] conducted an _a
posteriori_ study to show that deep learning models can be trained with
experimental data to generate closure models for chemical scalars in Reynolds-
averaged Navier–Stokes (RANS) simulations of turbulent jet flames. These
studies have demonstrated the ability of deep learning models in predicting
closure models, but tend to offer little physical understanding.

Computational studies of high-pressure non-premixed flames were pioneered by
Bray et. al. [2, 3, 18]. One work [18] examined the effect of different
Damköhler numbers (Da) on autoignition in high-pressure non-premixed flames
under decaying homogeneous isotropic turbulence. In the present work, we
perform DNS calculations of inert and reacting LOX/GCH4 non-premixed mixtures
in the presence of decaying turbulence, under different Da, in order to
evaluate algebraic and data-driven models for predicting unclosed SGS stresses
under high pressure. Within this context, this _a priori_ study has the
following objectives:

  * 1.

To identify and quantify limitations of conventional algebraic SGS models in
transcritical flows.

  * 2.

To evaluate the insight gained from the random forest, an interpretable
machine learning algorithm, in constructing a data-informed, interpretable SGS
model.

The mathematical models for simulating the turbulent transcritical flows are
presented in Section 2. Details regarding the DNS configuration are discussed
in Section 3. Section 4 describes the SGS models and data-driven methods
employed in the present work. Results from this _a priori_ study are discussed
in Section 5, before offering concluding remarks in Section 6.

##  2 Mathematical Models

The governing equations that are solved in the present study are the
conservation equations for mass, momentum, energy, and chemical species:

$\displaystyle\partial_{t}{\rho}+\nabla\cdot({\rho}{\boldsymbol{u}})$
$\displaystyle=0$ (1a)
$\displaystyle\partial_{t}({\rho}{\boldsymbol{u}})+\nabla\cdot({\rho}{\boldsymbol{u}}{\boldsymbol{u}})$
$\displaystyle=-\nabla\cdot({p}\boldsymbol{I})+\nabla\cdot{\boldsymbol{\tau}}_{v}$
(1b)
$\displaystyle\partial_{t}({\rho}{e_{t}})+\nabla\cdot[{\boldsymbol{u}}({\rho}{e_{t}}+{p})]$
$\displaystyle=-\nabla\cdot{\boldsymbol{q}}_{v}+\nabla\cdot[({\boldsymbol{\tau}}_{v})\cdot{\boldsymbol{u}}]$
(1c)
$\displaystyle\partial_{t}({\rho}Y_{k})+\nabla\cdot({\rho}{\boldsymbol{u}}Y_{k})$
$\displaystyle=-\nabla\cdot\boldsymbol{{j}}_{v}+{\dot{\omega}}_{k}\quad\textrm{where}\quad
k=1,2,...,N_{s}-1$ (1d)

with density $\rho$, velocity vector $\boldsymbol{u}$, pressure $p$, specific
total energy $e_{t}$, stress tensor $\boldsymbol{\tau}$, and heat flux
$\boldsymbol{q}$. $Y_{k}$, $\boldsymbol{j}_{k}$, and ${\dot{\omega}_{k}}$ are
the mass fraction, diffusion flux, and source term for the $k^{\textrm{th}}$
species, while subscript $v$ denotes viscous quantities.

In the a priori analysis carried out in this study, a top-hat filter $H$ with
a desired filter size $\overline{\Delta}$ is applied on a arbitrary quantity
$\phi$ from the DNS data through a volume integral:

$\displaystyle\overline{\phi(\boldsymbol{x})}$
$\displaystyle=\int_{V}\phi(\boldsymbol{x})H(\boldsymbol{x}-\boldsymbol{y},{\overline{\Delta}})d\boldsymbol{y}$
(2a) $\displaystyle H(\boldsymbol{x}-\boldsymbol{y},{\overline{\Delta}})$
$\displaystyle=\begin{cases}\frac{1}{\overline{\Delta}}&\textrm{for
}|\boldsymbol{x}-\boldsymbol{y}|\leq\overline{\Delta}/2,\\\ 0&\text{otherwise.
}\\\ \end{cases}$ (2b)

and Favre-averaged quantity:

$\widetilde{\phi}=\frac{\overline{\rho\phi}}{\overline{\rho}}$ (3)

where ${\overline{\,\cdot\,}}$ denotes a filtered quantity and
$\widetilde{\,\cdot\,}$ is a Favre-filtered quantity. After filtering, the
governing equations become:

$\displaystyle\partial_{t}{\overline{\rho}}+\nabla\cdot({\overline{\rho}}\widetilde{\boldsymbol{u}})$
$\displaystyle=0$ (4a)
$\displaystyle\partial_{t}({\overline{\rho}}\widetilde{\boldsymbol{u}})+\nabla\cdot({\overline{\rho}}\widetilde{\boldsymbol{u}}\widetilde{\boldsymbol{u}})$
$\displaystyle=-\nabla\cdot({\overline{p}}\boldsymbol{I})+\nabla\cdot({\overline{\boldsymbol{\tau}}}_{v}+{\boldsymbol{\tau}}^{sgs})$
(4b)
$\displaystyle\partial_{t}({\overline{\rho}}\widetilde{e}_{t})+\nabla\cdot[\widetilde{\boldsymbol{u}}({\overline{\rho}}\widetilde{e}_{t}+{\overline{p}})]$
$\displaystyle=-\nabla\cdot({\overline{\boldsymbol{q}}}_{v}+{\boldsymbol{q}}^{sgs})+\nabla\cdot[({\overline{\boldsymbol{\tau}}}_{v}+{\boldsymbol{\tau}}^{sgs})\cdot\widetilde{\boldsymbol{u}}]$
(4c)
$\displaystyle\partial_{t}({\overline{\rho}}\widetilde{{Y}}_{k})+\nabla\cdot({\overline{\rho}}\widetilde{\boldsymbol{u}}\widetilde{{Y}}_{k})$
$\displaystyle=-\nabla\cdot(\boldsymbol{{\overline{j}}}_{v}+\boldsymbol{{j}}^{sgs})+{\overline{\dot{{\omega}}}}_{k}\quad\textrm{where}\quad
k=1,2,...,N_{s}-1$ (4d)

with superscript $sgs$ denoting subgrid-scale quantities. Exact subgrid-scale
quantities can then be extracted directly from the filtered DNS or
approximated through the models described in Section 4.

The Peng-Robinson cubic state equation [19] is employed to model real-fluid
thermodynamics under transcritical conditions:

$p=\frac{\rho RT}{1-b\rho}-\frac{a\rho^{2}}{1+2b\rho-b^{2}\rho^{2}}$ (5)

with mixture-specific gas constant $R$. The coefficients $a$ and $b$ account
for effects of intermolecular forces and volumetric displacement, and are
dependent on temperature and composition [20]. Details regarding the
evaluation of specific heat capacity, internal energy, and partial enthalpy
from the Peng-Robinson state equation is described in Ma et al. [21]

In this study, the two-step five-species CH4-BFER mechanism [22] is employed,
which has been applied to investigate a supercritical gas-turbine combustor at
20 MPa [23]. In DNS of trans- and supercritical combustion, reduced chemical
mechanisms [24, 25] have been employed to circumvent large computational costs
incurred by solving non-ideal state equations. Takahashi’s high-pressure
correction [26] is used to evaluate the binary diffusion coefficients. Since
only two species are used in the inert simulations, the binary diffusion
coefficients are exact. Thermal conductivity and dynamic viscosity are
evaluated using Chung’s method with high-pressure correction [27]. For multi-
species mixtures in the reacting cases, Chung’s pressure correction is known
to produce oscillations, especially for dynamic viscosity. Hence, transport
properties of the mixture are evaluated through mole-fraction-averaging of
individual species. A similar approach has been applied in prior studies [5,
12].

Simulations are performed by employing an unstructured compressible finite-
volume solver [28, 21, 29]. A central scheme, which is 4th-order accurate on
uniform meshes, is used along with a 2nd-order ENO scheme. The ENO scheme is
activated only in regions of high local density variations using a threshold-
based sensor to describe sharp interfaces present in transcritical flows. Due
to the density gradients present at trans- and supercritical conditions, an
entropy-stable flux correction technique [21] is used to dampen non-linear
instabilities in the numerical scheme. The double-flux method by Ma et al.
[21] is used with a dynamic sensor to eliminate spurious pressure
oscillations. A Strang-splitting scheme is employed for time-advancement,
combining a strong stability preserving 3rd-order Runge-Kutta (SSP-RK3) scheme
for integrating the non-stiff operators with a semi-implicit scheme [30] for
advancing the chemical source terms.

##  3 DNS Configuration

Inert and reacting direct numerical simulations are performed on a three-
dimensional cubic domain, with length $L$, of LOX/GCH4 mixture shown in Fig.
1. In this setup, a spherical liquid oxygen core, with a radius $r=0.25L$, is
initialized in gaseous methane, where the radial profile of the initial
conditions is chosen to match one-dimensional calculations of inert and
reacting mixtures under the same conditions. Fuel and oxidizer temperatures,
$T_{\textrm{CH4}}$ and $T_{\textrm{O2}}$ have been set at 120 K and 300 K,
respectively, while the pressure is set at 10 MPa. These operating conditions
are chosen to match practical LOX/GCH4 combustors, and were investigated in
previous studies [31, 32]. Periodic boundary conditions are used for all
boundaries in the inert case. In the reacting cases, non-reflecting pressure
outlets are used in both boundaries in the $x$-direction, while remaining
boundaries are periodic.

Figure 1: DNS investigated at initial time $t=0$ and one eddy turnover time
$t=t_{I}$. Isosurface shows stoichiometric mixture fraction $Z=0.2$ for the
inert case.

The initial velocity profile was generated with a synthetic isotropic
turbulence generator by Saad et al. [33] with zero mean velocity, based on the
von Kármán-Pao energy spectrum:

$\displaystyle E(\kappa)$ $\displaystyle=\alpha\frac{u^{\prime
2}}{\kappa_{I}}\frac{(\kappa/\kappa_{I})^{4}}{[1+(\kappa/\kappa_{I}]^{17/6}}\exp\left[-2\left(\frac{\kappa}{\kappa_{\eta}}\right)^{2}\right],$
(6a) $\displaystyle\alpha$ $\displaystyle=1.453,$ (6b)
$\displaystyle\kappa_{I}$ $\displaystyle=0.746834/l_{I},$ (6c)

where $u^{\prime}$ is the RMS velocity, $\kappa$ is the wave number, and
$\kappa_{\eta}$ the Kolmogorov wave number. The chosen scaling constant
$\alpha$ and large-eddy wavenumber $\kappa_{I}$ are typical for isotropic
turbulence [34]. In all cases, the integral lengthscale $l_{I}$ and root-mean-
squared (RMS) velocity fluctuation $u^{\prime}$ have been chosen to produce a
turbulent Reynolds number Ret of 80, which has been computed with the averaged
kinematic viscosity of oxygen and methane at 120 K and 300 K, respectively.

In the reacting cases, two different Damköhler numbers, Da, of 780 and 10 are
investigated, corresponding to flamelet and unsteady regimes [35],
respectively. The Damköhler number is given by the ratio of physical timescale
$t_{phys}$ and chemical timescale $t_{chem}$:

$\textrm{Da}=\frac{t_{phys}}{t_{chem}},$ (7)

where $t_{chem}$ is approximated from the extinction strain rate of a one-
dimensional counterflow diffusion flame of a LOX/GCH4 mixture under similar
conditions, and physical time is evaluated from the eddy turnover time
$t_{phys}=t_{I}$. Figure 2(a) shows that the mean temperature $\langle
T\rangle$ is lower when Da = 10 than when Da = 780, due the presence of local
extinction. This is also reflected in Figure 2(b) where the consumption of CH4
is slower in the case Da = 10 than the case Da = 780. This decrease in
temperature and composition also results in a slower decay of turbulence, as
shown by the mean turbulent kinetic energy $\langle\textrm{TKE}\rangle$,
normalized by the initial TKE, shown in Fig. 2(c).

(a) Temperature.

(b) Mass fraction.

(c) Turbulent kinetic energy.

Figure 2:  Temporal evolution of global temperature $T$, mass fraction
$Y_{k}$, and normalized turbulent kinetic energy TKE for two reacting cases..

In this study, analysis is performed on all cases after
$t=\textrm{argmax}(t_{I},t_{chem})$, which is typically done for DNS of
combustion under decaying turbulence in order to ensure the flow fields are
independent of initialization [36]. Instantaneous flow fields for axial
velocity component $u_{1}$, mixture fraction $Z$, and mixture-fraction
conditioned temperature $T$ for the reacting cases at $t=0$ and $t=t_{I}$ are
shown in Fig. 3.

(a) Case Da = 10 at initial time $t=0$.

(b) Case Da = 10 at eddy turnover time $t=t_{I}$.

(c) Case Da = 780 at eddy turnover time $t=t_{I}$.

Figure 3:  Axial velocity $u_{1}$, mixture fraction $Z$, and conditional
temperature $T$ for the reacting cases at transverse location $z=0$.

Table 1 summarizes the DNS cases examined in this study. The domain lengths in
all direction were chosen to be eight times the size of the integral
lengthscale $l_{I}$ to minimize effects of the boundary conditions. Cell size
$\Delta$ is prescribed in the order of the Kolmogorov lengthscale, ensuring
that all lengthscales are resolved. In addition, a mesh refinement study was
performed, where the energy spectra of velocity was found to converge between
$128^{3}$ and $256^{3}$. Simulations for all three cases are advanced with a
convective CFL number of unity, corresponding to timesteps of 2.5 and 0.5 ns
for cases Da = 780 and Da = 10, respectively. The simulations were performed
using 960 Intel Xeon (E5-2698 v3) processors, and 2.3 $\mu$s and 0.6 $\mu$s of
physical time could be completed in about an hour wall clock time for cases Da
= 780 and Da = 10, respectively.

Table 1: Summary of DNS cases.

Case | $N_{x,y,z}$ |  $L_{x,y,z}$ [$\mu$m] | Ret |  $l_{I}$ [$\mu$m] |  $\eta_{k}$ [$\mu$m] |  $\Delta$ [$\mu$m] |  $t_{I}$ [$\mu$s]
---|---|---|---|---|---|---|---
Inert | $128$ | 500 | 80 | 62.5 | 2.32 | 3.91 | 286
Da = 780 | $128$ | 500 | 80 | 62.5 | 2.32 | 3.91 | 286
Da = 10 | $128$ | 60 | 80 | 7.50 | 0.278 | 0.469 | 4.12

##  4 Subgrid-scale models and data-driven methods

In this study, we employ conventional algebraic and data-driven methods for
predicting the subgrid-scale fluxes from the LES momentum equation (Eq. 4c):

$\tau_{ij}^{sgs}={\overline{\rho}}(\widetilde{u_{i}u_{j}}-\widetilde{u}_{i}\widetilde{u}_{j})$
(8)

Two algebraic SGS models, namely the Vreman and the gradient model, and random
forest regressors are evaluated. Additionally, we demonstrate the employment
of random forest feature importance scores for assisting the discovery of
algebraic SGS models by sparse regression.

###  4.1 Algebraic models

The Vreman SGS model [7] is derived from the eddy-viscosity hypothesis:

$\tau_{ij}^{sgs}\simeq-2{\overline{\rho}}\nu_{\text{SGS}}\widetilde{S}_{ij}+\frac{1}{3}\tau_{kk}\delta_{ij},$
(9)

where $S_{ij}$ is the velocity strain tensor, and $\delta_{ij}$ is the
Kronecker delta. The eddy viscosity $\nu_{\text{SGS}}$ is evaluated for a
filter width ${\overline{\Delta}}$ as follows:

$\displaystyle\nu_{\text{SGS}}$
$\displaystyle=C_{v}\sqrt{\frac{B}{a_{ij}a_{ij}}},$ (10a) $\displaystyle
a_{ij}$ $\displaystyle=\dfrac{\partial\widetilde{u}_{i}}{\partial x_{j}},$
(10b) $\displaystyle B$
$\displaystyle=\beta_{11}\beta_{22}-\beta_{12}^{2}+\beta_{11}\beta_{33}-\beta_{13}^{2}+\beta_{22}\beta_{33}-\beta_{23}^{2},$
(10c) $\displaystyle\beta_{ij}$
$\displaystyle={\overline{\Delta}}^{2}a_{ki}a_{kj},$ (10d)

where a Vreman coefficient $C_{v}$ of 0.07 is typically used in isotropic
turbulence.

The gradient model by Clark et al. [8] is extracted from the first term in the
Taylor series expansion of the filtering operation, and is given by:

$\tau_{ij}^{sgs}\approx{\overline{\rho}}C_{g}{\overline{\Delta}}^{2}\dfrac{\partial\widetilde{u}_{i}}{\partial
x_{k}}\dfrac{\partial\widetilde{u}_{j}}{\partial x_{k}},$ (11)

where a coefficient $C_{g}$ of 1/12 is typically used when a top-hat filter is
employed. In the present study, we wil evaluate both models and compare
results against DNS data and a data-driven approach.

###  4.2 Random forest regressor

In this study, we employ the random forest as our regression algorithm. Random
forests [37] consist of an ensemble of decorrelated Classification and
Regression Trees (CARTs) [38]. CARTs are a machine learning approach for
formulating prediction models from data by recursively partitioning the
inputted feature space, and fitting a simple prediction within each final
partition. The partitioning of the feature space can be represented as a
decision trees. Decision trees are supervised graph based model wherein the
tree consists of nodes and edges. The internal (or non-terminal) nodes of the
tree represent splits based on learned partitions of the feature space. Each
leaf (or terminal) node is associated with a numerical value for regression
trees (as opposed to categorical targets for classification trees).

During the training phase, the structure of the decision tree and the
partitions associated with each node are inferred. During each step of this
phase, exhaustive sets of splits over different input features are evaluated.
The split leading to maximal decrease in prediction variance is selected at
the associated node. The procedure continues to make recursive splits based on
dataset until it has reduced the overall variance below a given threshold or
upon reaching a given stopping parameter (for instance, upon reaching a
maximal depth of the tree).

During the prediction phase, a new sample is traversed down the tree from the
root node to a leaf node, wherein its path is determined based on the
partition at each internal node. Once a leaf node is reached, the numerical
value associated with the specific leaf node is outputted as the prediction
for the sample.

Decision trees are non-parametric and can model arbitrarily complex relations
without any _a priori_ assumptions, but are prone to overfitting due to the
greedy nature of the inference algorithm. Thus, decision trees have low bias
but high variance.

Random Forests are an ensemble learning algorithm, wherein the predictions of
ensembles of decorrelated decision trees are aggregated so as to give a final
meta-model with low bias and low variance. The decorrelation amongst the
individual trees in the ensemble is achieved using bootstrapping [37] in
conjunction with feature bagging [39] during the training of each decision
tree. The final prediction of the resulting ensemble model is by averaging the
predictions of all trained individual trees (or equivalently, via
aggregation).

While accurate predictions are an important goal for machine learning models,
in many fields of application it is just as important to derive understanding
from the model. At the basic level, this can be embodied via feature
importances, wherein the trained model also provides information regarding
importances of different input features to the final prediction. Such measures
of model interpretability provide insight into the underlying rationale
learned by the model during training and can lead to high confidence in the
model. Similarly, such interpretability measures can lead to data-driven
discovery of new relationships between the input features and targets. Random
forests provide the Mean Decrease Impurity (MDI) importance measure for all
the features in the input set [40]. Here, the importance of a feature is given
by aggregating the weighted decrease in variance for all the nodes where the
specific feature is used as the criterion for partitioning the feature space.

In the present investigation, the random forest regressor implementation from
the Scikit-learn library [41] is used. Here, a random forest consisting of
fifty decision trees is employed. The hyperparameters of the random forest are
determined using a random grid search approach with a 3-fold cross-validation
set. Training is performed once _a priori_ , and requires 88s of walltime with
8 CPUs, when trained on data coarsened for three different filter sizes from a
single timestep. Prediction time for a $64^{3}$ dataset requires 2.4 s on a
single CPU.

###  4.3 Sparse regression for model discovery

A linear model $\hat{f}_{i}$ for $m$ number of samples is typically expressed
as the weighted sum of independent quantities $X$:

$\hat{f}_{i}=\sum_{j=1}^{n}\beta_{j}X_{ij}\quad 0\leq i\leq m$ (12)

with $n$ number of model coefficients $\beta$ being employed.

When samples of the ground truth f are available, the model coefficients can
be found with the $l_{1}$-norm regularized least-squares or _lasso_ method
[42]:

$\displaystyle\min_{\boldsymbol{\beta}}\left\\{\frac{1}{m}||\mathbf{f}-\mathbf{X}\boldsymbol{\beta}||^{2}_{2}+\lambda||\boldsymbol{\beta}||_{1}\right\\},$
(13)

where $\lambda$ is a regularization parameter for controlling the tradeoff
between the least-squares fit and the $l_{1}$-norm. Since, the optimization
scheme in Eq. 13 also minimizes the $l_{1}$-norm of the model coefficients
$||\boldsymbol{\beta}||_{1}$, lasso encourages sparsity, i.e. reduces the
number of terms in the linear model, as zero-valued model coefficients are
preferred.

In the context of discovering subgrid-scale models, non-linearities can be
introduced by replacing X non-linear functions $G(\mathbf{X})$ of the original
variables. In this study, we construct a model with non-linear variables by
evaluating $d$-order polynomial functions:

$\displaystyle\hat{\boldsymbol{f}}$
$\displaystyle=G^{d}(\mathbf{X})\boldsymbol{\beta},$ (14a) $\displaystyle
G^{d}(\mathbf{X})$
$\displaystyle=\begin{bmatrix}1&X_{11}&X_{12}&\cdots&X_{1n}&X_{11}^{2}&X_{11}X_{12}&\cdots&X_{1n}^{d}\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\
1&X_{m1}&X_{mn}&\cdots&X_{mn}&X_{m1}^{2}&X_{m1}X_{m2}&\cdots&X_{mn}^{d}\end{bmatrix},$
(14b) $\displaystyle\boldsymbol{\beta}$
$\displaystyle=\begin{bmatrix}\beta_{0}&\beta_{1}&\cdots&\beta_{k}\end{bmatrix}^{T},\quad
k=\sum^{d}_{i=1}n^{i}$ (14c)

Equation 14 demonstrates that the dimensionality of this approach scales to
order of polynomial functions $\mathcal{O}(mn^{d})$. Hence, the number of
candidate variables must be reduced for this method to remain tractable. In
this work, we employ the random forest feature importance score to reduce the
number of candidate variables.

##  5 Results

###  5.1 Algebraic SGS models

_A priori_ analysis is performed by comparing SGS stresses $\tau_{ij}^{sgs}$
computed from filtered DNS, with SGS stress modeled by the Vreman (Eq. 9) and
gradient (Eq. 11) models. The performance of the SGS models is evaluated
through the Pearson correlation coefficient, which measures the linear
correlation between two variables. A Pearson correlation of 1.0 and $-1.0$
corresponds to perfectly positive and negative linear relationships,
respectively, whereas a correlation of 0.0 indicates a negligible linear
relationship.

Figure 4 presents the resulting Pearson correlation between exact and
algebraically modeled SGS stresses for three different filter widths
${\overline{\Delta}}$ for all three DNS cases specified in Table 1, at time
$t=1.5t_{I}$. For all three cases and filter sizes, negative correlations and
weak positive correlations ranging from approximately $-0.6$ to 0.4 are
observed for the Vreman model. Negative correlations suggest the presence of
counter-gradient diffusion, which causes the employment of eddy-viscosity
models such as the Vreman model to be ineffective. In all three cases and
three filter sizes, strong positive correlations, ranging from 0.5 to 0.95,
suggest that the gradient model is highly suitable for modeling SGS stresses
in transcritical inert and reacting flows.

(a) Inert case.

(b) Da $=780$.

(c) Da $=10$.

Figure 4: Pearson correlations between exact and algebraically modeled SGS
stresses for three different filter widths ${\overline{\Delta}}$.

The effectiveness of the Vreman and gradient models are further assessed by
examining the conditional Pearson correlation for $\tau^{sgs}_{1i}$ with
respect to the mixture fraction $\widetilde{Z}$ at filter size
${\overline{\Delta}}=2\Delta$. The mixture fraction for the reacting cases
have been evaluated using Bilger’s definition. Figure 5(a) shows that weak
correlations ranging from $-0.4$ to 0.5 are observed throughout the inert
case. In both reacting cases in Figures 5(b) and 5(c), the presence of
countergradient diffusion is much larger than the inert case, as denoted by
the presence of highly negative correlations $(-0.8)$ in the Vreman model. In
the inert case, the gradient model has the highest correlation of
approximately 1.0 in pure methane and pure oxygen, and the lowest correlation
of 0.6 when $Z=0.5$. For the case Da $=780$, the gradient model has the lowest
correlation (0.7) close to the pure oxygen stream, with the correlation
steadily increasing as the mixture approaches stoichiometry
($\widetilde{Z}_{st}=0.2$), after which the correlations remain high (0.85 to
1.0). For the case Da $=10$, the correlations for the gradient model are high
(0.8 to 1.0) throughout the entire mixture.

(a) Inert case.

(b) Da $=780$.

(c) Da $=10$.

Figure 5: Conditional Pearson correlations with respect to mixture fraction
$\widetilde{Z}$ between exact and algebraically modeled SGS stresses
$\tau^{sgs}_{1i}$ for a single filter width ${\overline{\Delta}}=2\Delta$.

The accuracy of the gradient model in predicting the magnitude of SGS stresses
is evaluated by examining the least squares fit between the exact and modeled
SGS stresses. A slope greater than unity indicates underprediction of the
modeled SGS stresses, while a slope less than unity indicates overprediction.
Figure 6 shows that the slopes from the gradient model range from 1 to 4.5.
The average of the slopes is 1.98, which suggests that the gradient model with
a constant coefficient should employ $C_{g}=1/6$, instead of the typical
$C_{g}=1/12$. However, since a wide range of coefficients are observed, a
dynamic gradient model scheme is likely more suited in _a posteriori_
simulations. This is further reinforced by good results from _a posteriori_
evaluations of the dynamic gradient model from transcritical inert DNS [43].

Figure 6: Slopes from a least squares fit of exact and gradient modeled SGS
stress for three different filter widths ${\overline{\Delta}}$.

###  5.2 Random forest SGS models

The _a priori_ analysis performed in Section 5.1 is repeated in this section
for the SGS stresses modeled by random forest regressors. Table 2 summarizes
the random forests employed in this study. All random forests are trained with
snapshots at one eddy turnover time $t=t_{I}$ and tested on the three cases at
$t=1.5t_{I}$. Two different sets of feature, or inputs, are employed to train
the random forests. One feature set corresponds to a domain-blind random
forest RF_BLIND, consisting only of velocity, and the first and second spatial
derivatives of velocity. The other set considers Galilean invariant basis
functions constructed from strain $\widetilde{S}_{ij}$ and rotation
$\widetilde{R}_{ij}$ tensors as features, shown to predict anisotropy well in
a previous study [44]. These Galilean invariant features are used to train the
random forest RF_INFORM. In order to investigate the generalizability of
random forests in the absence of a vast representative dataset, we evaluate
the predictive performance of three additional random forests RF_INERT,
RF_DA780, and RF_DA10, which are trained solely from the inert, Da $=780$, and
Da $=10$ cases, respectively.

Table 2: Random forests employed in this study.

Random forest | RF_BLIND | RF_INFORM | RF_INERT | RF_DA780 | RF_DA10
---|---|---|---|---|---
Features (Input) | $\widetilde{u}_{i}$, $\dfrac{\partial\widetilde{u}_{i}}{\partial x_{j}}$, $\dfrac{\partial^{2}\widetilde{u}_{i}}{\partial x_{j}x_{k}}$ |  $\widetilde{S}_{ij}$, $\widetilde{S}_{ik}\widetilde{S}_{kj}$, $\widetilde{R}_{ik}\widetilde{R}_{kj}$, | $\widetilde{u}_{i}$, $\dfrac{\partial\widetilde{u}_{i}}{\partial x_{j}}$, $\dfrac{\partial^{2}\widetilde{u}_{i}}{\partial x_{j}x_{k}}$ | $\widetilde{u}_{i}$, $\dfrac{\partial\widetilde{u}_{i}}{\partial x_{j}}$, $\dfrac{\partial^{2}\widetilde{u}_{i}}{\partial x_{j}x_{k}}$ | $\widetilde{u}_{i}$, $\dfrac{\partial\widetilde{u}_{i}}{\partial x_{j}}$, $\dfrac{\partial^{2}\widetilde{u}_{i}}{\partial x_{j}x_{k}}$
$\widetilde{S_{ik}}\widetilde{R_{kj}}-\widetilde{R_{ik}}\widetilde{S_{kj}}$
Output | $\tau_{ij}^{sgs}$
Training data ($t=t_{I}$) | All three cases | All three cases | Inert case | Da $=780$ case | Da $=10$ case
Testing data ($t=1.5t_{I}$) | All three cases

Figure 7 presents the Pearson correlation between exact SGS stresses and the
SGS stresses modeled by the random forests RF_BLIND and RF_INFORM. Figure 7(a)
shows that strong correlations (0.4 to 0.95) are observed when the random
forest is trained with an uninformed approach, which is similar to the
gradient model and higher than the Vreman model in Fig. 4. Figure 7(b)
demonstrates that the employment of invariant basis functions as features
decreases the range of correlations (0.35 to 0.9) by 0.05. This small decrease
is likely caused by the additional constraints placed on the random forest
when forming a hypothesis space.

(a) RF_BLIND.

(b) RF_INFORM.

Figure 7: Pearson correlation between exact and random forest modeled SGS
stresses for three different filter widths ${\overline{\Delta}}$.

Figure 8 presents the Pearson correlations between exact and random forest SGS
stresses $\tau^{sgs}_{1i}$ conditioned to mixture fraction $\widetilde{Z}$ at
${\overline{\Delta}}=2\Delta$. In the inert case, shown in Fig. 8(a), highest
correlation from RF_BLIND of approximately 0.95 is observed in pure methane
and pure oxygen, and lowest correlation of 0.5 when $Z=0.5$. For the case Da
$=780$ in Fig. 8(b), RF_BLIND possesses the lowest correlation (0.7) close to
the oxygen stream, with the correlation steadily increasing as the mixture
approaches stoichiometric conditions ($\widetilde{Z}_{st}=0.2$), after which
the correlations remain high (0.85 to 1.0). For the case Da $=10$, shown in
Fig. 8(c), the correlations for the gradient model are high (0.8 to 1.0)
throughout the entire mixture. The conditional Pearson correlation produced
from RF_BLIND in all three cases are similar qualitatively and quantitatively
to correlations from the gradient model in Fig. 5. This suggests that RF_BLIND
has approximated a function similar to the gradient model, even when trained
solely on exact SGS stresses and without any prior knowledge of the gradient
Model. The correlations from RF_INFORM share similar qualitative behaviors as
the correlations from RF_BLIND, but with up to a 0.2 lower values.

(a) Inert case.

(b) Da $=780$.

(c) Da $=10$.

Figure 8: Conditional Pearson correlations as a function of mixture fraction
$\widetilde{Z}$ between exact and random forest modeled SGS stresses
$\tau^{sgs}_{1i}$ for a single filter width ${\overline{\Delta}}=2\Delta$.

Figure 9 presents slopes from a least squares fit between the exact and the
random forest SGS stresses. Figure 9(a) shows that the slopes from RF_BLIND
range from 0.25 to 1.6, with an average slope of 0.96, which demonstrates
excellent agreement between modeled and exact magnitudes of SGS stresses. The
employment of invariant features leads to lower slopes (0.25 to 1.35), with an
average slope of 0.867, as presented in Fig. 9(b). The use of the invariant
feature set not only leads to lower correlations but also to an overprediction
in magnitudes of SGS stresses.

(a) RF_BLIND.

(b) RF_INFORM.

Figure 9: Slopes from a least squares fit of exact and random forest modeled
SGS stress for three different filter widths ${\overline{\Delta}}$.

Figure 10 compares the exact and modeled SGS stress
$\tau^{sgs}_{12}/{\overline{\rho}}$ at filter width
${\overline{\Delta}}=4\Delta$. In the inert case, both SGS stresses from the
gradient model and RF_BLIND are in good agreement with the exact term. For Da
= 780, the gradient model is in better agreement with the exact term than
RF_BLIND. This is further supported by the difference in Pearson correlation
for this particular case shown by the gradient model (0.9) and RF_BLIND (0.6)
in Figs. 4 and 7, respectively. For Da = 10, RF_BLIND predicts the magnitude
of the SGS stress better than the gradient model, especially near the domain
boundaries, which is also observed in the slopes shown by RF_BLIND (0.9) and
the gradient model (1.5) shown in Figs. 6 and 9.

Figure 10: Comparison of exact and modeled SGS stress
$\tau^{sgs}_{12}/{\overline{\rho}}$ [m2s-2] at filter width
${\overline{\Delta}}=4\Delta$ at axial location $x=0$.

Figure 11 presents Pearson correlations from examining the generalizability of
random forests in the presence of limited data. Random forest RF_INERT
demonstrates a similar range of correlations (0.5 to 0.85) to RF_ALL when
tested on the inert case with a filter size ${\overline{\Delta}}=2{\Delta}$.
However, lower ranges are observed for RF_INERT when tested on the cases Da
$=780$ (0.4 to 0.75) and Da $=10$ (0.5 to 0.9). RF_DA780 also possesses a
similar correlation as RF_BLIND when tested on case Da $=780$ (0.5 to 0.9),
but worse correlations when tested on the inert case (0.4 to 0.8) and case Da
$=10$ (0.8 to 0.9). Lastly, RF_DA10 performs similarly to RF_BLIND when tested
on Da $=10$ (0.85 to 0.95) but performs worse when tested on the inert (0.5 to
0.8) and Da $=780$ (0.55 to 0.8) cases. These three random forests perform as
well as RF_BLIND on a test set that is represented well by the training set.
However, the effectiveness of random forests decreases when modeling on out-
of-sample distributions. Nevertheless, these out-of-sample predictions are
more accurate than the Vreman model, thus demonstrating a appreciable degree
of generalizability.

Figure 11:  Pearson correlations between exact and random forest modeled SGS
stresses for a single filter width ${\overline{\Delta}}=2\Delta$.

###  5.3 Data-driven discovery of SGS models

In this section, we examine how the interpretability of random forests can be
employed as a tool for model discovery.

Figure 12 presents feature importance scores extracted from RF_BLIND for
$\tau^{sgs}_{1i}$. For all three SGS stresses $\tau^{sgs}_{1i}$ shown, the
highest scores are from $\dfrac{\partial\widetilde{u}_{1}}{\partial x_{k}}$
and $\dfrac{\partial\widetilde{u}_{i}}{\partial x_{k}}$ for three spatial
dimensions. We employ this observation to formulate a sparse regression
problem (see Eqs. 12, 13 and 14):

$\displaystyle\frac{\tau^{sgs}_{ij}}{{\overline{\rho}}u^{\prime
2}}=f_{ij}\left[G^{d=2}\left(\frac{{\overline{\Delta}}}{u^{\prime}}\dfrac{\partial\widetilde{u}_{i}}{\partial
x_{k}},\frac{{\overline{\Delta}}}{u^{\prime}}\dfrac{\partial\widetilde{u}_{j}}{\partial
x_{k}}\right)\right]$ (15)

where the independent variables consist of 2nd-order polynomial functions of
the non-dimensionalized selected features. Equation 15 is non-dimensionalized
by density, filter width and initial root-mean-squared velocity to ensure
dimensional consistency in the final model. This is essential for improving
the dimensionality of this sparse regression problem. Since the dimensionality
scales with $n^{d}$ for $n$ number of candidate variables, as discussed in
Section 4.3, the employment of the feature importance score for reducing 30
candidate variables to 6 candidate variables results in a 25-fold improvement
in dimensionality.

Figure 12:  Fifteen feature importance scores from RF_BLIND. The other fifteen
features, with importance scores less than 0.02, are not shown for brevity.

The following equations present the SGS model that resulted from applying
sparse regression:

$\displaystyle\tau_{11}^{sgs}$
$\displaystyle\simeq{\overline{\rho}}{\overline{\Delta}}^{2}\left(0.116\dfrac{\partial\widetilde{u}_{1}}{\partial
x_{1}}\dfrac{\partial\widetilde{u}_{1}}{\partial
x_{1}}+0.191\dfrac{\partial\widetilde{u}_{1}}{\partial
x_{2}}\dfrac{\partial\widetilde{u}_{1}}{\partial
x_{2}}+0.207\dfrac{\partial\widetilde{u}_{1}}{\partial
x_{3}}\dfrac{\partial\widetilde{u}_{1}}{\partial x_{3}}\right)$ (16a)
$\displaystyle\tau_{12}^{sgs}$
$\displaystyle\simeq{\overline{\rho}}{\overline{\Delta}}^{2}\left(0.113\dfrac{\partial\widetilde{u}_{1}}{\partial
x_{1}}\dfrac{\partial\widetilde{u}_{2}}{\partial
x_{1}}+0.102\dfrac{\partial\widetilde{u}_{1}}{\partial
x_{2}}\dfrac{\partial\widetilde{u}_{2}}{\partial
x_{2}}+0.134\dfrac{\partial\widetilde{u}_{1}}{\partial
x_{3}}\dfrac{\partial\widetilde{u}_{2}}{\partial x_{3}}\right)$ (16b)
$\displaystyle\tau_{13}^{sgs}$
$\displaystyle\simeq{\overline{\rho}}{\overline{\Delta}}^{2}\left(0.119\dfrac{\partial\widetilde{u}_{1}}{\partial
x_{1}}\dfrac{\partial\widetilde{u}_{3}}{\partial
x_{1}}+0.117\dfrac{\partial\widetilde{u}_{1}}{\partial
x_{2}}\dfrac{\partial\widetilde{u}_{3}}{\partial
x_{2}}+0.109\dfrac{\partial\widetilde{u}_{1}}{\partial
x_{3}}\dfrac{\partial\widetilde{u}_{3}}{\partial x_{3}}\right)$ (16c)
$\displaystyle\tau_{22}^{sgs}$
$\displaystyle\simeq{\overline{\rho}}{\overline{\Delta}}^{2}\left(0.215\dfrac{\partial\widetilde{u}_{2}}{\partial
x_{1}}\dfrac{\partial\widetilde{u}_{2}}{\partial
x_{1}}+0.135\dfrac{\partial\widetilde{u}_{2}}{\partial
x_{2}}\dfrac{\partial\widetilde{u}_{2}}{\partial
x_{2}}+0.164\dfrac{\partial\widetilde{u}_{2}}{\partial
x_{3}}\dfrac{\partial\widetilde{u}_{2}}{\partial x_{3}}\right)$ (16d)
$\displaystyle\tau_{23}^{sgs}$
$\displaystyle\simeq{\overline{\rho}}{\overline{\Delta}}^{2}\left(0.123\dfrac{\partial\widetilde{u}_{2}}{\partial
x_{1}}\dfrac{\partial\widetilde{u}_{3}}{\partial
x_{1}}+0.116\dfrac{\partial\widetilde{u}_{2}}{\partial
x_{2}}\dfrac{\partial\widetilde{u}_{3}}{\partial
x_{2}}+0.134\dfrac{\partial\widetilde{u}_{2}}{\partial
x_{3}}\dfrac{\partial\widetilde{u}_{3}}{\partial x_{3}}\right)$ (16e)
$\displaystyle\tau_{33}^{sgs}$
$\displaystyle\simeq{\overline{\rho}}{\overline{\Delta}}^{2}\left(0.251\dfrac{\partial\widetilde{u}_{3}}{\partial
x_{1}}\dfrac{\partial\widetilde{u}_{3}}{\partial
x_{1}}+0.177\dfrac{\partial\widetilde{u}_{3}}{\partial
x_{2}}\dfrac{\partial\widetilde{u}_{3}}{\partial
x_{2}}+0.124\dfrac{\partial\widetilde{u}_{3}}{\partial
x_{3}}\dfrac{\partial\widetilde{u}_{3}}{\partial x_{3}}\right)$ (16f)

The resulting model can be rewritten as:

$\displaystyle\tau_{ij}^{sgs}$
$\displaystyle\simeq{\overline{\rho}}{\overline{\Delta}}^{2}\left(C_{1}\dfrac{\partial\widetilde{u}_{i}}{\partial
x_{1}}\dfrac{\partial\widetilde{u}_{j}}{\partial
x_{1}}+C_{2}\dfrac{\partial\widetilde{u}_{i}}{\partial
x_{2}}\dfrac{\partial\widetilde{u}_{j}}{\partial
x_{2}}+C_{3}\dfrac{\partial\widetilde{u}_{i}}{\partial
x_{3}}\dfrac{\partial\widetilde{u}_{j}}{\partial x_{3}}\right)$ (17)

where the resulting model coefficients $C_{\\{1,2,3\\}}$ range from 0.102 to
0.251. Equation 17 is similar in form to the gradient model (Eq. 11), but
possesses three model coefficients instead of one. By observing that
$C_{\\{1,2,3\\}}$ are of the same order of magnitudes, and collapsing the
three coefficients by evaluating the average model coefficients, we recover
the gradient model:

$\displaystyle\tau_{ij}^{sgs}$
$\displaystyle\simeq{\overline{\rho}}C_{x}{\overline{\Delta}}^{2}\dfrac{\partial\widetilde{u}_{i}}{\partial
x_{k}}\dfrac{\partial\widetilde{u}_{j}}{\partial x_{k}}$ (18)

where the model coefficient $C_{x}=0.147$ is similar in value to the suggested
model coefficient of 0.167 from Section 5.1. This result demonstrates that the
employment of sparse regression, in conjunction with random forest feature
importance can be employed to discover an algebraic expression, similar to the
effective gradient model, for modeling subgrid scale stresses in transcritical
flows.

Since the present method relies on the random forest feature importance score,
a statistical test must be employed to test for the effects of significant
correlation amongst the features. If multiple features in the modeling basis
are significantly correlated, they act as exchangeable surrogates for each
other during the calculation of feature importance scores (this is similar to
the phenomenon of multicollinearity in classical statistics [45]). Under such
conditions, metrics such as the MDI are susceptible to correlation bias, and
can generate erroneous importance scores [46, 47]. As a note, almost all
algorithms for estimating feature importance, including SHAP (Shapley additive
explanations) [48] exhibit such correlation bias. As an alternative, Principal
Component Analysis may be utilized to engender orthogonal bases for new
features that are independent. However, these derived features are often
difficult to ascribe physical meanings to, obfuscating their utility toward
interpretability.

Multicollinearity is likely not an issue for RF_BLIND since we know from
finite differencing schemes that the spatial derivative of a quantity in a
cell cannot be evaluated locally, but must be evaluated from neighboring
cells. Nevertheless, we utilize the Spearman correlation as a statistical test
for evaluating the correlation amongst the features in the modeling basis.
While the Pearson correlation is a statistical tool used for evaluating linear
relationships, the Spearman correlation evaluates the monotonicity of
variables in both linear and non-linear functions, i.e., whether the
increasing or decreasing trend is being preserved. Spearman correlations of 1
and $-1$ correspond to a perfect monotonic relationship, while 0 corresponds
to a negligible monotonic relationship. Figure 13 shows that Spearman
correlations between different features from RF_BLIND are weak (between $-0.4$
and 0.4), which indicates weak correlation amongst the features.

Figure 13: Spearman correlation matrix for selected features from RF_BLIND.
Features with correlations less than 0.2 are not shown for brevity.

##  6 Conclusions

DNS of inert and reacting transcritical LOX/GCH4 non-premixed mixtures under
decaying turbulence were performed. Pressure and temperature were chosen to
correspond to conditions in rocket combustors. A priori analysis was conducted
by comparing exact subgrid-scale stresses from Favre-filtered DNS data with
algebraic and data-driven SGS models.

A priori analysis showed that the SGS stresses evaluated by Vreman SGS model
correlated poorly with the corresponding exact terms, likely caused by the
presence of countergradient diffusion. In contrast, good correlations are seen
from the gradient SGS model. Results demonstrated a wide range of magnitude
errors in the gradient model, which suggests that a dynamic gradient model
approach is suited in _a posteriori_ simulations. Random forests demonstrated
high correlations when trained on datasets which are representative of the
test sets, with reasonable predictions for the magnitude of subgrid-scale
stresses. However, correlations were shown to decrease significantly when
tested out-of-sample.

Sparse regression was performed to discover an algebraic expression for SGS
stresses from non-linear transformations of velocity and its derivatives. The
interpretability of random forests was demonstrated to reduce the
dimensionality of the sparse regression problem by 25 times, by employing the
feature importance score for variable selection. The derived algebraic
expression was shown to be similar to the gradient model.

Results demonstrate that random forests can perform as effectively as suitable
algebraic models when modeling subgrid stresses, if trained on a sufficiently
representative database. However, in the absence of such a database, this good
performance is not replicated. Nevertheless, the employment of random forests
can provide insight into the discovery of subgrid-scale models through the
feature importance score. Additional work on interpretable deep learning,
which take advantage of autoencoders, could be an exciting direction for
future work. In addition, the present study should be complemented with an a
posteriori study, and extended to other SGS closure terms that form chemical
source terms and non-linear state equations, to generate further insight.

## Acknowledgments

The authors gratefully acknowledge financial support from NASA with award No.
80NSS C18C0207, the Department of Energy, National Nuclear Security
Administration under award No. DE-NA0003968, and the German Research
Foundation (Deutsche Forschungsgemeinschaft – DFG) in the framework of the
Sonderforschungsbereich Transregio 40. Resources supporting this work are
provided by the National Energy Research Scientific Computing Center, a U.S.
Department of Energy Office of Science User Facility operated under contract
No. DE-AC02-05CH11231.

## References

  * Bellan [2020] J. Bellan, Future challenges in the modelling and simulations of high-pressure flows, Combust. Sci. Technol. 192 (2020) 1199–1218.
  * Chang et al. [1996] C. Chang, Y. Zhang, K. N. C. Bray, B. Rogg, Modelling and simulation of autoignition under simulated diesel-engine conditions, Combust. Sci. Technol. 113 (1996) 205–219.
  * Zhu et al. [1996] M. Zhu, K. N. C. Bray, B. Rogg, PDF modelling of spray autoignition in high pressure turbulent flows, Combust. Sci. Technol. 120 (1996) 357–379.
  * Ma et al. [2019] P. C. Ma, H. Wu, T. Jaravel, L. Bravo, M. Ihme, Large-eddy simulations of transcritical injection and auto-ignition using diffuse-interface method and finite-rate chemistry, Proc. Combust. Inst. 37 (2019) 3303–3310.
  * Oefelein [2019] J. C. Oefelein, Advances in modeling supercritical fluid behavior and combustion in high-pressure propulsion systems, AIAA Paper 2019-0634 (2019).
  * Chung et al. [2020] W. T. Chung, P. C. Ma, M. Ihme, Examination of diesel spray combustion in supercritical ambient fluid using large-eddy simulations, Int. J. Engine Res. 21 (2020) 122–133.
  * Vreman [2004] A. W. Vreman, An eddy-viscosity subgrid-scale model for turbulent shear flow: Algebraic theory and applications, Phys. Fluids 16 (2004) 3670–3681.
  * Clark et al. [1979] R. A. Clark, J. H. Ferziger, W. C. Reynolds, Evaluation of subgrid-scale models using an accurately simulated turbulent flow, J. Fluid Mech. 91 (1979) 1–16.
  * Selle et al. [2007] L. C. Selle, N. A. Okong’o, J. Bellan, K. G. Harstad, Modelling of subgrid-scale phenomena in supercritical transitional mixing layers: an a priori study, J. Fluid Mech. 593 (2007) 57–91.
  * Smagorinsky [1963] J. Smagorinsky, General circulation experiments with the primitive equations, Mon. Weather Rev. 91 (1963) 99–164.
  * Bardina et al. [1980] J. Bardina, J. Ferziger, W. C. Reynolds, Improved subgrid-scale models for large-eddy simulation, AIAA Paper 1980-1357 (1980).
  * Unnikrishnan et al. [2019] U. Unnikrishnan, J. C. Oefelein, V. Yang, A priori analysis of subfilter scalar covariance fields in turbulent reacting LOX-CH4 mixing layers, AIAA Paper 2019-1495 (2019).
  * Moin et al. [1991] P. Moin, K. Squires, W. Cabot, S. Lee, A dynamic subgrid-scale model for compressible turbulence and scalar transport, Phys. Fluids A 3 (1991) 2746–2757.
  * Lapeyre et al. [2019] C. J. Lapeyre, A. Misdariis, N. Cazard, D. Veynante, T. Poinsot, Training convolutional neural networks to estimate turbulent sub-grid scale reaction rates, Combust. Flame 203 (2019) 255–264.
  * Seltz et al. [2019] A. Seltz, P. Domingo, L. Vervisch, Z. M. Nikolaou, Direct mapping from LES resolved scales to filtered-flame generated manifolds using convolutional neural networks, Combust. Flame 210 (2019) 71–82.
  * Henry de Frahan et al. [2019] M. T. Henry de Frahan, S. Yellapantula, R. King, M. S. Day, R. W. Grout, Deep learning for presumed probability density function models, Combust. Flame 208 (2019) 436–450.
  * Ranade and Echekki [2019] R. Ranade, T. Echekki, A framework for data-based turbulent combustion closure: A posteriori validation, Combust. Flame 210 (2019) 279–291.
  * Lakshmisha et al. [1995] K. N. Lakshmisha, B. Rogg, K. N. C. Bray, PDF modelling of autoignition in nonpremixed turbulent flows, Combust. Sci. Technol. 105 (1995) 229–243.
  * Peng and Robinson [1976] D.-Y. Peng, D. B. Robinson, A new two-constant equation of state, Ind. Eng. Chem. Fundam. 15 (1976) 59–64.
  * Poling et al. [2001] B. E. Poling, J. M. Prausnitz, J. P. O’Connell, The Properties of Gases and Liquids, McGraw-Hill, 2001.
  * Ma et al. [2017] P. C. Ma, Y. Lv, M. Ihme, An entropy-stable hybrid scheme for simulations of transcritical real-fluid flows, J. Comput. Phys. 340 (2017) 330–357.
  * Franzelli et al. [2012] B. Franzelli, E. Riber, L. Y. Gicquel, T. Poinsot, Large eddy simulation of combustion instabilities in a lean partially premixed swirled flame, Combust. Flame 159 (2012) 621 – 637.
  * Chong et al. [2017] S. T. Chong, Y. Tang, M. Hassanaly, V. Raman, Turbulent mixing and combustion of supercritical jets, AIAA Paper 2017-0141 (2017).
  * Bellan [2017] J. Bellan, Direct numerical simulation of a high-pressure turbulent reacting temporal mixing layer, Combust. Flame 176 (2017) 245–262.
  * Bushe et al. [2020] W. K. Bushe, C. Devaud, J. Bellan, A priori evaluation of the double-conditioned conditional source-term estimation model for high-pressure heptane turbulent combustion using DNS data obtained with one-step chemistry, Combust. Flame 217 (2020) 131–151.
  * Takahashi [1975] S. Takahashi, Preparation of a generalized chart for the diffusion coefficients of gases at high pressures, J. Chem. Eng. Japan 7 (1975) 417–420.
  * Chung et al. [1988] T. H. Chung, M. Ajlan, L. L. Lee, K. E. Starling, Generalized multiparameter correlation for nonpolar and polar fluid transport properties, Ind. Eng. Chem. Res 27 (1988) 671–679.
  * Khalighi et al. [2011] Y. Khalighi, J. W. Nichols, F. Ham, S. K. Lele, P. Moin, Unstructured large eddy simulation for prediction of noise issued from turbulent jets in various configurations, AIAA Paper 2011-2886 (2011).
  * Chung et al. [2021] W. T. Chung, A. A. Mishra, N. Perakis, M. Ihme, Data-assisted combustion simulations with dynamic submodel assignment using random forests, Combust. Flame 227 (2021) 172–185.
  * Wu et al. [2019] H. Wu, P. C. Ma, M. Ihme, Efficient time-stepping techniques for simulating turbulent reactive flows with stiff chemistry, Comput. Phys. Commun. 243 (2019) 81–96.
  * Huo and Yang [2011] H. Huo, V. Yang, Supercritical LOX/methane combustion of a shear coaxial injector, AIAA Paper 2011-326 (2011).
  * Unnikrishnan et al. [2018] U. Unnikrishnan, J. C. Oefelein, V. Yang, Direct numerical simulation of a turbulent reacting liquid-oxygen/methane mixing layer at supercritical pressure, AIAA Paper 2018-4564 (2018).
  * Saad et al. [2017] T. Saad, D. Cline, R. Stoll, J. C. Sutherland, Scalable tools for generating synthetic isotropic turbulence with arbitrary spectra, AIAA J. 55 (2017) 327–331.
  * Bailly and Juves [1999] C. Bailly, D. Juves, A stochastic approach to compute subsonic-noise using linearized Euler’s equations, AIAA Paper 1999-1872 (1999).
  * Williams [2006] F. Williams, Descriptions of nonpremixed turbulent combustion, AIAA Paper 2006-1505 (2006).
  * Schoepplein et al. [2018] M. Schoepplein, J. Weatheritt, R. Sandberg, M. Talei, M. Klein, Application of an evolutionary algorithm to les modelling of turbulent transport in premixed flames, J. Comput. Phys. 374 (2018) 1166 – 1179.
  * Breiman [2001] L. Breiman, Random forests, Mach. Learn. 45 (2001) 5–32.
  * Breiman et al. [1984] L. Breiman, J. Friedman, R. Olshen, C. Stone, Classification and Regression Trees, Routledge, 1984\.
  * Amit et al. [1997] Y. Amit, D. Geman, K. Wilder, Joint induction of shape features and tree classifiers, IEEE Trans. Pattern Anal. and Mach. Intell. 19 (1997) 1300–1305.
  * Louppe et al. [2013] G. Louppe, L. Wehenkel, A. Sutera, P. Geurts, Understanding variable importances in forests of randomized trees, Adv. Neural Inf. Process. Syst. 26 (2013).
  * Pedregosa et al. [2011] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, E. Duchesnay, Scikit-learn: Machine learning in Python, J. Mach. Learn. Res. 12 (2011) 2825–2830.
  * Tibshirani [1996] R. Tibshirani, Regression shrinkage and selection via the lasso, J. R. Stat. Soc. Series B Stat. Methodol. 58 (1996) 267–288.
  * Taskinoglu and Bellan [2010] E. S. Taskinoglu, J. Bellan, A posteriori study using a DNS database describing fluid disintegration and binary-species mixing under supercritical pressure: heptane and nitrogen, J. Fluid Mech. 645 (2010) 211–254.
  * Ling et al. [2016] J. Ling, R. Jones, J. Templeton, Machine learning strategies for systems with invariance properties, J. Comput. Phys. 318 (2016) 22–35.
  * Kumar et al. [1975] T. K. Kumar, et al., Multicollinearity in regression analysis, Rev. Econ. Stat 57 (1975) 365–366.
  * Altmann et al. [2010] A. Altmann, L. Toloşi, O. Sander, T. Lengauer, Permutation importance: a corrected feature importance measure, Bioinformatics 26 (2010) 1340–1347.
  * Toloşi and Lengauer [2011] L. Toloşi, T. Lengauer, Classification with correlated features: unreliability of feature ranking and solutions, Bioinformatics 27 (2011) 1986–1994.
  * Lundberg and Lee [2017] S. M. Lundberg, S.-I. Lee, A unified approach to interpreting model predictions, Adv. Neural Inf. Process. Syst. 30 (2017) 4765–4774.