Physical interpretation of neural network-based nonlinear eddy viscosity
  models
Neural network-based turbulence modeling has gained significant success in
improving turbulence predictions by incorporating high--fidelity data. However,
the interpretability of the learned model is often not fully analyzed, which
has been one of the main criticism of neural network-based turbulence modeling.
Therefore, it is increasingly demanding to provide physical interpretation of
the trained model, which is of significant interest for guiding the development
of interpretable and unified turbulence models. The present work aims to
interpret the predictive improvement of turbulence flows based on the behavior
of the learned model, represented with tensor basis neural networks. The
ensemble Kalman method is used for model learning from sparse observation data
due to its ease of implementation and high training efficiency. Two cases,
i.e., flow over the S809 airfoil and flow in a square duct, are used to
demonstrate the physical interpretation of the ensemble-based turbulence
modeling. For the flow over the S809 airfoil, our results show that the
ensemble Kalman method learns an optimal linear eddy viscosity model, which
improves the prediction of the aerodynamic lift by reducing the eddy viscosity
in the upstream boundary layer and promoting the early onset of flow
separation. For the square duct case, the method provides a nonlinear eddy
viscosity model, which predicts well secondary flows by capturing the imbalance
of the Reynolds normal stresses. The flexibility of the ensemble-based method
is highlighted to capture characteristics of the flow separation and secondary
flow by adjusting the nonlinearity of the turbulence model.
Xin-Lei Zhang, Heng Xiao, Solkeun Jee, Guowei He
# Physical interpretation of neural network-based nonlinear eddy viscosity models

###### Abstract

Neural network-based turbulence modeling has gained significant success in improving turbulence predictions by incorporating high-fidelity data. However, the interpretability of the learned model is often not fully analyzed, which has been one of the main criticism of neural network-based turbulence modeling. Therefore, it is increasingly demanding to provide physical interpretation of the trained model, which is of significant interest for guiding the development of interpretable and unified turbulence models. The present work aims to interpret the predictive improvement of turbulence flows based on the behavior of the learned model, represented with tensor basis neural networks. The ensemble Kalman method is used for model learning from sparse observation data due to its ease of implementation and high training efficiency. Two cases, i.e., flow over the S809 airfoil and flow in a square duct, are used to demonstrate the physical interpretation of the ensemble-based turbulence modeling. For the flow over the S809 airfoil, our results show that the ensemble Kalman method learns an optimal linear eddy viscosity model, which improves the prediction of the aerodynamic lift by reducing the eddy viscosity in the upstream boundary layer and promoting the early onset of flow separation. For the square duct case, the method provides a nonlinear eddy viscosity model, which predicts well secondary flows by capturing the imbalance of the Reynolds normal stresses. The flexibility of the ensemble-based method is highlighted to capture characteristics of the flow separation and secondary flow by adjusting the nonlinearity of the turbulence model.

keywords: Machine learning, turbulence modeling, ensemble Kalman inversion, physical interpretability +
Footnote †: journal:

## 1 Introduction

Data-driven turbulence modeling has emerged as an important approach for predicting turbulent flows [1], which constructs functional mappings from mean velocity to the Reynolds stress by incorporating observa
tion data. Over the past few years, this paradigm of turbulence modeling has been pursued from various aspects, including choice of training data, development of training strategy, and representative form of the Reynolds stress. As for the training data, both the Reynolds stress and velocity data have been used for learning turbulence models. The velocity data become advocated for model learning as they are relatively straightforward to obtain in practical applications compared to the Reynolds stress data [2]. Regarding the training strategies, the conventional _a priori_ approach [3; 4; 5] trains neural network-based models without involving the RANS solver, which is pointed out [6] to have inconsistency issues in posterior tests. For this reason, the model-consistent training [7; 8; 9; 2; 10; 11] has been proposed to improve the predictive abilities of learned models by coupling the neural network and the RANS equation during the training process. Besides the two research lines mentioned above, the representative form of turbulence closure has also been investigated to empower the model with generalizability across different classes of flows. It is one critical step toward the ultimate goal of discovering unified turbulence models from data.

Various strategies have been proposed to represent the Reynolds stress, such as neural-network-based multiplicative correction [12], eigen perturbation method [13; 14; 15], symbolic expression [16], tensor basis neural network [4] and so on. Specifically, the neural-network-based multiplicative correction is introduced to modify turbulent production terms in turbulence transport equations, which can improve the velocity prediction of separated flows but is still under the Boussinesq assumption. The eigen perturbation method is proposed to present the Reynolds stress based on the eigen decomposition of the Reynolds stress tensor. The obtained eigenfunctions have physical interpretations to indicate the magnitude, shape, and orientation of the Reynolds stress tensor. This representation is a general form to represent the Reynolds stress but requires careful selection of the input features to ensure the Galilean invariance. To overcome these limitations, the nonlinear eddy viscosity model is often used as the base model, which is beyond the Boussinesq assumption and regards scalar invariants associated with velocity gradients as model inputs. Different techniques, including symbolic expression and neural networks, have been introduced to represent the Reynolds stress based on the nonlinear eddy viscosity model. In this work, we focus on the neural network-based representation, i.e., tensor basis neural network [4].

The tensor basis neural network is able to represent the anisotropy of the Reynolds stress flexibly due to its great expressive power. Neural networks have expressive power that increases exponentially with the depth of the network. Hence, it has the potential to achieve a universal or at least unified model to represent various flow characteristics. That is, one model form is applicable to multiple classes of flows, such as attached flows, separated flows, and corner flows, possibly with internal switching or branching. While such universality is not the objective of this work, it is appealing to have such possibilities in the future. However, the tensor basis neural network has intrinsic drawbacks due to the weak equilibrium assumption and the black-box feature. On the one hand, although the tensor basis neural network is the most general
nonlinear eddy viscosity model, it is still a local model under the weak equilibrium assumption. That is, the Reynolds stress anisotropy only depends on the local velocity gradient. To address this issue, the vector-cloud neural network [17] has been proposed to enforce the nonlocal dependence in the representative form. On the other hand, the trained neural network is still a black box and has encountered an interpretability crisis for neural network-based turbulence modeling. Therefore, it is of significant necessity to interpret the physical mechanism behind the learned neural network and guide the development of turbulence closures.

In this work, we aim to physically interpret the behavior of the learned turbulence model in terms of predictive improvement. Neural networks can represent complex functional relationships between physical quantities but have poor interpretability on the learned model behavior. In contrast, symbolic models are often assumed as interpretable since they can provide the causes and effects of the model behavior in the _a priori_ sense. It is noted that when the learned symbolic model provides a complicated expression that is highly composited or has many high-order terms, which would also be difficult to interpret. Some post-hoc approaches, such as the Shapley additive explanations (SHAP) method [18], have been proposed to interpret the black-box neural network models. These methods can indicate the importance value of each input feature on the neural network output [19] in the _a posteriori_ sense. However, they cannot provide physical insights into the mechanism of the learned model for improving the RANS prediction.

