Surfzone eddies enhance the dispersion and transport of contaminants, bacteria, and larvae across the nearshore, altering coastal water quality and ecosystem health. During directionally spread wave conditions, vertical vortices (horizontal eddies) are injected near the ends of breaking crests. Energy associated with these eddies may be transferred to larger-scale, low-frequency rotational motions through an inverse energy cascade, consistent with two-dimensional turbulence. However, our understanding of the relationships between the wave conditions and the dynamics and energetics of low-frequency surfzone eddies are largely based on numerical modeling. Here, we test these relationships with remotely sensed and *in situ* observations from large-scale directional wave basin experiments with varying wave conditions over alongshore-uniform barred bathymetry. Surface velocities derived with particle image velocimetry were employed to assess the spatial scales of low-frequency surfzone eddies and compute structure functions with alongshore velocities. Second-order structure functions for directionally spread waves ( $ \sigma \theta \u2265 10 \xb0$) are consistent with energy flux to larger or smaller length scales, while normally incident, unidirectional waves do not display this behavior. Third-order structure functions suggest that the surfzone flows exhibit a bidirectional energy cascade—a direct cascade to smaller and inverse cascade to larger length scales—during large directional spreads waves ( $ \sigma \theta \u2265 18 \xb0$). However, there is not decisive evidence of an inverse energy cascade for moderate directional spreads ( $ \sigma \theta = 10 \xb0$). Energy flux varies by cross-shore location and increases with increasing directional spread and wave height. Eddy decorrelation length scales weakly depend on wave directional spread. These findings advance our understanding of the dynamics linking wave breaking to large-scale rotational motions that enhance mixing and lead to rip currents, important conduits for cross-shore material exchange.

## I. INTRODUCTION

The surf zone extends from the shoreline to the seaward limit of depth-induced wave breaking. As the boundary between land and open ocean, this region is widely relied upon for beach recreation and tourism and plays a crucial role in coastal ecosystem health. Water quality within the surf zone is often degraded due to the anthropogenic influences, such as pathogens from point or non-point sources and excess nutrients from terrestrial runoff (Halpern , 2008; Boehm and Soller, 2013), prompting beach closures and harming coastal ecosystems (Boehm , 2002; Moret , 2005). Transport and dispersion of materials (i.e., pollutants, nutrients, and larvae) within the surf zone and between the surf zone and deeper waters is primarily controlled by wave-forced surfzone currents and eddies (Hally–Rosendahl , 2014; 2015; and Hally–Rosendahl and Feddersen, 2016). Over relatively short timescales (tens of minutes), surfzone eddies—horizontal rotational motions with lower frequencies than surface gravity waves—enhance lateral dispersion and lead to offshore ejections of material in transient rip currents (Johnson and Pattiaratchi, 2004; Spydell , 2007; Brown , 2009; Clark , 2010; Feddersen , 2011; Hally–Rosendahl , 2014; 2015; Suanda and Feddersen, 2015; Hally–Rosendahl and Feddersen, 2016; Kumar and Feddersen, 2017b; 2017a; Grimes , 2020a; 2020b; and Choi and Roh, 2021). Despite the ubiquity and critical role of surfzone eddies in determining the fate of material within the nearshore, the forcing, dynamics, and energetics of these eddies are still poorly understood.

