# Uncertainty-aware Consistency Learning for
Cold-Start Item Recommendation

Taichi Liu<EMAIL_ADDRESS>Tsinghua University ,  Chen Gao
<EMAIL_ADDRESS>Tsinghua University Huawei Noah’s Ark Lab ,  Zhenyu Wang
<EMAIL_ADDRESS>Tsinghua University ,  Dong Li
<EMAIL_ADDRESS>Huawei Noah Ark’s Lab ,  Jianye Hao<EMAIL_ADDRESS>Huawei Noah’s Ark Lab ,  Depeng Jin<EMAIL_ADDRESS>Tsinghua University
and  Yong Li<EMAIL_ADDRESS>Tsinghua University

(2023)

###### Abstract.

Graph Neural Network (GNN)-based models have become the mainstream approach
for recommender systems. Despite the effectiveness, they are still suffering
from the cold-start problem, i.e., recommend for few-interaction items.
Existing GNN-based recommendation models to address the cold-start problem
mainly focus on utilizing auxiliary features of users and items, leaving the
user-item interactions under-utilized. However, embeddings distributions of
cold and warm items are still largely different, since cold items’ embeddings
are learned from lower-popularity interactions, while warm items’ embeddings
are from higher-popularity interactions. Thus, there is a seesaw phenomenon,
where the recommendation performance for the cold and warm items cannot be
improved simultaneously. To this end, we proposed a Uncertainty-aware
Consistency learning framework for Cold-start item recommendation (shorten as
UCC) solely based on user-item interactions. Under this framework, we train
the teacher model (generator) and student model (recommender) with consistency
learning, to ensure the cold items with additionally generated low-uncertainty
interactions can have similar distribution with the warm items. Therefore, the
proposed framework improves the recommendation of cold and warm items at the
same time, without hurting any one of them. Extensive experiments on benchmark
datasets demonstrate that our proposed method significantly outperforms state-
of-the-art methods on both warm and cold items, with an average performance
improvement of 27.6%.

Recommender System; Cold-Start Item; Graph Neural Networks

††journalyear: 2023††copyright: acmlicensed††conference: Proceedings of the
46th International ACM SIGIR Conference on Research and Development in
Information Retrieval; July 23–27, 2023; Taipei, Taiwan††booktitle:
Proceedings of the 46th International ACM SIGIR Conference on Research and
Development in Information Retrieval (SIGIR ’23), July 23–27, 2023, Taipei,
Taiwan††price: 15.00††doi: 10.1145/3539618.3592078††isbn:
978-1-4503-9408-6/23/07††ccs: Information systems Recommender systems

##  1\. Introduction

Recommender systems have become a fundamental service for filtering
information by learning from collected user behavioral data and inferring
users’ personalized demands. Currently, GNN-based models, with their strengths
in learning from structured data and capturing high-order similarity, have
become state-of-the-art approaches in various recommendation tasks (Gao et
al., 2023). GNNs have a strong ability to incorporate information from the
neighboring nodes, helping to understand the local context of nodes. However,
it leads to the poor performance for the cold-start recommendation (Volkovs et
al., 2017; Zhang et al., 2019; Hu et al., 2019; Zhao et al., 2022;
Abdollahpouri et al., 2017; Cañamares and Castells, 2018; Yin et al., 2012),
which occurs when a new item node is added to the graph. Specifically, there
are not enough neighbors for the new item due to few user-item interactions.
The cold-start recommendation problem can be particularly challenging in
scenarios where the user-item graph rapidly evolves, such as short-video
recommendation (Ling et al., 2022).

To address this challenge, existing models (Volkovs et al., 2017; Zhu et al.,
2021; Du et al., 2022; Chen et al., 2022; Togashi et al., 2021) have proposed
various techniques for incorporating additional information into GNNs, such as
social networks, node attributes, and other metadata, which can partly improve
the performance for cold-start items. However, cold items’ embeddings are
still continuously learned from fewer interactions, while warm items’
embeddings are still mainly learned from adequate interactions (Abdollahpouri
et al., 2019, 2021), which leads to the issue that the distributions of
embeddings of cold and warm items are significantly different. That is, the
collected interaction of cold items cannot reflect the real item
characteristics due to less exposure, thus resulting in the seesaw phenomenon,
i.e., improves the recommendation of either the cold or warm items will hurt
the other (Chen et al., 2020a; Zhang et al., 2021).

To this end, we proposed a framework for bridging the distribution gap between
cold items and warm items in the cold-start recommendation, with two key
designs included: a) Uncertainty-aware Interaction Generation, b) Graph-based
Teacher-Student Consistency Learning. In particular, we first introduce low-
uncertainty interactions and take use of generated interactions for embedding
learning, to address the distribution difference which causes the seesaw
phenomenon. We then develop a teacher-student training paradigm that improves
the robustness (Park et al., 2006; Wei et al., 2021) in embedding learning for
cold and warm items. Specifically, we propose a contrastive loss to keep the
consistency of item embeddings before and after the interaction generation. We
further fine-tune the student model on the recommendation task, maintaining
the consistency between the teacher model and student model. Our contributions
can be summarized as follows:

  * •

We propose approach the cold-start problem from the novel perspective of
addressing the distribution difference between cold and warm items, by
proposing a general UCC framework.

  * •

Under the framework, we first address the distribution gap between cold and
warm items through our uncertainty-aware interaction generation, and keep both
item-level and model-level similarity with consistency learning, which can
improve model’s robustness.

  * •

We conduct extensive experiments to demonstrate the effectiveness of our
framework. The results shown an average improvement of 27.6% in ranking
metrics compared to the baseline method on two widely-used benchmark datasets.

Figure 1. The overall framework of UCC. A teacher model (the upper blue frame)
is first trained under weak augmentation regularization. Then interactions are
generated with uncertainty-aware selection process. Finally the generated
interactions serve for the student model training (the lower green frame).
Strong augmentation based consistency regularization and integrated teacher
representation increase the robustness of the student model.

##  2\. Problem Formulation

Notations. The common GNN-based model paradigm can be formulated as follows:
Let $\mathcal{U}=\\{u\\}$, $\mathcal{I}=\\{i\\}$ denote the set of users and
items respectively. The interactions of user and item can be regard as
$\mathcal{O}=\left\\{(u,i^{+})\mid
u\in\mathcal{U},i^{+}\in\mathcal{I}\right\\}$, where each pair represents each
observed feedback.

Input: Most GNN-based models treat user-item interactions as a bipartite
graph, where the node set is $\mathcal{V}=\mathcal{U}\cup\mathcal{I}$ and the
edge set is $\mathcal{G}=(\mathcal{V},\mathcal{O}^{+})$. In the training
process of graph representation learning, we use
$\textbf{E}_{u}=[e_{u_{1}},\cdots,e_{u_{M}}]\in\mathbb{R}^{M\times
D},\textbf{E}_{i}=[e_{i_{1}},\cdots,e_{i_{N}}]\in\mathbb{R}^{N\times D}$ to
denote the user embedding and item embedding, where D is the dimension size of
embedding, $M,N$ are the number of users and items respectively.

Output: Recommender models estimate the relations of unobserved user-item
pairs through the dot products of their embeddings (Grover and Leskovec, 2016;
He et al., 2017; Hsieh et al., 2017). Specifically, the scores $\\{s_{mn}\\}$
of a given user m and a given item n, where $s_{mn}$ is calculated as:

(1) $s_{mn}=e_{u_{m}}e_{i_{n}}^{T},$

then rank items by their scores - a high score means that the user prefers.
The top-k items from the ranking list are adopted into a candidate set and
recommended to the user.

##  3\. Methodology

###  3.1. Uncertainty-aware Interaction Generation

For the cold-start recommendation, one main reason for the limited performance
is the different distribution between cold items and warm items caused by
popularity bias. Hence, we propose to generate similar distributed
interactions for cold items as warm items, and we need to select accurate and
reliable interactions from unobserved interactions.

However, since totally unbiased and accurate interactions for the original
data used are unavailable, especially for cold items, it is hard to
discriminate whether an interaction is reliable. We thus propose to estimate
the uncertainty degree of each user-item interaction, and select low-
uncertainty interactions for each item. We estimate the uncertainty of the
particular interaction between the user $u_{m}$ and item $i_{n}$ by their
similarity. We apply the cosine distance to calculate the similarity:

(2)
$d_{mn}=\frac{\left|e_{u_{m}}e_{i_{n}}^{T}\right|}{\left\|e_{u_{m}}\right\|\left\|e_{i_{n}}\right\|},$

With the pre-trained recommender, we have the ranking scores
$\\{s_{mn}\\}_{n=1}^{N}$ of each item $i_{n}$ for all users. The average of
all ranking scores measures the overall uncertainty degree of the item $i_{n}$
as follows:

(3) $\bar{s}_{n}=\frac{1}{M}\sum\limits_{m=1}^{M}s_{mn},$

where $\bar{s}_{n}$ provides a description of the overall interaction
uncertainty of the item $i_{n}$. We aim to force the cold items with generated
interactions to have a similar distribution as the warm items to improve the
recommendation performance for both cold and warm items. A small $d_{nm}$
interactions are regarded as uncertain ones and are filtered out in the
generation phase. This selection strategy based on the interaction uncertainty
guides the generator on more certain interactions to mitigate the popularity
bias in the cold-start problem. The selection is as follows:

(4) $\widehat{\mathcal{O}}_{n}=I(d_{mn}>\alpha\bar{s}_{n}),$

where $\alpha$ is a pre-defined parameter and $I$ is the indicator function.

###  3.2. Teacher-Student Consistency Learning

The consistency learning for our framework consists of two parts:

####  3.2.1. Item-level consistency learning.

We adopt the contrastive loss for the item embeddings between pre-generation
and post-generation. Inspired by (Sohn et al., 2020), we adopt two
augmentations - strong augmentation and weak augmentation. Specifically, the
weak augmentation drops out the edges in the graph with a dropout ratio
$\rho$, which can be formulated as follows:

(5) $\mathcal{G}^{w}=(\mathcal{V},\mathcal{M}\cdot\mathcal{O}^{+}),$

where $\mathcal{M}\in\\{0,1\\}^{|\mathcal{O}^{+}|}$ is the masking vector. The
strong augmentation based on generated labels adds more edges to the graph,
which can be formulated as follows,

(6) $\mathcal{G}^{s}=(\mathcal{V},\mathcal{O}^{+}+\widehat{\mathcal{O}}),$

Two augmentation operations on the graph generates two different views of each
node. For the item node $i$, we denote its two different views as
$z_{i}^{\prime}$, $z_{i}^{\prime\prime}$ (on weak and strong graphs
respectively). We perform a consistency regularization for these different
views, to encourage the similarity between different views of the same node.
To implement the consistency regularization, the contrastive loss (Chen et
al., 2020b) is as follows,

(7) $\displaystyle\mathcal{L}_{\textbf{cr,
item}}=\sum_{i\in\mathcal{I}}-\log\frac{\exp({\rm
sim}(z_{i}^{\prime},z_{i}^{\prime\prime})/\tau)}{\sum_{j\in\mathcal{I}}\exp({\rm
sim}(z_{i}^{\prime},z_{j}^{\prime\prime})/\tau)},$

where $\mathcal{L}_{\textbf{cr, item}}$ denotes the item-side consistency
regularization loss for the teacher and the student model, respectively. Here
${\rm sim}(\cdot,\cdot)$ denotes the similarity function, for which we use
cosine similarity, and $\tau$ is a pre-defined hyper-parameter. The user-side
consistency regularization loss $\mathcal{L}_{\textbf{cr, user}}$ can be
similarly calculated, and we omit it due to the limit of space. The final
consistency loss $\mathcal{L}_{\textbf{cr}}=\mathcal{L}_{\textbf{cr,
item}}+\mathcal{L}_{\textbf{cr, user}}$ serves for consistency regularization.
For the recommendation loss we use Bayesian personalized ranking (BPR) (Rendle
et al., 2012) loss to optimize the representation learning, which can be
formulated as follows,

(8) $\mathcal{L}_{\textbf{rec}}=\sum_{(u,i^{+},i^{-})\in
O}-\ln\sigma\left(\hat{y}_{ui^{+}}-\hat{y}_{ui^{-}}\right)+\lambda\|\Theta\|_{2}^{2},$

where $(u,i^{+})$, $(u,i^{-})$ represent observed/unobserved interaction pairs
and $\|\Theta\|_{2}^{2}$ is L2-regularization of model’s parameters. Then our
final loss of the whole framework is:

(9)
$\mathcal{L}_{\textbf{total}}=\mathcal{L}_{\textbf{rec}}+\mu\mathcal{L}_{\textbf{cr}},$

where $\mu$ is a hyper-parameter.

####  3.2.2. Model-level consistency learning.

To maintain the consistency between the teacher model and student model, We
denote the embedding of the student model’s embedding as $\textbf{E}^{s}$.
After each iteration, we propose to accumulate the teacher embedding into the
student embedding:

(10) $\textbf{E}^{s}\leftarrow\gamma\textbf{E}^{s}+(1-\gamma)\textbf{E}^{t},$

where $\gamma$ denotes a momentum term that controls the impact of the teacher
model. Here the accumulation step includes both user embeddings and item
embeddings.

##  4\. Experiments

###  4.1. Experimental Setting

Datasets. We conduct experiments on benchmark datasets: Yelp and Amazon-Book,
following the same 10-core setting as (He et al., 2017, 2020). We split all
user-item interactions into training, validation, and testing set with the
ratio of 7:1:2, evaluating the top-$K$ recommendation performance with two
widely-used metrics Recall@K and NDCG@K where $K=20$, following (He et al.,
2020).

Baselines. We select LightGCN(He et al., 2020) as the backbone of our GNN-
based model. Since we mainly focus on modeling user-item interactions but not
the features of items, we compare our model with state-of-the-art
recommendation models which can be classified into two classes: generative
model NB(Saito, 2020) and denosing model IRBPR(Wang et al., 2021b), ADT(Wang
et al., 2021a), SGL(Wu et al., 2021) for user-item interactions.

Table 1. Performance comparison on benchmark datasets. Dataset | Yelp | Amazon-Book
---|---|---
Metric | Recall | NDCG | Recall | NDCG
LightGCN | 0.0639 | 0.0515 | 0.0410 | 0.0318
IRBPR | 0.0538 | 0.0438 | 0.0325 | 0.0246
NB | 0.0640 | 0.0526 | 0.0425 | 0.0328
ADT | 0.0691 | 0.0388 | 0.0335 | 0.0176
SGL | 0.0673 | 0.0554 | 0.0477 | 0.0378
UCC | 0.0713 | 0.0587 | 0.0523 | 0.0410
Improv. | +11.6% | +14.0% | +27.6% | +28.9%

(a) Yelp

(b) Amazon-book

Figure 2. Performance comparison over item groups.

###  4.2. Recommendation Performance Comparison

We evaluate the model from two aspects, overall recommendation performance and
cold-start recommendation performance.

####  4.2.1. Overall Performance.

The details of the comparison are listed in Table 1. The following points are
observed:

$\bullet$ Our UCC outperforms previous models significantly on all datasets.
On the Yelp dataset, our method outperforming recall of LightGCN by 11.6%.
Compared to the state-of-the-art graph based method SGL, our method also
obtains a higher recall of SGL by 6%. The improvement is consistent for
Amazon-Book. This demonstrates the effectiveness of our uncertainty-aware
consistency learning.

$\bullet$ Our method improves Recall@20 significantly more on the Amazon-Book
dataset - 27.6% compared to LightGCN and 9.6% higher than SGL. The reason is
that the interactions for cold items in amazon-book are more biased and
sparse. This validates that our method targets at the cold-start problem and
well alleviates the restriction of different distribution between cold items
and warm items.

####  4.2.2. Cold-Start Performance.

We further compare our methods with LightGCN on cold and warm items separately
to illustrate its effectiveness for cold-start recommendation. We split items
into ten groups based on the popularity, meanwhile keeping the total number of
interactions of each group the same. We treat the last two groups as cold-
start items. The larger the GroupID is, the warmer the items are. We evaluate
our method on different groups and plot the Recall@20 in Figure 2. From the
results, we notice that our method outperforms LightGCN on all groups, no
matter whether the items are warm or cold. In comparison, many current methods
improve the performance of cold items at the cost of warm item accuracy. Also,
our method improves the recall of LightGCN more prominently on cold items. On
the Yelp dataset, our method obtains nearly 7 times compared to recall of
LightGCN. The most significant improvement comes from the groupid 1 items of
Amazon-Book. In this group, the recall of amazon-book obtains is 0.0009. In
contrast, the recall of our method obtains 0.0045, improved by 400%. This
strongly demonstrates the effectiveness of our method for cold-start
recommendation and the problem of the seesaw phenomenon.

###  4.3. Ablation Study of UCC

In this section, we perform ablation studies to evaluate each component of
UCC, which consists of two parts: Uncertainty-aware Interaction Generation and
Teacher-Student Consistency Learning.

####  4.3.1. Uncertainty-aware interaction generation

We calculate the average number of generated interactions for different item
groups. We find the number of generated interactions is adaptive for different
items. For cold items, the average number of produced interactions is about
3.5, which is about 1.0 for warm items in comparison. Since the generated low-
uncertainty interactions contribute more for cold items, the uncertainty
guided strategy well alleviating the different distribution between warm items
and cold items. We also notice that the simple usage of item-side generated
interactions helps improves the performance. The performance increase
indicates that our item-side adaptive generated interactions can indeed
alleviate the cold-start problem. User-side generated interactions, in
comparison, cannot improve the recommendation ability of the model - even
lowering the Recall@20 from 0.0673 to 0.0671. This is because user-side
generated interactions are biased and exacerbate the distribution between cold
and warm items. This demonstrates the superiority of our uncertainty-aware
interaction generation.

####  4.3.2. Teacher-Student Consistency Learning

Figure 3. Study of teacher-student learning.