In this work, we investigate the physical interpretability of the learned turbulence model, represented with tensor basis neural networks [4]. The ensemble Kalman method is adopted to learn turbulence models from sparse observation data, including the lift force and velocity. We show that the behavior of the learned neural network is physically interpretable to improve flow predictions on two canonical flows, i.e., separated flow in the S809 airfoil and secondary flow in a square duct. The ensemble method can adjust the nonlinearity of the learned model to capture the different flow characteristics. Moreover, the capability of the ensemble Kalman method is shown in learning turbulence models from very sparse observation data. In addition, the normalization strategy is investigated to avoid feature clustering due to the stagnation point of airfoil flows. We note that the interpretability in this work refers to the model behavior of the trained neural network. It is different from the interpretability of neural networks in the machine learning community, which aims to present the features of neural networks in an understandable term, e.g., indicating the importance of input features with specific contribution values [18].

The rest of the paper is outlined as follows. The ensemble-based modeling methodology is elaborated in Section 2. The case setups and the training results are presented in Section 3 and 4, respectively. Finally, the paper is concluded in Section 5.
## 2 Methodology

For incompressible turbulent flows, the mean flow can be described by the RANS equation as

\[\begin{split}\nabla\cdot\mathbf{u}&=0\\ \mathbf{u}\cdot\nabla\mathbf{u}&=-\nabla p+\nu\nabla^{2}\mathbf{u }-\nabla\cdot\mathbf{\tau},\end{split} \tag{1}\]

where \(p\) is the mean pressure normalized by the flow density, \(\mathbf{u}\) is the velocity vector, \(\nu\) represents the molecular viscosity, and \(\mathbf{\tau}\) indicates the Reynolds stress1 to be modeled. Here we aim to construct neural-network-based turbulence models by incorporating available observations, such as lift force and velocity measurements. In the following, we introduce the Reynolds stress representation and the ensemble-based training method adopted in this work.

Footnote 1: Here we followed Pope’s convention [20] of defining Reynolds stress as the covariance of the velocity fluctuations i.e., \(\tau_{ij}=\left\langle u^{\prime}_{i}u^{\prime}_{j}\right\rangle\). We note that in the literature (e.g., [21]) it is more common to call \(-\left\langle u^{\prime}_{i}u^{\prime}_{j}\right\rangle\) the Reynolds stress because of its role in the RANS momentum equations.

### Neural-network-based turbulence closure

The tensor basis neural network [4] is used to represent the Reynolds stress due to the flexibility to represent the anisotropy of Reynolds stress. In the tensor basis neural network, the Reynolds stress \(\mathbf{\tau}\) is decomposed into a deviatoric part and an isotropic part, as

\[\begin{split}\mathbf{\tau}&=2k\sum_{\ell=1}^{10}g^{( \ell)}\mathbf{T}^{(\ell)}+\frac{2k}{3}\mathbf{I},\\ \text{with}\quad g^{(\ell)}&=g^{(\ell)}\left(\theta _{1},\dots,\theta_{5}\right),\end{split} \tag{2}\]

where \(k\) is the turbulent kinetic energy, \(\mathbf{T}\) is the tensor basis, \(g^{(\ell)}\) is the coefficient of the tensor basis to be determined, \(\mathbf{\theta}\) is the scalar invariants, and \(\mathbf{I}\) is the identity matrix. The \(g\) functions are represented with neural networks in this work which approximates functional mappings from the scalar invariants \(\mathbf{\theta}\) to the basis coefficients. There are ten independent tensor bases based on the Cayley-Hamilton theory [22] and five scalar invariants for incompressible flows. In the 2D scenario, only two scalar invariants and three tensor bases are remained [20]. Further, the third tensor basis can be incorporated in the pressure term for incompressible flows, leaving only two scalar invariants. The first four tensor bases can be written as

\[\begin{split}\mathbf{T}^{(1)}&=\hat{\mathbf{S}}, \qquad\mathbf{T}^{(2)}=\hat{\mathbf{S}}\hat{\mathbf{W}}-\hat{\mathbf{W}}\hat{ \mathbf{S}},\\ \mathbf{T}^{(3)}&=\hat{\mathbf{S}}^{2}-\frac{1}{3} \{\hat{\mathbf{S}}^{2}\}\mathbf{I},\qquad\mathbf{T}^{(4)}=\hat{\mathbf{W}}^{2 }-\frac{1}{3}\{\hat{\mathbf{W}}^{2}\}\mathbf{I}.\end{split} \tag{3}\]
In the formula above, \(\{\cdot\}\) denotes the trace operator, and the \(\hat{\mathbf{S}}\) and \(\hat{\mathbf{W}}\) are the normalized strain rate and the rotation rate based on the turbulence time scale \(\tau_{s}\), i.e.,

\[\begin{split}\hat{\mathbf{S}}&=\tau_{s}\mathbf{S} \quad\hat{\mathbf{W}}=\tau_{s}\mathbf{W}\\ \text{with}&\mathbf{S}&=\frac{1}{2}( \nabla\mathbf{u}+\nabla\mathbf{u}^{\top})\quad\text{and}\quad\mathbf{W}=\frac{1}{2}( \nabla\mathbf{u}-\nabla\mathbf{u}^{\top}).\end{split} \tag{4}\]

The time scale \(\tau_{s}\) can be estimated with the turbulent kinetic energy \(k\) and the dissipation rate \(\varepsilon\) or the specific dissipation rate \(\omega\). It is noted that the time scale becomes zero as we approach the wall. Hence one can bound the time scale with the Kolmogorov scale [23] as

\[\tau_{s}=\max\left(\frac{k}{\varepsilon},C_{\tau}\sqrt{\frac{\nu}{\varepsilon} }\right), \tag{5}\]

where \(C_{\tau}\) is constant and set as 6 in this work.

### Normalization of input features

The input features of the neural networks should be scaled within \([-1,1]\) to accelerate the training convergence. The min-max normalization is able to confine the input features within the range of \([0,1]\) (see e.g., Ref. [2]). A normalized feature \(\hat{\theta}\) can be formulated as \(\hat{\theta}=(\theta-\theta_{\min})/(\theta_{\max}-\theta_{\min})\), where the subscript'min' and'max' indicate the minimum and maximum value of a given feature \(\theta\). However, when there exist singular points with extremely large magnitudes in computational domains, this normalization strategy can lead to severe feature clustering. For instance, the velocity gradient near a stagnation point can have an extremely large value. Using the global maximum value to normalize entire input features will lead to most feature values clustering around 0, which would significantly affect the training performance.

In this work, the scalar invariants \(\hat{\mathbf{\theta}}\) are normalized with the local time scale \(\tau_{s}\) [e.g., 3; 15; 24] based on

