# Support-vector-machine with Bayesian optimization for lithofacies
classification using elastic properties

Yohei Nishitsuji
Department of Geoscience and Engineering
Delft University of Technology
Delft, The Netherlands
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
&Jalil Nasseri
Summit Exploration and Production Limited
London, UK
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>

###### Abstract

We investigate an applicability of Bayesian-optimization (BO) to optimize
hyperparameters associated with support-vector-machine (SVM) in order to
classify facies using elastic properties derived from well data in the East
Central Graben, UKCS. The cross-plot products of the field dataset appear to
be successfully classified with non-linear boundaries. Although there are a
few factors to be predetermined in the BO scheme such as an iteration number
to deal with a trade-off between the prediction accuracy and the computational
cost, this approach effectively reduces possible human subjectivity connected
to the architecture of the SVM. Our proposed workflow might be beneficial in
resource-exploration and development in terms of subsurface objective
technical evaluations.

_K_ eywords Deep learning $\cdot$ Bayesian optimization $\cdot$ Hyperparameter
$\cdot$ SVM $\cdot$ Oil&Gas $\cdot$ Elastic properties cross-plots $\cdot$
Quantitative interpretation

##  1 Introduction

Quantitative interpretation (QI) is of importance for the purpose of improving
insight from subsurface for exploration and development of natural resources.
Among many seismic attributes, seismic inversion products are the most
commonly used products for QI [e.g., 2, 6, 17]. The elastic properties cross-
plots are often used to assess different facies (e.g., sand and shale) and
presence of hydrocarbon in reservoirs. When the properties are linearly
separable in the original cross-plot space (acoustic impedance, $AI$ and
velocity ratio, $\nicefrac{{V_{P}}}{{V_{S}}}$), a linear projection can be
useful because of its simple and fast calculations to distinguish facies.
However, the cross-plot data in practice are often exhibited as complex
clusters, which essentially requires non-linear separations. To tackle this
problem, we come up with support-vector-machines [SVM; e.g., 15, 14], a
machine-learning algorithm.

The SVM is originally developed to solve a binary (two) classification problem
by a linear separation. Having called kernel-function, however, the SVM can
handle non-linear problems too. Moreover, not only the binary problems but
also multiclass-classification problems can be solved when an appropriate
aggregation strategy is implemented. A couple of hyperparameters require to be
optimized to perform the SVM effectively and efficiently for any given
problems. In terms of finding optimal hyperparameters, the simplest approach
can be an exhaustive-grid search [e.g., 16]. This method seeks every single
grid within a given searching range, literally one by one. Even though when
the number of the SVM’s hyperparameters to be optimized is small, this method
becomes an inefficient approach when the searching range becomes huge. A time-
efficient approach we would like to apply in this study is called Bayesian-
optimization [BO; e.g., 7].

The BO is a global optimization algorithm that finds hyperparameters of black-
box (unknown) functions in as few iterations as possible by searching a good
balance between so-called exploration and exploitation [e.g., 4]. There are
many industries and applications that come up with the BO in practice such as
social network, accommodation fare metasearch and travel route aggregator. The
exploration factor seeks where the variance of the current posterior
measurement is high for next round of evaluation (iteration), whereas the
exploitation factor focuses on where the means of the posterior is low. In
other words, the exploration part completely ignores what the posterior has
already obtained so far, whilst the exploitation part completely ignores what
the posterior has not obtained so far.

