2.1. Experimental set-up

The ISWs were generated by a jet-array wave maker (JAW) (Chu, Lynett & Luhar Reference Chu, Lynett and Luhar2025) in a $2.2$ m long, $0.2$ m wide, and $0.5$ m deep acrylic flume. Within the flume, a miscible two-layer system was set up as background stratification, separated by a thin pycnocline. The lower layer was filled with a salt solution of density $\rho _2 = 1020 \pm 2$ kg m⁻³. This was confirmed via multiple measurements using a salinity probe. Fresh water with density $\rho _1 = 998$ kg m⁻³ was then introduced into the flume until the total depth $H = h_1 + h_2$ reached approximately $11$ cm. The experiments were carried out with fixed nominal depth ratio $h_2/h_1 = 9/2$ , with nominal wave amplitudes $a = -1, -1.5, -2$ cm.

The JAW system is composed of two independent, piston-driven cylinders, as shown in figure 1(b). The cylinders containing fresh water and salt solution are connected to the upper and lower layers of the flume, respectively. The pistons were driven by servo motors (Kollmorgen AKM44J), which were controlled using MATLAB scripts interfaced with a National Instruments data acquisition system (NI 6009) to generate prescribed volumetric fluxes (inflows and outflows) within each layer. To ensure near-uniform flow within each layer, fluid was driven from each cylinder through four uniformly spaced pipes, which subsequently fed through honeycomb flow straighteners. For each trial, the signal to drive the JAW system was prescribed based on the closed-form extended Korteweg–de Vries (eKdV) solution using the nominal layer depths and amplitude to determine the required mass flux for the inlet boundary conditions. While the fully nonlinear theories, such as the Dubreil–Jacotin–Long (DJL) equation, may provide a more mathematically rigorous description of ISWs (Lamb Reference Lamb2023), there are minimal differences in the computed inlet time series and initial wave profiles predicted using the DJL solver (Dunphy, Subich & Stastna Reference Dunphy, Subich and Stastna2011) and the eKdV solution for the parameter ranges considered here. We recognise that the wave evolves and transforms during the wave–structure interaction, such that the term ‘internal solitary wave’ may not apply. However, for simplicity, we retain the term ISW to describe the wave throughout the paper. The influence of a slight variation between nominal and actual layer depths on wave generation using the JAW system is addressed in Chu et al. (Reference Chu, Lynett and Luhar2025). For the conditions considered in this paper, the actual layer depths differed from the nominal values by ${\lt } 0.2$ cm. This is expected to translate into wave amplitude differences of 10 % or less relative to the nominal prescribed values.

Surface canopy structures were installed approximately $100$ cm downstream of the inlet, positioned just beneath the free surface at vertical location $z_c = 10$ cm as shown in figure 1(a). The canopies were rigidly constrained and not allowed to drift with the waves. The canopies were three-dimensionally printed to span the full channel width, with constant height $h_c = 1$ cm. Each canopy consisted of uniform element width $d$ , equal streamwise and lateral spacing $s_x = s_y$ , and vertical spacing $s_z$ , as illustrated in figure 2(a). Two canopy lengths ( $l_c = 5h_1 = 10$ cm and $l_c = 10h_1 = 20$ cm) and two porosity levels ( $n = 0.648, 0.964$ ) were tested. For the range of ISWs generated in our experiments, the characteristic wavelengths remained at $\lambda \approx 55$ cm. Consequently, both canopy lengths were significantly shorter than the wavelength ( $l_c \lt 0.5 \lambda$ ). Thus transient edge effects from the surface canopy are expected during the interaction. To isolate these transient effects, a numerical study is presented in Appendix A using an extended canopy with canopy length $l_c \gt 2 \lambda$ . The printed dimensions were $(d, s_x=s_y, s_z) = (0.4, 0.6, 0.6)$ cm for the lower-porosity canopies $(n=0.648)$ , and $(d, s_x=s_y, s_z) = (0.2, 2.3, 0.8)$ cm for the higher-porosity canopies $(n=0.964)$ .

Figure 2.

(a) Cross-sectional view and nomenclature of the canopy structure. (b) Isometric views of the transitional canopy with porosity $n=0.964$ . (c) Isometric views of the dense canopy with porosity $n=0.648$ .