Horizontal surfzone eddies, or vertical vortices, can be generated by individual breaking crests (Peregrine, 1998). Eddies can also be generated by wave group forcing (Reniers , 2004; Long and Özkan Haller, 2009) or shear instabilities in mean alongshore currents (Bowen and Holman, 1989; Dodd , 1992; Oltman–Shay , 1989; Feddersen, 1998; and Özkan Haller and Kirby, 1999). Recent findings indicate that vertical vorticity injected by individual breaking waves may be the primary eddy-generation mechanism on alongshore-uniform beaches (Feddersen, 2014; Elgar and Raubenheimer, 2020), and of comparable importance to bathymetrically controlled eddies on alongshore-variable beaches (O'Dea , 2021; Baker , 2021). Under directionally spread wave conditions—when wave energy is distributed broadly around the mean angle—interference patterns result in alongshore varying wave amplitudes (Longuet-Higgins, 1956). In the surf zone, the highest amplitude sections of the wave break first, farther offshore, along finite regions, known as short-crested wave breaking. In a wave-averaged reference, these forcing terms are alongshore uniform, and other eddy-generation mechanisms are not predicted to occur for shore-normal, directionally spread waves (e.g., Long and Özkan Haller, 2005; 2009). On a wave-by-wave basis, the along-crest gradient in the breaking force, with magnitudes and length scales associated with breaking wave characteristics, leads to the injection of positive and negative vertical vorticity near crest ends [Fig. 1(a); Peregrine and Bokhove, 1998; Bühler and Jacobson, 2001; Johnson and Pattiaratchi, 2006; Sullivan 2007; Clark 2012; and Choi and Roh, 2021].

Energy associated with these breaking-wave-forced eddies may be nonlinearly transferred to larger scale rotational motions through an inverse energy cascade, consistent with a forced two-dimensional (2D) turbulence system [Fig. 1(b); Peregrine and Bokhove, 1998; Spydell and Feddersen, 2009; Feddersen, 2014; and Elgar and Raubenheimer, 2020]. Forced 2D turbulence dynamics have been observed in idealized laboratory experiments (Sommeria, 1986; Paret and Tabeling, 1997; Rutgers, 1998; Vorobieff , 1999; and Clercx and Van Heijst, 2009), in numerical simulations (Lilly, 1969; Siggia and Aref, 1981; Frisch and Sulem, 1984; and Herring and McWilliams, 1985), and in multiple natural systems that exhibit large lateral to vertical aspect ratios (Celani , 2010; Kelley and Ouellette, 2011; and Xia , 2011), including large-scale motions in the atmosphere (Leith, 1971; Gage and Nastrom, 1986; Lilly, 1989; and Lindborg, 1999) and ocean (Rhines, 1979; Stewart , 1996; Scott and Wang, 2005; and Balwada , 2022).

In the surf zone, an inverse energy cascade may provide an explanation for the enhanced very low-frequency eddy activity observed for directionally spread wave fields (e.g., at frequencies around an order of magnitude lower than the peak wave frequency; Spydell , 2007; Spydell , 2009; Spydell and Feddersen, 2009; Feddersen, 2014; Spydell, 2016; and Elgar and Raubenheimer, 2020) and for the broad distribution of low-frequency motions ubiquitous on alongshore-uniform beaches (Noyes , 2005; Spydell and Feddersen, 2009; and Hally–Rosendahl , 2015). Eddies may become quasi-2D in the outer edge of the surf zone, particularly for the highest frequency rotational motions (Lippmann and Bowen, 2016; Henderson , 2017; and Baker , 2021), where vortex tilting, vortex stretching, and vertical diffusion may lead to more rapid decay of eddies than expected for solely 2D eddies (Uchiyama , 2017; McWilliams , 2018).

Despite the relevance of an inverse energy cascade mechanism within the surf zone, there are a limited number of modeling (Spydell and Feddersen, 2009; Feddersen, 2014; and Marchesiello , 2021) and observational (Elgar and Raubenheimer, 2020; Elgar , 2023) studies assessing the dynamical relevance of a forced 2D turbulence system for surfzone processes. Field observations provide evidence that surfzone flows are consistent with 2D turbulence when waves are directionally spread (Elgar , 2023), but there is no observational evidence to address how the presence, energetics, and limiting length scales of an inverse energy cascade vary with wave forcing. In the field environment, it is difficult to sample the wide range of eddy spatial and temporal scales as well as isolate the wave-forced-eddy inverse energy cascade and its variation with wave conditions. Thus, we performed experiments in a laboratory wave basin with a directionally spread wave field on an alongshore-uniform barred beach (Baker , 2023) to address how nonlinear transfers of energy between breaking-wave-injected vorticity and larger scale rotational motions vary with the wave conditions. We employed particle image velocimetry (PIV) on optical imagery to characterize surfzone flow fields and assess their consistency with 2D turbulence.

In this paper, we provide a brief overview of 2D turbulence and its applications to surfzone flows in Sec. II. Laboratory setup and analyses are described in Sec. III, and results assessing the low-frequency eddy activity in the surf zone with *in situ* and remote sensing observations are in Sec. IV. Implications of our results and comparisons with previous literature are discussed in Sec. V, and a brief summary of findings is in Sec. VI.

## II. 2D TURBULENCE THEORETICAL OVERVIEW AND APPLICABILITY TO THE SURF ZONE

*g*is the gravitational acceleration,

*h*is the mean water depth, η is the sea level relative to a still water level, and $ h + \eta $ is the total water depth. The body forces are a source term representing the depth-averaged wave breaking external force ( $ F br$), and a sink term ( $ F s$) representing friction [e.g., $ \tau b / ( \rho h )$, the bottom stress normalized by

*h*and the water density, ρ] and viscous forces (e.g., $ \nu \u2207 2 u$, where ν is the kinematic viscosity). Other forcing terms such as wind stress are not considered here.

Consider these dynamics in a closed fluid system of constant depth in which vorticity is continually injected at an intermediate horizontal length scale (*l _{i}*) by external rotational forcing. In three-dimensional (3D) turbulence, vortex stretching and strain self-amplification result in a forward energy cascade from

*l*to smaller length scales (Pope, 2000), where energy is eventually dissipated by viscous forces. However, in a forced, steady-state, strictly 2D turbulence system, there is no axis perpendicular to the flow, and thus vortex stretching is absent (Kraichnan, 1967; Kellay and Goldburg, 2002; and Boffetta and Ecke, 2012). Instead, energy is transferred nonlinearly to larger scales through an inverse cascade, ultimately reaching a maximum length scale constrained by the horizontal spatial extent of the system. At the largest scales of the system, energy is dissipated (Batchelor, 1969; Paret and Tabeling, 1997; and Boffetta , 2000) or accumulates, leading to spectral condensation if limited dissipative forces are present (Hossain , 1983; Smith and Yakhot, 1993; Smithr and Yakhot, 1994; Chertkov , 2007; Xia , 2011; and Fang and Ouellette, 2021). Vorticity gradients, however, are not bounded in two-dimensions, enabling a direct cascade of enstrophy ( $ \Omega = 1 2 \u27e8 \omega 2 \u27e9$, spatial average of vorticity squared) from

_{i}*l*to smaller length scales

_{i}*l*, where the enstrophy is dissipated by viscous forces (Leith, 1968; Rutgers, 1998). In quasi-2D turbulence with a bidirectional energy cascade, the remaining energy that is not transferred to larger scales cascades to smaller scales through a 3D direct energy cascade (Smith , 1996; Celani , 2010; Shats , 2010; Byrne , 2011; Xia , 2011; Musacchio and Boffetta, 2017; Alexakis and Biferale, 2018; and Musacchio and Boffetta, 2019).

_{D}The surfzone geometry has some similarities and some substantial differences with the classical flat-bottomed 2D turbulence assumptions. The surf zone has spatially varying water depths, and thus, the assumptions leading from Eqs. (3) to (4) may not be satisfied. Expanding the second term on the left-hand side of Eq. (3), $ u \xb7 \u2207 [ \omega h + \eta ] = u \u2202 \u2202 x ( \omega h + \eta ) + v \u2202 \u2202 y ( \omega h + \eta )$. Considering alongshore-uniform bathymetry for simplicity, the term with an alongshore derivative is zero, and the first term can be expanded to $ u \u2202 \u2202 x ( \omega h + \eta ) = ( 1 h + \eta ) u \u2202 \omega \u2202 x \u2212 u \omega h 2 \u2202 h \u2202 x$. The second term $ u \omega h 2 \u2202 h \u2202 x$ is small when variations in the water depth ( $ \Delta h$) are small relative to the water depth (*h*). For typical field and laboratory surfzone bathymetries, this may be reasonably satisfied in the surf zone for shallow sloped beaches and in the relatively flat trough region of barred or terraced beaches. For example, for the laboratory barred bathymetry presented here, $ \Delta h / h$ is small for most of the surf zone region, and thus, it is reasonable to consider the similarities to eddy evolution in a flat-bottom system. Specifically, $ \Delta h / h \u2264$ 0.2 for ∼85% of distance between the bar crest and the shoreline. ( $ \Delta h / h$) is larger for the remaining ∼15% of the surf zone, including near the shoreline where water depths are small and swashzone processes dominate, and in the deepest part of the trough, where $ \Delta h$ relative to the bar crest approaches the same order as the average water depth.

Surfzone vorticity forcing may be reasonably analogous to a closed forced 2D turbulence system with some notable differences (Bühler and Jacobson, 2001; Terrile and Brocchini, 2007). Short-crested breaking waves inject vorticity with length scales *l _{i}* associated with the along-crest gradient in the breaking force [ $ \u2207 \xd7 F br$, Eq. (4), Fig. 1(a)], which may be maximum near crest ends at the transitions from a broken ( $ F br \u2260 0$) to unbroken water surface ( $ F br = 0$; Bühler and Jacobson, 2001; Johnson and Pattiaratchi, 2006; Clark , 2012; and Kirby and Derakhti, 2019). For short-crested wave breaking in deep water, a half vortex-ring between crest ends dissipates similar to 3D turbulence (Peregrine and Bokhove, 1999; Pizzo and Melville, 2013). In contrast, for short-crested wave breaking in shallow water, oppositely signed vertical vortices generated near crest ends have scales greater than the water depth [Fig. 1(a); Peregrine and Bokhove, 1998; Spydell , 2019; and Elgar and Raubenheimer, 2020]. In an irregular, directionally spread wave field, the strength and scales of injected vortices may depend on individual crest characteristics and may change as waves propagate shoreward over varying bathymetry and interact with currents. Further work is needed to understand to what extent the unsteady and localized aspects of the forcing lead to inhomogeneous eddy fields or other differences with a classic forced 2D turbulence system (Clark , 2012; Baker, 2023).

Surfzone dissipation and boundary conditions vary from classical 2D turbulence with some similarities. If surfzone flow is consistent with 2D turbulence, energy associated with injected surfzone vortices will be non-linearly transferred to larger scales (*l _{0}*) that may be constrained by the nature of dissipative processes such as bottom friction (Fang and Ouellette, 2021; Grianik , 2004) and the surfzone geometry [e.g., the surfzone width, $ W sz$, Fig. 1(b)]. In the surf zone, drag may be enhanced by wave orbital velocities and bottom roughness. Consequently, it is unclear if friction is small enough that a turbulent cascade can develop. In addition, the surf zone has an open boundary eddy-kinetic energy flux to the inner shelf, i.e., eddies can be ejected from the outer surf zone to the shelf in the form of a transient rip current [Kumar and Feddersen, 2017a; Grimes and Feddersen, 2021, Fig. 1(b)], where dynamics transition to a quasi-2D system (Lippmann and Bowen, 2016; Baker , 2021). The ratio of bottom friction to the wave forcing term can be estimated by approximating bottom stress with a quadratic drag $ \u221d C d u 0 u \xaf$, where

*u*

_{0}is the wave orbital velocity and $ u \xaf$ is the mean flow velocity (Feddersen , 2000), and approximating the wave forcing term for shore-normal waves, $ \u221d \u2202 \u2202 x 3 2 E$ where

*E*is the wave energy, assuming depth-limited breaking with threshold γ

_{b}(Longuet-Higgins and Stewart, 1964). While the spatial scales of short-crested breaking may be shorter than variations in the mean flow, we conservatively assume that these vary over similar lengthscales, such that the ratio of body forcing terms $ | F s | / | F br |$ is similar to the ratio of vorticity forcing terms $ | \u2207 \xd7 F s | / | \u2207 \xd7 F br |$. This ratio then scales as $ ( C d u 0 u \xaf ) ( 3 16 g \gamma b 2 \alpha h ) \u2212 1$, where α is the beach slope. For the laboratory conditions investigated here, we estimate this ratio to be of order $ 10 \u2212 2$, with a conservative range of $ 10 \u2212 3$ to $ 10 \u2212 1$, and thus, we expect that the forcing is dominant over friction and an inverse cascade could develop. For field conditions with shore-normal waves, the ratio is also typically small, however may be closer to order 1 for conditions with smaller slopes and water depths with larger friction and flows.

## III. LABORATORY OBSERVATIONS

### A. Experimental setup and observations

Laboratory experiments were performed in the Directional Wave Basin at the O.H. Hinsdale Wave Research Laboratory (HWRL), Oregon State University, USA (Fig. 2). The basin dimensions are 48.8 m in the cross-shore direction, 26.5 m in the alongshore direction, and up to 2 m deep. The origin of the basin coordinate system is centered along the width of the basin at the bottom of the wavemaker. The positive coordinate in the cross-shore (*x*) is to the west toward the beach, alongshore (*y*) is toward the south (right handed coordinate system), and elevation (*z*) is upward from the base of the basin. The water depth (*h*) of the bar crest was *h* = 0.19 m when the basin still water level was set to *h* = 1.07. The basin bathymetry was configured with a flat region with $ z = 0$ m from *x* = 0–22 m and a planar $ 1 : 10$ sloping beach from $ x = 22 \u2013 48.8$ m. An alongshore-uniform concrete bar with a triangular cross-section rested on the metal slope and extended from $ x = 25.4 \u2013 29.2$ m with a 0.30 m tall apex [Fig. 2(b)]. Although the 1:10 slope may be steeper than found on natural beaches, the surf similarity number (the ratio of wave steepness to beach slope; Battjes, 1975) resembles some field settings (e.g., *H _{s}* = 1 m,

*T*= 10 s, and a 2.5°–3° beach slope). See Baker (2023) for more discussion on this topic.

The piston-type wavemaker consisting of 29 segments can simulate multidirectional wave fields following the second-order multi-directional wavemaker theory (Schäffer and Steenberg, 2003). The wave field was generated using a JONSWAP spectrum given a significant wave height (*H _{s}*), peak wave period (

*T*), and a spectral width gamma parameter set to 3.3 (Hasselmann , 1973) and a $ \u2009 cos 2 s$ directional distribution from a mean wave angle (θ) and directional spread ( $ \sigma \theta $; Baker , 2023). Wave conditions offshore of the toe of the beach slope were measured with 15 surface-piercing wire resistance gauges (

_{p}*x*= 19.0, 20.4 m, Fig. 2, squares). Trials are reported in terms of the offshore wave conditions estimated from a sub-sampled array of the wire resistance gauges from minute 15–25 of each 45 min trial (Table I, grouped by wave heights less than or greater than 0.25 m). For the trials presented here, the offshore

*H*varied from 0.24 to 0.30 m, where the onset of breaking typically occurred just offshore of or near the bar crest ( $ x = 27.1 \u2212 27.6$ m), and

_{s}*T*was 2 s. Wave conditions were near shore-normal ( $ \theta \u2248 0 \xb0$) with directional spreads ranging from narrow to large ( $ \sigma \theta = 2 \xb0 \u2212 26 \xb0$).

_{p}Trial . | H (m)
. _{s} | T (s)
. _{p} | $ \theta ( \xb0 )$ . | $ \sigma \theta \u2009 ( \xb0 )$ . | $ x sz$ (m) . | λ (m)
. _{c} |
---|---|---|---|---|---|---|

G1a | 0.28 | 2.1 | −0.1 | 2.4 | 27.1 | 16.5 |

G1b | 0.27 | 2.0 | 0.4 | 16.8 | 27.3 | 7.7 |

G1c | 0.26 | 2.0 | −5.9 | 25.2 | 27.3 | 5.0 |

G1d | 0.27 | 2.0 | −3.3 | 26.1 | 27.3 | 4.5 |

G2a | 0.25 | 2.0 | 0.2 | 2.4 | 27.2 | 4.2 |

G2b | 0.23 | 2.0 | −1.0 | 9.6 | 27.4 | 17.8 |

G2c | 0.22 | 2.0 | 0.8 | 18.3 | 27.6 | 5.7 |

G2d | 0.24 | 2.0 | 4.8 | 24.2 | 27.6 | 4.5 |

G2e | 0.23 | 2.0 | 2.1 | 25.9 | 27.5 | 4.7 |

Trial . | H (m)
. _{s} | T (s)
. _{p} | $ \theta ( \xb0 )$ . | $ \sigma \theta \u2009 ( \xb0 )$ . | $ x sz$ (m) . | λ (m)
. _{c} |
---|---|---|---|---|---|---|

G1a | 0.28 | 2.1 | −0.1 | 2.4 | 27.1 | 16.5 |

G1b | 0.27 | 2.0 | 0.4 | 16.8 | 27.3 | 7.7 |

G1c | 0.26 | 2.0 | −5.9 | 25.2 | 27.3 | 5.0 |

G1d | 0.27 | 2.0 | −3.3 | 26.1 | 27.3 | 4.5 |

G2a | 0.25 | 2.0 | 0.2 | 2.4 | 27.2 | 4.2 |

G2b | 0.23 | 2.0 | −1.0 | 9.6 | 27.4 | 17.8 |

G2c | 0.22 | 2.0 | 0.8 | 18.3 | 27.6 | 5.7 |

G2d | 0.24 | 2.0 | 4.8 | 24.2 | 27.6 | 4.5 |

G2e | 0.23 | 2.0 | 2.1 | 25.9 | 27.5 | 4.7 |

On the beach slope, near-bed pressure and velocities were measured with a 12-sensor array of colocated pressure gauges and acoustic Doppler velocimeters (ADVs, Vectrino Profilers) mounted about 0.05 m above the bed. Colocated sensors were first deployed in two alongshore transects in the surf zone (*x* = 28.4 and 30.7 m, $ y = \u2212 8.2$ to 3.2 m, Fig. 2, circles and triangles) and then re-deployed for repeated wave conditions in the inner shelf (Fig. 2, diamonds). The inner-shelf array consisted of an alongshore transect spanning the width of the tank ( $ y = \u2212 9.8$ to 9.7 m) just offshore of the surf zone (*x* = 26.6 m) and three vertically stacked sensors farther offshore ( $ x , y = 24.5 , \u2212 0.1$ m). The ADV at station p9 ( $ x , y = 30.7 , \u2212 6.6$ m) did not record data during several trials (e.g., G2*b,d*). All *in situ* sensors synchronously sampled at 100 Hz. Additional details about the experimental setup, sensor quality control, and computation of offshore wave conditions are described in Baker (2023).

Our analysis considers images from the alongshore-centered camera (c2, $ x , y , z = 47.0 , \u2009 0.0 , \u2009 10.1$ m; Baker , 2023) that was oriented in the offshore direction with an oblique view of the water surface (268.8° azimuth, $ 22.8 \xb0$ tilt, $ 0.1 \xb0$ roll). Images were collected at 8 Hz with $ \u2264 1$ cm resolution in the surf zone and up to 2 cm resolution near the bridge. Images were rectified into HWRL coordinates (Bruder and Brodie, 2020; Palmsten and Brodie, 2022) with the intrinsic and extrinsic parameters from stereo processing (Baker , 2023). The z-elevation for rectifications was obtained from the stereo derived elevation maps in the surf zone and set to the still water elevation offshore. Additionally, a multibeam lidar ( $ x , y , z = 22.0 , \u2212 3.0 , 7.3$ m) scanned the water surface at 10 Hz with a footprint spanning farther offshore than the stereo reconstructions. The surfzone edge was computed from alongshore-averaged cross-shore profiles of the wave height estimated from lidar measurements (Baker , 2023).

### B. Remotely sensed surface velocities

In the surf zone, foam and bubbles generated by wave breaking result in bright features in visible-band imagery that contrast with darker non-breaking regions. The residual foam following an active breaker advects with surface currents, often imparting a rich signal to track mean alongshore currents (Chickadel , 2003; Perkovic , 2009; and Almar , 2016) and mean velocities (Anderson , 2021). However, remotely sensing velocities at low frequencies is challenging in the surf zone, partially because bright breaking waves propagating through the field of view often lead to an onshore bias. Here, we estimated highly resolved horizontal, surface velocities with a feature-based particle image velocimetry (PIV) technique developed by Chickadel (2011). This method tracks shifts in coherent surface texture (tracers such as bubbles generated during active breaking, residual foam, and white plastic, perforated, hollow golf balls) across a sequence of visible-band images, rectified in laboratory coordinates.

PIV was performed on the 8 Hz, 1 cm resolution rectified imagery from the camera c2 over the width of the tank ( $ y = \u2212 13.3$–13.3 m) and from the shoreline to just offshore of the surf zone ( $ x = 25 \u2212 33.5$ m). Following Chickadel (2011), the PIV technique uses a phase correlation between image pairs to better refine feature motion given large-scale background variation (Adrian, 2005). Sub-pixel refinement was achieved by computing the peak of a 2D Gaussian fit to the phase map (Eckstein , 2008). A multi-pass approach with decreasing interrogation windows and 50% overlap started with 256 × 256 pixels, followed by 128 × 128 pixels, 64 × 64 pixels, 32 × 32 pixels, and culminated with 16 × 16 pixels. This resulted in an 8 Hz timeseries of cross- (u) and alongshore (v) velocities with an effective spatial resolution of 8 cm in the cross- and alongshore direction. Velocity estimates were omitted if they exceeded known physical constraints (|*u*|>4 m/s, |*v*|>2 m/s) or if the signal to noise ratio fell below 0.68, indicative of insufficient image contrast for optimal PIV performance (Chickadel , 2011).

Surface PIV velocities were comprised of processes spanning frequency and spatial scales, including breaker celerity and low-frequency horizontal motions. We were interested in the low-frequency surface velocities ( $ f \u2264 0.12$ Hz) to study surfzone eddy processes. Therefore, we removed instantaneous surface velocities at locations of breaking crests to reduce contamination by surface gravity waves on the low-frequency velocity signal. The breaking crests were identified with the Remote Breaker Identification Scheme (RBIS; Baker , 2023) and enlarged by 20 cm on all sides. RBIS detects active breaking by selecting foamy pixels that exceed an intensity threshold and removing regions that are below a water-surface-elevation threshold and likely residual foam. The remaining surface velocities were filtered with a time-space moving median filter over 8 s and 24 cm in the cross- and alongshore direction (nine-point spatial median). This filtering technique was chosen to optimize the consistency between surface and in situ alongshore velocities while preserving velocities near vorticity injection timescales. Grid locations with PIV estimates less than 50% of the time for a given filter window were omitted. After point removal and filtering, surface velocities were retained $ > 95 %$ of the time for the inner surf zone and >85% of the time at the outer edge of the surf zone. For purposes of the turbulence statistics that require homogeneous flow, we removed measurements at locations where the standard deviation of the alongshore eddy velocity (σ_{v}) was less than 20% of the alongshore-averaged σ_{v}. These lower σ_{v} are due to stationary sub-surface features that were visible through the water surface (e.g., small creases between concrete bar sections and metal plates used to attach *in situ* sensor bodies). Vorticity was estimated from gridded low-frequency surface velocity components using a central-difference approach (Patankar, 1980). Additionally, *in situ* observations were median filtered over 8 s for comparisons between surface and *in situ* velocities (Sec. III B).

Despite removing velocities associated with breaking waves, uncertainties in low-frequency velocity estimates may persist due to incomplete removal of breaking crests as well as particles not following with the surface flow. Median filtered surface velocities were less consistent with *in situ* velocities when crests were not removed, resulting in twice the onshore bias for cross-shore velocities and larger root mean squared differences (RMSDs). When the distance around extracted breaking crests was expanded from 20 to 30 cm, surface velocities were slightly more consistent with *in situ* velocities in the outer surf zone (i.e., RMSDs decreased by $ < 5 %$), while PIV-*in situ* velocity comparisons did not significantly change in the inner surf zone. This very small improvement in PIV-*in situ* velocity comparisons did not compensate for the reduction in the number of PIV estimates. Lowering the threshold values for the RBIS-identified crests (i.e., selecting more crest features or larger regions) did not improve PIV-*in situ* velocity comparisons. RMSDs between surface PIV and *in situ* velocities were similar when longer median filters were applied (e.g., a 10 s median filter), while they were larger when shorter median filters were used (e.g., a 6 s median filter).

### C. Surfzone eddy statistics

We examined characteristics of low-frequency surfzone flows to assess how they vary with wave forcing characteristics. Here, low-frequencies ( $ f \u2264 0.12$ Hz) were defined as $ 4 \xd7 T p$, consistent with previous literature defining low frequencies (e.g., f = 0.04 Hz for $ T p \u2248 8$ s; Lippmann , 1999; MacMahan , 2010). *In situ* velocities were analyzed to assess frequencies over which flow is primarily rotational as well as the tank seiche modes (see the Appendix). Hereafter, we refer to the 8 s median filtered velocities as $ u = ( u , v )$, the 10 min time-average of **u** as $ u \xaf$, and the eddy velocities as $ u \u2032 = u \u2212 u \xaf$.

*S*) as a function of the alongshore position (

_{p}*y*),