In this study, we utilize the BO to optimize the SVM’s hyperparameters for the
multiclass classification problems associated with elastic impedance cross-
plot data obtained from well data in the East Central Graben in UKCS (Figure
1). We use 10 wells in this study including 22/9-5, 22/14b-5, 22/15-3,
22/15-4, 22/24c-11, 22/30a-16, 23/16d-6, 23/16b-9, 2316f-12, and 23/21-5 (UK
National Data Repository: https://ndr.ogauthority.co.uk). The target interval
in this study is the Paleocene sediments where the Forties member of Sele
Formation is the main reservoir sandstone.

Figure 1: The chosen 10 wells (see the introduction for more details) for the
study in the East Central Graben in the UK North Sea.

##  2 The Bayesian-Optimization Applied to The Support-Vector-Machine

###  2.1 Multiclass Classification of The Support-Vector-Machine

In this section, we briefly introduce the multiclass classifications by the
SVM. The SVM, which is based on statistical learning theory and structural
risk minimization, was initially developed for the purpose of separating
binary problems [e.g., 14]. We show a schematic of the binary problem using
the SVM in Figure 2. Considering of actual classifications in practice,
however, there are many cases that we need to deal with more than three
classes. For example, five chosen classes in this study are the geological
facies such as shale, tuff, brine-, oil-, gas-sandstone. A plenty number of
studies connected to aggregation strategy has been reported to deal with
multiclass classifications [e.g., 5]. To date, this is one of ongoing-research
topics in machine-learning community.

Figure 2: Schematic image how SVM classifies the binary problem non-linearly.
With the kernel trick, the original input dimension is implicitly converted to
the feature dimension where the decision boundary can linearly separate
classes. The estimated linear boundary in the feature dimension then converts
back to the original dimension being as the non-linear boundary.

We assume that we have a training dataset $T$ which has $p$-th classes for the
multiclass classifications by the SVM, such that:

$\displaystyle
T=\left\\{\left(\mathbf{x}_{1},y_{1}\right),\left(\mathbf{x}_{2},y_{2}\right),\cdots\left(\mathbf{x}_{n},y_{n}\right)\right\\},\,\mathbf{x}_{i}\in\mathbb{R^{\mathrm{\mathit{m}}}},\,y_{i}\in\\{1,2,\cdots,p\\},$
(1)

where $\mathbf{x}_{i}$ and $y_{i}$ denotes $i$-th feature vector whose feature
dimension is $m$ and their associated classes (labels). For example, $p=2$
means a binary-class problem. The SVM minimizes a difference between the
expected risk, which we cannot directly obtain due to the lack of
probabilistic-distribution information including unknown data which are not
observed yet, and empirical risk [e.g., 14]. The minimization allows us to
find a maximum margin which divides $\mathbf{x}_{i}$ into two classes if
$p=2$. The decision boundaries drawn by the margins can be either linear or
non-linear via hyperplane using so-called kernel function. Among several
kernel functions, the gaussian kernel is commonly used because of its better
performance than the others [e.g., 10]. Therefore, we adopt the gaussian
kernel in this study.

There are two hyperparameters require optimization for a possible gain in
performance of the SVM. The first parameter is called sigma, which presents in
a gaussian kernel. This parameter determines the complexity (non-linearity) of
the margins. The second parameter is penalization, which grants the rate of
misclassification. The penalization parameter can be derived from the
Lagrangian-dual-problem for the SVM (not shown here for the sake of brevity).
Moreover, BO determines the type of aggregation strategy of either one-versus-
one (OvO) or one-versus-all (OvA) for solving multi binary class problems.

###  2.2 The Gaussian-Process

We adapt the BO for the optimization framework and determine the
aforementioned two hyperparameters and the aggregation strategy for optimizing
SVM. With the regime of the BO, a probability distribution to an unknown-
target function $f$, in which hyperparameters $h_{i}$ are present, is assumed
to follow the gaussian-process prior:

$\displaystyle p\left(\begin{array}[]{c}f(h_{1})\\\ f(h_{2})\\\ \vdots\\\
f(h_{n})\end{array}\right)\sim\mathcal{N}\left\\{\left(\begin{array}[]{c}\mu(h_{1})\\\
\mu(h_{2})\\\ \vdots\\\
\mu(h_{n})\end{array}\right),\left(\begin{array}[]{cccc}K(h_{1},h_{1})&K(h_{1},h_{2})&\cdots&K(h_{1},h_{n})\\\
K(h_{2},h_{1})&K(h_{2},h_{2})&\cdots&K(h_{2},h_{n})\\\
\vdots&\vdots&\ddots&\vdots\\\
K(h_{n},h_{1})&K(h_{n},h_{2})&\cdots&K(h_{n},h_{n})\end{array}\right)\right\\},$
(14)

where $\mathcal{N}$ , $\mu$ and $K$ indicate normal (gaussian) distribution,
the mean function (usually set to be zero) and the covariance function which
can be derived from the utilized the gaussian kernel. As one can see from
equation 14, the gaussian process, which is completely defined by the mean
function as well as the covariance function, is a group of random variables
whose distributions are all based on the gaussian distributions. While the
covariance is independent of observations, the mean is a linear combination of
each observation. For the brevity, we simplify equation 14 using matrix-vector
notation:

$\displaystyle p(\mathbf{f})\sim\mathcal{N}(\mathbf{f};$
$\displaystyle\boldsymbol{\mu},\mathbf{K}),$ (15)

where $\mathbf{K}$ is a $n\times n$ matrix, whilst $\mathbf{f}$ and
$\boldsymbol{\mu}$ are $n\times 1$ column vectors. $p(\mathbf{f})$ is the
prior information that can be used in the Bayesian regression [e.g., 11],
which is given by:

$\displaystyle
p(\mathbf{f}\mid\mathbf{D})=\frac{p(\mathbf{f})p(\mathbf{D}\mid\mathbf{f})}{p(\mathbf{D})}$
, (16)

where $\mathbf{D}$ is a set of observation data (also often called evidence)
in the matrix-vector notation of $D_{1:n}\left\\{h_{1:n},f(h_{1:n})\right\\}$.
Given new evidence, a posterior probability, the left-hand side term in
equation 16, is updated. Note that we assume noise-free case in equation 16
for simplicity. By solving equation 16, we expect to choose where to observe
the function next. The following equation is adapted to make the choice to be
fully automated:

$\displaystyle\mathbf{h}_{n+1}=\underset{\mathbf{h}}{argmax}\,\alpha(\mathbf{h};\mathbf{D})$
, (17)

where $\alpha$ is an acquisition function, which contributes to the objective
function $\underset{\mathbf{h}}{argmax\,}\mathbf{f}.$ Although several
policies (tasks and criteria) associated with the acquisition function have
been ever proposed in the machine-learning community [e.g., 9, 7, 12], the
following policy explained hereinafter is used in this work.

###  2.3 The Expected-Improvement-Per-Second

The expected improvement, developed by Jones et al. [8] based on the pioneer
work by Mockus et al. [7], is commonly used within the number of the
acquisition functions. The policy of this function is to evaluate the expected
value of possible improvement of the objective function. A searching point
where to observe next is based on the largest expected improvement. The
acquisition function of the expected improvement can be written as:

$\displaystyle\alpha=\int_{-\infty}^{\infty}max(0,\mathbf{f}-\mathbf{f}_{best},)\mathcal{N}(\mathbf{f};\boldsymbol{\mu},\mathbf{K})d\mathbf{f}$
, (18)

in which the next searching point with the largest expected improvement is
automatically selected one after the other. It can be noted that the left-hand
side term in the equation 18 increases when the mean function and the
covariance (Kernel) function decreases and increases, respectively. Therefore,
this policy implicitly has a trade-off between the exploitation factor (driven
by low mean) and the exploration factor (driven by high covariance).
Tendentially, an explicit trade-off parameter exists in other policies such as
the upper-confidence based policy [e.g., 13] which we do not use here.

Snoek et al. [12] improved the performance of the expected improvement by
putting the time-weight on its evaluation time, which is called the expected-
improvement-per-second (EIPS). The EIPS seeks a point to evaluate the next
from where the value of the expected improvement per second is largest. The
EIPS finds the optimal parameters faster than using the expected improvement
by the Markov-chain-Monte-Carlo [e.g., 12]. In our study, we adopt the EIPS
policy for the BO to automatically optimize the hyperparameters of the SVM
with respect to solve the facies-classification problems associated with the
elastic properties cross-plot data.

##  3 Implementation Test

We carry out a synthetic test to check feasibility of our computational
implementation on an $AI$-$\nicefrac{{V_{P}}}{{V_{S}}}$ cross-plot prior to
the actual classification using the field dataset. By means of the feasibility
test we assure that there are no technical problems caused by our
implementation per se. The number of target facies defined by the well data
for classification is five in this study. We therefore synthesize five
classes. We generate 1,000 supervisors (training data) and 50 test data for
each class using a gaussian randomizer. After training the SVM with the BO, we
check the performance of the test classification. The learning process in
general can be optimized within 20 iterations by the BO [3], so we execute 20
iterations. The strategy of the multiclass classification is either the OvO or
the OvA depending on the BO. The classification results are shown in Figure 3
using the EIPS acquisition function.

Figure 3: Classification result of synthetic elastic properties ($AI$ and
$\nicefrac{{V_{P}}}{{V_{S}}}$) using the SVM after 20 BO iterations of the
hyperparameters. The filled circles are the training data, the hollow circles
are the test data, and the crosses are the ground truth. C, B, O, S and L in
the legend box stand for condensate, brine-bearing sand, oil-bearing sand,
shale and limestone facies. The five different background colours correspond
with the classified facies. The details of 20 BO iterations are shown in
Figure 4, and a confusion matrix of the test accuracy is shown in Figure 5.
Figure 4: Cross plots of cost values after 20th BO iteration using the
synthetic data (Figure 3) for: (a) the penalization parameter and sigma; (b)
the cost and the sigma; (c) the penalization parameter and the cost; (d) the
number of iterations and the minimum cost. The stars and arrow correspond to
where the optimal values are found.

The cost function in Figure 4 shows that the optimized hyperparameters are
consequently found in 17th iteration, but major improvement of the cost is
obtained in 8th iteration. From Figures 4a-c we find out that the BO searches
the hyperparameters from both aspects of the exploitation and the exploration
followed by equation 18. This means that several iterations (represented by
circles in Figure 4a-c) are focused on neighbouring locations (e.g., where the
penalization parameters are higher than about 2*100 in Figure 4c) by the
exploitation factor whilst some iterations are randomly realized by the
exploration factor. Based on this experiment, we do not observe major
technical issues caused by our implementation. In Figure 5, we plot a
confusion matrix for the test accuracy of 90.8 %. The numbers in each bin
indicate the test accuracy per target class when the sample numbers in the
identical bin are used. Note that the accuracy shown here stands for the total
accuracy of all kinds of used facies. In this synthetic example, two classes
of oil and gas bearing sands are well separated and distinguished from the
others as the main target. A prediction accuracy of 90.0 % is achieved with
only 20 iterations for the utilized BO.

Figure 5: Confusion matrix of the synthetic data classification in Figure 3.
The percentile in each bin is the test accuracy between the target- and
resulted classes. The integer shown in each bin is the used sample number.

##  4 Results and Discussions

The Figure 6 shows the cross-plot of the elastic parameters, $AI$ and
$\nicefrac{{V_{P}}}{{V_{S}}}$, for the Palaeocene formations derived from ten
chosen wells. There are always overlaps between different facies/classes and
we intend to improve and maximize the separation between the classes by
implementing the methods described in the sections above. The initial
classification results show the optimized parameters are consequently found in
5th iteration out of 20 iterations (Figure 7). The test (prediction) accuracy
after 20 iterations result in 66.2 % at this realization (Figure 8a). Note
that the test accuracy might differ slightly after running consequent
realizations due to the dynamic nature of the BO.

Figure 6: Facies classification result using elastic properties ($AI$ and
$\nicefrac{{V_{P}}}{{V_{S}}}$) from the chosen 10 wells through the SVM after
20th BO iteration of the hyperparameters. The legend and the background
colours are same as Figure 3. The details of 20 BO iterations are shown in
Figure 7, and a confusion matrix of the test accuracy is shown in Figure 8a.
Figure 7: Same as Figure 4, but for the field dataset (Figure 6). Figure 8:
Confusion matrices of the field data classification with: (a) 20 iterations in
Figures 6 and 7; (b) 40 iterations in Figures 9 and 10; (c) 200 iterations in
Figures 11 and 12.

The Figure 7 indicates the sigma are intensively exploited between 100and 10-2
in Figure 7b, where the cost is relatively low among all of the 20 iterations.
On the contrary, we also see that the BO tries to explore where there are
little observations available i.e., the sigma around 10-4 in Figure 7b.

The decision boundaries in the initial classification in general appear to be
good visually (Figure 6), however, there are rooms to improve the separation
between the facies. For instances, it looks like that an alternative boundary
between the brine sand and the limestone can improve the separation rate such
that the limestone class with the white background in Figure 6 could have been
more expanded towards left by exchanging with the brine-sand class with the
dark grey background. On the contrary, it is too hard to uphold the decision
boundaries between classes with significant overlaps such as the oil- and
condensate-bearing sands or brine sand and shale visually. However, any
further improvement in other scenarios can be examined quantitatively by the
confusion matrix, i.e., Figure 8a. A certain number of miss-classified plots
is inevitable because of faceis overlaps on $AI$-$\nicefrac{{V_{P}}}{{V_{S}}}$
cross-plot. When we come to think of such over-lapping cases in general, our
subjectivity might prone to be dictated unless automatic methods like the BO-
based-SVM are adopted.

We investigate further improvement in classification and alternative decision
boundaries by increasing the number of iteration; we doubled it to be 40
times. The results are shown in Figures 8b-10.

Figure 9: Same as Figure 6 but for 40 iterations. The details of 40 BO
iterations are shown in Figure 10, and a confusion matrix of the test accuracy
is shown in Figure 8b. Figure 10: Same as Figure 7 but for 40 iterations
(Figure 9).

Although the improvement is consequently found in the 38th iteration, the
large improvement of the cost function is rapidly gained within the first
several iterations (Figure 10d). The searching range of the hyperparameters in
Figure 10 gets broader for both small and large values than those in Figure 7.
A possible correlation between the sigma and the penalization parameter in
terms of the cost begins to appear in Figure 10a. The most of pairs with the
cost show lower than 0.4 to be likely emerged when both the sigma and the
penalization parameter have similar values in logarithmic scale. This is not
surprising when we consider of the hyperparameters associated with the SVM.
Recall that the sigma is responsible to the margins’ non-linearity whereas the
penalization parameter is related to the misclassification. Hence, it can be
interpreted that the BO automatically responds to the trade-off nature of
these two hyperparameters.

There are certain improvements in the new classification scenario in Figure 9
in comparison with Figure 6 such as visible changes in the boundary between
the brine sand and the limestone. Not only the total test accuracy but also
the accuracy of the oil-bearing sand in the new classification increases by 1
% (Figures 8a-b). Since we still visually see the potential improvement of the
aforementioned boundary between the brine sand and the limestone, we carry out
one more additional test with 200 iterations. The results are shown in Figures
8c, 11 and 12.

Figure 11: Same as Figure 6 but for 200 iterations. The details of 200 BO
iterations are shown in Figure 12, and a confusion matrix of the test accuracy
is shown in Figure 8c. Figure 12: Same as Figure 7, but for 200 iterations
(Figure 11).

Like other last two realizations by the 20- (Figure 7d) and the 40-iteration
(Figure 10d), the 200-iteration (Figure 12d) also shows a major improvement of
the cost function within the first several BO iterations. The possible
correlation discussed in Figure 10a becomes more prominent (Figure 12a). In
overall, the 200-iteration yields more than 1 % improvement of the test
accuracy. For the oil-bearing sand, the accuracy improves to 55.2 %, which is
about 15 % higher than when the 40-iteration is tested.

As demonstrated hitherto, the hyperparameters of the SVM can be automatically
optimized within the framework of the BO. When the cost function appears to be
even stable (e.g., Figure 7d), the test accuracy might be improved further by
simply increasing the iteration if the current classification is far from the
global minima. The iteration number might be determined by a trade-off between
the targeted test accuracy and the computational cost. For the field dataset
used in the study, there is a fundamental limitation for the classifications
relating to the overlapping facies or classes. Some additional elastic
properties as input might be considered in future investigations in order to
achieve better classification beyond the current limitation, which also would
cost more computations. With this regard, so-called multidimensional scaling
[e.g., 1] could be useful to be implemented.

The computational cost by the grid search tends to be enormous especially when
the searching range is large but utilizing the BO is a versatility choice for
geophysical applications as long as we use conventional computers rather than
quantum computers which are under development.

##  5 Conclusions

We presented an application of Bayesian optimization (BO) with support vector
machine (SVM) in order to classify cross-plot products of elastic properties
which could be used in seismic QI studies. The applications showed in this
study used the data from 10 wells in East Central Graben in the UKCS. We
achieved the classification results with visually feasible appearance having
non-linear decision boundaries without choosing the SVM hyperparameters
subjectively. The cost function related to hyperparameters optimization by the
BO exhibited to be stable before finishing the first dozen iterations, but the
test accuracy and the visual classification slightly improved by increasing
the total iteration number. However, it is difficult to predetermine which
iteration number is optimal because there is a trade-off between the
prediction accuracy and the computational cost. Nonetheless, the BO
application is potentially more efficient in optimizing the hyperparameters of
the SVM than a conventional exhaustive-grid search. So, the application of BO
with SVM can be useful for well and consequently seismic quantitative
interpretations and facies analysis in resource exploration and development
subsurface studies.

### Acknowledgments

For SVM applications and Bayesian optimization, LIBSVM and BayesOpt were used
in this study, respectively.

## References

  * Borg and Groenen [1997] I. Borg and P. J. F. Groenen.  _Modern Multidimensional Scaling_.  Springer, 1997.  URL https://link.springer.com/book/10.1007/978-1-4757-2711-1.
  * Buland et al. [1996] A. Buland, M. Landrø, M. Andersen, and T. Dahl.  AVO inversion of Troll field data.  _Geophysics_ , 61, 1996.  doi: 10.1190/1.1444078.
  * Czarnecki et al. [2015] W. M. Czarnecki, S. Podlewska, and A. J. Bojarski.  Robust optimization of SVM hyperparameters in the classification of bioactive compounds.  _Journal of Cheminformatics_ , 7, 2015.  doi: 10.1186/s13321-015-0088-0.
  * Desautels et al. [2014] T. Desautels, A. Krause, and J. W. Burdick.  Parallelizing exploration-exploitation tradeoffs in gaussian process bandit optimization.  _Journal of Machine Learning Research_ , 15, 2014.  URL https://dl.acm.org/doi/10.5555/2627435.2750368.
  * Galar et al. [2011] M. Galar, A. Fernández, E. Barrenechea, H. Bustince, and F. Herrera.  An overview of ensemble methods for binary classifiers in multi-class problems: experimental study on one-vs-one and one-vs-all schemes.  _Pattern Recognition_ , 44, 2011.  doi: 10.1016/j.patcog.2011.01.017.
  * Harrison and Urosevic [2012] C. B. Harrison and M. Urosevic.  Seismic processing, inversion, and avo for gold exploration — case study from western australia.  _Geophysics_ , 77, 2012.  doi: 10.1190/geo2011-0506.1.
  * Jonas et al. [1978] M. Jonas, T. Vytautas, and Z. Antanas.  _The Application of Bayesian Methods for Seeking the Extremum_.  Elsevier, 1978.  URL https://www.researchgate.net/publication/248818761_The_application_of_Bayesian_methods_for_seeking_the_extremum.
  * Jones et al. [1998] D. R. Jones, M. Schonlau, and W. J. Welch.  Efficient global optimization of expensive black-box functions.  _Journal of Global Optimization_ , 13, 1998.  doi: 10.1023/A:1008306431147.
  * Kushner [1964] H. J. Kushner.  A new method of locating the maximum point of an arbitrary multipeak curve in the presence of noise.  _Journal of Basic Engineering_ , 86, 1964.  doi: 10.1115/1.3653121.
  * Scholkopf et al. [1997] B. Scholkopf, S. Kah-Kay, C. Burges, F. Girosi, P. Niyogi, T. Poggio, and V. Vapnik.  Comparing support vector machines with gaussian kernels to radial basis function classifiers.  _IEEE Transactions on Signal Processing_ , 45, 1997.  doi: 10.1109/78.650102.
  * Shahriari et al. [2016] B. Shahriari, K. Swersky, Z. Wang, R. P. Adams, and N. de Freitas.  Taking the human out of the loop: a review of Bayesian optimization.  _Proceedings of the IEEE_ , 106, 2016.  doi: 10.1109/JPROC.2015.2494218.
  * Snoek et al. [2012] J. Snoek, H. Larochelle, and R. P. Adams.  Practical Bayesian optimization of machine learning algorithms.  In _Advances in Neural Information Processing Systems_ , page 2951–2959. NIPS, 2012.  URL https://dl.acm.org/doi/10.5555/2999325.2999464.
  * Srinivas et al. [2010] N. Srinivas, A. Krause, S. M. Kakade, and M. Seeger.  Gaussian process optimization in the bandit setting: no regret and experimental design.  In _Proceedings of the 27th International Conference on International Conference on Machine Learning_ , page 1015–1022. ICML, 2010.  URL https://dl.acm.org/doi/10.5555/3104322.3104451.
  * Vapnik [1995] V. Vapnik.  _The Nature of Statistical Learning Theory_.  Springer, 1995.  URL https://link.springer.com/book/10.1007/978-1-4757-2440-0.
  * Vapnik and Lerner [1963] V. Vapnik and A. Lerner.  Pattern recognition using generalized portrait method.  _Automation and Remote Control_ , 24, 1963.  URL https://ci.nii.ac.jp/naid/10020952249/.
  * Wang et al. [2014] G. Wang, T. R. Carr, Y. Ju, and C. Li.  Identifying organic-rich Marcellus shale lithofacies by support vector machine classifier in the Appalachian basin.  _Computers & Geosciences_, 64, 2014.  doi: 10.1016/j.cageo.2013.12.002.
  * Wang et al. [2022] H. Wang, Z. Wang, J. Ma, L. Li, Y. Wang, M. Tan, Y. Zhang, S. Cui, and Z. Qu.  Effective pressure prediction from 4D seismic AVO data during CO2-EOR and storage.  _International Journal of Greenhouse Gas Control_ , 113, 2022.  doi: 10.1016/j.ijggc.2021.103525.


11institutetext: 1 Shenzhen Key Laboratory of Smart Healthcare Engineering,
Department of Biomedical Engineering, Southern University of Science and
Technology
2 Department of Computer Science, Southern University of Science and
Technology, Shenzhen, 518055, China
11email<EMAIL_ADDRESS>
# Partial Least Square Regression via Three-factor SVD-type Manifold
Optimization for EEG Decoding

Wanguang Yin 11 Zhichao Liang 11 Jianguo Zhang 22 Quanying Liu 11

###### Abstract

Partial least square regression (PLSR) is a widely-used statistical model to
reveal the linear relationships of latent factors that comes from the
independent variables and dependent variables. However, traditional methods
to solve PLSR models are usually based on the Euclidean space, and easily
getting stuck into a local minimum. To this end, we propose a new method to
solve the partial least square regression, named PLSR via optimization on bi-
Grassmann manifold (PLSRbiGr). Specifically, we first leverage the three-
factor SVD-type decomposition of the cross-covariance matrix defined on the
bi-Grassmann manifold, converting the orthogonal constrained optimization
problem into an unconstrained optimization problem on bi-Grassmann manifold,
and then incorporate the Riemannian preconditioning of matrix scaling to
regulate the Riemannian metric in each iteration. PLSRbiGr is validated with a
variety of experiments for decoding EEG signals at motor imagery (MI) and
steady-state visual evoked potential (SSVEP) task. Experimental results
demonstrate that PLSRbiGr outperforms competing algorithms in multiple EEG
decoding tasks, which will greatly facilitate small sample data learning.

###### Keywords:

Partial least square regression, Bi-Grassmann manifold, Riemannian
preconditioning, matrix manifold optimization, EEG decoding.

##  1 Introduction

Extracting the latent factor is an essential procedure for discovering the
latent space of high-dimensional data; thereby proposed lots of regression
model for latent semantic or variable analysis. Among of them, partial least
squares regression (PLSR) is a well-established model for learning the latent
variables, by a sequential way to learn the latent sub-spaces while maximally
maintaining the correlations between the latent variables of independent
variables $X$ and dependent variables $Y$. Specifically, it can be described
as the projection of two variables $X$ and $Y$ onto the lower dimensional sub-
spaces $U$ and $V$. This application can be found in a various of fields,
including chemometrics [2, 8], chemical process control [7, 13], and
neuroscience [5, 9].

As a result, there emerges a lot of partial least square regression models,
while most of them are based on the Euclidean space to solve the latent
factors and identify the latent factors column by column in each iteration. A
drawback of this approach is that it easily converges to a spurious local
minimum and hardly obtain a globally optimal solution. Given an extension of
PLSR to the N-way tensors with rank-one decomposition, which can provide an
improvement on the intuitive interpretation of the model [3]. However, N-way
PLSR suffers from high computational complexity and slow convergence in
dealing with complex data.

To address problems above, we propose a three-factor SVD-type decomposition of
PLSR via optimization on bi-Grassmann manifold (PLSRbiGr). Hence, its
corresponding sub-spaces $U$ and $V$ associated with independent variable $X$
and dependent variable $Y$ can be solved by using a nonlinear manifold
optimization method. Moreover, we leverage the Riemannian preconditioning of
matrix scaling $S^{T}S$ to regulate the Riemannian metric in each iteration
and self-adapt to the changing of subspace [10]. To validate the performance
of PLSRbiGr, we conduct two EEG classification tasks on the motor imagery (MI)
and steady-state visual evoked potential (SSVEP) datasets. The results
demonstrate PLSRbiGr is superior to the inspired modification of PLSR
(SIMPLSR) [6], SIMPLSR with the generalized Grassmann manifold (PLSRGGr),
SIMPLSR with product manifold (PLSRGStO), as well as the sparse SIMPLSR via
optimization on generalized Stiefel manifold (SPLSRGSt) [4]. Our main
contributions can be summarized as follow:

  * •

We propose a novel method for solving partial least square regression (PLSR),
named PLSR via optimization on bi-Grassmann manifold (PLSRbiGr), which
decomposes the cross-covariance matrix ($X^{T}Y$) to an interpretable subspace
($U$ and $V$) and simultaneously learns the latent factors via optimization on
bi-Grassmann manifold (Sec 2.2).

  * •

We present Riemannian metric that equipped with Riemannian preconditioning, to
self-adapt to the changing of subspace. Fortunately, Riemannian
preconditioning can largely improves the algorithmic performance (i.e.
convergence speed and classification accuracy) (Sec 3.2).

  * •

The results of EEG decoding demonstrate that PLSRbiGr outperforms conventional
Euclidean-based methods and Riemannian-based methods for small sample data
learning (Sec 3.1).

##  2 Method

###  2.1 Review of Partial Least Squares Regression

Partial least square regression (PLSR) is a wide class of models for learning
the linear relationship between independent variables $X\in\mathbb{R}^{I\times
N}$ and dependent variables $Y\in\mathbb{R}^{I\times M}$ by means of latent
variables. We present the schematic illustration of partial least square
regression (PLSR) in Fig. 1.

To predict dependent variables $Y$ from independent variables $X$, PLSR finds
a set of latent variables (also called the latent vectors, score vectors, or
components) by projecting both $X$ and $Y$ onto a lower dimensional subspace
while at the same time maximizing the pairwise covariance between the latent
variables $T$ and $B$, which are presented in Eq. (1) & Eq. (2).

Figure 1: The PLSR decomposes the independent variables (EEG)
$X\in\mathbb{R}^{I\times N}$ and dependent variables (sample labels)
$Y\in\mathbb{R}^{I\times M}$ as a sum of rank-one matrices.

$\small X=TP^{T}+E=\sum_{r=1}^{R}t_{r}p_{r}^{T}+E$ (1)

$\small Y=BQ^{T}+F=\sum_{r=1}^{R}b_{r}q_{r}^{T}+F$ (2)

where $T=\left[t_{1},t_{2},\ldots,t_{R}\right]\in\mathbb{R}^{I\times R}$ is a
matrix of $R$ extracted latent variables from $X$, and
$B=\left[b_{1},b_{2},\ldots,b_{R}\right]\in\mathbb{R}^{I\times R}$ are latent
variables from $Y$, their columns have the maximum correlations between each
other. In addition, $P\in\mathbb{R}^{N\times R}$ and $Q\in\mathbb{R}^{M\times
R}$ are the loading matrices, and $E\in\mathbb{R}^{I\times N}$ and
$F\in\mathbb{R}^{I\times M}$ are the residuals with respect to $X$ and $Y$.

However, most of the current methods for solving PLSR are based on the
Euclidean space by performing the sum of a minimum number of rank-one
decomposition to jointly approximate the independent variables $X$ and
dependent variables $Y$. No existing method solves PLSR with bi-Grassmann
manifold optimization. To this end, we propose a novel method for solving PLSR
via optimization on the bi-Grassmann manifold.

####  2.1.1 Regression model

To predict the dependent variable given a new sample $X_{new}$, its regression
model is given by:

$Y^{new}=X_{new}C=X_{new}W\left(P^{T}W\right)^{-1}DQ^{T}$ (3)

where $Y^{new}$ is the predicted output, and
$C=U\left(P^{T}U\right)^{-1}DQ^{T}$ is the regression coefficient calculated
from the training data. Following to [4], we can obtain the learned sub-spaces
$U$ & $V$, and its corresponding latent factors $T=XU$ & $B=YV$ by solving the
SVD-type model via optimization on bi-Grassmann manifold.

###  2.2 PLSRbiGr: PLSR via Optimization on Bi-Grassmann Manifold

In the sequel, we present the SVD-type decomposition via optimization on bi-
Grassmann manifold. As shown in Fig. 2, we treat the cross-product matrix (or
tensor) of independent variables and dependent variables as the input data,
and then perform the three factor decomposition via optimization on bi-
Grassmann manifold.

Figure 2: SVD-type decomposition via optimization on bi-Grassmann manifold.

####  2.2.1 Three-factor SVD-type decomposition

The SVD-type decomposition of cross-product matrix (i.e. $Z$) can be
formulated as the following tensor-matrix product:

$\small
f\left(U,V,S\right)=\mathop{\mathrm{argmin}}_{S,U,V}\frac{1}{2}\left\|S\times_{1}U\times_{2}V-Z\left(k\right)\right\|_{F}^{2}$
(4)

where $k$ denotes the count of iterations. To apply the Riemannian version of
conjugate gradient descent or trust-region method, it needs to obtain the
partial derivatives of $f$ with respect to $U$, $V$ and $S$, that is given by

$\displaystyle\mathrm{grad}f\left(U,V,S\right)$ $\displaystyle=\frac{\partial
f\left(U,V,S\right)}{\partial U}$ (5)
$\displaystyle=\left(USV^{T}-Z\left(k\right)\right)VS^{T}$
$\displaystyle=Z\left(k+1\right)VS^{T}$

$\displaystyle\mathrm{grad}f\left(U,V,S\right)$ $\displaystyle=\frac{\partial
f\left(U,V,S\right)}{\partial V}$ (6)
$\displaystyle=\left(USV^{T}-Z\left(k\right)\right)^{T}US^{T}$
$\displaystyle=Z\left(k+1\right)^{T}US^{T}$

$\displaystyle\mathrm{grad}f\left(U,V,S\right)$ $\displaystyle=\frac{\partial
f\left(U,V,S\right)}{\partial S}$ (7)
$\displaystyle=\left(S\times_{1}U\times_{2}V-Z\left(k\right)\right)\times_{1}U^{T}\times_{2}V^{T}$
$\displaystyle=U^{T}\left(USV^{T}-Z\left(k\right)\right)V$
$\displaystyle=U^{T}Z\left(k+1\right)V$

Then, we project the Euclidean gradient (i.e. Eq.(5), Eq.(6), and Eq.(7)) onto
the Riemannian manifold space.

####  2.2.2 SVD-type decomposition via optimization on bi-Grassmann manifold

In this subsection, we derive the bi-Grassmann manifold optimization for the
SVD-type decomposition and $\mathcal{M}$ is a manifold equipped with
Riemannian metric. Recall that Stiefel manifold is the set of matrices whose
columns are orthogonal, that is denoted by

$\small St\left(N,R\right)=\left\\{U\in\mathbb{R}^{N\times
R}|U^{T}U=\mathrm{I}_{R}\right\\}$ (8)

For a Stiefel manifold $St\left(N,R\right)$, its related Grassmann manifold
$Gr\left(N,R\right)$ can be formulated as the quotient space of
$St\left(N,R\right)$, under the equivalence relation defined by the orthogonal
group,

$\small Gr\left(N,R\right)=St\left(N,R\right)/\mathcal{O}\left(R\right)$ (9)

here, $\mathcal{O}\left(R\right)$ is the orthogonal group defined by

$\small\mathcal{O}\left(R\right)=\left\\{W\in\mathbb{R}^{R\times
R}|W^{T}W=WW^{T}=\mathrm{I}_{R}\right\\}$ (10)

Moreover, for the SVD-type decomposition, the optimization sub-spaces can be
expressed as the following bi-Grassmann manifold:

$\small\mathcal{M}:=Gr\left(N,R\right)\times
Gr\left(M,R\right)\times\mathbb{R}^{R\times R}$ (11)

that is equipped with following Riemannian metric:

$\small\mathrm{g}_{Z}\left(\xi_{Z},\eta_{Z}\right)=tr\left(SS^{T}\xi_{U}^{T}\eta_{U}\right)+tr\left(S^{T}S\xi_{V}^{T}\eta_{V}\right)+tr\left(\xi_{S}^{T}\eta_{S}\right)$
(12)

In practice, the computational space (i.e. $Z\in\mathbb{R}^{N\times M}$) can
be first decomposed into orthogonal complementary sub-spaces (normal space
i.e. $\mathrm{N}_{Z}\mathcal{M}$, and tangent space i.e.
$\mathrm{T}_{Z}\mathcal{M}$), and then the tangent space can be further
decomposed into the other two orthogonal complementary sub-spaces (horizontal
space i.e. $\mathrm{H}_{Z}\mathcal{M}$, and vertical space i.e.
$\mathrm{V}_{Z}\mathcal{M}$), and eventually we project the Euclidean gradient
to the horizontal space defined by the equivalence relation of orthogonal
group [1].

####  2.2.3 Riemannian Gradient

To obtain the Riemannian gradient, it needs to project the Euclidean gradient
onto the Riemannian manifold space [1], that is

$\displaystyle\mathrm{Grad}f\left(U,V,S\right)$
$\displaystyle=\Pi_{U}\mathrm{grad}f\left(U,V,S\right)$ (13)
$\displaystyle=\left(\mathrm{grad}f\left(U\right)-UB_{U}\right)\left(SS^{T}\right)^{-1}$

$\displaystyle\mathrm{Grad}f\left(U,V,S\right)$
$\displaystyle=\Pi_{V}\mathrm{grad}f\left(U,V,S\right)$ (14)
$\displaystyle=\left(\mathrm{grad}f\left(V\right)-VB_{V}\right)\left(S^{T}S\right)^{-1}$

where $\left(SS^{T}\right)^{-1}$ and $\left(S^{T}S\right)^{-1}$ is the scaling
factors. $\Pi_{U}=\mathrm{P}_{U}^{h}\mathrm{P}_{U}^{v}$ is the project
operator that involves two steps of operation, one is a mapping from ambient
space to the tangent space $\left(\textit{i.e.,}\mathrm{P}_{U}^{v}\right)$,
and the other is a mapping from tangent space to the horizontal space
$\left(\textit{i.e.,}\mathrm{P}_{U}^{h}\right)$. Therefore, the computational
space of Stiefel manifold $St\left(N,M\right)$ can be decomposed into tangent
space and normal space, and the tangent space can be further decomposed into
two orthogonal complementary sub-spaces (i.e. horizontal space and vertical
space) [10]. Once the expression of Riemannian gradient and Riemannian Hessian
are obtained, we can conduct the Riemannian manifold optimization by using
Manopt toolbox [1].

##  3 Experiments and Results

To test the performance of our proposed algorithm, we conduct experiments on a
lot of EEG signal decoding, whose performance is compared to several well
known algorithms, including the statistically inspired modification of PLSR
(SIMPLSR), SIMPLSR with the generalized Grassmann manifold (PLSRGGr), sparse
SIMPLSR via optimization on generalized Stiefel manifold (SPLSRGSt) [4, 6],
and higher order partial least squares regression (HOPLSR) [12].

###  3.1 EEG Decoding

In this subsection, we test the efficiency and accuracy of our proposed
algorithm (PLSRbiGr) on the public PhysioNet MI dataset [11]. We compare
PLSRbiGr with other existing algorithms, including the statistically inspired
modification of PLSR (SIMPLSR), SIMPLSR with the generalized Grassmann
manifold (PLSRGGr), sparse SIMPLSR via optimization on generalized Stiefel
manifold (SPLSRGSt) [4, 6], and higher order partial least squares (HOPLSR)
[12]. To evaluate the performance of decoding algorithms, we use Accuracy
(Acc) as the evaluation metric to quantify the performance of comparison
algorithms.

Figure 3: Accuracy of 2-class MI classification task on PhysioNet MI dataset

In training, we set the 4-fold-cross-validation to obtain the averaged
classification accuracy in testing samples. As shown in Fig. 3, PLSRbiGr
generally achieves the best performance in comparison to the existing methods.
The used PhysioNet EEG MI dataset consists of 2-class MI tasks (i.e. runs 3,
4, 7, 8, 11, and 12, with imagine movements of left fist or right fist) [11],
which is recorded from 109 subjects with 64-channel EEG signals (sampling rate
equals to 160 Hz) during MI tasks. We randomly select 10 subjects from
PhysioNet MI dataset in our experiments. The EEG signals are filtered with a
band-pass filter (cutoff frequencies at $7\sim 35$ Hz) and a spatial filter
(i.e. Xdawn with 16 filters), therefore the resulting data is represented by
$trials\times channel\times time$.

Figure 4: Accuracy of 4-class SSVEP classification task on Macau SSVEP dataset

The used Macau SSVEP dataset contains 128-channel EEG recordings from 7
subjects sampled at 1000 Hz, which was recorded at University of Macau with
their ethical approval. There are four types visual stimulus with flashing
frequency at 10Hz, 12Hz, 15Hz, and 20Hz. To increase the signal-noise ratio
(SNR) and reduce the dimension of raw EEG signals, the EEG signals are
filtered with a band-pass filter (cutoff frequencies at $7\sim 35$ Hz) and a
spatial filter (i.e. Xdawn with 16 filters). The epochs of 1-second EEG
signals before every time point of the SSVEP data were extracted and down-
sampled to 200 Hz, therefore the resulting data is represented by $trial\times
channel\times time$. Fig. 4 presents the classification accuracy of all
comparison algorithms on the SSVEP dataset.

###  3.2 Effects of Riemannian Preconditioning

Furthermore, we test the effects of Riemannian preconditioning. In PLSRbiGr,
the scaling factors $S^{T}S$ and $SS^{T}$ of Riemannian preconditioning
provides an effective strategy to accelerate the convergence speed. As shown
in Table 1, the classification accuracy and running time of PLSRbiGr equipped
with Riemannian preconditioning are closely better than results of other
methods, such as PLSRGGr, PLSRGStO, and SPLSRGSt that have not taken into
account of Riemannian preconditioning.

Table 1: The effects of Riemannian preconditioning on the classification
accuracy and running time. The test experiments were conducted on Macau SSVEP
dataset by using PLSRbiGr.

ID Accuracy Running time (s) preconditioned non-preconditioned preconditioned
non-preconditioned S001 0.8487±0.0148 0.3063±0.0394 0.0571 0.0570 S002
0.7036±0.0123 0.1211±0.0241 0.0624 0.0544 S003 0.5903±0.0183 0.1854±0.0278
0.0576 0.0551 S004 0.7360±0.0231 0.1056±0.0258 0.0960 0.0591 S005
0.8026±0.0182 0.3187±0.0483 0.0692 0.0560 S006 0.7987±0.0075 0.1522±0.0285
0.0578 0.0574 S007 0.8375±0.0124 0.2135±0.0375 0.0593 0.0557 mean
0.7596±0.0152 0.2004±0.0330 0.0656 0.0563

##  4 Discussion and Conclusion

In this paper, we propose a novel method, named partial least square
regression via optimization on bi-Grassmann manifold (PLSRbiGr) for EEG signal
decoding. It features to find the objective solution via optimization on bi-
Grassmann manifold. Specifically, to relax the orthogonality constraints of
objective function, PLSRbiGr converts the constrained optimization problem in
Euclidean space to an optimization problem defined on bi-Grassmann manifold,
thereby its corresponding subspaces (i.e. $U$ and $V$) can be learned by using
Riemannian manifold optimization instead of the traditional methods that by
deflating the residuals in each iteration. In practice, PLSRbiGr can also be
used for image classification, and many other prediction tasks, which has no
such limitations. More importantly, extensive experiments on MI and SSVEP
datasets suggest that PLSRbiGr robustly outperforms other methods optimized in
Euclidean space and has a fast convergence than other algorithms that solved
in Riemannian manifold space (Table 1), and the scaling factor of Riemannian
preconditioning provides a good generalization ability between subjects and
robust to variances (Table 1).

A limitation of our method is that the column of cross-product matrix equals
to the class of training samples, thereby the low-rank nature of PLSRbiGr is a
certain value. To address such problem, several different directions can be
carried out in our feature work. For example, we only consider the case of
cross-product matrix, when cross covariance is a multi-way array (tensors), it
needs to further consider the product manifold optimization. For example, the
higher order partial least squares (HOPLSR) can also be defined on the
Riemannian manifold space [12], and directly optimized over the entire product
of $N$ manifolds, thereby the rank of cross-product tensor can be
automatically inferred from the optimization procedures. Another issue that
needs to be further investigated is how to scale the computation of Riemannian
preconditioning to the higher order partial least squares (HOPLSR).

## References

  * [1] P.-A. Absil, R. Mahony, and R. Sepulchre.  Optimization algorithms on matrix manifolds.  Princeton University Press, 2009.
  * [2] R. G. Brereton and G. R. Lloyd.  Partial least squares discriminant analysis for chemometrics and metabolomics: H ow scores, loadings, and weights differ according to two common algorithms.  Journal of Chemometrics, 32(4):e3028, 2018.
  * [3] R. Bro.  Multiway calibration. multilinear pls.  Journal of chemometrics, 10(1):47–61, 1996.
  * [4] H. Chen, Y. Sun, J. Gao, Y. Hu, and B. Yin.  Solving partial least squares regression via manifold optimization approaches.  IEEE transactions on neural networks and learning systems, 30(2):588–600, 2018.
  * [5] Y. Chu, X. Zhao, Y. Zou, W. Xu, G. Song, J. Han, and Y. Zhao.  Decoding multiclass motor imagery eeg from the same upper limb by combining riemannian geometry features and partial least squares regression.  Journal of Neural Engineering, 17(4):046029, 2020.
  * [6] S. De Jong.  Simpls: an alternative approach to partial least squares regression.  Chemometrics and intelligent laboratory systems, 18(3):251–263, 1993.
  * [7] Y. Dong and S. J. Qin.  Regression on dynamic pls structures for supervised learning of dynamic data.  Journal of Process Control, 68:64–72, 2018.
  * [8] K. Hasegawa, M. Arakawa, and K. Funatsu.  Rational choice of bioactive conformations through use of conformation analysis and 3-way partial least squares modeling.  Chemometrics and Intelligent Laboratory Systems, 50(2):253–261, 2000.
  * [9] D. A. Hoagey, J. R. Rieck, K. M. Rodrigue, and K. M. Kennedy.  Joint contributions of cortical morphometry and white matter microstructure in healthy brain aging: a partial least squares correlation analysis.  Human brain mapping, 40(18):5315–5329, 2019.
  * [10] H. Kasai and B. Mishra.  Low-rank tensor completion: a riemannian manifold preconditioning approach.  In International Conference on Machine Learning, pages 1012–1021. PMLR, 2016.
  * [11] G. Schalk, D. J. McFarland, T. Hinterberger, N. Birbaumer, and J. R. Wolpaw.  Bci2000: a general-purpose brain-computer interface (bci) system.  IEEE Transactions on biomedical engineering, 51(6):1034–1043, 2004\.
  * [12] Q. Zhao, C. F. Caiafa, D. P. Mandic, Z. C. Chao, Y. Nagasaka, N. Fujii, L. Zhang, and A. Cichocki.  Higher order partial least squares (hopls): a generalized multilinear regression method.  IEEE transactions on pattern analysis and machine intelligence, 35(7):1660–1673, 2012.
  * [13] J. Zheng and Z. Song.  Semisupervised learning for probabilistic partial least squares regression model and soft sensor application.  Journal of process control, 64:123–131, 2018.


,

# Effects of the equilibrium model on impurity transport in tokamaks

A Skyman, L Fazendeiro, D Tegnered, H Nordman, J Anderson and P Strand
Department of Earth and Space Sciences, Chalmers University of Technology,
SE-412 96 Gothenburg, Sweden<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Gyrokinetic simulations of ion temperature gradient mode and trapped electron
mode driven impurity transport in realistic tokamak geometry are presented and
compared with results using simplified geometries. The gyrokinetic results,
obtained with the GENE code in both linear and non-linear mode, are compared
with data and analysis for a dedicated impurity injection discharge at JET.
The impact of several factors on heat and particle transport are discussed,
lending special focus to tokamak geometry and rotational shear. To this end,
results using $s-\alpha$ and concentric circular equilibria are compared to
results with magnetic geometry from a JET-experiment. To further approach
experimental conditions, non-linear gyrokinetic simulations are performed with
collisions and a carbon background included.

The impurity peaking factors, computed by finding local density gradients
corresponding to zero particle flux, are discussed. The impurity peaking
factors are seen to be reduced by a factor of $\sim 2$ in realistic geometry
compared to the simplified geometries, due to a reduction of the convective
pinch. It is also seen that collisions reduce the peaking factor for low $Z$
impurities, while increasing it for high charge numbers, which is attributed
to a shift in the transport spectra towards higher wave-numbers with the
addition of collisions. With the addition of roto-diffusion, an over-all
reduction of the peaking factors is observed, but this decrease is not
sufficient to explain the flat carbon profiles seen at JET.

###### pacs:

28.52.Av, 52.25.Vy, 52.30.Ex, 52.30.Gz, 52.35.Ra, 52.55.Fa, 52.65.Tt

##  1 Introduction

Impurity transport is a matter of crucial relevance for Tokamak fusion plasmas
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21,
22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32] due to the contribution of
impurities to radiation losses and plasma fuel dilution. The impurities can
originate either from the sputtering of wall and divertor materials or from
deliberate impurity injection in order to reduce power load. The impurities in
a tokamak thus cover a large range in charge number $Z$, making it necessary
to study the scaling of impurity transport with $Z$. This is even more
relevant after the installation of the new ITER-like wall at JET [33], with
its beryllium wall and tungsten divertor. In addition, a typical plasma
discharge already has such a low particle density, that even a modest amount
of impurities can greatly dilute the plasma, thereby reducing the power
output. Since most impurity sources are located at the edge of the tokamak,
the resulting impurity profile in the core of a given discharge will depend on
the balance between diffusive and advective processes.