These two porosities were chosen to span the so-called dense and transitional regimes for canopy flows (Nepf Reference Nepf2012). Previous work shows that flow penetration into canopies from adjacent non-vegetated regions depends on the dimensionless frontal area ratio $\lambda_{\kern-1pt f}$ , which represents the frontal area of the canopy per unit of planar area in the streamwise direction. For $\lambda_{\kern-1pt f} \ll 0.1$ , the canopy can be considered sparse, and there is significant flow penetration into the canopy. For $\lambda_{\kern-1pt f} \gg 0.1$ , the canopy can be considered dense. In this ‘dense’ case, a region of strong shear develops at the interface between the canopy and the adjacent unobstructed flow. This canopy shear layer can generate large-scale turbulent rollers via mechanisms similar to the Kelvin–Helmholtz instability. Canopies with $\lambda_{\kern-1pt f} \approx 0.1$ exhibit transitional behaviour between the sparse and dense regimes. For the canopy geometries considered in this paper, the less porous canopy has frontal area ratio $\lambda_{\kern-1pt f} \approx 0.64$ , and the more porous canopy has frontal area ratio $\lambda_{\kern-1pt f} \approx 0.16$ . These canopies may therefore be classified as dense and transitional, respectively.

Table 1 summarises the parameters for the $12$ experiments conducted. Each case is named after the prescribed wave amplitude (A), canopy length (L) and canopy type – dense (D) or transitional (T) – corresponding to its porosity. In this table, $h_{1,m}$ and $h_{2,m}$ are the measured layer depths obtained using planar laser-induced fluorescence (PLIF), as described below. For the following figures, experimental data are represented by hollow markers as listed in table 1, while the simulation results are represented by the corresponding solid markers.

Table 1.

Theoretical and measured values used for the experiments.

2.2. Measurement techniques

Interactions between ISWs and canopies were captured by synchronised PLIF and two-dimensional two-component particle imaging velocimetry (PIV). Polyamide particles with mean diameter $55\,{\unicode{x03BC}}$ m were introduced via floating sponges during the establishment of the upper layer for particle tracking. While the particles are slightly denser than both the fresh water ( $\rho _1 = 998$ kg m⁻³) and the salt solution ( $\rho _2 \approx 1020$ kg m⁻³), with mean density $1030$ kg m⁻³, the settling rate is negligible over the time scale of the ISW passage. Additional details are provided in Chu et al. (Reference Chu, Lynett and Luhar2025). However, PIV measurements in prior work were subject to insufficient particle density in the uppermost region of the freshwater layer, which is the location of our surface canopy structure, due to particle settling effects (Chu et al. Reference Chu, Lynett and Luhar2025). To ensure sufficient seeding density around the canopy, a revised two-step filling method was employed here. First, the freshwater layer was built with sponges floating on the free surface until the water level reached the canopy height. Second, the sponges were repositioned on top of the canopy as illustrated in figure 3, with the seeded fresh water slowly dripping into the upper layer. Due to the sponge repositioning, the pycnocline thickness estimated from PLIF increased to approximately $1.5$ cm in the experiments considered here (as opposed to remaining below $1$ cm in the previous work). The methodology for characterising the pycnocline position ( $\zeta$ ) and thickness ( $\Delta \zeta$ ) from the PLIF data follows the framework established by Chu et al. (Reference Chu, Lynett and Luhar2025), in which a hyperbolic tangent function is fitted to the PLIF intensity profile. The inflection point of the fitted profile defines the vertical position of the pycnocline, and the pycnocline thickness is determined by the coefficient parametrising the steepness of the hyperbolic tangent function; half the pycnocline thickness is absorbed into each layer depth ( $h_{1,m}, h_{2,m}$ ).

Figure 3.

An illustration of the two-step filling method to ensure sufficient particle density in the upper layer.

Figure 4.

(a) Schematic of the synchronised PLIF and PIV measurement system. Observation window coverage for (bi) short canopy, (bii) long canopy conditions.

