# Inferring wall-bounded coherent structures from two-dimensional turbulent fields via SHAP analysis

Samuel Molina-Casino, Andrés Cremades, Sergio Hoyas, Jose I. Cardesa, François Chedevergne, Ricardo Vinuesa

► To cite this version:

Samuel Molina-Casino, Andrés Cremades, Sergio Hoyas, Jose I. Cardesa, François Chedevergne, et al. Inferring wall-bounded coherent structures from two-dimensional turbulent fields via SHAP analysis. European Journal of Mechanics - B/Fluids, 2025, 114, pp.204304. 10.1016/j.euromechflu.2025.204304. hal-05075817v3

HAL Id: hal-05075817

https://hal.science/hal-05075817v3

Submitted on 15 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Inferring wall-bounded coherent structures from two-dimensional turbulent fields via SHAP analysis

Samuel Molina-Casino<sup>a, ID,*</sup>, Andrés Cremades<sup>b, ID</sup>, Sergio Hoyas<sup>c, ID</sup>, José I. Cardesa<sup>a, ID</sup>,
François Chedevergne<sup>a, ID</sup>, Ricardo Vinuesa<sup>b, ID,*</sup>

## ARTICLE INFO

Dataset link: https://github.com/KTH-FlowAI

Keywords:

Turbulence

Explainable deep learning

Shapley values

Coherent structures

## ABSTRACT

The analysis of coherent structures in turbulent flows is essential for understanding and controlling turbulence. Explainable Deep Learning (XDL) has been successfully used to identify the most important flow structures in turbulent channels (Cremades et al., 2024), but its computational cost grows significantly as Reynolds number increases. This study adapts the XDL framework to assess whether observations from time-resolved volumetric data of direct numerical simulations (DNS) can be inferred using two-dimensional (2D) slices. Limiting the analysis to 2D subsets reduces the cost of three-dimensional (3D) post-processing and facilitates its application to planar particle image velocimetry (PIV) experimental data. We employ SHapley Additive exPlanations (SHAP), a game-theoretic method, to identify dynamically relevant flow regions by assigning importance scores. Using DNS data from a turbulent channel flow at a friction Reynolds number of 125, we compare SHAP-based structures derived from 3D and 2D data with classically studied structures, including Q events, streaks, and vortex clusters. Our findings show that SHAP-based structures from 2D data statistically agree with those from full 3D fields and previously reported in the literature. These results underscore the potential of SHAP analyses for studying coherent structures with 2D data, enabling experimental and DNS research at higher Reynolds number.

## 1. Introduction

Although wall-bounded turbulence has been studied for decades [1], it remains a dynamic field of research, continually drawing diverse methodologies to address its unresolved questions. The industrial implications of a breakthrough in this area from a technological and economic point of view are being caught up by the more recent and urgent need to decrease the impact of mankind on climate change [2]. However, to make a significant contribution to the advancement of modeling or controlling [3] wall-bounded turbulent flows, the primary challenge lies in deepening the understanding of the underlying flow physics. Achieving this requires both extensive data availability and the development of effective analysis techniques.

In terms of data availability, the increase of computational power and the development of highly efficient numerical schemes has enabled two types of scale-resolving simulations: direct numerical simulations (DNS) [4,5], where all the turbulence scales are resolved, and large-eddy simulations (LES) [6], where only the smallest scales are modeled. While these methods provide unprecedented insight into turbulence

dynamics, their prohibitive computational cost restricts application to Reynolds numbers below industrial relevance or limits geometric complexity to academic flow configurations. On the other hand, data provided by experimental techniques is much less limited in terms of Reynolds numbers. Hot-wire anemometry, laser-Doppler velocimetry, and particle image velocimetry (PIV) [7] are techniques that can be found in studies handling Reynolds numbers well in excess of those in the most ambitious DNSs [8]. Yet experimental methods often provide limited access to certain flow quantities or restrict measurements to small regions of the flow field [9,10].

One approach to study wall-bounded turbulence relies on the analysis of coherent structures [11], such as Q events [12], streaks [13], and vortices [14], which are associated with regions of high shear, wall-parallel velocity fluctuations, and rotation regions, respectively. However, the precise connection between these classical structures and their impact on the flow dynamics remains elusive. This is particularly true when it comes to ranking them in terms of dynamical importance.

<footnote><sup>a</sup> ONERA/DMPE, Université de Toulouse, Toulouse, F-31055, France</footnote>

<footnote><sup>b</sup> FLOW, Engineering Mechanics, KTH Royal Institute of Technology, Stockholm, SE-100 44, Sweden</footnote>

<footnote><sup>c</sup> Instituto Universitario de Matemática Pura y Aplicada, Universitat Politècnica de València, Valencia, 46022, Spain</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: samuel.molina_casino@onera.fr (S. Molina-Casino), rvinuesa@mech.kth.se (R. Vinuesa).</footnote>
Recent studies have leveraged Explainable Deep Learning (XDL) to analyze coherent structures in turbulence. XDL aims to make machine-learning models more transparent by identifying which features contribute most significantly to a given prediction. In particular, SHapley Additive exPlanations (SHAP) [15,16] offer a data-driven approach to ranking turbulence structures by their influence on the future evolution of the flow, providing a more comprehensive understanding of the dynamics. SHAP-based analyses by Cremades et al. [17,18] have shown that the structures identified as most important through SHAP values do not always coincide with the classical structures obtained by traditional methods, revealing a clear dependency of their importance for dynamics prediction with wall-normal distance. Furthermore, their findings revealed that the most influential structures are not necessarily always those that contribute the most to the Reynolds shear stress, notably highlighting the relevance of streamwise streaks in the buffer layer. Crucially, their work reveals that the classically studied coherent structures only paint a partial picture of wall-bounded turbulence.

Insights from XDL methodologies like SHAP extend beyond refining our understanding of coherent structures, now driving advancements in wider areas of turbulence research. These techniques show emerging potential to inform more targeted active flow control strategies; for example, by guiding Deep Reinforcement Learning (DRL) systems [3, 19] towards more effective policies for drag reduction or separation control, as exemplified in the recent work by Beneitez et al. [20]. Concurrently, in turbulence modeling, SHAP provides valuable tools for interpreting and potentially improving data-driven closures for Reynolds-averaged Navier-Stokes (RANS) equations. By highlighting flow features critical to model performance or discrepancies, SHAP can aid in developing more physically sound models [21,22]. The promise of such applications underscores the importance of developing versatile and computationally accessible SHAP-based analytical frameworks, which is a central motivation for the present study.

The analysis of the importance of flow structures based on SHAP values has been primarily applied to 3D flow fields [17], using point-by-point calculations to evaluate the relevance of each individual grid point. Nevertheless, studies based on SHAP values have so far been restricted to datasets with low friction Reynolds numbers due to the computational cost of this technique.

This situation naturally raises the question of whether meaningful insights into flow physics can be gained by analyzing two-dimensional subsets of volumetric flow data. Indeed, if the information excluded when analyzing data subsets does not hinder access to the conclusions inferred from the full 3D fields, then two compelling opportunities arise from the viewpoint of using data-driven techniques: analyzing experimental results at higher Reynolds numbers, or carrying out post-processing on larger existing 3D datasets of which only a strategically selected subset is kept.

