Deformation of ambient chemical gradients by sinking spheres

A sphere sinking through a chemical gradient drags fluid with it, deforming the gradient. The sphere leaves a trail of gradient enhancement that persists longer than the velocity disturbance in the Reynolds $10^{-2}\leqslant Re\leqslant 10^{2}$, Froude $10^{-1}\leqslant Fr\leqslant 10^{3}$ and Péclet $10^{2}<Pe\leqslant 10^{6}$ regime considered here. We quantify the enhancement of the gradient and the diffusive flux in the trail of disturbed chemical left by the passing sphere using a combination of numerical simulations and scaling analyses. When $Fr$ is large and buoyancy forces are negligible, dragged isosurfaces of chemical form a boundary layer of thickness $\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D70C}}$ around the sphere with diameter $l$. We derive the scaling $\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D70C}}/l\sim \mathit{Pe}^{-1/3}$ from the balance of advection and diffusion in the chemical boundary layer. The sphere displaces a single isosurface of chemical a maximum distance $\mathit{L}_{Def}$ that increases as $\mathit{L}_{Def}/l\sim l/\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D70C}}\sim \mathit{Pe}^{1/3}$. Increased flux through the chemical boundary layer moving with the sphere is described by a Sherwood number, $Sh\sim l/\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D70C}}\sim \mathit{Pe}^{1/3}$. The gradient enhancement trail extends much farther than $\mathit{L}_{Def}$ as displaced isosurfaces slowly return to their original positions through diffusion. In the reference frame of a chemical isosurface moving past the sphere, a new quantity describing the Lagrangian flux is found to scale as $\mathit{M}\sim (\mathit{L}_{Def}/l)^{2}\sim \mathit{Pe}^{2/3}$. The greater $\mathit{Pe}$ dependence of $\mathit{M}$ versus $Sh$ demonstrates the importance of the deformation trail for determining the total flux of chemical in the system. For $\mathit{Fr}\geqslant 10$, buoyancy forces are weak compared to the motion of the sphere and the preceding results are retained. Below $\mathit{Fr}=10$, an additional Froude dependence is found and $l/\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D70C}}\sim Sh\sim Re^{1/6}Fr^{-1/6}Pe^{1/3}$. Buoyancy forces suppress gradient deformation downstream, resulting in $\mathit{L}_{Def}/l\sim Re^{-1/3}Fr^{1/3}Pe^{1/3}$ and $\mathit{M}\sim Re^{-1/3}Fr^{1/3}Pe^{2/3}$. The productivity of marine plankton – and therefore global carbon and oxygen cycles – depends on the availability of microscale gradients of chemicals. Because most plankton exist in the fluids regime under consideration, this work describes a new mechanism by which sinking particles and plankton can stir weak ambient chemical gradients a distance $\mathit{L}_{Def}$ and increase chemical flux in the trail by a factor of $\mathit{M}$.

A sphere sinking through a chemical gradient drags fluid with it, deforming the gradient. The sphere leaves a trail of gradient enhancement that persists longer than the velocity disturbance in the Reynolds 10 −2 Re 10 2 , Froude 10 −1 Fr 10 3 and Péclet 10 2 < Pe 10 6 regime considered here. We quantify the enhancement of the gradient and the diffusive flux in the trail of disturbed chemical left by the passing sphere using a combination of numerical simulations and scaling analyses. When Fr is large and buoyancy forces are negligible, dragged isosurfaces of chemical form a boundary layer of thickness δ ρ around the sphere with diameter l. We derive the scaling δ ρ /l ∼ Pe −1/3 from the balance of advection and diffusion in the chemical boundary layer. The sphere displaces a single isosurface of chemical a maximum distance L Def that increases as L Def /l ∼ l/δ ρ ∼ Pe 1/3 . Increased flux through the chemical boundary layer moving with the sphere is described by a Sherwood number, Sh ∼ l/δ ρ ∼ Pe 1/3 . The gradient enhancement trail extends much farther than L Def as displaced isosurfaces slowly return to their original positions through diffusion. In the reference frame of a chemical isosurface moving past the sphere, a new quantity describing the Lagrangian flux is found to scale as M ∼ (L Def /l) 2 ∼ Pe 2/3 . The greater Pe dependence of M versus Sh demonstrates the importance of the deformation trail for determining the total flux of chemical in the system. For Fr 10, buoyancy forces are weak compared to the motion of the sphere and the preceding results are retained. Below Fr = 10, an additional Froude dependence is found and l/δ ρ ∼ Sh ∼ Re 1/6 Fr −1/6 Pe 1/3 . Buoyancy forces suppress gradient deformation downstream, resulting in L Def /l ∼ Re −1/3 Fr 1/3 Pe 1/3 and M ∼ Re −1/3 Fr 1/3 Pe 2/3 . The productivity of marine plankton -and therefore global carbon and oxygen cycles -depends on the availability of microscale gradients of chemicals. Because most plankton exist in the fluids regime under consideration, this work describes a new mechanism by which sinking particles and plankton can stir weak ambient chemical gradients a distance L Def and increase chemical flux in the trail by a factor of M.

