Exploring ship track spreading rates with a physics-informed Langevin particle parameterization

The rate at which aerosols spread from a point source injection, such as from a ship or other stationary pollution source, is critical for accurately representing subgrid plume spreading in a climate model. Such climate model results will guide future decisions regarding the feasibility and application of large-scale intentional marine cloud brightening (MCB). Prior modeling studies have shown that the rate at which ship plumes spread may be strongly dependent on meteorological conditions, such as precipitating versus non-precipitating boundary layers and shear. In this study, we apply a Lagrangian particle model (PM-ABL v1.0), governed by a Langevin stochastic differential equation, to create a simplified framework for predicting the rate of spreading from a ship-injected aerosol plume in sheared, precipitating, and non-precipitating boundary layers. The velocity and position of each stochastic particle is predicted with the acceleration of each particle being driven by the turbulent kinetic energy, dissipation rate, momentum variance, and mean wind. These inputs to the stochastic particle velocity equation are derived from high-fidelity large-eddy simulations (LES) equipped with a prognostic aerosol–cloud microphysics scheme (UW-SAM) to simulate an aerosol injection from a ship into a cloud-topped marine boundary layer. The resulting spreading rate from the reduced-order stochastic model is then compared to the spreading rate in the LES. The stochastic particle velocity representation is shown to reasonably reproduce spreading rates in sheared, precipitating, and non-precipitating cases using domain-averaged turbulent statistics from the LES.

1 Introduction

Lying beneath regions of large-scale subsidence and above cool subtropical eastern ocean basins, stratocumulus clouds exert a strong negative radiative forcing (cooling effect) on the climate, as these bright low clouds are able to efficiently reflect shortwave radiation back to space in comparison to the relatively dark ocean surface below (Hartmann and Short, 1980). Stratocumulus decks observed from space often reveal narrow bands of enhanced albedo (known as ship tracks) as a result of ship emissions (Conover, 1966). The underlying physical mechanism behind the observed cloud brightening in ship tracks is referred to as the Twomey effect, which describes the relationship between increasing the number of cloud condensation nuclei (CCN) and the corresponding higher albedo as the increased surface area of the resulting smaller and more numerous droplets act as a more reflective surface than a lower CCN environment (Twomey, 1974, 1977). However, the Twomey effect does not act in isolation, and various additional aerosol–cloud interactions may occur with an increase in CCN, such as reduced collision-coalescence efficiency and the suppression of precipitation (Albrecht, 1989).

Because the Twomey effect could be enhanced by increased liquid water path (LWP) and/or extended cloud lifetime associated with precipitation suppression, Latham (1990) postulated that the injection of sea salt aerosols into stratocumulus-topped boundary layers (marine cloud brightening; MCB) may be a viable method to offset a substantial portion of greenhouse warming and potentially circumvent climate tipping points, such as a collapse of the Arctic ice sheet (Rasch et al., 2009). However, the clouds do not always increase in brightness when aerosol concentrations are increased. For example, aerosol increases in non-precipitating clouds can lead to decreases in cloud fraction and liquid water path (Ackerman et al., 2004; Toll et al., 2017). In models, the smaller droplets (that result from the increase in CCN with added aerosols) more readily evaporate (Wang et al., 2003), in part because they remain near the cloud top for longer periods of time as a result of decreased sedimentation rates (Ackerman et al., 2004; Bretherton et al., 2007). The aforementioned properties of smaller droplets promote increased entrainment efficiency, which may decrease LWP and cloud fraction and act to darken the cloud (Wood, 2007). Observational estimates of the degree to which an aerosol perturbation may act to alter LWP shed little light on the aerosol–cloud response, with studies finding depleted LWP (Christensen et al., 2023; Sato et al., 2018; Diamond and Wood, 2020; Gryspeerdt et al., 2019; Segrin et al., 2007; Coakley and Walsh, 2002), augmented LWP (Gryspeerdt et al., 2019; McCoy et al., 2018), and background meteorological condition dependence (Bender and Sentelhas, 2018; Christensen and Stephens, 2011, 2012).

Given the muddled results from satellite retrievals and the grid spacing required to numerically resolve both ship tracks and mesoscale circulations, process studies regarding the response of low clouds to aerosol injections often employ large-eddy simulation (LES), capable of resolving fine-scale turbulent structures that are crucial for estimating local aerosol and microphysical process rates and cloud-top entrainment rates (Lewellen and Lewellen, 1998). Previous LES studies of ship tracks have shown the LWP response to be dependent on background aerosol concentrations, with clean boundary layers exhibiting larger LWP in the ship track region, while polluted boundary layers experience the opposite (Wang et al., 2011; Berner et al., 2015; Chun et al., 2023). In a LES study of an idealized summertime subtropical stratocumulus regime, Chun et al. (2023) found that regardless of background aerosol or free-tropospheric moisture, the Twomey effect remained larger than the cloud fraction and LWP adjustments in all cases and resulted in cloud brightening of varying magnitudes over a 2 d period. Considering these promising LES results and a greater understanding of the environments in which MCB would be most effective, global modeling efforts of aerosol injections are a critical component in determining the feasibility of large-scale deployment and illuminating potential downsides related to regional climate variability brought about by inhomogeneous radiative forcing (Latham et al., 2008; Jones et al., 2009; Rasch et al., 2009).

Global climate model (GCM) studies of regional aerosol and CCN perturbations of the eastern subtropical oceans have suggested that MCB may be able to offset much of (if not all) the warming from a doubling of CO₂ (Jones et al., 2009; Rasch et al., 2009; Hill and Ming, 2012; Ahlm et al., 2017). However, the nature of the regional and global climate responses may differ substantially across the perturbation strategies and the assumed size distributions of aerosols (Wood, 2021). One major complicating factor of using GCMs to probe MCB feasibility is that they suffer from insufficient low cloud cover over the eastern subtropical ocean regions in comparison to observations (Xie et al., 2018), primarily as a result of under-resolved vertical and horizontal processes (Lee et al., 2022). Global-scale high-resolution (2–5 km horizontal grid spacing) simulations have recently become possible; however, such configurations are only practical for the simulation of a few days to months (Khairoutdinov et al., 2022) and not the decades or centuries that are required of GCMs. Future GCM runs of MCB strategies that do not involve instantaneous perturbations of entire ocean basins will require information about how injected aerosol plumes spread within the grid, and such information will need to be relayed to the radiation and microphysics parameterizations to better capture the spatial heterogeneity at small scales. In addition to MCB, a computationally efficient model of aerosol spreading tied to turbulent dynamics may be broadly applicable to injected stratospheric aerosols and the spread of hazardous chemicals and aerosols.

Previous attempts to constrain ship track spreading rates from satellite images aimed to estimate average lateral plume-spreading behavior (Durkee et al., 2000; Patel and Shand, 2022), but assuming a constant rate of plume spreading may result in subgrid plume fraction imprecision that leads to compounding nonlinear errors in the resolved-scale properties, as prior LES modeling studies have shown that the rate at which ship plumes spread may be strongly dependent on meteorological conditions, such as precipitating versus non-precipitating boundary layers (Prabhakaran et al., 2024) or wind shear (Berner et al., 2015). Recent efforts to represent subgrid plumes, such as the plume-in-grid (PIG) method with adaptive grids (Sun et al., 2022), allow for time-dependent changes in the horizontal spreading rate as a function of wind shear but require grid refinements in the presence of plumes. The associated increase in computational demand may be a bottleneck for the assessment of large-scale injection strategies. Alternatively, by leveraging Lagrangian particle-based methods, subgrid plumes may be characterized by statistical descriptions of the flow field at singular points in space and time and not bound by standard Gaussian diffusion and dispersion processes that assume fixed plume dispersion rates (Pope, 2000). The utilization of scalable approaches necessitates an accurate and computationally efficient representation of subgrid particle trajectories.

In this study, we formulate a turbulence-driven Lagrangian particle model, governed by a Langevin stochastic differential equation, to create a simplified and computationally efficient framework for predicting the rate of horizontal spreading from a ship-injected aerosol plume. The velocity and position of each stochastic particle is controlled by the acceleration of the surrounding fluid, which depends on the turbulent kinetic energy (TKE), dissipation rate, momentum variance, and mean wind. These inputs to the stochastic particle velocity equation are derived from high-fidelity and large-domain (204.8 km × 25.6 km) large-eddy simulations (LES) equipped with a prognostic aerosol–cloud microphysics scheme (Berner et al., 2013) to simulate aerosol injection from a ship into a cloud-topped marine boundary layer. By minimizing the error between the Gaussian fits of the Lagrangian particle model and LES ship track widths, we constrain a free parameter within the particle model established in Pope (2000). Using the fully parameterized particle model, we then study horizontal ship track spreading rates across a range of plausible vertical shear magnitudes in the northeastern Pacific boundary layer under precipitating and non-precipitating conditions. Conditionally averaged turbulent statistics from both within the ship plume and across the entire domain are used to judge spreading rate sensitivity and behavior in the particle model.