*r*is the spatial lag between observations, $ \u27e8 \xb7 \u27e9$ is the temporal average, and

*p*is the structure function order (here, either second-order,

*S*

_{2}, or third-order,

*S*

_{3}. The second-order structure function can be expanded to

*R*, is a dimensionless function calculated as

_{vv}*R*decreases to zero at long lags where the flow is no longer correlated in homogeneous, isotropic turbulence.

_{vv}Turbulence statistics were computed from *in situ* [Fig. 2(a), circles and triangles] and surface $ v \u2032$ from 15 to 25 min. Surface-velocity *S*_{2} and *S*_{3} were computed from $ \u2212 8 < y < 8$ m at intervals of 8 cm in the cross-shore from $ x = 28.7 \u2212 31$ m [Fig. 2(a), region encompassed by the white dashed rectangle]. Constraining the structure function analysis region to the center of the tank reduced signals associated with the concrete side walls, such as near-wall pressure gradients that led to weak mean circulation patterns from $ \u2212 10 > y > 10$ m (Fig. 2) that could alter eddy processes. Moreover, the cross-shore extent was limited to offshore of the region dominated by infragravity wave processes and strong mean cross-shore velocities (*x* > 31 m) and shoreward of locations where sub-surface static features visible through the water column biased surface velocities low in a regularly spaced alongshore pattern (*x* < 28.7 m, see Sec. III B). For the 10 min averaged structure functions, the effective samples per lag (accounting for the median filter) for *r* = 0.24 to 16 m varied from 5025 to 75 for structure functions calculated at each cross-shore position and from 45 225 to 675 for the surfzone-average structure functions.

By dimensional arguments (Kraichnan, 1967), *S*_{2} scales with $ \u223c ( \epsilon r ) 2 / 3$ for the inverse energy cascade inertial region [Fig. 3(a)], where ε is proportional to the energy flux, $ \epsilon = \u2212 d E / d t$, with kinetic energy per unit mass, $E$ [Fig. 3(a)]. In the absence of mean circulation patterns or other processes altering large-scale rotational motions, *S*_{2} asymptotes to twice the velocity variance ( $ 2 \u27e8 v \u2032 2 \u27e9$), because the velocities are uncorrelated above the largest length scale of the inverse energy cascade. The surfzone-average *S*_{2} (average from $ x = 28.7 \u2212 31$ m) were fitted using the least squares method with $ C r \gamma $ from $ r = 1 \u2212 3$ m to evaluate the consistency of the surfzone flow with an inverse cascade, where the exponent (γ) and constant (*C*) vary freely per trial. Additionally, we estimated the energy flux for each trial as $\epsilon \u223c C 3 / 2$ assuming γ = 2/3. We emphasize that this energy flux should be considered the relative, rather than an absolute, energy flux, as we do not assume a value for the universal Kolmogorov constant [typically $ O$(1)]. The fitted lag range was selected to be longer than lags impacted by median filtering and within the inverse energy cascade inertial range for all trials (i.e., three times the surfzone-average water depth to less than the surfzone width, which may be limiting length scales of the inverse energy cascade).

Consequentially, *S*_{2} also scales with $ r 2 / 3$ for the direct energy cascade inertial range in 3D turbulence. Therefore, the energy cascade direction cannot be determined with S2. Conveniently, the sign of the third-order structure function [*S*_{3}, Fig. 3(b)] indicates if the turbulence is consistent with a 2D inverse energy cascade ( $ S 3 > 0$) or 3D direct energy cascade ( $ S 3 < 0$; Xia , 2008; Cerbus and Chakraborty, 2017; Alexakis and Biferale, 2018; and Xie and Bühler, 2019). In a 2D inverse energy cascade, *S*_{3} increases linearly with lag ( $ S 3 \u223c 3 / 2 \epsilon r$), while decreases linearly with lag (*S*_{3} ∼ $ \u2212 4 / 5 \epsilon r$) for a 3D direct energy cascade. For calculations herein, the surfzone-average *S*_{3} (average from $ x = 28.7 \u2212 31$ m) was fitted using the least squares methods with *Cr* over lags at which *S*_{3} is positive and linearly increasing from either *r* = 2 − 4 m or *r* = 3–5 m. Third-order structure functions have been used to oceanographic applications (Poje , 2017; Balwada , 2022).

*l*

_{0}, a measure of the extent of the region where velocities are correlated) was calculated as

