The horizontal dispersion of materials with a constant rising speed under the exclusive influence of ocean surface boundary layer (OSBL) flows is investigated using both three-dimensional turbulence-resolving Lagrangian particle trajectories and the classical theory of dispersion in bounded shear currents generalized for buoyant materials. Dispersion in the OSBL is caused by the vertical shear of mean horizontal currents and by the turbulent velocity fluctuations. It reaches a diffusive regime when the equilibrium vertical material distribution is established. Diffusivity from the classical shear dispersion theory agrees reasonably well with that diagnosed using three-dimensional particle trajectories. For weakly buoyant materials that can be mixed into the boundary layer, shear dispersion dominates turbulent dispersion. For strongly buoyant materials that stay at the ocean surface, shear dispersion is negligible compared to turbulent dispersion. The effective horizontal diffusivity due to shear dispersion is controlled by multiple factors, including wind speed, wave conditions, vertical diffusivity, mixed layer depth, latitude, and buoyant rising speed. With all other meteorological and hydrographic conditions being equal, the effective horizontal diffusivity is larger in wind-driven Ekman flows than in wave-driven Ekman–Stokes flows for weakly buoyant materials and is smaller in Ekman flows than in Ekman–Stokes flows for strongly buoyant materials. The effective horizontal diffusivity is further reduced when enhanced mixing by breaking waves is included. Dispersion by OSBL flows is weaker than that by submesoscale currents at a scale larger than 100 m. The analytic framework will improve subgrid-scale modeling in realistic particle trajectory models using currents from operational ocean models.
The transport and dispersion of dissolved and particulate materials in the ocean, such as salts, gases, marine organisms, and various marine pollutants, determine the pathways and the concentration of those materials, respectively. Some particulate materials, such as spilled oil, marine debris, and some marine organisms, are positively buoyant and are confined in the ocean surface boundary layer. Understanding the transport and dispersion of oceanic buoyant materials is therefore important to the marine environment, ocean ecology, and Earth’s climate.
Material transport and dispersion are more commonly and naturally understood in a Lagrangian reference frame than in an Eulerian reference frame. Material transport can easily be quantified using an ensemble of trajectories from a large number of individual particles, while material dispersion and mixing are usually understood using the ensemble of particle-pair displacements. Assuming that the separation of any two particles with trajectories X1 and X2 initially separated by D0 is , the relative dispersion of two particles initially separated by D0 is defined by (e.g., Özgökmen et al. 2012; Romero et al. 2013; Poje et al. 2014; Kamenkovich et al. 2015)
For anisotropic flows, the relative dispersion is often defined in two dimensions (i, j) as a tensor whose trace is the total relative dispersion . The relative diffusivity tensor is defined as
It is obvious from Eq. (2) that κ is a constant when D2 grows linearly with time and is otherwise a time dependent variable. Because of its importance in a range of environmental and geophysical problems in both the atmosphere and the ocean, material dispersion has been studied extensively using theoretical techniques, laboratory measurements, in situ observations, and computer simulations, after the seminal work by Richardson (1926). It is now understood that dispersion in natural environments is primarily caused by turbulent flows and by sheared mean currents. When the kinetic energy spectrum is steep with a slope smaller than or equal to −3, such as in the enstrophy cascade regime, dispersion is nonlocal as a result of strong turbulent strain fields and D2 grows exponentially with time (Bennett 1984). Within the inertial subrange with the famous −5/3 slope in the kinetic energy spectrum, D2 grows with t3 and is in the Richardson regime (Batchelor 1950). When particles are close so that their relative velocity is constant, D2 is proportional to t2 and is in the ballistic regime. When particles are far from each other so that the respective velocities of different particles are uncorrelated, dispersion is caused purely by random walk; D2 grows with t and is in the diffusive regime (e.g., Davidson 2006). For particles in a transversely sheared mean current, D2 grows with t3 when particles can extend in an infinite domain (e.g., Bennett 1987; LaCasce 2008). When particles are confined in a transversely sheared mean current, such as in a pipe or in an open channel, D2 grows with t (e.g., Taylor 1953; Fischer et al. 1979).
These theoretical principles have been used to understand and explain dispersion data from Lagrangian observing platforms in the ocean, such as drifters and floats (e.g., Davis 1991; Lumpkin et al. 2017) and from Lagrangian particle trajectories in realistic ocean model solutions. These studies cover the Atlantic Ocean (e.g., Lumpkin and Elipot 2010), the Pacific Ocean (e.g., Kirwan et al. 1978), and a number of marginal seas such as the Gulf of Mexico (e.g., LaCasce and Ohlmann 2003), the Adriatic Sea (Haza et al. 2008), the Gulf of La Spezia (e.g., Haza et al. 2010), the Nordic sea (Koszalka et al. 2009), and the Liguro–Provençal Basin (Schroeder et al. 2011). Those studies identify three distinct regimes for the evolution of particle dispersion: when particle-pair separation is below the deformation radius (viz., the enstrophy injection scale, typically from a few tens of kilometers to a hundred kilometers depending on latitudes; Chelton et al. 1998), dispersion is nonlocal and D2 grows exponentially with time. When particle separation is larger than the deformation radius and smaller than a few hundred kilometers, D2 is in the Richardson regime and grows with t3. When particle separation is larger than a few hundred kilometers, particle motions are uncorrelated and D2 grows linearly with time.
Horizontal particle dispersion at scales smaller than about 10 km was not studied until recently, partly because the physical processes at these scales (i.e., submesoscale and turbulent flows) were not well elucidated. Recent studies (e.g., Ohlmann et al. 2012; Romero et al. 2013; Poje et al. 2014; Beron-Vera and LaCasce 2016; D’Asaro et al. 2018) show that submesoscale motions are important in dispersion, although these studies still differ in whether the dispersion is local in the Richardson regime or is nonlocal with exponential growth in time.
A critical process missing in the abovementioned studies is three-dimensional small-scale boundary layer currents. While the effect of boundary layer turbulence is inherent to many dissolved and particulate oceanic substances, drifters and floats used for dispersion studies are usually advected by two-dimensional horizontal velocities at a fixed density surface or close to the air–sea interface. In numerical studies, existing studies of particle dispersion rely on particles advected by velocity fields from regional ocean models that do not resolve boundary layer turbulence and have a horizontal resolution from hundreds of meters to a few kilometers. Smaller-scale motions are represented by horizontally homogeneous and isotropic random motions. To the best knowledge of the authors, there are only a few papers studying the effect of small-scale oceanic processes on dispersion. Wang et al. (2018) examined the effect of internal waves on horizontal dispersion and found that D2 grows with t1.5. While studying the transport and dispersion by convective turbulence, Mensa et al. (2015) noted that a weak wind forcing can greatly enhance dispersion in convective turbulence. However, surface ocean boundary layer currents are dominated by winds and waves over much of the global ocean (e.g., Belcher et al. 2012; D’Asaro et al. 2014), and no existing study has focused on the effect of wind- and wave-driven boundary layer currents on horizontal material dispersion. One of the possible reasons is that boundary layer currents are usually studied with large-eddy simulation models that are configured on a computation domain of hundreds of meters in the horizontal directions (e.g., Skyllingstad and Denbo 1995; Tejada-Martínez and Grosch 2007) and therefore cannot be used to study material evolution over sufficiently large areas.
The purposes of this study are 1) to derive the analytic depth-averaged effective horizontal diffusivity for materials with a constant buoyant rising speed under the exclusive influence of ocean surface boundary layer flows; 2) to determine the respective importance of shear dispersion and turbulent dispersion in the ocean surface boundary layer; and 3) to study the effect of meteorological, geographic, and water-column conditions on the horizontal dispersion of oceanic buoyant materials. The remainder of the paper is organized as follows: section 2 presents a computational model and an analytic model for horizontal dispersion of buoyant materials with the detailed derivation of the analytic model in the appendix; section 3 presents results from both models and compares the two models; section 4 discusses the importance of the studied problem in the realistic ocean and the way the analytic model could be applied in realistic ocean models; and section 5 is a summary.
2. Model description
The study is based on two models: a large-eddy simulation model with Lagrangian particles that generates particle trajectories in ocean surface boundary layer turbulence and an analytic model of horizontal effective diffusivity for depth-averaged material concentration based on shear dispersion theory.
a. Large-eddy simulation model with Lagrangian particles
1) Model description
The modeling framework is based on the National Center for Atmospheric Research large-eddy simulation (NCAR-LES) model for the simulation of ocean surface boundary layer. The model equations (e.g., Craik and Leibovich 1976; McWilliams et al. 1997; Suzuki and Fox-Kemper 2016) and the numerical techniques (e.g., Sullivan et al. 1996; Sullivan and Patton 2011) are documented in previous publications and are not repeated here. The model has been used to study boundary layer turbulence driven by different surface meteorological and lateral large-scale flow forcings (e.g., Sullivan and McWilliams 2010; Sullivan et al. 2012; Hamlington et al. 2014; Sullivan and McWilliams 2018) and to study the evolution of different tracers (e.g., Liang et al. 2011, 2012, 2013; Brunner et al. 2015; Kukulka and Brunner 2015; Smith et al. 2016) and produces solutions that agree well with observations (e.g., Kukulka et al. 2009; Liang et al. 2017).
The location of buoyant Lagrangian particles is calculated by solving a set of ordinary differential equations for the trajectory of each individual particle:
where x is the Lagrangian trajectory, ur is the resolved velocity interpolated from the LES solutions of resolved velocity, ust is the wave-phase-averaged Stokes drift, wb is the buoyant rising speed, k is the unit vector in the vertical direction, and is the displacement due to subgrid motions. In the current study, the subgrid-scale displacement is calculated following the random displacement model (RDM) by Wilson (2015), which has been shown to be consistent with the Eulerian diffusion equation for material concentrations:
where ks,i is the subgrid-scale tracer diffusivity from the LES model, dξi is a Gaussian white noise with variance dt, and subscript i indicates spatial dimension. The contribution from subgrid-scale motion is negligible except close to the surface where turbulence is poorly resolved. The Lagrangian model embedded in the parallelized framework of the NCAR-LES model has been successfully applied to simulate the evolution of gas bubbles in the upper ocean using 8 million Lagrangian particles (Liang et al. 2017). Previously, Lagrangian particle models have also been used to study dispersion (e.g., Weil et al. 2004; Mueller and Veron 2009) and the modulation of momentum flux by sea spray (e.g., Richter and Sullivan 2013) in the atmosphere and boundary layer turbulence (Harcourt and D’Asaro 2010) and vertical mixing in the ocean (Kukulka et al. 2012a).
2) Lateral boundary conditions for Lagrangian particles
In most turbulence-resolving simulations of geophysical turbulence, domain periodicity has to be assumed in both horizontal dimensions. Since turbulence requires time and distance to develop, domain periodicity allows the computation of turbulence at relatively low computational cost and permits the study of turbulence evolution over a relatively long time period. While the computational domain has a finite size, the doubly periodic boundary implies an infinitely large horizontal domain and horizontal homogeneity at the scale of the computational domain for both the flow and tracers. Numerically, periodic boundary conditions facilitate the use of spectral methods that are more accurate than finite volume or finite different methods. While turbulence-resolving simulations in a doubly periodic domain have achieved considerable success in geophysical turbulence, particularly in the study of vertical processes in boundary layers, the finite-size domain, usually from hundreds of meters in the ocean to a few kilometers in the atmosphere, is a major obstacle to the study of material dispersion in the ocean that usually requires a large domain and a long integration.
To circumvent the finite-size domain limitation, Matheou and Bowman (2016) developed a recycling method for tracer dispersion modeling in atmospheric boundary layer. The method was later generalized to simulate tracer evolution in oceanic boundary layer by Chen et al. (2016). In this method, the periodic boundary conditions are retained for turbulence computation. Tracers are allowed to evolve out of the finite domain and are advected by periodic (recycling) turbulence fields. While the method does not require a large computation domain for turbulence, the computation domain for the tracer is determined by the spatial extents the simulated material occupies. In this study, the recycling method is applied to Lagrangian particles. Compared to the two previous studies, the application of the recycling method in Lagrangian particles does not require additional computing time when particles spread and therefore offers a more powerful approach to studying material dispersion.
3) Model configuration
The model is configured in a domain of 300 m × 1200 m × 150 m with 160 × 640 × 96 grid points. Wind speed at 10-m height U10 is set to 10 m s−1 in the positive x direction for the two primary experiments. Additional experiments with U10 = 5 and 15 m s−1 will also be presented for comparison. The mixed layer is initially 30 m deep. Two simulations, with and without nonbreaking surface gravity wave (Stokes drift) forcing, are performed. The Stokes drift due to surface gravity waves is calculated as (e.g., Sullivan et al. 2007; Harcourt and D’Asaro 2008) , where F(ω, θ) is the equilibrium wave directional spectrum proposed by Donelan et al. (1985) and later modified by Alves et al. (2003). For each set of forcing conditions, a spinup simulation without Lagrangian particles was first carried out for more than 1 day, followed by a run with Lagrangian particles for about 10 days. In each simulation, there are eight groups of Lagrangian particles, each group with different buoyancy rising speeds, (i.e., 1, 2, 4, 8, 10, 20, 40, and 80 mm s−1).
b. Analytic model for horizontal effective diffusivity by sheared boundary layer currents
We seek an analytic expression for asymptotic effective horizontal diffusivities in the ocean surface boundary layer flows. The study of one-dimensional (1-D) horizontal dispersion by shear current originates from Taylor (1953), who first derived an analytic solution for dispersion in a pipe. The analysis was extended to dispersion for atmospheric boundary layer flows (Saffman 1962), for open-channel flows (e.g., Fischer et al. 1979), and for a variety of steady and unsteady oceanic currents (e.g., Young et al. 1982; Young and Jones 1991). The effective horizontal dispersion coefficient for depth-averaged material concentration in vertically bounded flows with two horizontal dimensions can be written as (Esler and Ramli 2017)
where , the angle bracket is the depth average defined as , and kυ and kH are the vertical and horizontal diffusivity, respectively, which may include the effects from both molecular Brownian motions and turbulent fluctuation in the respective directions. Here, depth-averaged material concentration is advected by current () and is dispersed following the effective diffusivity in the horizontal directions. In Eq. (5), dispersion due to shear and that due to turbulent fluctuation are in the first and the second terms of the diagonal elements, respectively. According to Eq. (5), strong vertical shear of horizontal velocity and small effective vertical diffusivity are favorable conditions for strong shear dispersion. While Eq. (5) can be applied to neutrally buoyant materials in the ocean surface boundary layer, we here generalize the theory for oceanic materials with constant buoyant rising speeds.
1) Effective horizontal diffusivity for buoyant particles
Using the assumption for a random displacement model (e.g., Rodean 1996), the governing equation for the probability density function of particle concentration at location x and time t is given by the Fokker–Planck equation in two horizontal dimensions (Rodean 1996; Esler and Ramli 2017):
where . Here, the primary difference from Esler and Ramli (2017) is the addition of the fourth term on the left-hand side of the equation representing the buoyant rising effect. Horizontal turbulent diffusivity kH in Esler and Ramli (2017) is also generalized to , , and because horizontal velocity fluctuation in ocean boundary layer turbulence is anisotropic. The effective horizontal diffusivity for depth-averaged material concentration based on Eq. (6) can be derived using multiscale expansions and the method of moments (e.g., Saffman 1962) as (see appendix)
where the angle bracket is the depth average. In each component of , the two terms represent the effect of shear dispersion and turbulent dispersion, respectively. The F, M, and N are depth-dependent shape functions for mean concentration, and the perturbation concentrations in x and y directions are given by
where C is a constant to ensure , , and . Overbars in the above equations indicate averages weighted by F:
The mathematical derivation of Eqs. (7) and (8) is presented in the appendix. For neutrally buoyant tracers (wb = 0), F is a constant with depth, and Eqs. (7) and (8) reduce to Eq. (5) (see appendix). Equations (7)–(9) are valid when . The can equal any positive value when the expression for is more generally written as Eq. (A10).
Equations (7)–(9) form a theoretical framework to describe the evolution of buoyant materials under the exclusive influence of boundary layer currents. According to the framework, a patch of buoyant materials with a constant buoyant rising speed wb move at velocities , determined by Eq. (9) and spread at a rate calculated with Eqs. (7) and (8). The depth-averaged effective horizontal diffusivity [Eq. (7)] and velocity [Eq. (9)] are only accurate after the vertical profile of material concentration follows the equilibrated profile F(z) [Eq. (8a)]. The time scale for the establishment of F(z) can be estimated as , where is the smaller of the boundary layer depth and . When is smaller than the boundary layer depth , the time scale is proportional to and decreases with increasing . When is larger than or equal to , the time scale for the establishment of F(z) is independent of . Assuming that = 30 m and m2 s−1, the time scale is 9000 s when and is smaller than 9000 s when .
2) Diagnosis of shear-induced horizontal effective diffusivity using a 1-D ocean model
In this study, we use Eq. (7) to assess the effect of meteorological, geographic, and water-column conditions on horizontal dispersion of buoyant materials due to the vertical shear of boundary layer currents. When considering only the effective horizontal diffusivity due to the vertical shear of horizontal boundary layer current, Eq. (7) reduces to
The remaining unknowns in Eq. (10) are mean horizontal velocities (u and υ). In this study, we will present results using mean horizontal boundary layer currents from the LES model described in section 2a and from a 1-D model with a parameterized turbulent mixing effect. In a 1-D model, horizontal velocities in the ocean surface boundary layer are calculated as
Here, f is the Coriolis parameter (e.g., Cushman-Roisin and Beckers 2011). When the wave effect is not considered, and vertical diffusivity kυ is calculated by the K-profile parameterization (KPP) (e.g., Large et al. 1994):
where c1 = 0.4 is the von Kármán constant, is the ocean-side friction velocity, is the boundary layer depth that is chosen to be the smaller of the Ekman depth () and the mixed layer depth, and is the shape function with . According to Eq. (13), vertical diffusivity is controlled by friction velocity and boundary layer depth.
The is an enhancement factor proposed by McWilliams and Sullivan (2000) to account for the effect of Langmuir turbulence with , which is the turbulent Langmuir number (McWilliams et al. 1997); is a shape function proposed by McWilliams and Huckle (2006) to account for enhanced mixing near the ocean surface due to breaking waves, with and H the Heaviside step function. Equation (14) reduces to Eq. (13) when and . Six experiments are carried out to explore the effect of buoyant rising speed, wind, wave, stratification, and vertical diffusivity on effective horizontal diffusivity, and the parameters of the experiments are listed in Table 1.
Equations (8)–(11) are discretized in the vertical direction z into a set of linear equations. Horizontal velocities from Eq. (11) are plugged into Eqs. (8) and (9). Numerical solutions for F, M, and N can be obtained efficiently using either interactive commercial software such as MATLAB or the Linear Algebra Package (LAPACK) available in both Fortran and C++.
a. LES of Lagrangian particles with a constant rising speed
Upper-ocean turbulent flows have been extensively studied in the past few decades using observational and numerical techniques (e.g., Plueddemann et al. 1996; Smith 1998; Sullivan and McWilliams 2010; D’Asaro 2014, and the references therein). Here, we show only profiles of mean current and turbulent velocity variance (Fig. 1) to facilitate the discussion of horizontal material dispersion. In both simulations (Ekman turbulence and Langmuir turbulence), the surface current is to the right of the wind. The horizontal current rotates clockwise and weakens with increasing water depth (Figs. 1a and 1b). The vertical shear of the downwind current is weaker in the upper few meters in Ekman turbulence than in Langmuir turbulence, primarily because of the fast decay of the wave-phase-averaged Stokes drift in the upper few meters when the wave effect is included. In the rest of the boundary layer, vertical shear of both downwind and crosswind currents is stronger in Ekman turbulence than in Langmuir turbulence, since vertical mixing by locally shear-driven turbulence is weaker than by wave-driven nonlocal Langmuir cells. Horizontal velocity variances are also distinctly different between Ekman turbulence and Langmuir turbulence, consistent with previous studies (e.g., Li et al. 2005). Downwind velocity variance is larger than crosswind velocity variance in Ekman turbulence (Fig. 1c). In Langmuir turbulence, crosswind velocity variance dominates downwind velocity variance in much of the boundary layer because of the counterrotating downwind coherent vortices (Fig. 1d). Vertical velocity variance is larger in Langmuir turbulence than in Ekman turbulence. The mean material concentration decreases with depth owing to buoyant rising, and the material penetration decreases with increasing buoyant rising speed (Figs. 1e and 1f). For material of the same buoyant rising speed, its concentration decays faster in Ekman turbulence than in Langmuir turbulence.
The horizontal locations of particles three days after they are released at a region of 300 m × 300 m centered at the origin in an Ekman boundary layer are shown in Fig. 2. Particles of different buoyant rising speeds are depicted in different colors. Lagrangian particles are able to travel for tens of kilometers in horizontal directions, much larger than the size of a computational domain without incurring additional computational resources. This demonstrates the capability of the model to study horizontal upper-ocean dispersion. Particles of all buoyant rising speeds are advected in both the downwind and crosswind directions to the right of the wind. Both the downwind and crosswind movement are faster for particles of larger buoyant rising speed. This is because particles with stronger preferred upward motion spend more time in the upper portion of the boundary layer where both downwind and crosswind velocity is stronger. The dependence of particle transport is consistent with previous large-eddy simulation studies with Eulerian concentration models (e.g., Yang et al. 2015).
Horizontal particle dispersion in the Ekman boundary layer is strongly anisotropic for all buoyant rising speeds. However, the anisotropy decreases when wb > 4 mm s−1. The major axis of dispersion is about 30° to the right of the wind and rotates toward the downwind direction. For the four selected buoyant rising speeds, dispersion decreases with increasing wb. The dispersion is uniform horizontally for weakly buoyant particles (the left inset of Fig. 2). When wb increases, the effect of convergence and divergence by the two-dimensional horizontal turbulent currents becomes evident (the middle and right insets of Fig. 2) and particles aggregate into patches and streaks, though the streaks and patches are not as evident as in Langmuir turbulence. This is consistent with Sundermeyer et al. (2014), which shows surface material is aligned in streaks without significant wave forcing.
In a wave-driven Langmuir boundary layer, particles also move along and to the right of the wind (Fig. 3). Displacement in both the downwind and the crosswind directions is smaller than in Ekman currents for weakly buoyant particles ( mm s−1) that can be mixed into the boundary layer. The mean horizontal current is weaker in the wave-driven boundary layer, except close to the surface when Stokes drift dominates (Figs. 1a and 1b). Strongly buoyant particles (wb ≥ 40 mm s−1) move farther in the downwind direction than in an Ekman boundary layer because they stay at the surface and are advected primarily by the Stokes drift.
The dispersion of particles with similar buoyant rising speeds in a Langmuir boundary layer is visually less anisotropic than in an Ekman boundary layer. For strongly buoyant particles (wb = 40 mm s−1), dispersion along the major axis is close to that along the minor axis. The difference in spreading rate among particles of different buoyant rising speeds is not visually evident in the instantaneous snapshots and will be quantified in the next subsection. The major axis of dispersion is also about 30° to the right of the wind and rotates counterclockwise toward the direction of the wind. Three dispersion patterns—the diluted plume, the blurred plume, and the fingered plume—have been identified by Yang et al. (2014), who also developed a regime diagram matching the drift-to-buoyancy parameter () with the dispersion patterns. The three patterns are shown in the insets of Fig. 3 with parameter Db equal to 215.6, 21.6, and 5.4, respectively. The relation between parameter Db and the dispersion pattern is in agreement with the study by Yang et al. (2014). Weakly buoyant particles (wb = 5 mm s−1; Db = 215.6) are primarily under the influence of three-dimensional boundary layer currents, and they evenly disperse horizontally (the left inset of Fig. 3). Moderately buoyant particles (wb = 10 mm s−1; Db = 21.6) are partially influenced by two-dimensional horizontal boundary layer currents (the middle inset of Fig. 3), displaying a blurred pattern. Strongly buoyant particles (wb = 40 mm s−1; Db = 5.4) are aggregated in the convergence streaks of two-dimensional horizontal currents (the right inset of Fig. 3). The streaks are stronger and horizontally more separated in the crosswind direction in a Langmuir boundary layer than those in an Ekman boundary layer.
2) Horizontal relative dispersion and absolute diffusivity
To quantify the dispersion, the time evolution of particle-pair relative displacement and the asymptotic horizontal diffusivity are studied (Figs. 4 and 5). Particle pairs are tracked if their initial separation is less than 1 m. Figure 4a plots the time evolution of relative displacement square D2 in major and minor axes, respectively, for 80 000 particle pairs in an Ekman boundary layer. Relative dispersion D2 undergoes an initial exponential growth regime for the first tens of seconds. After approximately 100 s, D2 grows with t2, implying a ballistic dispersion regime. For weakly buoyant particles (wb ≤ 4 mm s−1), D2 grows with t1.0 after about 104 seconds. For particles of larger buoyant rising speeds (wb = 10 and 40 mm s−1), the diffusive dispersion regime starts earlier, at around 3000 s. An examination of the vertical profile time series shows that the asymptotic diffusive regime is reached when the vertical material profile is equilibrated.
To study the relative importance of the mean horizontal current and the turbulent velocity fluctuation, the displacement due to horizontally averaged velocity and turbulent velocity are respectively calculated. For the weakly buoyant particles (wb ≤ 4 mm s−1), D2 (dominated by D2major) under the sole influence of horizontal currents grows approximately with t2 from 102 to about 104 seconds after the initial exponential growth. During this period, particles are not bounded by the mixed layer, and the ballistic dispersion regime is due to unbounded shear dispersion. When the equilibrium vertical distribution is reached (t > 104 seconds), D2 growth is diffusive, consistent with the theory of bounded shear dispersion. The transition to the diffusive regime is earlier for particles with larger buoyant speed. This is qualitatively consistent with the discussion in section 2b(2) that the theoretical time scale required to establish an equilibrated vertical profile decreases with increasing buoyant rising speed. Relative dispersion D2 due to the vertical shear of mean currents changes significantly with wb and spans a few orders of magnitude (Fig. 4b). Similar to the time evolution of D2 due to mean current, D2 due to horizontal turbulent velocity also undergoes an initial exponential growth and then a ballistic growth before reaching the asymptotic diffusive regime in which turbulent velocity for different particles is uncorrelated (Fig. 4c).
In the diffusive regime, the absolute diffusivity (Fig. 4d) is half of the relative diffusivity given by Eq. (2) (e.g., LaCasce 2008). The diffusivity in major (Kmajor) and minor (Kminor) directions and the direction of the major axis of dispersion θ can be calculated as (e.g., Oh et al. 2000)
The absolute diffusivity in the major axis Kmajor is about 6.2 m2 s−1 and is more than 15 at wb = 2 mm s−1. It is the largest among all wb considered (solid black line in Fig. 4d). The absolute diffusivity in the minor axis Kminor is about one order of magnitude smaller than Kmajor. The ratio between Kmajor and Kminor is 68 when wb = 4 mm s−1, the largest among the selected wb, and is around 12 at wb = 40 mm s−1. The direction between the major axis of dispersion and wind θmajor is approximately 26.5° for wb = 1 mm s−1 and decreases to about 5° for large wb (≥10 mm s−1) (Fig. 4e, black line). The diffusivity calculated with Eq. (10) using mean currents and diagnosed vertical diffusivity from the LES-particle simulation (green lines in Figs. 4d and 4e) agrees reasonably well with those calculated with the particle statistics from the LES simulations (red lines in Figs. 4d and 4e), showing the ability of a reduced-physics model in accurately predicting horizontal dispersion due to vertical shear of boundary layer currents. The variability of Kmajor, Kminor, and θmajor with wb due to turbulent horizontal velocity is relatively small compared to the variability due to mean current shear. The angle between Kmajor and U10 (θmajor) is small (~5° to the right of the wind). When wb is smaller than about 10 mm s−1, Kmajor by the mean current is larger than that by turbulent horizontal velocity. For strongly buoyant particles, Kmajor, Kmajor, and θmajor are dominated by horizontal turbulent velocity.
Figures 5a to 5c plot the time evolution of relative dispersion for 80 000 particle pairs in a Langmuir boundary layer. Similar to in an Ekman boundary layer, dispersion in a Langmuir boundary layer reaches the asymptotic diffusive regime at around 103–104 seconds. Before the diffusive regime, relative dispersion undergoes an exponential growth stage and a power law growth stage with a power close to 2.
The asymptotic absolute diffusivity Kmajor increases with wb and reaches its maximum of 2.2 m2 s−1 (or about 6) when wb = 20 mm s−1. Both the absolute and the normalized maximum horizontal diffusivity is smaller in a Langmuir boundary layer than in an Ekman boundary layer. Dispersion is less anisotropic in a Langmuir boundary layer than in an Ekman boundary layer. The largest Kmajor/Kminor is less than 15 when wb = 20 mm s−1 (see Fig. 7d below). When wb = 40 mm s−1, Kmajor/Kminor is less than 1.5. The ratio increases toward wb = 80 mm s−1. The angle between Kmajor and U10 also is to the right of the wind when wb ≤ 20 mm s−1 and is to the left when wb > 20 mm s−1. When wb = 80 mm s−1, θmajor is about 80° to the left of the wind (Fig. 7e). Similar to the total relative dispersion, the relative dispersion due to the mean horizontal current also reaches the asymptotic diffusive regime after 103–104 seconds. The total dispersion and diffusivity are dominated by the shear-induced diffusivity when wb ≤ 20 mm s−1. When wb = 40 mm s−1, the contribution from the shear of mean current and that from turbulent velocity fluctuation are comparable. Particles with wb = 80 mm s−1 cannot be mixed into the boundary layer ,and dispersion is dominated by turbulent velocity fluctuation.
While Eq. (10) reasonably predicts the effective horizontal diffusivity due to the vertical shear of mean boundary layer currents, the relation between diffusivity due to turbulent fluctuating velocity and the horizontal effective diffusivity due to turbulence is still undecided. Therefore, four additional experiments were carried out using the LES-particle model driven by wind (U10) of 5 and 15 m s−1, with and without Stokes drift forcing. Other initial and boundary layer conditions are the same as the two experiments with U10 = 10 m s−1. Similar to when U10 = 10 m s−1, the horizontal effective diffusivity due to horizontal velocities is small compared to diffusivity due to shear (not shown). In an Ekman boundary layer, is scaled by with the ocean-side friction velocity and the unit length. When wb is large (wb > 20 mm s−1) and shear dispersion is negligible, the mean of , , and are (Figs. 6a and 6b)
In a Langmuir boundary layer, the horizontal effective diffusivity due to horizontal turbulent velocities is scaled by . When wb > 20 mm s−1, the mean of , , and are (Figs. 6c and 6d)
b. Horizontal effective diffusivity from a 1-D column model
1) Horizontal dispersion of buoyant materials in an Ekman layer
We first consider an ocean surface Ekman layer driven by a steady wind located at 10 m above the ocean surface of 10 m s−1 (U10 = 10 m s−1) at 45°N (simulation S1 in Table 1). For an Ekman layer bounded by rotation (unstratified; Fig. 7), the vertical shear of downwind current is much stronger near the surface than within the boundary layer, while the vertical shear of the crosswind current is relatively uniform within the boundary layer (Fig. 7a). Concentration decreases with water depth for positively buoyant materials, and the concentration decays faster for materials of larger buoyant rising speeds (Fig. 7b). Horizontal dispersion is highly anisotropic for materials of all buoyant rising speeds wb: the ratio between the diffusivity in the major axis and that in the minor axis is about 20 when mm s−1 and increases monotonically to more than 1000 when mm s−1. The diffusivity in the minor axis decreases monotonically with increasing buoyant rising speed. The diffusivity in the major axis increases with increasing buoyant rising speed and reaches the maximum, close to 12 m2 s−1, at mm s−1 (Fig. 7c). At this buoyant rising speed, materials concentrate at the depths where the vertical shear of the downwind speed is the largest (Fig. 1b). When mm s−1, decreases with increasing wb. The major axis of dispersion is about 45° to the right of the wind for the smallest rising speeds (wb = 0.5 mm s−1) and rotates toward the direction of the wind with increasing buoyant rising speed (Fig. 7d). The influence of the vertical shear of crosswind current υ decreases when wb increases, since the shear in υ is relatively uniform in the boundary layer.
2) Effect of water-column stratification
For an Ekman layer bounded by stratification (simulation S2 in Table 1 and Fig. 8), vertical viscosity and diffusivity kυ is smaller [Eq. (13)] than in simulation S1; therefore, velocity shear is stronger and material concentration decreases faster than in unstratified Ekman layer (Figs. 8a and 8b). Similar to in an unstratified Ekman layer, dispersion is strongly anisotropic (Fig. 8c) and Kmajor/Kminor is about twice as large as in unstratified Ekman layer. The effective horizontal diffusivity in both the major and minor axis is smaller than in an unstratified Ekman layer; however, the diffusivity normalized by is slightly larger than in an unstratified Ekman layer. For wb > 0.5 mm s−1, diffusivity decreases monotonically with increasing wb. The major axis of dispersion is also to the right of the wind (Fig. 8d). The misalignment of the major axis of dispersion with wind is smaller than in unstratified Ekman layer, implying a relatively larger influence by downwind current shear than by crosswind current shear.
3) Effect of vertical diffusivity
Vertical diffusivity kυ is a major uncertainty in ocean models, and solutions from sensitivity experiments with a doubled kυ by setting c1 = 0.8 in Eq. (13) are also presented (simulation S3 in Table 1 and Fig. 9). The vertical shear of both downwind and crosswind horizontal currents with a doubled kυ is weaker than with the standard kυ value (Fig. 9a vs Fig. 7a). For material of the same wb, concentration decays slower when kυ is larger (Fig. 9b vs Fig. 7b). Similar to the previous two experiments with standard kυ values (simulations S1 and S2), dispersion is strongly anisotropic: Kmajor/Kminor is about 25 when mm s−1. Although Kmajor/Kminor also increases with increasing wb, the increase is slower in simulation S3 than in S1 and S2. The Kmajor/Kminor becomes larger than 1000 for mm s−1. For weakly buoyant materials (wb < 8 mm s−1), when they can be mixed into the boundary layer in both simulations S1 and S3, effective horizontal diffusivity with doubled kυ is smaller than with the standard kυ value because of weaker shear and larger vertical diffusivity. Maximum horizontal diffusivity in the major axis is smaller than 2 m2 s−1, much smaller than that with the standard kυ value (>11 m2 s−1). Materials with wb ≥ 8 mm s−1 can be mixed into the boundary layer by the larger kυ in simulation S3 but not in simulation S1; horizontal diffusivity is larger in simulation S3 than in simulation S1.
4) Effect of latitude
The sensitivity of horizontal diffusivity to latitude is investigated through comparing two solutions at 15°N (simulation S4 and Fig. 10) and at 45°N (Fig. 7) with all other conditions the same. Here we set the mixed layer depth at 15°N the same as the Ekman depth at 45°N for the ease of comparison because the Ekman depth is deeper at 15° than at 45°N. Vertical shear of horizontal currents in both downwind and crosswind directions is evidently stronger at 15° than at 45°N (Fig. 10a vs Fig. 7a). The vertical profiles for materials of the same buoyant rising speed at the two different latitudes are the same (Fig. 10b vs Fig. 1b) because vertical diffusivity is the same. The effective horizontal diffusivity is larger at 15° than at 45°N in both the major and minor directions (Fig. 10c vs Fig. 1c) because of the stronger current shear at 15°N. The Kmajor/Kminor is larger at 15° than at 45°N for all buoyant rising speeds The misalignment of the major axis of dispersion with wind is smaller at 15° than at 45°N, implying a relatively larger influence by downwind current shear than by crosswind current shear.
5) Effect of nonbreaking surface gravity waves
Nonbreaking surface gravity waves directly modify horizontal current through the Stokes–Coriolis term (e.g., Polton et al. 2005). When with a substantial along-wind component, they also interact with the currents and generate Langmuir circulation that penetrates the boundary layer (e.g., D’Asaro 2014), resulting in vertical diffusivity larger than in shear-driven Ekman turbulence through the boundary layer (e.g., McWilliams and Sullivan 2000). Breaking waves impact the vertical profiles of horizontal currents and vertical diffusivity, mainly close to the surface (e.g., McWilliams et al. 2012). The effects of nonbreaking surface gravity waves are represented by the Stokes drift [ust in Eq. (11)] and by the enhanced vertical diffusivity [ε(Lat) in Eq. (14)]. We here consider only one scenario when the wave is in equilibrium with the wind (simulation S5 in Table 1). In the real ocean, waves can be of any magnitude and in any direction (e.g., van Roekel et al. 2012; McWilliams et al. 2014). Those scenarios can also be studied with the analytic framework.
Under the influence of nonbreaking surface gravity waves, the Lagrangian downwind current decays faster in the upper few meters and is more uniform inside the boundary layer than in the Ekman layer (Fig. 11a vs Fig. 7a). The vertical shear of the crosswind current is weaker throughout the boundary layer with waves than without. Materials of the same buoyant rising speed decay more slowly with depth with waves than without (Fig. 11b vs Fig. 7b), since vertical diffusivity is larger in a wave-driven boundary layer current than in an Ekman layer. Similar to in an unstratified Ekman layer, horizontal dispersion in a wave-driven boundary layer is highly anisotropic (Fig. 11c vs Fig. 7c). The Kmajor/Kminor is about 24 when wb = 0.5 mm s−1 and increases with increasing wb. The increase in Kmajor/Kminor with wb, however, is slower than without waves. Horizontal diffusivity in the major axis increases with increasing wb for weakly buoyant materials. The buoyant rising speed when Kmajor reaches maximum is about 13 mm s−1 with waves, larger than that in an Ekman layer, since vertical diffusivity is larger in wave-driven turbulence. The maximum horizontal diffusivity is smaller than in an Ekman layer. The misalignment of the major axis of dispersion with wind is smaller in a Langmuir boundary layer than in an Ekman boundary layer (Fig. 11d vs Fig. 7d), again implying a relatively larger influence by downwind current shear than by crosswind current shear.
6) Effect of breaking surface gravity waves
When the enhanced near-surface mixing due to breaking waves is included with nonbreaking waves in the model (simulation S6), the near-surface vertical shear of horizontal current reduces (Fig. 12a vs Fig. 11a). Buoyant materials penetrate deeper into the boundary layer with the inclusion of breaking waves than without the breaking wave effect (Fig. 12f vs Fig. 11b). Horizontal diffusivity is smaller with breaking waves than without breaking waves (Figs. 12c and 2c) since the vertical shear of horizontal currents is weaker and vertical diffusivity is larger. While Kmajor/Kminor is larger with breaking waves than without when wb = 5 mm s−1, the increase in Kmajor/Kminor is slower than without breaking waves. When wb is greater than 3 mm s−1, Kmajor/Kminor is smaller with breaking waves included with nonbreaking waves than without breaking waves, indicating weaker dominance of vertical shear of downwind current with breaking waves. The misalignment between the major axis of dispersion with wind is slightly larger with breaking waves included for all buoyant rising speeds considered (Fig. 12d vs Fig. 2d).
a. Importance of shear dispersion by ocean surface boundary layer currents to spilled oil and ocean plastic
Dispersion by the vertical shear of horizontal boundary layer currents is an order of magnitude larger than that by turbulent velocity fluctuations in wind- and wave-driven boundary layers, but it is only important for materials that are able to be mixed into the boundary layer by turbulence [Eq. (8) and results in section 3b]. That is when is comparable to the depth of the boundary layer . For an Ekman layer, scales with , with being the ocean-side friction velocity (e.g., Large et al. 1994). The ratio between and becomes . Friction velocity is 6 mm s−1 when U10 = 5 m s−1, following the drag law of Large and Pond (1981), implying that shear dispersion by the boundary layer current is significant when wb is about 6 mm s−1 or smaller at this wind speed. When the wind is stronger, materials with a larger wb can suspend in the boundary layer and be dispersed by the sheared current. There are two important pollutants in the ocean, including marine plastic and spilled oil, which usually have buoyant rising speed of a few millmeters per second. Many marine plastic debris are weakly buoyant with a rise speed of that range (e.g., Reisser et al. 2015). They are observed to be distributed well into the ocean surface boundary layer (e.g., Kukulka et al. 2012b; Reisser et al. 2015). After an oil spill event, chemical dispersants are usually applied to break oil slicks into small droplets of a hundred micrometers or less in diameter (e.g., Li and Garrett 1998), so that the small oil droplets can suspend in the water column and be degraded by bacteria (e.g., Joye et al. 2016). For oil droplets with density ρoil of 850 kg m−3 and a diameter doil of 100 μm, wb is 0.95 mm s−1 using Stokes’s law . Shear dispersion in the ocean surface boundary layer is therefore important for both ocean plastic and oil droplets.
b. Importance of horizontal dispersion by ocean surface boundary layer flows
Simulations in this study show that the horizontal diffusivity due to boundary layer processes is on the order of 100–101 m2 s−1 for a wind speed of 10 m s−1 located at 10-m height, with the variability influenced by latitude, boundary layer depth, bottom stratification, wave conditions, and particle buoyant rising speed. Dispersion due to geostrophic mesoscale ocean processes undergoes exponential growth when particle separation is smaller than the size of mesoscale eddies (from a few kilometers to tens of kilometers); that is, so that with an exponent λ on the order of 2 × 10−7–10−6 (Young et al. 1982). Therefore, the scale of diffusivity due to the boundary layer current is larger than that due to mesoscale processes at a scale of 1–10 km. A recent observational study (Poje et al. 2014) shows that dispersion by submesoscale currents can be an order of magnitude larger than that by mesoscale flows. The horizontal effective diffusivity increases with horizontal scale of the material and is larger than 10 m2 s−1 when the scale is larger than a hundred meters or so (Fig. 6 in Poje et al. 2014). The horizontal dispersion due to boundary layer currents is likely dominated by submesoscale and mesoscale currents at a scale larger than 100 m.
c. Parameterization of horizontal dispersion in the ocean surface boundary layer
With advances in computing technology and computational methods, current predictive ocean models usually resolve submesoscale currents (e.g., Haza et al. 2016; Barkan et al. 2017) but not boundary layer processes. Lagrangian particle–tracking models (e.g., North et al. 2011; Paris et al. 2013; Thyng and Hetland 2014) driven by velocity fields from predictive ocean models usually use a subgrid-scale model that does not include boundary layer currents and is independent to buoyant rising speed. The analytic model in section 2b(1) can be used to mechanistically determine subgrid-scale velocity in those models: the effective horizontal diffusivity is first calculated using Eqs. (7)–(9) and Eqs. (14)–(16), and the subgrid-scale displacement is then related to the diffusivity using a random displacement model (e.g., Griffa 1996), such as Eq. (4).
Horizontal diffusivity due to turbulent velocity fluctuation are important for strongly buoyant materials, although they are an order of magnitude smaller than that due to shear dispersion for weakly buoyant materials. It is given for Ekman turbulence () and for a Langmuir turbulence–driven wave in equilibrium with wind (). In the realistic ocean, the characteristic of turbulence is influenced by forcing not included here, such as another combination of wind and wave forcing, surface buoyancy flux (e.g., Li et al. 2005; Belcher et al. 2012), the angle between wind and wave (e.g., van Roekel et al. 2012), and wind direction (e.g., Liu et al. 2018). A relationship more generic than Eqs. (15) and (16) will only be possible through a systematic study that involves a comprehensive survey of LES solutions under different meteorological and hydrographic conditions. Future studies will also be carried out to investigate how boundary layer currents and submesoscale currents together modulate horizontal material dispersion.
d. Dispersion of buoyant materials with different rising speed wb by ocean surface boundary layer flows
This study focuses on the dispersion of particles with the same buoyant rising speed. Particles with different buoyant rising speeds move at different advection velocities [Eq. (9) and Figs. 2 and 3] because they occupy different portion of the water column. The differential advection separates particles with different buoyant rising speeds more than dispersion does.
Two different approaches are used to understand and quantify the effect of ocean surface boundary layer flows on the dispersion of materials with a constant upward buoyant speed. The first approach employs a three-dimensional Lagrangian particle–tracking model driven by turbulent velocity fields from a large-eddy simulation model. Model solutions are used to calculate time-dependent relative dispersion and asymptotic horizontal diffusivity. The second approach generalizes the classical theory on dispersion in a bounded shear flow to materials that experience a constant upward buoyant velocity and computes the asymptotic depth-averaged effective horizontal diffusivity for any combination of horizontal currents, vertical diffusivity, and buoyant rising speed.
The particle trajectory solutions show that particles of different buoyant rising speeds move and spread in different directions at different rates because they occupy different depths of the water column. Effective horizontal diffusivity is larger in an Ekman boundary layer than in a Langmuir boundary layer for weakly buoyant particles that can be mixed into a boundary layer by both types of turbulence. For buoyant particles that could be mixed into the boundary layer by Langmuir turbulence but not by Ekman turbulence, horizontal dispersion is stronger in a Langmuir boundary layer than in an Ekman boundary layer. Dispersion is anisotropic. In an Ekman boundary layer, the major axis of dispersion is about 30° to the right of the wind for weakly buoyant particles and rotates toward the direction of the wind with increasing buoyant rising speed. In a Langmuir boundary layer, the major axis of dispersion is also about 20° to the right for weakly buoyant particles. However, it is approximately 80° to the left of the wind for strongly buoyant particles because of the strong crosswind velocity fluctuation associated with Langmuir cells. Dispersion experiences an exponential growth and a ballistic regime, before transitioning to the asymptotic diffusive regime. The asymptotic regime is established when material distribution in the vertical direction reaches equilibrium. The horizontal effective diffusivity from shear dispersion theory agrees reasonably well with diffusivity calculated using 3D particle trajectories. For weakly buoyant materials that could be mixed into the boundary layer, shear dispersion dominates turbulent dispersion. For strongly buoyant materials that cannot be mixed into the boundary layer, turbulent dispersion dominates.
Using the shear dispersion theory with a Reynolds-averaged ocean model, we demonstrated that the variability of shear-driven effective horizontal diffusivities to buoyant rising speed, wind, wave, stratification, latitude, vertical mixing, and breaking waves. For particles with small buoyant rising speed that can be mixed into the boundary layer by turbulence, a deep boundary layer, weak vertical mixing, and low latitude favors horizontal dispersion by the vertical shear of boundary layer currents. For particles that are able to be mixed into the boundary layer by both shear-driven Ekman turbulence and wave-driven Langmuir turbulence, horizontal dispersion is stronger in Ekman turbulence than in Langmuir turbulence. For particles that can be mixed into the boundary layer only by Langmuir turbulence but not by Ekman turbulence, horizontal dispersion is stronger in Langmuir turbulence than in Ekman turbulence. Horizontal dispersion is also sensitive to latitude and vertical mixing. Vertical mixing determines horizontal dispersion by controlling both the vertical shear of horizontal velocity and the amount of materials that could be mixed into the boundary layer. While the KPP parameterization for vertical diffusivity kυ is used in section 3b(2), other parameterizations for kυ (e.g., Kantha and Clayson 1994; Umlauf and Burchard 2003; Harcourt 2013) can also be applied with Eqs. (10)–(12) to calculate the effective horizontal diffusivity. Recent efforts in improving the parameterization of vertical mixing (e.g., Harcourt 2015; Sinha et al. 2015; Reichl et al. 2016; Li and Fox-Kemper 2017; Taylor 2018) also improve the prediction of horizontal dispersion in ocean models.
The authors thank two anonymous reviewers for their comments and suggestions. Numerical solutions used for the plots are publicly available through the Gulf of Mexico Research Initiative Information and Data Cooperative (GRIIDC) (https://data.gulfresearchinitiative.org; https://doi.org/10.7266/N7CV4G9S). The study received support from National Key Research Program of China (2017YFA0604100). JHL is funded by the National Science Foundation (NSF) through Grants OCE-1521018 and OCE-1558317. JHL, KAR, and JCM are funded through a grant from the Gulf of Mexico Research Initiative. XW is funded by NSF Grant DMS-1620026 and AFOSR Grant FA9550-15-1-0051. PPS and JCM are funded by the Office of Naval Research through Grant N00014-12-1-0105. Computations were performed on supercomputing facilities at Louisiana State University, through the Louisiana Optical Network Infrastructure (LONI), and at the National Center for Atmospheric Research (NCAR). NCAR is sponsored by NSF.
Derivation of Eq. (7)
We start from the Fokker–Planck equation for the probability density function of particle concentration at location x and time t [Eq. (4)]:
subject to no-flux boundary conditions ( at ). We then perform the following coordinate transformation:
With respect to the new coordinates, Eq. (A1) can be written as
where is the differential of a variable in horizontal dimensions. Let , , and , where ε ≪ 1 is the ratio between the vertical and horizontal characteristic scales (Esler and Ramli 2017). Equation (A2) can be rescaled as
Since the physical meaning of () is the same as that of (), we will use () for the rescaled equation; that is,
Substituting the ansatz into Eq. (A3), we obtain a power series in terms of . The equation of order is
subject to boundary conditions . We can assume that , where satisfies the following equation:
subject to the boundary conditions
The equation of order is
We can assume that and obtain the following two equations for M and N, respectively:
subject to no-flux boundary conditions, similar to Eq. (A5). By the definition that and , , which is consistent with the no-flux boundary conditions, and the above equations are solvable.
The equation of order is
Averaging this equation from −h to 0 along the z direction, we have
where C is a constant chosen to ensure , with , and with . Using the fact that , , and , the effective Eq. (A9) can be written in the original coordinate system as
which is the governing equation for depth-averaged material concentration advected by current () and dispersed following a horizontal effective diffusivity .
According to Eq. (A4), F is a constant with z, and the solution of the above equation is