We show the effectiveness of our teacher-student learning approach. We design
two comparing experiments: 1) Finetune learning: We train the student
representations with generated interactions directly from the teacher model
embeddings without accumulating teacher representations. 2) Continuous
learning: We train the student embeddings with accumulating with teacher
embeddings without using generated interactions. We plot the results in Figure
3. We observe that our teacher-student learning approach obtains a higher
performance. The superior comes from two aspects. Compared with finetune
learning, our student model maintains more information from the teacher model.
The information within the teacher embeddings alleviates the negative effect
of noise from generated interactions. Compared with continuous learning, we
provide more interactions for cold items, thus mitigating the cold-start
problem. Through the proposed teacher-student learning framework, similar
distributions of information and data are accumulated together, thus helping
obtain a better recommendation performance.

##  5\. Conclusion

In this paper, we identify that the main limitation comes from the
interaction-distribution difference between cold and warm items. To address
it, we propose an uncertainty-aware consistency learning framework. The
extraordinary performance on benchmark datasets and its easy-to-use property
make our framework practical in both academia and industry.

##  6\. ACKNOWLEDGEMENT

This work was supported by the National Key Research and Development Program
of China under 2022YFB3104702, the National Natural Science Foundation of
China under 62272262, 61972223, U1936217, and U20B2060, the Fellowship of
China Postdoctoral Science Foundation under 2021TQ0027 and 2022M710006.

## References

  * (1)
  * Abdollahpouri et al. (2017) Himan Abdollahpouri, Robin Burke, and Bamshad Mobasher. 2017\.  Controlling popularity bias in learning-to-rank recommendation. In _Proceedings of the eleventh ACM conference on recommender systems_. 42–46.
  * Abdollahpouri et al. (2019) Himan Abdollahpouri, Robin Burke, and Bamshad Mobasher. 2019\.  Managing popularity bias in recommender systems with personalized re-ranking.  _arXiv preprint arXiv:1901.07555_ (2019).
  * Abdollahpouri et al. (2021) Himan Abdollahpouri, Masoud Mansoury, Robin Burke, Bamshad Mobasher, and Edward Malthouse. 2021.  User-centered evaluation of popularity bias in recommender systems. In _Proceedings of the 29th ACM Conference on User Modeling, Adaptation and Personalization_. 119–129.
  * Cañamares and Castells (2018) Rocío Cañamares and Pablo Castells. 2018.  Should I follow the crowd? A probabilistic analysis of the effectiveness of popularity in recommender systems. In _The 41st International ACM SIGIR Conference on Research & Development in Information Retrieval_. 415–424.
  * Chen et al. (2022) Hao Chen, Zefan Wang, Feiran Huang, Xiao Huang, Yue Xu, Yishi Lin, Peng He, and Zhoujun Li. 2022.  Generative adversarial framework for cold-start item recommendation. In _Proceedings of the 45th International ACM SIGIR Conference on Research and Development in Information Retrieval_. 2565–2571.
  * Chen et al. (2020a) Jiawei Chen, Hande Dong, Xiang Wang, Fuli Feng, Meng Wang, and Xiangnan He. 2020a.  Bias and debias in recommender system: A survey and future directions.  _arXiv preprint arXiv:2010.03240_ (2020).
  * Chen et al. (2020b) Ting Chen, Simon Kornblith, Mohammad Norouzi, and Geoffrey Hinton. 2020b.  A simple framework for contrastive learning of visual representations. In _ICML_.
  * Du et al. (2022) Jing Du, Zesheng Ye, Lina Yao, Bin Guo, and Zhiwen Yu. 2022.  Socially-aware Dual Contrastive Learning for Cold-Start Recommendation. In _Proceedings of the 45th International ACM SIGIR Conference on Research and Development in Information Retrieval_. 1927–1932.
  * Gao et al. (2023) Chen Gao, Yu Zheng, Nian Li, Yinfeng Li, Yingrong Qin, Jinghua Piao, Yuhan Quan, Jianxin Chang, Depeng Jin, Xiangnan He, et al. 2023\.  A survey of graph neural networks for recommender systems: challenges, methods, and directions.  _ACM Transactions on Recommender Systems_ 1, 1 (2023), 1–51.
  * Grover and Leskovec (2016) Aditya Grover and Jure Leskovec. 2016.  node2vec: Scalable feature learning for networks. In _Proceedings of the 22nd ACM SIGKDD international conference on Knowledge discovery and data mining_. 855–864.
  * He et al. (2020) Xiangnan He, Kuan Deng, Xiang Wang, Yan Li, Yongdong Zhang, and Meng Wang. 2020\.  Lightgcn: Simplifying and powering graph convolution network for recommendation. In _SIGIR_.
  * He et al. (2017) Xiangnan He, Lizi Liao, Hanwang Zhang, Liqiang Nie, Xia Hu, and Tat-Seng Chua. 2017\.  Neural collaborative filtering. In _WWW_.
  * Hsieh et al. (2017) Cheng-Kang Hsieh, Longqi Yang, Yin Cui, Tsung-Yi Lin, Serge Belongie, and Deborah Estrin. 2017\.  Collaborative metric learning. In _Proceedings of the 26th international conference on world wide web_. 193–201.
  * Hu et al. (2019) Liang Hu, Songlei Jian, Longbing Cao, Zhiping Gu, Qingkui Chen, and Artak Amirbekyan. 2019\.  Hers: Modeling influential contexts with heterogeneous relations for sparse and cold-start recommendation. In _Proceedings of the AAAI Conference on Artificial Intelligence_ , Vol. 33. 3830–3837.
  * Ling et al. (2022) Chen Ling, Jeremy Blackburn, Emiliano De Cristofaro, and Gianluca Stringhini. 2022.  Slapping Cats, Bopping Heads, and Oreo Shakes: Understanding Indicators of Virality in TikTok Short Videos. In _14th ACM Web Science Conference 2022_. 164–173.
  * Park et al. (2006) Seung-Taek Park, David Pennock, Omid Madani, Nathan Good, and Dennis DeCoste. 2006.  Naïve filterbots for robust cold-start recommendations. In _Proceedings of the 12th ACM SIGKDD international conference on Knowledge discovery and data mining_. 699–705.
  * Rendle et al. (2012) Steffen Rendle, Christoph Freudenthaler, Zeno Gantner, and Lars Schmidt-Thieme. 2012.  BPR: Bayesian personalized ranking from implicit feedback.  _arXiv preprint arXiv:1205.2618_ (2012).
  * Saito (2020) Yuta Saito. 2020\.  Asymmetric tri-training for debiasing missing-not-at-random explicit feedback. In _SIGIR_.
  * Sohn et al. (2020) Kihyuk Sohn, David Berthelot, Nicholas Carlini, Zizhao Zhang, Han Zhang, Colin A Raffel, Ekin Dogus Cubuk, Alexey Kurakin, and Chun-Liang Li. 2020.  Fixmatch: Simplifying semi-supervised learning with consistency and confidence.  _NeurIPS_ (2020).
  * Togashi et al. (2021) Riku Togashi, Mayu Otani, and Shin’ichi Satoh. 2021.  Alleviating cold-start problems in recommendation through pseudo-labelling over knowledge graph. In _Proceedings of the 14th ACM international conference on web search and data mining_. 931–939.
  * Volkovs et al. (2017) Maksims Volkovs, Guangwei Yu, and Tomi Poutanen. 2017.  Dropoutnet: Addressing cold start in recommender systems.  _NeurIPS_ (2017).
  * Wang et al. (2021a) Wenjie Wang, Fuli Feng, Xiangnan He, Liqiang Nie, and Tat-Seng Chua. 2021a.  Denoising implicit feedback for recommendation. In _Proceedings of the 14th ACM international conference on web search and data mining_. 373–381.
  * Wang et al. (2021b) Zitai Wang, Qianqian Xu, Zhiyong Yang, Xiaochun Cao, and Qingming Huang. 2021b.  Implicit Feedbacks are Not Always Favorable: Iterative Relabeled One-Class Collaborative Filtering against Noisy Interactions. In _Proceedings of the 29th ACM International Conference on Multimedia_. 3070–3078.
  * Wei et al. (2021) Yinwei Wei, Xiang Wang, Qi Li, Liqiang Nie, Yan Li, Xuanping Li, and Tat-Seng Chua. 2021\.  Contrastive learning for cold-start recommendation. In _Proceedings of the 29th ACM International Conference on Multimedia_. 5382–5390.
  * Wu et al. (2021) Jiancan Wu, Xiang Wang, Fuli Feng, Xiangnan He, Liang Chen, Jianxun Lian, and Xing Xie. 2021.  Self-supervised graph learning for recommendation. In _SIGIR_.
  * Yin et al. (2012) Hongzhi Yin, Bin Cui, Jing Li, Junjie Yao, and Chen Chen. 2012.  Challenging the long tail recommendation.  _arXiv preprint arXiv:1205.6700_ (2012).
  * Zhang et al. (2019) Jiani Zhang, Xingjian Shi, Shenglin Zhao, and Irwin King. 2019\.  Star-gcn: Stacked and reconstructed graph convolutional networks for recommender systems.  _arXiv preprint arXiv:1905.13129_ (2019).
  * Zhang et al. (2021) Yang Zhang, Fuli Feng, Xiangnan He, Tianxin Wei, Chonggang Song, Guohui Ling, and Yongdong Zhang. 2021.  Causal intervention for leveraging popularity bias in recommendation. In _Proceedings of the 44th International ACM SIGIR Conference on Research and Development in Information Retrieval_. 11–20.
  * Zhao et al. (2022) Minghao Zhao, Le Wu, Yile Liang, Lei Chen, Jian Zhang, Qilin Deng, Kai Wang, Xudong Shen, Tangjie Lv, and Runze Wu. 2022\.  Investigating Accuracy-Novelty Performance for Graph-based Collaborative Filtering. In _Proceedings of the 45th International ACM SIGIR Conference on Research and Development in Information Retrieval_. 50–59.
  * Zhu et al. (2021) Yongchun Zhu, Ruobing Xie, Fuzhen Zhuang, Kaikai Ge, Ying Sun, Xu Zhang, Leyu Lin, and Juan Cao. 2021.  Learning to warm up cold item embeddings for cold-start recommendation with meta scaling and shifting networks. In _Proceedings of the 44th International ACM SIGIR Conference on Research and Development in Information Retrieval_. 1167–1176.