\[\begin{split}\hat{\theta}_{1}&=\{\tilde{\mathbf{S} }^{2}\},\qquad\hat{\theta}_{2}=\{\tilde{\mathbf{W}}^{2}\},\\ \tilde{\mathbf{S}}&=\frac{\mathbf{S}}{\|\mathbf{S} \|+1/\tau_{s}},\quad\text{and}\qquad\tilde{\mathbf{W}}=\frac{\mathbf{W}}{\| \mathbf{W}\|+1/\tau_{s}}.\end{split} \tag{6}\]

With this specific normalization, the scalar invariants can be scaled within \([-1,1]\) to avoid feature clustering along certain directions. The normalized scalar invariants \(\mathbf{\theta}\) are used as the neural network inputs, and the coefficients \(g\) of the tensor bases are regarded as the outputs. Further, the neural network outputs \(g\) are combined with the tensor bases \(\mathbf{T}\) to form the anisotropic part of the Reynolds stress. The obtained Reynolds stress is used to predict the velocity and pressure fields by solving the RANS equations. Moreover, the constructed Reynolds stress \(\mathbf{\tau}\) is used to compute the turbulence production term in the turbulent kinetic energy and dissipation rate transport equations. Further, the neural network weights are optimized by incorporating observation data based on the ensemble Kalman method, which will be illustrated in the following subsection.
### Model-consistent training with ensemble Kalman method

Model-consistent training couples a neural network and a CFD solver during the training process. By doing this, it can ensure consistency between the training and prediction environments, thereby alleviating the ill-conditioning of the RANS model operator [25]. Moreover, this strategy can leverage sparse observation data, e.g., velocity measurements, to train the neural network-based model. This is in contrast to the prior training where the model is often trained with the full field data of the Reynolds stress and has poor generalizability due to the inconsistency issue [10]. The model-consistent training amounts to finding the optimal weights of neural networks that lead to the best fit with the sparse observation data.

Various training methods can be used to perform the model-consistent training, including the adjoint method [7], the ensemble method [10], and the genetic programming method [8]. We use the ensemble Kalman method for model training due to its non-derivative nature and good training efficiency. The ensemble method is a statistical inference method that uses an ensemble of samples to guide the optimization [26], which has been used for the physical modeling of subsurface flows [27] and turbulent flows with high Reynolds numbers [11; 24]. We use this method to train the turbulence model represented with the tensor basis neural network. The update scheme of the ensemble Kalman method can be formulated as

\[\begin{split}&\mathbf{w}_{j}^{i+1}=\mathbf{w}_{j}^{i}+\mathsf{K}(\mathsf{ y}_{j}-\mathsf{H}\mathbf{w}_{j}^{i})\\ &\text{with}\quad\mathsf{K}=\mathsf{PH}^{\top}(\mathsf{H}\mathsf{ PH}^{\top}+\mathsf{R})^{-1}.\end{split} \tag{7}\]

Herein \(\mathsf{H}\) is the local gradient of the model prediction \(\mathcal{H}[\mathbf{w}]\) with respect to the weights of neural networks \(\mathbf{w}\), \(\mathsf{P}\) is the model error covariance, \(\mathsf{R}\) is the observation error covariance, \(\mathsf{y}\) is the observation data, and \(i\) and \(j\) represent the index of optimization iteration and sample, respectively. The model operator \(\mathsf{H}\) is often avoided to be computed by reformulating the Kalman gain matrix as

\[\mathsf{K}=\mathsf{S}_{w}\mathsf{S}_{y}^{\top}(\mathsf{S}_{y}\mathsf{S}_{y}^{ \top}+\mathsf{R})^{-1}.\]

The square-root matrices \(\mathsf{S}_{w}\) and \(\mathsf{S}_{y}\) are defined as

\[\mathsf{S}_{w}^{i} =\frac{1}{\sqrt{N_{e}-1}}\left[\mathbf{w}_{1}^{i}-\overline{\mathbf{w}}^ {i},\mathbf{w}_{2}^{i}-\overline{\mathbf{w}}^{i},\cdots,\mathbf{w}_{N_{e}}^{i}-\overline{ \mathbf{w}}^{i}\right], \tag{8a}\] \[\mathsf{S}_{y}^{i} =\frac{1}{\sqrt{N_{e}-1}}\left[\mathcal{H}[\mathbf{w}_{1}^{i}]- \mathcal{H}[\overline{\mathbf{w}}^{i}],\mathcal{H}[\mathbf{w}_{2}^{i}]-\mathcal{H}[ \overline{\mathbf{w}}^{i}],\cdots,\mathcal{H}[\mathbf{w}_{N_{e}}^{i}]-\mathcal{H}[ \overline{\mathbf{w}}^{i}]\right],\] (8b) \[\overline{\mathbf{w}}^{i} =\frac{1}{N_{e}}\sum_{j=1}^{N_{e}}\mathbf{w}_{j}^{i}, \tag{8c}\]

which are estimated from the samples at every iteration. In this work, we use the ensemble-based Kalman update scheme for learning turbulence models in a model-consistent manner. As pointed out in Ref. [10], in scenarios having large data sets, e.g., time-dependent three-dimensional flow fields, the present algorithm would be computationally expensive and need to incorporate reduced-order techniques such as the truncated singular value decomposition [28].
Note that the ensemble Kalman method can train the neural network-based model with multiple observation data, including the measurements at various flow conditions. Specifically, we can incorporate the observation data at different flow conditions sequentially. It is achieved by training neural networks with each observation data in several inner loops. The maximum iteration number of the inner loop is set as 3 in this work based on our sensitivity study. Moreover, the observation data are shuffled randomly before training, which allows escaping from local minima similar to the stochastic gradient descent method [29]. Further, the Kalman update scheme is used to incorporate the observation data in the shuffled order till the entire data sets are traversed. After that, the observation data will be reshuffled and continue to be incorporated with the ensemble Kalman method. The practical implementation of the ensemble method is presented in A. One can also augment the observation with data from different flow conditions. However, this may drop into local minima and lead to unsatisfactory predictive accuracy in certain cases since the ensemble method aims to reduce the L2 norm of the total data misfit. By shuffling the training data, the method is able to find the global minimum and provide more accurate turbulence models based on our numerical tests.

## 3 Case setup

We use two cases to demonstrate the physical interpretation of the ensemble-based turbulence modeling, i.e., the flow over the S809 airfoil and the flow in a square duct. The two cases represent canonical separated flows and secondary flows, respectively. Both are challenging for conventional linear eddy viscosity models. The distinct flow characteristics are able to examine the flexibility of the ensemble-based method in learning interpretable models from partial observation. The details of the case setup are described in the following subsections.

### Flow over S809 airfoil