Introduction
A spherical particle sinking at intermediate Reynolds number through a chemical gradient draws fluid downward, disturbing the local composition field without generating turbulence. If the diffusivity of the chemical is small compared to the particle's size and speed, the particle leaves behind a trail in which the chemical gradient is intensified in some regions and weakened in others. This gradient deformation is evident in figure 1, modified from Yick et al. (2009), which shows a microscale schlieren image of a sphere sinking at low Reynolds number in a salinity-stratified tank. The schlieren technique visualizes changes in the background linear density gradient; here the gradient is enhanced around the sphere and in two tracks that extend 20 sphere diameters downstream. Variation in the gradient across chemical isosurfaces must also result in changes in the diffusive flux in the vicinity of the trail. Yick et al. (2009) focused on the increased drag on a sphere due to stratification and did not characterize gradients in the trail or the diffusive flux.
A number of studies have depicted disturbance of stratification by falling spheres in order to describe increased drag (Eames & Hunt 1997;Srdić-Mitrović, Mohamed & Fernando 1999;Yick et al. 2009;Zhang, Mercier & Magnaudet 2019;Magnaudet & Mercier 2020), settling speed anomalies (Abaid et al. 2004;Camassa et al. 2009Camassa et al. , 2010Camassa et al. , 2013Doostmohammadi, Dabiri & Ardekani 2014), internal wave generation (Mowbray & Rarity 1967) and buoyant jets behind the sphere (Ochoa & Van Woert 1977;Torres et al. 2000;Hanazaki, Kashimoto & Okamura 2009a;Hanazaki, Konishi & Okamura 2009b;Hanazaki, Nakamura & Yoshikawa 2015). These studies indicate that perturbations of the compositional field by a sphere depend on viscosity, diffusivity and buoyancy, which we parameterize by the Reynolds number Re = Wl/ν, the Péclet number Pe = ReSc = Wl/κ and the Froude number Fr = W/Na, where W is the velocity of the sphere, l is the sphere diameter, a is the sphere radius, ν is the kinematic viscosity, κ is the diffusivity, N is the Brunt-Väisälä frequency and Sc = ν/κ is the Schmidt number. The purpose of this work is to quantify gradient deformation and diffusive flux enhancement by a sphere descending through a general chemical gradient in terms of the Reynolds, Froude and Péclet numbers. At low Fr, disturbance of the compositional gradient generates buoyancy forces, but when Fr is sufficiently large, buoyancy forces are negligible and the chemical behaves as a passive, diffusive substance advected by the sphere. Both cases are examined in this work.
Alteration of diffusive fluxes by gradient deformation is relevant to any system with objects moving through inhomogeneities at low to intermediate Reynolds numbers, such as particles sinking in the ocean and atmosphere. Sinking detritus and swimming plankton drag fluid and may generate small-scale fluxes of nutrients and other chemicals important for microscale ocean ecology (Katija & Dabiri 2009;Noss & Lorke 2014;Wang & Ardekani 2015;Houghton & Dabiri 2019), although there is some debate as to whether plankton mix vertical stratification at large scales (Dewar et al. 2006;Visser 2007;Houghton et al. 2018). Dust and bacteria are ubiquitous in the atmosphere and act as a substrate for ice nucleation and cloud formation (DeMott et al. 2003;Bowers et al. 2009). Microscale stirring by these particles could influence nucleation rates and factor into the ecology of airborne bacteria (Grover & Pruppacher 1985;Lighthart 1997). Another potential application concerns the precipitation of heavy solids at the top of terrestrial planetary cores. Falling iron generates convective motions that have been linked to the generation of the ancient magnetic fields of the Moon (Laneuville et al. 2014) and Mars (Davies & Pommier 2018) and the present magnetic field of Ganymede (Rückriemen, Breuer FIGURE 1. Microscale schlieren image of a sphere sinking through salinity stratification for Re ≈ 2, Fr ≈ 10 and Pe ≈ 10 3 (reproduced from Yick et al. (2009)). The greyscale intensity indicates the absolute value of the magnitude of the gradient of the salinity anomaly. Superimposed curves are isosurfaces of salinity from our simulation using equivalent parameters. The scale bar is 3 mm.
& Spohn 2015). The influence of this 'iron snow' on background stratification in the planetary core is not known but may affect the nature and vigour of convection.
The factors that characterize the extent of gradient deformation are the thickness of the chemical boundary layer δ ρ , the distance over which chemical isosurfaces are dragged L Def , the chemical flux through the boundary layer Sh and a new quantity M describing the Lagrangian flux.
Deformation of a plane of dye by a sphere was first described by Darwin (1953), who derived the position of dye particles in potential flow after the sphere travelled an infinite distance. Lighthill (1956) found asymptotic solutions for the location of the dye at finite times during the passage of the sphere. For a sphere moving through a linear gradient of dye, the positions of a single dye isosurface at time increments t are equivalent to isosurfaces of the dye field at concentration increments ρ, in the frame of reference of the sphere. Eames & Hunt (1997) used this equivalence to calculate the composition and flow fields for a sphere moving in weak stratification at high Re. In the absence of diffusion, isosurfaces pile up at the stagnation point upstream of the sphere as they are stretched infinitely far. The composition and its gradient are then singular on the surface of the sphere, which they predicted would generate a downstream jet.
Numerical simulations conducted by Torres et al. (2000) for Re > 25 demonstrated that stratification collapses the standing vortices that normally appear in the wake of a sphere, and the displaced isosurfaces generate a buoyant jet on their return. In addition to confirming the prediction from Eames & Hunt (1997), they showed that diffusion allows the sphere to pass through isosurfaces so that they are only dragged a finite distance, allowing steady solutions to be found. The maximum vertical distance an isosurface can be dragged, L Def , is correlated with the extent of gradient 892 A33-4 B. G. Inman, C. J. Davies, C. R. Torres and P. J. S. Franks deformation, as we will show later. Hanazaki et al. (2015) employed a dimensional analysis of the buoyancy time scale to estimate this distance as lπFr for Re 200, 1 Fr 10 and Pe > 10 5 . Yick et al. (2009) combined experiments and numerical results to empirically estimate the dragged distance as l √ Fr for Re < 1, Fr < 1 and Pe < 10 3 . Both they and Hanazaki et al. (2009b) suggested that L Def should increase with decreasing diffusivity, but neither provided Sc or Pe scalings.
Compositional isosurfaces that are dragged downward by the sphere constitute a chemical boundary layer near the surface that is analogous to the classic problem of mass transfer from a sphere. For Re → ∞ and Sc → 0, the thickness of the chemical boundary layer is described by δ ρ /l ∼ Pe −1/2 (Schlichting 1968;Yih 1969). Previous studies of spheres sinking through salinity gradients have used this relation to justify the grid spacing in their numerical models but noted that it underestimates the chemical boundary layer thickness for finite Re and Sc (Torres et al. 2000;Hanazaki et al. 2009aHanazaki et al. , 2015Doostmohammadi et al. 2014). In this work, we derive a new scaling for δ ρ in the 10 −2 Re 10 2 and Pe > 10 regime. We show that this scaling is critical for determining L Def , Sh and M.
Here we present a combined analytical and numerical study of gradient deformation and diffusive flux enhancement by a sphere descending through a linear chemical gradient. Methods are detailed in § 2. In § 3.1, we describe the general process of gradient deformation with a representative case and define the quantities δ ρ , L Def , Sh and M that are used in the subsequent analysis to characterize the influence of diffusion ( § 3.2), viscosity ( § 3.3) and buoyancy ( § 3.4). Previous work on spheres sinking through salinity gradients has focused on Re 200 or Re < 10, Fr < 10 and Sc = 7, 700 (Torres et al. 2000;Yick et al. 2009;Hanazaki et al. 2015). The parameter range considered here is 10 −2 Re 10 2 , 10 −1 Fr 10 3 and 10 2 < Pe 10 6 . This regime is relevant to particles and plankton in the ocean. The potential for sinking ocean particles to deform ambient chemical gradients is discussed in § 4.