# Closures for Multi-Component Reacting Flows
based on Dispersion Analysis

Omkar B. Shende<EMAIL_ADDRESS>Ali Mani<EMAIL_ADDRESS>Department of Mechanical Engineering, Stanford University, Stanford, CA 94035,
USA

###### Abstract

This work presents algebraic closure models associated with advective
transport and nonlinear reactions in a Reynolds-averaged Navier-Stokes context
for a system of species subject to binary reactions and transport by advection
and diffusion. Expanding upon analysis originally developed for non-reactive
transport in the context of Taylor dispersion of scalars, this work extends
the modified gradient diffusion model explicated by Peters (Turbulent
Combustion, 2000) and based on work by Corrsin (JFM, vol. 11, p.407-416)
beyond single-component transport phenomena and involving nonlinear reactions.
The presented model forms, from this weakly-nonlinear extension of the
original dispersion theory, lead to an analytic expression for the eddy
diffusivity matrix that explicitly captures the influence of the reaction
kinetics on the closure operators. Furthermore, we demonstrate that the
derived model form directly translates between flow topologies through a
priori and a posteriori testing of a binary species system subject to
homogeneous isotropic turbulence. Using two- and three-dimensional direct
numerical simulations involving laminar and turbulent flows, it is shown that
this framework improves prediction of mean quantities compared to previous
results. Lastly, the presented model form, collapses to the earlier gradient
diffusion and its modified version derived by Corrsin in the limits of non-
reactive species and linear reactions, respectively.

##  I Introduction

Building reduced-order models for turbulent reacting flows is theoretically
and computationally challenging as the underlying chemical and transport
processes are individually complex and a thorough understanding of the coupled
effects of these phenomena remains elusive. However, deeper insight into the
mechanisms by which turbulent transport and reaction dynamics influence each
other is essential for the future design of efficient systems.

In considering the evolution of the concentration of passive scalars in an
incompressible flow, a one-way coupling between the fluid momentum and the
scalars is typically assumed in writing the transport equations. In
particular, if considering the concentrations of two scalars involved in a
binary reaction system, the concentrations being denoted $C_{1}$ and $C_{2}$,
in an imposed solenoidal velocity field, $U$, the governing equations can be
formulated as

$\frac{\partial C_{i}}{\partial
t}+\nabla\cdot\left(UC_{i}\right)=D_{m}\nabla^{2}C_{i}-AC_{i}C_{j\neq i},$ (1)

where $A$ is a reaction coefficient and $D_{m}$ is a scalar diffusivity which
is assumed to be the same for both indexed scalars in this work. This equation
is valid in the dilute limit, where the heat release from the reaction is so
small as to affect neither the density nor the reaction kinetics. [1]

While this form of the equation incorporates all effects in the system, in
practice, the Reynolds-averaged version of equation 1, written as

$\frac{\partial\overline{C_{i}}}{\partial
t}+\nabla\cdot\left(\overline{U}_{\,}\overline{C}_{i}\right)+\nabla\cdot\left(\overline{U^{\prime}C_{i}^{\prime}}\right)=D_{m}\nabla^{2}\overline{C_{i}}-A\overline{C}_{i}\overline{C}_{j\neq
i}-A\overline{C^{\prime}_{i}C^{\prime}_{j\neq i}},$ (2)

has more utility due to its lower computational cost. In Eqn. 2 primed terms
denote fluctuations about the mean quantities, which are denoted by over-bars.
In general, for a stochastic system, this averaging is done over multiple
ensembles of flow realizations, and the corresponding set of partial
differential equations for the velocity field are the Reynolds-averaged Navier
Stokes (RANS) equations. Finding models for the two types of unclosed terms in
Eqn. 2, $\overline{U^{\prime}C_{i}^{\prime}}$ and
$\overline{C^{\prime}_{i}C^{\prime}_{j\neq i}}$, is an important problem in
predictive science and for applications involving combustion [2], atmospheric
and marine pollution [3], electrolyte solutions in electrochemical
applications [4], and other fields. In this work, such a model will be derived
using a weakly-nonlinear extension of dispersion analysis.

###  I.1 Taylor Dispersion

For non-reacting scalar contaminants, many models exist for tracking the
dispersion of passive scalars. The most common analytic model invoked for the
unclosed transport term in the non-reactive equations, which is also used for
some reactive cases, follows the generalized gradient diffusion hypothesis,
for which the unresolved flux of a single passive scalar is of the form

$\overline{u_{i}^{\prime}C_{1}^{\prime}}=-D_{eff}\frac{\partial\overline{C_{1}}}{\partial
x_{i}},$ (3)

where $D_{eff}$ is a generic effective eddy diffusivity that aggregates
turbulent transport effects. When there is strong scale separation between the
mean fields and the underlying fluctuating velocity fields, this local model
is exact, with the caveat that $D_{eff}$ may be a non-isotropic tensor when
the underlying fluctuating flow is not statistically isotropic. [5, 6]

In the framework of this model form, there is a rich body of literature on
analytic models that find this effective diffusivity in canonical flows. In
this study we denote the effective diffusivity in the absence of reaction as
$D^{0}$, and refer to it as the ”standard” eddy diffusivity. Correspondingly,
we denote $D_{eff}=D^{0}$ as referring to the standard gradient diffusion
model, which is agnostic to the reactivity of scalars.

There are many ways to measure this defining parameter. In particular, for the
case of parallel flow in a pipe, [7] shows that a passive scalar experiences
enhanced diffusive transport in the longitudinal direction. We can define a
Péclet number as $Pe\sim uL/D_{m}$ such that, at high values of the Péclet
number, $D^{0}$ dominates the standard molecular diffusivity, $D_{m}$. In
parallel flows, it is generally shown that this quantity scales as

$\frac{D_{0}}{D_{m}}\sim\frac{u^{2}L^{2}}{D_{m}^{2}},$ (4)

where $u$ is the characteristic convective speed, and $L$ is some
characteristic spanwise length scale for the flow. This idea can be expanded
to applications outside of pipe flow, and other studies have expanded the
range of problems and fields to ones that Taylor’s original analysis does not
address. [8, 9, 10, 11]

###  I.2 Other Scalar Transport Models

Thus far, we have only looked at non-reactive models. However, many tracers of
interest are not passive; to reach a multi-component framework, the first step
is the simplest possible reaction involving two species, where the
concentration of the second species is orders of magnitude larger than the
first everywhere in the domain. In some applications, this allows an
assumption that the second species is maintained in the domain at a constant
concentration and unaffected by the reaction dynamics, leading to a first-
order linear reaction term for the first species. In the presence of such a
first-order linear reaction, [2] used work by [12] to propose a model that
incorporates the chemistry effects for the unclosed scalar flux term. It can
be written as

$\overline{u_{i}^{\prime}C_{1}^{\prime}}=-\frac{D^{0}}{1+A\overline{C_{2}}\tau_{mix}}\frac{\partial\overline{C_{1}}}{\partial
x_{i}}.$ (5)

In this model, $\tau_{mix}$ is a characteristic time scale for mixing by the
flow. In the denominator, one can form a Damköhler number, relating the
timescale of the reaction to the turbulence timescale, and can therefore write
the flux as the gradient times an effective diffusivity that scales as

$D_{eff}\sim\frac{D_{0}}{1+Da},$ (6)

where $D^{0}$ is the aforementioned ”standard” turbulent diffusivity. When the
reaction rate is sufficiently slow or when examining the non-reacting limit
where $Da$ vanishes, this recovers the gradient diffusion model exactly. While
$D^{0}$ still needs to be determined from other means, the advantage of Eqn. 6
is that it analytically captures the impact of reaction on $D_{eff}$.

Equation 6 contributes critical understanding by explicitly revealing the
dependence of effective diffusivity (or eddy diffusivity) on the reaction
kinetics and, specifically, the reaction coefficient. Although transport and
kinetics are controlled by independent terms at the microscopic continuum
level, when behavior at a macroscopic level, governed by ensemble-averaged
equations, is sought, kinetics influence transport.