It is well known empirically that shaping of the plasma has a beneficial
effect on energy confinement, the Troyon limit and the Greenwald density
limit. For example, elongation ($\kappa$) is one of the factors considered in
empirical scaling laws (e.g. [34]) of the energy confinement time. Only in
recent years, however, has the impact of the geometric configuration on micro-
turbulence been rigorously ascertained, e.g. from simulations of the
gyrokinetic equations [35, 36, 37, 38, 39].

It has been observed that the heat-flux scales with elongation at fixed
average temperature gradients, and that shaping effects enhance residual zonal
flows, thus increasing the nonlinear critical temperature gradient [36].
Elongation has been seen to be the dominant effect on micro-turbulence on ion
temperature gradient (ITG) and trapped electron (TE) drift-wave instabilities,
compared with e.g. triangularity, up–down asymmetry and Dimits shift [38].
However, in spite of these results, the effects of the equilibrium model on
impurity transport has not been studied in any detail.

Previous work in this area includes [21, 22, 23], in which the impurity
transport in a dedicated impurity injection experiment at JET was studied
using $s-\alpha$ geometry. It was shown that the profile steepness parameter
of the injected impurities was qualitatively reproduced by the simulations. On
the other hand, a discrepancy was found for the background carbon profile
which was too steep at mid radius compared to the experimental profile.

