2.1. Generation of high-fidelity database

The observational datasets for DA were obtained through synchronized time-resolved PIV/PLIF measurements. Additionally, delayed detached eddy simulations (DDESs) were conducted to quantitatively benchmark the reconstruction capabilities of the DA procedure. The experiment set-up, depicted in figure 1, involves a water tunnel submerged within a large glass tank. To ensure a stable mainstream flow, a honeycomb and a contraction section were positioned upstream of the 600-mm-long test section, which features a cross-section of 80 mm (width) × 50 mm (height). The free stream turbulence level is roughly 5 %. During the experiments, the cross-flow velocity at the inlet of the test section was maintained at ${U_0} = 0.38\;\textrm{m}\;{\textrm{s}^{ - 1}}$, corresponding to a Reynolds number of 150 000 based on the distance between the leading edge of the test plate and the coolant injection holes. Two different velocity ratios ($VR = {U_j}/{U_0} = 0.4$ and $1.2$, where ${U_j}$ is the mean jet velocity) were examined by adjusting the jet flow rate. These ratios represent both attached and detached flow regimes, resulting in jet Reynolds numbers (based on the tube diameter) of 1200 and 3600, respectively. Both the cross-flow and jet flow Reynolds number falls in the scope of the typical engine condition (Ye et al. Reference Ye, Liu, Zhu and Luo2019). A constant-head tank was employed to maintain a consistent water level, while a jet chamber was used to stabilize the jet before introduction into the mainstream. Multi-plane measurements ($z/D = 0,\;0.3$ and $0.5$) were conducted to comprehend the three-dimensional flow features in JICF. The jet originated from a circular tube with dimensions: diameter ($D$) of 8 mm, length ($L$) of 40 mm, and an inclination angle ($\alpha $) of 30°. This configuration was selected for its representativeness of film cooling applications, highlighting the highly non-uniform and secondary flow characteristics within the hole.

Figure 1. Experimental set-ups for the simultaneous PIV/PLIF measurements.

A continuous 532-nm laser illuminated PIV particles and served as the excitation source for rhodamine solution in PLIF measurements. Synchronized measurements were performed using two high-speed cameras (Dimax HS, PCO) to capture instantaneous velocity and concentration field. The velocity fields were determined using the adaptive cross-correlation PIV algorithm with 32 × 32 pixels interrogation window and 50 % overlap, resulting in a grid spacing of $0.05D$. Rhodamine dye, with a concentration of ${C_0} = \; 0.05\;\textrm{mg}\;{\textrm{L}^{ - 1}}$, was employed as the jet fluid source. To eliminate the excitation source, a long-wave-pass filter (570 nm) was positioned before the camera lens. The sampling frequency was 1000 Hz and PLIF images had a spatial resolution of $0.02D$.

The computational domain mirrors the experimental set-up depicted in figure 1, with the mainstream inlet boundary located $4.0D$ upstream of the jet exit centre. The numerical simulation incorporates the water chamber, consistent with established film cooling practices, to resolve the recirculation occurring inside the hole before encountering the cross-flow (Bodart, Coletti & Eaton Reference Bodart, Coletti and Eaton2013). In figure 2(c), time-averaged velocity profiles extracted at $x/D ={-} 1.5$ in three measurement planes are shown. Notably, the velocity profiles exhibit good agreement, supporting the notion of treating the mainstream inlet velocity as uniformly distributed along the spanwise direction. To replicate the experimental environment in the numerical simulation, the inflow statistics are derived from time-resolved PIV measurements. The spectral synthesizer, originally proposed by Kraichnan (Reference Kraichnan1970) and modified by Smirnov, Shi & Celik (Reference Smirnov, Shi and Celik2001), is employed to impose the desired flow statistics. Since PIV data for the Reynolds stresses out of the $xy$-plane are unavailable, we resemble these stresses based on the distribution from the high-fidelity simulation of a turbulent inclined JICF by Bodart et al. (Reference Bodart, Coletti and Eaton2013) with a mean turbulence intensity of 5 %. These inflow statistics, extracted at $x/D ={-} 1.5$ from PIV measurements, are applied to the DDES with the mainstream inlet boundary located $4.0D$ upstream of the jet exit centre (figure 2b). The cross-flow boundary layer characteristics $1.5D$ upstream of the jet orifice are then verified using the DDES statistics for velocity and Reynolds stress behaviour. Despite some discrepancies in turbulence statistics, the mean velocity agrees well with the experiments, as shown in figure 2(d–g). For the Spalart–Allmaras (SA) model, we impose the turbulent mean velocity profile and the actual eddy viscosity determined from PIV, as depicted in figure 2(c). The eddy viscosity ratio is determined by the linear relationship of the Boussinesq hypothesis, as suggested by Weatheritt et al. (Reference Weatheritt, Zhao, Sandberg, Mizukami and Tanimoto2020). At the chamber inflow, velocity and concentration are imposed without any superimposed turbulence. Figure 2(b) depicts the numerical meshes of the JICF configuration. The computational grids are block-structured and locally refined with 100 inflation boundary layers in the region of interest ($- 1.5 < x/D < 5$, $0 < y/D < 2$ and $- 1.5 < z/D < 1.5$) to better resolve the gradients near the shear layer and solid surfaces. The grid resolution ensures that the first cell centre is located at ${y^ + } < 1.5$ for all walls, except for the upper wall where a slip boundary condition is applied. This set-up adequately captures the blockage effect (Li, He & Liu Reference Li, He and Liu2023). An in-house dynamic DDES code (He, Liu & Yavuzkurt Reference He, Liu and Yavuzkurt2017) is employed for the simulation on a refined mesh comprising either 9.1 million ($VR = 0.4$) or 10.3 million ($VR = 1.2$) elements. The time step size, $\Delta t_{DDES}^ + = \Delta {t_{DDES}}{U_0}/D = 0.01$, ensures that the maximum Courant–Friedrichs–Lewy number is always less than 1. A statistically steady state is reached after 100 time units, where $T = D/{U_0}$. Turbulent quantities are then averaged over a period of 200T, which leads to an acceptable level of convergence for this fully non-homogeneous configuration in which no spatial averaging can be employed to accelerate the convergence.

Figure 2. (a) Schematic of the computational domain, (b) computational mesh with every 10 grid nodes displayed, (c) profiles of the streamwise velocity and the eddy viscosity ratio ${v_t}/v$ determined through PIV measurements and (d–g) comparison of velocity and turbulence statistics extracted from $x/D ={-} 1.5$.

2.2. Model validation