*a*= 6 m [where $ R v v ( r ) < 0.15$ for all trials] because the velocity correlation function [ $ R v v ( r )$] did not simply asymptote to zero for all trials. Trends were similar when

*l*

_{0}was calculated as the lag when $ R v v ( r )$ decreases to below 0.2.

We computed *S*_{2} and *S*_{3} in the alongshore with low-frequency alongshore (longitudinal) velocities, as opposed to cross-shore (transverse) velocities, consistent with Spydell and Feddersen (2009), Elgar and Raubenheimer (2020), and Elgar (2023). Structure functions were computed in the alongshore to ensure that longer lags may be resolved, which were not feasible in the cross-shore due to the insufficient foam to estimate surface velocities offshore of the surf zone and inhomogeneity of flows in the inner shelf relative to the surf zone. Structure functions computed from alongshore velocities minimized any influence of cross-shore seiching patterns on the velocity signal. Regardless, cross-shore seiching patterns at low-frequencies were small for directionally spread waves, and there was no evidence of alongshore seiching patterns (the Appendix). Finally, alongshore velocity wavenumber spectra are not presented, because the maximum length scales at which the flows were correlated were similar to alongshore length of the analysis region, thus inhibiting repeated signals that are necessary to adequately resolve the alongshore wavenumbers of the eddy field.

## IV. RESULTS

### A. *In situ* and surface velocity temporal variability

Low-frequency (8 s median filtered) surface velocities are compared with *in situ* velocities. These comparisons are for higher-frequency surface velocity estimates than previous surfzone studies (i.e., at timescales on the order of 4 waves; Chickadel , 2003; Anderson , 2021; and Elgar , 2023). Alongshore surface velocities agree reasonably well with near-bottom *in situ* velocities [Figs. 4(b), 4(d), 5(b), and 5(d)], while cross-shore surface velocities are less consistent with near-bottom *in situ* velocities [Figs. 4(a), 4(c), 5(a), and 5(c)]. Generally, surface and *in situ* velocities are better correlated for directionally spread waves relative to unidirectional waves as well as in the inner surf zone relative to the outer surf zone (Table II). Surface and *in situ* 1 min averaged velocities are more correlated (*r*) than the 8 s median filtered velocities (e.g., 1-min averaged velocity correlations for the inner surf zone: $ r u =$ 0.69, $ r v =$ 0.93 and outer surf zone: $ r u =$ 0.48, $ r v =$ 0.72 for all *in situ* sensors in trial *G2d* and Fig. 5, black circles) with little bias for the alongshore velocities (e.g., 1 min averaged velocity bias: <0.01 m/s). *In situ* and surface alongshore velocity auto-spectra are not statistically different, while *in situ* and surface cross-shore velocity auto-spectra are only statistically different at $ f \u2248 0.04$ Hz for inner surfzone velocities and *f* < 0.015 Hz for outer surfzone velocities (not shown).

Measurement . | Statistic . | G1a
. | G1c
. | G1d
. | G2d
. |
---|---|---|---|---|---|