The paper is organized as follows: Sect. 2 describes the theoretical background of the Langevin particle model and details the numerical implementation of the approach. Section 3 outlines the LES configurations, Langevin particle model input parameters, and the plume width calculation used to compare the LES to the particle model. Section 4 compares the LES and particle model plume widths across the different shear cases and discusses the input parameters needed for optimal performance. Section 5 summarizes and discusses the findings of this work.

2 Description of plume model

This work aims to efficiently model the turbulent dispersion of atmospheric aerosols in the marine boundary layer for potential use as a subgrid plume parameterization in a GCM. In service of that effort, we consider a simplified but appropriately parameterized Langevin model and employ notation that is largely based on that used in Pope (2000). Specifically, we consider the unbounded (d= one-, two-, three-dimensional) turbulent dispersion of a passive conserved (no sources or sinks) scalar tracer, i.e., ϕ(x,t) with arbitrary units (e.g., temperature, concentration) and a known initial condition, i.e., ϕ₀(x). Being passive, this tracer has no effect on the material properties of the air or flow field through which it is transported (i.e., density, ρ; kinematic viscosity, ν; molecular diffusion, Γ) and thus has no effect on its own transport mechanism.

In the following sections, we will lay out the equations that govern our atmospheric plume model. We will begin with the Eulerian formulation representative of the LES framework and from there work towards the Lagrangian formulation that corresponds to the numerical particle model we introduce in Sect. 2.3.

2.1 Governing equations

The Eulerian conservation equation that governs the evolution of the scalar tracer field is the advection–diffusion equation

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial \mathit{\varphi }}{\partial t}}+\mathrm{\nabla }\cdot \left(\mathit{U}\mathit{\varphi }\right)=\mathrm{\Gamma }{\mathrm{\nabla }}^{\mathrm{2}}\mathit{\varphi },\end{array}$$

with the initial condition at t₀=0

$$\begin{array}{}\text{(2)}& \mathit{\varphi }(\mathit{x},t_{\mathrm{0}})={\mathit{\varphi }}_{\mathrm{0}}\left(\mathit{x}\right).\end{array}$$

The unbounded domain assumption, along with the diffusive nature of the spreading process, implies the following decay condition and long-time asymptotic solution

$$\begin{array}{}\text{(3)}& {\displaystyle }\underset{x_i\to \mathrm{\infty }}{lim}\mathit{\varphi }(\mathit{x},t)& {\displaystyle }=\mathrm{0},t\ge \mathrm{0},i\in \mathit{\{}\mathrm{1},\mathrm{2},\mathrm{3}\mathit{\}},\text{(4)}& {\displaystyle }\underset{x_i\to -\mathrm{\infty }}{lim}\mathit{\varphi }(\mathit{x},t)& {\displaystyle }=\mathrm{0},t\ge \mathrm{0},i\in \mathit{\{}\mathrm{1},\mathrm{2}\mathit{\}},\text{(5)}& {\displaystyle }\underset{x_{\mathrm{3}}\to \mathrm{0}}{lim}\mathit{\varphi }(\mathit{x},t)& {\displaystyle }=f(\mathit{y},t),t\ge \mathrm{0},\mathit{y}\in {\mathbb{R}}^{\mathrm{2}},\text{(6)}& {\displaystyle }\underset{t\to \mathrm{\infty }}{lim}\mathit{\varphi }(\mathit{x},t)& {\displaystyle }=\mathrm{0},\mathit{x}\in {\mathbb{R}}^{\mathrm{3}}.\end{array}$$

Here, $\mathit{U}(\mathit{x},t)\in {\mathbb{R}}^{\mathrm{3}}\left[{\text{m\hspace{0.17em}s}}^{-\mathrm{1}}\right]$ is the fluid velocity, $\mathrm{\Gamma }\in \mathbb{R}\left[{\text{m}}^{\mathrm{2}}{\text{s}}^{-\mathrm{1}}\right]$ is the constant diffusion coefficient (with the exact formulation depending on the definition of ϕ), and f(y,t) is an arbitrary boundary condition at the earth's surface (x₃=0). Note that we formulate this for d=3, although formulations for $d=\mathrm{1},\mathrm{2}$ may be written similarly.

The evolution of the fluid velocity field, U(x,t) is governed by the Navier–Stokes equations

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{\mathrm{D}\mathit{U}}{\mathrm{D}t}}=-{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}\mathrm{\nabla }p+\mathit{\nu }{\mathrm{\nabla }}^{\mathrm{2}}\mathit{U}+F_s,\end{array}$$

where we make the assumption that the base-state fluid density only varies in the vertical, ρ₀(z) [kg m⁻³] (anelastic approximation). Within Eq. (7), the kinematic viscosity is defined as $\mathit{\nu }:=\mathit{\mu }/\mathit{\rho }\left[{\text{m}}^{\mathrm{2}}{\text{s}}^{-\mathrm{1}}\right]$, wherein μ is the coefficient of viscosity that is assumed to be a constant in a Newtonian fluid. The mean fluid pressure is denoted by $p\left[\text{kg}{\text{m}}^{-\mathrm{1}}{\text{s}}^{-\mathrm{2}}\right]$, and any additional sources or sinks of momentum (e.g., apparent forces such as the Coriolis force, centrifugal force, or the buoyancy contributions in the vertical velocity component) are captured by the forcing term, F_s. Finally, the $\mathrm{D}/\mathrm{D}t$ operator is the material derivative defined as

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{\mathrm{D}}{\mathrm{D}t}}:={\displaystyle \frac{\partial }{\partial t}}+\mathit{U}\cdot \mathrm{\nabla }.\end{array}$$

By introducing a filtering operation that separates resolved fluid motion (U_i) from unresolved subgrid motion (u_i), Eq. (7) becomes the Reynolds-averaged Navier–Stokes (RANS) equation

$$\begin{array}{}\text{(9)}& {\displaystyle \frac{\partial U_i}{\partial t}}+U_j{\displaystyle \frac{\partial U_i}{\partial x_j}}=-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\displaystyle \frac{\partial p}{\partial x_i}}+\mathit{\nu }{\displaystyle \frac{{\partial }^{\mathrm{2}}{U}_i}{\partial x_j\partial x_j}}-{\displaystyle \frac{\partial {\mathit{\tau }}_{ij}}{\partial x_j}}+F_s,\end{array}$$

where the subgrid Reynolds stresses, ${\mathit{\tau }}_{ij}:=\stackrel{\mathrm{‾}}{u_i{u}_j}\left[{\text{m}}^{\mathrm{2}}{\text{s}}^{-\mathrm{2}}\right]$, represent a sink of momentum brought about by unresolved (subgrid-scale) fluctuations in velocity. A turbulence closure for τ_ij is necessary to close the system of equations given in Eq. (9). Once a relationship between the mean flow and subgrid flow has been established and an equation for pressure and conservation of energy have been solved, the mean fluid velocity can be updated in time and space, but as is the case with LES, this requires substantial computing power at fine spatial and temporal resolution, and we instead wish to describe the flow field that dictates tracer transport in a simplified manner.

2.2 Lagrangian governing equations

Given a turbulent flow field, defined in terms of its mean velocity $\stackrel{\mathrm{‾}}{\mathit{U}}(\mathcal{X},t)$, over the region 𝓧⊂ℝ³, Reynolds stresses τ_ij, and scalar dissipation rate ϵ [m² s⁻³], we are concerned with modeling the mean field of our aerosol plume $\stackrel{\mathrm{‾}}{\mathit{\varphi }}(\mathit{x},t)$ based on its initial condition ϕ₀(x). The Lagrangian nature of our model indicates that rather than describing our system with continuous fields in time and space, we will instead follow parcels of fluid that are characterized by their position and velocity, and we deploy the notation ${\mathit{X}}^{\left(n\right)}\left(t\right),{\mathit{U}}^{\left(n\right)}\left(t\right),n=\mathrm{1},\mathrm{\dots },N_{\mathrm{p}}$. Under this notational convention, N_p is the number of particles in the system, and each carries an equal portion (m_p) of the total mass of fluid ℳ contained in a fixed volume 𝒱 such that

$$\begin{array}{}\text{(10)}& m_p:={\displaystyle \frac{\mathcal{M}}{N_{\mathrm{p}}}}\equiv {\displaystyle \frac{\mathit{\rho }\mathcal{V}}{N_{\mathrm{p}}}}.\end{array}$$

We note that the assumption of equal mass for all particles is convenient but not required, as many Lagrangian models allow for unequal particle masses that may also vary in time (Monaghan, 2012; Tartakovsky et al., 2016; Avesani et al., 2015; Cherfils et al., 2012; Bosler et al., 2017; Schmidt et al., 2020, 2019; Engdahl et al., 2019; Jiao et al., 2022).