The mean turbulent flow, as determined by DDES and SA models, was validated against PIV/PLIF measurements in the mid-plane (figures 3 and 4). The scalar transport in SA model employs a straightforward, isotropic GDH, where the turbulent diffusivity is prescribed with a fixed turbulent Schmidt number, denoted as $S{c_t} = 0.85$ (Kays Reference Kays1994). Regarding the $x$-direction velocity, all profiles exhibit favourable agreement between experimental data and DDES results. However, the SA model fails to accurately capture jet spreading and diffusion, leading to noticeable discrepancies. Two significant distinctions between the default SA model and PIV/DDES results are observed from the contours in figure 3(a,c). First, the pattern of the jet core, as indicated by the $x$-direction velocity, varies noticeably. The experimental and DDES results reveal a broader jet core that extends deeper into the mainstream, whereas the default SA model predicts a narrower core squeezed into the upstream region of the jet, leading to increased $y$-direction velocity. Second, the $x$-direction velocity downstream of the injection, simulated using the default SA model, is notably lower compared with the reference data. This discrepancy is particularly pronounced in the velocity profiles in figure 4(a–d). There are two local maxima for each velocity ratio, representing the upstream and downstream shear layers, respectively. The SA model consistently underestimates $x$-direction velocity for both shear layers across all velocity ratios. Additionally, it is important to emphasize that the cross-flow incoming boundary layer is thicker than the other two, indicating an overestimation of eddy viscosity in front of the jet. This discrepancy stems from the deficiency of the SA model in handling strong adverse pressure gradient (APG) flows. The SA model assumes equilibrium conditions in the boundary layer, but under strong APG conditions, this assumption fails. The destruction term in the SA model decays turbulence too quickly, leading to overestimated eddy viscosity and a thicker boundary layer.

Figure 3. Contours of the PIV/PLIF, SA and DDES results at the mid-plane. (a,b) $VR = 0.4$, (c,d) $VR = 1.2$.

Figure 4. Comparison of the PIV/PLIF, SA and DDES results at four stations.

Regarding scalar distribution, the streamwise decay of the scalar field is accurately predicted by DDES, and the extent of the jet's shear layer aligns well with the PLIF data. However, as illustrated in figures 3(b,d) and 4(e–h), the deviation trends in the scalar field predicted by the SA model are entirely different from those in the velocity fields. The isotropic model broadly underestimates spreading, resulting in a longer downstream transport of scalar towards the wall after injection and a slower decay of the jet. Thus, the jet scalar penetrates deeper into the mainstream for SA than for DDES and experimental results. This penetration is influenced by both the convective and turbulent scalar flux transport. The convective transport is affected by the Reynolds stress, while the turbulent transport is affected by the turbulent scalar diffusion. In the default SA model, the linear eddy viscosity model (LEVM) leads to limited dissipation, and the gradient diffusion hypothesis with uniform Schmidt number results in insufficient turbulent diffusion. This suggests the necessity for careful consideration of both Reynolds stress and turbulent scalar flux in JICF for accurate scalar mixing.

3.1. Identification of counter-gradient transportation

Figure 5 presents the evolution of concentration and vorticity fields at two $VR\textrm{s}$. Each case shows two snapshots separated by a dimensionless time interval equivalent to approximately 0.01 in terms of $T = D/{U_0}$. To facilitate comparison between the vorticity and concentration fields, the jet boundary is outlined in the vorticity field using a threshold of $C/{C_0} = 0.8$. These vivid results aid in intuitively understanding the evolution of concentration and vorticity fields and their dynamic relationship. The snapshots reveal mixing from the Kelvin–Helmholtz (K–H) mode and the downstream growth of flow structures. Figure 5(b,d) displays clockwise (blue) and anticlockwise (red) K–H vortices along the jet leading edge and from the separated wake region in the tube. Anticlockwise vortices with positive vorticity originate from the tube boundary layer, while clockwise vortices shed from the vorticity feed associated with the horseshoe vortex. The interaction between the anticlockwise vortex and the upstream horseshoe vortex leads to noticeable oscillations during vortex shedding.

Figure 5. Simultaneous snapshots of (a,c) concentration and (b,d) vorticity field with a dimensionless time interval of 0.01 at two $VR\textrm{s}$. (a,b) VR = 0.4, (c,d) VR = 1.2.

The mean velocity gradient plays a dominant role in the shear layer dynamics. In the low-velocity-ratio scenario (figure 5a,b), the positive mean velocity gradient in the upper shear layer promotes the formation of clockwise vortices. Anticlockwise vortices, however, undergo vorticity cancellation and decay as they are convected downstream of the jet orifice. Observations at locations A and B reveal that the spatiotemporal evolution of concentration structures is influenced by the corresponding vortical structures: local clockwise vortical structures induce a clockwise curl in the local concentration field. As the velocity ratio increases to 1.2, the inclined jet fully detaches from the wall. This intensifies vortical structures, leading to a denser and stronger mixing process and more fragmented concentration structures compared with at $VR = 0.4$ (figure 5c,d). In this higher velocity ratio scenario, the negative mean velocity gradient in the upper shear layer promotes the formation of counterclockwise vortices. Despite this shift, a similar correlation persists between the spatiotemporal evolution of concentration and vorticity fields at locations A and B: the anticlockwise curl of local concentration field arises from the influence of an anticlockwise vortical structure.

The synchronous correlation between concentration and vorticity fields offers valuable insights into the turbulent scalar flux distribution in the mixing process of JICF. In the PIV/PLIF measurement, turbulent scalar flux comprises two components: the streamwise term ($\overline {{u^{\prime}_1}c^{\prime}} /{U_0}C$) and the wall-normal term ($\overline {{u_{\textrm{2}}^{\prime}}c^{\prime}} /{U_0}C$). Figure 6 illustrates the contours of normalized (a,c) streamwise and (b,d) wall-normal components of turbulent scalar flux at two velocity ratios. At $VR = 0.4$, the streamwise turbulent scalar flux shows significant magnitude in the upstream shear layer with a negative value, whereas it turns positive at $VR = 1.2$. As for the wall-normal turbulent scalar flux $\overline {u_{\textrm{2}}^{\prime}c^{\prime}} $, the distribution pattern is quite different. In both examined velocity ratios, the wall-normal term exhibits intensity in the windward shear layer with positive values but is considerable weaker in the downstream shear layer.

Figure 6. (a,c) Streamwise and (b,d) wall-normal turbulent scalar flux determined by synchronized PIV/PLIF. (a,b) $VR = 0.4$, (c,d) $VR = 1.2$.