Onshore u | In situ σ (m/s) _{u} | 0.08 | 0.07 | 0.07 | 0.06 |

Surface σ (m/s) _{u} | 0.07 | 0.11 | 0.11 | 0.10 | |

RMSD (m/s) | 0.10 | 0.11 | 0.10 | 0.08 | |

r | 0.73 | 0.43 | 0.52 | 0.58 | |

s | 0.61 | 0.64 | 0.81 | 0.87 | |

Bias (m/s) | 0.08 | 0.00 | 0.03 | 0.03 | |

Onshore v | In situ σ (m/s) _{v} | 0.03 | 0.06 | 0.07 | 0.07 |

Surface σ _{v} | 0.02 | 0.10 | 0.11 | 0.08 | |

RMSD (m/s) | 0.04 | 0.08 | 0.08 | 0.05 | |

r | -0.07 | 0.61 | 0.62 | 0.82 | |

s | -0.03 | 0.91 | 0.96 | 0.97 | |

Bias (m/s) | 0.00 | 0.02 | 0.00 | 0.00 | |

Offshore u | In situ σ (m/s) _{u} | 0.08 | 0.10 | 0.10 | 0.09 |

Surface σ (m/s) _{u} | 0.08 | 0.09 | 0.08 | 0.06 | |

RMSD (m/s) | 0.11 | 0.11 | 0.11 | 0.11 | |

r | 0.33 | 0.50 | 0.42 | 0.41 | |

s | 0.33 | 0.45 | 0.34 | 0.28 | |

Bias (m/s) | 0.06 | 0.05 | 0.04 | 0.07 | |

Offshore v | In situ σ (m/s) _{v} | 0.02 | 0.04 | 0.04 | 0.03 |

Surface σ (m/s) _{v} | 0.02 | 0.06 | 0.06 | 0.04 | |

RMSD (m/s) | 0.03 | 0.05 | 0.05 | 0.04 | |

r | 0.19 | 0.53 | 0.56 | 0.53 | |

s | 0.21 | 0.80 | 0.90 | 0.58 | |

Bias (m/s) | 0.00 | 0.00 | 0.00 | 0.00 |

Measurement . | Statistic . | G1a
. | G1c
. | G1d
. | G2d
. |
---|---|---|---|---|---|

Onshore u | In situ σ (m/s) _{u} | 0.08 | 0.07 | 0.07 | 0.06 |

Surface σ (m/s) _{u} | 0.07 | 0.11 | 0.11 | 0.10 | |

RMSD (m/s) | 0.10 | 0.11 | 0.10 | 0.08 | |

r | 0.73 | 0.43 | 0.52 | 0.58 | |

s | 0.61 | 0.64 | 0.81 | 0.87 | |

Bias (m/s) | 0.08 | 0.00 | 0.03 | 0.03 | |

Onshore v | In situ σ (m/s) _{v} | 0.03 | 0.06 | 0.07 | 0.07 |

Surface σ _{v} | 0.02 | 0.10 | 0.11 | 0.08 | |

RMSD (m/s) | 0.04 | 0.08 | 0.08 | 0.05 | |

r | -0.07 | 0.61 | 0.62 | 0.82 | |

s | -0.03 | 0.91 | 0.96 | 0.97 | |

Bias (m/s) | 0.00 | 0.02 | 0.00 | 0.00 | |

Offshore u | In situ σ (m/s) _{u} | 0.08 | 0.10 | 0.10 | 0.09 |

Surface σ (m/s) _{u} | 0.08 | 0.09 | 0.08 | 0.06 | |

RMSD (m/s) | 0.11 | 0.11 | 0.11 | 0.11 | |

r | 0.33 | 0.50 | 0.42 | 0.41 | |

s | 0.33 | 0.45 | 0.34 | 0.28 | |

Bias (m/s) | 0.06 | 0.05 | 0.04 | 0.07 | |

Offshore v | In situ σ (m/s) _{v} | 0.02 | 0.04 | 0.04 | 0.03 |

Surface σ (m/s) _{v} | 0.02 | 0.06 | 0.06 | 0.04 | |

RMSD (m/s) | 0.03 | 0.05 | 0.05 | 0.04 | |

r | 0.19 | 0.53 | 0.56 | 0.53 | |

s | 0.21 | 0.80 | 0.90 | 0.58 | |

Bias (m/s) | 0.00 | 0.00 | 0.00 | 0.00 |

The standard deviation of surface eddy velocities ( $ \sigma u , \sigma v$, standard deviation of 8 s median filtered velocities) are on average ∼0.02 m/s faster than near-bottom eddy velocities in the inner surf zone [Figs. 6(a), 6(c), and 6(f)]. Mean surface alongshore velocities are similar to near-bottom velocities, while mean cross-shore velocities are faster near-bottom than at the surface [Figs. 6(b), 6(d), and 6(g)].

Discrepancies between surface velocities and near-bottom *in situ* velocities may not necessarily be indicative of PIV estimate errors. More consistent surface and near-bottom *in situ* alongshore velocities than cross-shore velocities [Table II, compare Figs. 5(a) and 5(c) vs 5(b) and 5(d)] may result from an onshore bias in cross-shore velocities introduced by bright breaking waves or a sheared water column in the cross-shore direction. Sheared cross-shore flows, for example, due to onshore mass flux near the surface balanced by return flow at the bottom, may explain the stronger offshore-directed flows measured with near-bottom *in situ* sensors compared to surface velocities, particularly at the offshore sensors near the bar crest [e.g., *u* slope = 0.28 for *G2d*, Fig. 5(c)]. Higher correlations in the inner surf zone compared to the outer surf zone [compare Figs. 5(a) and 5(b) vs 5(c) and 5(d)] may be due to the more dense, temporally persistent foam in the inner surf zone and higher signal-to-noise ratios [e.g., stronger alongshore fluctuations in Figs. 5(b) vs 5(d)]. This is further supported by the more similar surface and *in situ* velocities in the outer surf zone for larger waves (e.g., *G1d*) where breaking (and thus, foam generation) starts farther offshore. Due to lack of persistent foam (inhibiting tracking), correlations between *in situ* and surface estimates are lower for unidirectional waves. Alongshore velocity fluctuations during unidirectional waves are small ( $ \sigma v = \xb1 0.05$ m/s) and on the order of PIV-estimate uncertainties. Additionally, offshore of the bar crest, surface velocities may be biased offshore due to intermittent estimates in wave troughs only feasible during large waves ( $ x \u2264 28$ m, $ | \u27e8 u \u27e9 | \u2264 0.1$ m/s, std. $ u \u2248 0.14 \u2212 0.20$ m/s). Structure functions are not presented for these cross-shore locations ( $ 31 < x < 28.7$ m).

### B. Surfzone flow spatial patterns

Mean surface currents in the surf zone [ $ 28.7 < x < 31$ m, away from side walls, Figs. 6(b) and 6(d)] for directionally spread waves ( $ \sigma \theta \u2265 10 \xb0$) are much weaker than eddy velocities [Figs. 6(a) and 6(c)], particularly for the alongshore velocity component. Mean alongshore currents are similar to alongshore eddy velocities for unidirectional waves [ $ v \xaf \u2248 \sigma v \u2248 0.02$ m/s, Figs. 6(a) and 6(b)]. Mean alongshore currents are small and similar across the surf zone [ $ v \xaf < 0.03$ m/s, Fig. 6(b)], except for near the swash zone ( $ v \xaf < 0.05$ m/s). Mean cross-shore velocities are offshore directed in the regions with the strongest breaking (near the bar crest, $ 27 < x < 29$ m, and beach, x > 30.5 m) and onshore directed within the trough [ $ 29 < x < 30.5$ m, Fig. 6(d)]. In addition to small, primarily cross-shore, mean currents in the center of the tank, an onshore directed mean current is evident near side walls in surface velocities (e.g., $ u \xaf \u2264 0.1$ m/s for G2d, Fig. 6). Additionally, cross-shore seiche modes at low frequencies were small for directionally spread waves, and there was no evidence of alongshore seiche modes (see more in the Appendix).

Low-frequency velocity fields are primarily rotational for directionally spread waves (the Appendix). Alongshore eddy velocities increase shoreward from $ \sigma v \u2248 0.05$ m/s in the outer surf zone to $ \sigma v \u2248 0.15$ m/s near the swash zone [Fig. 6(a)]. Cross-shore eddy velocities are fairly uniform across the surf zone [ $ \sigma u \u2248 0.10 \u2212 0.15$ m/s, Fig. 6(c)]. The large σ_{u} for unidirectional waves without the presence of rotational flows may be due to a lack of foam abundance (necessary for PIV) and irrotational cross-shore velocity fluctuations at wave set timescales. Eddy velocities are fairly homogeneous along the surf zone [e.g., Fig. 6(g)], where the surfzone-average σ_{v} varies in the alongshore by <0.020 m/s for trials *G1a-d* and <0.035 m/s for trials *G2a-e*. On average across the surf zone, the alongshore variability in σ_{v} is about twice that of $ v \xaf$.