In the present study, an experimental equilibrium is used and compared with
results using simpler $s-\alpha$ and circular equilibria [37]. Since the
impurity profile is a result of a balance between diffusion and advection,
where the main advective term is the inward curvature pinch, a strong
dependence on the equilibrium model may be expected. In addition, a more
realistic physics description, including the effects of rotational shear on
the instabilities and impurity transport (roto-diffusion), is investigated
using gyrokinetic simulations. Realistic non-linear gyrokinetic simulations
are performed taking into account collisionality, and the presence of a $2\%$C
background, consistent with the conditions applying for an impurity injection
experiment at JET.

The paper is organised as follows. Section 2 describes the gyrokinetic
simulations performed, and the fundamental concepts of impurity transport.
Section 2.2 describes the particular JET discharge considered. In Section 3,
results for the the eigenvalues, turbulent fluctuation levels and transport
are presented, along with a discussion on the scaling of the impurity peaking
factor with impurity charge $Z$ for the various models considered. A
conclusion then follows in Section 4.

##  2 Background

###  2.1 Gyrokinetic model

The gyrokinetic simulations in this study were performed using the gyrokinetic
turbulence code GENE111http://www.ipp.mpg.de/~fsj/gene/ [40, 41] where a
flux–tube domain was considered. This is an Eulerian–type code (i.e.,
employing a fixed grid in phase space), which allows for arbitrarily long
simulation times. Both quasi- and nonlinear simulations including kinetic ions
and electrons, and at least one impurity species were performed. All
impurities were treated as fully kinetic species with low concentrations. For
studying the effects of collisions, we considered the GENE implementation of
the Landau–Boltzmann type collision operator [41].