The synchronised PLIF and PIV system consisted of two 8-bit CMOS monochrome cameras (Mako-U-130B) of $1024 \times 1280$ pixel resolution operating at $50$ Hz for $40$ s. The measurement duration was chosen to ensure that the entirety of the ISW was captured. A 5 W continuous laser with emission wavelength 532 nm and built-in sheet generation optics was used to illuminate the flow field from below (figure 4 a). For the PLIF measurements, Rhodamine 6G dye with peak excitation wavelength 525 nm and peak emission wavelength 548 nm was pre-mixed with the fresh water. An optical filter with a passband between 536 nm to 564 nm was therefore used to isolate the fluorescent freshwater layer. Assuming single-pixel resolution for the PLIF measurements, the dye field is resolved to better than 0.2 mm.

The cameras shared an overlapping field of view measuring $12\times15\ \text{cm}^2$ . For the short canopy case $(l_c/h_1 = 5)$ , the field of view extended from upstream to downstream of the canopy, enabling full observation of the interaction process, as illustrated in figure 4(bi). However, for the long canopy case $(l_c/h_1 = 10)$ , the field of view could not cover the entire canopy length. As such, the field of view was positioned to focus on the downstream half of the canopy. Specifically, the field of view encompassed slightly more than the downstream half of the canopy with an additional $2$ cm beyond, as illustrated in figure 4(bii). The MATLAB PIVLab package was used to post-process the images using a multi-pass algorithm. The initial pass involved a $64\times64\ \text{pixel}^2$ interrogation window. Two subsequent passes with $32\times32\ \text{pixel}^2$ interrogation windows with $50$ % overlap yielded $52\times32$ velocity vectors per snapshot.

3.1. Model set-up

Numerical simulations of ISW interactions with surface canopies were conducted using ANSYS Fluent. A two-dimensional computational domain was designed to match the experimental set-up illustrated in the previous section, consisting of total depth $H = 0.11$ m and flume length $L=2.2$ m. The domain was discretised using a uniform hexahedral mesh with resolution $\Delta x = 2\times 10^{-3}$ m and $\Delta z = 1\times 10^{-3}$ m. Grid resolution in the vertical direction $(z)$ was selected to be no greater than 10 % of the smallest wave amplitude. Grid resolution in the streamwise $(x)$ direction was controlled by the grid aspect ratio ( $\Delta x / \Delta z$ ) to avoid amplified numerical dissipation. To replicate the experimental conditions, the initial stratification was imposed using a hyperbolic tangent density profile with fixed pycnocline thickness 1.5 cm. The upper freshwater layer was assigned density $\rho _1 = 998$ kg m⁻³ and dynamic viscosity $\mu _1 = 1.00 \times 10^{-3}$ Pa s, while the lower salt solution layer was defined by $\rho _2 = 1020$ kg m⁻³ and dynamic viscosity $\mu _2 = 1.09 \times 10^{-3}$ Pa s.

To reproduce the wave-generation mechanism of the JAW system, mass flow inlets were applied at the upstream boundary of both layers, as illustrated in figure 5(a). The upper channel allowed only fresh water ( $\dot {m_1}$ ), while the lower channel allowed only ocean water ( $\dot {m_2}$ ). The mass flow rates were calculated by integrating the horizontal velocity profiles from the eKdV solution across the respective layer thicknesses, multiplied by the corresponding fluid density. No-slip boundary conditions were imposed along the channel walls and flume bed. The top surface was modelled as a no-shear boundary to approximate the free surface.