Governing equations
We consider a sphere sinking through a linear, vertical gradient of a diffusive chemical, such as the nutrient nitrate in the upper ocean. At low Froude numbers, disturbance of the compositional gradient induces buoyancy forces. When Fr 10 3 , buoyancy forces are negligible and the sphere can then be envisioned as travelling in the direction of a gradient of passive, dilute chemical. Both cases are examined in this work and explicit reference to the Froude number is made to distinguish them. In the stationary reference frame of a sphere falling with speed W, the dimensional forms of the governing equations are and where ρ is the density of the fluid, t is time, u is the velocity vector, p is the pressure, g is the acceleration due to gravity, µ is the dynamic viscosity of the fluid and κ is the diffusivity coefficient of the chemical. The Boussinesq approximation is used and it is assumed that variation in ρ is due to differences in chemical concentration alone. Therefore spatial gradients of ρ are referred to as chemical gradients. The flow is assumed to be axisymmetric and soẑ is the unit vector in the direction of the vertical axis of symmetry and u = (u, w), with u the radial and w the vertical components of velocity in cylindrical coordinates. Laboratory experiments with spheres descending through both stratified and unstratified fluids have exhibited axisymmetric flow patterns when Re < 200 (Taneda 1956;Ochoa & Van Woert 1977;Yick et al. 2009;Hanazaki et al. 2009a;Okino, Akiyama & Hanazaki 2017). Numerical simulations are conducted in the reference frame of the sphere such that isosurfaces of chemical appear to move upward past the sphere at constant speed W. The background, unperturbed chemical field and hydrostatic pressure are functions of z and t:ρ where ∂ρ/∂z is the undisturbed vertical gradient of chemical and ρ 0 is a reference value of the total fluid density at z = 0 and t = 0. The total density and total pressure are decomposed into the background states and respective perturbations ρ and p : (2.6) and p =p(z, t) + p (r, z, t). (2.7) Substituting the ρ and p decompositions in (2.6) and (2.7) into (2.1), subtracting the background state and applying the Boussinesq approximation produce where ν = µ/ρ 0 is the kinematic viscosity. A similar substitution of (2.6) into (2.2) gives Dρ Dt = − ∂ρ ∂z (w − W) + κ∇ 2 ρ . (2.9) The non-dimensional forms of (2.8) and (2.9) are obtained by scaling distances by the diameter of the sphere l = 2a, velocities by the background flow W, pressure perturbation by ρ 0 W 2 and chemical perturbation by −l(∂ρ/∂z) as follows: (2.11) The non-dimensional parameters in (2.10) and (2.11) are the Froude number Fr = W/Na, Reynolds number Re = Wl/ν and Péclet number Pe = Wl/κ, where N is the Brunt-Väisälä frequency defined by When Fr is small, ρ enters the momentum equation (2.10) through the buoyancy term, and disturbance of the chemical gradient generates buoyancy forces. When Fr 10 3 , the buoyancy term is negligible compared to the other terms for the parameter range studied and ρ no longer has an influence on the momentum equation (i.e. (2.10) and (2.11) are decoupled). This is confirmed by comparing Fr = 10 3 simulations to test cases where the buoyancy term is set to zero. Therefore, ρ signifies a passive, dilute substance at Fr = 10 3 . In this case, ρ diffuses away over time through the last term in (2.11) because w − 1 → 0 away from the sphere and ρ = 0 is the background state (2.6) which is enforced on the outer boundary as described below. In other words, diffusion allows the chemical field to assume its background state following a perturbation, in the absence of buoyancy forces.
Because (2.3), (2.10) and (2.11) do not explicitly include the time development of pressure, we employ a diagnostic Poisson equation for pressure. Taking the divergence of (2.10) gives where D = ∇ · u. The time discretization of the diagnostic pressure equation satisfies the divergence-free requirement in (2.3), as described in the next section. Boundary conditions on the surface of the sphere (z 2 + r 2 = a 2 ) are no-slip and no chemical flux: u = (u, w) = (0, 0) (2.13) and ∇ρ ·n = 0, (2.14) wheren is the unit vector normal to the sphere surface. The non-dimensional form of the no-flux boundary condition (2.14) is  Torres et al. (2000) and Yick et al. (2009). Solutions for the pressure, velocities and chemical concentration are calculated each time step on a computational grid using finite differences with the marker and cell (MAC) method (Harlow & Welch 1965). The MAC method has been modified for incompressible stratified fluids (Hanazaki 1988) and density diffusion (Torres et al. 2000), and is adjusted here for Re < 200 and chemical gradients with diminishing buoyancy forces.
The time evolution of the flow is discretized using an explicit Euler method. Dropping the primes that denote non-dimensional variables for clarity in this section only, equations (2.10)-(2.12) are rewritten as where n is the current time step corresponding to an integration time of t = n t. The initial condition is an impulsive start in which ρ = 0, u = 0 and w = W everywhere except at the boundary of the sphere where w = 0 as per (2.13). At time step n, the Poisson equation for pressure (2.23) is solved with given values for u n and ρ n using the successive over-relaxation method. In accordance with the MAC method, D n+1 is set equal to zero to preserve the incompressibility condition (2.3) and D n is retained to prevent the accumulation of numerical round-off errors. Inserting p n into (2.21), u and ρ are then calculated for the next time point n + 1. This process is repeated for each subsequent time step until a steady-state criterion is reached (see § 2.4).
2.3. Curvilinear grid Equations (2.21)-(2.23) are transformed from the physical plane (r, z) -represented by the curvilinear grid in figure 2 -to a rectilinear computational domain (ξ , η) as described in Torres et al. (2000): A novel curvilinear grid was developed for this work in order to resolve the chemical gradient both near the sphere as well as hundreds of sphere diameters downstream (figure 2). Curvilinear grid rays and arcs correspond to computational coordinates ξ and η, respectively. As in Torres et al. (2000) and Hanazaki et al. (2015), grid arcs (constant η) are clustered towards the sphere to better resolve the thin chemical boundary layer. Using Pe −1/2 to estimate the thickness of the boundary layer (Schlichting 1968;Yih 1969), the number of grid points in this layer is three for Pe = 10 6 . The Pe −1/2 estimate is overly conservative for high Péclet numbers as will be demonstrated in § 3.2. This is comparable to the grid used in Hanazaki et al. (2015) for Re = 200 and Sc = 700. Near the sphere, grid rays (constant ξ ) are normal to its surface to accurately implement the no-flux and no-slip boundary conditions on . The curvilinear grid used in this study. Every fifth arc and ray are drawn, for clarity. Arcs are clustered towards the sphere to resolve the boundary layer. Rays are clustered towards the downstream axis of symmetry, the positive z axis. Rays are perpendicular close to the surface of the sphere, and then curve slightly towards the z axis to preserve resolution at great distances downstream.
a curve. Rays are clustered towards the forward (upstream) end of the sphere where chemical isosurfaces collect before being punctured by the sphere. As in Hanazaki et al. (2015), rays are also clustered towards the rear (downstream) stagnation point to improve resolution of wake structures.
In the large-Péclet-number regime of interest for this study, variation of chemical perturbations is primarily transverse to and concentrated near the z axis far downstream of the sphere (figure 1). The goal, then, in generating a grid is to maintain high ξ resolution transverse to the downstream z axis over a very large domain while remaining computationally tractable. Straight rays cause a loss in resolution near the axis of symmetry with increasing z; instead, at some distance from the sphere the rays are curved from their initial sphere-normal trajectory to become parallel with the z axis using a quadratic Bézier scheme. This preserves both ξ and η resolution near the sphere and provides excellent ξ resolution at great distances downstream.
The computational grid used in this work is 720 ξ by 1226 η grid points and extends 1200 sphere diameters (l) to the outer boundary. The smallest grid spacing normal to the sphere is η = 4.39 × 10 −4 l. For the extreme case Pe = 10 6 , there are 25 grid points within a chemical boundary layer of thickness 1.89 × 10 −2 l, calculated using the method displayed in figure 4. At the downstream edge of the domain ray curvature provides a smallest ξ grid spacing of ξ = 6 × 10 −3 l and 11 grid points within a sphere radius (a) of the z axis. The simulation can run up to a time Wt/2a = 1200 before isosurfaces that started at the sphere reach the boundary. The grid-independence of results is demonstrated in the Appendix.