For the simulation domain, all three geometries used a flux tube with periodic
boundary conditions in the perpendicular plane. The nonlinear simulations were
performed using a $96\times 96\times 32$ grid in the normal, bi-normal, and
parallel spatial directions respectively; in the perpendicular and parallel
momentum directions, a $12\times 48$ grid was used. For the linear and
quasilinear computations, a typical resolution was $12\times 32$ grid points
in the normal and parallel directions, with $12\times 64$ grid points in
momentum space. To validate the choice of resolutions, convergence tests for
mode structure and spectra were performed.

The nonlinear simulations were typically run up to $t=600\,a/c_{\text{s}}$ for
the experimental geometry scenario, where $a$ represents minor radius and
$c_{\text{s}}=\sqrt{T_{e}/m_{i}}$. For the “full scenario,” in which a $2\%$C
background and collisions were both present, the simulation time was instead
$t\approx 450\,a/c_{\text{s}}$. This shorter time span is due to the much more
computationally intensive nature of the second scenario, by a factor of $\sim
2$–$3$. The increased numerical cost was mainly due to collisions, though in
part also due to the addition of the $2\%$C background.

###  2.2 Experimental profiles and parameters

The physical parameters used in the simulations are presented in table. 1 and
were chosen so as to be consistent with the experimental values taken from the
JET database for discharge #67330. This was previously analysed in [21], where
the $s-\alpha$ geometry model was used. The discharge was part of an impurity-
dedicated set of discharges with low MHD activity [32]. In these, extrinsic
impurities (Ne, Ar and Ni) were injected via laser ablation and gas injection.
Diffusivity $D_{Z}$ and convective velocity $V_{Z}$ were determined by
matching spectroscopic data with results obtained from predictive transport
codes (see [21] and in particular Fig. 4). Profiles of density and temperature
for the discharge are shown in figure 1a.

(a) profiles of $T_{i,e}$ and $n_{i}$

(b) profiles of $v_{\text{tor}}$, $q$ and $s$

Figure 1: (colour online) Radial profiles of the background parameters for JET
discharge #67730 at $t=47.5\,\mathrm{s}$. The simulations were performed at
mid radius ($r/a=0.5$; indicated).