A mixture model was employed to simulate the multiphase configuration of the two-layer stratified system. Turbulence was modelled using the realisable $k\text{-}\epsilon$ Reynolds-averaged Navier–Stokes model. While the $k\text{-}\epsilon$ model has known limitations in handling strong adverse pressure gradients, recent comparisons with field measurements (Trowbridge et al. Reference Trowbridge, Helfrich, Reeder, Medley, Chang, Jan, Ramp and Yang2025) demonstrate its predictive accuracy in the bottom boundary layer beneath ISWs during free-stream acceleration and early deceleration. However, the model is unable to reproduce observations during the late deceleration and within the wake. In the present two-dimensional study, the $k\text{-}\epsilon$ model is employed primarily to characterise the macroscopic ISW–canopy interaction. The realisable formulation was specifically selected for its improved performance in capturing flow separation and boundary layers with variable eddy viscosity. To ensure accurate reproduction of wave amplitude evolution, the model coefficients within the $k\text{-}\epsilon$ formulation were calibrated to avoid overprediction of numerical dissipation relative to experimental measurements. Specifically, the coefficients were set as $c_{2,\epsilon } = 1.6$ , $\sigma _k =0.8$ and $\sigma _\epsilon =0.9$ . Numerical solution procedures employed the coupled pressure–velocity scheme combined with the PRESTO pressure discretisation. Second-order implicit upwind schemes were applied to all transport equations, including momentum, turbulent kinetic energy, dissipation rate and transient formulation. A fixed time step $ \Delta t = 5 \times 10^{-3}$ s was employed throughout all simulations. This value was selected to satisfy the Courant–Friedrichs–Lewy (CFL) stability criteria. The CFL number was estimated based on the dominant streamwise propagation of an ISW of $a = 2$ cm by using the nominal phase speed $c$ as the upper limit for the horizontal particle speed. This estimation resulted in a maximum CFL number $c\, \Delta t/ \Delta x \approx 0.18$ .

Figure 5(b) presents a snapshot of density, horizontal velocity and vertical velocity contours of a simulated ISW corresponding to $(h_1, h_2, a) = (2,9,2)$ cm, evaluated $1$ m downstream from the inlet. For consistency in the subsequent analysis, the streamwise coordinate origin $(x=0)$ is defined at $1$ m downstream from the inlet, where the leading edge of the surface canopies would otherwise be. The wave profile of the generated ISW in our numerical model is closely aligned with the eKdV solution for the same layer depths and wave amplitude.

Figure 5.

(a) Schematic showing the numerical model. (b) Density contour for a simulated ISW. The black dashed curve represents the analytical eKdV profile used as the inlet boundary condition.

3.2. Canopy zone characterisation

The surface canopy structure was modelled as a homogeneous porous zone within the computational domain, as illustrated in figure 5. The porous region was defined by a specified porosity $n$ , and a viscous resistance parameter $1/K$ , representing the inverse of permeability ( $K$ ). While the porosity was directly calculated as listed in table 1, accurate determination of viscous resistance was more challenging in experiments due to the small measurable pressure gradient across our canopy dimensions.

Previous experimental studies by Vijay & Luhar (Reference Vijay and Luhar2024) characterised discrete viscous resistance values for both isotropic and anisotropic porous substrates with cubic lattice geometries similar to those considered in this paper. However, no general predictive relationship was established for estimating the viscous resistance of these porous structures. In the present study, the following Kozeny–Carman (KC) equation was adopted as a base model for estimating permeability (Carman Reference Carman1939):

(3.1) \begin{align} K_{\textit{KC}} = \frac {\beta }{\sigma ^2} \frac {n^3}{(1-n)^2}. \end{align}

In this model, $\sigma$ represents the specific surface area for the structure, and $\beta$ is the so-called Kozeny constant, which depends on the specific porous geometry. The KC model was initially derived under the assumption of non-interconnected, parallel tube bundles as a simplified representation of complicated porous formations. More recent findings by Schulz et al. (Reference Schulz, Ray, Zech, Rupp and Knabner2019) indicate that a value $\beta =0.2$ yields the best agreement with experimental data for irregular porous structures. The specific surface is defined as the surface area $S_c$ divided by the solid volume $V_c$ of a unit canopy cell. For the geometries considered in this paper, the unit canopy cell consists of a central cube with side length $d$ , and three square columns oriented along the principal axes as shown in figure 6. As such, the specific surface area can be estimated using the equation

(3.2) \begin{align} \sigma = \frac {S_c}{V_c} = \frac {6d^2 + 4 s_x d + 4 s_y d + 4 s_z d}{d^3 + d^2s_x + d^2 s_y + d^2 s_z}. \end{align}

Figure 6.

Unit cell for the calculation of a specific surface.