The particles in this model are all independently and identically distributed; i.e., two particles beginning at the same position and under the same fluid conditions (or at identical times) will be governed by equivalent position and velocity densities and thus have the same underlying probability density function (PDF). Further, their corresponding trajectories do not depend on one another. For this reason, we only need consider the dynamics of a single, arbitrary particle to describe this model, and we will denote its properties as ${\mathit{X}}^{*}\left(t\right),{\mathit{U}}^{*}\left(t\right)$. If we accept the argument of Taylor (1921) that the molecular diffusion provides a negligible contribution to the transport of the aerosol plume for the considered problem of atmospheric transport with a relatively high Reynolds number (Re) compared to the mean fluid flow and turbulent motions. In the case where molecular diffusion is neglected, Eq. (6) is no longer valid, given that the plume will fail to spread indefinitely in a flow with no turbulence. As such, the tracer ϕ is conserved along the path of a fluid parcel (particle), and the evolution of the mean aerosol plume field $\stackrel{\mathrm{‾}}{\mathit{\varphi }}$ is fully specified by the statistical properties of the fluid particles in motion. It is then convenient that the particle velocity PDF $f^{*}(\mathcal{U}\mathrm{|}\mathit{x};t)$ is equivalent to the fluid velocity PDF $f(\mathcal{U};\mathit{x},t)$. Note here that 𝓤 is the independent, or dummy, variable for the PDF of the velocity random variable, U^*. For example, when employing the notation more common to probability theory, the PDF could also be written as $f_{{\mathit{U}}^{*}}({\mathit{u}}^{*}\mathrm{|}\mathit{x};t)$. In addition, note that f^* is a density for the velocity U^* conditioned on the particle being located at position x, which is also a function of the time t at which the velocity is sampled, whereas f is a density for U but is not conditional and is strictly a function of x and t. This directly implies that the first and second moments of the particle velocity are equal to the mean fluid velocity and Reynolds stresses, as given in Table 1.

Table 1 Relationship between velocity moments and mean velocity and Reynolds stresses.

Due to the properties presented above, as well as other correspondences between the fluid and particle systems that are given in Pope (2000, Table 12.1), we obtain the governing equations for the particles. The position of a particle evolves according to

$$\begin{array}{}\text{(11)}& {\displaystyle \frac{\mathrm{d}{\mathit{X}}^{*}\left(t\right)}{\partial t}}={\mathit{U}}^{*}\left(t\right),\end{array}$$

with the particle velocity obeying an Ornstein–Uhlenbeck (diffusion) process defined by the stochastic differential equation (SDE)

$$\begin{array}{}\text{(12)}& \mathrm{d}{\mathit{U}}^{*}\left(t\right)=\mathit{a}\left({\mathit{U}}^{*}\left(t\right),{\mathit{X}}^{*}\left(t\right),t\right)\mathrm{d}t+b\left({\mathit{X}}^{*}\left(t\right),t\right)\mathrm{d}\mathit{W}\left(t\right).\end{array}$$

Here, $\mathit{a}\left(\mathcal{U},\mathit{x},t\right)$ and b(x,t) are generically formulated drift and diffusion functions, respectively, and W(t) denotes the d-dimensional (standard) Brownian motion. The Fokker–Planck equation associated with Eqs. (11) and (12) governing the behavior of the particle velocity–position joint PDF, $f_P^{*}(\mathcal{U},\mathit{x};t)$, is

$$\begin{array}{}\text{(13)}& \begin{array}{rl}& {\displaystyle \frac{\partial f_P^{*}}{\partial t}}+{\mathcal{U}}_i{\displaystyle \frac{\partial f_P^{*}}{\partial x_i}}=\\ & -{\displaystyle \frac{\partial }{\partial {\mathcal{U}}_i}}\left[f_P^{*}{a}_i\left(\mathcal{U},\mathit{x},t\right)\right]+{\displaystyle \frac{{\left[b\left(\mathit{x},t\right)\right]}^{\mathrm{2}}}{\mathrm{2}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{f}_P^{*}}{\partial {\mathcal{U}}_i\partial {\mathcal{U}}_i}}.\end{array}\end{array}$$

Note that in Eq. (13) (and going forward) we apply Einstein notation to represent sums over the spatial coordinate indices $\mathit{\{}i,j,k\mathit{\}}$. Finally, after traversing a small jungle of mathematical manipulations (Pope, 2000), we arrive at the generalized Langevin model (GLM) for the particle velocity stochastic process

$$\begin{array}{}\text{(14)}& \begin{array}{rl}\mathrm{d}{U}_i^{*}\left(t\right)& =-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\displaystyle \frac{\partial P}{\partial x_i}}\mathrm{d}t+{\mathit{G}}_{ij}\left[U_j^{*}\left(t\right)-\stackrel{\mathrm{‾}}{U_j^{*}\mathrm{|}{\mathit{X}}^{*}\left(t\right)}\right]\mathrm{d}t\\ & +\sqrt{C_{\mathrm{0}}\mathit{ϵ}\left({\mathit{X}}^{*}\right(t),t)}\mathrm{d}{W}_i\left(t\right).\end{array}\end{array}$$

This SDE contains some previously undefined terms, including the particle pressure field P(x,t) that directly corresponds to the mean fluid pressure $\stackrel{\mathrm{‾}}{p}(\mathit{x},t)$. Additionally, we introduce the constant C₀ and the time-variable drift coefficient ${\mathit{G}}_{ij}\left({\mathit{X}}^{*}\left(t\right),t\right)$, where G_ij depends on local values of Reynolds stresses ($\stackrel{\mathrm{‾}}{u_i{u}_j}$), ϵ (a function of X^* and time), and cross-derivatives of the mean velocity $\left(\partial \stackrel{\mathrm{‾}}{U_i}/\partial x_j\right)$. The reason Eq. (14) is referred to as generalized is because it defines a class of models that are specified according to choices for G_ij and C₀.

2.2.1 Simplified Langevin model

There are some attractive properties of boundary-layer flows that allow us to eliminate terms from the GLM and arrive at a modification of what is referred to as the simplified Langevin model (Pope, 2000). First we assume that the large-scale (mean) pressure gradients are small enough to be negligible in the boundary layer, giving

$$\begin{array}{}\text{(15)}& {\displaystyle \frac{\partial P}{\partial {\mathit{x}}_i}}\equiv \mathrm{0}.\end{array}$$

Second, we impose an isotropic weight coefficient for the drift term, namely

$$\begin{array}{}\text{(16)}& {\mathit{G}}_{ij}:=-{\displaystyle \frac{\mathrm{3}}{\mathrm{4}}}{C}_{\mathrm{0}}{\displaystyle \frac{\mathit{ϵ}}{k}}{\mathit{\delta }}_{ij},\end{array}$$

where ϵ and k [m² s⁻²] are the scalar dissipation rate and turbulent kinetic energy (TKE), respectively.

Equivalently, Eq. (16) can be formulated as a constraint on the evolution of kinetic energy in homogeneous turbulence,

$$\begin{array}{}\text{(17)}& {\displaystyle \frac{\mathrm{3}}{\mathrm{2}}}{C}_{\mathrm{0}}\mathit{ϵ}+{\mathit{G}}_{ij}\left(\stackrel{\mathrm{‾}}{u_i{u}_j}\right)=\mathrm{0}.\end{array}$$

This constraint also serves to define the Lagrangian integral timescale, T_L, that can be related to the system's TKE, dissipation rate, and the scalar Langevin isotropic drift coefficient, 𝒢, as follows:

$$\begin{array}{}\text{(18)}& T_{\mathrm{L}}={\displaystyle \frac{k}{\frac{\mathrm{3}}{\mathrm{4}}{C}_{\mathrm{0}}\mathit{ϵ}}}=-{\displaystyle \frac{{\mathit{\delta }}_{ij}}{{\mathit{G}}_{ij}}}:=-{\mathcal{G}}^{-\mathrm{1}}.\end{array}$$

T_L is also referred to as the “relaxation timescale” for turbulent mixing and spreading because it characterizes the timescale over which turbulent fluctuating velocity reverts to the background mean velocity. Lastly, in isotropic turbulence, the velocity variance σ² can be related to the turbulent kinetic energy k as

$$\begin{array}{}\text{(19)}& {\mathit{\sigma }}^{\mathrm{2}}={\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}k,\end{array}$$

which is equivalently stated as

$$\begin{array}{}\text{(20)}& {\displaystyle \frac{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}{T_{\mathrm{L}}}}=C_{\mathrm{0}}\mathit{ϵ}.\end{array}$$

Taken together and imposing a scalar isotropic G≡𝒢δ_ij, Eqs. (15)–(18) result in the Langevin model we consider in this work, namely

$$\begin{array}{}\text{(21)}& {\displaystyle }\textcolor{red}{}\text{(22)}& {\displaystyle }={\displaystyle \frac{\stackrel{\mathrm{‾}}{U_i^{*}}-U_i^{*}}{T_{\mathrm{L}}}}\mathrm{d}t+\sqrt{{\displaystyle \frac{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}{T_{\mathrm{L}}}}}\mathrm{d}{W}_i\left(t\right),\end{array}$$

or in a more compact vectorized notation,