For the present study, the parameters used, including the MHD equilibrium
(from now on referred to as “experimental geometry”), were taken extracted at
mid radius ($r/a=0.5$) at $48\,\mathrm{s}$. The main quantities of interest
were elongation $\kappa=1.37$, triangularity $\delta=0.044$, and the fraction
of trapped particles $f_{t}=0.55$. The toroidal magnetic field intensity on
the axis was $B=3\,\mathrm{T}$, the major radius of the tokamak
$R=3\,\mathrm{m}$, electron temperature $T_{e}=1.55\,\mathrm{k}\,\mathrm{eV}$
and plasma $\beta=8\pi n_{e}T_{e}/B^{2}=1.28\times 10^{-3}$, and the
normalised gradient of $\beta$
$\alpha=-q^{2}R\mathrm{d}\beta/\mathrm{d}r=0.1026$. The safety factor in the
discharge was $q=2.2$ at the considered radius, obtained from EFIT. The
collisionality was $\nu_{ei}=0.07\,c_{\text{s}}/R$ ($\nu_{c}=0.28\times
10^{-3}$ in GENE). Effects of purely toroidal rotation were included through
the $\boldsymbol{E}\times\boldsymbol{B}$ shearing rate. The normalised
gradient scale lengths are defined as
$R/L_{n_{j}}=-(R/n_{j})(\mathrm{d}n_{j}/\mathrm{d}r)$ and
$R/L_{T_{j}}=-(R/T_{j})(\mathrm{d}T_{j}/\mathrm{d}r)$ where $R$ is the major
radius of the tokamak. For the realistic equilibria, the gradients were
rescaled to the minor radius with
$a=\sqrt{2\Phi_{\text{s}}/B}=1.26\,\mathrm{m}$, where $\Phi_{\text{s}}$ is the
toroidal flux at the last closed flux surface. In the considered discharge,
the Mach number was $M=0.21$, leading to
$\gamma_{\boldsymbol{E}\times\boldsymbol{B}}=-\frac{r}{q}\frac{1}{a}\frac{\partial
v_{\text{tor}}}{\partial r}=0.1\,c_{\text{s}}/R$, with the toroidal velocity
profile ($v_{\text{tor}}(r)$) as in figure 1b, evaluated at mid radius. Hence,
we only consider the effect of flow shear in the limit where the flow is
small, neglecting effects of centrifugal and Coriolis forces. These may,
however, be important for heavier impurities [9].

Table 1: Parameters for gyrokinetic simulations $T_{i}/T_{e}$: | $1.02$
---|---
$\beta$: | $1.28\times 10^{-3}$
$s$: | $0.75$
$q$: | $2.20$
$\varepsilon=r/R$: | $0.17$
$\alpha$: | $0.126$
$\kappa$: | $1.37$
$\delta$: | $0.044$
$k_{\theta}\rho_{\text{s}}$: |  $0.2$–$0.6$
$n_{e}$, $n_{i}+Z\,n_{Z}$: | $1.0\cdot 10^{19}\,\mathrm{m^{-}3}$
$n_{Z}$ _(trace)_ : | $10^{-6}\cdot 10^{19}\,\mathrm{m^{-}3}$
$Z$: |  $2$–$74$
$R/L_{n_{i,e}}$: | $2.7$
$R/L_{T_{i}},R/L_{T_{Z}}$: | $5.6$
$R/L_{T_{e}}$: | $5.6$
$\nu_{c}$: |  $0$, $0.28\times 10^{-3}$
$\gamma_{\boldsymbol{E}\times\boldsymbol{B}}$: |  $0$–$0.6$
$T_{e}$: | $1.55\,\mathrm{k}\,\mathrm{eV}$
$B$: | $3.04\,\mathrm{T}$

###  2.3 Magnetic equilibrium

The gyrokinetic simulations were performed in three different geometries:
$s-\alpha$, concentric circular [37] and an MHD equilibrium calculated for the
JET-discharge #67330 from EFIT data using the the TRACER code (see [35] for
details on the geometry implementation). A cross section of the experimental
magnetic equilibrium is shown in figure 2. For the parameters under
consideration, $\alpha$ is fairly small ($\alpha=0.1026$), meaning that the
$s-\alpha$ and circular model should be near equivalent. The $s-\alpha$ model,
however, suffers from inconsistencies of the order $\varepsilon=a/R$ in its
standard implementation [35, 37, 38]. Supplementary simulations using circular
geometry were therefore performed, in order to better differentiate between
effects of magnetic geometry and effects due to the inconsistent $s-\alpha$
metric.

Figure 2: Cross section of the magnetic equilibrium, for JET discharge #67730.

###  2.4 Impurity transport

For trace impurities, equation (4) can be uniquely written [21] as a linear
function of $\nabla n_{Z}$, offset by a convective velocity or “pinch”
$V_{Z}$:

$\Gamma_{Z}=-D_{Z}\nabla
n_{Z}+n_{Z}V_{Z}\Leftrightarrow\frac{R\Gamma_{Z}}{n_{Z}}=D_{Z}\frac{R}{L_{n_{Z}}}+RV_{Z},$
(1)

where $D_{Z}$ is the impurity diffusion coefficient, and $V_{Z}$ is the
impurity convective velocity, which includes roto-diffusion in the gyrokinetic
treatment. Both $D_{Z}$ and $V_{Z}$ are independent of $\nabla n_{Z}$ in the
trace impurity limit [16].

In the core of a steady-state plasma with fuelling from the edge (i.e. in the
absence of particle sinks or sources in the core), the impurity flux
$\Gamma_{Z}$ will go to zero. Since the gradient is a measure of how peaked
the impurity density profile is, the gradient of zero impurity flux is
referred to as the impurity peaking factor ($PF$). Inserting $R/L_{n_{Z}}=PF$
into the linearised equation (1), it can be seen that the peaking factor
quantifies the balance between convective and diffusive impurity
transport:222this definition is closely related to the _Péclet number_ [42,
43]

$PF=-\frac{R\,V_{Z}}{D_{Z}}.$ (2)

Specifically, the sign of the peaking factor is determined by the sign of the
pinch, meaning that $PF>0$ is indicative of a net inward impurity pinch,
giving a peaked impurity profile. Conversely, if $PF<0$ the net impurity pinch
is outward, leading to a hollow impurity profile (also called _flux
reversal_). For trace impurities, equation (1) is linear in $R/L_{n_{Z}}$, and
the peaking factor is found by calculating the impurity particle flux
($\Gamma_{Z}$) for a range of $R/L_{n_{Z}}$ and solving equation (1) for zero
flux. For non-trace species, the $D_{Z}$ and $V_{Z}$ depend on $R/L_{n_{Z}}$,
and the peaking factor has to be found by explicitly seeking the gradient
yielding zero particle flux.

Much of the observed difference between the TE and ITG mode dominated cases
can be understood from the convective velocity $V_{Z}$ in equation (1). To
lowest order in $Z^{-1}$, the pinch contains two terms that depend on $Z$
[16]:

  * •

thermal diffusion (thermopinch):

    * –

$V_{\nabla T_{Z}}\sim\frac{1}{Z}\frac{R}{L_{T_{Z}}}\omega_{r}$,

    * –

inward for TE mode $\left(V_{\nabla T_{Z}}<0\right)$, outward for ITG mode
$\left(V_{\nabla T_{Z}}>0\right)$,

  * •

parallel impurity compression:

    * –

$V_{\parallel,Z}\sim\frac{Z}{A_{Z}}k_{\parallel}^{2}\sim\frac{Z}{A_{Z}q^{2}}\approx\frac{1}{2q^{2}}$,

    * –

outward for TE mode $\left(V_{\parallel_{Z}}>0\right)$, inward for ITG mode
$\left(V_{\parallel_{Z}}<0\right)$.

Here $1/k_{\parallel}$ represents the wavelength of the parallel structure of
the turbulence. Due to the ballooning character of the modes considered, this
is proportional to the safety factor ($q$). In addition, the convective
velocity contains the curvature pinch, which to lowest order is independent of
the impurity charge [21].

The $Z$ dependence in the parallel impurity compression is expected to be
weak, since the mass number is approximately $A_{Z}\approx 2Z$ for an impurity
species with charge $Z$, and for the high $q$ considered, this term is
expected to be small compared to the other pinch contributions. The
thermodiffusive contribution, however, can dominate the transport for low $Z$
impurities (such as the Helium ash). For lower charge numbers, the second
order correction to the thermal pinch can become important:

$V_{T}\sim\left(\omega_{r}\frac{T_{Z}}{T_{e}Z}-\frac{7}{4}\left(\frac{T_{Z}}{T_{e}Z}\right)^{2}\right)\frac{R}{L_{T_{Z}}},$
(3)

where the second order term is inward, independent of the mode direction, and
finite Larmor-radius effects have been neglected [44, 45, 16, 21].

The direction of the contributions to the pinch are governed mainly by the
considered mode’s direction of rotation, which in the considered cases has a
different sign for TE and ITG modes [45]. Expanding the convective velocity in
(1) into thermal diffusion and “pure convection”,

$\Gamma_{Z}=-D_{Z}\nabla n_{Z}+n_{Z}\left(D_{T,Z}\nabla T_{Z}+V_{p,Z}\right),$
(4)

their relative contributions to the total peaking factor can be uniquely
determined for trace impurities, using the method described in e.g. [46]. Here
$V_{p,Z}$ includes the parallel compression and curvature pinches. In the
presence of sheared toroidal rotation ($v_{\phi}$), these terms are joined by
a roto-diffusive term, which is proportional to
$\mathrm{d}v_{\phi}/\mathrm{d}r$.

##  3 Results and discussion

###  3.1 Background turbulence and transport

####  3.1.1 Linear eigenvalue spectra

Using GENE, the linear eigenvalues in the studied geometries were calculated
at mid-radius of JET L-mode discharge #67330 for a range of different wave
numbers. The resulting spectra are compared in figure 3. The results in figure
3a show an overall destabilisation and shift toward higher
$k_{\theta}\rho_{\text{s}}$ when moving from $s-\alpha$ to more realistic
geometries, consistent with the results reported in [37, 47]. For
triangularity of the same order as the experimental value ($\delta=0.044$),
the effect on the eigenvalues of $\delta$ was found to be minute, and it is
therefore surmised that the elongation is the main factor behind the the
observed results.

Adding more realism in the form of collisions and a $2\%$ carbon background,
in accordance with the “full scenario” considered in the nonlinear runs, both
have a stabilising effect, as can be seen in figure 3b. The figure also shows
that the stabilising effect is larger for the addition of collisions than for
the carbon background, in particular for lower values of
$k_{\theta}\rho_{\text{s}}$, where most of the transport normally occurs. In
the second figure, the spectra for the ITG mode is supplemented with those of
the sub-dominant TE mode, except in the full scenario, where it is completely
stabilised, with $\gamma<0$ for the wave-numbers considered.

(a) eigenvalue spectra for different geometry models

(b) eigenvalue spectra in experimental geometry with added effects

Figure 3: (colour online) Growthrate spectra for $s-\alpha$, circular and
experimental magnetic equilibrium for an ITG mode dominated case with JET-like
parameters (3a), and for the experimental equilibrium with added degrees of
realism (3b). For the second figure, _expr._ represents full geometry case
with no added effects, in _expr.+C_ a background of $2\%C$ was added, and in
_full_ collisions are added. The last case corresponds to the “full-scenario”
NL runs in figure 7b. Parameters are presented in table. 1.

The effects of sheared toroidal rotation on the ITG growth rate in the
circular and experimental equilibrium are shown in figure 4. Both the
stabilising perpendicular velocity shear and the destabilising parallel
shearing rate are included. For the value of $q$ considered, the mode is
destabilised with increasing $\gamma_{\boldsymbol{E}\times\boldsymbol{B}}$.
This is because the stabilising component is proportional to $1/q$, wherefore
the destabilising effect dominates in the present case [48]. The shearing rate
in the considered JET experiment (marked in the figure) is, however, too small
to have a significant impact on the mode growth, and is left out of the
nonlinear treatment below.

Figure 4: (colour online) Scaling of ITG growthrate with
$\boldsymbol{E}\times\boldsymbol{B}$ shearing rate. The experimental shearing
rate has been indicated (vertical dashed line).

####  3.1.2 Nonlinear transport and fluctuation levels

By comparing the transport levels for the different geometries, the difference
in the turbulence can be assessed qualitatively. Time series and wave-number
spectra of heat and particle flux ($Q$ and $\Gamma$ respectively) for the main
ions were produced from nonlinear GENE simulations. The results are presented
in figure 5.

(a) timeseries of particle fluxes ($\Gamma$) with mean flux indicated