In this study, we develop a model in the Reynolds-averaged context that can
provide algebraic closures to the scalar evolution equations for a binary
reactant setup. We do this using the spirit of Taylor’s dispersion analysis.
By considering the case of a parallel flow as an analytical prototype, we
extend Taylor’s analysis to discover model forms for dispersion of scalars
undergoing binary reactions. Given the nonlinearity of this system, the
closure operators are expected to be dependent on the transported quantity
itself. We derive these dependencies through a weakly nonlinear extension of
Taylor’s original analysis. We then hypothesize that reactive mixing by
turbulent flows should involve model forms with similar nonlinearities, but
with different model coefficients. This is analogous to the situation of
nonreactive scalars, where mixing by both classes of flows can be reasonably
captured by gradient diffusion models. However, the diffusivity coefficients
differ based on the underlying flow topology, with $D^{0}/D_{M}$ scaling as
$Pe^{2}$ in parallel flows and as $Pe^{1}$ in turbulent flows.

Our weakly nonlinear model offers improvements to perturbation expansions
based on directly linearizing the equations performed by works like [13]. It
also recovers the scaling relationship between diffusivity and Damköhler
number derived in [14] and [15], although those works use nonlinear first-
order reactions. It also indicates that the effective eddy diffusivity is
linked to reaction, as has been shown in experiments. [16] This relationship
shows that the impacts of reaction chemistry and turbulent mixing are
fundamentally linked.

In what follows, we will first develop a model problem based on laminar flows
that can capture the relevant multi-physics at play for a linear reaction
setup and work to find a solution for the single-scale and the multi-scale
problems. This insight will guide generalization to a binary reaction problem,
which we can fully expand to general turbulent flows. The outcome of our
analysis is a unifying reduced-order model (ROM) in the RANS context that
captures the analytic impact of chemical kinetics on transport and reaction
closures while leveraging existing literature on non-reactive scalar transport
modeling.

##  II Model Problem

Figure 1: A schematic of the model problem, illustrating dimensions and
boundary conditions. The flow is depicted in black as $u^{\prime}(y)$, with a
representation of a reaction front given at the mid-plane. Axial will refer to
the $x$-direction, while spanwise will refer to the $y$-direction

In order to develop insights into a model form that captures closure terms in
a binary mixture, we propose the following illustrative setup that can be
treated semi-analytically. Consider a flow field that is parallel, steady, and
two-dimensional, contained in a domain that is an elongated box. We place
periodic boundary conditions on the two opposing elongated sides, as seen in
Fig. 1. By elongated, we mean that $L_{1}>>L_{2}$.

In such a two-dimensional domain, we define an averaging operator as

$\overline{f(x)}=\frac{1}{L_{2}}\int_{0}^{L_{2}}f(x,y)dy,$ (7)

which differs from the ensemble averaged used earlier. In this context, let us
use a periodic structure for the flow-field which imposes no mean velocity.
Namely, $u(x,y)=\overline{u}+u^{\prime}=U_{0}sin(ky)$, where $k=2\pi/L_{2}$.
At $x=-L1/2$ and $x=L_{1}/2$, the two domain boundaries, we prescribe
Dirichlet conditions for the scalar concentrations.

In addition, we consider a reaction of the form $C_{1}+C_{2}\rightarrow
C_{3}$, which is considered irreversible and fully activated, but follows a
finite rate law of mass action. While simpler than a fully three-dimensional
turbulent flow with reacting scalars, this model problem captures the
essential competing physics of mixing, transport, and reactions that are
present in the more realistic case. However, we can now formulate solutions in
this sandbox context and propose closures for each sub-problem that each
highlight different aspects of the overall ROM.

###  II.1 Linear Reaction, Single-Scale Flow

Let us first examine this problem in the context of a linear reaction. In this
setup, the concentration of one of the scalars, $C_{2}$, is held constant in
the entire domain at unity and does not evolve, while the maximum value of
$C_{1}$, at the left side of the domain, is held to be at least an order of
magnitude less than that of $C_{2}$. As $C_{2}$ is a constant, we will denote
the reaction rate as $A_{L}=AC_{2}$, denoting the product of a standard binary
species reaction coefficient times the concentration of the constant species.
The fluctuations for $C_{1}$ are governed by a transport equation as

$\frac{\partial C_{1}^{\prime}}{\partial
t}+\frac{\partial(u^{\prime}C_{1}^{\prime})^{\prime}}{\partial
x_{1}}+u^{\prime}(y)\frac{\partial\overline{C_{1}}}{\partial
x_{1}}=D_{m}\frac{\partial^{2}C_{1}^{\prime}}{\partial x_{i}\partial
x_{i}}-A_{L}C_{1}^{\prime},$ (8)

where $C_{1}=\overline{C_{1}}+C_{1}^{\prime}$ and so on. This equation can be
derived by finding the classical transport equation for $C_{1}$ and
subtracting from it the evolution equation for $\overline{C_{1}}$. Note that
we have implicitly assumed here that the molecular diffusivity, $D_{m}$, is
the same for all species, is constant, and is isotropic.

This equation is similar to the equation governing evolution of fluctuating
fields in the work of [7], with the primary distinction being the reaction
term on the right hand side. While Taylor was not explicitly considering this
type of periodic parallel flow in his work, the steps he follows allow us to
also simplify this equation to arrive at a more analytically approachable
differential equation.

First, we note that we are considering a domain where $L_{1}>>L_{2}$, from
which we can conclude that diffusion in the axial direction, $x$, can be
neglected when compared to diffusion in the spanwise direction, $y$. This can
be shown by performing a scaling analysis and seeing the associated term tends
to zero as the relevant Péclet number tends to infinity. This allows us to
neglect axial diffusion in Eqn. 8.

In addition, when $L_{1}>>L_{2}$, the time scale of mixing in the axial
direction is much larger than the time scale of mixing in the spanwise
direction. Therefore, one can decompose the concentration field as
$C_{1}(x,y)=\overline{C_{1}}(x)+C_{1}^{\prime}(x,y)$, allowing for evolution
of the mean field in the axial direction but considering the underlying field
to be ”well mixed” in the spanwise direction. In particular, the presence of
fast mixing in the spanwise direction implies that
$|C_{1}^{\prime}|<<|\Delta\overline{C_{1}}|$, where $|\Delta\overline{C_{1}}|$
is a characteristic change in the averaged concentration in the axial
direction.

This final conclusion allows us to see that the second term on the left-hand
side of Eqn. 8 must be much less significant than the advection of the mean
scalar concentration, represented by the third term on the left-hand side.
Therefore, for a leading-order approximation, we can neglect the former,
double primed term.

All of these physical assumptions imply that the flow is in an advection-
driven limit. The strong mixing in the spanwise direction that such a flow
implies also allows us to invoke a quasi-steady approximation for the
evolution of the fluctuating concentration. In particular, since we are most
interested in the long-time response, not the initial transient one, this
assumption does not unreasonably limit our analysis.

If the previous simplifications are applied to Eqn. 8, the dominant remaining
terms give us an approximate evolution equation that can be written as

$u^{\prime}(y)\frac{\partial\overline{C_{1}}}{\partial
x}=D_{m}\frac{\partial^{2}C_{1}^{\prime}}{\partial
y^{2}}-A_{L}C_{1}^{\prime}.$ (9)

Figure 2: Predictions of closures for the linear reaction problem subject to a
turbulent flow, $X=x/L_{1}$. The Linear Reaction model of Eqn. 5 is identical
to the ROM represented by Eqn. 12, and the former is overlaid by the latter in
the plot.

This gives us the dispersion of the fluctuating quantity as driven by the
flow. A further ansatz we make here is that the form of the fluctuating
components follows the velocity, as
$C_{1}(x,y)=\overline{C_{1}}(x)+f_{1}(x)sin(ky)$. That is, the fluctuating
field has some magnitude determined purely by axial position. This is
motivated by examining solely the leading-order term in the harmonic expansion
of the concentration fluctuations in the $y$-direction. Substituting this into
the fluctuation equation yields

$D_{m}k^{2}f_{1}(x)+A_{L}f_{1}(x)=-U_{0}\frac{\partial\overline{C_{1}}}{\partial
x}.$ (10)

Solving this equation for our unknown axial concentration magnitude yields

$f_{i}(x)=-\frac{U_{0}}{k^{2}D_{m}+A_{L}}\frac{\partial\overline{C_{i}}}{\partial
x}.$ (11)

The only unclosed term in Eqn. 8 is the scalar flux, and it can now be solved
as

$\overline{u^{\prime}C_{1}^{\prime}}=\overline{U_{0}f_{1}(x)sin^{2}(ky)}=-D_{eff}\frac{\partial\overline{C_{1}}}{\partial
x}.$ (12)

Here, $D_{eff}$ is given by

$D_{eff}=\frac{D_{0}}{1+A_{L}\tau_{mix}},$ (13)

where the standard eddy diffusivity, $D^{0}$, is determined from non-reactive
Taylor-type analysis. For this flow, it is given by

$D^{0}=\frac{U_{0}^{2}}{2D_{m}k^{2}}.$ (14)

In addition, $\tau_{mix}$ is a characteristic timescale for mixing given by

$\tau_{mix}=\frac{1}{D_{m}k^{2}}=\frac{D_{0}}{u_{rms}^{2}},$ (15)