$$\begin{array}{}\text{(23)}& \mathrm{d}{\mathit{U}}_t^{*}=\underset{\text{deterministic drift term}}{\underbrace{{\displaystyle \frac{\stackrel{\mathrm{‾}}{{\mathit{U}}_t^{*}}-{\mathit{U}}_t^{*}}{T_{\mathrm{L}}}}\mathrm{d}t}}+\underset{\text{Brownian motion term}}{\underbrace{\sqrt{{\displaystyle \frac{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}{T_{\mathrm{L}}}}}\mathrm{d}{\mathit{W}}_t}}.\end{array}$$

Here, the first bracketed term captures the memory effects of the relaxation timescale over which the Lagrangian velocity returns to the mean, while the second bracketed term is the Brownian motion contribution. Conceptually, the deterministic drift term is related to the representation of larger-scale flow features, while the Brownian motion term attempts to represent the random, smaller-scale turbulent fluctuations. The structure of Eq. (23) permits a broader range of behaviors than a traditional purely diffusive Gaussian plume model that represents turbulent mixing as a constant eddy diffusivity and is restricted to $\sqrt{t}$ growth (as discussed in Sect. 4.1). While the above equation is able to represent three-dimensional (d=3) flows, we restrict the following analysis to the horizontal dimensions (in particular, the x dimension) by vertically averaging boundary-layer quantities. Vertical averaging simplifies the particle model and focuses on horizontal plume spreading, which is the most GCM-relevant variable given the vertical spreading of the plume throughout the boundary-layer depth occurs on much shorter spatial and temporal scales. The vertical dimension reduction contains the implicit assumption that the boundary layer remains in an approximately well-mixed state, which is oftentimes the case in shallow, cloud-topped marine boundary layers.

2.3 Numerical implementation of Langevin particle model

The particle model we consider is composed of a collection of N_p particles that approximate the initial condition (IC) of the passive, conserved tracer Eq. (2) as a field composed of weighted kernel functions k(x,y). First, we note that hereafter we abandon the single arbitrary particle analysis applied in Sect. 2.2 and instead consider a full ensemble of particles identified by subscript index $i\in [\mathrm{1},N_{\mathrm{p}}]$. We formulate the model for the previously described boundary-layer flows in d=2 spatial dimensions in the horizontal or x–y plane. The IC may be formulated to be

$$\begin{array}{}\text{(24)}& {\displaystyle }{\mathit{\varphi }}_{\mathrm{0}}\left(\mathit{x}\right)& {\displaystyle }=\sum _{i=\mathrm{1}}^{N_{\mathrm{p}}}\underset{{\mathbb{R}}^{\mathrm{2}}}{\int }{m}_{i}k(\mathit{x}-\mathit{z})\mathit{\delta }(\mathit{z}-{\mathit{X}}_i^{*}(\mathrm{0}\left)\right)\mathrm{d}\mathit{z}\text{(25)}& {\displaystyle }& {\displaystyle }=\sum _{i=\mathrm{1}}^{N_{\mathrm{p}}}{m}_{i}k(\mathit{x}-{\mathit{X}}_i^{*}(\mathrm{0}\left)\right).\end{array}$$

We specify here that these kernels possess the standard properties that they are symmetric, translation invariant, non-negative, and integrate to unity; i.e.,

$$\begin{array}{}\text{(26)}& {\displaystyle }k(\mathit{x}-\mathit{y})\equiv k(\mathit{y}-\mathit{x}):=k\left(\mathit{z}\right),\text{(27)}& {\displaystyle }k\left(\mathit{z}\right)\ge \mathrm{0},\forall \mathit{z}\in {\mathbb{R}}^{\mathrm{2}},\text{(28)}& {\displaystyle }\underset{{\mathbb{R}}^{\mathrm{2}}}{\int }k\left(\mathit{z}\right)\mathrm{d}\mathit{z}=\mathrm{1}.\end{array}$$

The evolution of the tracer field is driven by the fluctuating Lagrangian velocities that follow Eq. (23), and the particle positions change according to

$$\begin{array}{}\text{(29)}& {\displaystyle \frac{\mathrm{d}{\mathit{X}}_i^{*}\left(t\right)}{\mathrm{d}t}}={\mathit{U}}_i^{*}\left(t\right).\end{array}$$

To solve the Langevin model (LM) governed by Eq. (23), we integrate the Langevin equation for velocity over a time step of length Δt employing the Euler–Maruyama method

$$\begin{array}{}\text{(30)}& \begin{array}{rl}{\mathit{U}}_i^{*}\left(t\right)& ={\mathit{U}}_i^{*}(t-\mathrm{\Delta }t)+{\displaystyle \frac{\stackrel{\mathrm{‾}}{{\mathit{U}}_i^{*}(t-\mathrm{\Delta }t)}-{\mathit{U}}_i^{*}(t-\mathrm{\Delta }t)}{T_{\mathrm{L}}}}\mathrm{\Delta }t\\ & +\sqrt{{\displaystyle \frac{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}\mathrm{\Delta }t}{T_{\mathrm{L}}}}}{\mathit{\xi }}_i\left(t\right).\end{array}\end{array}$$

Above, $\sqrt{\mathrm{\Delta }t}{\mathit{\xi }}_i\left(t\right)=\mathrm{\Delta }{\mathit{W}}_i\approx \mathrm{d}{\mathit{W}}_i\left(t\right)$ and ξ_i are a two-dimensional vector with independent and identically distributed entries drawn from a standard normal distribution, i.e., ${\mathit{\xi }}_i\left(t\right)\sim \mathcal{N}(\mathbf{0},\mathit{I})$. Finally, we use this updated Lagrangian velocity to update particle positions, again using forward Euler, given by

$$\begin{array}{}\text{(31)}& {\mathit{X}}_i^{*}(t+\mathrm{\Delta }t)={\mathit{X}}_i^{*}\left(t\right)+{\mathit{U}}_i^{*}\left(t\right)\mathrm{\Delta }t.\end{array}$$

3 Large-eddy simulation configuration used to inform particle model

The LES configuration in this study largely follows that of Chun et al. (2023), using version 6.10.9 of the System for Atmospheric Modeling (SAM; Khairoutdinov and Randall, 2003) with additional capabilities for representing the aerosol accumulation mode developed at the University of Washington (UW-SAM; Berner et al., 2013) to simulate an evolving ship plume with varying background aerosol conditions. The LES forcing profiles for the CONTROL simulation are derived from averaging shallow coastal stratocumulus boundary layers (400–800 m deep) in the northeastern Pacific during July 2003 using ECMWF Interim Reanalysis (Zhang et al., 2012; Blossey et al., 2013), and the mean forcings are constant in time for the duration of the simulations. Surface fluxes of latent and sensible heat respond to local conditions according to Monin–Obukhov similarity theory, and radiative transfer calculations were performed every 15 s (every five model time steps) using the Rapid Radiative Transfer Model for GCM applications (RRTMG; Mlawer et al., 1997) with a diurnally varying zenith angle. The evolution of liquid hydrometeors (cloud and rain droplets) was handled by the two-moment Morrison microphysics parameterization (Morrison and Grabowski, 2008), and cloud water was diagnosed via saturation adjustment. The subgrid turbulence was represented using a 1.5-order TKE closure, which allowed for anisotropic diffusivity (Deardoff, 1980). Scalar advection was calculated using the fifth-order ULTIMATE-MACHO scheme (Yamaguchi et al., 2011), and momentum advection was calculated using second-order finite differencing. All simulations were run for 12 h to fully develop turbulence overnight (spin-up period) and for an additional 39.5 h post-ship injection, with the initial injection occurring just before sunrise. The injection time was chosen to maximize the impact of the aerosol–cloud interaction, as sunrise is associated with a diurnal peak in the precipitation rate (Wood, 2012), and injecting during the nighttime hours would result in unchanged shortwave radiative forcing. The LES domain is turned so that the mean wind blows from north to south and is translated along with the mean wind at −10.5 m s⁻¹, which is the strength of the mean wind at 315 m altitude. The ship traverses the domain at the east–west center point (102.4 km) over a 50 min period with a domain-relative speed of 10.5 m s⁻¹ and an aerosol injection rate of 10¹⁶ particles per second.

The vertical grid contains 144 levels with variable grid spacing from 15 m near the surface to 5 m in the cloud layer and in the vicinity of the inversion layer. Grid spacing in the free troposphere increases to near 70 m at the domain top (1.555 km). In this study, our focus is on the zonal (cross-wind) spreading of an injected ship plume, and as such high-aspect-ratio (bowling alley) domain geometries are necessary to maximize the amount of time in which it takes the injected plume to fill the entire domain and to limit east–west plume-edge interactions as a result of periodic lateral boundary conditions. Simulations in Chun et al. (2023) use a 96 km × 9.6 km bowling alley domain with a uniform grid spacing of 50 m. Initial test simulations for the CONTROL run were done on a 102.4 km × 25.6 km horizontal grid (50 m grid spacing), requiring nearly 450 000 core hours, which translates to roughly 9 d of computational time on 2048 processors. After Gaussian curve fitting of the average boundary-layer aerosol concentration, the estimated 2σ ship plume width was approaching the domain size by hour 20 of the simulation, necessitating a wider x dimension to prevent plume-edge overlap. Modeling a 204.8 km × 25.6 km domain at 50 m grid spacing is computationally cost prohibitive and demands a larger horizontal grid spacing (100 or 200 m). Tests of larger grid spacing on the 102.4 km × 25.6 km domain revealed grid spacing sensitivities to coarsening, most notably at 200 m grid spacing (Figs. 1a, 2); however, characteristic mesoscale cell sizes as determined from LWP power spectra show minimal sensitivity to grid spacing (Fig. 1b).