In a previous study by Cremades et al. [18], a 2D analysis based on PIV data was conducted showing good agreement with the 3D-based conclusions included in the same paper. This first inspection relied on a pre-segmented field of Q events, which limited the exploration of the interplay between different structures, such as streaks and vortices.

In the present study, we investigate extending the point-by-point SHAP method to 2D analyses, comparing the insights obtained from a 2D framework with those from 3D analyses when trying to identify regions of high dynamic importance. Furthermore, the investigation also examines whether the correlation between SHAP-identified structures and classical coherent features (Q events, streaks, and vortices) remains robust when based solely on 2D data. Finally, the role of spanwise velocity in accurately segmenting structures by importance is explored. Any significant differences in the probability distributions of the obtained structures between the two approaches are considered, along with potential interpretations and justifications for such discrepancies.

The paper is structured as follows. Section 2 details the numerical setup of the simulation used to generate the high-fidelity 3D database, followed by a description of the methodology leading to SHAP value

<figure><img src="image_2.png" /><figcaption>Fig. 1. Numerical domain of the turbulent channel. The flow, with streamwise and spanwise periodicities, is from left to right. The purple and green planes mark the slices at $z = 0$ and $z = \pi/2$ of the domain used for the 2D temporal database. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

estimates and an outline of the coherent structure eduction process. Section 3 discusses the results and analyzes the impact of dimensionality reduction by comparing SHAP-based structures derived from both 2D and 3D implementations. The final section contains the conclusion of this work.

## 2. Methods

### 2.1. Numerical simulation and database

A canonical channel flow as depicted in Fig. 1 is used in the study. The DNS simulation and main database was previously obtained and described by Cremades et al. [18] using the LISO code [23]. The channel comprises two parallel planes with dimensions $L_x = 2\pi h$ and $L_z = \pi h$, separated by a vertical distance of $L_y = 2h$. Here, x, y, and z refer to the streamwise, wall-normal, and spanwise coordinates, respectively. Following the Reynolds decomposition, the instantaneous velocity fields, noted with upper-case $U = (U, V, W)$, are decomposed as $U(x, y, z, t) = \bar{U}(y) + u(x, y, z, t)$, where the bar represents averaging in x, z, and time, and the lower-case $u$ is the instantaneous turbulent fluctuation vector, composed by the three components, $u = (u, v, w)$. The simulation was performed at a friction Reynolds number of $Re_x = u_x h / \nu = 125$, where $\nu$ is the kinematic viscosity and $u_x = \sqrt{\tau_{xx}} / \rho$ is the friction velocity formed with the wall-shear stress $\tau_{xx}$ and the fluid density $\rho$. Quantities normalized by viscous units are denoted with the superscript + when scaled in terms of the friction velocity $u_x$ or the viscous length $\ell^* = \nu / u_x$.

The computational domain size was chosen to be sufficiently large to ensure that one-point statistics converge to their large-scale [24]. The associated grid comprises $n_x \times n_y \times n_z = 192 \times 201 \times 96$ points, yielding resolutions of $\Delta x^+ \approx 4$, $\Delta z^+ \approx 4$, and $\Delta y^+ = 1.5\eta$, being $\eta$ the Kolmogorov scale [25].

The database comprises 9000 instantaneous velocity fields sampled at regular intervals of $\Delta t^+ = 5$. This sampling rate has been demonstrated to be sufficient for capturing the flow dynamics at the specified $Re_x$ [18]. From this extensive original dataset, three distinct databases are generated for subsequent analysis:

* 3D database (3D). This dataset contains the full three-dimensional information. It provides a comprehensive description of the turbulent flow field and serves as the reference for fully 3D analyses.

* 2D database (2D-uv). This dataset is constructed by extracting planar $xy$ slices from the 3D database and retaining only the in-plane velocity components ($U$ and $V$).

* 2D database (2D-uvw). Similar to the 2D-uv dataset, this database involves extracting $xy$ slices but retains all three velocity components.
Classic two-component 2D particle image velocimetry (2C2D PIV), widely used in laboratory experiments, captures only two velocity components due to its single-camera setup, which is limited to in-plane motion measurements [26,27]. This configuration is mimicked by the 2D-uv database. However, as three-component 2D PIV (3C2D PIV, also known as stereoscopic PIV) [28] becomes more widespread, having separate 2D-uv and 2D-uvw datasets facilitates an incremental comparison between the analyses allowed for by the different planar PIV variants.

The 2D databases are assembled using two xy cuts of the full 3D computational field. These planes are selected at the largest distance possible, z = π/2, to avoid spurious correlation. Hence, by choosing slices at z = 0 and at z = π/2 (center of the channel) the amount of redundant information available during training is minimized, based on two-point-correlation analysis.

Measurements near the wall present significant challenges in experimental fluid dynamics. In many PIV studies of wall-bounded turbulence, reliable velocity data can typically only be obtained for wall-normal distances greater than approximately y⁺ ≈ 20. For instance, Kähler et al. [29] highlight that high velocity gradients, limited optical access and seeding issues in the near-wall region prevent accurate measurements of instantaneous velocity data very close to the wall. These challenges are intensified by reflections and image distortions from the wall surface, which further degrade measurement quality. PIV techniques also face significant limitations in resolving small-scale turbulence structures near walls due to insufficient spatial resolution and high uncertainties in regions with steep velocity gradients [26]. Achieving reliable near-wall turbulence characterization often requires advanced techniques like micro-PIV or volumetric Particle Tracking Velocimetry (PTV), which are resource-intensive and prone to errors in this critical region [30]. For a survey of additional techniques to measure the near-wall flow, the reader is referred to Ref. [31].

For all these reasons, the analysis will be restricted to y⁺ > 20, thus minimizing uncertainties associated with near-wall measurements and ensuring that the findings relate to the most frequently sampled region in experiments.

## 2.2. SHAP framework

SHAP (SHapley Additive exPlanations) belongs to a class of additive-feature-attribution methods that provide a framework for the interpretation of complex machine-learning models. These methods quantify how each input feature contributes to the model's predictions. In this study, features are defined as the individual grid points within the velocity field.

SHAP values derive from cooperative game theory and are based on Shapley values [15]. This method approximates highly non-linear predictive models with simpler, additive representations, making interpretation more accessible [32]. Conceptually, when multiple features contribute to a prediction, their individual importance is determined by measuring how much they enhance the model's output when added to different feature combinations. For any given instance, the overall outcome is decomposed into a baseline (what the model would predict without any features), $\phi_0$, and the incremental contributions (SHAP values) of features when included, denoted as $\phi_i$. Within this binary framework, each feature ($z'_i$) is either present or absent, and its SHAP value indicates how much its inclusion shifts the prediction from the baseline. Mathematically, the explainable model is formulated as:

$$g(z') = \phi_0 + \sum_{i=1}^{N} \phi_i z'_i; \quad z'_i \in \{0, 1\} \qquad (1)$$

where the $g(z')$ represents the simplified linear model as a function of the presence of each individual feature.

Although Shapley values offer a rigorous framework for explainability, their exact computation becomes impractical for high-dimensional

models, such as those used in turbulence simulations, due to their exponential scaling with feature count. To address this limitation, various approximation methods such as kernel-SHAP [32], deep-SHAP [16], tree-SHAP [33], and gradient-SHAP [34] have been developed [16, 22]. These techniques leverage model-specific properties or sampling strategies to efficiently estimate SHAP values while reducing the computational burden.

In this work, the gradient-SHAP framework is adopted. This method estimates the contribution of each feature using the already trained predicting model's gradients and computes expected gradients (EG) to assess how deviations from an uninformative baseline field affect the model output. The core idea of expected gradients [34] is to approximate the SHAP values by integrating the model's output gradients with respect to the input features along a path from the baseline input to the actual input. This path integration is typically performed by sampling m points along a straight line between the baseline and the actual input and averaging the gradients at these interpolated points. This process essentially measures the local sensitivity of the model's output to each feature. Unlike kernel-SHAP, which is model-agnostic, gradient-SHAP exploits the model's gradient information for computational efficiency. The resulting EG directly correspond to the SHAP values $\phi_i$ in Eq. (1), quantifying the additive contribution of each grid point's velocity components to the overall prediction error. A more detailed overview of the gradient-SHAP method is provided in Appendix A.

Fig. 2 summarizes the workflow of the methodology. The 3D temporal DNS dataset serves as ground truth, from which 2D temporal slices along the z-direction are extracted to form the 2D-uv and 2D-uvw datasets. Velocity fluctuations at time t are computed and used to train a U-Net model to predict the velocity field at t⁺ + Δt⁺. Once the predictive model is obtained, it is adapted to output the mean squared error (MSE) between the predicted and ground-truth velocity fields. Using the non-informative zero-fluctuation velocity field as the baseline reference, each grid point's contribution is evaluated by effectively setting it to a zero-fluctuation state. The gradients of the MSE with respect to the input field are then computed, yielding SHAP values that provide a measure of how small changes in each grid point influence the prediction error. Finally, a percolation-based thresholding method segments the SHAP values into distinct structures, delineating fluid regions of high predictive importance, analogous to classical coherent structure analysis.

By applying this methodology to each grid point of the input field, a field of SHAP values, $\phi(x, y, z, t)$, is generated. This field provides the importance of the streamwise ($\phi_u$), wall-normal ($\phi_l$), and spanwise velocities ($\phi_w$) when included. The absolute value of the SHAP vector is used as a combined metric to identify regions of high importance, which will be analyzed in detail in the next section.

The described workflow relies on the presence of a predictive model. To achieve this, a U-Net deep-learning architecture [35] is employed since this type of model is widely used for pattern recognition and segmentation tasks. The objective of the network is to forecast the velocity field Δt⁺ units in the future. The chosen topology consists of a contracting (downsampling) and an expansive (upsampling) path with four levels, as illustrated in Fig. 3. Table 1 provides a summary of its main parameters.

During the downsampling path, two convolutional layers [36] per level are applied, with an increasing number of 3 × 3 filters. The number of filters doubles at each level, starting with 24 filters in the first level for the 3D case and 32 for the 2D cases. Downsampling is performed using 2 × 2 average pooling layers. The upsampling path follows a symmetric structure, consisting of a transposed convolution followed by a standard convolution. At each level, the output from the corresponding downsampling path is incorporated, ensuring preservation of fine-scale spatial details. The Rectified Linear Unit (ReLU) is used as activation function in all layers except for the last one, where a sigmoid function is applied to generate the final velocity field prediction. The network is trained using the mean-squared-error (MSE)
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Workflow diagram summarizing the methodology. The 3D DNS database is used to extract 2D slices, which are then employed to train a U-Net model. The model's predictions are analyzed with gradient-SHAP to evaluate each grid point's contribution using their zero-fluctuation state as non-informative baseline. Finally, a thresholding method segments important fluid regions into coherent structures.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Schematic of the U-Net architecture used for the prediction model. The parameter $n_c$ refers to the number of channels in the input and output tensors, that is, the velocity components of the velocity ($n_c = 2$ for the 2D-uv case, and $n_c = 3$ for the 2D-uvw and 3D cases). The parameter $n_f^2$ represents the number of filters after each convolution operation, doubling at each level of the network. The number of filters used in the first level is 24 for the 3D case and 32 for the 2D cases.</figcaption></figure>

loss function, which minimizes the difference between the predicted field at $t^+ + \Delta t^+$ and the true field obtained from the DNS simulation.

The inputs of the U-net are the normalized instantaneous velocity fields obtained from the DNS simulation. The input tensor has a shape of $192 \times 201 \times 96 \times n_u$, and $192 \times 201 \times n_u$ for the 3D and 2D cases, respectively, where $n_u$ indicates the number of velocity components used as input. From the complete dataset, 6400 fields are used for training and validation (with an 80/20 split), while the remaining fields are reserved for testing and analysis.

As shown in Table 1, the prediction errors for all cases remain within the same order of magnitude. The slightly higher errors observed for the 2D cases are likely due to the lack of transverse continuity in the provided inputs, as well as the reduced information available in the 2D slices compared to the 3D volume. Despite this limitation, the observed error remains well within the acceptable range to proceed with the analysis.

<caption>**Table 1**</caption>

Summary of prediction networks and their average error ($\epsilon$). The error $\epsilon$ represents, for each velocity component, the volume-weighted average difference between predicted and DNS velocity fields normalized by the maximum velocity of the database and computed across all testing snapshots.

<table><thead><tr><th>Case</th><th>3D</th><th>2D-uv</th><th>2D-uvw</th></tr></thead><tbody><tr><td>Input field</td><td>u, v, w</td><td>u, v</td><td>u, v, w</td></tr><tr><td>Input tensor</td><td>192 x 201 x 96 x 3</td><td>192 x 201 x 2</td><td>192 x 201 x 3</td></tr><tr><td># train fields</td><td>5120</td><td>10 240</td><td>10 240</td></tr><tr><td># validation fields</td><td>1280</td><td>2560</td><td>2560</td></tr><tr><td># filters</td><td>24, 48, 96, 192</td><td>32, 64, 128, 256</td><td>32, 64, 128, 256</td></tr><tr><td># parameters</td><td>3.29 M</td><td>1.95 M</td><td>1.95 M</td></tr><tr><td>ε<sub>u</sub></td><td>1.00%</td><td>1.86%</td><td>1.68%</td></tr><tr><td>ε<sub>v</sub></td><td>0.77%</td><td>1.95%</td><td>1.86%</td></tr><tr><td>ε<sub>w</sub></td><td>0.85%</td><td>-</td><td>2.10%</td></tr></tbody></table>

## 2.3. Definition of coherent structures

Coherent structures are defined as connected regions of the flow satisfying a condition based on a threshold that is applied to a scalar
property of the flow. The algorithm used to identify connectivity between points is defined in terms of the six orthogonal neighbors in the cartesian mesh of the DNS [37]. The threshold definition varies from one structure type to another, and they will be defined in what follows.

The SHAP value at each point, denoted as $\phi(x, y, z, t)$, quantifies the importance of each velocity component in predicting the evolution of the flow at the next time step. To identify coherent structures based on these importance scores, a threshold criterion is applied to compare the absolute magnitude of $\phi$ against its scaled root mean square (RMS) value as follows:

$$ \sqrt{\phi_u^2(x, y, z, t) + \phi_v^2(x, y, z, t) + \phi_{uv}^2(x, y, z, t)} > H \sqrt{\phi_u^2(y) + \phi_v^2(y) + \phi_{uv}^2(y)}. \quad (2) $$

$H$ is the hyperbolic-hole size [38], or percolation index. This parameter determines the segmentation process and will be discussed in detail in the subsequent subsection.

Moving on towards the more classical coherent structures that will be compared against SHAP structures, Q events are regions of intense Reynolds shear stress characterized by a strong correlation between horizontal and vertical velocity fluctuations [12]. These structures are associated with momentum transfer and the production of turbulent kinetic energy. They are defined by:

$$ |u(x, y, z, t)v(x, y, z, t)| > \beta u'(y)v'(y), \quad (3) $$

where $\beta$ is the corresponding hyperbolic-hole size and the primes denote RMS quantities, i.e., $u' = \sqrt{u^2}$. Q events are further categorized based on quadrant analysis [39] into four types depending on the signs of the streamwise and wall-normal velocity fluctuations. Outward interactions ($u > 0, v > 0$) are fast structures moving upwards; ejections ($u < 0, v > 0$) are slow lumps of fluid rising from the wall towards the bulk; inward interactions ($u < 0, v < 0$) have slow horizontal velocity moving downwards; and finally sweeps ($u > 0, v < 0$) are fast fluid lumps moving towards the wall. Ejections and sweeps have been previously shown to be the most dominant structures for momentum transport [18].

The second type of classical structures under consideration are the so-called streamwise streaks, characterized by alternating low and high streamwise velocity regions aligned with the primary flow direction. Their definition is given by:

$$ \sqrt{u^2(x, y, z, t) + w^2(x, y, z, t)} > \alpha u_r, \quad (4) $$

where $\alpha$ is again a hyperbolic-hole size. Streamwise streaks can be categorized into high ($u > 0$) and low ($u < 0$) streamwise velocity regions. Low-speed streaks result from low-momentum fluid being lifted by streamwise vortices [40], playing a crucial role in the near-wall turbulent self-sustaining cycle. High-speed streaks, often appearing in a staggered distribution with low-speed streaks, are essential for generating and transporting turbulent kinetic energy in the boundary layer [13].

Finally, vortex clusters, or vortices, are regions where fluid rotation exceeds shear, indicating dominance of the rotational component ($\Omega$) over the strain component ($S$) in the instantaneous velocity-gradient tensor, $\nabla u = S + \Omega$. Several criteria to identify these clusters exist, however, this work employs the $\Delta$ criterion defined by Chong et al. [14]. This criterion is based on the second ($Q = 1/2 (\|\Omega\| - \|S\|)$) and third ($R = -\det(\nabla u)$) invariants of $\nabla u$, which for incompressible flows results in:

$$ \Delta(x, y, z, t) = Q^3 + \frac{27}{4} R^2 > \gamma \Delta'(y), \quad (5) $$

with $\gamma$ being the associated hyperbolic-hole size.

Unlike SHAP-based structures, classical streaks and vortices require the three-dimensional field to be computed. Consequently, these structures cannot be recovered from the 2D-$uv$ dataset. While a single component of vorticity can be resolved with planar PIV techniques,

the three-dimensional orientation cannot be assumed a priori to align such that most vorticity is contained in the resolved component. Hence, studying vortical structures in the 2D analyses is not feasible. For this reason, to assess the consistency between 2D SHAP-based observations and classical structures, the latter are computed from the full 3D dataset, as will be presented in Section 3. The comparison is then performed by extracting the same spanwise slices used to generate the 2D subsets.

## 2.4. Coherent structure eduction

A critical aspect of clustering domain points into connected regions like coherent structures is analyzing the sensitivity of the resulting regions to the threshold value. For SHAP-based structures, clustering is performed by comparing the local SHAP value to a threshold defined as the product of its standard deviation and $H$. Given that domain statistics vary with the wall-normal coordinate $y$, the threshold includes a $y$-dependent contribution stemming from the RMS of the SHAP values, $|\phi'|$ (Eq. (2)).

Regarding the constant component of the threshold, the hyperbolic-hole value $H$ is usually selected to maximize the number of structures in the domain [38]. Following this approach, $H$ is set to 1.62 for the 3D dataset based on the 3D percolation analysis. For 2D datasets, however, a more thorough examination is necessary. Fig. 4 shows the percolation analysis for both 3D and 2D-$uv$ cases, revealing a clear shift towards lower $H$ values in the 2D dataset.

The two main factors likely to contribute to this shift are the dimensional discrepancy and the ambiguity in structure identification. The first one comes from the fact that in 2D a structure's relative area does not directly correspond to its relative volume in 3D. For instance, structures aligned with the $z$-plane may appear large in a 2D slice but much smaller in the 3D volume if their spanwise dimension is limited. This can result in higher $V/V_{tot}$ ratios (volume with respect to the total identified structures), leading to artificially inflated or deflated 2D values.

The second factor is that a single 3D structure with significant spanwise extent may intersect a given $z$-plane at multiple spatially separated lobes. In a 2D slice, it becomes impossible to determine whether these separate intersecting areas belong to the same 3D structure or to different ones. At lower $H$ values, where segmented structures are larger, more intersections may occur. This increases the apparent number of structures in 2D, shifting the 2D $N/N_{tot}$ curve towards lower $H$ values.

Therefore, selecting the $H$ value in both 2D and 3D as that which maximizes the number of structures will lead to disagreeing thresholds applied to the same scalar quantity, thus blurring the intended analysis that the present study proposes geared towards 2D vs. 3D comparisons. Rather than choosing the same threshold as that in 3D, or maximizing the number of structures in 2D, the value of $H$ for the 2D datasets is chosen so that it leads to an overall area of SHAP structures that is equal to that of 2D slices of the 3D field segmented at the (3D) optimal threshold. This approach yields very similar fields of SHAP coherent structures when comparing the structures obtained in 2D against the spanwise slices of the 3D ones, as shown in Fig. 5. Although the areas of the fields are similar, it is important to note that several identified structures may belong to the same object in the 3D labeling, whereas in the 2D labeling they are treated as distinct structures. Nevertheless, this distinction does not affect the point-by-point analysis that will be carried out subsequently.

## 3. Results and discussion

After segmenting the flow domain into discrete regions, their physical and geometrical characteristics are analyzed to enable qualitative and quantitative comparison. Fig. 6 presents the joint probability density function (JPDF) of inner-scaled streamwise velocity fluctuations
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Comparison of the percolation curve for the 3D case (left) and the 2D-*uv* case (right). The *H* value to maximize the number of structures is 1.62 for the 3D case and 0.56 for the 2D. Two curves are represented: *V*/*V*<sub>tot</sub>, the ratio of the volume of the largest structure to the volume of all identified objects (areas in 2D); and *N*/*N*<sub>tot</sub>, ratio of the number of identified structures to the maximum number of them.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Comparison of the segmented binary field of SHAP values for an instantaneous snapshot. (Left) Results for a spanwise slice from the 3D dataset using the optimal percolation index (*H* = 1.62) determined via percolation analysis in 3D (maximizing the number of structures). (Right) Segmentation for the 2D-*uv* dataset using the *H* value derived from its own optimal percolation analysis (■), together segmented field using the percolation index that matches the overall total area of the structures between the 2D-*uv* fields and the 3D fields (■). Note that the spanwise slice for the 3D figure corresponds to the same spatial location as that of the 2D datasets. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

and wall-normal distance. To ensure fair comparison across datasets, the identified 3D structures have been sliced at the same spanwise positions used in the 2D cases. Similarly, Q events, streaks, and vortex structures are analyzed following the same approach as the SHAP-based ones. Preliminary checks on the planar data confirmed that the JPDFs of Q events and streaks computed from the 2D datasets (with streaks only accessible for 2D-*uvw*) aligned with those obtained from the 3D dataset. However, as mentioned in Section 2.3, the results shown for classical structures are based on the full 3D domain to ensure a comprehensive comparison.

**Fig. 7** quantifies the spatial overlap between SHAP-based structures and classical coherent structures for each dataset. The top-left panel of Fig. 7 illustrates an example of these overlapping regions for a small section of the domain, where black areas denote the shared volumes used in the coincidence calculations. This coincidence calculation is purely geometrical: for each pair of compared structures (e.g., SHAP and streaks), their coincidence is calculated as the ratio between their intersection (number of shared grid points, $\phi \cap$ streaks in the example) and their volume. For short, it will be expressed as $\phi \cap * / \phi$ when obtained relative to the volume of the SHAP-based structures, and $\phi \cap * / *$ when relative to the classic structure's volume, the symbol * representing any given classical structure.

When interpreting the coincidence results presented in Fig. 7, it is crucial to recall the primary aim of the SHAP-based methodology: to obtain an unbiased identification of dynamically important regions based directly on the flow data, rather than relying on or attempting to conform to a-priori definitions of classical coherent structures. Consequently, the coincidence analysis serves an evaluative purpose: to understand the relationship and overlap between SHAP-identified regions and these established classical structures, rather than being a prescriptive tool for empirically deriving the SHAP structures themselves.

A prior study by Cremades et al. [17] examined a larger channel ($8\pi h \times 3\pi h \times 2h$) at the same $Re_\tau$, considering only the 3D domain. Their

analysis demonstrated that SHAP-based structures offer a complete representation of turbulence dynamics. Classical structures, on the other hand, capture only part of the picture, exhibiting a degree of overlap with SHAP structures that depends on wall-normal distance. While Q events (mainly sweeps and ejections) and SHAP events were found to feature a high level of coincidence throughout the studied range of $y^+$, streaks represented most of the dynamics captured by the SHAP-based structures in the range $10 \le y^+ \le 40$, encompassing the lower end of the $y^+$ range studied in the present paper. Meanwhile, vortices were found to have limited overlap with SHAP-based structures, contributing at most 20% at the furthest distances from the wall.

The trends observed in the previous work [17] based on a 3D analysis are found to be consistent with current observations relying on a 2D approach, as shown by the JPDFs in Fig. 6. In all cases, there is a predominant probability of low streamwise velocity fluctuation structures, which is in very good agreement with the ejections shown in the Q events. Previous results demonstrated similar importance (per unit volume) for sweeps and ejections [18]. However, sweeps were found to be dynamically more important closer to the wall, while ejections took over for $y^+ > 15$ [17]. The results obtained in 2D support these findings, highlighting the role of ejections in the $y^+ > 20$ region. The high coincidence values for these structures also align with previous turbulence studies [11,12,41], which identified sweeps and ejections as primary mechanisms for momentum transport near the wall.

**Fig. 7** (bottom-left) shows that between 50% and 75% of SHAP-based structures coincide with Q events, mainly ejections. However, near $y^+ = 20$, streaks account for nearly 100% of the SHAP-based volume. These streaks primarily consist of low-speed structures, as shown in the JPDFs of Fig. 6. When calculating coincidence relative to the volume of streaks (Fig. 7, bottom-right), the intersection represents 25%-45% at those $y^+$ ranges, demonstrating that streaks occupy greater volume than the intersection. A similar conclusion can be drawn about
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Joint probability density functions (JPDFs) of various structures as a function of the wall-normal distance and the inner-scaled streamwise velocity fluctuation of the segmented domains. The results span a $y^+$ range between 20 and 125. (Top) JPDFs for SHAP-based structures in the 3D, 2D-*uv* and 2D-*uvw* datasets (left to right); (bottom) classical coherent structures: Q events, streaks and vortex clusters (left to right). To ensure consistency, the 3D results include only the slices corresponding to those used in the 2D cases. JPDFs of the wall-normal and spanwise velocity components are available in Appendix B.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Comparison of spatial coincidence between SHAP-based and classical structures for the 3D, 2D-*uv*, and 2D-*uvw* datasets. (Top-left) Example of a segmented domain, with dark regions (■) indicating intersections between SHAP-based and classical structures, ■ for Q events, ■ for streaks, and ■ for vortices. (Top-right) Spatial overlap between 3D SHAP-based structures and those obtained in the 2D cases. (Bottom-left) Coincidence between classical turbulence structures (Q events, streaks, and vortices) and SHAP-based structures in 3D, 2D-*uv*, and 2D-*uvw* cases relative to the volume of the SHAP-based structures. (Bottom-right) Coincidence relative to the volume of the classical structures. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

the volume of Q events, where the intersection only represents about 25% of their volume across all wall-normal distances. Both ejections and low-speed streaks represent fluid regions slower than the bulk

flow, and they appear to have a stronger impact on the predicted flow evolution. This trend is consistent across the 3D and the two 2D datasets.
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Comparison of SHAP-based and classical structures for 3D and 2D-y datasets, computed on XZ planes at $y^+$ values of 20 (left) and 100 (right). (Top) JPDF for all structures at each height, with continuous and dashed lines for $y^+ = 20$ and 100, respectively. (Middle) Spatial coincidence between structures, with bars normalized by SHAP (solid) or classical (hatched) volumes. (Bottom) Segmented structures on wall-normal planes, highlighting the spatial overlap (dark color). Color scheme matches that in Fig. 7. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

The limited overlap between vortices and SHAP-based structures in 3D suggests that SHAP methodologies primarily account for momentum transport mechanisms rather than rotational structures as the main predictors of flow evolution. This trend is similarly observed in the 2D datasets, although with slightly reduced coincidence values. It is crucial to reiterate that SHAP-based structures and classical structures are fundamentally distinct [18]; they are defined by different criteria to capture different facets of flow physics. Therefore, the selection of the hyperbolic-hole parameter ($H$) for SHAP structures is not aimed at maximizing their coincidence with classical structures by tuning $H$. Instead, as detailed in Section 2.4, the $H$-selection methodology was specifically designed to ensure consistency between the 2D and 3D SHAP analysis. This approach is pivotal for a meaningful comparison between the 2D and 3D SHAP analyses themselves. Consequently, it ensures that any subsequent comparisons of these SHAP structures with classical structures are conducted on an impartial basis, rather than reflecting an attempt to artificially force agreement.

Fig. 7 (top-right) shows that spatial coincidence between SHAP-based structures in 3D and 2D lies between 50% and 70%. This range likely stems from the point-by-point method used to quantify the spatial coincidence between structures. In particular, the irregular, non-smooth boundaries of SHAP-based structures (compared to the smoother perimeters of classical structures), pose challenges for a direct comparison. Even if the 2D and 3D SHAP structures have similar location and dimensions, they may still yield different point distributions, providing lower overlap values. Another possible factor is the percolation issue discussed in Section 2.4, which affect how structures are identified and compared across dimensions. Nevertheless, results indicate that the conclusions regarding the regions of dynamical importance and their common features with classical structures remain

aligned with those deduced from the 3D analysis, revealing that dimensionality reduction is not reducing the inferring capabilities of the framework.

A somewhat unexpected result in Fig. 7 is that including the spanwise velocity component in the 2D-$ww$ case does not significantly improve agreement with 3D results compared to the 2D-$wv$ case. This is a positive finding as it can significantly ease the experimental data acquisition process required to obtain data on which to apply the present methodology.

To further investigate the impact of the spanwise coordinate for various wall-normal distances, a supplementary dataset called 2D-y is added. This dataset comprises wall-normal slices (plane $y$) at $y^+$ values of roughly 20 and 100, where all three velocity components ($u$, $v$, and $w$) are considered. Experimentally, such data could be equally acquired using planar PIV, hence its relevance to the present study. Fig. 8 shows the results obtained on such wall-parallel planes. The structures identified at each height (bottom row) showcase the dominant presence of streaks at $y^+ = 20$ and their absence at $y^+ = 100$, supporting the relevant role of their volume when interpreting the coincidence curves shown in Fig. 7.

Once again, the PDFs of SHAP-based structures, both 3D and 2D-y, exhibit similar trends (top row of Fig. 8). At both heights the PDFs peak at positive and negative fluctuation velocities featuring a slight asymmetry in favor of higher probabilities for negative fluctuations. An increase in wall-normal distance appears to reduce the fluctuation velocity (positive or negative) at which the probabilities peak, with an even stronger emphasis on negative fluctuations than observed closer to the wall. Remarkably, the PDF of SHAP values at $y^+ = 100$ resembles a combination of the different classical structures at various ranges of velocity fluctuations, confirming that their relative importance depends not only on wall-normal distance but also on
the magnitude of the horizontal velocity fluctuation. This relationship diminishes closer to the wall (left), where SHAP structures tend to exhibit dominant negative horizontal fluctuation, with positive fluctuations being slightly less prevalent. In terms of relative topological coincidence (middle row), streaks again represent the highest volume portion of the SHAP-based structures at height $y^+ = 20$. This proportion drastically diminishes further from the wall, where Q events assume greater relative coincidence.

## 4. Conclusions

This study evaluates the effectiveness of SHAP methodologies in 2D datasets for providing accurate data-driven interpretations of turbulence dynamics. A U-Net model trained on DNS data of turbulent channel flow was employed to analyze SHAP-based importance fields in both 2D (with and without spanwise velocity) and 3D configurations. The comparative analysis between 2D and 3D approaches reveals a good agreement for $y^+ > 20$, confirming that 2D SHAP structures capture key momentum-transport mechanisms identified in 3D analyses even with reduced data dimensionality. Notably, Q events (mainly ejections at these heights) and streamwise streaks exhibit a high degree of spatial coincidence with SHAP-based structures, validating the framework's ability to identify dynamically significant regions in turbulence.

Nevertheless, certain aspects require further investigation. Future studies should assess the sensitivity of our findings to variations in Reynolds number, particularly at higher values where turbulence exhibits increased complexity. Another key aspect is the development of improved methods for determining the optimal hyperbolic-hole, H, to ensure robust structure clustering based solely on 2D flow fields. Establishing a systematic approach to selecting H will be crucial for enhancing the reliability of 2D-based analyses. Additional data across different flow configurations and regimes will be essential for establishing a comprehensive guideline in this regard.

Overall, the present study establishes 2D SHAP analysis as a reliable tool for wall-bounded turbulence investigations, offering a computationally efficient alternative to full 3D analyses and the possibility to use experimental data as provided by planar PIV techniques. Additionally, the 2D-uvw configuration, which includes spanwise velocity, produces results comparable to the 2D-uv case, suggesting the simplest experimental set-up can be chosen.

### CRediT authorship contribution statement

**Samuel Molina-Casino:** Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Conceptualization. **Andrés Cremades:** Writing – review & editing, Software, Methodology, Investigation, Data curation. **Sergio Hoyas:** Writing – review & editing, Resources, Funding acquisition, Data curation. **José I. Cardesa:** Writing – review & editing, Supervision. **François Chedevergne:** Writing – review & editing, Supervision. **Ricardo Vinuesa:** Writing – review & editing, Project administration, Methodology, Funding acquisition, Conceptualization.

### Code availability

The codes used for this work will be made available as open access at: https://github.com/KTH-FlowAI.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgments

The data has been obtained with support of grant PID2021-128676OB-I00 funded by MCIN/AEI/10.13039/ 501100011033 and by "ERDF A way of making Europe", by the European Union (SH). RV acknowledges the financial support from ERC grant no. 2021-CoG-101043998, DEEPCONTROL. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them.

## Appendix A. Gradient-SHAP methodology

For deep-learning models like the U-Net architecture used in this work, gradient-SHAP methods offer a computationally efficient way to estimate SHAP values using the so-called expected gradients (EG) [34], which build upon the concept of integrated gradients [42].

The core idea is to attribute the change in the model's output, relative to a baseline output, to the changes in its input features, relative to a baseline input. Let $f(x)$ be the output of the modified U-Net model (in our case, the MSE of the velocity prediction) for an input velocity field $x$. Let $x'$ be a baseline input (e.g., a zero-fluctuation field, as used in this study). The contribution of the $i$th feature (e.g., the velocity component at a specific grid point) to the difference $f(x)-f(x')$ is calculated by integrating the partial derivative of $f$ with respect to the $i$th feature along a straight-line path from $x'$ to $x$.

The EG value for the $i$th feature is thus approximated by:

$$EG_i(f, x) \approx (x_i - x'_i) \times \int_{\alpha=0}^{1} \frac{\partial f(x' + \alpha(x - x'))}{\partial x_i} d\alpha \quad (A.1)$$

In practice, this integral is often approximated numerically by summing the gradients at discrete steps along the path (effectively using linear interpolation). A common approximation for EG involves averaging the gradients obtained by perturbing the input features against multiple baseline samples, or by sampling points along the path from a single baseline to the input instance:

$$EG_i(f, x) \approx (x_i - x'_i) \times \sum_{k=1}^{m} \frac{\partial f(x' + \frac{k}{m}(x - x'))}{\partial x_i} \times \frac{1}{m} \quad (A.2)$$

where $m$ is the number of steps (or samples) taken along the path from the baseline $x'$ to the input $x$.

The obtained EG values are equivalent to the SHAP values $\phi_i$ from Eq. (1). SHAP values $\phi_u$, $\phi_v$, $\phi_w$ computed at each grid point in this study represent the individual contributions of the $u, v, w$ velocity components at that point to the model's prediction error (MSE), relative to the zero-fluctuation baseline. These values allow for the identification of spatially localized regions in the flow field that are most influential for the U-Net's predictive performance, as discussed in Section 2.2.

## Appendix B. Wall-normal and spanwise JPDFs

To ensure completeness, Figs. B.9 and B.10 present the joint probability functions of the wall-normal and spanwise velocity components for each structure type analyzed in this study.

### Data availability

The minimum representative data used in this study will be made available open access at: https://github.com/KTH-FlowAI. For the complete database, please contact the authors.
<figure><img src="image_10.png" /><figcaption>**Fig. B.9.** Joint probability density functions (JPDFs) of various structures as a function of the wall-normal distance and the inner-scaled wall-normal velocity fluctuation of the segmented domains. The results span a $y^+$ range between 20 and 125. (Top) JPDFs for SHAP-based structures in the 3D, 2D-*uv* and 2D-*uvw* datasets (left to right); (bottom) classical coherent structures: Q events, streaks and vortex clusters (left to right). To ensure consistency, the 3D results include only the slices corresponding to those used in the 2D cases.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. B.10.** Joint probability density functions (JPDFs) of various structures as a function of the wall-normal distance and the inner-scaled spanwise velocity fluctuation of the segmented domains. The results span a $y^+$ range between 20 and 125. (Top) JPDFs for SHAP-based structures in the 3D and 2D-*uvw* datasets (left to right); (bottom) classical coherent structures: Q events, streaks and vortex clusters (left to right). To ensure consistency, the 3D results include only the slices corresponding to those used in the 2D cases.</figcaption></figure>

## References

[1] O. Reynolds, An experimental investigation of the circumstances which determine whether the motion of water shall be direct or sinuous, and of the law of resistance in parallel channels, Philos. Trans. R. Soc. Lond. 174 (1883) 935–982, http://dx.doi.org/10.1098/rstl.1883.0029.

[2] E. Lagemann, S.L. Brunton, W. Schröder, C. Lagemann, Towards extending the aircraft flight envelope by mitigating transonic airfoil buffet, Nat. Commun. 15 (1) (2024) 5020, http://dx.doi.org/10.1038/s41467-024-49361-3.

[3] P. Varela, P. Suárez, F. Alcántara-Ávila, A. Miró, J. Rabault, B. Font, L.M. García-Cuevas, O. Lehmkuhl, R. Vinuesa, Deep reinforcement learning for flow control exploits different physics for increasing Reynolds number regimes, Actuators 11 (12) (2022) 359, http://dx.doi.org/10.3390/act11120359.

[4] J. Kim, P. Moin, R. Moser, Turbulence statistics in fully developed channel flow at low Reynolds number, J. Fluid Mech. 177 (1987) 133–166, http://dx.doi.org/10.1017/S0022112087000892.

[5] S. Hoyas, M. Oberlack, F. Alcántara-Ávila, S.V. Kraheberger, J. Laux, Wall turbulence at high friction Reynolds numbers, Phys. Rev. Fluids 7 (1) (2022) 014602, http://dx.doi.org/10.1103/PhysRevFluids.7.014602.
[6] R. Vinuesa, P. Negi, M. Atzori, A. Hanifi, D. Henningson, P. Schlatter, Turbulent boundary layers around wing sections up to Rec=1,000,000, Int. J. Heat Fluid Flow 72 (2018) 86-99, http://dx.doi.org/10.1016/j.ijheatfluidflow.2018.04.017.

[7] F. Scarano, M.C. Jacob, X. Carbonneau, E.R. Gowree, On the turbulent boundary layer over a flat plate at moderate Reynolds numbers, Phys. Fluids 34 (11) (2022) 115150, http://dx.doi.org/10.1063/5.0124498.

[8] D. Squire, C. Morrill-Winter, N. Hutchins, I. Marusic, M. Schultz, J. Klewicki, Smooth- and rough-wall boundary layer structure from high spatial range particle image velocimetry, Phys. Rev. Fluids 1 (6) (2016) 064402, http://dx.doi.org/10.1103/PhysRevFluids.1.064402.

[9] C. Tropea, A. Yarin, J. Foss, Chapter 12: Wall-bounded Flows. Handbook of Experimental Fluid Mechanics, in: Handbook of Experimental Fluid Mechanics, Springer, 2007.

[10] N. Hutchins, T. Nickels, I. Marusic, M. Chong, Hot-wire spatial resolution issues in wall-bounded turbulence, J. Fluid Mech. 635 (2009) 103-136, http://dx.doi.org/10.1017/S0022112009007721.

[11] J. Jiménez, Coherent structures in wall-bounded turbulence, J. Fluid Mech. 842 (2018) P1, http://dx.doi.org/10.1017/jfm.2018.144.

[12] A. Lozano-Durán, O. Flores, J. Jiménez, The three-dimensional structure of momentum transfer in turbulent channels, J. Fluid Mech. 694 (2012) 100-130, http://dx.doi.org/10.1017/jfm.2011.524.

[13] S. Kline, W. Reynolds, F. Schraub, P. Runstadler, The structure of turbulent boundary layers, J. Fluid Mech. 30 (4) (1967) 741-773, http://dx.doi.org/10.1017/S0022112067001740.

[14] M. Chong, A. Perry, B. Cantwell, A general classification of three-dimensional flow fields, Phys. Fluids A Fluid Dyn. 2 (5) (1990) 765-777, http://dx.doi.org/10.1063/1.857730.

[15] L. Shapley, A value for N-person games, in: H.W. Kuhn, A.W. Tucker (Eds.), in: Contributions To the Theory of Games, vol. II, Princeton University Press, 1953, pp. 307-318, http://dx.doi.org/10.1515/9781400881970-018.

[16] S. Lundberg, S. Lee, A Unified Approach to Interpreting Model Predictions, in: Advances in Neural Information Processing Systems, vol. 30, Curran Associates, Inc., 2017.

[17] A. Cremades, S. Hoyas, R. Vinuesa, Classically studied coherent structures only paint a partial picture of wall-bounded turbulence, 2024, http://dx.doi.org/10.48550/arXiv.2410.23189, arXiv:2410.23189.

[18] A. Cremades, S. Hoyas, R. Deshpande, P. Quintero, M. Lellep, W. Lee, J. Monty, N. Hutchins, M. Linkmann, I. Marusic, R. Vinuesa, Identifying regions of importance in wall-bounded turbulence through explainable deep learning, Nat. Commun. 15 (1) (2024) 3864, http://dx.doi.org/10.1038/s41467-024-47954-6.

[19] P. Suárez, F. Álcantara-Ávila, A. Miró, J. Rabault, B. Font, O. Lehmkuhl, R. Vinuesa, Active flow controls for drag reduction through multi-agent reinforcement learning on a turbulent cylinder at $SRe_D=3900$, 2025, http://dx.doi.org/10.48550/arXiv.2405.17655, arXiv:2405.17655.

[20] M. Beneitez, A. Cremades, L. Guastoni, R. Vinuesa, Improving turbulence control through explainable deep learning, 2025, http://dx.doi.org/10.48550/arXiv.2504.02354, arXiv:2504.02354.

[21] X. He, J. Tan, G. Rigas, M. Vahdati, On the explainability of machine-learning-assisted turbulence modeling for transonic flows, Int. J. Heat Fluid Flow 97 (2022) 109038, http://dx.doi.org/10.1016/j.ijheatfluidflow.2022.109038.

[22] A. Cremades, S. Hoyas, R. Vinuesa, Additive-feature-attribution methods: A review on explainable artificial intelligence for fluid dynamics and heat transfer, Int. J. Heat Fluid Flow 112 (2025) 109662, http://dx.doi.org/10.1016/j.ijheatfluidflow.2024.109662.

[23] F. Lluesma-Rodríguez, F. Álcantara-Ávila, M.J. Pérez-Quiles, S. Hoyas, A code for simulating heat transfer in turbulent channel flow, Mathematics 9 (7) (2021) 756, http://dx.doi.org/10.3390/math9070756.

[24] A. Lozano-Durán, J. Jiménez, Effect of the computational domain on direct simulations of turbulent channels up to $Re_x = 4200$, Phys. Fluids 26 (1) (2014) 011702, http://dx.doi.org/10.1063/1.4862918.

[25] A.N. Kolmogorov, The local structure of turbulence in incompressible viscous fluid for very large Reynolds numbers, Proc. R. Soc. Lond. Ser. A Math. Phys. Sci. 434 (1890) (1941) 9-13, http://dx.doi.org/10.1098/rspa.1991.0075.

[26] R.J. Adrian, Particle-imaging techniques for experimental fluid mechanics, Annu. Rev. Fluid Mech. 23 (1991) 261-304, http://dx.doi.org/10.1146/annurev.fl.23.010191.001401.

[27] M. Raffel, C.E. Willert, S.T. Wereley, J. Kompenhans, PIV recording techniques, in: M. Raffel, C.E. Willert, S.T. Wereley, J. Kompenhans (Eds.), Particle Image Velocimetry: A Practical Guide, Springer, Berlin, Heidelberg, 2007, pp. 97-121, http://dx.doi.org/10.1007/978-3-540-72308-0-4.

[28] G. Elsinga, F. Scarano, B. Wieneke, B. van Oudheusden, Tomographic particle image velocimetry, Exp. Fluids 41 (6) (2006) 933-947, http://dx.doi.org/10.1007/s00348-006-0212-z.

[29] C.J. Kähler, S. Scharnowski, C. Cierpka, Highly resolved experimental results of the separated flow in a channel with streamwise periodic constrictions, J. Fluid Mech. 796 (2016) 257-284, http://dx.doi.org/10.1017/jfm.2016.250.

[30] A. Schröder, D. Schanz, R. Geisler, S. Gesemann, C. Willert, Near-wall turbulence characterization using 4D-PTV shake-the-box, in: Measurement Science and Technology, IOP Science, Santa Barbara, California, 2015.

[31] R. Örlü, R. Vinuesa, Instantaneous wall-shear-stress measurements: Advances and application to near-wall extreme events, Meas. Sci. Technol. 31 (11) (2020) 112001, http://dx.doi.org/10.1088/1361-6501/aba06f.

[32] M.T. Ribeiro, S. Singh, C. Guestrin, "Why should I trust you?": explaining the predictions of any classifier, in: Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD '16, Association for Computing Machinery, New York, NY, USA, 2016, pp. 1135-1144, http://dx.doi.org/10.1145/2939672.2939778.

[33] S.M. Lundberg, G.G. Erion, S.-I. Lee, Consistent individualized feature attribution for tree ensembles, 2018, http://dx.doi.org/10.48550/arXiv.1802.03888, arXiv: 1802.03888.

[34] G. Erion, J.D. Janizek, P. Sturmfels, S.M. Lundberg, S.-I. Lee, Improving performance of deep learning models with axiomatic attribution priors and expected gradients, Nat. Mach. Intell. 3 (7) (2021) 620-631, http://dx.doi.org/10.1038/s42256-021-00343-w.

[35] O. Ronneberger, P. Fischer, T. Brox, U-net: convolutional networks for biomedical image segmentation, in: N. Navab, J. Hornegger, W. Wells, A. Frangi (Eds.), Medical Image Computing and Computer-Assisted Intervention – MICCAI 2015, Springer International Publishing, Cham, 2015, pp. 234-241, http://dx.doi.org/10.1007/978-3-319-24574-4_28.

[36] Y. LeCun, B. Boser, J. Denker, D. Henderson, R. Howard, W. Hubbard, L. Jackel, Handwritten Digit Recognition with a Back-Propagation Network, in: Advances in Neural Information Processing Systems, vol. 2, Morgan-Kaufmann, 1989.

[37] S. Gandía-Barberá, A. Cremades, R. Vinuesa, S. Hoyas, M.J. Pérez-Quiles, Sequential and parallel algorithms to compute turbulent coherent structures, Mathematics 12 (21) (2024) 3325, http://dx.doi.org/10.3390/math12213325.

[38] J.C. del Álamo, J. Jiménez, P. Zandonade, R.D. Moser, Self-similar Vortex Clusters in the Turbulent Logarithmic Region, J. Fluid Mech. 561 (2006) 329-358, http://dx.doi.org/10.1017/S0022112006000814.

[39] J.M. Wallace, H. Eckelmann, R.S. Brodkey, The wall region in Turbulent Shear flow, J. Fluid Mech. 54 (1) (1972) 39-48, http://dx.doi.org/10.1017/S0022112072000515.

[40] W. Schoppa, F. Hussain, Coherent structure generation in near-wall turbulence, J. Fluid Mech. 453 (2002) 57-108, http://dx.doi.org/10.1017/S002211200100667X.

[41] K. Osawa, J. Jiménez, Causal features in turbulent channel flow, J. Fluid Mech. 1000 (2024) A4, http://dx.doi.org/10.1017/jfm.2024.838.

[42] M. Sundararajan, A. Taly, Q. Yan, Axiomatic attribution for deep networks, in: Proceedings of the 34th International Conference on Machine Learning, PMLR, 2017, pp. 3319-3328.