where $u_{rms}$ is the root-mean-squared velocity of the underlying flow,
equal to $U_{0}/\sqrt{2}$ in this derivation. For this section, we will refer
to Eqn. 5 as the Linear Reaction model and Eqn. 13 as a special case of our
Reduced-Order Model (ROM) for this first-order reaction context.

Interestingly, for this problem, the ROM is identical to Eqn. 5, even though
it is derived in a laminar context, while the latter was created in the
framework of turbulent flows. It again has an effective Damköhler number in
the denominator, and shows that effective diffusivity decreases as the
reaction rate increases relative to the flow timescale. This lends confidence
that the choice of model problem should not hinder the relevance of the
conclusions here in application to turbulence.

In particular, the model form is an explicit function of $A_{L}$, with the
other two quantities, $D^{0}$ and $u_{rms}$, defined purely in terms of flow
parameters. This means that they can be derived independent of the reaction
kinetics, in non-reacting flow, and applied directly to this model for
reacting flows.

To test this model, we next examine its predictions against 2D DNS data as
well as the gradient diffusion model that ignores impact of reactions on the
advective closure. For this purpose, we consider a setting with
$Pe\equiv(U_{0}L_{2})/(2\pi D_{m})=100$, and $Da\equiv A_{L}\tau_{mix}=1$. For
the axial boundary conditions, we used Dirichlet conditions $x=-L_{1}/2$ and
$x=L_{1}/s$ by respectively enforcing $C_{1}$ to be equal to $0.1$ and $0$.

Figure 2 presents quantitatively a comparison between different results. For
this specific setting, the mean concentration assumes an exponential profile,
as seen in the figure. It can be seen that ROM produces results much closure
to the DNS than the gradient diffusion model. The gap between DNS and gradient
diffusion model is solely due to over prediction of $D_{eff}$ by the model.
For settings at higher $Da$ the correction offered by ROM becomes even more
crucial as this gap will be even larger.

###  II.2 Linear Reaction, Multi-Scale Flow

The previous section only considered a velocity field with a single relevant
length scale. However, a realistic flow has a spectrum of relevant scales and
the model problem setup allows us to consider a straightforward superposition
of multiple velocity fields. The total velocity now can be represented by
$u(x,y)=\sum U_{n}sin(nky)$, where $k=2\pi/L_{2}$ and $n$ refers to each scale
as corresponding to an integer multiple of the fundamental scale set by $k$.
Following the same procedure as the previous section, the unresolved scalar
flux can be expressed as

$\overline{u^{\prime}C_{1}^{\prime}}=-\sum_{n=1}^{\infty}\frac{U_{n}^{2}/2}{n^{2}k^{2}D_{m}+A_{L}}\frac{\partial\overline{C_{1}}}{\partial
x}.$ (16)

In this context, the dependence of $D_{eff}$ on $A_{L}$ is more complex.
Nevertheless, we still seek to express the dependence of $D_{eff}$ on $A$ in
terms of a simple expression that involves low order statistics, specifically
$D^{0}$ and $u_{rms}$. In order to do so, we first examine the expression
above in the asymptotic limits of large and small $A_{L}$. In the limit of
$A_{L}$ equaling zero, the above expression results in $D_{eff}$ equal to the
$D^{0}$ of the multi-scale flow field, recovering the non-reactive solution
from pure Taylor-type analysis. In the limit of large $A_{L}$, however, we
obtain that $D_{eff}=u_{rms}^{2}/A_{L}$.

Figure 3: Effective diffusivity as a function of $Da$. Each marker refers to a
different velocity field, with the ”Multiscale” field given by Eqn. 20. The
solid line represents Eqn. 17.

Interestingly, examining Eqn. 13 and updating values for $D^{0}$ and
$\tau_{mix}$ satisfies both limits exactly.

$D_{eff}=\frac{D^{0}}{1+A_{L}\tau_{mix}},$ (17)

where now

$D^{0}=\sum_{n=1}^{\infty}\frac{U_{n}^{2}/2}{n^{2}k^{2}D_{m}}$ (18)

and $\tau_{mix}$, the mixing time, is $D^{0}/u_{rms}^{2}$, where

$u^{2}_{rms}=\sum_{n=1}^{\infty}U_{n}^{2}.$ (19)

We reiterate that both $D^{0}$ and $u_{rms}$, and therefore $\tau_{mix}$, are
fundamental measures of the flow that are independent of the presence of the
reaction, and there is no $A_{L}$ dependence in these quantities.

One can examine this model by comparing the $D_{eff}$ it provides against
those from Eqn. 16, as shown in Fig.3. Among the various velocity profiles
tested, we considered the the following velocity profile, denoted as “Multi-
scale” in Figure 3, expressed as

$u^{\prime}(y)=\sum_{n=1}^{10}n^{-5/6}cos(nky),$ (20)

to mimic the multi-scale effects found in turbulence. The specific power,
$-5/6$, is chosen such that the eddy diffusivity at each wavenumber matches
that of Kolmogorov turbulence.

In Fig. 3, the simple model presented in Eqn. 17 is plotted across values of
the reaction rate. We see that it captures the precise behavior of $D_{eff}$
in the limits of small and large $A_{L}$ across a variety of imposed velocity
fields. Additionally, the model still fits the data reasonably in the
intermediate regimes where $A_{L}\tau_{mix}\sim O(1)$.

###  II.3 Binary Reaction

Having some confidence with the linear problem, we now turn back to examining
Eqn. 2 in its entirety, and note that the complete problem involves a full
binary system, where both $C_{1}$ and $C_{2}$ are free to evolve.

Before proceeding with derivation of a closure model, let us first test
whether we can naively modify Eqn. 5 to perform reasonably in the case of
binary reactions. Specifically, one may interpret that for evolution of each
species, $C_{i}$, the other species might be assumed ”locally constant.” This
allows interpreting an effective $A_{L}$ for each species as
$A_{L}=A\overline{C}_{j\neq i}$. With this simple modification, the previous
model form can be extended to a multi-scale model as

$\overline{u^{\prime}C_{i}^{\prime}}=-\frac{D}{1+A\tau\overline{C_{j}}}\frac{\partial\overline{C_{i}}}{\partial
x}.$ (21)

where $A_{L}$ is now replaced with the $A\overline{C_{2}}$ and
$A\overline{C_{1}}$ for mean advective fluxes associated with $C_{1}$ and
$C_{2}$, respectively.

We next examine the performance of this model form by comparing its prediction
against direct numerical simulation. In this case we consider a single-scale
velocity profile where $u(y)=U_{0}sin(ky)$, but where both reactants are now
free to vary. We choose $Pe\equiv(U_{0}L_{2})/(2\pi D_{m})=100$ and $Da\equiv
A\tau_{mix}C_{ref}=100$, where $C_{ref}$ is concentration imposed by Dirichlet
boundary conditions on each scalar. More specifically, the boundary condition
on $C_{1}=C_{ref}$ is imposed at $x=-L_{1}/2$ and $C_{2}=C_{ref}$ is imposed
at $x=L_{1}/2$. The other boundary condition for each scalar is a homogeneous
Dirichlet condition.

In Fig. 4, we see that the Linear Reaction model given by Equation 21 incurs
great errors near the middle of the domain, where $C_{1}$ and $C_{2}$ exist in
near-stoichiometric ratios. Surprisingly, the standard Gradient Diffusion
model, which carries no chemistry information, now outperforms this more
physics-based model.

Figure 4: Predictions of closures for the binary reaction problem subject to a
laminar flow, $X=x/L_{1}$.

Clearly, the Linear Reaction model is not sufficient to describe this system,
so let us return to Eqn. 2 and examine the full system of equations again. In
this framework, the governing equations for the fluctuations control the
evolution of the unclosed terms of interest. For $C_{1}$, this looks like

$\frac{\partial C_{1}^{\prime}}{\partial
t}+\frac{\partial(u^{\prime}C_{1}^{\prime})^{\prime}}{\partial
x_{1}}+u^{\prime}(y)\frac{\partial\overline{C_{1}}}{\partial
x_{1}}=D_{m}\frac{\partial^{2}C_{1}^{\prime}}{\partial x_{i}\partial x_{i}}\\\
-A\overline{C_{1}}C_{2}^{\prime}-A\overline{C_{2}}C_{1}^{\prime}\\\
+A(C_{1}^{\prime}C_{2}^{\prime})^{\prime}.$ (22)

Again following the work of [7], we can add one more assumption for the binary
system that was not relevant for the linear case, which is that the
fluctuations of the product of fluctuating quantities will be negligible when
compared to the mean field quantities.

If those simplifications are applied to Eqn. 22 for $C_{1}$, it can be
reformulated as

$u^{\prime}(y)\frac{\partial\overline{C_{1}}}{\partial
x}=D_{m}\frac{\partial^{2}C_{1}^{\prime}}{\partial
y^{2}}-A\overline{C_{1}}C_{2}^{\prime}-A\overline{C_{2}}C_{1}^{\prime},$ (23)

and an equivalent equation can be derived for $C_{2}$.