Figure 1 (a) LWP cross-sections at the end of the spin-up period (hour 12) in 100 km × 25 km simulations with 50 m horizontal grid spacing (top), 100 m grid spacing (middle), and 200 m grid spacing (bottom). (b) The smoothed LWP power spectra for the three different horizontal grid spacings. The sloped dashed line represents the $-\mathrm{5}/\mathrm{3}$ power law from satellite observations (Wood and Hartmann, 2006), and the vertical dashed line denotes the characteristic wavelength at the spectral peak (≈ 8 km cell size).

At $\mathrm{\Delta }x=\mathrm{\Delta }y=\mathrm{100}$ and 200 m, the LWP at the end of the spin-up period remains consistent with that for a 50 m grid spacing (Fig. 2a), but coarser grid spacing results in larger boundary-layer aerosol concentrations, weaker surface precipitation, and slightly deeper boundary layers (Fig. 2). Given that large-scale prescribed vertical motion remains the same between runs, an increase in boundary-layer depth in comparison to the CONTROL run indicates higher entrainment rates. In an attempt to reconcile the coarse grid simulations with the 50 m configuration, we apply “hyperdiffusion” to reduce entrainment rates and dampen grid-scale noise (Wyant et al., 2018) in the momentum equations using

$$\begin{array}{}\text{(32)}& {\displaystyle \frac{\mathrm{d}\mathit{U}}{\mathrm{d}t}}=-{\displaystyle \frac{(\mathrm{\Delta }x\mathrm{\Delta }y{)}^{\mathrm{2}}}{{\mathit{\tau }}_{\text{hyper}}}}{\mathrm{\nabla }}^{\mathrm{4}}\mathit{U},\end{array}$$

where τ_hyper is the diffusion timescale and has a value of 1200 and 60 s for $\mathrm{\Delta }x=\mathrm{\Delta }y=\mathrm{100}$ and 200 m, respectively. For the 200 m simulation, strong hyperdiffusion effectively mutes entrainment, resulting in excessive LWP and a shallow boundary-layer depth, all while continuing to underestimate the surface precipitation in comparison to the 50 m run. The inability of the 200 m simulation to match 50 m surface precipitation despite larger LWP and in-line aerosol concentrations suggests that there may be grid sensitivity to precipitation rates independent of entrainment. Using 100 m horizontal grid spacing with τ_hyper=1200 s results in boundary-layer depths, aerosol concentrations, rain rates, and LWP that are in agreement with the 50 m baseline simulation (Fig. 2) and we consequently make use of 100 m grid spacing with hyperdiffusion in order to use sufficiently large domains (204.8 km × 25.6 km).

Figure 2 Spin-up period (first 12 h of the simulation) time series for (a) LWP, (b) boundary-layer-averaged aerosol concentration, (c) surface rain rate, and (d) inversion height. Dashed lines indicate that the simulation included hyperdiffusion.

To explore the sensitivity of the particle model to input parameters from various regions within LES domain, the TKE, variances, and dissipation are conditionally sampled. The ship track (plume) is present if the grid-space-weighted aerosol concentration in the lowest 30 model grid levels is 3 times larger than the maximum column deviation. This stringent criterion ensures that non-ship columns are not incorrectly identified while potentially underestimating the fraction of the grid covered by the ship plume. Figure 3 shows the four conditional statistical categories identified in the LES, with the main conditional average of interest in this study being the cloudy region with ship track present (STcloud or in-plume). Statistics were output every 15 min, and 3-D files were output every 30 min.

Given the previously established spreading rate dependence on both wind shear (Berner et al., 2015) and precipitation (Prabhakaran et al., 2024), our sensitivity studies attempt to span a realizable range of zonal wind shear magnitudes and both precipitating and non-precipitating cases to assess the ability of the particle model to represent a broad range of environmental conditions.

Figure 3 CONTROL run snapshot at 12 h after ship injection (18:00 PDT, Pacific daylight time, hereafter referred to as local time, LT) illustrating the different regions being conditionally sampled in the LES: STcloud is ship track and cloudy, STclear is ship track and not cloudy, NSTcloud is no ship track and cloudy, and NSTclear is no ship track and not cloudy.

3.1 Shear sensitivity tests

Previous modeling studies of stratocumulus have generally focused on the impact of shear across the inversion layer that acts to deplete liquid water through enhanced entrainment and reduce overall TKE (Wang et al., 2008; McMichael et al., 2019; Zapata et al., 2021), while in this study we wish to explore wind shear “within” the boundary layer. The span of zonal shear (cross-plume wind) magnitudes in this study was achieved by altering the geostrophic wind profile to land at a post-spin-up shear that does not dramatically alter vertical shear near the cloud top and restricts the majority of the wind shear to the subcloud layer. In addition to adjusting the forcing wind profile, the weak shear case (WEAK) neglects the apparent Coriolis force. To better constrain zonal shear magnitudes in our simulations, we analyze 2208 Lagrangian trajectories from six different source regions in the northeastern Pacific during June–August in 2018–2021 using ECMWF Reanalysis v5 (ERA5) (Eastman and Wood, 2016; Mohrmann et al., 2019; Erfani et al., 2022). Restricting the analysis to grid levels in ERA5 may result in substantial errors, and we instead linearly interpolate the zonal wind to the estimated boundary-layer and cloud-base heights. The surface wind speed is approximated by the 10 m wind speed. If the cloud base is found to be below the estimated boundary-layer height and the cloud depth exceeds 50 m, we compute the boundary-layer, cloud-layer, and subcloud-layer zonal shear magnitudes as a bulk wind (vector) difference between the top and bottom of the layer. Median maximum boundary-layer zonal shear magnitudes in the ERA5 trajectories are ≈1.7 m s⁻¹, with median subcloud layer shear being almost double that of the cloud layer (Fig. 4).

The CONTROL case maximum zonal shear magnitude in the NSTcloud region is 0.9 m s⁻¹, which is near the 25th percentile of shear magnitudes from ERA5 trajectories, with shear magnitudes of up to 1.5 m s⁻¹ in STcloud during the first evening period when weak decoupling occurs (Fig. 4). For the strong shear case (STRONG), maximum zonal shear magnitudes in NSTcloud are 2.5 m s⁻¹ (75th percentile) with STcloud shear magnitudes exceeding 3.5 m s⁻¹ (> 90th percentile). During the overnight period, the STRONG case zonal shear in STcloud becomes weaker than the non-plume environment and remains weaker for the duration of the run (Fig. 4). The WEAK case maintains less than 0.1 m s⁻¹ of zonal shear, representing an extreme case of low wind shear not seen in the ERA5 trajectories (Fig. 4).

Figure 4 (a) Time series of bulk boundary-layer zonal shear magnitude in the cloudy ship track region (STcloud; dashed lines) and the cloudy non-ship environment (NSTclear; solid lines). The time of 06:00 LT corresponds to the initial injection time. (b) Maximum zonal shear magnitudes estimated from ERA5 Lagrangian trajectories in the northeastern Pacific.

Figure 5 Liquid water path (LWP) evolution of the four large-domain LESs from a bird's-eye view (x–y plane) beginning 1 h after ship injection (07:00 LT) and at 12 h intervals thereafter. The zonal shear vector in CONTROL and STRONG is from right to left.

3.1.1 Macroscale evolution of shear simulations

The varying shear magnitudes result in broadly similar x–y LWP evolution with no conspicuous signal of a ship perturbation 1 h after the initial injection (Fig. 5). All shear simulations exhibit cloud clearing near the ship-plume edge (with the plume edge being easily identifiable in Fig. 6) during the first evening period (hour 13) with the clear region recovering by the following morning (hour 25) (Fig. 5). This cloud-clearing feature arises from a buoyancy anomaly in the ship track region that induces a mesoscale circulation and results in subsidence near the plume edge that creates an inhospitable environment for cloud development (Chun et al., 2023; Wang et al., 2011; Prabhakaran et al., 2024). The CONTROL and STRONG cases initially contain homogeneous, closed-cell convection, which transitions to more broken cloud conditions in the non-ship region, with the ship track region maintaining near-overcast conditions for the first 25 h (Fig. 5). The overcast ship track region begins to break up during the second evening period, with the STRONG and WEAK cases experiencing more severe fragmentation than the CONTROL (Fig. 5). Non-ship-region mesoscale LWP structure is notably different in the WEAK case, with narrow bands of cloud from the first evening onward, indicative of open-cell convection (Fig. 5). Precipitation suppression brought about by the increase in CCN is evident in all shear simulations (Fig. 7). At 1 h after injection, the WEAK case has the most intense precipitating cells (Fig. 7), with a decrease in precipitation intensity later in the day in all cases, which is consistent with the typical diurnal cycle of stratocumulus precipitation (Wood, 2012). Local precipitation enhancement occurs on the down-shear side of the plume edge in the CONTROL and STRONG cases, becoming especially prominent during the second daytime period (Fig. 7). The background aerosol concentration is modified by the precipitation rate through scavenging, with lower environmental aerosol concentrations during the early morning hours, coinciding with the strongest precipitation (Figs. 6, 7). Background aerosol is able to recover during the second evening as entrainment and surface aerosol sources are larger than the aerosol losses from precipitation (Fig. 6).