The distribution of turbulent scalar flux in the jet shear layer is interpretable. When $VR = 0.4$, the streamwise component of jet velocity is lower than that of the mainstream velocity, yet the jet concentration is higher than that of the cross-flow. During turbulent fluctuation, a location in the windward shear layer is dominated by the jet and the cross-flow alternately due to K–H instability, as indicated in figure 5(a,b). When the cross-flow dominates at a certain time, the streamwise component of velocity is likely to be larger than the time-averaged streamwise velocity (${u^{\prime}_\textrm{1}} > 0$), while the local concentration is lower than the time-averaged concentration ($c^{\prime} < 0$). Similarly, when the jet dominates, $c^{\prime} > 0$ and ${u^{\prime}_\textrm{1}}$ is likely to be negative. Therefore, the streamwise turbulent scalar flux $\overline {{u_{\textrm{1}}^{\prime}}c^{\prime}} $ is negative in the windward shear layer. As for the wall-normal component of turbulent scalar flux $\overline {{u_{\textrm{2}}^{\prime}}c^{\prime}} $, the measurement results shown in figure 6 can also be interpreted in the same way. Moreover, the wide region and relatively large magnitude of positive $\overline {{u_{\textrm{2}}^{\prime}}c^{\prime}} $ in the windward shear layer for both cases reveal a significant relationship between the wall-normal velocity and concentration field in an inclined JICF.

Regarding the downstream shear layer, $\overline {{u_{\textrm{1}}^{\prime}}c^{\prime}} $ is also positive, as depicted in figure 6(a). It is logical to explain the case of $VR = 1.2$ by the alternating dominance of jet and cross-flow. However, the scenario seems counterintuitive for $VR = 0.4$. When considering the wall effect, even at low velocity ratio, the jet velocity in the downstream shear layer exceeds the cross-flow velocity, as shown in figure 3(a). Consequently, the streamwise turbulent scalar flux remains positive even at low $VR$. In the wall-normal direction, negligible correlation is observed within the downstream shear layer at $VR = 0.4$, and only weak negative correlations are observed at $VR = 1.2$.

Most RANS models reduce the turbulent scalar flux to a diffusion term of the form:

(3.1)\begin{equation}\overline {{u^{\prime}_i}c^{\prime}} ={-} \frac{{{\upsilon _t}}}{{S{c_t}}}{\partial _i}C.\end{equation}

In this scenario, turbulent diffusivity primarily relies on the turbulent viscosity (${\upsilon _t}$) in the momentum equation using the Reynolds analogy and on the turbulent Schmidt number ($S{c_t}$) using the GDH. In most turbulent regions, the mean concentration gradient and the mean scalar flux transport by turbulence retain the same sign across the flow. This indicates that scalar transport consistently occurs from regions of higher mean scalar concentration to those of lower concentration. However, the alignment between the two vectors, $\overline {{u^{\prime}_i}c^{\prime}} $ and ${\partial _i}C$, may be notably skewed. This can be identified by examining the sign of the product between each component of $\overline {{u^{\prime}_i}c^{\prime}} $ and ${\partial _i}C$, as depicted in figure 7, revealing instances of counter-gradient diffusion. Notably, two prominent features stand out: (a) in the upstream shear layer, a positive correlation between $\overline {{u^{\prime}_i}c^{\prime}} $ and ${\partial _1}C$ is observed for both $VR\textrm{s}$. Specifically, for $VR = 1.2$, a sizable positive region immediately follows the jet release, while for $VR = 0.4$, the correlation shifts to positive after $2D$ downstream of the jet orifice; (b) in the downstream shear layer near the wall, a positive correlation between $\overline {u_{2}^{\prime}c^{\prime}} $ and ${\partial _2}C$ is observed for both $VR\textrm{s}$, with a notably enhanced positive correlation for $VR = 1.2$. This correlation can be readily deduced by considering the distribution of ${\partial _i}C$, as evident in figure 3(b,d), in conjunction with the distribution of turbulent scalar flux illustrated in figure 6.

Figure 7. Identification of counter-gradient transport regions in (a,c) streamwise and (b,d) wall-normal directions determined by synchronized PIV/PLIF. (a,b) $VR = 0.4$, (c,d) $VR = 1.2$. Terms are non-dimensionalized by ${U_0}$ and D.

3.2. Local and non-local contribution

Counter-gradient transport is commonly observed in combustion problems (Pfadler et al. Reference Pfadler, Kerl, Beyrau, Leipertz, Sadiki, Scheuerlein and Dinkelacker2009), where it is generally understood to arise from thermal dilatation due to chemical reactions. However, the previous subsection reveals CGT in different flow regions without any chemical reactions, consistent with reports from authors such as Sau & Mahesh (Reference Sau and Mahesh2008), Schreivogel et al. (Reference Schreivogel, Abram, Fond, Straußwald, Beyrau and Pfitzner2016) and Milani, Ling & Eaton (Reference Milani, Ling and Eaton2020). This section focuses on clarifying the existence of CGT in these regions, identifying the dominant flow structures responsible and pinpointing their locations.

To elucidate the presence of CGT in various regions, it is useful to examine the scalar flux transport equation, presented in (3.2). Term ${{\boldsymbol{\mathcal P}}_{\boldsymbol{i}}}$ is the production term, showing how mean velocity and scalar gradients locally generate $\overline {{u^{\prime}_i}c^{\prime}} $. Term ${{\boldsymbol{\mathcal D}}_{\boldsymbol{i}}}$ represents redistribution. The source term, ${\boldsymbol{\varphi }_{\boldsymbol{i}}}$, arises from the pressure-scalar correlation, which is inherently non-local because pressure fluctuations are governed by an elliptic Poisson equation. Term ${{\boldsymbol{\mathcal E}}_{\boldsymbol{i}}}$ represents viscous destruction associated with small scales.

(3.2)\begin{equation}{D_t}\overline {{u^{\prime}_i}c^{\prime}} = \overbrace{{ - \overline {{u^{\prime}_j}c^{\prime}} {\partial _j}{U_i} - \overline {{u^{\prime}_j}{u^{\prime}_i}} {\partial _j}C}}^{{{{\boldsymbol{\mathcal P}}_{\boldsymbol{i}}}}} + {{\boldsymbol{\mathcal D}}_{\boldsymbol{i}}} + {\boldsymbol{\varphi }_{\boldsymbol{i}}} - {{\boldsymbol{\mathcal E}}_{\boldsymbol{i}}}.\end{equation}

Assuming that the experimental and numerical results have reached statistical steady-state values and considering physical symmetry,

(3.3)\begin{equation}\left. {\begin{array}{c@{}} {{\partial_3}\phi \equiv 0,}\\ {{\partial_j}{U_3} \equiv 0,} \end{array}} \right\}\end{equation}