Flow over the S809 airfoil has been widely used for numerical validation of turbulence models as well as their data-driven counterparts [12; 7]. Such flows are challenging for linear eddy viscosity models at large angles of attack due to the flow separation. Conventional RANS models cannot accurately predict the massive flow separation, which further leads to the overestimation of the lift force beyond the stall angle [30]. Here we aim to interpret the behavior of the neural network-based model learned from lift force measurements with the ensemble Kalman method.

The Reynolds number is \(Re_{c}=2\times 10^{6}\) based on the inflow velocity and chord length. The angle of attack \(\alpha\) varies from 1\({}^{\circ}\)-18\({}^{\circ}\). At large angles of attack, conventional turbulence models underestimate the separation zones [30], which leads to large discrepancies in the predictions of the lift force [30]. The unstructured mesh with around 78000 cells is used to discretize the computational domain. The mesh grid
from the work [7] is adopted in the work as shown in Fig. 1. The no-slip condition is employed on the airfoil surface. The height of the first cell in the normal direction corresponds to \(y^{+}\approx 1\).

For the flow around the S809 airfoil, the available measurement data is the lift force, which is the integral type data source. Such limited observation would increase the ill-posedness of the inverse problems [31; 32; 33]. That is, different model functions can provide similar lift forces. Further, the learned model could have poor predictive accuracy and robustness due to the ill-posedness issue. To alleviate the issue, we use observations at two angles of attack, i.e., \(8.2^{\circ}\) and \(14.24^{\circ}\) to train the model. The former is attached flows, while the latter is separated flows. Learning from both the attached and separated flows can provide a model with better predictive ability.

As for the setup of the ensemble-based learning algorithm, the number of samples is taken as 50. The initial relative variance of the samples is set as 0.1, which is used to draw the random samples. The measurements of the lift force [34] are used as training data. The relative observation error is set as 0.01. The first two scalar invariants \(\theta_{1}\) and \(\theta_{2}\), and the first two tensor coefficients \(g^{(1)}\) and \(g^{(2)}\) are used as the inputs and the outputs of the neural network, respectively. We use the \(k\)-\(\omega\) model [21] as the baseline model, which is suitable for complex boundary layer flows with adverse pressure gradient compared to standard \(k\)-\(\varepsilon\) model [35].

### Flow in a square duct

The secondary flow in a square duct is mainly driven by the imbalance of the Reynolds normal stresses \(\tau_{yy}-\tau_{zz}\)[36]. The linear eddy viscosity model is not able to predict the secondary flow since it cannot well estimate the anisotropy of the Reynolds stress. We use this case to demonstrate the flexibility of the ensemble

Figure 1: Computational domain and mesh grid for computations of flows over the S809 airfoil
method in building interpretable nonlinear models from sparse observation data of secondary flows.

The Reynolds number based on the bulk velocity and half of the duct height is \(Re_{h}=3500\) for this case. Only one quadrant of the physical domain is simulated, considering the symmetry of the flow to the centerlines along \(y-\) and \(z-\)axes. The mesh with \(50\times 50\) is used to discretize the domain. The non-slip condition is imposed on the wall, and the symmetry condition is imposed at the symmetry boundary.

As for the setup of the ensemble-based learning algorithm, the number of samples is set as 50. The initial variance of the weights is set as 0.1. The DNS data [37] are used to train the neural network-based model. The velocity profiles at \(y/H=0.25,0.5,0.75,1.0\) are regarded as the observation data. The total number of observation data points is 200. The relative observation error is set as 0.01. For this case, we use the first two scalar invariants \(\theta_{1}\) and \(\theta_{2}\) as the input features, and the first four tensor coefficients \(g^{(1-4)}\) as the model outputs. Compared to the linear eddy viscosity model, the tensor bases of \(\mathbf{T}^{(2)}\), \(\mathbf{T}^{(3)}\) and \(\mathbf{T}^{(4)}\) are introduced to capture the anisotropy of the Reynolds normal stress, thereby producing the secondary flow. The \(k\)-\(\varepsilon\) model [35] is used as the baseline model in this case, since this model is often taken as the base of nonlinear eddy viscosity model for secondary flows [38; 39].

In this work, the open source CFD library OpenFOAM [40] is used to solve the RANS equations with turbulence models. Specifically, the built-in solver _simpleFOAM_ is used to solve the RANS equations, given the Reynolds stress fields. The Reynolds stresses are constructed with the neural networks, and the scalar invariants from the RANS computation are taken as inputs of the networks. The weights of the neural network are updated with the observation data based on the ensemble Kalman method. The TensorFlow [41] library is used to construct the neural network, and the DAFI code [42] is used to implement the ensemble Kalman method. The test cases and the weights of the learned neural network are publicly available [43] for reproducibility.

Figure 2: Computational domain and mesh grid for the fully–developed square duct case.
## 4 Results

### Flows over S809 airfoil

#### 4.1.1 Training performance

The learned model can improve the predictions of the lift force compared to the baseline \(k\)-\(\omega\) model. The predicted lift with the learned and baseline models at two chosen angles of attack is listed in Table 2. The baseline model predicts the aerodynamic lift \(C_{l}=1.25\) at the angle of attack \(14.24^{\circ}\), which deviates much from the experimental observation \(C_{l}=1.05\)[34] due to the massive flow separation. In contrast, the learned model provides \(C_{l}=1.07\), which is in good agreement with the experimental data. At \(\alpha=8.2^{\circ}\), the boundary layer is attached, and the baseline model provides good prediction with \(C_{l}=0.97\). The learned model predicts \(C_{l}=1.0\), which is slightly deviated from the observation \(C_{l}=0.95\). That is because the learning method decreases the data misfit at the two flow conditions simultaneously. The significant decrease of lift force at the angle \(\alpha=14.24^{\circ}\) is achieved with the sacrifice of a slight discrepancy at \(\alpha=8.2^{\circ}\). In general, the learned model can provide good predictions close to the experimental measurements, which is not surprising since the experimental data are used to train the model function.

The prediction of the pressure coefficient \(C_{p}\) is improved with the learned model compared to the baseline model. Figure 3 shows the predicted \(C_{p}\) with comparison among the learned model, the baseline \(k\)-\(\omega\) model, and the experimental data. It can be seen that the baseline \(k\)-\(\omega\) model can predict well the pressure

\begin{table}
\begin{tabular}{c|c|c} \hline Cases & S809 airfoil & Square duct \\ \hline mesh counts & \(\approx 78000\) & 2500 \\ Reynolds number & \(Re_{c}=2\times 10^{6}\) & \(Re_{h}=3500\) \\ data & \(C_{l}\) (at \(\alpha=8.2^{\circ}\) and \(14.24^{\circ}\)) & \(\mathbf{u}\) (at \(y/H=0.25,0.5,0.75,1\)) \\ baseline model & \(k\)–\(\omega\) & \(k\)–\(\varepsilon\) \\ initial relative variance & 0.1 & 0.1 \\ relative observation error & 0.01 & 0.01 \\ sample size & 50 & 50 \\ \hline \end{tabular}
\end{table}
Table 1: Computational parameters used in the flow around the S809 airfoil and the flow in a square duct