The CONTROL and STRONG cases have comparable domain-averaged LWP, inversion heights, entrainment rates, surface precipitation rates, and surface fluxes during the first 15 h of the simulations, with the strongest divergence during the first overnight period as the STRONG case experiences less domain-averaged entrainment and an attendant lower average inversion height (Fig. 8). Surface precipitation rates in the WEAK case are more than twice as large as the CONTROL, which depletes LWP, reduces the entrainment rate, and results in a boundary-layer depth that is ≈15 % shallower than the CONTROL. Surface sensible heat fluxes remain similar between the shear cases, although surface latent heat fluxes in the WEAK case are ≈ 20 % smaller than the CONTROL (Fig. 8e, f).

Figure 6 Average dry aerosol number concentration of the bottom 30 grid levels for the four large-domain LESs from a bird's-eye view (x–y plane) beginning 1 h after ship injection (07:00 LT) and at 12 h intervals thereafter.

Figure 7 Surface precipitation rate evolution for the four large-domain LESs from a bird's-eye view (x–y plane) beginning 1 h after ship injection (07:00 LT) and at 12 h intervals thereafter. The zonal shear vector in CONTROL and STRONG is from right to left.

Figure 8 Domain-averaged time series of (a) liquid water path, (b) inversion height, (c) surface precipitation rate, (d) entrainment rate, (e) surface latent heat flux, and (f) surface sensible heat flux. Here, 06:00 LT corresponds to the initial injection time. The gray-shaded region is the time period with no shortwave radiation present.

3.2 Background aerosol sensitivity test

For the CONTROL case, the initial boundary-layer aerosol concentration is 20 # mg⁻¹, while the free-tropospheric aerosol concentration is 50 # mg⁻¹. For the POLLUTED case, the initial boundary-layer aerosol concentration is 130 # mg⁻¹, while the free-tropospheric aerosol concentration is 100 # mg⁻¹. The injection rate in the POLLUTED case is increased to 3.25 × 10¹⁶ s⁻¹ to generate an aerosol perturbation that is roughly the same size as the CONTROL on a percentage basis (Chun et al., 2023). The POLLUTED case large-scale subsidence is increased by 50 % in comparison to the CONTROL as a means to attain a similar boundary-layer depth in the presence of much larger LWP and entrainment rates (Fig. 8).

3.2.1 Macroscale evolution of polluted simulation

The POLLUTED case has a markedly different mesoscale structure than the three shear cases, with cloud fraction near unity and convective cells that remain closed for the duration of the simulation, along with no discernible ship track region in the LWP field (Fig. 5). While background aerosol concentrations in the shear cases decrease during the overnight period as a result of precipitation scavenging, the decrease in average aerosol concentration in the POLLUTED case is a result of the free troposphere being cleaner than the boundary layer and the entrainment acting as a sink instead of a source (Fig. 6). The increased aerosol concentration in the POLLUTED case is sufficient to shut down precipitation production for the duration of the simulation (Fig. 7) and provides an opportunity to explore the ability of the particle model to represent precipitation- or aerosol-concentration-induced responses.

3.3 Langevin particle model input parameters

3.3.1 Domain-averaged input parameters

Zonal variances modestly fluctuate diurnally, with a ≈20 % reduction in zonal variance from the overnight period into the evening period as the buoyancy production of TKE driven by radiative cooling near the cloud top is stunted by solar absorption (Fig. 9a). The meridional variance remains nearly constant for the entire simulation. The diurnal signal is most evident in the vertical velocity variance, with the variance decreasing by nearly 50 % by 16:00 LT of the first day (Fig. 9c). Vertical velocity variance ramps up after 18:00 LT as solar insolation decreases (Fig. 9c). The dissipation rate generally mirrors the TKE, with a minimum in dissipation rate in the evening leading to greater dissipation overnight associated with invigorated turbulence (Fig. 9d, e); however, the relationship between TKE and dissipation can be complex and vary in space and time, depending on local flow geometry, turbulent length scales, and stratification. The TKE divided by dissipation, $k/\mathit{ϵ}$, yields a timescale that is related to the relaxation of a particle's velocity to the mean velocity, and longer relaxation timescales (T_L) can be interpreted as a longer leash on any individual particle, allowing the particle to deviate farther from the mean wind profile before being tugged back strongly. The relaxation timescale for the CONTROL run peaks around 15:00 LT and is smallest and nearly constant during the night (Fig. 9f).

Figure 9 Turbulence-related input parameter time series relevant to the Langevin particle model: (a) zonal (east–west) velocity variance, (b) meridional (north–south) velocity variance, (c) vertical velocity variance, (d) the resolved TKE, (e) the resolved dissipation rate (ϵ), and (f) the relaxation timescale (T_L) without the C₀ correction. All quantities are averaged over the entire boundary-layer depth. Solid lines correspond to domain-averaged time series, while dashed lines correspond to in-plume-averaged (STcloud) time series.

The STRONG case domain-averaged zonal and meridional velocity variance shows little evidence of a diurnal cycle with nearly constant variance for the entire period (Fig. 9a, b). Vertical velocity variance does have a diurnal cycle and is almost indistinguishable from that of the CONTROL simulation (Fig. 9c). The overall TKE remains similar between the two cases, despite the STRONG case having a slightly dampened diurnal cycle (Fig. 9d). From 15:00 LT until sunset on the first day, the nearly equal TKE between the CONTROL and STRONG case coupled with a STRONG case dissipation rate that is marginally larger gives rise to a smaller relaxation timescale in the STRONG case (Fig. 9d, e, f). The WEAK case has an analogous diurnal cycle in zonal variance to the CONTROL, but the magnitude of all components of the variance is lower (Fig. 9a, b, c). The evening restrengthening of vertical velocity variance is less pronounced in the WEAK case, but the relaxation timescale is in line with the CONTROL and STRONG case.

The diurnal cycle of zonal and vertical velocity variance is much larger in the POLLUTED case (Fig. 9a, c), with zonal variances that plummet during the afternoon to values seen in the WEAK case and a rebound in zonal variance overnight that is nearly double that of the CONTROL by the early morning hours of the second day (Fig. 9a). The POLLUTED case relaxation timescale remains lower than all of the shear cases during the first day, with comparable relaxation timescales overnight (Fig. 9f). During the second day, the POLLUTED case relaxation timescale is considerably larger than the three other cases, suggesting a different relationship between TKE and dissipation from the first to the second day (Fig. 9f).

3.3.2 In-plume-averaged input parameters

Now focusing on the conditionally sampled ship track region (STcloud), the CONTROL run experiences enhanced TKE in comparison to the domain average (Fig. 9d). The diurnal cycle of zonal velocity variance is not as apparent in the ship plume, while the diurnal cycle of vertical velocity variance is magnified with large increases in vertical variance overnight (Fig. 9a, c). This is consistent with the findings of Chun et al. (2023), where plume turbulence intensification through the suppression of drizzle dominated over other potential factors, such as increased entrainment efficiency in the ship track region. Although the TKE in the ship track region is larger than the domain-averaged TKE, the dissipation is smaller during the evening of the first day (Fig. 9d, e), suggesting that the energy cascade at small scales is fundamentally different. There are several factors that may contribute to less efficient dissipation and an altered energy cascade. We speculate that the mesoscale circulation that develops in precipitating cases may cause reduced dissipation rates as a result of changing flow geometry, as larger coherent structures within the circulation have smaller internal velocity gradients. The dissipation rate is also dependent on stability, as more stable environments dissipate energy to smaller scales at a slower rate. In comparison to the non-ship track region, the ship track region is more stable and experiences larger negative buoyancy fluxes near the cloud base and stronger daytime decoupling, which also leads to longer ship track relaxation timescales as a result of less efficient dissipation (smaller dissipation rate) given a fixed amount of TKE. Additionally, organized horizontal TKE transport in the ship region may be causing the TKE to be dissipated in the non-ship region. The relative importance of each of these dissipation-modulating mechanisms is currently unknown. The relaxation timescale in the plume is ≈25 % longer than the domain average during the first daytime period but converges to the domain-averaged relaxation timescale during the overnight period, pointing to a return to more isotropic, homogeneous turbulence. The presence of decoupling challenges the well-mixed assumption; however, LES spreading rates in the subcloud and cloud layers differ negligibly from the boundary-layer-averaged spreading rate.