(3.4)\begin{equation}{{\boldsymbol{\mathcal P}}_{\boldsymbol{i}}} = \left\{ {\begin{array}{@{}c} { - \overline {{u^{\prime}_1}c^{\prime}} {\partial_1}{U_1} - \overline {{u^{\prime}_1}{u^{\prime}_1}} {\partial_1}C - \overline {{u^{\prime}_2}c^{\prime}} {\partial_2}{U_1} - \overline {{u^{\prime}_2}{u^{\prime}_1}} {\partial_2}C,}\\ { - \overline {{u^{\prime}_1}c^{\prime}} {\partial_1}{U_2} - \overline {{u^{\prime}_1}{u^{\prime}_2}} {\partial_1}C - \overline {{u^{\prime}_2}c^{\prime}} {\partial_2}{U_2} - \overline {{u^{\prime}_2}{u^{\prime}_2}} {\partial_2}C,}\\ { - \overline {{u^{\prime}_1}{u^{\prime}_3}} {\partial_1}C - \overline {{u^{\prime}_2}{u^{\prime}_3}} {\partial_2}C.} \end{array}} \right.\end{equation}

Next, considering the scalar field analogously, the production term ${{\boldsymbol{\mathcal P}}_{\boldsymbol{c}}}$ for scalar fluctuation variance can be expressed as

(3.5)\begin{equation}{{\boldsymbol{\mathcal P}}_{\boldsymbol{c}}} ={-} \overline {{u_{\textrm{1}}^{\prime}}c^{\prime}} {\partial _1}C - \overline {{u_{\textrm{2}}^{\prime}}c^{\prime}} {\partial _2}C.\end{equation}

The production term for scalar fluctuations reflects the transfer of scalar intensity from the mean field to the fluctuating field, akin to the turbulent kinetic energy cascade. A negative ${{\boldsymbol{\mathcal P}}_{\boldsymbol{c}}}$ signifies that scalar intensity is transferred from smaller scales to mean flow, indicating non-local behaviour.

Figure 8 shows the production terms for turbulent scalar flux and scalar fluctuation variance, as determined by synchronized PIV/PLIF. Building on the insights from figures 6 and 7, figure 8 identifies two distinct types of CGT regions. The first type, evident in the windward shear layer (figure 7a,c), is associated with the streamwise component of the turbulent scalar flux. Here, the production term aligns with the sign of the streamwise turbulent scalar flux, indicating primarily local CGT behaviour. This local CGT is mainly driven by the cross-gradient transport of turbulent scalar flux. In the grey regions depicted in figure 9, where ${\partial _2}C < 0$ and ${\partial _2}U > 0$, turbulent eddies bringing fluid from above (${u^{\prime}_\textrm{2}} < 0$) tend to carry a weaker concentration ($c^{\prime} < 0$) and higher streamwise velocity (${u^{\prime}_\textrm{1}} > 0$). Similarly, eddies that bring fluid from below (${u^{\prime}_\textrm{2}} > 0$) link $c^{\prime} < 0$ with ${u^{\prime}_\textrm{1}} < 0$. This pattern is consistent in jets with higher velocity ratios, where vertical turbulent transport produces positive $\overline {u_{\textrm{1}}^{\prime}c^{\prime}} $ in that region. These behaviours align with the alternating dominance of jet and cross-flow correlations discussed in § 3.1. To further elucidate, figure 9 depicts various components of the $\overline {u_{\textrm{1}}^{\prime}c^{\prime}}$ production term. Note that the GDH considers only the mean concentration gradient in the $x$-direction on $\overline {u_{\textrm{1}}^{\prime}c^{\prime}}$, reflected in the production term $- \overline {u_{\textrm{1}}^{\prime}{u^{\prime}_\textrm{1}}} {\partial _1}C$. However, in regions exhibiting CGT, the production term is dominated by vertical gradient components, specifically $- \overline {u_{\textrm{2}}^{\prime}c^{\prime}} {\partial _2}{U_1}$ and $- \overline {u_{\textrm{2}}^{\prime}{u^{\prime}_\textrm{1}}} {\partial _2}C$. These terms, shown as negative in figure 9(a) and positive in figure 9(b), correspond with the observed sign of $\overline {u_{\textrm{1}}^{\prime}c^{\prime}}$ in the grey regions.

Figure 8. Production terms of (a,b,d,e) turbulent scalar flux and (c, f) scalar fluctuation variance determined by synchronized PIV/PLIF. (a–c) $VR = 0.4$, (d–f) $VR = 1.2$. Terms are non-dimensionalized by ${U_0}$ and D.

Figure 9. Vertical profiles of different component of the $\overline {u_{\textrm{1}}^{\prime}c^{\prime}}$ production term, mean streamwise velocity and concentration. Profiles are extracted (a) at $x/D = 2.8$ and $z/D = 0$ for $VR = 0.4$ and (b) at $x/D = 0$ and $z/D = 0$ for $VR = 1.2$ using synchronized PIV/PLIF. The grey band represents the region of CGT. Terms are non-dimensionalized by ${U_0}$ and D.

To identify the flow structure responsible for the local CGT behaviour, the co-spectra ${{{\mathcal F}}_{{u^{\prime}_i}c^{\prime}}}$ of the turbulence transport term are presented. Figure 11(a,b) shows the first momentum of co-spectra $f{{{\mathcal F}}_{{u^{\prime}_i}c^{\prime}}}$, normalized by the sign of ${\partial _i}C$ to clearly depict the CGT contribution (Sen et al. Reference Sen, Chuangxin, Bing and Yingzheng2023). The Strouhal number, $St$, is defined as $St = f\kern0.07em D/{U_0}$. Figure 11(a,b) illustrates the co-spectrum for the point where CGT occurs. It is evident that ${u^{\prime}_2}c^{\prime}$ consistently portrays gradient transport behaviour across all investigated frequencies, with the large-scale structure of the windward shear-layer playing a dominant role. This large-scale structure, identifiable by its characteristic Strouhal number – corresponding to shear layer vortices at high VR and hairpin vortices at low VR – has been extensively documented in the literature (Mahesh Reference Mahesh2013). Regarding the streamwise component, regions with opposite sign effectively define the contributions from the low- and high-frequency components of the motion. The low-frequency components, related to the largest eddies, are inferred to be the main mechanism for the CGT of scalar flux. The change in sign occurs at approximately $St = 1.3$, corresponding to a length scale of $l/D = {U_c}/f\kern0.07em D \approx 0.3$, where ${U_c}$ is the convection velocity of the vortices (Sen et al. Reference Sen, Chuangxin, Bing and Yingzheng2023). Similarly, the length scale determined from the characteristic Strouhal number is in the range of $0.5\unicode{x2013} 1.0$. A local transport model requires that the characteristic scale of the transport mechanism is comparable to or smaller than the distance over which the mean gradient of the transported property changes appreciably (Hamba Reference Hamba2022). This suggests that the CGT behaviour aligns with the local turbulent effects. Evidence that cross-gradient effects are responsible for the streamwise CGT behaviour includes the energetic scales for both ${u^{\prime}_1}c^{\prime}$ and ${u^{\prime}_2}c^{\prime}$ being in the same frequency band, and the vertical scalar flux showing a strong correlation. This corroborates that vertical large-scale motion drives the local CGT behaviour. As a result, the scalar flux transport equation redistributes energy from the vertical to the streamwise component, with the CGT behaviour in the streamwise turbulence scalar flux governed by vertical turbulence motion.

A second region where counter gradient transport is present is near the wall, right after injection. In this case, a negative production term is observed for both streamwise and vertical turbulent scalar flux, and for scalar fluctuation variance. This indicates that local effects act as a sink and, therefore, favours a negative value of $\overline {u_{\textrm{2}}^{\prime}c^{\prime}}$. Since the resulting scalar flux is positive, other effects are overwhelming the production. To help explain this, figure 10 shows a complete budget of $\overline {{u^{\prime}_2}c^{\prime}} $, with vertical profiles of all terms of (3.2). The viscous destruction term is omitted since a first-order isotropic tensor does not exist, making the dissipation rate zero (Durbin & Reif Reference Durbin and Reif2010). Since we cannot directly obtain all budgets from PIV/PLIF measurements, only DDES budgets are depicted, with experimental results of the production term shown for comparison. The DDES results closely replicate the experimental data. It can be clearly seen that, since the resulting scalar flux is positive, the pressure-scalar correlation effects are overwhelming the production. This is especially true close to the wall, where pressure fluctuation is dominant (Durbin Reference Durbin2018). This suggests that non-local turbulent effects, primarily through fluctuating pressure, are responsible for generating a positive correlation between ${u^{\prime}_\textrm{2}}$ and $c^{\prime}$ in this region.

Figure 10. Streamwise turbulence scalar flux budget determined from PIV (scatters) and DDES (lines). Profiles are extracted at $x/D = 2.1$ and $z/D = 0$ (a) for $VR = 0.4$ and (b) for $VR = 1.2$. The grey band represents the region of CGT. Terms are non-dimensionalized by ${U_0}$ and D.

To investigate the physical mechanism and coherent structure of the flow in the non-local CGT region, the co-spectra ${{{\mathcal F}}_{{u^{\prime}_i}c^{\prime}}}$ of the turbulence transport term are presented in figure 11(c,d). The results clearly show that both components of the scalar flux exhibit CGT behaviour, and the energetic structure responsible for the non-local CGT behaviour occurs at a much lower frequency compared with the local one. Since low-frequency components dominate the scalar flux, it can be inferred that these components, linked to the largest eddies in the flow, are the primary mechanism for the negative production of scalar fluctuation intensity seen in figure 8(c, f). The non-local nature can be further elucidated by the scale separation between mixing length theory and the large-scale structure. The large-scale vortical structures with a frequency of $St = 0.2$ and their high frequency harmonics dominate the unsteady behaviour near the wall. This frequency corresponds to a length scale of 3. According to Bodart et al. (Reference Bodart, Coletti and Eaton2013), these turbulent scales are much larger than the length scales over which ${\partial _i}C$ varies (less than 0.1 according to the mixing-length theory), meaning they can induce turbulent fluctuations that cannot be explained by local information alone. This indicated that the non-local nature of the fluxes – depending on the gradient not just at the point in question but over a broader region – should be considered.

Figure 11. Co-spectra of velocity and concentration at different locations determined using synchronized PIV/PLIF. (a,b) Local CGT region; (c,d) non-local CGT region. The grey band represents the CGT region. (a) $VR = 0.4$, locations P1 (2.8, 0.43) and P2 (2.8, 0.48); (b) $VR = 1.2$, locations P1 (0, 0.23) and P2 (0, 0.28); (c) $VR = 0.4$, locations P1 (2.1, 0.13) and P2 (2.1, 0.18); (d) $VR = 1.2$, locations P1 (2.1, 0.13) and P2 (2.1, 0.18). The notation $\textrm{sgn}({\cdot} )$ represents the sign function.

Another interesting phenomenon is that, unlike the dominant frequencies determined by the shear layer properties related to jet velocity, the dominant frequency in the non-local region remains constant for two different velocity ratios. This Strouhal number suggests wake vortex shedding, as reported by Moussa, Trischka & Eskinazi (Reference Moussa, Trischka and Eskinazi1977), who found that vortex shedding from JICF has a Strouhal number of 0.2 that is generally independent of VR. To support these observations, consider figure 12, which shows the non-local CGT regions identified using DDES and multiplane synchronized PIV/PLIF. Additionally, a snapshot of the Q criterion, coloured by ${u^{\prime}_\textrm{2}}c^{\prime}$, is depicted from a downward view. The distribution of the CGT region in the lateral direction resembles a classical wake shear layer, similar to those in the von Kármán street. The vortical structures further corroborate this observation. Thus, based on spectral characteristics, there is substantial evidence suggesting that the non-local CGT behaviour in the JICF may be associated with vortex shedding across the jet column.

Figure 12. Identification of non-local counter-gradient transport regions and visualization of the Q criterion. (a,b) Non-local counter-gradient transport regions identified using DDES and synchronized PIV/PLIF. (c,d) Snapshot of Q criterion (normalized $Q = 0.62$) coloured by ${u^{\prime}_2}c^{\prime}$ from a downward view, as determined by DDES. (a,c) $VR = 0.4$; (b,d) $VR = 1.2$.

3.3. Modelling implications

In the previous subsections, we identified locations exhibiting CGT behaviour and explored its origins. However, while insightful, the implication for turbulence modelling remains unclear. To delve into counter-gradient transport in turbulent scalar flux, we examine the angle ($\theta $) between the vectors $\overline {{u^{\prime}_i}c^{\prime}} $ and ${\partial _i}C$, expressed as

(3.6)\begin{equation}\theta = \textrm{arccos}\left( {\frac{{ - \overline {{u^{\prime}_i}c^{\prime}} {\partial_i}C}}{{|{\overline {{u^{\prime}_i}c^{\prime}} } |\cdot |{{\partial_i}C} |}}} \right).\end{equation}

In regions where the simple GDH holds precisely, $\theta $ should be $0^\circ $, indicating normal diffusion, whereas in areas exhibiting counter-gradient transport, $\theta $ exceeds $90^\circ $. Figure 13 shows contours of the angle $\theta $ at two different $VR\textrm{s}$, revealing key differences between previously identified local and non-local regions. Generally, $\theta $ is highest near the wall, where its presence induces strong anisotropy into the turbulent mixing. In the windward shear layer, where vertical turbulent motion drives local transport, $\theta $ typically ranges between $40^\circ $ and $70^\circ $, as the vertical component governed by the GDH is dominant. Conversely, regions of non-local transport, shown in figure 13(a,b), exhibit significant misalignment between $- \overline {{u^{\prime}_i}c^{\prime}} $ and ${\partial _i}C$ due to both components of the scalar flux depicting counter-gradient transport behaviour.

Figure 13. Counter-gradient transport at the mid-plane determined by synchronized PIV/PLIF. (a) $VR = 0.4$; (b) $VR = 1.2$.

From a modelling standpoint, regions where $\theta $ is significantly above zero cannot be well captured by the simple GDH described in (3.1), regardless of the diffusivity chosen. Introducing a spatially distributed $S{c_t}$ can capture some counter-gradient transport of turbulent scalar flux and would therefore be expected to model turbulent mixing more effectively. However, singular values of $S{c_t}$ could occur when $\theta > 90^\circ $ and the scalar transport equation may suffer from extreme numerical instability. Local and non-local behaviours require models of different complexity. Local behaviours can be modelled using local formulations but need a diffusivity matrix to capture cross-gradient effects (Milani et al. Reference Milani, Ling and Eaton2020). Non-local behaviours cannot be adequately modelled by any diffusivity-based model and may require solving separate transport equations for the scalar flux components. However, from a spectrum perspective, both the local and non-local behaviours correspond to large scales of the flow. The scalar flux can thus be decomposed into a gradient-type term representing the gradient diffusion and a term accounting for the stirring by large eddies. This approach extends the double-structure concept originally proposed by Townsend (Reference Townsend1976):

(3.7)\begin{equation}- \overline {{u^{\prime}_i}c^{\prime}} = \frac{{{\upsilon _t}}}{{S{c_t}}}{\partial _i}C + {( - \overline {{u^{\prime}_i}c^{\prime}} )^\ast },\end{equation}

where ${( - \overline {u^{\prime}_{i}c^{\prime}})^\ast }$ denotes the turbulent scalar flux due to large eddies.

4.1. Model deficiency representation

The classic RANS model attempts to estimate the unknown nonlinear source term in the mean equations to acquire the mean flow quantities. For incompressible mean flow, we first define a state vector $\boldsymbol{q}$ as

(4.1)\begin{equation}\boldsymbol{q} = \left\{ {\begin{array}{@{}c@{}} p\\ {{U_i}}\\ C \end{array}} \right\},\end{equation}

where p is the static pressure of the fluid; ${U_i}$ is the mean velocity vector, $i = 1,\;2,\;3$ in a three-dimensional space; and C is the mean scalar. We may then denote ${{\mathcal R}}$ as the conserved form of the incompressible RANS models, that is,

(4.2) \begin{equation}{{\mathcal R}}(\boldsymbol{q}) = \left\{ {\begin{array}{@{}c@{}} {{\partial_j}{U_j}}\\ {\partial_j}({U_i}{U_j}) + \dfrac{1}{\rho }{\partial_i}p - \upsilon {\partial_j}{\partial_j}{U_i} + {\partial_j}\overline {u_{i}^{\prime}{u_{j}^{\prime}} }\\ {{\partial_j}({U_j}C) - {\partial_j}\left( {\dfrac{\upsilon }{{Sc}}{\partial_j}C} \right) + {\partial_j}\overline {u_{j}^{\prime}c^{\prime}} } \end{array} } \right\} = \boldsymbol{0},\end{equation}

where $Sc$ is the Schmidt number and $- \overline {{u^{\prime}_i}{u^{\prime}_j}} $ is the Reynolds stress tensor. The commonly adopted closure approximation for the second-order correlation of fluctuation states follows a form based on the classical analogies of Reynolds and GDH (Pope Reference Pope2000), assumed as

(4.3)\begin{equation}\left. {\begin{array}{c@{}} { - \overline {{u^{\prime}_i}{u^{\prime}_j}} = {\upsilon_t}{\partial_j}{U_i} - \dfrac{2}{3}k{\delta_{ij}},}\\ { - \overline {{u^{\prime}_j}c^{\prime}} = \dfrac{{{\upsilon_t}}}{{S{c_t}}}{\partial_j}C,} \end{array}} \right\}\end{equation}