(b) timeseries of heat fluxes ($Q$) with mean flux indicated

(c) spectra of particle fluxes ($\Gamma$)

(d) spectra of heat fluxes ($Q$)

Figure 5: (colour online) Timeseries and spectra of particle ($\Gamma$) and
heat ($Q$) fluxes for nonlinear GENE simulations for $s-\alpha$ and circular
geometry, the experimental equilibrium, and the fully realistic case, with
parameters as in table. 1. Spectra show average over radial wavenumbers. The
timeseries and spectra have been normalised to the maximum of the
corresponding $s-\alpha$ entry. The spectra and time series show a decrease in
amplitude consistent with that expected from the linear eigenvalue spectra in
figure 3.

As can be seen in figures 5a and 5b, the turbulent fluctuation level is
increased when moving from $s-\alpha$ to realistic magnetic geometry. The same
trend holds true for the spectra in figures 5c and 5d. Simulations using
circular geometry [35, 37] place both the spectra and time series between the
$s-\alpha$ and experimental case without collisions. This is consistent with
the linear eigenvalues in figure 3, and with the results reported in [37].
With the addition of collisions and background carbon, the fluctuation levels
are brought down, due to a general reduction of the turbulence, as seen in the
the eigenvalue spectra figure 3b. The particle flux ($\Gamma$) is more
sensitive to collisions than the heat flux ($Q$), and therefore retains a
higher fluctuation level, despite the lower growthrates.

In the case with full realism, including collisions and carbon background, the
fluctuation levels are lowered to an intermediate level, consistent with
higher dissipation from collisions. This is also consistent with the linear
eigenvalues presented in figure 3, in which the growthrates for the “full
scenario” have an intermediate value between the circular (or $s-\alpha$) and
the case where only effects of realistic geometry is considered. As can be
seen in figures 5a and 5b, the full case was simulated for a shorter time
span, due to its larger computational cost, mainly due to the inclusion of
collisions, but also in part due to the $2\%$C background; see Section 2.1.

We note that the global confinement time in both L-mode and H-mode plasmas are
scaling favourably with elongation, opposite to the trend found here for the
core transport. This is likely a result of edge physics not included in the
present study [49].

In addition, the zonal flow activity was investigated for each case. However,
no significant differences were observed between the considered cases. The NL
results can thus be explained qualitatively by the linear physics.

###  3.2 Impurity peaking factors

Simulations of impurity transport using a _JET_ -like experimental magnetic
equilibrium have been compared to simulations with $s-\alpha$ and circular
geometry, using the gyrokinetic code GENE.

First, the effects sheared toroidal rotation in the different geometries is
investigated The effects on the impurity peaking factors are displayed in
figure 6. The peaking factors for several impurity species are shown as a
function of $\gamma_{\boldsymbol{E}\times\boldsymbol{B}}$ in the experimental
equilibrium. As observed, the peaking factors are reduced by the sheared
rotation. This leads to a reversal of the impurity pinch for
$\gamma_{\boldsymbol{E}\times\boldsymbol{B}}\gtrsim 0.23$ for all but the
lightest impurities ($Z\lesssim 4$) in the experimental equilibrium. In the
circular geometry the trends are the same, but the flux reversal for Tungsten
($Z=74$) is shifted to higher rotational sharing rate
($\gamma_{\boldsymbol{E}\times\boldsymbol{B}}\gtrsim 0.28$). This effect, due
to roto-diffusion, has been found to be a critical ingredient to include in
order to reproduce the Boron profile in ASDEX U [50]. For the considered
discharge, however, the shearing rate is small and hence it is neglected in
the non-linear analysis.

Figure 6: (colour online) Scaling of PF with
$\boldsymbol{E}\times\boldsymbol{B}$ shearing rate for different impurity
species in circular and shaped equilibrium for wavenumbers near the peak of
the corresponding growthrate spectra ($k_{\theta}\rho_{s}=0.3$ and
$k_{\theta}\rho_{s}=0.4$ respectively). The experimental shearing rate has
been indicated (vertical dashed line).

Next, the effects of the geometry on the scaling with impurity charge are
investigated. In figure 7a the scaling of the impurity peaking factor ($PF$)
with impurity charge ($Z$) is shown, for $s-\alpha$, circular and experimental
geometry. Non-linear (NL) gyrokinetic results are compared in figure 7b for
added degrees of realism. Error-bars in figure 7b correspond to a conservative
standard error estimate of $\pm\sigma$. This was calculated from the NL flux
data, where an effective sample size was gauged for the time series and used
to estimate the error for the mean flux. This estimate was then propagated
through to an error estimate for the peaking factor. The effects on the
$PF$-scaling of adding collisions and $2\%$C background, consistent with the
considered JET discharge, are shown in figure 7c.

For moderate to high impurity charge, weaker scalings of $PF$ with $Z$ were
consistently observed when departing from the $s-\alpha$ equilibrium, and the
level at which $PF$ saturates for high $Z$ was reduced in the experimental
geometry case. The lower $PF$ levels compared to $s-\alpha$ geometry can be
attributed to a reduction of the convective pinch due to effects of shaping,
but also due to effects of finite inverse aspect ratio ($\varepsilon$)
effects, as discussed in Section 2.2. To separate the effects of shaping from
the inconsistent $\varepsilon$ effects in the $s-\alpha$ model, simulations
using GENE’s circular geometry [37] were performed. Since $\alpha\ll 1$ in the
studied discharge, the difference between the results using circular and
experimental equilibrium in figure 7a can attributed to shaping effects,
mainly due to elongation. Further, a large increase in $PF$ was observed for
low $Z$ impurities, most notably He impurities, in circular and experimental
geometry. This is due to a reduction of the outward thermopinch in the more
realistic equilibrium models, as discussed in Section 3.3.

In figure 7b non-linear results using $s-\alpha$ geometry from [21] are
compared with NL gyrokinetic simulations using the experimental equilibrium
for two different sets of parameters. The first NL set of runs only uses the
experimental geometry, whereas the full case also includes a carbon background
and collisions, but neglects effects of sheared toroidal rotation (_full_ in
figure 3). With the introduction of these two effects, $PF$ is lowered for low
$Z$ and increased for high $Z$. As seen in the quasilinear runs, the NL GK
results are lower than those predicted by the $s-\alpha$ model, though a
qualitative agreement is reached in the full scenario, particularly for high
$Z$. When comparing QL and NL results, the former show a more dramatic scaling
than the latter and the QL results tend to over-estimate $PF$ for high $Z$, as
can be seen in figure 7c. In all models there is, as expected, a saturation of
$PF$ in the high $Z$ limit [16, 21, 22], and the observed NL and QL impurity
pinches qualitatively agree with the results in [21, 22].

Figure 7c shows the effect of adding collisions and a $2\%$ carbon background
to the $PF$ scaling with $Z$, using QL gyrokinetics. The “base case”
corresponds to the experimental MHD equilibrium, without any added effects
(_expr._ in figure 3). In all cases $k_{\theta}\rho_{\text{s}}=0.4$ was used,
since that is the approximate value for the maximum linear growth rate, as
shown in figure 3. As with the eigenvalue spectra (figure 3b), it is clear
that the addition of collisions to the model has a larger impact than the
addition of the $2\%$C background. The addition of the carbon background
slightly raises $PF$ for low $Z$, but the addition of collisions lowers it
considerably, as well as raising $PF$ for high charge numbers, as seen for the
NL runs. These results are consistent with previous observations of the
importance of impurity–main ion collisions for core impurity transport [51,
52]. We note that, although a reduction is seen in the peaking for carbon
impurities in the realistic case, the flat or hollow profile seen in
experiments is not reproduced. Though a further reduction is expected from
roto-diffusion, nonlinear simulations of the full scenario with sheared
toroidal rotation still show an inward nett transport of the background
carbon.

(a) comparison of $s-\alpha$, circular and experimental geometries for QL GENE

(b) comparison of NL GENE base case and full case (_expr._ and  _full_ in
figure 3) with results for $s-\alpha$ geometry from [21]; error-bars indicate
a standard error of $\pm\sigma$.

(c) effects of added realism on QL $Z$-scaling for $k_{\theta}\rho_{s}=0.4$

Figure 7: (colour online) Scalings of impurity peaking factor ($PF$) with
impurity charge number ($Z$) from JET discharge #67730 for $s-\alpha$,
circular and experimental equilibrium. Parameters as in table. 1.

###  3.3 Contributions to the impurity pinch

In order to gain more insight into the results presented in the last section,
the contributions to the peaking factor were calculated using the method
outlined in [46].

Figure 8 shows the effect of shaping on the thermodiffusive and convective
contributions to the impurity peaking factor (equation (4)). It is seen, that
the increase in peaking factor observed for lighter impurities is mainly due
to the thermopinch, where the inward second order contribution to the thermal
pinch comes into dominance over the outward term (equation (3)). This effect
becomes more pronounced for lower wave-numbers – especially in the shaped
equilibrium – due to the lower mode frequency (figure 3). For higher $Z$,
however, the peaking is determined by the balance of the inward convective
pinch and the diffusion. This contribution is only weakly dependent on the
impurity charge and wave-number, and is reduced substantially in the shaped
equilibrium.

In figure 9 diffusivities and pinches obtained from NL GENE are compared with
with data from [21]. From these results we note that, in contrast to the
diffusivity ($D_{Z}$), the convective velocity ($|V_{Z}|$) is lower in the
realistic equilibrium as compared to the $s-\alpha$ case, despite the
corresponding increase in fluctuation levels (figure 5b). This confirms the
conclusion that the main reason for the reduced peaking factors obtained in
the experimental equilibrium is a reduction of the convective velocity. When
comparing the collisionless shaped equilibrium with the full scenario, we note
that the diffusivity is decreased, as expected from the observed decrease in
fluctuation levels with the addition of collisions. For high $Z$, the
reduction in convective velocity is similar in both cases, compared to the
$s-\alpha$ equilibrium, leading to higher peaking factors in the full
scenario.

When comparing the scalings for low ($k_{\theta}\rho_{s}=0.2$) and high
($k_{\theta}\rho_{s}=0.5$) wavenumbers in the case with shaping, it is seen
that the thermal diffusion is dominated completely by the higher order terms
($\sim Z^{-2}$) in the low wavenumber case, turning from mostly outward to
inward. In the full case, the NL spectrum was seen to shift toward higher
wavenumbers, compared to the case with shaping but without collisions and
carbon background (figure 5). The difference between the non-linear
gyrokinetic results in figure 7b is therefore consistent with the
interpretation that higher order terms dominate the thermal diffusion for
lower wave numbers, thereby determining the shape of the peaking factor
scaling.

(a) contributions to the pinch for circular equilibrium

(b) contributions to the pinch for experimental equilibrium

Figure 8: (colour online) Contributions to the peaking factor
($PF_{k_{\theta}\rho_{s}}^{\text{tot}}$) from thermodiffusion
($PF_{k_{\theta}\rho_{s}}^{\nabla T}$) and pure convection
($PF_{k_{\theta}\rho_{s}}^{\text{conv}}$) as a function of impurity charge
($Z$). Results from QL GENE for wavenumbers as in figure 7, with and without
effects of shaping ($\kappa=1.0$, $\kappa=1.37$).  Figure 9: (colour online)
Scaling of the impurity diffusivity ($D_{Z}$) and total pinch ($RV_{Z}$) with
impurity charge for $s-\alpha$ geometry and experimental geometry, with and
without collisions and $2\%$C background. Data from NL GENE simulations with
standard errors $~{}10$–$20\%$ (omitted for clarity).

