Mean flow structure and velocity–bed shear stress maxima phase difference in smooth wall, transitionally turbulent oscillatory boundary layers: direct numerical simulations

Abstract Direct numerical simulations of oscillatory boundary-layer flows in the transitional regime were performed to explain discrepancies in the literature regarding the phase difference ${\rm \Delta} \phi$ between the bed-shear stress and free-stream velocity maxima. Recent experimental observations in smooth bed oscillatory boundary-layer (OBL) flows, showed a significant change in the widely used ${\rm \Delta} \phi$ diagram (Mier et al., J. Fluid Mech., vol. 922, 2021, A29). However, the limitations of the point-wise measurement technique did not allow us to associate this finding with the turbulent kinetic energy budget and to detect the approach to a ‘near-equilibrium’ condition, defined in a narrow sense herein. Direct numerical simulation results suggest that a phase lag occurs as the result of a delayed and incomplete transition of OBL flows to a stage that mimics the fully turbulent regime. Data from the literature were also used to support the presence of the phase lag and propose a new ${\rm \Delta} \phi$ diagram. Simulations performed for ${\textit {Re}}_{\delta }=671$ confirmed the sensitivity in the development of self-sustained turbulence on the background disturbances ($\textit{Re}_{\delta}=U_{o}\delta/\nu$, where $\delta=[2\nu/\omega]^{1/2}$ is the Stokes' length, $U_{o}$ is the maximum free stream velocity of the oscillation, $\nu$ is the kinematic viscosity and $\omega=2{\rm \pi}/T$ is the angular velocity based on the period of the oscillation T). Variations of the mean velocity slope and intersect values for oscillatory flows are also explained in terms of the proximity to near-equilibrium conditions. Relaminarization and transition effects can significantly delay the development of OBL flows, resulting in an incomplete transition. The shape and defect factors are examined as diagnostic parameters for conditions that allow the formation of a logarithmic profile with the universal von Kármán constant and intersect. These findings are of relevance for environmental fluid mechanics and coastal morphodynamics/engineering applications.


Introduction
In Mier, Fytanidis & García (2021), the experimental observations of mean flow structure and bed-shear stress/free-stream velocity maximum phase difference were presented for the case of intermittently turbulent oscillatory boundary-layer (OBL) flows over smooth walls. A revision of the bed-shear stress/free-stream velocity maxima phase difference diagram was proposed and a threshold value of Re δ = 763 was identified as a critical Re δ value for which phase lag starts being observed (Re δ = U o δ/ν, where δ = [2ν/ω] 1/2 is the Stokes length, U o is the maximum free-stream velocity of the oscillation, ν is the kinematic viscosity and ω = 2π/T is the angular velocity based on the period of the oscillation T). This new diagram explains inconsistencies in the literature regarding the instance when the maximum of the bed-shear stress was predicted with respect to the instance of free-stream velocity maximum. In the present work, flows in the same regime as those of Mier et al. (2021) will be examined in an effort to analyse their characteristics. These results are of relevance for environmental fluid mechanics applications, and coastal engineering and morphodynamics (Sleath 1984;Fredsøe & Deigaard 1992;Nielsen 1992;Garcia 2008;Sumer 2014).
OBL flows can be categorized into different regimes (Akhavan, Kamm & Shapiro 1991a;Pedocchi, Cantero & García 2011;Ozdemir, Hsu & Balachandar 2014), namely: (i) the laminar regime (Re δ < Re δ cr1 ), for which analytical solutions exist for the velocity and shear stress profiles (Batchelor 1967); (ii) the disturbed laminar regime (Re δ cr1 < Re δ < Re δ cr2 ), in which small perturbations are superimposed on the laminar profiles without altering the mean characteristics of the flow such as the mean velocity or shear stress profiles (Carstensen, Sumer & Fredsøe 2010); (iii) the intermittently turbulent regime (Re δ cr2 < Re δ < Re δ cr3 ), for which the flow tends to remain in a quasi-laminar state for part of the acceleration phase until turbulent bursts are observed later during the period (starting at the beginning of the deceleration phase and moving closer to the end of the acceleration phase as Re δ increases), altering both the mean flow velocity profiles and the bed-shear stress signature of the flow (Merkli & Thomann 1975;Hino et al. 1983;Akhavan et al. 1991a;Akhavan, Kamm & Shapiro 1991b); and (iv) the fully turbulent regime (Re δ > Re δ cr3 ) in which high turbulence levels are observed during the whole cycle of the oscillation and the logarithmic layer is valid for most of the time during the oscillation cycle, excluding a period close to the flow reversal (Jensen, Sumer & Fredsøe 1989). Different flow regimes have been identified to alter significantly the temporal variation of mean flow characteristics (Hino, Sawamoto & Takasu 1976;Jensen et al. 1989;Akhavan et al. 1991a;Mier et al. 2021) and bed-shear stress Mier et al. 2021) over the period. Even in the early works of Kajiura (1964), Kamphuis (1975) and Sarpkaya (1993), these researchers had recognized the effect of different flow regimes on the friction coefficient f w (defined as f w = 2τ max /2U 2 o , where τ max is the maximum of the ensemble-average bed-shear stress).
The values of the critical Reynolds numbers Re δ cr1 , Re δ cr2 and Re δ cr3 have been the subject of many studies. Reviews of these efforts can be found in works by Akhavan et al. (1991a,b), Sarpkaya (1993), Blondeaux & Vittori (1994), Blennerhassett & Bassom (2002), Ozdemir et al. (2014) and Thomas et al. (2015). For Re δ cr1 , a value of 85 is commonly accepted (Blondeaux & Seminara 1979;Akhavan et al. 1991a). The reported values for Re δ cr2 in the literature vary. Values of 500-550 are reported experimentally and numerically by Hino et al. (1976) and Jensen et al. (1989). A value of 600 is reported by Vittori & Verzicco (1998). Blennerhassett & Bassom (2002) proposed an Re δ cr2 value of ∼700. However, the actual value for the transition to intermittently turbulent regime seems to be affected by wall imperfections (Blondeaux & Vittori 1994;Vittori & Verzicco 1998), background disturbances (Ozdemir et al. 2014) and high frequency perturbations (Thomas et al. 2015). This possibly can explain the wide variability between the reported values in different experimental facilities (e.g. Merkli & Thomann 1975;Hino et al. 1976;Fishler & Brodkey 1991). Finally, Jensen et al. (1989) reported a value of Re δ cr3 = 3460. For this Reynolds number, Jensen et al. (1989) observed the existence of the logarithmic layer for 90 % of the period. The existence of the logarithmic layer is the criterion used by coastal engineers to determine the fully turbulent regime (e.g. Fredsøe 1984;Fredsøe & Deigaard 1992;Ozdemir et al. 2014). The analysis presented herein is not focused on the estimation of the critical Reynolds number for the intermittent turbulent regime. Instead, the presence of negative phase difference is examined in an effort to explain inconsistencies in the literature regarding the phase difference values. It is important to note that the negative phase difference occurs for Reynolds numbers higher than the threshold value Re δ cr2 whose estimation is challenging and seems to be affected by many parameters. Jensen et al. (1989) in their pioneering work presented results of bed-shear stress variations over time for a wide range of flows (Re δ between 257 and 3464). Of interest is the 'phase lead' diagram proposed by Jensen et al. (1989) which can be used to predict the phase difference between the instance when the maximum of the bed-shear stress happens with respect to the free-stream velocity. Such diagrams are commonly included in coastal engineering text books (e.g. Fredsøe & Deigaard 1992). According to this classic diagram, a bed-shear stress phase lead Δφ of π/4 can be predicted by the analytical solution for the laminar regime (for Re δ values up to ∼300), while for the fully turbulent regime the semi-empirical formula by Fredsøe (1984) can be applied (for Re δ larger than ∼1450). It is also accepted that, as Re δ increases, the phase lead will approach zero. This is expected to happen with a slow rate that scales with one over the logarithm of the Reynolds number ∼1/ log[Re δ ] (Spalart & Baldwin 1989). Although experimental and numerical results in the literature were contradictory (see Mier et al. (2021) for a review of these works), it is usually assumed that these two behaviours are smoothly connected, with the phase difference Δφ being reduced from π/4 to values of just below π/18 (10 • ) for Re δ = 1450 (see Fredsøe (1984), p. 1110, table 2). Mier et al. (2021) showed that this was not the actual behaviour, highlighting the presence of negative phase differences Δφ (phase lag) in a portion of the transitional regime. In addition, the variation of the mean flow structure over time was presented for a wide range of different Re δ . High turbulence levels during the deceleration phase were associated with the occurrence of the bed-shear stress phase lag with respect to the free-stream velocity maximum. However, due to the applied point-wise measurement technique in their analysis (laser Doppler velocimetry, LDV), it was really challenging to associate the presence of phase lag with the turbulent kinetic energy budget and also the presence of quasi-equilibrium conditions as presented in this work. Herein, the association of phase lag with the flow structure and the reason behind the presence of a phase lag are elucidated using direct numerical simulation (DNS) results. Hino et al. (1983) examined the flow structure of an oscillatory boundary layer at Re δ = 876. They presented turbulent statistics and showed that a logarithmic mean velocity profile with slope and intersect values close to those of equilibrium unidirectional boundary-layer flows exist for parts of the period. It is also worth noting that, after a close examination of the results by Hino et al. (1983) (figure 10 of their work), it can be seen 928 A33-3 that the phase difference for the case they presented should have a negative value. Jensen et al. (1989) studied the mean velocity and turbulent structure of the OBL for a much higher Re δ and identified a Re δ value of 3460 as the threshold condition Re δ cr3 for the fully turbulent OBL to occur, for which a logarithmic profile exists for almost all the parts of the period. However, similarly to many previous studies in the published literature, which are summarized later in the text as well as in Mier et al. (2021), both Hino et al. (1983) and Jensen et al. (1989) did not associate the turbulence characteristics with the presence of a phase lag in the wall-shear stress with respect to the free-stream velocity. Probably, this consistent neglect of the phase lag in the literature can be attributed partially to the difficulties associated with the measurement of wall-shear stresses in oscillatory flows (see § 3.1.1) but also due to the fact that the range of flow conditions studied by previous authors had either one or no Re δ cases that, according to the analysis by Mier et al. (2021), would exhibit the presence of a phase lag between wall-shear stress and free-stream velocity maxima.
Of interest to the analysis of the phase lag is the second burst of sediment entrainment which is commonly observed in oscillatory sheet flows (Ribberink et al. 2000(Ribberink et al. , 2008. This second sediment entrainment event during the deceleration phase has been puzzling the coastal engineering community. While the presence of enhanced Reynolds stresses during the deceleration phase has been recognized by coastal modellers (e.g. in the one dimension vertical 1DV analysis by Guizien, Dohmen-Janssen & Vittori 2003) as a characteristic of transitional OBL flows that enhances sediment entrainment, no detailed analysis is available in the literature that shows how the negative phase difference Δφ changes with Re δ . In fact, although the results from these Reynolds-averaged Navier-Stokes (RANS) simulations are promising, the predictions do not agree with the experimental observations (e.g. figure 7 in their work shows that the modified k−ω model predicts a secondary peak during the deceleration for Re δ = 565, Re w = 1.6 × 10 5 in their work's notation). This disagrees with the measurements by Jensen et al. (1989) and Mier et al. (2021). Mier et al. (2021) presented results of mean flow and turbulence statistics for a range of flows between 254 ≤ Re δ ≤ 1315. Their observations show that a logarithmic profile starts to exist in the middle of the deceleration phase for Re δ = 763 and as Re δ increases, the velocity profiles approach the log law for a longer part of the period and for a larger region of the boundary layer. These findings are in agreement with previous work by Hino et al. (1983, Akhavan et al. (1991a) and others, who have observed that a logarithmic profile starts appearing during the deceleration phase. However, the von Kármán constant κ and the intersect A s of those profiles do not seem to agree with those of the equilibrium logarithmic profile of the unidirectional zero-pressure gradient boundary-layer flows (ZPGBL) (Krug, Philip & Marusic 2017;Jiménez 2018). This is also consistent with some of the observations in recent DNS studies by Ebadi et al. (2019) for oscillatory channel flows and the large eddy simulation (LES) analysis by Kaptein et al. (2020) for oscillatory boundary layers, both of whom reported a range of slopes and intersects for the mean velocity profiles for oscillatory boundary layers that differ from those of equilibrium ZPGBL. We believe that this lack of equilibrium conditions is also relevant to the inconsistencies found in the literature regarding the presence of a phase lag between the bed-shear stress and the free-stream velocity. In addition, simplified models (e.g. Guizien et al. 2003;Blondeaux, Vittori & Porcile 2018), fail to accurately predict the underlying physics related to the turbulent flow-bed interaction and the presence of a phase lag between the bed-shear stress and the free-stream velocity in OBL, highlighting the need for further research of the bed-shear stress and turbulence characteristics in the transitional regime of OBL flows.
In the present analysis, oscillatory flows over smooth walls are considered. Oscillatory flows in real coastal applications occur over complex, rough, porous and moving beds (e.g. Pedocchi & Garcia 2009;Mazzuoli et al. 2020). More analysis is needed to examine the effects of the bed characteristics on the phase difference. However, the results for the analysis of the canonical/smooth-wall case may be relevant for more complex OBL flows.
The present work focuses on the theoretical and numerical analysis of turbulence characteristics of smooth-wall OBL of the same Re δ range as the one in Mier et al. (2021). Specifically, the available experimental data combined with DNS results are analysed to explain the mechanisms behind the presence of phase lags in transitional OBL flows and to elucidate the structure of unsteady boundary layers. The approximation of a 'quasi-equilibrium' state is discussed and the results are compared with other transitional flow data. Relaminarization and retransition effects are also discussed in the examined regime. Additionally, some turbulence characteristics are presented in comparison with the presence of phase lag in the bed-shear stress and with Re δ in the transitional regime.

Mathematical formulation and simulation details
DNS simulations were conducted to investigate the turbulence characteristics in OBL flow in the intermittently turbulent regime. The non-dimensional conservation equations governing the flow are where the non-dimensional quantities are denoted as tilded andũ = u/U o is the normalized instantaneous velocity, U o is the maximum free-stream velocity of the oscillation, δ is the Stokes' boundary thickness which is used to set dimensionless the spatial coordinates and derivatives (X = X /δ), P is the normalized pressureP = P/ρU 2 o andt is the normalized timet = tU o /δ. The last term of the right-hand side of (2.2) is the driving force, which was considered as a body force acting only in the streamwise direction e x . The governing equations were solved using the highly scalable spectral element method (SEM) based solver Nek5000 (Fischer, Lottes & Kerkemeier 2008). The P N − P N formulation was adopted for our simulations (Deville, Fischer & Mund 2002). In the SEM, the functions are represented as tensor-product polynomials of degree p within each element Ω e , e = 1, . . . , E. The domain Ω = ∪Ω e , where it is assumed that the elements do not overlap. The polynomial basis comprises Lagrange interpolating polynomials on Gauss-Lobatto-Legendre quadrature points, which ensures stability and allows efficient point-wise quadrature operations inside each e element (Deville et al. 2002). The incompressible Navier-Stokes algorithms implemented in Nek5000 ensure rapid (exponential) convergence in space and third-order accuracy in time. A de-aliasing procedure was used in our simulations following the 3/2 rule for the over-integration of the advection operation (Deville et al. 2002). The simulation results in the present analysis were conducted using seventh-order elements (p = 7) to maximize both spatial accuracy and computational efficiency for fast convergence of the simulations. All the linear terms were treated implicitly with pressure/velocity decoupling while the nonlinear advection terms were treated explicitly using third-order time integration and extrapolation schemes (third order backward differentiation and extrapolation schemes BDF3/EXT3). The simulations were conducted in a rectangular domain (figure 1). To simulate the development of the boundary layer under purely sinusoidal oscillatory forcing, periodic boundary conditions were imposed in the streamwise (x or x 1 ) and spanwise directions (z or x 3 ), while in the vertical non-homogeneous direction (y or x 2 ), a no-slip boundary condition was adopted at the bottom of the domain (y min ) and a stress-free boundary condition was imposed at the top of the domain (y max ).
Preliminary simulations were performed in a small domain (domain A, L × H × W = 160δ × 30δ × 40δ) which had the same width W and height H as in the works by Spalart & Baldwin (1989) and Ozdemir et al. (2014) for channel flows and double the length L of their domain. Most of the simulation results presented herein were conducted using a computational domain with L = 160δ length, H = 50δ height and W = 80δ width (domain B). This domain is significantly larger (double the size) in the streamwise and spanwise directions than the corresponding DNS channel simulations of similar flow regimes by Spalart & Baldwin (1989) and by Ozdemir et al. (2014). In the vertical direction the height was also larger than the corresponding half-channel height of the aforementioned studies. Finally, a simulation was conducted using a significantly longer domain (domain C, L × H × W = 250δ × 50δ × 80δ) and was used to validate that the results obtained were independent of the chosen domain size. In addition, a higher-order polynomial (eighth order) was adopted for this larger simulation to ensure that the results are also independent of the spatial resolution. This latter domain has similar dimensions to the domain used by Mazzuoli, Vittori & Blondeaux (2011) to study turbulent spot formation in OBL at similar Re δ . While no significant differences were noticed between the simulations of different domain size, it was observed that the statistics converge faster for the medium and large size domains compared with the smaller domain. This rate of convergence was evaluated after ignoring the first two periods, which were affected by the initial transient of each simulation. No significant differences were observed in the computation of friction factor and phase difference of bed-shear stress/free-stream velocity maxima, something that ensures the independence of the results from the size of the domain.
A summary of the examined cases as well as the grid resolution used for each run is presented in Here, L is the streamwise length, H is the height and W is the width of the domain; E x , E y and E z are the number of elements used in each direction, n x , n y and n z are the number of computational points/degrees of freedom in each direction (n i = E i p), n p is the total number of computational points for each simulation; dx + and dz + are resolutions in the streamwise and spanwise directions normalized using the maximum shear velocity u * max over the period and the kinematic viscosity ν. dy + min is the first point location in wall units; n T is the number of periods used to compute statistics. Letters A, B and C correspond to the computations using the small, medium and large domain sizes, respectively. kept uniform and is also included in table 1 in wall units (dx + = u * max dx/ν and dz + = u * max dz/ν, where u * max is the maximum ensemble-average shear velocity over the period and ν is the kinematic viscosity). In the vertical direction a Chebyshev distribution was adopted between y/δ = 0 and y/δ = 30 for all the domains. Uniform elements were adopted between y/δ = 30 and 50 for the medium and large domains (domains B and C). The present resolutions were chosen based on a series of tests, summarized in Appendix A, and were found to be sufficient for the target Reynolds numbers.
The initialization of the flow was conducted in a similar way as in the work by Ozdemir et al. (2014). Specifically, the laminar solution for ωt = π/2 was superimposed by an initial fluctuating velocity vector that followed a uniform probability distribution between −1 and 1 multiplied by 2.5 % and was scaled accordingly as a function of the vertical distance from the bed. This was used as the initial condition. A discussion about the need of special initialization for case 671B is included in paragraph 3.1. To avoid initialization effects, the first 2 periods of the flow were ignored for the computation of flow statistics. Thus, all the results presented herein were averaged over the remaining n T periods. The n T used from each case are also summarized in table 1. Analysis of the statistics showed that the results presented herein converged. Some characteristic plots that demonstrate the convergence of the statistics are shown in Appendix A. To ensure that enough statistics are collected, the symmetry between ωt and ωt + π was also exploited (see Appendix A).
All the simulations were conducted in the Petascale Computing Facility Blue Waters of the National Center for Supercomputing Applications of the University of Illinois at Urbana-Champaign over three years (2017)(2018)(2019)(2020). A typical simulation was performed using 512 computational nodes and required on average approximately 6000 node hours per period for simulations conducted using domain A, 12 000 node hours per period for simulations conducted using domain B and 50 000 node hours per period for simulations conducted using domain C. The total computational resource required for the present study is approximately 1 300 000 node hours.
The volume integrals of the turbulent kinetic energy (TKE), the TKE production and TKE dissipation rates were monitored at the beginning of each simulation to ensure that the initial conditions and the initial transient until the flow reaches equilibrium state did not affect the flow statistics. This same approach was used by Ozdemir et al. (2014). Similarly, to their observation, the flows considered in our study reached equilibrium after a couple of periods (in most of the cases a single period was enough). Thus, the two periods were discarded from each simulation before statistics were collected. In Vittori & Verzicco (1998) a different initialization was used based on analytical fields from Vittori (1992). They also introduced wall imperfections as a mechanism for triggering transition to turbulence in the Stokes layer similarly to the works by Blondeaux & Vittori (1994), Verzicco & Vittori (1996) and Costamagna, Vittori & Blondeaux (2003). In their analysis, convergence was obtained at the end of the first cycle for Re δ < 550 and Re δ > 600. For these flows, similarly to the present work, Vittori & Verzicco (1998) discarded the first two periods before collecting statistics. It is worth noting that Vittori & Verzicco (1998) observed that a higher number of periods was required to reach equilibrium at a Reynolds number close to 550-600. Later, Ozdemir et al. (2014) reported that the transition to turbulence at a Re δ close to 600 may depend on the levels of 'background disturbances'. In our analysis, the case of Re δ = 679 also showed signs of dependence on initialization regarding whether self-sustained turbulence would be achieved. As discussed later, similarly to the work by Ozdemir et al. (2014), self-sustained turbulence was maintained when the flow was initialized using a fully developed turbulent field. This is also reasonably close to the predictions from the stability analysis by Blennerhassett & Bassom (2002), who propose a critical Reynolds number of ∼709 for the transition. For the case of Re δ = 679 considered herein, the flow reached equilibrium after two periods.

Results and discussion
While all the simulations performed herein were initialized at the instance when the maximum free-stream velocity occurs, for the presentation of results, the typical convention used in the coastal engineering literature is adopted. According to this convention the results are presented in intervals within ωt = (0, π) for which the free-stream velocity is zero at ωt = 0, ωt = π/2 corresponds to the maximum free-stream velocity, the free stream velocity changes direction at ωt = π etc.

Mean velocity profiles and root-mean-square velocity fluctuations -Reynolds number effect
The presence of the universal logarithmic velocity profile has been associated with the equilibrium flow conditions when the mean/ensemble quantities become constant and independent of the streamwise location. The average velocity profile obeys the well-known log-law equation for a smooth boundary is assumed to be valid close to the wall (typically ∼20 % of the water depth, Nezu & Nakagawa 1993), while far from the wall the wake effects can become significant (Krug et al. 2017;Jiménez 2018). Equilibrium boundary layer (BL) flows defined e.g. in Clauser (1954), Rotta (1962) and Townsend (1980) for zero-pressure gradient flows and extended for boundary layers under pressure gradients (Clauser 1954;Townsend 1956;Coles 1957;Mellor & Gibson 1966) are flows for which the proportional contribution of each term in the flow equations remains constant with respect to the streamwise direction. The theoretical analyses by Townsend (1956) and Mellor & Gibson (1966) showed that, to reach a near-equilibrium state in spatially accelerating flows, a power-law relation ( is required for the free-stream velocity variation. The analyses of these flows have become the focus of many theoretical, numerical and experimental works (reviews can be found in Bobke et al. 2017;Kitsios et al. 2017;Vila et al. 2020). Recently, there have been efforts to collapse the velocity profiles of self-similar adverse pressure gradient BL flows (Kitsios et al. 2016;Bobke et al. 2017). Similar efforts to identify the proper velocity and length scales as well as the conditions that will collapse the velocity profiles of non-equilibrium BL flows, when those reach a near-equilibrium state, can be found in the spatially varying boundary-layer literature Castillo & Wang 2004, and more).
Here, the canonical purely sinusoidal OBL flows were examined, which belong to the general family of unsteady and non-equilibrium flows. However, even from the early works of Hino et al. (1976) and Jensen et al. (1989), the presence of a self-similar logarithmic profile had been observed. The present analysis focuses on describing the behaviour and identifying the characteristics of the BL properties and how these are associated with the mean velocity profile slope and intersect as well as the TKE budgets. There exist parts of the period that reach a 'near-equilibrium' state which will be loosely called 'quasi-equilibrium' here. For such stages, the flow characteristic and BL properties will be analysed in the following paragraphs, in an effort to identify the necessary and sufficient conditions that will result in the log-law profiles.
Clauser's parameter, β = (δ * /τ b )(dP/dx), where δ * is the displacement thickness and dP/dx is the mean streamwise pressure gradient, has been used in the past to examine spatially accelerating BL. This parameter is negative β < 0 for accelerating flows, positive β > 0 for decelerating flows, β = 0 for zero-pressure gradient flows and β → ∞ for flows that experience separation. In the context of temporally accelerating flows, Clauser's parameter can be written as β = (−ρδ * /τ b )(dU/dt). Thus, the equivalent variation of β in space (as a function of streamwise direction x, β(x)) can be easily expressed with respect to time/phase (β(ωt)). The structure of the boundary layer experiencing positive values of β will be significantly different from the one in a BL under negative β values (Perry, Marusic & Jones 2002). In addition, experimental and numerical results from the OBL literature had shown that flow relaminarizes for part of the acceleration phase.
Relaminarization is also known to occur when a severely strong favourable pressure gradient interacts with an initially turbulent boundary layer, causing the flow to turn fully laminar. The precursor to relaminarization, 'laminarescence' is the stage when the flow loses part of its turbulent behaviour without becoming fully laminar. Herein, the terms 'laminarization' and 'laminarization phase' are loosely used to describe the portion of the acceleration phase where the flow loses energy. Interested readers can refer to the reviews by Narasimha & Sreenivasan (1979) and Sreenivasan (1982), on the relaminarization of unidirectional flows. In oscillatory flows of the intermittently turbulent regime (and for Re δ larger than ∼600), as shown by Mier et al. (2021) and the previous works by Merkli & Thomann (1975), Hino & Sawamoto (1975), Akhavan et al. (1991a), Akhavan et al. (1991b), Ozdemir et al. (2014) and Carstensen et al. (2010), the turbulence statistics increase during the deceleration phase.
The ensemble-averaged velocity profiles were calculated after averaging the instantaneous fields with respect to the phase of the period and the homogeneous streamwise and spanwise directions. The ensemble-average value of a quantity φ is denoted asφ and is computed asφ( y, where S A is the horizontal homogeneous direction area and n T is the number of periods over which the ensemble averaging is conducted. The local fluctuation values are defined as φ = φ −φ. Due to the symmetry between the positive and negative parts of the oscillation, the results presented are only for half the period (ωt = 0 − π). No significant differences are observed between the positive and negative parts of the oscillation. Ozdemir et al. (2014) observed that, for a Re δ of 600, initializing the flow with the results of a higher Re δ (Re δ = 1000) can lead to a self-sustained transitional flow. When a similar initialization was attempted for the case 552A (initialized using the fluctuation results of case 1036A), the turbulence was not sustained and the initially turbulent flow lost its energy after a period. This was not the case herein when initializing the flow for the Re δ = 671. For this flow, similarly to the observations by Ozdemir et al. (2014), initializing the flow with the fluctuation field from case 1036B leads to a self-sustained transitional flow. When the standard initialization described earlier was used, the flow was turning laminar after a period. The results presented here for Re δ = 671 are those that correspond to the self-sustained turbulence case. For the other cases, no sensitivity to the initial conditions was observed, as self-sustained turbulence could be maintained and the ensemble-averaged profiles were independent of the initial conditions.
The mean flow/ensemble-averaged velocity profiles for all the examined flows of table 1 are shown in figure 2. The numerical results are plotted using wall units in figures 2(a), 2(b) and 2(c), for which streamwise velocity was normalized using the corresponding 1315  π 2 3π 4 10 0 10 1 10 2 10 3 10 0 10 1 10 2 10 3 10 0 10 1 10 2 10 3 10 0 10 1 10 2 10 3 10 0 10 1 10 2 10 3 10 0 10 1 10 2 10 3 10 0 10 1 10 2 10 3 10 0 10 1 10 2 10 3 y + y + y + 0 1 2 0 0.5 ensemble-average shear velocity u * , and are plotted vs y + coordinates (where y + = u * y/ν) for ωt = π/4, π/2 and 3π/4. The ensemble-average velocity defect profile normalized with free-stream velocity U ∞ and Stokes length δ are plotted in figure 2(d) for ωt = π/4, (e) for ωt = π/2 and ( f ) for ωt = 3π/4. Finally, the ensemble-average velocity defect profiles normalized with shear velocity u * vs y + are shown in figure 2(h) for ωt = π/4, (i) for ωt = π/2 and (g) for ωt = 3π/4. It is shown in these plots that, as Re δ increases, the velocity profiles for ωt = π/2 and 3π/4 approach a logarithmic behaviour. However, as is shown by the orange dashed lines, which correspond to a logarithmic fit of the data, these logarithmic profiles are not always self-similar and thus the slope and intersect values do not have physical meaning. No clear logarithmic profiles are observed for ωt = π/4, which is in agreement with the experimental observations of Mier et al. (2021), who observed that flow looses part of its energy and behaves similarly to a relaminarized flow during the parts of the period which experience severe favourable pressure gradient; namely, during the relaminarization phase. Hino et al. (1983) argued that OBL flows can show a behaviour similar to the fully turbulent state for Re δ > 800. Sarpkaya (1993) and Jensen et al. (1989) have shown that OBL can exhibit behaviour that mimics the fully turbulent regime for Re δ ≥ 1750, which was similar to the LES observations by Salon, Armenio & Crise (2007). Figure 3 shows results of the normal components of the Reynolds Stresses tensor in wall units for the examined flows. The streamwise, spanwise and vertical profiles of the ensemble-average u 2 + , w 2 + and v 2 + are plotted for ωt = π/4, π/2 and 3π/4, respectively. For the streamwise component, it is shown in figures 3(b) and 3(c) that there is a trend towards an overlapping profile for Re δ 1036 and 1315. For the examined Re δ range, the ensemble spanwise (figure 3e, f ) and vertical (figure 3h,i) components do not seem to converge to a profile for ω = tπ/2 and 3π/4, although there is indeed a trend towards the high Re δ values. Figures 3(a), 3(d) and 3(g) show that the Reynolds stresses grow faster for higher Re δ , which is a sign that, as expected, more turbulent flows will start developing earlier during the period. Still, no significant overlap towards similar profiles are observed for ωt = π/4.

Friction coefficient f w
The friction coefficient f w , defined as f w = 2τ max /ρU 2 o is a measure of the maximum ensemble-average bed-shear stress τ max over the period. The estimation of the values of f w and the development of graphs for its prediction for various wave conditions have been the focus of the early works by Kajiura (1964), Yalin & Russell (1966), Jonsson (1966), Riedel, Kamphuis & Brebner (1973) and Kamphuis (1975). These studies highlighted the importance of the flow regime on the friction coefficient (Kajiura 1964;Jensen et al. 1989;Sarpkaya 1993). Figure 4 summarizes the data from the literature combined with the results of the present analysis. Note that the f w data are plotted both with respect to Re δ and Re w , which is a different Reynolds number, defined using half of the oscillation amplitude (note the relationship between the two Reynolds numbers, Re w = Re 2 δ /2). The analytical solution and the theoretical solution by Fredsøe (1984) are also plotted. In figure 4, a plot using a linear axis is also included for the range between Re δ = 200 and Re δ = 2000. This plot shows the spread of the reported f w values in the transitional regime. Differences as high as 100 % are observed in some cases. These discrepancies may be relevant to the difficulties that experimentalists face when trying to measure the bed-shear stress in unsteady flows as well as the crisis that seem to take place near the critical Reynolds number. The DNS results of the present study are in good agreement with the previous experimental results by Hino et al. (1983), Jensen et al. (1989), Sarpkaya (1993) and Carstensen et al. (2010) and the DNS results of Spalart & Baldwin (1989), as well as the experimental observations of Mier et al. (2021). A detailed comparison between the DNS results and the experimental observations is included in table 2 of Appendix B. Data from Kamphuis (1975) underestimate the friction by a factor of 20%. In general, for the data in the literature a deviation from the laminar prediction starts being observed for Re δ > 600. For values Re δ ≥ 1036 the present data are in good agreement with the observations of Sarpkaya (1993) and Carstensen et al. (2010). The theoretical prediction by Fredsøe (1984) seems to predict well the experimental results for Re δ ≥ 3000, which is reasonably close to the Re δ cr3 value of 3460 proposed by Jensen et al. (1989) as the threshold value for the fully turbulent regime. 10 -1 f w Spalart and Baldwin (1989) Laminar solution Fredsøe's (1984)  The variation of shear stress over the period can be seen in figure 5. For Re δ = 552, it is shown that the bed-shear stress of the numerical results agree well with the prediction of the analytical solution of the laminar regime, τ b /ρU 2 o = √ 2/Re δ sin(ωt + π/4) (Batchelor 1967). Indeed, a π/4 radians phase lead Δφ of the bed-shear stress maximum with respect to the free-stream velocity has been observed. For Re δ = 671, two peaks can be observed, one during the acceleration phase, which can be linked to the laminar flow regime, and one weaker peak during the deceleration phase. The first peak happens closer to ωt = π/2 compared with the Δφ value of Re δ = 552 (Δφ = 37 • or 0.65 radians). Mier et al. (2021) associated the second peak with the transition to turbulence. As Re δ increases further, this second peak increases in magnitude and eventually the first, 'laminar' peak, vanishes due to the enhanced effect of the 'turbulent' peak for Re δ ≥ 1036. It is also important to note here that, for Re δ = 763, the bed-shear stress maximum happens during the deceleration phase, i.e. 'lagged' with respect to the free-stream velocity maximum. A similar behaviour has been observed for Re δ = 819 although the peak of the bed-shear stress seems to take place earlier during the deceleration phase compared with Re δ = 763. This is consistent with the experimental observations by Mier et al. (2021). After a close inspection of the ensemble-averaged bed-shear stress measurement by Hino et al. (1976) it is easy to conclude that they had also observed negative phase differences for a flow of Re δ = 876 (see p. 373, figure 10 in their work). Similarly, observations have been made in the literature both experimentally Carstensen et al. 2010) and numerically (Spalart & Baldwin 1989;Ozdemir et al. 2014;Bettencourt & Dias 2018;Ebadi et al. 2019). A more extended presentation and summary of these works can be found in Mier et al. (2021). However, all these studies have not stressed the need to revise the phase difference diagram that is included in the classic papers by Jensen et al. (1989), Sarpkaya (1993) and Carstensen et al. (2010), among others.
A revised version of the classic phase lead diagram (here plotted as phase difference) is plotted in figure 6. In addition to the data of the present study and the experimental measurements by Mier et al. (2021), the measurements done by Hino et al. (1976), Jensen et al. (1989), Sarpkaya (1993) and Carstensen et al. (2010) and the DNS results by Spalart & Baldwin (1989), Ozdemir et al. (2014) and Bettencourt & Dias (2018) are also plotted for comparison. Note that the phase difference Δφ values that are plotted here (e.g. for the case of Hino et al. (1983) and Jensen et al. (1989)) are those as read after processing the bed-shear stress signals rather than the positive values reported earlier. From figure 6 it can be observed that there is a clear threshold value Re δ ∼750 after which the data of the current study and the literature show negative phase difference. Small differences and discrepancies regarding the absolute value of the Δφ can be attributed to the challenges associated with the measurement of bed-shear stress under oscillatory conditions in general and the different temporal resolutions at which the bed-shear stresses were measured. Some dependencies of the numerical results on the height of the domain, which have been reported by Kaptein et al. (2019), may also explain some of the variations of the numerical results. Kaptein et al. (2019) used LES to examined the effect of h/δ ratio (where h is the height of the domain). They concluded that, for h/δ ≥ 40, the velocity, the turbulence characteristics and bed-shear stress results converge to those for h/δ → ∞, which is the focus of the present study. The data presented here show no significant effect on the height of the domain. The final results, however, were extracted in a domain with h/δ = 50, which is above the height-limited conditions reported by Kaptein et al. (2019). Finally, the consistent presence of negative phase differences in the literature shows in a convincing way that phase lag does exist in the intermittently turbulent regime for 750 ≤ Re δ ≤ 1000.
The DNS results and the Δφ measurements by Mier et al. (2021) are also summarized in table 2 (Appendix B). To ensure that the bed-shear stress signal will be captured in maximum detail, the bed-shear stress was computed while running the simulations at every time step, which had a significantly smaller size than 1 • and was determined in a way that the Courant-Friedrichs-Lewy number was less than 0.45 for all portions of the cycle. The results presented here are the ensemble-average results of these extracted values. Note that the symmetry between ωt and ωt + π was also used to enhance the collected statistics.

Flow structures
In figure 7, time variation of some characteristic dimensionless Reynolds numbers are reported, which will allow for the comparison of the examined OBL flows with previous studies from the literature. In figure 7 a, the shear Reynolds number Re δ * = u * δ/ν as introduced by Sarpkaya (1993) is presented. In figures 7(b) and 7(c), the classic Reynolds numbers that are usually used in studies of developing boundary layers are plotted, namely Re θ = U ∞ θ/ν, defined using the momentum thickness θ and the mean free-stream velocity U ∞ (with U ∞ (ωt) = U o sin(ωt)), and Re τ = u * y max /ν defined using the ensemble-average u * (which is also varying with ωt) and the thickness of the boundary layer y max defined in agreement with the conventions by Jensen et al. (1989) and Mier et al. (2021) as the height for which the shear stress becomes zeroτ = 0. The displacement thickness and the momentum thickness are computed as Note that, for the case of oscillatory flows, the definitions of (3.2) and (3.3) can lead to negative δ * and θ values during the acceleration phase. This is the result of the characteristic near-bed overshoots of the velocity that can be observed in oscillatory flows. The parts of the flow when overshoots occur are excluded from the present analysis since there is no clear physical meaning in negative momentum and displacement thickness values.
From figure 7 it can be observed that Re θ and Re τ continue to grow during the deceleration phase, despite the fact that both U ∞ and the shear velocity may decrease. This is explained later in figure 12(b), which shows that both y max and θ continue growing for a portion of the deceleration phase.
In figure 8, the coherent structures visualized using the λ 2 criterion introduced by Jeong & Hussain (1995) are shown for Re δ = 763 and at ωt = π/2, 7π/12 and 2π/3. The isosurfaces of λ 2 = −0.0025 are plotted, coloured using the velocity magnitude u/U o . In figure 8(a), the generation of hairpin packets that are randomly distributed in space is shown. These packets, generated towards the end of the acceleration phase, later grow and form what looks like the TS (figure 8b) described by Wu et al. (2017), until these TS merge later, during the deceleration phase (figure 8c). The fact that TS keep growing and merging for as long as 45 • (π/4 radians) during the deceleration phase is another sign that the OBL keeps developing during parts of the deceleration phase until near-bed flow reversal occurs. The final state of what looks like the final developed stage of the flow is shown in figure 9 for ωt = 3π/2, during which hairpin vortices have occupied most of the simulation domain. It is also instructional to notice the variation of bed-shear stress during this period (also plotted in figures 8 and 9). Indeed, as the flow structures start developing, the bed-shear stress continues to increase during the deceleration phase until ωt = 2π/3. Then, the bed-shear stress starts decreasing again under the effect of deceleration.
This late transition of the flow to an increased energy state is the reason for the phase lag in the bed-shear stress signal that is shown in figures 5 and 6.
In agreement with the work of Wu et al. (2017), the TS are generated randomly in space and time. However, it is clear that, as Re δ increases, the TS are generated earlier during the period (Carstensen et al. 2010). These OBL flow features have been examined in the past by Carstensen et al. (2010), who observed their presence for as low as Re δ = 550 (Re w = 1.5 × 10 5 in the original work). In the present analysis, TS started appearing for Re δ = 671, since Re δ = 552 behaves in a laminar fashion. Sarpkaya (1993) had identified a critical Re δ * value between 20 and 28 for which the high-low velocity streaks appear and disappear periodically. This threshold value of (Re δ * ∼24) is plotted with an orange dotted line in figure 7(a). This threshold value is in good qualitative agreement with our observations, which show Re δ * > 24 for part of the period of the intermittently turbulent flows examined here. Park et al. (2012) also observed that TS may appear as 'individual' for Re θ ∼300 and as 'merged' for Re θ ∼500 (plotted in figure 7(b) using blue and orange dashed lines, respectively). For the flow presented here (Re δ = 763), Re θ reaches a value of 436 for the state which is shown in figure 9. This value, although close enough to the Re θ ∼500 by Park et al. (2012), represents a state when the flow characteristics mimic those of fully turbulent flow. As shown in figure 10, this is the part of the flow when a velocity profile that resembles the characteristics of the logarithmic profile can be observed. Later, we will show that the actual slope and intersect of this velocity profile differ from those of the equilibrium BL. Figure 10 also includes the laminar solution and the logarithmic profile by VanDriest (1956) where κ = 0.41 and A v = 26. This profile agrees well with (3.1) for y + ≥ 30. A logarithmic fit to the data is also plotted with orange dashed lines. The LDV data by Mier et al. (2021) for the same Re δ , data by Hino et al. (1976) for Re δ = 876 and Spalart & Baldwin (1989) and Jensen et al. (1989) for Re δ = 1000 are plotted for comparison. It is shown that a logarithmic profile is observed for ωt = 2π/3 and 3π/4. These are the instances during the period when τ b and Re θ (also plotted in figure 10) are also increased. In addition, the shape factor H, defined as the ratio between displacement and momentum thickness (H = δ * /θ) approaches a constant value. This value is slightly larger than the H∼1.4 that was observed by Mier et al. (2021) as the converged value at the fully turbulent stage. As shown later, this happens due to an incomplete transition, in a narrow sense. Figure 10 shows a general good agreement between the DNS results and the experimental observations, with the exception of results for ωt = 5π/6, which is an instant really close to the near-bed flow reversal when τ b = 0. Thus, these deviations which are magnified due to the low u * values are attributed to the challenges associated with the measurement of τ b values that are close to zero. Also, the comparison with the results by Jensen et al. (1989) and Spalart & Baldwin (1989) shows that a logarithmic profile occurs earlier for higher Re δ values. Finally, the figure shows some deviations between the results presented here and the observations by Hino et al. (1976). These deviations had been reported earlier by Jensen (1989), who also observed that his velocity profiles differ from those by Hino et al. (1983) and highlighted the challenges of experimental works to estimate near-bed velocity profiles in OBL flows. Figure 11 show the corresponding results for Re δ = 1315. Again, an excellent agreement between the numerical and experimental results is observed (with the same deviations at ωt∼5π/6 when τ b → 0). As Re δ increases, the part of the circle for which a logarithmic  (Bobke et al. 2017). Mellor & Gibson (1966) showed that G values become constant in near-equilibrium conditions for flows with constant Clauser parameter β. Here, despite the fact that β values are not constant, we will evaluate the values of H and G that correspond to the part of the flows for which a self-similar velocity profile is observed. Figure 12 summarizes the results of the examined flows for the Re δ in the intermittently turbulent regime. The temporal variations of the skin-friction coefficient c f , normalized boundary-layer thickness y max /δ, displacement thickness δ * /δ and momentum thickness used across the period solely as diagnostic parameters for defining the part of the period that is 'near equilibrium' (when κ and A s get the standard values). This same approach has been used in the previous works by Hino et al. (1983) Figure 12 shows the path towards a transition to a state that mimics well the characteristics of the fully turbulent regime. This is illustrated for the examined flow range of Re δ 671 and 1315. Here, Re δ = 671 experiences a late transition that never actually gets to a stage that qualitatively looks similar to turbulent. The slope and intersect of the velocity profile deviate significantly from these of equilibrium for the whole period. Also, diagnostic parameters H and G do not get close to the threshold values of is the threshold value for phase lag to exist (figure 6), the quantities start approaching the behaviour of the fully turbulent regime during the deceleration phase; a logarithmic velocity profile is observed (see also figure 10), H and G values get close to the equilibrium values and κ and A s values advance towards the equilibrium values. In addition, the characteristic increase of the skin-friction coefficient c f takes place during the deceleration phase. As shown later, the absolute values of κ and A s do not actually reach the 0.41 and 5.1 values but they start getting close. A similar behaviour is observed for Re δ = 819. The transition starts earlier, but again, κ and A s values do not converge to the equilibrium values, although the near-equilibrium conditions allow the formation of a logarithmic profile during the deceleration phase. For Re δ = 1036, the behaviour starts changing, with the diagnostic parameters H and G reaching the equilibrium values earlier and the equilibrium κ and A s values start appearing for parts of the circle. This becomes even more convincing for Re δ = 1315. However, as the flow experiences a severe adverse pressure gradient (i.e. β ≥ 1) the flow starts deviating from these equilibrium conditions. This finding expands slightly the Re δ ≥ 1750 limit for OBL behaviour that mimics the fully turbulent regime observed by Sarpkaya (1993) and Jensen et al. (1989) and Salon et al. (2007). It seems that it is possible to observe near-equilibrium conditions that match the universal velocity slope and intersect (κ = 0.41 and A s = 5.1) for as low as Re δ ∼1000.

Turbulence characteristics, quasi-equilibrium and logarithmic profile
The TKE budgets were computed for the examined flows. The process adopted here is similar to that of Pedocchi et al. (2011), who conducted DNS for Re δ = 1679. The TKE budget can be written as where k is the TKE, Pr is the production of TKE, is the TKE dissipation rate, Π is the velocity pressure gradient term, Π s is the pressure-strain term, D is the viscous diffusion of TKE term and T is the turbulent transport of TKE term. All the terms are computed as presented in Pedocchi et al. (2011) and Vinuesa et al. (2017). For the analysis of the TKE budget both the profiles of each term and the integrals of production and dissipation inside the boundary layers are used. These integrals are defined as In figures 13 and 14 the TKE budgets over the period for Re δ 763 (case 763B) and 1315 (case 1315B) are plotted. All the terms are normalized using δ/U 3 o . In addition, the ratio of the integral of the production to the integral of dissipation I Pr /I inside the boundary layer is plotted using blue dashed lines. Also, the normalized integral of the TKE inside the boundary-layer thickness y max (I k = y max 0 k dy) is plotted divided by the maximum value |I k | max (red dashed line). The normalized bed-shear stress τ b /τ b max is plotted in green. Figure 13 shows that TKE terms start increasing during the acceleration phase. Production integral I Pr becomes larger than the absolute magnitude of the dissipation rate integral |I | for ωt ≥ π/6 and continues increasing until ωt∼π/2. After that, the I Pr /I ratio drops and approaches a value close to ∼1. This part of the period when I Pr /I ∼1 is the same part in which a logarithmic profile exists. Also, it is the same instance when the bed-shear stress and the integral of the k over the boundary-layer thickness reach their maximum values. Later, when y max decreases, I Pr /I and I k /|I k | max also decrease.
As Re δ increases, turbulent statistics increase earlier during the acceleration phase. In figure 14 it is shown that I Pr /I ratio becomes larger than 1 at ωt ∼π/12. The I k integrals continue to get their maximum values during the deceleration phase as y max continues to increase, until eventually becoming nearly zero at the near-bed reversal. Finally, I Pr /I approached a value of 1 at the parts of the circle when the logarithmic velocity profile exists.
It is also important to note that for the parts of the period when I Pr /I is nearly one, the TKE budget terms remain similarly constant when they are normalized in wall units using ν/u 4 * . Indeed, when the results are plotted with the shear-velocity/viscous normalization ( figure 15 shows the results for Re δ = 1315), this similarity becomes obvious. Pedocchi et al. (2011) performed a similar analysis for Re δ = 1679 (originally reported using Re w = 1.41 × 10 6 ). A peak production value at y + ∼10 was observed in their DNS results. A similar peak can be seen in figure 15. Further, close to the wall, where the dissipation rate ν/u 4 * is balanced by the viscous diffusion Dν/u 4 * and both obtain a value of ∼0.27, Π and Π s contributions are in general smaller across most of the period. This can be noticed in both figures 14 and 15. Some of their contributions may seem to become important when expressed in wall units at ωt∼11π/12, which should be expected due to the enhanced velocity-pressure gradient term Π . However, these are normalized with a Re θ = 111 ωt 0 π/2 π ωt 0 π/2 π ωt 0 π/2 π ωt 0 π/2 π ωt 0 π/2 π ωt 0 π/2 π ωt 0 π/2 π ωt 0 π/2 π ωt 0 π/2 π ωt 0 π/2 π ωt 0 π/2 π ωt I Pr /|I | Figure 13. TKE budget for Re δ = 763. Terms are normalized using δ/U 3 o . Also plotted are the time variation of the ratio between the integral of production and the integral of the dissipation rate I Pr /I , the normalized bed-shear stress τ b /τ bmax and the normalized TKE integral I k /I kmax .
shear velocity that is very close to zero and are still small when the absolute magnitude of these contributions is concerned. Finally, the variation of turbulent transport contributions remains small far from the wall (y + > 30), which is the region where logarithmic profile is in general observed. This is a region that all the terms except Pr and become practically zero. Close to the wall, the turbulent transport term switches sign between negative at y + ∼10 to positive for y + ∼4 until it becomes zero at the wall.
From the TKE budget analysis above it becomes obvious that TKE budget is dominated by the production and dissipation terms far from the wall (y + > 30). The presence of a logarithmic profile is usually associated with local equilibrium between the production and dissipation rate of TKE (Pr ≈ ). In figure 16 the ratio between TKE production Pr and the absolute value of dissipation rate | | is plotted for Re δ between 671 and 1315 for ωt 0 π/2 π ωt 0 π/2 π ωt 0 π/2 π ωt 0 π/2 π ωt 0 π/2 π ωt 0 π/2 π ωt 0 π/2 π ωt 0 π/2 π ωt 0 π/2 π ωt 0 π/2 π ωt 0 π/2 π ωt I Pr /|I | Figure 14. TKE budget for Re δ = 1315. Terms are normalized using δ/U 3 o . Also plotted are the time variation of the ratio between the integral of production and the integral of the dissipation rate I Pr /I , the normalized bed-shear stress τ b /τ bmax and the normalized TKE integral I k /I kmax .
various ωt values. It can be observed that the ratio Pr/| | does approach 1.0 (marked with an orange dashed line in figure 16) in regions where a logarithmic profile exists. This is the near-equilibrium behaviour described above. As Re δ increases, the region and part of the period for which local quasi-equilibrium exists increases. This is in agreement with the experimental observations of Mier et al. (2021), who observed that as Re δ increases the logarithmic profile of the mean velocity profile becomes valid in a more extended region and for a more extended part of the period.
(2017) simulated a developing boundary layer for 80 < Re θ < 3000 and there are a lot of similarities between the structure of the flow (TS and developing BL) in the examined range and their data. In their work, they evaluated the accuracy of their results and the establishment of equilibrium conditions by comparing the normalized dissipation rate + = /(u 4 * /ν) with the following estimation for the equilibrium conditions by Tennekes & Lumley (1972) and Balint, Wallace & Vukoslavčević (1991): where κ = 0.41. As Re θ values increase, their results converge to the prediction of (3.9). For example their results for Re θ = 150, 199, 334, 670 and 1000 are plotted in figure 17.  (b) (a) Figure 16. TKE production to absolute dissipation rate ratio Pr/| | for the examined flows.
The boundary layer grows and the near-wall dissipation rate values start growing until eventually the dissipation rate profiles converge to the prediction for Re θ = 670 after the characteristic overshoot at Re θ = 334. Figure 17 shows the results of the dissipation rate profile at the time instance when the profile is closest to the equilibrium profiles for Re δ 819. Figure 17 shows that this final stage is still far from the equilibrium results by Wu et al. (2017), but still, as the DNS results for + get closer to the prediction by (3.9), bed-shear stress starts exhibiting the characteristic secondary peak during the deceleration phase, which is in agreement with the observed behaviour by Mier et al. (2021). The Re θ value at the moment when second peak exists is 347. A movie in the supplementary material available at https://doi.org/10. 1017/jfm.2021.827 shows the variation of dissipation rate profiles over the period for Re δ 819 (movie 1). Although Re θ continues growing, the dissipation rate profile departs from the equilibrium profiles under the effect of the pressure gradient. The case in figure 17 corresponds to the bed-shear stress phase lag case. This late transition that takes place during the deceleration phase, is not actually complete in the sense that neither the slope κ and the intersect A s velocity logarithmic profiles reach the equilibrium values (0.41 and 5.1 respectively), nor does the dissipation rate profiles reach those of the equilibrium conditions. This incomplete and late transition is also observed for Re δ = 763, which is the threshold Re δ value for a phase lag to exist. However, the + profiles and the logarithmic velocity slope and intersect do get closer to the asymptotic values presented here as the Reynolds number increases.
At higher Re δ , the + profiles start to match those of the equilibrium condition by Wu et al. (2017) and the predictions by (3.9). Figures 18 and 19 show the profiles in the range of the period that matches the equilibrium data (the variation of the profiles over time can be found in movies 2 and 3). At the same time, the slope and intersect of the logarithmic  profile start to match their equilibrium values (figure 12) and the shape factor H and G values become ∼1.4 and ∼7.4, respectively. For Re δ = 1036 this happens for the part of the period between ωt = 85π/180 and 11π/18. This range increases to ωt = 7π/18 and 11π/18 for Re δ = 1315, which is also consistent with the diagnostic parameters of figures 12 and 16. In other words, a logarithmic velocity profile may be observed, when the ratio of Pr/| | becomes ∼1, and other diagnostic parameters examined here, such as the shape factors H and G get values similar to the fully turbulent regime values. However, the slope and intersect of these logarithms will not be those of the equilibrium profile until Re δ ∼1000. As shown in the dissipation profiles of figure for 17, this is explained due the fact that the transition remains incomplete, in the narrow sense of dissipation profiles that do not match those of the equilibrium conditions and could not be predicted by (3.9). Flows in this range (600 < Re δ < 1000) do not develop quickly enough to complete the transition before the flow starts being affected by the deceleration/adverse pressure gradient. When Re δ is larger than approximately ∼1000, the equilibrium slope and intersect are achieved and OBL flows start mimicking the behaviour of the fully developed ZPGBL after the completion of the retransition process. For these flows transition starts earlier during the period which allow its completion to happen.

Conclusions
DNSs were conducted for oscillatory flow conditions over a flat smooth wall in the transitionally turbulent regime. The range of the examined Re δ values was the same as in the experimental study by Mier et al. (2021). In such a range of Re δ , experimental and numerical observations in the literature show inconsistencies regarding the phase difference Δφ between the instance of the period when the maximum in bed-shear stress occurs with respect to the instance when the free-stream velocity maximum occurs. Analyses of the mean flow structure and TKE budgets were also conducted and the near-equilibrium state of an unsteady OBL flow was defined. The key results of the numerical and theoretical analyses presented here are summarized below: (i) The present analysis confirms the large-scale experimental observations by Mier et al. (2021) who observed that a phase lag between the bed-shear stress and free-stream velocity maxima exists for intermittently turbulent OBL flows within the range of 763 < Re δ < 1000; thus, a modification of the widely used 'phase lead' diagram found in the literature (e.g. Jensen et al. 1989) and in many coastal engineering handbooks (e.g. Fredsøe & Deigaard 1992) is required. Figure 6 summarizes data of the present work together with experimental and DNS results found in the literature and can be used to estimate the phase difference Δφ between the bed-shear stress and the free-stream velocity maxima. The typical phase lead of 0.79 radians (45 • ) is observed until Re δ ∼550. After that, the phase lead decreases until it reaches a value of ∼0.61 rads (∼35 • ) for Re δ < 750. For Re δ ∼750, a phase lag (Δφ < 0) of 0.46 rads (26.5 • ) occurs. For higher Re δ the phase lag becomes smaller until it becomes zero at approximately Re δ = 1000; Δφ has positive values up to approximately ∼π/18 (∼10 • ) for Re δ = 1450. Finally, as Re δ increases further, Δφ decreases again following the theoretical solution of Fredsøe (1984) (agreement seems excellent for Re δ ≥ 3000).
(ii) The present analysis concludes that the presence of a phase lag between the bed-shear stress and free-stream velocity maxima is the result of a late transition to a stage that mimics the characteristics of quasi-equilibrium conditions during the deceleration phase. However, this transition remains incomplete as neither the ensemble-average velocity profile (slope and intersect) nor the diagnostics for the quasi-equilibrium condition, namely the shape factor H and the defect shape factor G, reach their equilibrium values. Nevertheless, in these parts of the period, the TKE production to dissipation ratio Pr/ becomes ∼1, which would seem as a necessary but not sufficient condition for a logarithmic profile to exist in OBL flows.
(iii) Finally, DNSs performed for Re δ = 671 confirmed the sensitivity of flow behaviour to background disturbances. Specifically, two different initialization techniques were applied: one by using a 2.5 % disturbance following the approach by Ozdemir et al. (2014) and one by initializing the simulation using a fully turbulent field. In the first case, the flow turned laminar after a couple of periods. On the other hand, when the flow was initialized using a fully turbulent field, self-sustained turbulence was maintained. This confirms the findings by Ozdemir et al. (2014), who observed a crisis near the critical Reynolds number for intermittent turbulent OBL flows. This observation is also reasonably close to the predictions based on the stability analysis by Blennerhassett & Bassom (2002), who propose a critical Reynolds number of ∼709 for the transition. of the results for Re δ = 763, which was the most sensitive case that required the longest domain to become developed. In figure 22, the mean velocity profiles for three different polynomial orders, sixth to eighth, are presented. It is shown that there is practically no difference in the results as we increase or decrease the resolution. This is to be expected given the resolution reported in table 1. Finally, preliminary simulations using 5 periods in domain B were compared against the results from the final number of periods that was summarized in table 1. In the latter, the symmetry between ωt and ωt was also used to enhance the flow statistics. A typical comparison for the TKE statistics over 5 and 10 periods (20 realizations) is shown in figure 23. No significant difference was observed (error typically less than 1.5 %-3.0 %) when additional periods were considered for both the TKE characteristics and the bed-shear stress for all the examined cases. This may be expected given the significantly large area over which the spatial averaging is applied.