\begin{table}
\begin{tabular}{c|c|c|c|c} \hline  & \(\alpha\) & baseline \(k\)–\(\omega\) & learned model & experiment [34] \\ \hline \(C_{l}\) & \(14.24^{\circ}\) & 1.25 & 1.07 & 1.05 \\ \(C_{l}\) & \(8.2^{\circ}\) & 0.97 & 1.00 & 0.95 \\ \hline \end{tabular}
\end{table}
Table 2: Summary of the prediction in the aerodynamic lift \(C_{l}\) with the learned and baseline \(k\)–\(\omega\) models compared to the experimental data for the S809 airfoil
distribution on the surface of the S809 airfoil at \(\alpha=8.2^{\circ}\). However, at \(\alpha=14.24^{\circ}\) the baseline model underestimates the suction pressure on the upper surface of the airfoil. In contrast, the learned model with the ensemble method is able to predict \(C_{p}\) in better agreement with the experimental data for both angles of attack. The results demonstrate that the learned model can leverage integral data, i.e., lift force \(C_{l}\), to improve the prediction of wall pressure distribution.

#### 4.1.2 Physical interpretation of the model behavior

The learned model can accurately predict the aerodynamic lift and the wall pressure distribution beyond the stall angle compared to the baseline model. Such improvements can be interpreted by analyzing the model behavior in terms of the friction coefficient, flow separation, and modeled quantities. Therefore, we further provide comparisons between the learned model and the baseline model in the following.

The friction coefficient \(C_{f}\) on the airfoil is investigated to interpret the reduction of \(C_{l}\) with the learned model at \(\alpha=14.24^{\circ}\). Figure 4 shows \(C_{f}\) from the baseline and the learned models. The learned model leads to the early onset of flow separation, while the baseline model delays the separation. The early separation would lead to an enlarged re-circulation region. Therefore, the early occurrence of the flow separation is responsible for the modification of \(C_{p}\) on the suction side as shown in Fig. 3. Also, it is observed that the friction coefficient is reduced on the pressure side. That is because the flow around the airfoil changes globally with the learned model due to the upstream separation point on the upper surface. This global modification of the flow also affects the flow on the lower surface.

The enlarged re-circulation region with the learned model can be clearly seen in Figure 5, which presents streamlines around the airfoil. The baseline \(k\)-\(\omega\) model predicts a relatively small separation bubble compared to that with the learned model, which is consistent with the lift prediction listed in Table 2. The learned

Figure 3: Wall pressure coefficient \(C_{p}\) at \(\alpha=8.2^{\circ}\) and \(14.24^{\circ}\) with the learned model and the baseline \(k\)–\(\omega\) model compared to the experimental data [34] for the S809 airfoil
model produces a sufficient massive separation region, which leads to the improvement of \(C_{l}\).

We investigate the model function \(g\) to interpret the reason for the early onset of the flow separation. The learned model for the S809 airfoil is almost a linear eddy viscosity model. It can be seen from Figure 6 which presents the learned \(g\) functions on \(\theta_{1}\) at fixed planes of \(\theta_{2}/\theta_{\max}=0.25\) and \(0.75\). The magnitude of the learned \(g^{(1)}\) function decreases to around \(0.05\) and \(0.075\) at the plane of \(\theta_{2}/\theta_{\max}=0.25\) and \(0.75\), respectively, while the magnitude for the baseline model is constant at \(0.09\). In contrast, the learned \(g^{(2)}\) function is almost zero with the order of magnitude of \(10^{-6}\), which is similar to the baseline model, i.e., \(g^{(2)}=0\). Therefore, the learned model can be considered a linear eddy viscosity model, which is capable of capturing the flow separation on the S809 airfoil.

Recall that the Reynolds stress anisotropy is the linear combination of the learned function \(g(\mathbf{\theta})\) and the tensor bases \(\mathbf{T}\), i.e., \(\mathbf{b}=\sum g^{(\ell)}\mathbf{T}^{(\ell)}\). Fig. 7 presents the contour plots of the magnitude of each tensor component. The first term is a linear term on the strain rate (see Eqs. (2) and (3)). For the airfoil case, the nonlinear term is almost zero as shown in Fig. 7(b), which further confirms that the learned model can be considered a linear model under the Boussinesq assumption. Such linear models can also achieve good predictions in the lift force for the 2D airfoil case. This is consistent with the work of Singh et al. [30; 12],

Figure 4: The comparison of the friction coefficient between the learned model and the baseline model at the angle of attack of \(14.24^{\circ}\). The round circles in the right panel indicate the separation locations.

Figure 5: Predicted separation bubbles at angles of attack \(14.24^{\circ}\) with the learned model and the \(k\)–\(\omega\) model for the S809 airfoil case.
where a multiplicative correction is added in the turbulence transport equation to modify the eddy viscosity and improve the lift prediction beyond the stall angle. The linear eddy viscosity assumption is sufficient in this study of the S809 airfoil.

Since the learned model can be regarded as a linear eddy viscosity model, we further investigate the effects of the learned eddy viscosity on the model prediction. Figure 8 shows the predicted eddy viscosity with the learned model compared to the baseline \(k\)-\(\omega\) model. The eddy viscosity is computed based on the \(g^{(1)}\) function, which can be formulated as

\[\nu_{t}=-\frac{g^{(1)}k}{C_{\mu}\omega}.\]

The model constant is \(C_{\mu}=0.09\). It can be seen clearly that the eddy viscosity is reduced, particularly around the upstream boundary layer and the separated region. The eddy viscosity can transfer the energy

Figure 6: Plots of the learned mapping between the scalar invariants \(\mathbf{\theta}\) and the tensor coefficient \(\mathbf{g}\), compared to the baseline for the S809 airfoil case. For the learned model, the plots indicate the learned function at \(\theta_{2}/\theta_{\max}=0.25\) and \(0.75\).

Figure 7: Learned tensor components for the S809 airfoil case. The arrow indicates the direction of incoming flow (\(14.24^{\circ}\)).
from the outer flow to the boundary layer, which is able to restrain momentum reduction and further flow separation. The reduced eddy viscosity would weaken the energy transfer from the outer flow and enable the boundary layer to be less resistant to the adverse pressure gradient. As such, the reduced eddy viscosity in the upstream can induce a relatively early onset of the flow separation and further a large separation region. We note that the reduced eddy viscosity at the upstream boundary layer is also observed [44] from the model learned with the adjoint-based method [12], which further confirms the physical interpolation for the improved model prediction. In addition, the eddy viscosity at the lower surface is reduced as well compared to the baseline \(k\)-\(\omega\) model. This explains the friction coefficient reduction on the pressure side as presented in Fig. 4, because the wall shear stress (or the velocity gradient) is sensitive to a small change near the wall - here the eddy viscosity is slightly reduced.