The low-frequency velocity fields during directionally spread waves include rotational, coherent eddy structures with varying spatiotemporal persistence [Figs. 7(d), 7(f), and 8(b)–8(d)], while coherent structures are largely absent for unidirectional waves [Fig. 7(b)]. Eddies exhibit differing length scales for moderate ( $ \sigma \theta = 10 \xb0$) and large ( $ \sigma \theta = 24 \xb0$) directional spreads, and can advect along the surfzone, occasionally interacting with other coherent structures. White perforated, hollow golf balls and foam that collect in convergence zones were occasionally ejected offshore in a transient rip current [e.g., Figs. 7(c) and 7(d) at *y* = −7 m, 7(e) and 7(f) at *y* = −4 m, 8(a) and 8(b) at *y* = −8 m].

### C. Surfzone flow consistency with two-dimensional turbulence

Second-order structure functions (*S*_{2}) were computed from low-frequency surface and *in situ* alongshore velocities to quantify the spatial scales of coherent eddy structures and evaluate the consistency of surfzone flow with an energy cascades to smaller or larger length scales (Fig. 9, Sec. III C). A power law scaling ( $ C r \gamma $) was fit to *S*_{2} from lags of $ r =$ 1–3 m, where larger C suggests higher energy flux (ε) and larger γ indicates more rapid spatial decorrelation. For an inverse energy cascade, $ C \u223c \epsilon 2 / 3$ and $ \gamma = 2 / 3$ (Sec. III C). The structure function is nearly independent of cross-shore position for unidirectional waves [*G2a*, Fig. 9(a), circles] and for directionally spread cases at *r* < 1.5 m [Figs. 9(b) and 9(c)]. For directionally spread waves, γ and ε are larger in the inner surf zone than the outer surf zone at *r* > 1.5 [ $ \sigma \theta \u2265 10 \xb0$, Figs. 9(b) and 9(c), differences between blue and light green circles].

The shape of the surfzone-average *S*_{2} from $ x = 28.7 \u2212 31$ m vary with directional spread (Fig. 9, black curves) and are similar to *S*_{2} computed over the trough region ( $ x = 28.7 \u2212 31$ m, blue circles). The surfzone-average *S*_{2} is nearly flat ( $ \gamma = 1 / 20$) at r > 0.6 m for unidirectional waves, suggesting that velocities are uncorrelated at longer lags. In contrast, the surfzone-average *S*_{2} is correlated from $ r = 1 \u2212 3$ m (positive slope) for directionally spread waves ( $ \sigma \theta = 10 \xb0 , \u2009 24 \xb0$). For these directional spreads, *S*_{2} is nearly consistent with the energy cascade self-similar power law of $ ( r ) 2 / 3$ [Figs. 9(b) and 9(c), red lines]. Based on the best fit curves, *S*_{2} scales with $ r 0.68$ and $ r 0.69$ for trials with $ \sigma \theta = 10 \xb0$ and $ 24 \xb0$, respectively. The scaling exponent (γ) varies by ±0.1 depending on the spatial and temporal scales of the surface velocity filter (e.g., removing the spatial median filter or increasing the temporal median filter to 10 s) and the range of lags over which the power laws were fit (e.g., 0.5–4 m). Additionally, the velocity decorrelation length scale (*l*_{0}, the characteristic size of large eddies) decreases from around 3 to 2.5 m for $ \sigma = 10 \xb0$– $ 24 \xb0$. Based on testing with a white noise signal, the spatial median applied to surface velocity estimates increases the correlation at small lags, altering the power law at $ r \u2264 0.24$ m. The temporal filter likely also increases the correlation at small lags (e.g., for a 0.1 m/s flow, an 8 s median filter may increase the correlation at *r* < 0.8 m). Thus, we are unable to make definitive conclusions regarding the minimum length scales of an inverse cascade with *S*_{2}.

On average, *S*_{2} calculated from *in situ* sensors in the inner surf zone (*x* = 30.7 m) are three times larger than *S*_{2} from *in situ* sensors in the outer surf zone (*x* = 28.4 m). For directionally spread waves, the magnitude of *in situ* *S*_{2} is smaller at most lags relative to the surface velocity *S*_{2}, particularly in the outer surf zone (Fig. 9, compare triangles and circles). Despite their magnitude differences, the *in situ* *S*_{2} from the inner surf zone (*x* = 30.7 m) are similarly shaped to the surfzone-average surface velocity *S*_{2} (Fig. 9, compare green triangles with black curve). For moderate directional spreads ( $ \sigma \theta = 10 \xb0$), the inner surf zone *in situ* *S*_{2} follows $ \u223c r 0.8$ between $ r = 3.2 \u2212 6.4$ m, similar to the surfzone-average surface velocity *S*_{2} between these lags. In contrast, for large directional spreads ( $ \sigma \theta = 24 \xb0$), the inner surf zone *in situ* velocity *S*_{2} between $ r = 3.2 \u2212 6.4$ m is uncorrelated (nearly flat), similar to the surfzone-average surface velocity *S*_{2} at *r* > 4 m. The *in situ* velocity *S*_{2} from the outer surf zone (*x* = 28.4 m) are similarly shaped to the outer-most surface velocity *S*_{2} (Fig. 9, compare blue triangles with dark blue circles).

The surfzone-average *S*_{2} from surface velocities are compared for all trials to assess the consistency of the patterns observed in trials *G2a,b,d* for other wave heights ( $ H s \u2248 0.24$, 0.27 m) and directional spreads [ $ \sigma \theta = 2 \xb0 \u2212 26 \xb0$, Table I, Fig. 10(a)]. For unidirectional waves with both wave heights, *S*_{2} is uncorrelated at *r* > 0.6 m (*G1a*, *G2a*, Fig. 10, black curves). For trials with $ \sigma \theta \u2265 10 \xb0$, *S*_{2} approximately satisfy the energy cascade power law ( $ \gamma = 0.66 \xb1 0.1$, red line, Fig. 10) and the energy flux (ε) is slightly lower for smaller waves, where differences between curves are significant ( $ > 95 %$ confidence interval). Although still consistent with $ r 2 / 3$, *S*_{2} from trial *G2b* with $ \sigma \theta = 10 \xb0$ has overall lower energy flux, inferred from the magnitude of *S*_{2} (Fig. 10, dashed purple curve) than trials with larger directional spreads from $ r = 0.6 \u2013 5$ m. Additionally, the surfzone flow fields generally remain more correlated at longer lags for trials with smaller spreads [Fig. 10(b)], except for unidirectional waves. For example, trial G2*b* ( $ \sigma \theta = 10 \xb0$) does not decrease to *R _{vv}* = 0.5 until

*r*= 3.2 m, relative to the trial G2

*d*, where

*R*= 0.5 at

_{vv}*r*= 1.4 m. Most of the trials decrease to $ R v v \u2248 0.1$ by

*r*= 7.5 m, except for trial

*G*1

*d*, that remains correlated ( $ R v v \u2248 0.15$) at all lags. Second-order structure functions computed from transverse velocities exhibit similar trends.

The third-order structure function from alongshore velocities (*S*_{3}) vary with cross-shore position and directional spread (Fig. 11). Consistent with *S*_{2}, the energy flux (slope of *S*_{3}) is larger in the inner than outer surf zone for directionally spread waves [Figs. 11(b) and 11(c), compare light green vs dark blue]. For unidirectional waves, *S*_{3} is approximately zero at all lags, suggesting no energy flux is present, consistent with *S*_{2}. The surfzone-averaged *S*_{3} is negative for trial G2b ( $ \sigma \theta = 10 \xb0$) at nearly all lags suggesting that energy is primarily transferred to smaller scales. In contrast, for trial G2d, the surfzone-averaged *S*_{3} transitions from negative over $ r = 1 \u2212 2.6$ m to positive over $ r = 2.6 \u2212 5$ m with an approximately linear slope.

These trends in *S*_{3} are similar for different wave heights (Fig. 12). The surfzone-averaged *S*_{3} is near zero for unidirectional wave cases, while the surfzone-averaged *S*_{3} for trials with large directional spread ( $ \sigma \theta > 10 \xb0$) are often negative at *r* < 1 m and are positive at *r* > 2 m. For large directional spreads, the rate of linear increase in *S*_{3} is greater for trials with larger wave heights as expected due to the dependence of *S*_{3} on energy flux (i.e., $ 2 3 \epsilon r$). Third-order structure functions computed from both the transverse and longitudinal velocities [i.e., $ S 3 = \u27e8 v ( r ) ( v ( r ) 2 + u ( r ) 2 ) \u27e9$] also include a mixed forward ( $ S 3 < 0$) and inverse ( $ S 3 > 0$) cascade for large directional spread cases. The only wave conditions where the surfzone-averaged *S*_{3} is negative over all lags is for moderate directional spread ( $ \sigma \theta = 10 \xb0$).

## V. DISCUSSION