##  4 Conclusions

The effects of the choice of equilibrium model on impurity transport in ITG/TE
mode driven turbulence were studied using gyrokinetic simulations. Results
were obtained and contrasted for the experimental MHD equilibrium and the
simpler $s-\alpha$ and concentric circular geometries. These results were
extended by adding degrees of realism in the physics description, such as a 2%
carbon background, collisions, and sheared toridal rotation. The gyrokinetic
results, obtained with the GENE code in both quasi- and non-linear mode, were
compared with results from a previous study, as well as with expectations from
a dedicated impurity injection L-mode discharge at JET. It is found that the
different equilibria give qualitatively similar results, but with significant
quantitative differences.

Linearly, a destabilisation and shift of the growthrate spectrum to higher
wave numbers was seen when departing from the simpler geometries, mainly due
to elongation. This resulted in larger heat and particle transport, as seen in
the nonlinear simulations. The addition of collisions was seen to stabilise
the spectrum and reduce transport.

The effect of sheared rotation on the mode stability and the impurity
transport was studied. A weak destabilisation of the ITG mode with
$\boldsymbol{E}\times\boldsymbol{B}$ shear was found. However, this was seen
to be a minor contribution for experimentally feasible values of the
rotational shearing rate for the considered discharge.

The impurity peaking factors ($PF$s), computed by finding local density
gradients corresponding to zero particle flux, were derived. The peaking
factor was observed to saturate at levels far below neo-classical predictions
for high impurity charge numbers. The level of the saturation was considerably
lowered when the experimental equilibrium was used, typically by a factor of
$\sim 2$. Comparing the diffusivity and convective pinch velocity in the
different geometries, it was found that the reduction of $PF$ was mainly due
to a reduction of the convective pinch.

By decomposing the contributions to the peaking factor into thermodiffusion
and pure convection, it was seen that second order contributions to the
thermal diffusion became important for low impurity charge. Whereas the first
order term is proportional to $\omega_{r}$, and therefore outward for ion
temperature gradient modes, the second order contribution to the pinch is
always inward, and therefore leads to higher $PF$s. This explains the more
pronounced increase in peaking observed for low charge numbers in the shaped
equilibria, where the real frequency is lower.

Collisions were seen to affect the impurity peaking factors considerably,
lowering $PF$ for light impurities, while increasing it for heavy impurities.
It was seen, that the addition of collisions and background carbon shifted the
nonlinear spectra towards higher wavenumbers, assosciated with higher real
frequencies. Therefore, the the difference between the nonlinear results in
the experimental equilibrium could be also understood through the second order
term in the thermopinch, which was concluded to dominate this pinch
contribution for the non-collisional case.

However, even combining the effect of collisions with the effect of geometry
and rotation in the most realistic simulations, the observed peaking factor
for carbon was still too large to explain the rather flat C profile observed
in the experiment. One explanation for this descrepancy could lie in the
sensitivity of the thermopinch to main ion and impurity temperature, and to
the gradients of thereof. A moderate variation in these may bring the
background carbon peaking factor down to zero, however, such a treatment is
left for future studies. We also note that the treatment of rotation employed
in this study neglects centrifugal and Coriolis forces, which may lower the
impurity peaking further.

## Acknowledgements

The main simulations were performed on resources provided on the
Lindgren333See http://www.pdc.kth.se/resources/computers/lindgren/ for details
on Lindgren high performance computer, by the Swedish National Infrastructure
for Computing (SNIC) at Paralleldatorcentrum (PDC). Additional computations
were carried out on resources at Chalmers Centre for Computational Science and
Engineering (C3SE)444See http://www.c3se.chalmers.se, also provided by SNIC.

The authors would like to thank F Jenko, T Görler, F Merz, MJ Püschel, D Told,
and the rest of the GENE team at IPP–Garching for their valuable support and
input.

## References

  * [1] Fröjdh M, Liljeström M and Nordman H 1992 Nucl. Fusion 32 419
  * [2] Basu R, Jessen T, Naulin V and Rasmussen J J 2003 Phys. Plasmas 10 2696
  * [3] Estrada-Mila C, Candy J and Waltz R W 2005 Phys. Plasmas 12 022305
  * [4] Naulin V 2005 Phys. Rev. E 71 015402
  * [5] Priego M, Garcia O E, Naulin V and Rasmussen J J 2005 Phys. Plasmas 12 062312
  * [6] Fülöp T and Weiland J 2006 Phys. Plasmas 13 112504
  * [7] Bourdelle C, Garbet X, Imbeaux F, Casati A, Dubuit N, Guirlet R and Parisot T 2007 Phys. Plasmas 14 112501
  * [8] Dubuit N, Garbet X, Parisot T, Guirlet R and Bourdelle C 2007 Phys. Plasmas 14 042301
  * [9] Camenen Y, Peeters A G, Angioni C, Casson F J, Hornsby W A, Snodin A P and Strintzi D 2009 Phys. Plasmas 16 012503
  * [10] Fülöp T, Braun S and Pusztai I 2010 Phys. Plasmas 17 062501
  * [11] Futatani S, Garbet X, Benkadda S and Dubuit N 2010 Phys. Rev. Lett. 104 015003
  * [12] Hein T and Angioni C 2010 Phys. Plasmas 17 012307
  * [13] Moradi S, Tokar M Z and Weyssow B 2010 Phys. Plasmas 17 012101
  * [14] Fülöp T and Moradi S 2011 Phys. Plasmas 18 030703
  * [15] Nordman H, Singh R, Fülöp T, Eriksson L G, Dumont R, Anderson J, Kaw P, Strand P, Tokar M and Weiland J 2008 Phys. Plasmas 15 042316
  * [16] Angioni C and Peeters A G 2006 Phys. Rev. Lett. 96 095003
  * [17] Angioni C, Carraro L, Dannert T, Dubuit N, Dux R, Fuchs C, Garbet X, Garzotti L, Giroud C, Guirlet R, Jenko F, Kardaun O J W F, Lauro-Taroni L, Mantica P, Maslov M, Naulin V, Neu R, Peeters A G, Pereverzev G, Puiatti M E, Pütterich T, Stober J, Valovič M, Valisa M, Weisen H, Zablotsky A, ASDEX Upgrade team and JET–EFDA contributors 2007 Phys. Plasmas 14 055905
  * [18] Nordman H, Fülöp T, Candy J, Strand P and Weiland J 2007 Phys. Plasmas 14 052303
  * [19] Angioni C, Peeters A G, Pereverzev G V, Bottino A, Candy J, Dux R, Fable E, Hein T and Waltz R E 2009 Nucl. Fusion 49 055013
  * [20] Fülöp T and Nordman H 2009 Phys. Plasmas 16 032306
  * [21] Nordman H, Skyman A, Strand P, Giroud C, Jenko F, Merz F, Naulin V, Tala T and the JET–EFDA contributors 2011 Plasma Phys. Contr. F. 53 105005
  * [22] Skyman A, Nordman H and Strand P 2012 Phys. Plasmas 19 032313
  * [23] Skyman A, Nordman H and Strand P 2012 Nucl. Fusion 52 114015
  * [24] Mollén A, Pusztai I, Fülöp T, Kazakov Y O and Moradi S 2012 Phys. Plasmas 19 052307
  * [25] Moradi S, Pusztai I, Mollén A and Fülöp T 2012 Phys. Plasmas 19 032301
  * [26] Kazakov Y O, Fülöp T and Van Eester D 2013 Nucl. Fusion 53 053014
  * [27] Mollén A, Pusztai I, Fülöp T and Moradi S 2013 Phys. Plasmas 20 032310
  * [28] Dux R, Neu R, Peeters A G, Pereverzev G, Mück A, Ryter F and Stober J 2003 Plasma Phys. Contr. F. 45 1815
  * [29] Puiatti M E, Valisa M, Mattioli M, Bolzonella T, Bortolon A, Coffey I, Dux R, von Hellermann M, Monier-Garbet P, Nave M F F, Ongena J and contributors to the EFDA–JET workprogramme 2003 Plasma Phys. Contr. F. 45 2011
  * [30] Puiatti M E, Valisa M, Angioni C, Garzotti L, Mattioli M, Carraro L, Coffey I and Sozzi C 2006 Phys. Plasmas 13 042501
  * [31] Giroud C, Barnsley R, Buratti P, Coffey I H, von Hellermann M, Jupén C, Lawson K D, Meigs A, O’Mullane M, Whiteford A D, Zastrow K D and the JET–EFDA contributors 2007 Nucl. Fusion 47 313
  * [32] Giroud C et al 2009 Impurity characterisation in JET operational scenarios 12th International Workshop on “H-mode Physics and Transport Barriers” Princeton, USA URL http://iopscience.iop.org/0029-5515/50/6
  * [33] Matthews G F, Edwards P, Greuner H, Loving A, Maier H, Martens P, Philipps V, Riccardo V, Rubel M, Ruset C, Scmidt A and Villedieu E 2009 Phys. Scr. 2009 014030
  * [34] Doyle E J et al 2007 Plasma confinement and transport Progress in the ITER Physics Basis vol 47 (Nucl. Fusion) chap 2, p S18
  * [35] Xanthopolous P and Jenko F 2006 Phys. Plasmas 13 092301
  * [36] Belli E A, Hammett G W and Dorland W 2008 Phys. Plasmas 15 092303
  * [37] Lapillonne X, Brunner S, Dannert T, Jolliet S, Marioni A, Villard L, Görler T, Jenko F and Merz F 2009 Phys. Plasmas 16 032308
  * [38] Burckel A, Sauter O, Angioni C, Candy J, Fable E and Lapillonne X 2010 J. Phys.: Conf. Ser. 260 012006
  * [39] Told D and Jenko F 2010 Phys. Plasmas 17 042302
  * [40] Jenko F, Dorland W, Kotschenreuther M and Rogers B N 2000 Phys. Plasmas 7 1904
  * [41] Merz F 2008 Gyrokinetic Simulation of Multimode Plasma Turbulence Ph.d. thesis (monography) Westfälischen Wilhelms-Universität Münster
  * [42] Arter W 1995 Rep. Prog. Phys. 58 1
  * [43] Bakunin O G 2003 Plasma Phys. Contr. F. 45 1909
  * [44] Nilsson J, Liljeström M and Weiland J 1990 Phys. Fluids B 2 2568
  * [45] Weiland J 2000 Collective Modes in Inhomogeneous Plasmas (Bristol, UK: IoP Publishing)
  * [46] Casson F J, Peeters A G, Angioni C, Camenen Y, Hornsby W A, Snodin A P and Szepesi G 2010 Phys. Plasmas 17 102305
  * [47] Anderson J, Nordman H and Weiland J 2000 Plasma Phys. Contr. F. 42 545
  * [48] Barnes M, Parra F I, Highcock E G, Schekochihin A A, Cowley S C and Roach C M 2011 Phys. Rev. Lett. 106 175004
  * [49] ITER Physics Expert Group on Confinement and Transport, ITER Physics Expert Group on Confinement Modelling and Database and ITER Physics Basis Editors 1999 Plasma confinement and transport ITER Physics Basis vol 39 (Nucl. Fusion) chap 2, p 2175
  * [50] Angioni C, McDermott R, Fable E, Fischer R, Pütterich T, Ryter F, Tardini G and the ASDEX Upgrade Team 2011 Nucl. Fusion 51 023006
  * [51] Angioni C, Peeters A G, Jenko F and Dannert T 2005 Phys. Plasmas 12 112310
  * [52] Moradi S, Tokar M Z, Singh R and Weyssow B 2009 Nucl. Fusion 49 085007