The physical interpretation can empower the learned model with good predictive ability. Figure 9 presents the predicted lift force at different angles of attack from \(1^{\circ}\) to \(18^{\circ}\). The results show that the learned model can be well generalized to other flow conditions at different angles of attack. It can be seen that the baseline model has significant discrepancies in the lift coefficients for angles of attacks larger than around \(7.5^{\circ}\). In contrast, the learned model improves the prediction on \(C_{l}\) across the angle \(\alpha\). The reason for the improved prediction in other angles is due to the appropriate estimation of flow separation, as shown at the top of Figure 9. Specifically, at the small angle of attack, e.g., \(\alpha=1^{\circ}\), the baseline and learned models produce similar attached flow around the airfoil, and hence both predict the lift force in good agreement with the experiment. However, at the large angle of attack \(\alpha=11^{\circ}\), the baseline model still expects attached flow, which leads to the overestimation in \(C_{l}\) compared to the experimental data. In contrast, the learned model captures the flow separation near the trailing edge, which reduces the lift force and provides a good agreement with the lift force data. Additionally, at \(\alpha=18^{\circ}\), the baseline model predicts the flow separation but still underestimates the separation bubble size, which leads to the lift force being larger than the experimental measurement. The learned model leads to the early onset of the flow separation and provides

Figure 8: Comparison of the eddy viscosity between the learned and baseline \(k\)–\(\omega\) models for the S809 airfoil. Note that the eddy viscosity is normalized by the molecular viscosity as \(\nu_{t}/\nu\). The arrow indicates the direction of incoming flow (\(14.24^{\circ}\)).
a larger separation bubble compared to the baseline model, thereby improving the lift force prediction. The predictive performance of the learned model at additional angles \(\alpha=11^{\circ}\) and \(18^{\circ}\) can be found in Appendix B. Further generalization tests with different geometries are beyond the scope of the present work and will be conducted in the near future. We note that the adjoint-based learning method has been used in the S809 airfoil, demonstrating that the learned model can be generalized well for different configurations such as the S805 and S814 airfoils [12]. Here we use the ensemble Kalman method which is comparable to the adjoint method in model learning as demonstrated in Ref. [10]. Hence the current learned model could be generalized to other cases as the adjoint method does.

The S809 airfoil case has been used in various works [12; 7; 44; 11], including Singh et al. (2017) [12], where the learned model suppresses the turbulent production at the upstream boundary layer, leading to early flow separation. The difference between the present work and previous studies mainly lies in two aspects. First, the current modeling framework and the training method are different from the previous works. Specifically, the nonlinear eddy viscosity model is used in this study since it is flexible to capture both separated and secondary flows, while previous studies including Singh et al. (2017) [12] use a linear eddy viscosity model which is not able to predict secondary flows. Moreover, the ensemble Kalman method is used to train the neural networks in this work, while previous studies [12; 7; 44] often use the adjoint-based method for model inference. Second, the case in this work highlights the consistency of the learned model behaviors independent of the training techniques, in addition to the discussion on the capability of the data-driven method in prediction improvement as was only done in previous studies. That is, different

Figure 9: Tests on various angles of attack with comparison among the baseline model, the learned model, and the experiment [34] for the S809 airfoil case. The training cases are also indicated in the plot.
model representations and training methods lead to similar predictive improvements and model behaviors.

### Flow in a square duct

#### 4.2.1 Training performance

The ensemble Kalman method can learn a neural network-based turbulence model with improved velocity prediction for the square duct case. It can be seen from Figure 10 where the vector and contour plots of velocities are presented. The vector plots are presented in the first column of Figure 10, where the isolines indicate the levels of \(u_{y}=0.5,1.0\) and \(1.2\). It shows that the baseline model cannot predict the in-plane secondary flow, while the ensemble-based learned method can estimate the in-plane velocity vectors in a similar pattern as the DNS data. The contour plots of \(u_{x}\) and \(u_{y}\) are presented in the last two columns of Figure 10. The plots of \(u_{z}\) are omitted for brevity since it is symmetric to the vertical velocity \(u_{y}\). The axial velocity \(u_{x}\) with the baseline model and the learned model both have good agreement with the DNS data. The vertical velocity \(u_{y}\) is not captured at all by the baseline model, while the learned model can capture similar patterns as the DNS results.

Figure 10: Velocity \(u_{x}\) and \(u_{y}\) predicted from the learned models (center row) and baseline model (bottom row), compared with the ground truth (top row), for the square duct case. The velocity vectors are plotted along with contours of the streamwise velocity \(u_{x}\).
The learned model improves the prediction of the in-plane velocity by capturing the Reynolds stress imbalance and Reynolds shear stress. It is supported by Figure 11 where the Reynolds stress components and the imbalance of Reynolds normal stress are presented. The in-plane velocity is driven by the Reynolds stress imbalance \(\tau_{yy}-\tau_{zz}\) and the Reynolds shear stress \(\tau_{yz}\) based on the axial vorticity transport equation [36]:

\[u_{y}\frac{\partial\omega_{x}}{\partial y}+u_{z}\frac{\partial\omega_{x}}{ \partial z}-\nu\nabla^{2}\omega_{x}+\frac{\partial^{2}}{\partial y\partial z} \left(\mathbf{\tau}_{zz}-\mathbf{\tau}_{yy}\right)+\left(\frac{\partial^{2}}{\partial y \partial y}-\frac{\partial^{2}}{\partial z\partial z}\right)\mathbf{\tau}_{yz}=0. \tag{9}\]

For this reason, capturing the in-plane velocity requires well estimating Reynolds normal stress imbalance \(\tau_{yy}-\tau_{zz}\) and Reynolds shear stress \(\tau_{yz}\). From Figure 11, the learned model shows significant improvements in the prediction of \(\tau_{yz}\) and \(\tau_{yy}-\tau_{zz}\) compared to the baseline. Although the model still has discrepancies with the DNS data near the duct center, the noticeable improvement in \(\tau_{yz}\) and \(\tau_{yy}-\tau_{zz}\) allow us to obtain good agreement to the DNS data in the in-plane velocity prediction. Specifically, the baseline model estimates almost zero for both the Reynolds shear stress \(\tau_{yz}\) and the imbalance of Reynolds normal stresses \(\tau_{yy}-\tau_{zz}\) in the entire computational domain. In contrast, the learned model can predict them similarly obtained in the DNS, which significantly improves the in-plane velocity as shown in Fig. 10. For the Reynolds normal stresses \(\tau_{xx}\) and \(\tau_{yy}\) and the Reynolds shear stress \(\tau_{xy}\), the baseline model and the learned model give similar predictions since the in-plane velocity can not guide the training in these Reynolds stress components.