In a forced 2D turbulence system, energy is non-linearly transferred from vorticity injection length scales (*l _{i}*) to larger scales (

*l*

_{0}, see Sec. II). An inverse energy cascade from short-crested breaking wave injected vorticity may explain the presence of large-scale rotational motions found in numerical models and observations (Peregrine, 1999; Bonneton , 2010; Spydell and Feddersen, 2009; Feddersen, 2014; Elgar and Raubenheimer, 2020; and Marchesiello , 2022). Remotely sensed field observations demonstrated that third-order structure functions (

*S*

_{3}) were consistent with an inverse energy cascade for a narrow range of directional spreads ( $ \sigma \theta = 24 \xb0 \u2212 29 \xb0$; Elgar , 2023). Additionally, modeled wave number spectral fluxes were consistent with energy transferred to larger scales for directionally spread wave conditions ( $ \sigma \theta = 30 \xb0$; Marchesiello , 2022). Here, we add to existing literature by revealing how the inverse energy cascade assessed with observations varies with wave directional spread, which alters the breaker crest lengths and the number of crest ends that generate vorticity (Baker , 2023).

We find that the surfzone-average, second-order structure function (*S*_{2}) generally conforms with the expected power law for energy cascades to larger/smaller scales ( $ \gamma = 2 / 3$) for directionally spread wave conditions from $ \sigma \theta = 10 \xb0 \u2013 26 \xb0$ (Fig. 10). Conversely, *S*_{2} for unidirectional waves ( $ \sigma \theta = 2 \xb0$) does not scale with $ \gamma = 2 / 3$. Unidirectional waves lack random wave-generated vorticity by alongshore-variable wave breaking and resulted in significantly less eddy activity than trials with directionally spread waves [e.g., Fig. 13(b)]. Notably, γ is slightly larger than 2/3 for trials with larger wave heights and wider surf zones. The γ of *S*_{2} increases shoreward for all directionally spread trials. Thus, the larger γ for trials with larger wave heights may be due to the surfzone-averaged *S*_{2} not accounting for changes in the surfzone width. These results are generally in agreement with results from the inner surf zone of phase-resolved, depth-averaged numerical simulations of shore-normal waves over alongshore-uniform bathymetry, where surfzone flow was consistent with an energy cascade for $ \sigma \theta = 14 \xb0$ and $ 20 \xb0$ (Spydell and Feddersen, 2009).

Energy flux (*ε*), estimated from fitting $ ( \epsilon r ) 2 / 3$ to *S*_{2}, is larger in the inner surf zone than in the outer surf zone (Fig. 9, compare light vs dark) and increases with increasing directional spread and wave height [Fig. 13(e)]. Larger directional spreads result in more crest ends (Baker , 2023) where vorticity may be injected (Peregrine, 1998; Clark , 2012). Additionally, due to the dependence of the breaking force on wave dissipation ( $ F br \u2009 \u221d \u2009 D$; Clark , 2012), larger wave heights with the same directional spread (thus, the same average along-crest length) may increase injection rates per crest ( $ \u2207 \xd7 F br = \u2202 \omega / \u2202 t$; Spydell and Feddersen, 2009; Feddersen, 2014). The increase in ε shoreward and with directional spread is consistent with previous depth-averaged numerical simulations (Spydell and Feddersen, 2009).

The alongshore eddy velocities, characterized as the surfzone-averaged standard deviation of v [ $ \u27e8 \sigma v \u27e9 sz$, Fig. 13(b)], are independent of directional spread and increase with wave height, except for unidirectional waves. Conversely, the decorrelation length scales (*l*_{0}, characteristic eddy length scale) weakly decrease with increasing wave directional spread (from $ 2.9 \u2212 2.4$ for $ \sigma \theta = 10 \xb0 \u2212 26 \xb0$) and do not exhibit strong dependence on wave height [Fig. 13(c)]. The decorrelation length scales are smaller than the surfzone width (*W*_{sz}). The weak length-scale dependence on directional spread could be due to longer waves forcing larger-scale vortical motions (Baker, 2023). Based on field observations, the maximum length scales may be set by the surfzone width (4.4–4.9 m in the lab; Elgar , 2023) as low-frequency eddy activity is significantly smaller just outside the surf zone than within it (Kennedy, 2005; Spydell , 2007; Spydell and Feddersen, 2009; and Feddersen, 2014).

We assessed the direction of the energy cascade with third-order structure functions since this cannot be determined from second-order structure functions unless the forcing length scales are well established. Third-order structure functions (*S*_{3}) suggest that energy fluxes can be bidirectional, a mixture of forward ( $ S 3 < 0$) and inverse energy cascades ( $ S 3 > 0$), over lags where $ S 2 \u223c r 2 / 3$. These energy flux directions depend on wave conditions and cross-shore position (Figs. 11 and 12). For large directional spread waves, the surfzone-averaged *S*_{3} are positive and linearly increase often starting at $ r \u223c 2$ m, consistent with an inverse energy cascade to large scales. These results are generally consistent with alongshore third-order structure functions in field observations (Elgar , 2023). Energy flux is estimated by fitting *S*_{3} with $ 2 / 3 \epsilon r$ to lags where *S*_{3} is positive and linearly increasing (from either $ r = 2 \u2013 4$ m or $ r = 3 \u2013 5$ m). Therefore, ε from *S*_{3} is only related to energy flux to larger scales, while ε from *S*_{2} may include energy flux in both directions. Similar to energy flux estimated from *S*_{2}, ε increases with wave height and is about twice the magnitude of ε from *S*_{2} for trials with $ H s \u2248 0.24$ m. The exception to this is for $ \sigma \theta = 10 \xb0$ as *S*_{3} was negative at all lags (fit from $ r = 2 \u2013 4$ m). Additionally, ε from *S*_{3} for larger wave heights (trials *G1a-d*) does not increase with wave height. This may suggest that energy transfer to larger scales is limited for a given wave height. If there is too much random forcing in a two-dimensional turbulence system, formation of large-scale coherent eddies may be disrupted, limiting the likelihood of an inverse energy cascade (Kellay and Goldburg, 2002) at larger wave heights. Furthermore, similar to the less consistent *S*_{2} scaling for larger waves, the analysis does not compensate for the changes in surfzone width, and therefore, the high-energy inner surf zone may occupy a relatively greater portion of the fixed analysis region for trials with larger waves (*G1a–d*).

Unlike larger directional spreads, the surfzone-averaged *S*_{3} for the moderately directional spread waves ( $ \sigma \theta = 10 \xb0$) is negative at all lags, but the sign varies by cross-shore position. This may suggest energy is not primarily transferred to longer length scales as inferred from *S*_{3} for trials with larger directional spreads. Eddy forcing may occur at large length scales for the moderate directional spread waves, leading to the larger decorrelation length scale [Fig. 13(c)].

In an irregular, directionally spread wave field, vorticity injection may occur at a broad distribution of length scales (Fig. 3, yellow region; Feddersen, 2014) that vary with wave field characteristics. Thus, we anticipate that the vorticity injected by random short-crested breaking may vary with directional spread (and therefore, crest length and ends; Baker , 2023), wave height, and period based on previous numerical modeling and field observations (Clark , 2012; Suanda and Feddersen, 2015). Eddy injection at longer length scales may occur if the input scales are proportional to the breaking crest length (Baker, 2023). If the maximum forcing scale is $ \lambda c / 2$ and the upper limit of an inverse energy cascade is the surfzone width ( $ W sz$), then, when $ \lambda c / W sz \u2248 2$, forcing may occur up to decorrelation length scale. This may possibly explain the lack of a clear inverse energy cascade for $ \sigma \theta = 10 \xb0$, where $ \lambda c / W sz = 1.6$. However, future work could investigate the relationships between the wave conditions, vorticity injection length scales, and energy fluxes at moderate directional spreads.

There are several limitations to assessing the surfzone flow with structure function power laws. The energy flux and Kolmogorov constant may not be constant (not accounted for here, Landau and Lifshitz, 2013). The power law scaling neglects local fluctuations (intermittency), which may be ubiquitous within the surf zone. The power law scaling is only valid in local inertial ranges far from the forcing and dissipation scales, which may be violated for moderate directional spread cases. Additionally, the surf zone can be non-homogeneous and non-isotropic in the cross-shore (Fig. 6). Furthermore, theory for third-order structure functions was developed for unidirectional energy transfer (Cerbus and Chakraborty, 2017; Xie and Bühler, 2019), rather than bidirectional cascades relevant in the surf zone.

In the framework described here, we assume that eddies are primarily 2D with horizontal scales much larger than water depth (Fig. 10). However, the inferred depth-dependence of the structure function computed from surface velocities vs near bottom *in situ* velocities suggests that low-frequency motions may exhibit vertical structure, especially in the outer surf zone where the water column was sheared (Fig. 6, compare surface and *in situ* velocities). The depth-variable low-frequency fluctuations are consistent with observations and numerical simulations that suggested vertical variability in very low-frequency currents is more significant near prominent bar-trough features in the outer surf zone (Lippmann and Bowen, 2016; Henderson , 2017; and Baker , 2021). Additionally, recent numerical modeling studies suggest that the breaking force ( $ F br$) is better represented as an applied surface stress, than a depth-uniform body force, leading to vertical variability in the wave-forced vorticity (Kirby and Derakhti, 2019). Depth-varying wave-forced vorticity may alter the efficiency of an inverse energy cascade based on findings from numerical simulations that resolved 3D processes (Uchiyama , 2017; McWilliams , 2018; and Marchesiello , 2022).