Steady state
As the flow develops, the perturbation variables u , w , p and ρ can asymptote towards a steady state after an amount of time that depends on the non-dimensional parameters chosen. All results presented in this work are steady state in the frame of reference of the sphere, meaning that the spatial structure of perturbation variables u , w , ρ and gradients of total variables such as ∇ρ do not change with time. Because we are interested in the structure of chemical gradients that can propagate hundreds of sphere diameters downstream, simulations are run until a steady-state criterion of |f t+1 This criterion ensures that 99 % of the gradient structure is steady while simultaneously satisfying |g n+1 − g n | max /(g n max − g n min ) 10 −6 , where g represents u, p and ρ (cf. Torres et al. 2000). Simulations with parameters Re 10 and Fr = 10 −1 display unsteady, periodic behaviour. These simulations fail the steady-state criterion and are excluded from the results.
Note that the basic chemical fieldρ(z, t) and total fluid density ρ increase linearly with time (equations (2.4) and (2.6)), so that the shapes of all isosurfaces ρ + (∂ρ/∂z)Wt are constant, but the value associated with each advances monotonically. In this way isosurfaces separated by ρ = 1 also represent the time development of a single isosurface with time intervals of t = 1 as it moves past the sphere (Hanazaki et al. 2015). When reference is made to the shape of an isosurface at a given time, the initial time t = Wt/2a = 0 is defined as the time when the portion of the isosurface far from the sphere is located at z = 0.