The profiles of velocity and the Reynolds stress at \(y/h=0.25,0.5,0.75,1\) are provided in Fig. 12. It can be seen that the streamwise velocity \(u_{x}\) is similar between the baseline model and the learned model, and both can have good agreement with the DNS data. As for the in-plane velocity \(u_{y}\), the baseline model is not able to predict the in-plane velocity and provide \(u_{y}=0\) at the entire domain. In contrast, the learned model significantly improves the prediction of \(u_{y}\) in better agreement with the DNS data. The plots of the Reynolds stress show that the learned model provides better predictions in the imbalance of the Reynolds normal stress \(\tau_{yy}-\tau_{zz}\) than the baseline model. As for the Reynolds shear stress \(\tau_{yz}\), both the learned and baseline models have noticeable discrepancies from the DNS data. Also, it is observed that the learned model has larger discrepancies near the diagonal line of the computational domain compared to the baseline model, which is consistent with the plots in Fig. 11. Additional results of the velocity and the Reynolds stresses at \(y/h=0.2,0.4,0.6,0.8\) are presented in B.

#### 4.2.2 Physical interpretation of model behavior

The behavior of the learned model can be interpreted based on the learned tensor coefficients \(g\). In the secondary flow, the axial velocity \(u_{x}\) is orders of magnitude larger than the in-plane velocity. Also, only four Reynolds stress components, i.e., Reynolds shear stress \(\tau_{xy}\), \(\tau_{xz}\), \(\tau_{yz}\), and Reynolds normal stress imbalance \(\tau_{yy}-\tau_{zz}\), affect the velocity [45]. The former two components affect the axial velocity, and the latter components of \(\tau_{yz}\) and \(\tau_{yy}-\tau_{zz}\) affect in-plane velocity. It can be further derived [45] that only the
coefficient \(g^{(1)}\) and the combination \(g^{(2)}-0.5g^{(3)}+0.5g^{(4)}\) can be learned with velocity data in the scenario of only first four tensor bases. Moreover, there is only one independent scalar invariant since \(\theta_{1}\approx-\theta_{2}\)[2]. Therefore, we investigate the functional mapping from the scalar invariant \(\theta_{1}\) to the coefficient \(g^{(1)}\) and the combination \(g^{(2)}-0.5g^{(3)}+0.5g^{(4)}\).

The coefficient \(g^{(1)}\) and the combination of \(g^{(2-4)}\) are shown in Figure 13. The learned function of \(g^{(1)}\) can be seen from Figure 13(a). Note that the coefficient \(g^{(1)}\) is equivalent to the \(-C_{\mu}\) of the \(k\)-\(\varepsilon\) model. The difference lies in that the coefficient has dependencies on local scalar invariants in this work rather than a constant, i.e., \(-0.09\). The \(g^{(1)}\) function with the learned model varies slightly from \(-0.87\) to \(-0.78\). For the small scalar invariants that are located around the duct center, the magnitude of the \(g^{(1)}\) function is less than \(0.08\). As the scalar invariant increases, the magnitude increases to around \(0.087\), which is slightly less than the baseline value (i.e., \(0.09\)). The baseline model provides the combination of \(g^{(2-4)}\) at almost zero, which cannot capture the in-plane velocity. In contrast, the learned model increases the magnitude of the

Figure 11: Reynolds normal stresses \(\tau_{yy}\) and \(\tau_{zz}\), Reynolds shear stresses \(\tau_{yz}\), and imbalance of Reynolds normal stresses \(\tau_{yy}-\tau_{zz}\) predicted from the learned model (center row) and the baseline model (bottom row), compared with the ground truth DNS (top row), for the square duct case.
combination at the range of nearly \([0.0025,0.01]\). This leads to nonlinear functional mappings between the Reynolds stress and the strain rate. Such nonlinear models capture the Reynolds shear stress \(\tau_{yz}\) and the Reynolds normal stress imbalance \(\tau_{yy}-\tau_{zz}\), which further improve the prediction of the in-plane velocity.

The ensemble-based model-consistent training is flexible to provide interpretable models based on sparse observations. It can be supported by the results of the tensor components as shown in Figure 14. It shows that the linear tensor component \(g^{(1)}\mathbf{T}^{(1)}\) from the learned model is larger than other nonlinear components, i.e., \(g^{(2)}\mathbf{T}^{(2)}\), \(g^{(3)}\mathbf{T}^{(3)}\), and \(g^{(4)}\mathbf{T}^{(4)}\), but at similar magnitudes. This is in contrast to the S809 airfoil case, where the linear tensor component is larger than the nonlinear components by several orders of magnitude, as shown in Fig. 7. The relatively large magnitude of the nonlinear tensors in this case is due to the secondary

Figure 12: Prediction of velocity and Reynolds stress along profiles at \(y/H=0.25,0.5,0.75,1\) with comparison among the learned model, the baseline model, and the experimental data, for the square duct case
flow characteristics that are driven by the imbalance of the Reynolds normal stress. The linear tensor \(g^{(1)}\mathbf{T}^{(1)}\) cannot capture the anisotropy of the Reynolds stress, and the nonlinear components play dominant roles in predicting in-plane velocities. Hence, for the square duct case, the ensemble-based training leads to a nonlinear model with considerable magnitude for the nonlinear terms.

In general, the ensemble-based method provides an interpretable turbulence model with appropriate nonlinearity according to limited observation data. For instance, in the scenario of separated flows over airfoils, the optimization of linear eddy viscosity is able to remedy the deficiency on the adverse pressure gradient as shown in Figure 7. In contrast, for the secondary flow, the nonlinear terms are required to accurately estimate the imbalance of the Reynolds normal stress as shown in Figure 14, which is the driving force for the spanwise vorticity. We emphasize that the available observations are often sparse in practical applications, e.g., lift force and sparse velocity measurements, as used in this work. Such severe ill-posedness poses challenges to the training method in learning dominant physical mechanisms with various flow characteristics. Hence, the flexibility of ensemble-based training is demonstrated in discovering interpretable models from sparse data.

The square duct case has been used in Zhang et al. (2022) [10], which is a proof of concept for the ensemble-based learning method. In contrast, the current study aims to demonstrate the flexibility of the ensemble method in capturing separated and secondary flows by adjusting the nonlinearity of the turbulence model. Specifically, it is observed here that the ensemble method can learn a linear eddy viscosity model for the separated flow and a nonlinear eddy viscosity model for the secondary flow. This is different from the previous work [10], which is a proof of concept for the ensemble-based learning method. Moreover, here we

Figure 13: Comparison of the model function \(g^{(1)}\) and the combination \(g^{(2)}-0.5g^{(3)}+0.5g^{(4)}\) between the learned and the baseline models for the square duct case
use sparse DNS data to train neural network models, which shows the capability of the ensemble Kalman method to handle sparse data in realistic applications. In contrast, in the previous work [10], full-field flow data approximated with the quadratic model of Shih (1993) [38] is used as synthetic truth, which is not identical to the DNS data.

## 5 Conclusions