Our results provide the first direct evidence of the dependence of an inverse energy cascade within the surf zone on directional spread. These contributions are key in aiding our understanding of energy transfers across spatial scales within the surf zone, which impacts material transport and dispersion. Building on these observations, future work is needed to assess how findings from experiments translate to environments with varying bathymetric profiles, including planar profiles as studied in previous work (Kumar and Feddersen, 2017a; Suanda and Feddersen, 2015) and alongshore-variable bathymetry with fixed vorticity injection locations (Suarez , 2023). Additionally, future modeling and observational studies could assess the dependence of vorticity dynamics on directional spread outside of the constraints of a finite-width laboratory (see more discussion in Baker , 2023), link the vorticity injection process to the inverse energy cascade, and further explore the three-dimensionality of these processes.

## VI. CONCLUSION

In a forced quasi-two-dimensional turbulence system, vorticity is injected at an intermediate length scale and transferred nonlinearly to both smaller and larger scales. Our results provide evidence that energy associated with vorticity generated by short-crested breaking waves can be transferred to larger and smaller length scales through a bidirectional cascade (inverse and direct energy cascade). We find that the presence of an inverse energy cascade as well as the energy flux and decorrelation length scales vary with wave directional spread and height.

Specifically, we resolve vorticity dynamics at 8 cm resolution across the surf zone of a large-scale directional wave basin with velocities from surface PIV techniques applied to imagery and *in situ* sensors. Alongshore surface velocity fluctuations, estimated at higher frequencies relative to previous surfzone studies (e.g., Chickadel , 2003; Anderson , 2021; and Elgar , 2023), are similar to near-bottom *in situ* velocities. Our findings reveal that the eddy activity (outside of frequencies where cross-shore seiching may alter the signal) is a much larger signal than mean circulation in the laboratory and includes a sheared water column, particularly near the bar crest. We establish that second-order alongshore velocity structure functions are consistent with the expected power law for an energy cascade when waves are directionally spread ( $ \sigma \theta \u2265 10 \xb0$), but not when waves are near unidirectional, where the vorticity injection mechanism is absent. Based on third-order structure functions, we show that surfzone flow is consistent with an inverse energy cascade from intermediate to larger length scales for large directional spreads, but not for moderate directional spreads. Both the second- and third-order structure functions demonstrate that the energy flux varies by surfzone cross-shore region. We also find that the decorrelation length scales weakly decrease, while the energy flux increases, with increasing directional spread. Furthermore, energy flux and eddy velocities vary with wave height.

Our findings contribute to the understanding of mechanisms facilitating energy exchange between scales within the surf zone. Future work could examine the sensitivity of these processes to varying bathymetric profiles. A fundamental understanding of surfzone vorticity dynamics can aid predictions of the dispersion and transport of pollutants, nutrients, and larvae in the nearshore, as these eddies can exit the surf zone, transporting material into deeper waters.

## ACKNOWLEDGMENTS

We thank Nirnimesh Kumar, who both inspired this line of scientific inquiry and co-lead the laboratory campaign. The authors thank the staff of the O.H. Hinsdale Wave Research Laboratory at Oregon State University for helping to run experimental trials and obtain, collect, and quality-control *in situ* measurements. Also, the authors thank Jim Thomson at the University of Washington Applied Physics Laboratory and Steve Elgar at Woods Hole Oceanographic Institute for their input on the data analysis. Support was provided by the National Science Foundation, United States, a National Defense Science and Engineering Graduate Fellowship, and a Burges Fellowship. M.M. was also supported by the ONR Young Investigator Program. M.L.P was funded by the U.S. Naval Research Laboratory Base Funding and the U.S. Geological Survey Coastal-Marine Hazards and Resources Program. K.B. was funded by the U.S. Army Corps of Engineers Coastal and Ocean Data Systems Program and by the U.S. Army Engineer Research and Development Center's Military Engineering Basic Research Program from the Assistant Secretary of the Army for Acquisition, Logistics, and Technology. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Christine M. Baker:** Conceptualization (equal); Data curation (lead); Formal analysis (lead); Funding acquisition (supporting); Investigation (equal); Methodology (lead); Software (supporting); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (lead). **Melissa Moulton:** Conceptualization (equal); Formal analysis (supporting); Funding acquisition (lead); Investigation (equal); Methodology (supporting); Validation (supporting); Visualization (supporting); Writing – review & editing (supporting). **C. Chris Chickadel:** Methodology (supporting); Software (lead); Validation (supporting); Writing – review & editing (supporting). **Emma S. Nuss:** Methodology (supporting); Writing – review & editing (supporting). **Margaret Palmsten:** Investigation (supporting); Methodology (supporting); Writing – review & editing (supporting). **Katherine Brodie:** Investigation (supporting); Methodology (supporting); Writing – review & editing (supporting).

## DATA AVAILABILITY

The data and PIV output that support the findings of this study are openly available in the DesignSafe-CI Data Depot (Baker 2023).

### APPENDIX: ROTATIONAL AND IRROTATIONAL FLOWS IN THE SURF ZONE

The relative contribution of velocity fluctuations explained by irrotational (surface gravity waves, seiche patterns) and rotational (eddies) motion as well as tank seiching modes (standing wave oscillations) was assessed by comparing the magnitudes of the velocity spectrum ( $ S u u + S v v$) and sea-surface elevation spectral levels converted to equivalent velocity using the linear dispersion relationship ( $ S \eta \eta $ EV). The spectra were computed using a Hanning window period of 256 s with an overlap period of 128 s over a 35 min time series (5–40 min, DOF = 29). Velocities are attributed to irrotational motions if $ S u u + S v v \u2248 S \eta \eta $ EV, while velocities are associated with rotational motions or possibly at a seiche node if $ S u u + S v v > S \eta \eta $ EV (Lippmann , 1999; MacMahan , 2010; Feddersen , 2011; and Elgar , 2019).

Here, a representative example is described from a 35 min time series (5–40 min) from a sensor in the outer surf zone ( $ x , y = 28.4 , \u2212 6.2$ m) for unidirectional and directionally spread wave conditions (trials *G2a*,*d*: $ \sigma \theta = 2 \xb0 , \u2009 24 \xb0$, Fig. 14). $ S u u + S v v$ is 2–4 times greater than $ S \eta \eta $ EV at $ 0.10 < f < 0.03$ Hz and more than 10 times larger at $ f \u2264 0.02$ [Fig. 14(a), compare dashed curves]. This suggests that surfzone flow is primarily explained by rotational (eddy-like) motions at low-frequencies for directionally spread waves. These trends are consistent across sensor locations and for trials with smaller directional spreads (e.g., $ \sigma \theta = 10 \xb0$). In contrast, for unidirectional waves, $ S u u + S v v$ is more similar to $ S \eta \eta $ EV at *f* < 0.03 Hz [Fig. 14(a), compare solid curves], suggesting that flow is irrotational (driven by sea-surface elevation fluctuations) and may be attributed to cross-shore seiche patterns.

Cross-shore seiching modes at low-frequencies (*f* < 0.12 Hz) decrease with increasing directional spread [Fig. 14(b), white region] and decreasing wave period (e.g., from *T _{p}* = 4 to 2 s, not shown). Directionally spread waves with

*T*= 2 s do not exhibit significant cross-shore seiching patterns [Fig. 14(b), dashed yellow curves]. For example, the most significant cross-shore seiching mode at

_{p}*f*= 0.076 Hz, apparent for unidirectional waves, is nearly absent for directionally spread waves. Despite this, structure functions were conservatively calculated solely with low-frequency alongshore velocities that show little evidence of significant alongshore seiching modes [Fig. 14(b), red curves].

Observed cross-shore seiche modes at 0.043, 0.076, and 0.116 Hz [Fig. 14(b), vertical gray lines] are at slightly lower frequencies than estimates with the Merian's formula for a rectangular basin with uniform depth (i.e., 0.048, 0.095, and 0.143 Hz for $ n = 1 \u2212 4$, Rabinovich, 2009). As there is a slope starting at $ x =$ 22 m from the wavemaker, Merian's formula may not produce the exact seiching periods (e.g., estimates for a triangular basin with akin geometry have slightly lower frequency seiching modes, Rabinovich, 2009). Notably, the difference between $ S u u + S v v$ and $ S \eta \eta $ EV at *f* = 0.076 Hz for unidirectional waves [Fig. 14(a), compare solid curves] is much smaller for sensors in the inner surf zone (*x* = 30.7 m). This suggests that there may be a cross-shore seiche node near *x* = 28.4 m at *f* = 0.076 Hz and anti-node near *x* = 30.7 m that may alter breaking patterns and currents. Additionally, there was no evidence that seiching patterns grew over time based on Empirical Orthogonal Function (EOF) analysis of *in situ* sensors computed with short time increments over the 45 min trial (not shown).

## REFERENCES

*Environmental Toxicology: Selected Entries from the Encyclopedia of Sustainability Science and Technology*

*Fluid Mechanics: Landau and Lifshitz: Course of Theoretical Physics*

*Numerical Heat Transfer and Fluid Flow*

*Coastal Engineering 1998*

*Handbook of Coastal and Ocean Engineering*