where k is the turbulent kinetic energy and ${\delta _{ij}}$ is the Kronecker symbol. A constant value for $S{c_t}$ is often adopted ($S{c_t} = 0.85$ in all of the following cases). As demonstrated in §§ 2 and 3, the LEVM and GDH pose great challenges to time-averaged predictions. The solution space impedes the capability of the DA method to reproduce the second-order correlation of fluctuations. To address modelling deficiencies, we introduce a spatially distributed body force correction ${F_i}$ and ${F_C}$ into the Reynolds stress and turbulent scalar flux, as follows:

(4.4)\begin{equation}\left. {\begin{array}{c@{}} { - {\partial_j}\overline {{u^{\prime}_i}{u^{\prime}_j}} = {\partial_j}\left( {{\upsilon_t}{\partial_j}{U_i} - \dfrac{2}{3}k{\delta_{ij}}} \right) + {F_i},}\\ { - {\partial_j}\overline {{u^{\prime}_j}c^{\prime}} = {\partial_j}\left( {\dfrac{{{\upsilon_t}}}{{S{c_t}}}{\partial_j}C} \right) + {F_C}.} \end{array}} \right\}\end{equation}

The inclusion of the LEVM and GDH in (4.4) serves to enhance the effective diffusion in the RANS equations, thereby ensuring numerical stability. Indeed, DA can also be performed to optimize ${\upsilon _t}$ and $S{c_t}$. However, as demonstrated by Ling et al. (Reference Ling, Ruiz, Lacaze and Oefelein2017), singular values would occur in strongly anisotropic flow in a JICF. For non-local behaviour, the optimized coefficient of the Laplace operator is negative, generating counter-gradient transport that concentrates low-concentration regions into high-concentration regions. This can cause local increases in concentration, leading to oscillations or even blow-ups in the numerical solution. The direct assimilation of ${F_C}$ ensures robustness and convergence of the iteration in the DA process. Overall, the DA scheme demonstrates stability and convergence similar to the baseline RANS model, with instabilities occurring only when the baseline RANS computation becomes unstable. To ensure strictly divergence-free estimated force correction, we introduce the Stokes–Helmholtz decomposition as follows:

(4.5)\begin{equation}{F_i} ={-} \partial \varphi /\partial {x_i} - {\varepsilon _{ijk}}\partial {\psi _k}/\partial {x_j},\end{equation}

with the scalar potential $\varphi $, the vector potential $\psi $ and the Levi–Civita symbol ${\varepsilon _{ijk}}$. For more details on the scalar and vector potential, please refer to our previous work (Li et al. Reference Li, He and Liu2022). By substituting (4.4)–(4.5) into (4.2) and integrating the scalar potential $\varphi $ into pressure p, we have

(4.6)\begin{equation}{{\mathcal R}}(\boldsymbol{q}) = \left\{ {\begin{array}{@{}c@{}} {{\partial_j}{U_j}}\\ {{\partial_j}({U_i}{U_j}) + \dfrac{1}{\rho }{\partial_i}\bar{p} - {\partial_j}[(\upsilon + {\upsilon_t}){\partial_j}{U_i}] - {f_i}}\\ {{\partial_j}({U_j}C) - {\partial_j}\left[ {\left( {\dfrac{\upsilon }{{Sc}} + \dfrac{{{\upsilon_t}}}{{S{c_t}}}} \right){\partial_j}C} \right] - {F_C}} \end{array}} \right\} = \boldsymbol{0}.\end{equation}