Similar to the CONTROL case, the WEAK and STRONG cases both show stronger turbulence in STcloud in comparison to the domain averages (Fig. 9). The STRONG case has a peak in zonal velocity variance around 16:00 LT, while the WEAK case zonal variance peaks during the overnight hours (Fig. 9a). In-plume relaxation timescales in the STRONG and WEAK cases are lower than their respective domain averages after 18:00 LT of the first day, indicating a more efficient dissipation of energy from then on.

As a consequence of being non-precipitating, the POLLUTED case does not exhibit large differences between in-plume- and domain-averaged turbulent statistics over the first 15 h (Wang et al., 2011; Chun et al., 2023; Prabhakaran et al., 2024). In-plume data points after hour 15 are questionable given that conditionally sampled plume fraction underwent a rapid decrease as the aerosol concentration criteria established before ship injection became too restrictive for the POLLUTED case.

3.4 Plume width calculation for particle model comparison

The particle model lends itself to easily identifying the 1σ plume width simply by computing the standard deviation of the particle position PDF and multiplying by 2 to obtain the full width. The LES plume width is complicated by the presence of the background aerosol variability, and we employ a Gaussian curve fitting procedure on the boundary-layer-average aerosol concentration to define the 1σ plume width. We use the following curve fitting equation

$$\begin{array}{}\text{(33)}& f\left(x\right)=A\exp \left[-{\displaystyle \frac{(x-B{)}^{\mathrm{2}}}{\mathrm{2}{C}^{\mathrm{2}}}}\right]+D,\end{array}$$