We again make the ansatz that
$C_{i}(x,y)=\overline{C_{i}}(x)+f_{i}(x)sin(ky)$, following the velocity.
Substituting this into the fluctuation equations yields the following system
of two equations:

$D_{m}k^{2}f_{1}(x)+A\overline{C_{1}}f_{2}(x)+A\overline{C_{2}}f_{1}(x)=-U_{0}\frac{\partial\overline{C_{1}}}{\partial
x}$ (24)

and

$D_{m}k^{2}f_{2}(x)+A\overline{C_{2}}f_{1}(x)+A\overline{C_{1}}f_{2}(x)=-U_{0}\frac{\partial\overline{C_{2}}}{\partial
x}.$ (25)

This is a coupled, but linear, set of equations for $f_{1}$ and $f_{2}$ that
can be solved to find that

$f_{1}(x)=\frac{-U_{0}(k^{2}D_{m}+A\overline{C_{1}})}{k^{2}D_{m}(k^{2}D_{m}+A\overline{C_{1}}+A\overline{C_{2}})}\frac{\partial\overline{C_{1}}}{\partial
x}+\frac{U_{0}(A\overline{C_{1}})}{k^{2}D_{m}(k^{2}D_{m}+A\overline{C_{1}}+A\overline{C_{2}})}\frac{\partial\overline{C_{2}}}{\partial
x}$ (26)

and

$f_{2}(x)=\frac{-U_{0}(k^{2}D_{m}+A\overline{C_{2}})}{k^{2}D_{m}(k^{2}D_{m}+A\overline{C_{1}}+A\overline{C_{2}})}\frac{\partial\overline{C_{2}}}{\partial
x}+\frac{U_{0}(A\overline{C_{2}})}{k^{2}D_{m}(k^{2}D_{m}+A\overline{C_{1}}+A\overline{C_{2}})}\frac{\partial\overline{C_{1}}}{\partial
x}.$ (27)

We can also write these terms, the result of using a weakly-nonlinear
extension of dispersion analysis, as the superposition of two gradients
multiplied by a prefactor representing an effective diffusivity. The equations
for the scalar flux now look like

$\begin{bmatrix}\overline{u^{\prime}C_{1}^{\prime}}\\\
\overline{u^{\prime}C_{2}^{\prime}}\end{bmatrix}=-\begin{bmatrix}D_{11}&D_{12}\\\
D_{21}&D_{22}\end{bmatrix}\frac{\partial}{\partial
x}\begin{bmatrix}\overline{C_{1}}\\\ \overline{C_{2}}\end{bmatrix},$ (28)

where $D_{kl}$ represents a diffusivity associated with the flux of the $k$-th
species due to gradient in the species. The square matrix formed by these
coefficients is also positive definite so long as the mean scalar fields
maintain their positivity. The coefficients are written as,

$\displaystyle\begin{split}&D_{11}=\frac{D^{0}(1+A\overline{C_{1}}\tau_{mix})}{1+A\tau_{mix}\overline{C_{1}}+A\tau_{mix}\overline{C_{2}}}\;\;\;\;\;\;D_{12}=-\frac{D^{0}A\tau_{mix}\overline{C_{1}}}{1+A\tau_{mix}\overline{C_{1}}+A\tau_{mix}\overline{C_{2}}}\\\
&D_{21}=-\frac{D^{0}A\tau_{mix}\overline{C_{2}}}{1+A\tau_{mix}\overline{C_{1}}+A\tau_{mix}\overline{C_{2}}}\;\;\;\;\;\;D_{22}=\frac{D^{0}(1+A\overline{C_{2}}\tau_{mix})}{1+A\tau_{mix}\overline{C_{1}}+A\tau_{mix}\overline{C_{2}}}\end{split}$
(29)

using the notation from Section II.1.

Perhaps most novel, however, is that determining the actual local form of the
fluctuations means it is now possible to close not only the scalar transport
term, but also the unclosed reaction terms in the standard RANS equations.
This closure equation looks like

$\overline{AC_{1}^{\prime}C_{2}^{\prime}}=\frac{A}{u_{rms}^{2}}\left(D_{11}\frac{\partial\overline{C_{1}}}{\partial
x}+D_{12}\frac{\partial\overline{C_{2}}}{\partial
x}\right)\left(D_{21}\frac{\partial\overline{C_{1}}}{\partial
x}+D_{22}\frac{\partial\overline{C_{2}}}{\partial x}\right).$ (30)

This provides our Reduced-Order Model (ROM) in its most complete form. It
solves a reduced-degree-of-freedom system resulting from the RANS equations.
The three algebraic components describe an entire model form, and require only
prescriptions of parameters of the underlying flow, specifically $D^{0}$ and
$u_{rms}$. It is important to note that this model, under the appropriate
limits, captures the Gradient Diffusion model and the Linear Reaction as its
special cases. The Gradient Diffusion model is realized in the limit of
$A\rightarrow 0$. The Linear Reaction model, for finite $A_{L}$, can be
realized in the limit of $\overline{C_{2}}\rightarrow\infty$ and $A\rightarrow
0$.

Comparing this model against the DNS data presented in Fig.4, we see that the
new model recovers from the errors of the Linear Reaction model, and
outperforms the Gradient Diffusion model.

The multi-scale analysis performed in Section II.2 is omitted here as it would
yield a similar summation form as in that section, but as seen in the linear
regime, a model form from the single-scale case is still valid.

##  III Application to Turbulent Flows

Now that we have built up the full ROM, we can turn to realistic turbulent
flows. A great strength of the approach used in the previous section is that
we can borrow the model form from the laminar analysis directly and use it in
an analogous domain of homogeneous, isotropic turbulence. For this flow
configuration, as a recapitulation, the ROM is given by Eqns. 28 through 30,
with the only difference that flow parameters $D^{0}$ and $\tau_{mix}$ are now
properties of a turbulent flow. As the ROM captures reaction-dependent effects
explicitly, these parameters are agnostic to the presence of a reacting system
and non-reactive flow measurements provide the values needed.

Specifically, [6] explicates the process of finding these values for non-
reactive turbulent flows using the Macroscopic Forcing Method (MFM). In
particular, [17] used MFM to measure the true non-local form of the
macroscopic eddy diffusivity, and additionally found values for $D^{0}$ in the
purely local and isotropic limit applicable for this work.

###  III.1 A priori analysis

(a)

(b)

Figure 5: Instantaneous visualization of a binary reaction system subject to
homogeneous, isotropic turbulence (a) An isosurface of the instantaneous $Q$
criterion that indicates eddy cores (b) An isosurface of the instantaneous
reaction term, $AC_{1}C_{2}$, colored by axial coordinate. In both figures,
axes are equally scaled.

In particular, we replace the steady, parallel flow examined in the previous
section with three dimensional (3D) homogeneous isotropic turbulence (HIT) in
an elongated domain to generalize the model problem. The code of [18] in an
incompressible mode was adapted for this work to simulate HIT in a 3D domain
of size $(2\pi)^{3}$. The resulting flow field was periodically extended in
the $x$-direction to generate a computational domain for scalar transport with
dimensions $20\pi\times 2\pi\times 2\pi$. The scalar transport equations are
solved with boundary conditions identical to the 2D binary reaction problem.
This provides a realistic reaction zone in the middle of the domain, far from
the boundaries, as seen in Fig. 5.

For our flow simulations, we adopted parameters from [17] for a case of
$Re_{\lambda}=26$, which was obtained by setting the HIT box size of $L=2\pi$,
the turbulent forcing parameter following the prescription of [19] equal to
$0.2792$, and the kinematic viscosity to $\nu=0.0263$.

This setup allows adoption of the values for non-reactive eddy diffusivity.
For this specific case, [17] reported that
$D^{0}=0.86u_{rms}^{4}/\epsilon=0.96$ and
$\tau_{mix}=0.86u_{rms}^{2}/\epsilon=1.02$, where $\epsilon=0.79$ is the
dissipation rate of turbulent kinetic energy and $u_{rms}=0.97$ is the single-
component root-mean-squared velocity.

For scalar transport, we consider molecular diffusivity $D_{m}=0.0263$,
matching $\nu$, $C_{ref}=1$ and a reaction rate $A=100$. This leads to
$Pe\equiv D^{0}/D_{m}=37$ and $Da\equiv A\tau_{mix}C_{ref}=102$.

Figure 5 shows imagery captured near the middle of the computational domain.
Figure 5a shows vigorous turbulent structures that mix the flow and the scalar
fields. In response, reaction fronts are not planar; instead, one can observe
highly stretched and distorted structures. Similar to the laminar problem,
however, the ensemble-averaged fields will be smooth and one-dimensional.

Figure 6: Predictions of closures for the binary reaction problem subject to a
turbulent flow. $X=x/L_{1}$ and DNS results are presented with 95% confidence
intervals accounting for statistical fluctuations.

In Fig.6, the closure terms for the ensemble-averaged binary reaction problem
are examined in the defined turbulent context. Similar to the laminar binary
reaction problem, the full ROM recovers from the errors in capturing the true
transport closures that are incurred by the Linear Reaction model of Eqn. 21.
For these specific closure terms, the standard Gradient Diffusion model
appears to match the DNS data more closely than the ROM. However, an overall
assessment that considers both reaction and transport closures reveals the
advantage of the ROM, as the Gradient Diffusion transport closure model offers
no answers to the equally vital reaction closure question. Quantifying this
advantage can be accomplished by performing a posteriori analysis and
comparing 1D model predictions against DNS data.