This work investigates the physical interpretation of neural-network-based turbulence modeling with the ensemble Kalman method. The observation data, including aerodynamic lift and velocity measurements, are used to train the turbulence model represented with a tensor-basis neural network. The method is applied to the flow around the S809 airfoil and the flow in a square duct. Both cases show that the learned model significantly improves the flow predictions, and the model improvement can be interpreted from a physical viewpoint. In the S809 airfoil, the learned model reduces the eddy viscosity around the upstream boundary layer and captures the appropriate onset of the flow separation, which improves the prediction of the lift force compared to the baseline \(k\)-\(\omega\) model. The learned model can be well generalized to different angles of attack. In the square duct case, the learned model produces a nonlinear eddy viscosity model, which captures the imbalance of the Reynolds normal stress and the in-plane velocity. The ensemble Kalman method can provide appropriate turbulence models based on limited observation data. For the flow over the S809 airfoil, the training method provides an optimized linear eddy viscosity model based on the lift force measurements, which is able to capture the flow separation. In contrast, for the flow in a square duct, the training method provides a nonlinear eddy viscosity model to estimate the anisotropy of Reynolds stress and capture the in-plane secondary flows.

## Appendix A Practical implementation

The practical implementation of the ensemble-based turbulence modeling framework is detailed in this subsection. Given the observation error \(\mathsf{R}\), the data set \(\mathsf{y}\), and the sample variance \(\sigma\), the training procedure

Figure 14: Learned tensor components for the square duct case
is summarized briefly below.

1. Pre-training: To obtain the initial weight \(\mathbf{w}^{0}\) of the neural network, we pre-train the network to be equivalent to a linear eddy viscosity model such that \(g^{(1)}=-0.09\) and \(g^{(2-10)}=0\). The obtained weights \(\mathbf{w}^{0}\) are set as the initial value for model training [2].
2. Initial sampling: We assume that the weights are independent and identically distributed (i.i.d.) Gaussian random variables with mean \(\mathbf{w}^{0}\) and variance \(\sigma^{2}\). As such, we draw random samples of the weights through the formula \(\mathbf{w}_{j}=\mathbf{w}^{0}+\mathbf{\epsilon}_{j}\), where \(\mathbf{\epsilon}\sim\mathcal{N}(0,\sigma^{2})\).
3. Feature extraction: The velocity field \(\mathbf{u}\) and turbulence time scale \(\tau_{s}\) are used to compute the scalar invariants \(\mathbf{\theta}\) and the tensor bases \(\mathbf{T}\) based on the equations (3) and (6). The scalar invariants are normalized and then adopted as the inputs of the neural network. Further, the tensor bases are employed to construct the Reynolds stress by combining with the outputs of the neural network as illustrated in step 4.
4. Evaluation of Reynolds stress: The input features \(\mathbf{\theta}\) are propagated to the basis coefficient \(\mathbf{g}\) with each realization of the weights \(\mathbf{w}\). Then the Reynolds stress can be constructed by combining the coefficient \(g\) and the tensor basis \(\mathbf{T}\) based on Eq. (2).
5. Propagation to mean flow fields: The mean velocity is obtained by solving the RANS equations for each constructed Reynolds stress. Moreover, the turbulence kinetic energy and the dissipation rate are obtained by solving the turbulence transport equations.
6. Update weights of neural networks: The iterative ensemble Kalman method is used to update the weights of the neural network based on Eq. (7). In the scenario of multiple observations, e.g., the S809 airfoil case in this work, the data sets are randomly shuffled and then incorporated sequentially. Specifically, In the S809 airfoil case, the data from two different flow conditions are shuffled to generate a data set with random ordering. Then the observation data is incorporated sequentially in the shuffled order. The observation is reshuffled once the entire data sets are traversed. Besides, for each data, the Kalman update is iterated in an inner loop, and the maximum of the iteration step is set as 3 based on our sensitivity study. The random data ordering can escape from local minima [29] that provide good predictions for one case but inferior results for other cases based on our numerical tests.

If the ensemble variance is smaller than the observation error or the total iteration maximum is reached, the training is considered converged; otherwise, continue to Step 3 until the convergence criterion is met.
## Appendix B Predictive performance for unseen data

We show additional prediction results of the learned model on unseen data for both the S809 airfoil case and the square duct case in this section.

In the S809 airfoil case, the lift force measurements from angles of attack of \(8^{\circ}\) and \(14^{\circ}\) are used for training, which can improve predictions of aerodynamic lift at unseen angles of attack. Here we show that wall pressure prediction can be also improved with the learned model at unseen angles of attack. The wall pressure predictions with the learned model in angles of attack of \(11^{\circ}\) and \(18^{\circ}\) are shown in Figure 15, with comparison to the experimental data and prediction of the baseline \(k\)-\(\omega\) model. It can be seen that the baseline model underestimates the surface pressure on the suction side of the airfoil, which leads to large discrepancies in the predicted aerodynamic lift as presented in Fig. 9. In contrast, the learned model significantly improves the prediction of the wall pressure distribution and eventually the lift coefficient at both the angles, compared to the baseline \(k\)-\(\omega\) model.

In the square duct case, the velocity along profiles of \(y/h=0.25,0.5,0.75\), and \(1.0\) are used for training and lead to local predictive improvement in both the velocity and the Reynolds stress. Here we provide the model prediction at four unseen locations, i.e., \(y/h=0.2,0.4,0.6\), and \(0.8\). The results are shown in Fig. 16, with a comparison to the DNS data and the baseline \(k\)-\(\varepsilon\) model. Apparently, at these unobserved locations, the learned model predicts well the velocity component \(u_{y}\) and the difference of the imbalance of Reynolds normal stress \(\tau_{yy}-\tau_{zz}\). The learned model also yields a non-zero shear component \(\tau_{yz}\), while the baseline \(k\)-\(\varepsilon\) model yields zero shear. The latter is qualitatively incorrect based on the DNS data.

Figure 15: Prediction of wall pressure coefficient \(C_{p}\) at \(\alpha=11^{\circ}\) and \(18^{\circ}\) with the learned model and the baseline \(k\)-\(\omega\) model compared to the experimental data [34] for the S809 airfoil case
## Acknowledgment

XLZ and GH are supported by the NSFC Basic Science Center Program for "Multiscale Problems in Nonlinear Mechanics" (No. 11988102). XLZ also acknowledges support from the National Natural Science Foundation of China (No. 12102435) and the China Postdoctoral Science Foundation (No. 2021M690154). HX acknowledges the support from the National Research Foundation of Korea (No. NRF-2021H1D3A2A01096296) during his sabbatical visit to Gwangju Institute of Science and Technology, where this work was performed.

Figure B.16: Prediction of velocity and Reynolds stress along profiles at \(y/H=0.2,0.4,0.6,0.8\) with comparison among the learned model, the baseline model, and the experimental data, for the square duct case