where A represents the amplitude of the aerosol perturbation (# mg⁻¹) in the ship plume, x is the position along the x axis (km), B is the location of the ship plume center (km), C is the standard deviation (km), and D is the background aerosol concentration (# mg⁻¹). The algorithm must be provided with initial guesses of $A,B,C,D$, which were $[\mathrm{130.0},\mathrm{102.4},\mathrm{8.0},\mathrm{20.0}]$, respectively, for the CONTROL and shear cases and $[\mathrm{350.0},\mathrm{102.4},\mathrm{8.0},\mathrm{130.0}]$, respectively, for the POLLUTED case. Gaussian curve fitting along the x dimension was performed at each y location and then averaged along the entire y dimension to create a single 1σ plume width estimate, which is then multiplied by 2 to get an actual plume width.

4 Results

4.1 Large-eddy simulation plume width results

For an initial delta function injection and purely diffusive plume spreading, we define the time-varying width of the plume as

$$\begin{array}{}\text{(34)}& \mathcal{W}\left(t\right)=\mathrm{2}\sqrt{\mathrm{2}\mathcal{D}t},\end{array}$$

where 𝒟 is the coefficient for eddy diffusivity and t is time.

When analyzing the eddy diffusivity output from the 1.5-order TKE closure in the LES, maximum subgrid diffusivity values are near 0.75 m² s⁻¹, which results in a plume width of 500 m at hour 10 after injection, while the CONTROL 1σ plume width at hour 10 is 24 km (Fig. 10). It is then immediately obvious that the spreading of the aerosol is not solely related to mixing done by the smallest scales. As an alternative calculation, 𝒟 can be estimated as the product of a characteristic zonal eddy velocity (u_c) and a characteristic length scale (l) on which the transport occurs (𝒟=u_cl). In the CONTROL run, the characteristic zonal eddy velocity is ≈ 0.3 m s⁻¹, while the characteristic length scale may potentially be driven by mesoscale eddies that are initially ≈8 km wide. The resulting constant Gaussian diffusion equation estimates a width of 26.3 km at hour 10, which is much closer to the LES plume width (Fig. 10). While assuming the larger scales are doing a majority of the transport related to the ship track spreading does improve the simple Gaussian diffusion model performance, it fails to capture the accelerating growth (growth rates grow from 2 km h⁻¹ in the first 9 h to near 3 km h⁻¹ during hours 9–13) in the first 15 h after injection in the CONTROL run (Fig. 10). Given the inaccuracies of the Gaussian diffusion assumption in the first 15 h and uncertainties regarding the characteristic velocity and time-dependent characteristic length scales, an approach tied to the turbulence properties in the LES is desired and motivates the model developed in Sect. 2.

Figure 10 The single standard deviation (multiplied by 2 to get plume width) width estimated from the Gaussian curve fitting procedure applied to the LES aerosol concentrations at 30 min intervals. Hour 0 corresponds to 06:00 LT. The dotted–dashed line represents a constant plume growth rate of 2 km h⁻¹ (Heffter, 1965; Durkee et al., 2000). The dashed line represents a Gaussian diffusion curve (Eq. 34) where u_c=0.3 m s⁻¹ and l=8 km. The shaded region is the nighttime period.

The CONTROL run plume grows at 2–3 km h⁻¹ during the first 13 h before decreasing around sunset and remaining near 1 km h⁻¹ during the night. Further reductions in plume growth (< 1 km h⁻¹) during the second day lead to a final plume width of 48.8 km (Fig. 10). The STRONG case plume growth remains similar to the CONTROL for the first 5 h after injection before a burst of 3 km h⁻¹ growth (Fig. 10), and by hour 10 the STRONG case plume is > 3 km wider than the CONTROL. Similarly to the CONTROL, the STRONG case plume growth rate lessens during the first evening and remains nearly constant at 1 km h⁻¹ throughout the night and into the second day, resulting in a final plume width of 60.8 km (Fig. 10). WEAK case growth rates are initially slower than the CONTROL case with a plume width of 18 km at hour 10 (6 km narrower than CONTROL), but overnight spreading rates are consistent with the STRONG and CONTROL cases (Fig. 10). The WEAK case plume spreading during the second day persists at ≈ 1 km h⁻¹, and the final plume width is 44 km. The POLLUTED case experiences rapid growth in the first few hours after injection, but average spreading rates between hours 5–15 are slower in comparison to the three other cases (Fig. 10) and consistent with previous modeling studies indicating slower growth in non-precipitating boundary layers (Prabhakaran et al., 2024). Again, overnight and second-day spreading rates for the POLLUTED case are ≈ 1 km h⁻¹, suggesting overnight horizontal spreading rates may be independent of zonal variance intensity. In the sensitivity cases examined here, the difference in LES plume width at hour 15 between the STRONG and POLLUTED cases is 14 km, emphasizing the importance of accounting for different background conditions when estimating plume growth rates.

4.2 Particle model results for sheared cases

The particle model requires initial conditions for both particle position and particle velocity, with the standard deviation of particle position being initialized as half of the first available plume width in the LES and the standard deviation of particle velocity being set to 0. Additionally, the particle model contains a free parameter, C₀, that is found by spanning a range of C₀ values using 20 000 particle simulations with 2 min time steps and determining the value associated with the minimum cumulative least-squares error compared to the LES. Since we have two different input parameter categories (domain average and in-plume) we find two distinct C₀ values for the CONTROL, which are C₀=0.37 for domain-averaged statistics and C₀=0.69 for the in-plume statistics. These constants are then applied to the sensitivity cases in hopes of being physically representative of broad-ranging turbulence behavior. By running numerous ensemble members with a low particle number, the number of total particles modeled can be reduced to near 1000 (50 ensemble members, 20 particles each) with minimal deterioration in model performance and consistency. The results shown here use 50 ensemble members of 100 particles each, with each separate realization of the ensemble mean being stable and in agreement with simulations with 20 000 particles (or greater).

Figure 11 Ensemble particle model results for the CONTROL run using the original relaxation timescale formulation (T_L) laid out in Eq. (18). (a) Particle model results using the domain-averaged input parameters, (b) particle model results using the in-plume (STcloud) averages, and (c) the error between the particle model and the LES plume width for domain-averaged and in-plume input parameters. Individual dashed lines represent each ensemble member.

The domain-averaged input parameters are able to largely capture the plume spreading throughout the entire simulation, with errors at any point not exceeding 2.5 km (Fig. 11). While the particle model broadly captures the spread, subtleties such as the downward inflection in the spreading rate near hour 13 are only captured when the model is forced with turbulence data from within the ship plume (Fig. 11b). Neither the domain averages nor in-plume averages capture the second-day slowing of plume growth shortly after sunrise (near hour 25) (Fig. 11). When the particle model is applied to the two shear sensitivity cases, the performance becomes case dependent. Both the domain-averaged and in-plume particle model forcings are able to time the deceleration in spreading rate adeptly; however, plume widths in the STRONG case are underestimated substantially (Fig. 12).

Figure 12 Ensemble particle model results (only ensemble averages are shown) from the particle model are dashed lines, while solid lines are LES plume widths. (a) Particle model forced with domain-averaged statistics and (b) particle model forced with in-plume (STcloud) statistics. C₀ is 0.37 and 0.69 for the domain-averaged and in-plume runs, respectively. The original relaxation timescale formulation (T_L) was used.

It is desirable for the particle model to capture the clear divergence in LES spreading rates between the STRONG case and CONTROL between hours 7–10 after injection. Zonal variance in the STRONG case is larger during this time frame, and on its own it would result in faster spreading; however, the STRONG case relaxation timescale is equal to or less than that of the CONTROL, leading to particle model spreading rates that are nearly identical. The original formulation of the relaxation timescale assumes isotropic turbulence, but we are instead interested in the variance along the dimension of the spread, which in this case is the zonal variance. We introduce a modified relaxation timescale calculation (T_m) that not only considers the spread dimension variance but also maintains an assumption of isotropic dissipation

$$\begin{array}{}\text{(35)}& T_{\mathrm{m}}={\displaystyle \frac{\frac{\mathrm{1}}{\mathrm{2}}{u^{\prime }}^{\mathrm{2}}}{\frac{\mathrm{3}}{\mathrm{4}}{C}_{\mathrm{m}}\mathit{ϵ}}},\end{array}$$

where ${u^{\prime }}^{\mathrm{2}}$ is the zonal variance and C_m is a new constant that must be found through another round of CONTROL run optimization. The C_m constant is 0.15 for the domain-averaged forcing and 0.29 for the in-plume forcing. The new C_m values account for the magnitude of T_m being smaller in comparison to T_L, but the diurnal cycle of T_m is more amplified in the STRONG case (Fig. 13).

Figure 13 (a) The original relaxation timescale formulation (T_L) and (b) the modified version of the relaxation timescale (T_m), which focuses only on the zonal variance. Solid lines correspond to domain-averaged time series, while dotted–dashed lines correspond to the in-plume-averaged (STcloud) time series.

Figure 14 Ensemble particle model results (only ensemble averages are shown) from the particle model using the modified relaxation timescale (T_m). Particle model results are dashed lines, while solid lines are LES plume widths. (a) Particle model forced with domain-averaged statistics. (b) Particle model forced with in-plume (STcloud) statistics. C_m is 0.15 and 0.29 for the domain-averaged and in-plume runs, respectively.

Figure 15 Ensemble average particle model result for the POLLUTED case. Particle model results are dashed lines, and the solid line is the LES plume width. PM (C₀=0.37) is using the original relaxation timescale formulation T_L and C₀.The PM of modified TS (C_m=0.15) uses T_m and the CONTROL-optimized C_m. The PM of modified TS (C_np=0.50) uses T_m and a new optimized constant C_np that applies to non-precipitating cases.

The use of T_m resolves the relatively slow spread in the STRONG case during the 7–10 h period, with final plume widths in line with the LES when using domain-averaged forcings. Changing T_m does not appreciably change the CONTROL or WEAK case results using domain-average or in-plume statistics (Fig. 14). Using T_m with in-plume forcing results in good performance through hour 15, followed by poor performance during the nighttime and second-day periods (Fig. 14b). By applying T_m and C_m=0.15, the CONTROL, STRONG, and WEAK case particle model results all agree well with LES plume widths, with changes in diurnal-cycle-related spreading rates being captured using only domain-averaged turbulent statistics.

Figure 16 The x–z cross-sections of zonal velocity with pink going left to right and green going right to left. The CONTROL run containing a strong mesoscale circulation in the vicinity of the ship is on the left. The POLLUTED case containing no mesoscale circulation is shown on the right.

4.3 Particle model results for the POLLUTED case

The POLLUTED case forced with domain-averaged statistics using both T_L and T_m performs poorly in comparison to LES plume width, with growth that is much too extreme in the first 10 h of the particle simulation (Fig. 15), signifying that C₀ and C_m are not translating as “universal” constants to the non-precipitating POLLUTED case. All precipitating cases (CONTROL, STRONG, WEAK) develop mesoscale circulations shortly after ship injection that intensify until sunset, when stronger boundary-layer turbulence interferes with the circulation (Fig. 16). These mesoscale circulations are believed to arise from the suppression of precipitation and the associated buoyancy anomaly in the ship track region (Prabhakaran et al., 2024; Chun et al., 2023; Wang et al., 2011), but what mechanisms sustain, intensify, or destroy them is currently unknown. Even in the absence of a complete understanding of mesoscale circulation dynamics, it remains a critical deviation from homogeneous, isotropic turbulence for which the particle model was developed. In this sense, C₀ and C_m are specifically optimized to represent cases with mesoscale circulations present. Increased spreading rates are achieved through a reduction in isotropic dissipation in Eq. (35) through multiplication with (C₀,C_m), allowing for more rapid spread than would otherwise be possible given the domain-averaged forcings. Accounting for anisotropy with the Langevin model is commonly done by adding a correction term dependent on spatial derivatives of the variance (Legg and Raupach, 1982; Dehbi, 2008); however, such a correction is not desirable given that the purpose of this simplified model is to be used as a subgrid parameterization wherein variance gradients are not accessible. Therefore, we optimize C₀ for the POLLUTED case separately and propose the use of two constants for domain-averaged statistics using T_m: one for precipitating cases where mesoscale circulations develop (C_p=0.15) and one for non-precipitating cases where turbulence remains nearly isotropic (C_np=0.50). By correcting the POLLUTED case to not implicitly represent a mesoscale circulation, the particle model is aligned with the LES width.

5 Conclusions

Estimating subgrid plume fraction in a climate model grid box necessitates a method of approximating the rate at which the plumes grow. Neglecting such growth or assuming constant or Gaussian-diffusion-based growth may cause nonlinear errors that lead to unreliable or unrealistic responses to MCB. We employ a Lagrangian particle model, driven by large-domain LES output of ship tracks, and assess the reduced-order model's ability to represent plume spreading in environments with different shear magnitudes and environments with or without precipitation.

Using only 5000 particles (50 ensemble members with 100 particles) in the naturally parallel stochastic particle model, both the domain average and conditionally sampled plume TKE, variances, and dissipation rates result in good agreement with the LES CONTROL plume width. Extending the CONTROL-optimized free parameter C₀ to the STRONG case resulted in poor performance because the original relaxation timescale formulation T_L assumes that each component of the variance is equally contributing to the zonal spread. Instead, we apply a modified relaxation timescale T_m that only focuses on the variance in the direction of the spread (zonal direction in this study). This modified formulation results in appropriately larger spreading rates in the STRONG case and minimally impacts the CONTROL and WEAK case spreading rates. The particle model formulation that performs best in comparison to the LES shear sensitivity case widths uses domain-averaged input parameters; the modified relaxation timescale T_m; and a re-optimized C₀, named C_m, which accounts for the new relaxation timescale formulation.

In-plume turbulent statistics perform better than domain-averaged quantities during the first 15 h after injection, but as nocturnal turbulence disrupts the mesoscale circulation, the daytime relationship between the plume-optimized turbulence constant (C_m) and the dissipation rate breaks down and results in larger errors thereafter. As the sun sets, the domain-averaged statistics continue to represent spreading rates well during the night and into the second day, potentially as a result of the domain-averaged C_m being less sensitive to the termination of the mesoscale circulation.

Applying the best-performing particle model conditions (domain-averaged input parameters, T_m, C_m) to the POLLUTED case generates excessive first-day spreading rates as the CONTROL case that C_m was optimized for contained a large mesoscale circulation that aided in plume dispersion and no mesoscale circulation exists in the POLLUTED case. Using C_m in the POLLUTED case artificially decreases dissipation, although no such anisotropy is present. When C_m is optimized for the POLLUTED non-precipitating environment C_np, we find that the particle model is again able to accurately recreate the spreading rate geometry in the LES.

The particle model is able to capture the impacts of an anisotropic, spread-accelerating mesoscale feature in precipitating cases using only domain-averaged input parameters through a change in the free parameter C₀. Using C_m=0.15 for precipitating environments that behave anisotropically during the daytime and C_np=0.50 for non-precipitating environments is one potential way of easily dealing with anisotropic drift without the addition of an extra term that would require information not available to the subgrid parameterization, such as the spatial gradient of variance. The Langevin particle model is able to represent spreading rates better than traditional methods of constant Gaussian diffusion, all while using domain-averaged turbulence statistics. Considering that domain-averaged information from the turbulence parameterization within a climate grid box is the only available input into a would-be subgrid particle model, the particle model performance in this study is promising, suggesting such an approach may be a viable method of cheaply and accurately modeling plume spreading in different environments. Having access to subgrid horizontal variances to calculate T_m within a GCM grid box would require a higher-order closure such as CLUBB (Larson and Golaz, 2005; Guo et al., 2015) or perhaps an additional parameterization based on resolved wind shear and boundary-layer-integrated radiative cooling.

While we have examined a broad range of shear magnitudes in this paper, only a small portion of the variance and relaxation timescale parameter space has been explored. Deeper boundary-layer cases associated with stratocumulus-to-cumulus transitions can have accelerating spreading rates overnight (Prabhakaran et al., 2024) with spreading rates in excess of 4 km h⁻¹. Future plans are to apply the particle model to a much broader range of conditions. It is also worth noting that with a sufficiently large library of high-resolution LES, it may be possible to machine learn a time-dependent 𝒟 in Eq. (34); however, the computational cost of multi-day ship track simulations currently precludes the creation of such a library.