Here, the modified pressure $\bar{p}$ and solenoidal part of ${F_i}$ are expressed as follows:

(4.7)\begin{equation}\left. {\begin{array}{c@{}} {\bar{p} = p + {\textstyle{2 \over 3}}k + \rho \varphi }\\ {{f_i} ={-} {\varepsilon_{ijk}}\partial {\psi_k}/\partial {x_j}} \end{array}} \right\}.\end{equation}

4.2. Continuous adjoint formulation

The DA procedure is to find the optimal forcing ${f_i}$ and ${F_C}$ by minimizing the cost function ${{\mathcal J}}$, i.e. the discrepancy between the state variables obtained by the corrected model (4.6) and the limited observations, subject to the governing equations. This is expressed as follows:

(4.8)\begin{align}\left. {\begin{array}{*{20}{c}} {{{{\mathcal J}}_{{U_i}}} = \xi \int_{{\mathcal M}} {{{{\mathcal F}}_{{x_i}}}{{\left( {\dfrac{{{U_i} - {U_{Obs,i}}}}{{{U_0}}}} \right)}^2}\textrm{d}{{\mathcal M}}} \; + \dfrac{\zeta }{2}\int_{{\mathcal M}} {{{(\,{f_i})}^2}\,\textrm{d}{{\mathcal M}}\; } }\\ {{{{\mathcal J}}_C} = {\xi_C}\int_{{\mathcal M}} {{{{\mathcal F}}_{{x_i}}}{{\left( {\dfrac{{C - {C_{Obs}}}}{{{C_0}}}} \right)}^2}\textrm{d}{{\mathcal M}}} \; + \dfrac{\zeta }{2}\int_{{\mathcal M}} {{{({F_C})}^2}\,\textrm{d}{{\mathcal M}}\; } } \end{array}}\!\!\! \right\}\quad \!\!\!\textrm{s}\textrm{.t}\textrm{.}\;{{\mathcal R}}({{{\mathcal R}}_{NS}},{{{\mathcal R}}_p},{{{\mathcal R}}_C}) = 0.\end{align}

Here, ${{\mathcal M}}$ is the computational domain; ξ and ${\xi _C}$ are dimension converter with dimensions $[{L^2} \cdot {t^{ - 3}}]$ and $[{M^2} \cdot {L^{ - 6}} \cdot {t^{ - 1}}]$, respectively, and possess a value of unity. This ensures dimensional consistency. Term ${{{\mathcal F}}_{{x_i}}}$ is a masking function that is defined to specify the region where the measurements are made. If observation data are available, ${{{\mathcal F}}_{{x_i}}} = 1$ holds; otherwise, ${{{\mathcal F}}_{{x_i}}} = 0$. Here, we incorporate knowledge on the correction force by noting that the Reynolds stress and turbulent scalar flux usually exhibit large-scale structures that vary across the length scales of the mean flow (Franceschini, Sipp & Marquet Reference Franceschini, Sipp and Marquet2020). This term penalizes undesired oscillations resulting from sparse observations. The hyperparameter $\zeta $ is introduced to balance the measurement discrepancy and smoothness of the correction. Here, ${{\mathcal R}} = {{\mathcal R}}({{{\mathcal R}}_{NS}},{{{\mathcal R}}_p},{{{\mathcal R}}_C})$ are the momentum equations, the continuity equation, and the scalar transport equation presented in (4.6), respectively.

In a variational approach, the linear equations for the small perturbations of the mean flow are used to derive the adjoint equations. We first define the adjoint state vector ${\boldsymbol{q}^\ast }$, as follows:

(4.9)\begin{equation}{\boldsymbol{q}^\ast } = \left\{ {\begin{array}{@{}c@{}} {{p^\ast }}\\ {U_i^\ast }\\ {{C^\ast }} \end{array}} \right\}.\end{equation}

To minimize the energy of the cost function, the following Lagrangian function is used:

(4.10)\begin{equation}{{\mathcal L}}(\boldsymbol{q},{\boldsymbol{q}^{\boldsymbol{\ast }}},{f_i},{F_C}) = {{\mathcal J}} + \int_{{\mathcal M}} {(U_i^\ast {{{\mathcal R}}_{NS}} + {p^\ast }{{{\mathcal R}}_p} + {C^\ast }{{{\mathcal R}}_C})\,\textrm{d}{{\mathcal M}}} .\end{equation}

The Lagrangian multipliers $U_i^\ast $, ${p^\ast }$ and ${C^\ast }$ are the adjoint variables corresponding to (4.6), and are referred to as the adjoint velocity, adjoint pressure and adjoint scalar, respectively. The optimal distribution of ${f_i}$ and ${F_C}$ can be determined by obtaining the sensitivities of the Lagrange operator ${{\mathcal L}}$ with respect to ${f_i}$ and ${F_C}$, while setting the total variations of ${{\mathcal L}}$ with respect to the other state variables to vanish, according to