To evaluate the applicability of the KC model for the cubic lattice geometries considered in this paper, the predicted resistance values were compared against the measurements reported by Vijay & Luhar (Reference Vijay and Luhar2024), as shown in table 2. The experiments of Vijay & Luhar (Reference Vijay and Luhar2024) considered cubic cell geometries with fixed stem thickness $d= 4\times 10^{-3}$ m and varying spacings in the principal directions. Each specimen was labelled with a three-letter code representing the spacings along the principal directions. Here, we consider the isotropic cases HHH, MMM and LLL, denoting high, medium and low spacings (i.e. going from more porous to less porous). As shown in table 2, the KC model provides a reasonable estimate of the measured viscous resistances, with predictions agreeing with the measurements within $\pm 15\,\%$ for the HHH and MMM cases. For the densest LLL case considered by Vijay & Luhar (Reference Vijay and Luhar2024), the predicted viscous resistance is a factor of 2–3 times higher. Following this limited validation exercise, additional numerical sensitivity analyses were carried out to identify the most appropriate resistance values for the canopies considered in the experiments.

Table 2.

Viscous resistance (permeability) estimates obtained using the KC model, and comparison with measurements by Vijay & Luhar (Reference Vijay and Luhar2024) for similar geometries.

Table 3.

Sensitivity tests evaluating the impact of the chosen viscous resistance on the integrated velocity error $e$ from (3.3).

As a sensitivity test, numerical simulations were performed for the dense and transitional surface canopies subjected to ISWs using the exact viscous resistance values predicted by the KC model, and values multiplied by a factor of 5 lower and higher, i.e. corresponding to viscous resistances $5/K_{\textit{KC}}$ , $1/K_{\textit{KC}}$ and $1/5K_{\textit{KC}}$ . These tests compared the predicted and measured horizontal velocity profiles at the primary wave crest for ISWs of $a = 2$ cm under the long transitional and long dense canopies (A2L2T and A2L2D). Velocity profiles were compared at the streamwise location $x/h_1 = 10.5$ immediately downstream of the canopy, where experimental measurements were first available across the full water depth following the interaction between the ISW and the surface canopy. The difference was quantified using the $l_1$ norm, normalised by the nominal phase speed of the ISW, $c$ , defined as

(3.3) \begin{align} e = \frac {1}{N} \sum _{i=1}^N\left |\frac {u_{n, i}-u_{m, i}}{c}\right |, \end{align}

where $u_n$ and $u_m$ are the numerical and measured horizontal velocity components, and $N$ is the total number of data points across the entire water depth.

Figure 7.

Comparison of horizontal velocity profiles between experiments and simulations with varying viscous resistances for cases (a) A2L2T and (b) A2L2D. Red hollow markers represent experimental profiles extracted from PIV measurements. Solid curves indicate $1/K_{\textit{KC}}$ , dashed curves indicate $5/K_{\textit{KC}}$ , and dotted curves indicate $1/5K_{\textit{KC}}$ .

Viscous resistance values scaled by a factor of $5$ from the KC model yielded better agreement with experimental results (table 3). For both the transitional and dense canopies, this adjustment more accurately captured the acceleration within the upper-layer fluid and the deceleration downstream of the canopy structure, as demonstrated by the velocity profiles in figure 7. Two factors may contribute to the consistently higher viscous resistance from the KC model. First, the unit cell used to compute the specific surface area (figure 6) neglects the contribution from the bottom stem surfaces. This exclusion may have a negligible impact on multi-layer canopy structures such as those used in the experiments by Vijay & Luhar (Reference Vijay and Luhar2024). However, it becomes increasingly significant in single-layer structures such as the canopies employed in the present study. In particular, the fully exposed lower interface of the canopy to the surrounding fluid likely enhances the impact of the bottom layer on overall resistance. Second, the KC model strictly applies only for the Darcy (creeping) flow regime, characterised by a low pore-scale Reynolds number $Re_p \lt 1$ . In the present study, an estimated pore-scale Reynolds number $Re_p = u d/ \nu \gt O(10)$ – using the stem width $d$ , the horizontal velocity component within the upper layer fluid $u$ , and the kinematic viscosity of water $\nu =10^{-6}$ m $^2$ s⁻¹ – places the flow in the transitional regime. The elevated $Re_p$ implies enhanced inertial effects and stronger resistance compared to the conditions assumed in the original KC formulation. Despite these limitations, the numerical simulations capture experimental trends reasonably well using viscous resistance values of $5/K_{\textit{KC}}$ .