Deformation of chemical gradients by sinking spheres
We begin by describing deformation of a linear chemical gradient by a sinking sphere and defining the δ ρ , L Def , Sh and M variables that will be used to quantify the influence of diffusivity, viscosity and buoyancy in subsequent sections. For now, we will treat a representative case where the chemical gradient induces negligible buoyancy forces: Re = 1, Fr = 10 3 and Sc = Pe = 10 3 . As explained in § 2.1, equation (2.11) is decoupled from (2.10) at large Fr and ρ behaves as a passive chemical subject to diffusion. A dimensional realization of the representative case is a 1 mm diameter particle sinking at 1 mm s −1 through a weak salinity gradient ∂ρ/∂z ≈ 10 −4 kg m −4 with diffusivity 10 −9 m 2 s −1 in a tank of water with viscosity 10 −6 m 2 s −1 . The buoyancy frequency is N = 10 −3 s −1 with gravity g ≈ 9.8 m s −2 and reference density ρ 0 = 1025 kg m −3 .
A sinking sphere sets the surrounding fluid into motion, dragging isosurfaces of chemical concentration downward ( figure 3). An isosurface that began upstream will stretch and wrap around the sphere until the divergence of the concentration gradient reaches a critical threshold ∇ 2 ρ = Pe (Hanazaki et al. 2015). This threshold is derived from the steady-state form of (2.11) applied at the sphere surface where (u, w) = 0. Molecular diffusivity then creates a hole in the isosurface at the leading edge of the sphere, which occurs just after t = 5 in this example (figure 3b). As the sphere punches through, the attached portion of the isosurface migrates rearward on the sphere and stretches vertically (figure 3c).
Although there is no flux of material through the surface of the sphere, wrapped isosurfaces form a thin boundary layer next to the sphere (figure 4a). This resembles the concentration boundary layer in the classic problem of mass transfer from a sphere, an analogy that will prove useful for explaining scalings in the next section. The thickness of the chemical boundary layer, δ ρ , is calculated as shown in figure 4(a). Because the gradient normal to the surface is zero, a line of best fit is drawn at the inflection point of the ρ profile. The intersection of the fitted line and the maximum value of ρ gives the local thickness of δ ρ (Verzicco & Camussi 1999;Breuer et al. 2004;Zhou & Xia 2010). Local values ±30 • of the equator are averaged to calculate δ ρ , avoiding the region of flow separation when Re 20. We also calculated the boundary layer thickness using methods based on percentage of outer ρ (Schlichting 1968), the inverse of the Sherwood number (Kiørboe, Ploug & Thygesen 2001) and the second moment of ρ (Weyburne 2006), and found them all to give similar results. An empirical viscous layer thickness δ u is calculated in a similar manner to the chemical boundary layer, as shown in figure 4(b).
From initial upstream contact to detachment downstream, a ρ isosurface is dragged and deformed by the sphere from its initial unperturbed position to a maximum vertical length L Def . In this case, this length is L Def ≈ 15 sphere diameters for the isosurface with value ρ Def = ρ + (∂ρ/∂z)Wt ≈ −15 (dashed line in figure 5a). After L Def is reached, the isosurface either detaches from the rear of the sphere, or pinches off downstream near the sphere and the hole in the isosurface closes. Reformed isosurfaces return to their unperturbed (flat) state at some distance downstream determined by a combination of molecular diffusion and buoyancy, the latter of which is negligible in this case. The slow relaxation of isosurfaces by diffusion is evident 200 sphere diameters downstream where isosurfaces are still vertically deformed by 7 sphere diameters ( figure 5, top panel). The entire trail of chemical disturbance is therefore much longer than the maximum deformed distance L Def .
The magnitude of the chemical gradient |∇ρ| normal to isosurfaces increases above the background state where contours come together and decreases where they spread apart during deformation by the sphere. The deviation of the magnitude of the gradient from the background state is here defined as this case, G is enhanced by a factor of ∼40 over the background gradient. On the z axis, the gradient is reduced from the background state downstream; here isosurfaces pinch off from the sphere and begin to return to their unperturbed state. In three dimensions, the tracks of deformed gradient form a cylindrical plume of gradient enhancement that surrounds a core of gradient reduction. These features stretch for hundreds of sphere diameters downstream until they are erased by diffusion (figure 5, top panel). This gradient deformation trail constitutes a standing structure like a wake, except that it persists to a greater degree in time and space than the velocity perturbations due to large Sc. Gradient enhancement is not balanced by gradient reduction over the course of the deformation event depicted in figure 5. This results in a net increase in the diffusive flux of chemical across isosurfaces that can be evaluated in two ways. Near the sphere, flux through the boundary layer is measured by the Sherwood number: where integration is performed over the surface of the sphere A, ∂ρ /∂n is the gradient normal to the sphere surface in the chemical boundary layer (see figure 4a), ρ ∞ = 0 is the value far from the sphere and ρ 0 is the value at the surface of the sphere (Friedlander 1957(Friedlander , 1961Sherwood, Pigford & Wilke 1975;Clift, Grace & Weber 1978;Karp-Boss, Boss & Jumars 1996). The numerator is the total flux through the thin chemical boundary layer where the velocity becomes vanishingly small at high Sc. Since there is no flux at the surface, we must calculate the values of ρ 0 from the simulation. The denominator is the steady-state diffusive flux in the absence of motion.
For the case shown in figure 5, the boundary layer flux is increased by a factor of Sh ≈ 4 compared to a stationary sphere. The total flux including the trail can be calculated in the Lagrangian frame of reference following an isosurface S that starts upstream, is deformed by the sphere and then returns to an unperturbed state downstream. Note that S is the chemical isosurface displayed in figure 3 at different times t . Because the gradient vector of magnitude G is always normal to an isosurface, and because the diffusivity is constant, G is proportional to the local net increase or decrease of chemical flux across an isosurface. Integrating G over isosurface S gives the flow rate F of material across the isosurface at a given time, F = S G dS. A time course of F following a single isosurface is shown in figure 6(a). The flow rate increases until t = L Def /W, after which the isosurface is no longer deformed by the sphere and diffusion slowly flattens the isosurface for another ∼450 time units. The total Lagrangian flux is summarized by where G is integrated over isosurface S out to the horizontal edge of the domain and over the time T for the isosurface to move from the sphere to the downstream boundary. Because G → 0 away from the deformation trail, the varying width of the curvilinear domain does not influence the calculation. The cumulative outer (time) integral of M is shown in figure 6(b). Deformations at times t < L Def /W contribute only 12 % of the total M, which emphasizes the importance of the long trail to the total flux in the system. For the representative case, M is an order of magnitude higher than the background, unperturbed state.
To separate the influences of viscosity, diffusivity and buoyancy forces on gradient deformation, we will compare simulations covering combinations of the three-dimensional parameter space 10 −2 Re 10 2 , 10 −1 Fr 10 3 and 10 Pe 10 6 . The Péclet number was chosen instead of the Schmidt number in order to treat viscous and diffusive effects separately. Over the full parameter space, the variables δ ρ , L Def , Sh and M show the greatest dependence on the Péclet number. As such, we address the influence of diffusivity on gradient deformation before moving on to viscosity and buoyancy.

The Péclet number dependence of gradient deformation
In the absence of buoyancy forces (Fr 10 3 ), equation (2.11) is decoupled from (2.10) and ρ behaves as a passive chemical. The Péclet number then controls the deformation of chemical isosurfaces and gradients based on the balance of terms in (2.11). We posit that the chemical boundary layer surrounding a sphere sinking through a linear gradient of chemical is analogous to the compositional boundary layer in the classic problem of mass transfer from a falling sphere. The Pe dependence of δ ρ , L Def , Sh and M can be explained using a boundary layer analysis similar to that of the mass transfer problem. As Pe increases, the chemical boundary layer near the equator thins (figure 7a). At low Pe, the w − 1 term of (2.11) is balanced by diffusion and the chemical boundary layer is weak or absent. For Pe > 10 2 , diffusion is small and the convective term u · ∇ρ dominates, strengthening the gradient normal to the sphere and decreasing δ ρ . The thickness of the boundary layer is estimated well by Pe −1/3 in the high-Pe regime (figure 7a). This scaling can be rationalized through dimensional analysis of the chemical boundary layer equation. At steady state, the chemical distribution very close to the sphere is determined by the standard boundary layer equation with an additional source term: where u is the velocity in the x direction parallel to the surface, v is the velocity in the n direction perpendicular to the surface and the source term γ is the dimensional form of w − 1 transformed into the (x, n) coordinate system (Landau & Lifshitz 1959;Levich 1962;Acrivos & Goddard 1965;Grossmann & Lohse 2000). This is in the frame of reference of the chemical boundary layer, close to the surface of the sphere. The advection, source and diffusion terms are all of the same order of magnitude in the chemical boundary layer for Pe 10 2 (figure 8a). We begin our analysis by balancing the advection and diffusion terms. The velocity parallel to the surface scales as v ∼ uδ ρ /l at large Pe, which confirms that ∂u/∂x ∼ Uδ ρ /δ u l is of the same order of magnitude as ∂v/∂n ∼ v/δ ρ in the chemical boundary layer. The line represents one-to-one correspondence.
When Sc is large, δ ρ < δ u and the parallel velocity is estimated by u ∼ Uδ ρ /δ u , where δ u is the thickness of the viscous boundary layer and U ∼ W is the velocity parallel to the surface at δ u (Levich 1962;Grossmann & Lohse 2000). The incompressibility condition (2.3) requires that ∂u/∂x ∼ Uδ ρ /δ u l is of the same order of magnitude as ∂v/∂n ∼ v/δ ρ in the chemical boundary layer. This implies that v ∼ uδ ρ /l, which is borne out by the simulations at large Pe in figure 8(b). The left-hand components of (3.4), u∂ρ /∂x ∼ Uδ ρ ρ /δ u l and v∂ρ /∂n ∼ vρ /δ ρ , are then of the same order of magnitude and only one is needed for the analysis (Levich 1962). The advection and diffusion terms of (3.4) are then estimated by At high Re, the viscous boundary layer scales according to the classic relation δ u /l ∼ Re −1/2 , and therefore (3.5) implies that δ ρ /l ∼ Re −1/2 Sc −1/3 (Lighthill 1950;Morgan & Warner 1956;Levich 1962;Grossmann & Lohse 2000). In the intermediate regime 10 −2 Re 10 2 addressed here, the viscous boundary layer is weak or absent and does not hold to the classic scaling (figure 7b). This is because the advective terms in the momentum equation (2.10) are subdominant and the pressure gradient is everywhere balanced by the viscous term (figure 7c). At the upper end of our regime (Re = 10 2 ), the advective terms are comparable to the other terms but Re is still not large enough for the Re −1/2 scaling to apply. From this we conclude that δ u is of the same order of magnitude as l at intermediate Re and the velocity near the sphere is estimated by u ∼ Uδ ρ /l. Equation This scaling agrees with the analyses of Acrivos & Taylor (1962) and Acrivos & Goddard (1965) for mass transfer from a sphere at small Re and large Pe in which they assumed δ ρ /l ∼ Sh −1 ∼ Pe −1/3 . The maximum dragged length L Def increases as the chemical boundary layer thins with increasing Pe (figure 9). At high Pe, this length scales as Pe 1/3 . The inverse relationship between δ ρ and L Def can be understood as follows. When a chemical isosurface approaches z = 0, the portion near the surface of the sphere enters the boundary layer while the rest continues to advect past the sphere. The portion near the sphere passes through the boundary layer until it pinches off behind the sphere at exactly the time that the rest of the surface reaches L Def (see figure 5a). In this situation the time scales for an isosurface to diffuse through the boundary layer and to traverse a distance L Def are comparable. The advective time scale for the isosurface to reach L Def is t Def ∼ L Def /W. The passage of an isosurface through the boundary layer occurs on the diffusion time scale t δ ∼ δ ρ 2 /κ since the velocity is small near the sphere. Setting t Def = t δ , L Def /W = δ ρ 2 /κ can be rewritten as (L Def /l)(l/δ ρ ) 2 = Pe. Substituting (3.7), L Def l ∼ Pe 1/3 . (3.8) The passage of the sphere causes the flux to increase across chemical isosurfaces in the boundary layer. The ratio of the total flux through the chemical boundary layer to the purely diffusive flux in the absence of motion is expected to scale to first order as Sh ∼ l/δ ρ at high Pe because the chemical flux near the sphere is confined to the boundary layer (Lévêque 1928;Levich 1962;Acrivos & Goddard 1965). This relation is confirmed in figure 10(b) and implies that the boundary layer Sherwood number scales as Sh ∼ Pe 1/3 (3.9) Chemical gradient deformation by sinking spheres . At the chemical boundary layer, the chemical anomaly scales as ρ ∼ Pe(u − 1)δ ρ 2 at large Pe, Fr = 10 3 and 10 −2 Re < 10 2 . The dashed line is the best fit.
at large Pe (figure 10a). The analogy does not hold at low Pe because there is no boundary layer or corresponding chemical flux and so the simulations deviate from the classic result Sh = 2 as Pe → 0. After pinching-off near the rear of the sphere, isosurfaces return to their original positions through diffusion over a time scale much greater than the boundary layer time scale, t δ . The flux of chemical across these isosurfaces is elevated compared to the background throughout the deformation trail. At large Pe, the total flux of material due to enhanced gradients in the trail, M, exhibits slopes that fit Pe 2/3 (figure 10c). Here Sh differs from M because the former is a property of the chemical boundary layer in the frame of reference of the sphere (Eulerian) and the latter is calculated moving with the isosurfaces (Lagrangian). Following an isosurface, diffusion must occur over a distance L Def to flatten that isosurface when it reaches the end of the trail at time t M ∼ L Def 2 /κ. Since the integrand F in (3.3) is essentially constant between L Def and the end of the trail (see figure 6a), the definite integral of M can be approximated The analogy between mass flux from a sphere and a sphere sinking through a gradient provides a rationale for the scalings l/δ ρ ∼ Sh ∼ L Def /l ∼ M 1/2 observed in the simulations. This approach is successful at high Pe but does not hold at low Pe where the chemical boundary layer is absent.
Returning to (3.4), we now consider the balance of the source and diffusion terms γ ∼ κ(∂ 2 ρ /∂n 2 ) (Zhang et al. 2019). In non-dimensional form, this balance is estimated by u − 1 ∼ ρ /(Peδ ρ 2 ), where u(x, n) ∼ w(r, z) is the velocity parallel to the surface near the equator of the sphere. The resulting scaling ρ ∼ Pe(u − 1)δ ρ 2 agrees with the numerical simulations as shown in figure 11.

Reynolds number flow regimes and gradient deformation
Fluid motion around a sphere between Re = 10 −2 and 10 2 is marked by the well-known transition at Re ≈ 24 from laminar, viscosity-dominated flow to the formation of axisymmetric standing vortices in the wake (e.g. Taneda 1956;Magnaudet, Rivero & Fabre 1995). This change in flow regime is evident in the shape of chemical isosurfaces and gradients in figure 12, in which Pe and Fr are both held at 10 3 . The shape of the toroidal vortex in our simulations agrees with the experimental measurements of Taneda (1956). Recirculation in the standing vortices leads to a steady-state pattern of gradient deformation that is markedly different from lower-Re simulations.
Instead of migrating all the way around the sphere, attached isosurfaces accumulate near the separation streamline (dashed red line in figure 12a). Once an isosurface reaches L Def = 19.7, it pinches off at z = 4 and the unattached portion continues to advance downstream in a similar fashion to the lower-Reynolds-number cases. The portion of the isosurface that remains attached to the sphere is then partially The change in gradient deformation over 10 −2 Re 10 2 holding Pe constant is minor (figure 13a). However, L Def , Sh and M exhibit higher Pe scalings at Re = 10 2 than explained in the previous section (figure 13b). The standing vortices delay 'pinching-off' of contours at the rear of the sphere, increasing L Def . Enhanced chemical gradients along the edge of the vortices contribute to larger Sh and M. The change in Pe dependency indicates that the scaling behaviour observed at Re = 10 2 lies in between the δ ρ /l ∼ Pe −1/3 = Re −1/3 Sc −1/3 relation described here and the classic δ ρ /l ∼ Re −1/2 Sc −1/3 result found at larger Re. We expect the δ u /l ∼ Re −1/2 scaling at higher Re, but this remains to be tested. Therefore we posit that Re = 10 2 is intermediate and suggest a continuous (rather than sharp) transition of l/δ ρ ∼ Sh ∼ L Def /l ∼ M 1/2 between the asymptotic limits Re 1 and Re 1.

Suppression of downstream chemical gradients by buoyancy forces
For Fr < 10 3 , equations (2.10) and (2.11) are recoupled through the buoyancy term and displaced isosurfaces influence the velocity field. Experiments by Yick et al. (2009) showed that buoyancy forces reduce the distance that isosurfaces are dragged by the sphere at low Re. This is evident comparing Fr = 10 3 in figure 5(a) to Fr = 1 in figure 14(a). Buoyancy constrains the negative velocity anomaly to a narrow region close to the sphere (figure 14c). Attached isosurfaces then migrate around the sphere more quickly and detach sooner, reducing L Def . At the same time, buoyant shear near the equator slightly thins the chemical boundary layer. The buoyancy force accelerates the return of isosurfaces and generates a positive w − 1 anomaly along the z axis, i.e. faster than the background flow W (figure 14c). The trail disappears by z ≈ 20, which is an order of magnitude shorter than the trail observed for Fr = 10 3 where isosurfaces flatten only through diffusion ( figure 5). The influences of buoyancy forces on the deformation variables for Re = 1 and Pe = 10 3 are shown in figure 15. Above Fr = 10, buoyancy forces become insignificant and Pe dominates as described in § 3.2. For Fr < 10, a thinner chemical boundary layer results in a slightly larger Sh by definition, while a smaller L Def and a shorter trail cause M to decrease substantially. Yick et al. (2009) provided an empirical estimate of L Def /l ∼ Fr 1/2 for Re 1, Fr 1 and Pe < 10 3 , which fits the case shown in figure 15. The boundary layer argument we develop in § 3.2 can be generalized for scaling the deformation variables in the intermediate Re, Fr < 10 and large Pe regime addressed here. First, we assume that the viscous boundary layer thickness is represented by the natural length scale of a viscous and buoyant flow, δ u /l ∼ (ν/N) 1/2 /l = (Fr/Re) 1/2 (Gargett 1988;Yick et al. 2009). Figure 16(a) shows that this scaling is borne out by the numerical results. The boundary layer equation (3.4) then has dimensions of and δ ρ l ∼ Fr 1/6 Re 1/6 Pe 1/3 . (3.12) We note that (3.12) differs from Hanazaki et al. (2015) because they define δ ρ as the half-width at half-maximum of ∇ 2 ρ /Pe in the downstream jet when Re > 10 2 . Since the diffusive boundary layer time scale is equal to the advective L Def time scale (see § 3.2), L Def /W = δ ρ 2 /κ and  Fr 1/3 Pe 2/3 Re 1/3 TABLE 1. Scalings for the variables δ ρ , L Def , Sh and M for Pe > 10 2 and 10 −2 Re 10 2 . All scalings are analytical except M for Fr < 10, which is empirical.

Discussion
Deformation of chemical gradients across the three-dimensional parameter space 10 −2 Re 10 2 , 10 Fr 10 3 and 10 2 < Pe 10 6 is primarily controlled by the size and speed of the sphere versus the diffusivity of the chemical. Employing a boundary layer analysis based on the classic problem of mass transfer from a sphere, we derive the new relation l/δ ρ ∼ Sh ∼ L Def /l ∼ M 1/2 ∼ Pe 1/3 for spheres sinking through chemical gradients at large Pe in the absence of buoyancy forces (table 1). The stronger Pe dependence of the new Lagrangian flux quantity, M, versus the familiar Sh indicates that the trail of deformed gradients is an important consideration when determining the total chemical flux in the system. An increased Pe dependence is observed at Re = 10 2 due to strengthened gradients along the edges of standing vortices.
Below Fr = 10 3 , the momentum and composition equations are coupled and variation in ρ contributes to buoyancy forces. Buoyancy forces in the 10 Fr < 10 3 regime do not influence gradient deformation. Applying dimensional analysis of the boundary layer equation for large Pe and Fr < 10, we find an additional Fr dependence for the deformation variables while preserving the Pe dependence from the case of no buoyancy forces. Holding Pe constant and Fr < 10, our simulations agree with L Def /l ∼ Fr 1/2 empirically determined by Yick et al. (2009) for Re < 1 and L Def /l ∼ πFr proposed by Hanazaki et al. (2015) for Re 10 2 . Our work provides an analytical basis for a Pe dependence of L Def /l ∼ Re −1/3 Fr 1/3 Pe 1/3 that expands the regime to large Pe, 1 < Fr < 10 and intermediate Re. While Sh increases with buoyancy forces because of boundary layer thinning, both L Def and M are suppressed by buoyancy forces downstream. In a recent study by Zhang et al. (2019), the authors derive a buoyancy length scale s in the asymptotic Re 1 and Re 1 regimes in order to isolate the origin of the drag force on a sphere settling through density stratification. As s is distinct from the assumed δ u and derived δ ρ length scales at low Fr in § 3.4, we consider our scalings complementary. Additionally, the scaling laws in this work present a comprehensive description of a system that can act as a reference for more complex systems, such as non-spherical particles and unsteady settling.
These combined results imply that significant deformation of chemical gradients by moving particles occurs in natural environments when Pe > 100 and Fr > 10. Oceanic particles, for example, span an immense range of sizes (from 1 µm to 1 cm) and sinking rates, which correspond to 10 −2 < W/l < 10 s −1 (Stemmann, Jackson & Ianson 2004;McDonnell & Buesseler 2010). The buoyancy frequency of the upper ocean is typically 10 −3 < N < 10 −2 s −1 (Emery, Lee & Magaard 1984;Houry et al. 1987) and particle Froude numbers fall within the range 1 < Fr < 10 4 . Stratification is unlikely to suppress deformation of chemical gradients by particles with a sinking speed to size ratio of W/a > 10 −1 s −1 . The Péclet number for fluid transported vertically by a particle in stratification depends on the particle size and speed and the diffusivity of the stratifying agent. The diffusivity coefficient for salinity stratification is κ S ≈ 10 −9 m 2 s −1 (Poisson & Papaud 1983). The particle size and speed must be Wl > 10 −7 m 2 s −1 for significant deformation of the salinity gradient to occur. The diffusivity coefficient for thermal stratification is κ T ≈ 10 −7 m 2 s −1 and then Wl > 10 −5 m 2 s −1 . Although statistics for particle size and speed vary significantly throughout the oceans, these calculations suggest that many marine particles can deform salinity stratification and, to a lesser extent, thermal stratification.
Gradients of inorganic nutrients and dissolved organic matter are important for microscale ecology in the oceans (Azam 1998). These molecules are larger than sodium and chloride and so their diffusivities are smaller than 10 −9 m 2 s −1 . The diffusion coefficients of nutrients such as phosphate, sulphate and iron are ∼10 −10 m 2 s −1 (Li & Gregory 1974) and polysaccharides can have diffusivities as low as 10 −11 m 2 s −1 (Callaghan & Lelievre 1986). Large-scale vertical gradients of these compounds exist across the pycnocline in most of the ocean. Therefore particles sinking through weak stratification are even more likely to deform weak gradients of nutrients and organic compounds, enhancing the local diffusive flux by a factor of M to the potential benefit of nearby microbes. While swimming motions were not investigated in this work, models of motile microorganisms 'squirming' through weak stratification appear to cause patterns of chemical deformation that differ from those of a passively sinking sphere (e.g. Doostmohammadi, Stocker & Ardekani 2012). The potential for plankton to modify local chemical gradients and fluxes is an intriguing subject for future work.

Conclusions
We have presented a mechanism for deformation of an ambient chemical gradient by a sphere and the resulting trail of gradient enhancement that was imaged in previous laboratory experiments but never described (Yick et al. 2009). Using the analogy of mass flux from a sphere, we compare dimensional analysis of the boundary layer equation and numerical simulations to quantify gradient deformation in the regime of large Pe and intermediate Re. The chemical boundary layer thickness, Sherwood number and maximum dragged distance of an isosurface scale as l/δ ρ ∼ Sh ∼ L Def /l ∼ Pe 1/3 in the absence of buoyancy forces. The relation δ ρ /l ∼ Pe −1/3 arises from the classic balance between advection and diffusion within the chemical boundary layer surrounding the sphere. The time scale for an isosurface to diffuse through the boundary layer is equal to the advective time scale for an isosurface to travel a distance L Def , which leads to the scaling L Def /l ∼ Pe 1/3 . The Lagrangian flux following a chemical isosurface is M ∼ Pe 2/3 > Sh due to enhanced gradients in the trail that persist over a diffusive time scale proportional to (L Def /l) 2 . At Fr < 10 3 , disturbance of the chemical gradient causes buoyant velocities. For 10 Fr < 10 3 , buoyancy forces have little effect on gradient deformation. For Fr < 10 and large Pe, a similar analysis yields l/δ ρ ∼ Sh ∼ Re −1/6 Fr 1/6 Pe −1/3 and L Def /l ∼ Re −1/3 Fr 1/3 Pe 1/3 . The trail is greatly diminished by strong buoyancy forces and an empirical relation for the Lagrangian flux is M ∼ Re −1/3 Fr 1/3 Pe 2/3 . Gradient and diffusive flux enhancement by sinking particles may contribute to the microscale ecology of the oceans. were conducted on the Comet resource located at the San Diego Supercomputer Center through the generous support of NSF XSEDE grant TG-OCE160016. C.J.D. acknowledges a Natural Environment Research Council personal fellowship (reference NE/L011328/1). We are grateful to three reviewers and the editor for their detailed and constructive comments, which helped to improve the manuscript.

Declaration of interests
The authors report no conflict of interest.

Appendix. Validation of grid-independence
To ensure that the results described in this work are grid-independent, simulations of what is assumed to be the most demanding case (Re = 1, Fr = 10 3 , Pe = 10 6 ) were computed on two additional grids with different resolutions but the same physical domain size as the grid detailed in § 2.3. Grid A comprises 360 ξ by 613 η grid points with the smallest grid spacing normal to the sphere η = 8.79 × 10 −4 l. Grid B consists of 1440 ξ by 2451 η grid points with smallest spacing η = 2.20 × 10 −4 l. Grids A and B have roughly half and double the spatial resolution, respectively, of the original grid and use the same clustering algorithm. Simulations reach steady state ( § 2.4) at approximately the same non-dimensional time t for all grids, whereas the computational time was 2 days for grid A, 1 week for the original grid and 1 month for grid B. Similar to figure 4, chemical and velocity boundary layer profiles and thickness calculations are displayed for the extreme case in figure 17. While there are slight differences between the chemical profiles near the sphere surface in figure 17(a), the calculated δ ρ /l differs by 2 % between the different grids. With negligible differences between δ ρ and δ u calculated at different resolutions, the original grid sufficiently resolves the boundary layer and the grid-independence of the results is verified.