(4.11)\begin{equation}\delta {{\mathcal L}} \equiv \underbrace{{{\delta _{\boldsymbol{q}}}{{\mathcal L}}}}_{0} + {\delta _{{\boldsymbol{q}^\ast }}}{{\mathcal L}} + {\delta _{{f_i}}}{{\mathcal L}} + {\delta _{{F_T}}}{{\mathcal L}}.\end{equation}

The adjoint equations can be obtained by enforcing the zero constraint in (4.11), resulting in

(4.12)\begin{align}{{{\mathcal R}}^\ast }({\boldsymbol{q}^{\boldsymbol{\ast }}}) = \left\{ {\begin{array}{@{}c@{}} {{\partial_j}U_j^\ast }\\ {U_j^\ast {\partial_i}{U_j} - {U_j}{\partial_j}U_j^\ast+ \dfrac{1}{\rho }{\partial_i}{p^\ast } - {\partial_j}[(\upsilon + {\upsilon_t}){\partial_j}{U_i}] + 2\xi {{{\mathcal F}}_{{x_i}}}\dfrac{{{U_i} - {U_{Obs,i}}}}{{{U_0}}}}\\ { - {U_j}{\partial_j}{C^\ast } - {\partial_j}\left[ {\left( {\dfrac{\upsilon }{{Sc}} + \dfrac{{{\upsilon_t}}}{{S{c_t}}}} \right){\partial_j}{C^\ast }} \right] + 2{\xi_C}{{{\mathcal F}}_{{x_i}}}\dfrac{{C - {C_{Obs}}}}{{{C_0}}}} \end{array}} \right\} = \boldsymbol{0}.\end{align}

The boundary conditions are derived as follows. For the inflow and the wall, where the primary state variable ${U_i}$, C is specified, the boundary conditions are as follows:

(4.13)\begin{equation}\left. {\begin{array}{c@{}} {U_i^\ast= \boldsymbol{0},}\\ {{n_i}{\partial_i}{p^\ast } = 0,}\\ {{C^\ast } = 0,} \end{array}} \right\}\end{equation}

where ${n_i}$ is unit surface normal vector. For the outflow boundary, where the Neumann condition is used for the primary variables, the adjoint boundary conditions are as follows:

(4.14)\begin{equation}\left. {\begin{array}{c@{}} {{U_n} \cdot U_\tau^\ast{+} (\upsilon + {\upsilon_t})({n_i}{\partial_i}{U_\tau }) = \boldsymbol{0},}\\ {{U_n} \cdot U_n^\ast{+} (\upsilon + {\upsilon_t})({n_i}{\partial_i}{U_n}) = {p^\ast },}\\ {{U_n} \cdot {C^\ast } + \left( {\dfrac{\upsilon }{{Sc}} + \dfrac{{{\upsilon_t}}}{{S{c_t}}}} \right)({n_i}{\partial_i}{C^\ast }) = 0.} \end{array}} \right\}.\end{equation}

The subscripts n and $\tau $ denote the normal and tangential components of the variables, respectively. When $\boldsymbol{q}$ and ${\boldsymbol{q}^{\boldsymbol{\ast }}}$ satisfy (4.6) and (4.12)–(4.14), the sensitivity of the Lagrange function ${{\mathcal L}}$ becomes independent of the change in the adjoint variables, such that the gradient of the cost function can be computed, as follows:

(4.15)\begin{equation}\left. {\begin{array}{c@{}} {\dfrac{{\textrm{d}{{{\mathcal J}}_{{U_i}}}}}{{\textrm{d}{f_i}}} = \dfrac{{\partial {{\mathcal L}}}}{{\partial f}} ={-} U_i^\ast{+} \zeta {F_i},}\\ {\dfrac{{\textrm{d}{{{\mathcal J}}_C}}}{{\textrm{d}{F_C}}} = \dfrac{{\partial {{\mathcal L}}}}{{\partial {F_C}}} ={-} {C^\ast } + \zeta {F_C}.} \end{array}} \right\}\end{equation}

The gradient is determined for the entire computational domain and then used to update the force correction in an optimal sense, as

(4.16)\begin{equation}\left. {\begin{array}{c@{}} {f_i^{n + 1} = f_i^n - \beta \dfrac{{\partial {{\mathcal L}}}}{{\partial {f_i}}},}\\ {F_C^{n + 1} = F_C^n - \beta \dfrac{{\partial {{\mathcal L}}}}{{\partial {F_C}}},} \end{array}} \right\}\end{equation}

where $\beta $ is the step length. The iterative framework for flow reconstruction from mean field is as follows. It begins with an initial guess for ${f_i}$ and ${F_C}$, followed by an initialization procedure where (4.6) is solved without any corrections to achieve a converged base flow for the optimization process. The force correction is then used in (4.6) to obtain the primary variable distributions. Subsequently, the adjoint equations are solved using the boundary condition (4.13)–(4.14) computed by the cost function. The adjoint solution is then used to determine the sensitivity of the Lagrangian function (4.15), which is then employed to update the correction term (4.16). Finally, the iterative loop is repeated until a convergence criterion is met.

4.3. Data assimilation set-up

The benchmark flow configuration used for DA mirrors the experiment (figure 2a), with the mainstream inlet velocity profile derived from PIV results. Nevertheless, as demonstrated in our prior work (Li et al. Reference Li, He and Liu2022), the inlet velocity profile can also be accurately computed through a precursor RANS simulation. At the outflow boundary, a Dirichlet condition is specified for pressure and a Neumann condition for velocity. The concentration field is obtained via a scalar transport model, with rhodamine's mass fraction set to $5 \times {10^{ - 8}}$ at the jet inlet and 0 at the mainstream inlet, mirroring experimental conditions. Adjoint boundary conditions align with those in (4.13)–(4.14). The convection terms in the primary-adjoint system are discretized using a linear-upwind scheme to ensure second-order accuracy. To reduce the computational mesh size whilst maintaining reasonable flow blockage, slip-wall conditions are applied to neglect boundary layer effects on the side and upper walls (Li et al. Reference Li, He and Liu2023). A structured mesh with 1.05 million grids is used for the JICF configuration after a grid independence test. Observations for assimilation are extracted from the PIV/PLIF data at the $Z/D = 0$ plane (table 1), with exclusion of the near-wall region ($d/D < 0.2$) to mimic the large uncertainty. The regularization parameter $\zeta $ is fixed at 0.001 to mitigate any spurious fluctuations in correction determination. Further insights into the impact of regularization on convergence studies can be found in our earlier work (Li et al. Reference Li, Zhang, Zhou, He and Liu2024). Validation of velocity reconstruction is presented in § 5.1, followed by the application of DA for scalar field correction in § 5.2.

Table 1. Observation set-up and regularization for various cases