###  III.2 A posteriori analysis

(a)

(b)

Figure 7: Results of _a posteriori_ analysis, with DNS quantities indicated
with overbars, $X=x/L_{1}$. (a) A comparison of the predicted mean scalar
profiles in the reaction zone for the ROM and the Gradient Diffusion models
(b) A comparison of the predicted mean scalar profile for $C_{1}$ with
different aspects of the ROM active.

In this section, we solve Eqn. 2 directly by invoking some of the closure
models heretofore presented. The problem setup and parameters used in this
section are identical to those for the _a priori_ analysis. The only exception
is that for the RANS equation, instead of Dirichlet boundary conditions for
the scalars, Neumann boundary conditions with slopes matching the DNS profiles
are used. This is done because the DNS develops axial boundary layers near the
Dirichlet conditions due to local outflow advection. For the DNS, we have
ensured these artificial effects do not pollute the reaction zone results by
ensuring the boundaries are far from the reaction zone. These boundary layers
are absent in the RANS case as the mean velocity is zero. Appropriate matching
of RANS solutions to the DNS ones should consider concentration profiles
outside of these artificial boundary layers. We have done so by matching the
slopes of the RANS concentration profiles to those of the DNS outside of the
boundary layers, but far from the reaction zones.

In Fig. 7(a), we compare the predicted mean profiles of $C_{1}$ and $C_{2}$ to
the results derived from the DNS described in the previous section. It is
obvious that when all closures are considered, the ROM outperforms the
standard Gradient Diffusion model.

To demonstrate the relative importance of each of the closure terms, Fig.7(b)
shows the effects of toggling each of the closure terms on the overall
accuracy of scalar concentration predictions.

##  IV Discussion and Conclusions

In this work, we have introduced a framework to extend the analysis of the
dispersion of a single passive tracer to analysis of system of species
undergoing binary reactions. The resulting model form derived herein
introduces closures to both advective flux and reaction terms with nonlinear
interactions between the mean state of all species, even though the primary
set of equations being solved undergo a linearization procedure. The proposed
framework also captures the non-reactive and linear reaction limits as special
cases of the binary reaction problem. In contrast to [13], the derived model
leads to a positive-definite diffusion operator across all Damköhler numbers,
so numerical stability is maintained when the introduced ROM is used in a one-
dimensional scalar equation solver.

By considering four different test cases, we showed that the standard Gradient
Diffusion and Linear Reaction model do not suitably capture the behaviour of
mean scalar fields across all problems. The ROM, however, can address all
these different regimes consistently. Perhaps the most novel aspect of this
work is the direct translation of the obtained model form between the laminar
prototype problem and the test problem involving turbulent flows. Most
critically, no tuning of parameters was used in this analysis, as methods like
MFM permit the measurement of $D^{0}$ directly for turbulent flows. However,
it is unclear whether omitting a constant prefactor in
$\tau_{mix}=D^{0}/u_{rms}^{2}$ is appropriate for turbulent settings. In the
laminar case, this definition with no need for an additional prefactor was
derived analytically. In the turbulent case, we simply adopted the laminar
definition without retuning its prefactor that could otherwise account for the
change in flow topology.

An interesting feature of the developed ROM is that advective fluxes of any
one reactant depend on mean gradients of other species involved in the
reaction. This leads to an effective diffusivity matrix that shows cross-
diffusion coefficients. In this sense, our results are consistent with the
recent work of [20] in which a two-dimensional system undergoing a slightly
different reaction mechanism was considered. However, the additional advantage
of the recent ROM is that it consistently models all multiphysics effects by
offering closure expressions to both reaction terms as well as transport
terms. Furthermore, the present ROM provides analytic model forms that capture
explicitly the influence of reaction constants on each closure operator. To
this end, only other inputs needed are the non-reactive eddy diffusivities and
the mixing time scale to form the appropriate Damköhler numbers.

It should be noted that this analysis is not identical to reacting flow
formulations that consider just a mixture fraction or related quantity, as in
[21], as their scalar. The mixture fraction maps both concentration fields in
a binary system to a single scalar that is conserved in space. This is in
contrast to the partial differential equations considered here.

Another promising avenue of further inquiry is applying the proposed ROM to a
Large-Eddy Simulation context as a subgrid-scale modes, where mean space
variables are replaced with filtered variables and quantities like $u_{rms}$
are evaluated at the grid scale. Additionally, future work may be focused on
capturing non-locality, as in [17], which should account for the discrepancy
of the ROM values and the DNS values.

##  V Acknowledgements

Support for this work was provided by the National Science Foundation Graduate
Research Fellowship Program, the Stanford Graduate Fellowships in Science and
Engineering, and National Science Foundation Extreme Science and Engineering
Discovery Environment resources.

## References

  * Avrin [1997] J. D. Avrin, Incompressible reacting flows, Transactions of the American Mathematical Society 349, 3875 (1997).
  * Peters [2000] N. Peters,  _Turbulent Combustion_, Cambridge Monographs on Mechanics (Cambridge University Press, 2000).
  * Seinfeld and Pandis [2016] J. H. Seinfeld and S. N. Pandis,  _Atmospheric chemistry and physics : from air pollution to climate change_ , third edition. ed. (John Wiley and Sons, Hoboken, New Jersey, 2016).
  * Rosales and Nava [2017] M. Rosales and J. L. Nava, Simulations of turbulent flow, mass transport, and tertiary current distribution on the cathode of a rotating cylinder electrode reactor in continuous operation mode during silver deposition, Journal of The Electrochemical Society 164, E3345 (2017).
  * Kraichnan [1987] R. H. Kraichnan, Eddy viscosity and diffusivity: exact formulas and approximations, Complex Systems 1, 805 (1987).
  * Mani and Park [2021] A. Mani and D. Park, Macroscopic forcing method: A tool for turbulence modeling and analysis of closures, Phys. Rev. Fluids 6, 054607 (2021).
  * Taylor [1953] G. I. Taylor, Dispersion of soluble matter in solvent flowing slowly through a tube, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 219, 186 (1953).
  * Taylor [1954] G. I. Taylor, The dispersion of matter in turbulent flow through a pipe, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 223, 446 (1954).
  * Aris [1956] R. Aris, On the dispersion of a solute in a fluid flowing through a tube, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 235, 67 (1956).
  * Lighthill [1966] M. J. Lighthill, Initial Development of Diffusion in Poiseuille Flow, IMA Journal of Applied Mathematics 2, 97 (1966).
  * Shapiro and Brenner [1986] M. Shapiro and H. Brenner, Taylor dispersion of chemically reactive species: Irreversible first-order reactions in bulk and on boundaries, Chemical Engineering Science 41, 1417 (1986).
  * Corrsin [1961] S. Corrsin, The reactant concentration spectrum in turbulent mixing with a first-order reaction, Journal of Fluid Mechanics 11, 407–416 (1961).
  * Battiato and Tartakovsky [2011] I. Battiato and D. Tartakovsky, Applicability regimes for macroscopic models of reactive transport in porous media, Journal of Contaminant Hydrology 120-121, 18 (2011).
  * Elperin  _et al._ [2014] T. Elperin, N. Kleeorin, M. Liberman, and I. Rogachevskii, Turbulent diffusion of chemically reacting gaseous admixtures, Phys. Rev. E 90, 053001 (2014).
  * Elperin  _et al._ [2017] T. Elperin, N. Kleeorin, M. Liberman, A. N. Lipatnikov, I. Rogachevskii, and R. Yu, Turbulent diffusion of chemically reacting flows: Theory and numerical simulations, Phys. Rev. E 96, 053111 (2017).
  * Watanabe  _et al._ [2014] T. Watanabe, Y. Sakai, K. Nagata, and O. Terashima, Turbulent schmidt number and eddy diffusivity change with a chemical reaction, Journal of Fluid Mechanics 754, 98–121 (2014).
  * Shirian and Mani [2019] Y. Shirian and A. Mani, Eddy diffusivity in homogeneous isotropic turbulence (2019), arXiv:1905.08379 [physics.flu-dyn] .
  * Pouransari  _et al._ [2016] H. Pouransari, M. Mortazavi, and A. Mani, Parallel variable-density particle-laden turbulence simulation (2016), arXiv:1601.05448 .
  * Rosales and Meneveau [2005] C. Rosales and C. Meneveau, Linear forcing in numerical simulations of isotropic turbulence: Physical space implementations and convergence properties, Physics of Fluids 17, 095106 (2005).
  * Prend  _et al._ [2021] C. J. Prend, G. R. Flierl, K. M. Smith, and A. K. Kaminski, Parameterizing eddy transport of biogeochemical tracers, Geophysical Research Letters 48, e2021GL094405 (2021).
  * Wang  _et al._ [2015] H. Wang, M. Sun, Y. Yang, and N. Qin, A passive scalar-based method for numerical combustion, International Journal of Hydrogen Energy 40, 10658 (2015).