1. Introduction
Settling of porous particles in stratified fluids is a ubiquitous process in many environments, for example in the ocean (Asper et al. Reference Asper, Deuser, Knauer and Lohrenz1992; Li, Yuan & Wang Reference Li, Yuan and Wang2003). Porous marine particles, ‘marine aggregates’, play a fundamental role in oceanic biogeochemical cycles as they withdraw nutrients and especially carbon from surface waters into deeper regions of the ocean. The sinking of marine aggregates and their remineralisation during their descent eventually determine how efficiently the biological carbon pump transfers photosynthetically fixed carbon to depth. The accurate representation of the biological carbon pump in Earth system models is thus a key challenge to improve projections of global biogeochemical cycles and, in particular, the oceanic carbon sink (Ilyina & Friedlingstein Reference Ilyina and Friedlingstein2016; Maerz et al. Reference Maerz, Six, Stemmler, Ahmerkamp and Ilyina2020). Since ocean stratification is projected to intensify under ongoing climate change (Bopp et al. Reference Bopp2013), insights into effects of stratification on settling dynamics of marine aggregates are required to enable quantifying its potential effect on the future biological carbon pump.
In situ observation of marine aggregates settling at low to moderate Reynolds numbers 
$\textit {O}(0.01)\leq \textit {Re}\leq \textit {O}(10)$ showed increased retention times in stratification (Asper Reference Asper1985; Alldredge & Gotschalk Reference Alldredge and Gotschalk1988; Riebesell Reference Riebesell1992; Alldredge & Crocker Reference Alldredge and Crocker1995; MacIntyre, Alldredge & Gotschalk Reference MacIntyre, Alldredge and Gotschalk1995; Alldredge et al. Reference Alldredge, Cowles, McIntyre, Rins, Donaghay, Greenlaw, Holliday, Dekshenieks, Sullivan and Zaneveld2002). We here define the Reynolds number as 
$\textit {Re}=aW/\nu$ for sphere radius 
$a$, stationary sinking velocity 
$W$ and kinematic molecular viscosity 
$\nu$. Stratification, expressed in terms of the Brunt–Väisälä frequency, 
$N$, typically ranges from 
$N=\textit {O}(0.001\,\textrm {s}^{-1})$–
$\textit {O}(0.1\,\textrm {s}^{-1})$ under oceanic conditions up to 
$N\lesssim \textit {O}(1\,\textrm {s}^{-1})$ in estuaries and fjords (Boehrer & Schultze Reference Boehrer and Schultze2008; Geyer, Scully & Ralston Reference Geyer, Scully and Ralston2008). The Brunt–Väisälä frequency is defined as 
$N=\sqrt {-(g/\rho _0)\gamma }$, where 
$g$ is the gravitational acceleration, 
$\rho _0$ a reference density and 
$\gamma =\partial \rho /\partial z$ is the vertical density gradient of the ambient water. Studies of in situ marine aggregates suggested cessation of settling due to reaching neutral buoyancy (e.g. Riebesell Reference Riebesell1992; Alldredge & Crocker Reference Alldredge and Crocker1995), or the buoyancy of the interstitial fluid carried by the porous particles to cause increased retention during settling in stratification (Kindler, Khalili & Stocker Reference Kindler, Khalili and Stocker2010). Such retention is particularly associated with pycnocline layers, which are defined through strong vertical density gradients (MacIntyre et al. Reference MacIntyre, Alldredge and Gotschalk1995; Alldredge et al. Reference Alldredge, Cowles, McIntyre, Rins, Donaghay, Greenlaw, Holliday, Dekshenieks, Sullivan and Zaneveld2002).
Settling of porous particles across sharp pycnocline layers has thus far been studied in models and experiments to understand the underlying fluid dynamics (Kindler et al. Reference Kindler, Khalili and Stocker2010; Camassa et al. Reference Camassa, Khatri, McLaughlin, Prairie, White and Yu2013; Prairie et al. Reference Prairie, Ziervogel, Arnosti, Camassa, Falcon, Khatri, McLaughlin, White and Yu2013, Reference Prairie, Ziervogel, Camassa, McLaughlin, White, Dewald and Arnosti2015; Panah, Blanchette & Khatri Reference Panah, Blanchette and Khatri2017). These sharp pycnocline layers extend over a thickness of order of magnitude similar to the particle radius and separate an upper, lighter fluid phase from a lower, denser one. For the settling dynamics across these sharp pycnocline layers, a number of cases can be distinguished depending on the density difference ratio 
$\xi =\delta \rho /\Delta \rho _0$ (comparing the density increment over the pycnocline layer, 
$\delta \rho$, with a particle's excess density in the lighter phase, 
$\Delta \rho _0$) and the particle permeability to flow 
$k$ (Kindler et al. Reference Kindler, Khalili and Stocker2010; Camassa et al. Reference Camassa, Khatri, McLaughlin, Prairie, White and Yu2013; Prairie et al. Reference Prairie, Ziervogel, Arnosti, Camassa, Falcon, Khatri, McLaughlin, White and Yu2013, Reference Prairie, Ziervogel, Camassa, McLaughlin, White, Dewald and Arnosti2015; Panah et al. Reference Panah, Blanchette and Khatri2017; Prairie et al. Reference Prairie, Montgommery, Proctor and Ghioroso2019). In the case of impermeability to flow (which is a fair approximation for most marine aggregates; Moradi et al. Reference Moradi, Liu, Iversen, Kuypers, Ploug and Khalili2018) the external flow and density fields can be expected to be largely similar to those of solid particles. Solid spheres settling in linear stratification at low Reynolds number (
$\textit {Re}\leq 1$) and large Schmidt number (
$\textit {Sc}=\nu /D$, representative of salt as a stratifying agent) entrain lighter fluid by distorting isopycnals at the sphere surface which leads to the formation of a density boundary layer around a particle and a buoyant wake (figure 1; e.g. Srdić-Mitrović, Mohamed & Fernando Reference Srdić-Mitrović, Mohamed and Fernando1999; Higginson, Dalziel & Linden Reference Higginson, Dalziel and Linden2003; Yick et al. Reference Yick, Stocker, Peakock and Torres2009; Zhang, Mercier & Magnaudet Reference Zhang, Mercier and Magnaudet2019). Stratification effects also cause extra stresses at the sphere's surface by the formation of toroidal vortices due to baroclinic torque (List Reference List1971; Ardekani & Stocker Reference Ardekani and Stocker2010; Zhang et al. Reference Zhang, Mercier and Magnaudet2019) which has been shown to be the dominant contribution to drag enhancement at low to intermediate 
$\textit {Re}$ (Zhang et al. Reference Zhang, Mercier and Magnaudet2019). The interplay of vorticity generation and the buoyancy by lighter fluid dragged by the particles leads to distinct stratification, i.e. settling regimes separated by the relative length scales of stratification, viscosity and diffusivity, in which drag enhancement can be described by scaling laws based on 
$\textit {Re}$, 
$\textit {Sc}$ and the Froude number, 
$\textit {Fr}=W/(aN)$ (Zhang et al. Reference Zhang, Mercier and Magnaudet2019).
Figure 1. Conceptual representation of stationary, stratified flow around a porous and impermeable sphere of radius 
$a$ and density 
$\rho _p$. At constant velocity 
$W$, the impact of isopycnal deflection on the drag exerted on the sphere is modelled as a density boundary layer of thickness 
$\delta$ and density 
$\rho _\delta$. The red dot indicates the stagnation point.
In addition to these external effects commonly referred to as stratification drag enhancement, the excess density of porous particles changes when settling in stratification via an adjustment of interstitial and ambient fluid through the exchange of stratifying agent. This mass adjustment effect is most pronounced when particles are impermeable to flow, i.e. the exchange of interstitial pore water and ambient fluid is governed only by diffusion. If the particle excess density is large compared to that of the lower fluid phase, 
$\xi <1$, the influence of the interstitial fluid mass adaptation is negligible and only slight retention is observed due to viscous entrainment of fluid of lower density into the denser phase forming a buoyant wake which decelerates the particle descent below a sharp pycnocline (Srdić-Mitrović et al. Reference Srdić-Mitrović, Mohamed and Fernando1999; Camassa et al. Reference Camassa, Falcon, Lin, McLaughlin and Mykins2010). At moderate Reynolds numbers (
$5<\textit {Re}<15$), the wake pinches off rapidly and retreats upwards, while in the Stokes regime (
$\textit {Re}\ll 1$), the stem of entrained fluid dissolves due to diffusion and, hence, retardation increases (Camassa et al. Reference Camassa, Falcon, Lin, McLaughlin and Mykins2010). If, however, a particle is impermeable to flow and the density increment exceeds the excess density of the particle, hence 
$\xi >1$, settling itself becomes diffusion-limited (Li et al. Reference Li, Yuan and Wang2003; Kindler et al. Reference Kindler, Khalili and Stocker2010; Panah et al. Reference Panah, Blanchette and Khatri2017), i.e. particles settle in response to the diffusive mass adaptation of the interstitial fluid.
As stratification in the ocean typically extends over several to tens of metres (Li et al. Reference Li, Cheng, Zhu, Trenberth, Mann and Abraham2020), which is much larger than typical sizes of marine aggregates (
$\textit {O}(1\,\mathrm {\mu }\textrm {m})\leq a\leq \textit {O}(1\,\textrm {cm})$), viscous and diffusive effects on settling velocities of marine aggregates may not be as separable as suggested by models for sharp pycnocline layers. On the contrary, it is reasonable to assume that in extended, linear gradients a porous sphere reaches a stationary state (
$W=\textrm {const.}$) in which both drag enhancement and delayed interstitial mass adaptation contribute to reducing settling velocities as compared to settling in unstratified fluids. As an additional effect, the density boundary layer can be expected to modulate the interstitial mass adaptation with the ambient fluid.
Understanding the consequences of stratification-dependent sinking of marine aggregates for biogeochemical cycles is highly desirable, but thus far has been impeded by limitations inherent to in situ observations and large uncertainties of key parameters such as aggregate excess density, which demands model-based investigations. Ocean biogeochemical cycles are, however, typically studied in regional or Earth system models which are limited by computational performance and thus require parameterisation of subgrid-resolution processes (e.g. the Max Planck Institute's Earth system model 1.2-LR features a nominal horizontal resolution of 150 km and a minimum vertical layer thickness of 10 m; Mauritsen et al. Reference Mauritsen2019). A parameterisation of net effects of stratification on marine aggregate sinking is thus essential to enable Earth system models to investigate this effect on future biogeochemical cycles and carbon fluxes under ongoing climate change. Therefore, we here focus on highly porous, impermeable particles settling in linear stratification and aim at (i) detailed insights into the involved forces and (ii) providing approximate parameterisations that enable one to represent and study linear stratification effects on settling of marine aggregates in large-scale, spatially explicit regional to global ocean biogeochemistry models.
To determine the relationship between stratification drag enhancement and mass adaptation, we examined the settling dynamics of porous (
$\epsilon =0.95$, ratio of void to total volume) and impermeable spheres (here defined as 
$k\lesssim 10^{-12}\,\textrm {m}$; see also figure 14) in linear stratification for low Reynolds numbers and a wide range of Froude numbers as commonly found in marine environments. Steady-state lattice Boltzmann simulations and schlieren experiments were performed in which we focused on buoyancy effects of the interstitial and boundary layer fluid, as well as the impact of diffusive exchange between interstitial, boundary layer and ambient fluid. To simplify the quantification of increased retention of particles settling in stratification, we derive scaling laws and present their application to several test cases. This article is structured as follows. In § 2 the experiment method, basic equations and the simulation method are described along with some validation results. Subsequently, in § 3 flow and density fields of porous particles are characterised and referenced to solid-sphere behaviour with particular emphasis on the density boundary layer and the interstitial fluid properties as well as their impact on settling velocities. Scaling laws for both stratification drag enhancement and particle buoyancy are derived. Finally, in § 4 scaling analysis is applied to illustrate and discuss the retention of porous and solid particles in linear stratification under oceanic conditions.
2. Methods
2.1. Settling velocity measurements and schlieren visualisation
In order to simultaneously determine settling velocity and visualise the density perturbations of a porous particle settling through a linear stratified fluid, we conducted schlieren experiments in a settling chamber. The settling chamber had a square base of 
$d_c=22.4\,\textrm {cm}$ with height 
$h_c=38\,\textrm {cm}$ (figure 2). The tank was filled to a depth of 37 cm using a double bucket system to generate a linear stratification based on two fluids (Oster Reference Oster1965). The fluids in the buckets were aqueous glycerol mixtures with contents that ranged from 
$60$ to 
$90\,\text {wt}{\%}$ glycerol. The glycerol content as well as stratification strengths were varied to adjust Reynolds and Froude numbers of the particles. The mean densities throughout all experiments varied between 
$1150$ and 
$1227\,\textrm {kg}\,\textrm {m}^{-3}$. Note that, within one experiment, density changes were typically below 
$1\,\%$. To avoid evaporation, the fluid surface was covered with shading balls. Experiments were performed in an isolated basement room where temperature variations did not exceed 
$\Delta T=\pm 0.5\,^{\circ }$C.
Figure 2. Experimental method for simultaneous measurement of settling velocity and schlieren visualisation.
Highly porous, but impermeable particles were produced using fibres as described in Dörgens et al. (Reference Dörgens, Ahmerkamp, Müssig, Stocker, Kuypers, Khalili and Kindler2015). Briefly, in a semi-manual process, polyester fibres with 
$20\,\mathrm {\mu }\textrm {m}$ diameter were agglomerated to nearly perfect spheres. The porosities of 10 replicates were in the range 
$\epsilon = 92\,\%$ to 94 % with mean radius 
$a =0.68\,\textrm {cm}$. Particles were released into the flow tank using a cone with a prolonged inlet to avoid lateral movements and rotations. The settling was imaged using a pco1600 camera (PCO) fitted with an AF-S Nikkor zoom lens (
$f=16$ to 
$88\,\textrm {mm}$, 
$f$ number 
$=3.522$) positioned 20 cm in front of the tank. The rear surface of the tank was coloured black. The field of view was 
$20\,\textrm {cm}\times 5\,\textrm {cm}$, at a resolution of 
$27\,\mathrm {\mu }\textrm {m}\,\textrm {px}^{-1}$. A sequence of 7 to 35 images was recorded at 2 Hz and processed using Simulink Matlab (2011b). Prior to analysis, a reference image (
${\tilde {I}_{ref}}$) was subtracted from each image including the aggregate (
${I}$) and used to normalise: 
${\tilde {I} = (I-I_{ref})/I_{ref}}$. In the reference image 
${\tilde {I}}$, the largest connected white area was dissected by tracing the boundaries. The centre of the particle was identified as the centre of gravity of an ellipse fitted to the boundary. The sinking velocity was determined from a sequence of images by calculating the five-point central difference of the particle's position over time.
The schlieren method allows for the reconstruction of density perturbations based on refractive index changes (e.g. Yick, Stocker & Peakock Reference Yick, Stocker and Peakock2007). These refractive index changes were visualised using a pattern consisting of randomly distributed dots with varying brightness. The pattern was printed onto waterproof transparent paper and attached to the inner wall of the flow tank, and laterally illuminated by two diffuse 150 W halogen lamps. The pattern was imaged at 1 Hz using an SLR camera (Nikon D90) fitted with an AF-S Nikkor lens (
${f}=45\,\textrm {mm}$, 
$f$ number 
$=1.4$ to 16) installed 
$L_c=20\,\textrm {cm}$ in front of the tank. The field of view was 
$10\,\textrm {cm}\times 5\,\textrm {cm}$, at a resolution of 
$49\,\mathrm {\mu }\textrm {m}\,\textrm {px}^{-1}$. As reference, 10 images were obtained and averaged in an undisturbed situation. Subsequently, a particle was released and the density-disturbed image was recorded and cross-correlated with the reference image using PIV View 2C (Pivtec, Göttingen, Germany) to determine the displacement of the random pattern. The interrogation windows for the cross-correlation consisted of 
$24\,\textrm {px} \times 24\,\textrm {px}$ with 50 % overlap, resulting in a 
$12\,\textrm {px} \times 12\,\textrm {px}$ grid. In this set-up, the maximal traceable shifts were 
${\sim }10\,\mathrm {\mu }\textrm {m}$ at a spatial resolution of 
$588\,\mathrm {\mu }\textrm {m}$. Based on the displacement image, the density field was reconstructed using a tomographic algorithm (Appendix A). Experiments were performed in 13 density stratifications of different strengths with 10 replicate particles and were optimised to resemble the Reynolds and Froude numbers of the numerical model. Schmidt numbers were above 700 and of 
$\textit {O}(1000 - 10\,000)$. Permeabilities of the fibre particles were measured to be 
$k=3.87\times 10^{-12}\,\textrm {m}$ resulting in Darcy numbers of 
$\textit {O}(10^{-6})$ (Dörgens et al. Reference Dörgens, Ahmerkamp, Müssig, Stocker, Kuypers, Khalili and Kindler2015); see also figure 13 and table 1 for an overview of the experiment parameters.
Table 1. Parameters of the schlieren and settling experiments, where 
$\nu$ is the kinematic viscosity, 
$\rho$ the mean density, 
$N$ the buoyancy frequency and 
$W$ the settling velocity. Parameter 
$R_H$ is the theoretical Reynolds number in the absence of stratification. Each of the experiments was performed for 3 to 5 individual particles.
2.2. Model formulation
High-resolution numerical simulations were performed to obtain the density and flow fields. The model was adapted from an earlier formulation used for the investigation of stratification drag enhancement of non-porous spheres (Torres et al. Reference Torres, Hanazaki, Ochoa, Castillo and van Woert2000; Larrazabal, Torres & Castillo Reference Larrazabal, Torres and Castillo2003; Yick et al. Reference Yick, Stocker, Peakock and Torres2009; Liu, Kindler & Khalili Reference Liu, Kindler and Khalili2012). We considered the flow of a linearly stratified fluid at constant velocity 
$W$ past a stationary sphere:
$$\begin{gather} \frac{\partial\boldsymbol{u}}{\partial t}+\frac{1}{\epsilon}\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}={-}\epsilon \boldsymbol{\nabla}p-\frac{\rho'}{\textit{Fr}^2}\boldsymbol{j}+\frac{1}{\textit{Re}}\nabla^2 \boldsymbol{u}-A\frac{\epsilon}{\textit{Da}\textit{Re}} \boldsymbol{u}, \end{gather}$$
$$\begin{gather}\frac{\partial{\rho'}}{\partial t}+\frac{1}{\epsilon}\boldsymbol{u}\boldsymbol{\cdot}{\rho'}= w-1+ \frac{1}{\textit{Re}\textit{Sc}}\nabla^2\rho', \end{gather}$$
where length is scaled by 
$a$, velocity by the undisturbed velocity 
$W$ (resembling the constant settling velocity), pressure perturbations by 
$\rho _0 W^2$ and density perturbations 
$\rho '$ by 
$-a\gamma$. The density perturbation represents the density contrast induced by the particle in comparison to an undisturbed linearly stratified fluid: 
$\rho '= (\rho -\rho _B) / (-a\gamma ) = \Delta \rho / (-a\gamma$), where 
$\rho _B$ is the density of the undisturbed ambient fluid in linear stratification and 
$-a\gamma$ is the vertical density change along one particle radius. Therefore the density contrast 
$\Delta \rho / (-a\gamma )$ can be interpreted as the deflection distance of the pycnoclines from their equilibrium position.
Here, 
$\boldsymbol {u}=(u,w)$ is the fluid velocity vector with radial and vertical components, 
$p$ is the pressure, 
$\boldsymbol {j}$ is the vertical unit vector, positive upwards, 
$\epsilon$ is the porosity of the particle and 
$\textit {Da}=k/a^2$ is the Darcy number, with 
$k$ the permeability. The Darcy number was adjusted to 
$\textit {Da}=10^{-12}$, which is comparable to permeabilities of 
$k=10^{-18}$ to 
$10^{-16}\,\textrm {m}$ (for particle radii in the range 
$a=10^{-4}$ to 
$10^{-3}\,\textrm {m}$) and, therefore, resembling that for typical impermeable marine particles (Kiørboe et al. Reference Kiørboe, Grossart, Ploug and Tang2002) (see also figure 14 for tests on permeability effects). The fluid domain is defined by 
$\epsilon =1, A=0$, where 
$A$ is a binary switch, and the particle by 
$\epsilon =0.95$, 
$A=1$ (Basu & Khalili Reference Basu and Khalili1999). For the numerical implementation a quasi-two-dimensional model was employed using a lattice Boltzmann method for axial symmetric flows (see figure B for numerical implementation, grid system and boundary conditions). Model convergence was ensured through the drag coefficient and interstitial density. When changes were less than 
$10^{-7}$ between 1000 time steps, the model was considered to be converged (figure 17). Simulations were performed for both porous and solid spheres for the parameter regime spanned by 
$0.05\leq \textit {Re}\leq 10$ and 
$0.1\leq \textit {Fr}\leq 100$ at Schmidt number 
$\textit {Sc}=\nu / D=700$, where 
$\nu$ is the kinematic viscosity and 
$D$ the diffusion coefficient.
The momentum-exchange and stress-integral methods were employed to evaluate the force on the sphere. Namely, the drag coefficient 
$C_D^S$ was computed as the sum of the pressure and viscous drag coefficients, 
$C_P^S$ and 
$C_V^S$, respectively:
$$\begin{gather} C_P^S={-}\frac{1}{\dfrac{1}{2}\rho W^2{\rm \pi} a^2}\int_Sp\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{j}\,\textrm{d} S, \end{gather}$$
$$\begin{gather}C_V^S=\frac{1}{\dfrac{1}{2}\rho W^2{\rm \pi} a^2}\int_S\mu\boldsymbol{n}\boldsymbol{\cdot}((\boldsymbol{\nabla}\boldsymbol{u})+(\boldsymbol{\nabla}\boldsymbol{u})^{\textrm{T}})\boldsymbol{\cdot}\boldsymbol{j}\,\textrm{d} S, \end{gather}$$
where 
$\boldsymbol {n}$ is the unit vector normal to the sphere's surface 
$S$ and 
${\mu }$ the dynamic viscosity. Henceforth, we express the stratification drag force normalised by the drag force in the absence of stratification (
$C_D^H$): 
$C_D^N=C_D^S/C_D^H$.
In order to explore the mechanisms governing the drag enhancement by stratification, a force decomposition scheme was applied following Zhang et al. (Reference Zhang, Mercier and Magnaudet2019). The velocity and pressure fields were decomposed into the form 
$\boldsymbol {u} =\boldsymbol {u}_H + \boldsymbol {u}_{\rho }$ and 
$p = p_H + p_{\rho }$ with
$$\begin{gather} \frac{\partial\boldsymbol{u}_H}{\partial t}+\frac{1}{\epsilon}\boldsymbol{u}_H\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u}_H={-}\epsilon\boldsymbol{\nabla}p_H+\frac{1}{\textit{Re}}\nabla^2 \boldsymbol{u}_H-A\frac{\epsilon}{\textit{Da}\textit{Re}} \boldsymbol{u}_H, \end{gather}$$
$$\begin{gather}\frac{\partial\boldsymbol{u_{\rho}}}{\partial t}+\frac{1}{\epsilon}(\boldsymbol{u_{\rho}}+\boldsymbol{u}_H) \boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}_{\rho} + \frac{1}{\epsilon}\boldsymbol{u}_{\rho}\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u}_H ={-}\epsilon\boldsymbol{\nabla}p_{\rho}+\frac{1}{\textit{Re}}\nabla^2 \boldsymbol{u}_{\rho}-\frac{\rho'}{\textit{Fr}^2}\boldsymbol{j}-A\frac{\epsilon}{\textit{Da}\textit{Re}} \boldsymbol{u}_{\rho}, \end{gather}$$
where (
$\boldsymbol {u}_H$, 
$p_H$) are the velocity and pressure in the homogeneous fluid and (
$\boldsymbol {u}_{\rho }$, 
$p_{\rho }$) are density-induced contribution in stratified fluid. The buoyancy-induced pressure disturbance in the ambient fluid can further be decomposed into the form 
$p_{\rho } = p_{\rho \rho } + p_{\rho u} + p_{\rho \omega }$, where 
$p_{\rho \omega }$ obeys a Laplace equation, while the remaining two contributions satisfy
$$\begin{gather} \nabla^2 p_{\rho\rho} ={-}\frac{1}{\textit{Fr}^2}\boldsymbol{\nabla}\rho'\boldsymbol{\cdot}\boldsymbol{j}, \end{gather}$$
$$\begin{gather}\nabla^2 p_{\rho u} ={-} \boldsymbol{\nabla}\boldsymbol{\cdot}\{(\boldsymbol{u}_{\rho} + \boldsymbol{u}_H)\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}_{\rho} +\boldsymbol{u}_H\boldsymbol{\cdot}\boldsymbol{\nabla}u_{\rho} + \boldsymbol{u_{\rho}} \boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}_H\}, \end{gather}$$
with 
$p_{\rho \rho }$ and 
$p_{\rho u}$ vanishing in the far field. These two Poisson equations were solved by a lattice Boltzmann method (Chai & Shi Reference Chai and Shi2008). The external stratification-induced forces can then be expressed as
\begin{align} C^N_D - 1 &= \underbrace{-\frac{1}{\dfrac{1}{2}\rho W^2{\rm \pi} a^2 C^H_D}\int_S p_{\rho\rho}\boldsymbol{j}\,\mathrm{d}S}_{F_{\rho\rho}} \underbrace{-\frac{1}{\frac{1}{2}\rho W^2{\rm \pi} a^2 C^H_D}\int_S p_{\rho u}\boldsymbol{j}\,\mathrm{d}S}_{F_{\rho u}} \ldots \nonumber\\ &\quad {}\underbrace{+\frac{1}{\dfrac{1}{2}\rho W^2{\rm \pi} a^2 C^H_D}\int_S \boldsymbol{T}_{\rho\omega}\boldsymbol{\cdot}\boldsymbol{j}\,\mathrm{d}S}_{F_{\rho\omega}}, \end{align}
where 
$\boldsymbol {T}_{\rho \omega }=-p_{\rho \omega }+ \mu ((\boldsymbol {\nabla }\boldsymbol {u}_p)+(\boldsymbol {\nabla }\boldsymbol {u}_p)^{\textrm {T}})$ is the vorticity-induced stress tensor. The contribution 
$F_{\rho \rho }$ is an additional Archimedes-like force due to the deflection of the isopycnals, 
$F_{\rho u}$ is an inertial force associated with 
$\boldsymbol {u}_p$ and 
$F_{\rho \omega }$ results from the vorticity-induced baroclinic torque (vortical force). The contributions of 
$F_{\rho u}$ and 
$F_{\rho \rho }$ are directly calculated by solving the Poisson equations (2.7) and (2.8); 
$F_{\rho \omega }$ is subsequently determined as 
$F_{\rho \omega } = C_D^N-1-F_{\rho u}-F_{\rho \rho }$.
The numerical method was verified in three steps. The first was confirming that the homogeneous drag coefficients (in the absence of stratification) of solid as well as porous and impermeable spheres converged towards common values (figures 3a and 17b). The homogeneous drag coefficients were found to closely resemble the empirical relation (White Reference White2005)
\begin{equation} C_D^H=\frac{12}{\textit{Re}}+\frac{6}{1+\sqrt{2\textit{Re}}}+0.4. \end{equation}
Second, the implementation of stratification was tested by comparing our results for solid spheres with earlier quantitative results on stratification drag enhancement (Zhang et al. Reference Zhang, Mercier and Magnaudet2019) (figure 3b) showing excellent agreement for the normalised stratification drag 
$C^N_D$, as well as the two force components 
$F_{\rho \rho }$ and 
$F_{\rho \omega }$ at 
$\textit {Re}=0.05$. Results of additional tests on the convergence criteria for various 
$\textit {Re}$ and 
$\textit {Fr}$ are shown in figure 17.
Figure 3. Model validation based on comparison with existing literature. (a) Reynolds number dependence of the homogeneous drag coefficient 
$C_D^H$ of solid and porous spheres as compared with the empirical relation (2.10) of White (Reference White2005) (see figure 17 for logarithmic representation). (b) The normalised drag coefficient 
$C_D^N$ of solid spheres settling in linear stratification as a function of the Froude number 
$\textit {Fr}$ in comparison with direct numerical simulations obtained by Zhang et al. (Reference Zhang, Mercier and Magnaudet2019) for 
$Re = 0.05$. Here 
$F$ refers to stratification-induced forces, and 
$F_{\rho \rho }$, 
$F_{\rho \omega }$ and 
$F_{\rho u}$ are the force components of 
$C_D^N$. See text for further information.
Third, the geometries of the modelled and experimentally determined density perturbations were compared (figure 15a,b). The overall pattern qualitatively confirms earlier results based on solid spheres (Yick et al. Reference Yick, Stocker, Peakock and Torres2009; Zhang et al. Reference Zhang, Mercier and Magnaudet2019). In the vicinity of the particle, lighter fluid is entrained, forming a density boundary layer with strongly compressed isopycnals. At the lee side, a wake is formed which extends to a distance of several radii. For 
$\textit {Re}=0.66$ and 
$\textit {Fr}=1.37$, the boundary layer is narrow, resulting from the dominance of buoyant forces, while for 
$\textit {Re}=0.1$ and 
$\textit {Fr}=0.88$ viscous forces dominate, resulting in a wide wake and wide boundary layer. Overall, the numerical model represents the observed geometry well for non-dimensional density contrasts (compared with the undisturbed linear stratification) above values of 0.1, i.e. densities displaced by more than one-tenth of the particle radius. For weaker density contrasts, deviations become visible for more strongly stratified fluids, which can be attributed to variations in the Schmidt number. As the Schmidt numbers of the experiments were exceeding 
$\textit {Sc}=700$, we subsequently focus on the numerical model and refer to experimental validations whenever possible.
3. Results and discussion
3.1. Flow and concentration fields
The mean adjusted flow fields around porous and solid spheres are characterised by a vertical succession of recirculation regions (figure 4). Stratification brings about a form of vertical confinement of the flow, as the relative buoyancy of a fluid parcel restricts its vertical deflection. The entrainment of lighter fluid yields a horizontal density gradient in the wake of the particle which translates into baroclinic torque. As a result, toroidal vortices are formed, in agreement with analytical solutions for low-Reynolds-number point forces in a stratified fluid (List Reference List1971; Ardekani & Stocker Reference Ardekani and Stocker2010) and the recent analysis of solid-particle settling in density gradients (Zhang et al. Reference Zhang, Mercier and Magnaudet2019).
Figure 4. Lattice Boltzmann simulations presented as streamlines and scaled vertical velocity 
$w/W-1$ (left-hand side of each panel) and isopycnals combined with the scaled density contrast 
$\Delta \rho /(-a\gamma )$ (right-hand side). (a) Solid sphere at 
$\textit {Re}=0.05$, 
$\textit {Fr}=1$. Porous sphere at (b) 
$\textit {Re}=0.05$, 
$\textit {Fr}=1$ (regime R2; cf. Zhang et al. Reference Zhang, Mercier and Magnaudet2019); (c) 
$\textit {Re}=0.5$, 
$\textit {Fr}=1$ (R2); (d) 
$\textit {Re}=0.5$, 
$\textit {Fr}=10$ (R3). The red dots in (a–d) indicate the position of the rear stagnation points. (e–h) Zoomed depictions of (a–d).
The structure of the velocity field is found to be largely the same for porous and impermeable particles. The flow field in the vicinity of the sphere features two distinct regions separated by a rear stagnation point (red dot, figure 4). Entrained fluid surrounding the sphere effectively travels with the sphere, implying the deformation of the isopycnals, while further downstream the isopycnals are detached and relax towards their equilibrium position. In accordance with Zhang et al. (Reference Zhang, Mercier and Magnaudet2019), the flow structure shows a strong dependency on stratification strength and viscosity for large Péclet numbers 
$\textit {Pe}=aW/D=\textit {Sc} \textit {Re}\geq 1$.
The velocity fields contain a first indication of two different regimes which can be identified with two of the stratification regimes R2 and R3 as described by Zhang et al. (Reference Zhang, Mercier and Magnaudet2019) (cf. § 3.4). For faster settling and stronger stratification, the stagnation point approaches the sphere surface (figure 4c), which is consistent with the formation of a rear buoyant jet at larger Reynolds numbers 
$\textit {Re}\geq \textit {O}(10)$ (Torres et al. Reference Torres, Hanazaki, Ochoa, Castillo and van Woert2000; Hanazaki, Kashimoto & Okamura Reference Hanazaki, Kashimoto and Okamura2009; Okino, Akiyama & Hanazaki Reference Okino, Akiyama and Hanazaki2017) corresponding to regime 2, R2, for 
$\textit {Sc}^{-1/2}\ll \textit {Fr} \ll \textit {Re}^{-1}$ (retaining the notation of Zhang et al. (Reference Zhang, Mercier and Magnaudet2019)). For the opposing case of weaker stratification, the stagnation point increasingly recedes from the sphere (figure 4d) and the asymmetry of the velocity field vanishes at 
$\textit {Fr}\geq \textit {O}(100)$ when buoyant forces become negligible and the flow field is largely determined by viscous forces, corresponding to regime 3, R3, 
$\textit {Fr}\gg \textit {Re}^{-1}$ (Zhang et al. Reference Zhang, Mercier and Magnaudet2019).
In the proximity of the sphere surface, distinct differences between porous and solid cases were observed. For the porous case, the stagnation point is consistently closer to the surface for small 
$\textit {Fr}$ and large 
$\textit {Re}$, while as the stagnation point recedes from the surface (
$\textit {Fr}\gg \textit {Re}^{-1}$) these differences vanish (figure 5a). The position of the stagnation point is associated with the vertical velocities in the wake which indicate enhanced ascending velocities for stronger stratification, i.e. small 
$\textit {Fr}$, and large 
$\textit {Re}$.
Figure 5. (a) Modelled vertical velocity profiles in the centre of the wake (
$r/a=0$; solid and dashed lines depict porous and solid sphere values, respectively). (b) Density perturbation in the centre of the wake (
$r=0$) and (c) the corresponding density perturbation in the boundary layer at the equator (
$z/a=0$, cf. figure 4a–c) and inside the porous particle. The shaded areas indicate the particle interior.
The differences in velocities in the wake can be associated with an enhanced density contrast of the boundary layer fluid (figure 5b). Diffusive exchange with the interstitial pore fluid at the particle surface alters the density contrast of the density boundary layer while the boundary layer thickness remains largely unchanged when compared to solid particles (figure 5c). The thickness of the density boundary layer shows a dependency on 
$\textit {Fr}$ and 
$\textit {Re}$. The equatorial concentration profiles indicate substantial advective exchange of boundary layer fluid visible as slight kinks in the concentration profile at the sphere surface in figure 5(c) which implies lower-density fluid being fed into the wake altering the density contrast of the wake fluid, too.
The diffusive exchange between interstitial and boundary layer fluid controls the mass adaptation of the particle. Thus, the diffusive exchange also impacts the external density and velocity field, altering the stratification drag enhancement. The overall mass adaptation of the particle itself is limited by the viscous turnover of boundary layer fluid. To rationalise the effect of the diffusive exchange of interstitial and boundary layer fluid, we derive scaling laws for the density boundary layer thickness, the density contrast of the density boundary layer and the interstitial density in the following sections.
3.2. The density boundary layer
In the case of porous spheres, the diffusive exchange of interstitial and ambient fluid potentially alters the fluid density in the direct vicinity of a particle. We consider this altered layer as a density boundary layer, which can enhance the Archimedes-like contribution to stratification drag enhancement. The thickness of the density boundary layer, 
$\delta$, buffers the diffusive exchange between the external density field and the interstitial liquid.
Here, 
$\delta$ was estimated as the azimuthal stoss-side average of the distance between the porous sphere surface and the point at which the normalised density contrast increased to 95 % of that of the ambient fluid. The 
$95\,\%$ threshold is an operational definition and changing this threshold will affect the determined volume. For example, a 
$99\,\%$ threshold results in an increase of 
$\delta$ by 
$60\,\%$ and an 
$80\,\%$ threshold will decrease 
$\delta$ by 
$30\,\%$. However, while the magnitude is sensitive to the threshold, we did not observe an effect on the scaling slopes. In accordance with the flow field, we observed a regime separation at 
$\textit {Fr} \approx \textit {Re}^{-1}$ (figure 6a). The best fit to our results (figure 6b,c) yields
\begin{equation} \frac{\delta_{R2}}{a}=0.45\textit{Re}^{{-}0.66}\textit{Fr}^{0.61}\end{equation}and
\begin{equation} \frac{\delta_{R3}}{a}=1.1\textit{Re}^{{-}0.71}\textit{Fr}^{0.14}. \end{equation}
Figure 6. Scaling of the boundary layer thickness 
$\delta /a$ for 
$Re = 0.5$ and 
$Re = 0.05$ (a) and as a generalised empirical relation for regime R2 (b) and regime R3 (c). The boundary layer density contrast 
$\Delta \rho _\delta$ for 
$Re = 0.5$ and 
$Re = 0.05$ (d) and as a generalised empirical relation for regime R2 (e) and regime R3 (f). (g) The scaling relationships for 
$\Delta \rho _p/(-a\gamma )$ within the two regimes for 
$Re=0.5$ and 
$Re=0.05$ in log–log representation and the generalised empirical relation for the density contrast as a function of 
$\textit {Re}$ and 
$\textit {Fr}$ for regime R2 (h) and regime R3 (i). Blue and orange coloured symbols refer to regime 2 and regime 3, respectively (cf. figure 10).
In regime R2, when stratification is strong, the impact of diffusion vanishes and the density boundary layer thickness can be rationalised by the balance of viscous and buoyant forces through the natural length scale 
$\delta \sim (\nu /N)^{1/2}$, which implies 
$\delta /a\sim (\textit {Fr}/\textit {Re})^{1/2}$ (Yick et al. Reference Yick, Stocker, Peakock and Torres2009). In regime R3 the influence of 
$\textit {Fr}$, i.e. buoyant force, becomes negligible and the boundary layer thickness is largely defined by the interaction of inertial and viscous forces, represented through 
$\textit {Re}$. Overall, experiments and numerical results match well in both regimes. The density boundary layer fluid volume is found to be largely independent of diffusive exchange with the particle pore fluid – as indicated by isopycnals for porous and solid spheres (figure 5c,d) as well as the scaling in (3.1) and (3.2).
The density contrast of the boundary layer is, however, determined by the interaction of the external forcing through the ambient fluid and diffusive equilibration with the interstitial fluid. The density contrast of the boundary layer fluid with respect to an undisturbed stratified fluid 
$\Delta \rho _\delta / (-a \gamma )$ was evaluated as the average within the shell with width 
$\delta / a$ surrounding the particle (figure 6d–f). The best collapse (in terms of least square errors) was found for
\begin{equation} \frac{\Delta\rho_{\delta,R2}}{-a\gamma}=\frac{(\rho_\delta-\rho_B)}{-a\gamma }={-}39\textit{Re}^{0.43} \textit{Fr}^{0.26}\end{equation}and
\begin{equation} \frac{\Delta\rho_{\delta,R3}}{-a\gamma}={-}40\textit{Re}^{0.45}\textit{Fr}^{0.06}, \end{equation}
where 
$\rho _\delta$ represents the density of the density boundary layer and 
$\Delta \rho _\delta$ its excess density with respect to the ambient fluid. The simulations seem to overestimate the experimentally derived density contrast of the density boundary layer. The mismatch likely results from not completely negligible permeability in the experiments (cf. figure 14) and the fact that the experiments did not reach perfect stationarity. Further away from the sphere, the concentration fields demonstrate little effect of porosity when compared with those of solid spheres, where the boundary layer fluid passes into a slender concentration wake structure downstream from the sphere, extending up to O(
$10a$) in length (figure 4b,c).
3.3. Interstitial mass adaptation
The interstitial mass adaptation is controlled by the density boundary layer, i.e. the mass adaptation of the particle generally depends on the viscous and buoyant forces in the external flow field. The pore fluid excess density 
$\Delta \rho _p$ was evaluated as the spatial average of the interstitial density contrast:
\begin{equation} \frac{\Delta\rho_p}{-a\gamma}={-}\frac{3}{2a^4\gamma}\int_{v} \Delta\rho\,\textrm{d}{v},\end{equation}
where 
$\Delta \rho =\rho _p-\rho _B$. The interstitial pore water density was found to increase as 
$\Delta \rho _p\sim \textit {Pe}$ (data not shown). The remaining Froude number dependence can be ascribed to the fact that viscous exchange in the density boundary layer mitigates the exchange between the interstitial and the ambient fluid. The boundary layer density 
$\rho _\delta$ is reduced with respect to an undisturbed density field 
$\rho$.
Following Fick's first law, the dimensional diffusive exchange between the particle and the surrounding density boundary layer can be approximated via a shell model:
\begin{equation} \frac{\partial\rho_p}{\partial t}\sim 3\frac{D}{\epsilon a^2}(\rho_{\delta}-\rho_p), \end{equation}
where 
$\rho _{\delta }$ can be exchanged by our scaling relationships (3.3) and (3.4) for the two regimes (figure 6d–f):
\begin{equation} \frac{\partial\rho_p}{\partial t} \sim 3\frac{D}{\epsilon a^2} \left(\alpha\textit{Re}^{n}\textit{Fr}^{m} - \rho_p-\rho_B \right),\end{equation}
with 
$\alpha = 39 a \gamma$, 
$n=0.43$, 
$m=0.26$ in regime R2 and 
$\alpha = 40 a \gamma$, 
$n=0.45$, 
$m\approx 0$ in R3.
Settling is steady when the rate of change of the boundary layer density is equivalent to that of the external field experienced by the particle:
\begin{equation} \frac{\partial\rho_p}{\partial t}=\frac{\partial\rho_\delta}{\partial t}=\gamma W. \end{equation} Under these conditions a direct relationship between 
${\Delta \rho _p}/{-a\gamma }$ and 
${\Delta \rho _\delta }/{-a\gamma }$ can be derived. The transport equation (2.2) inside the porous particle reduces to the non-dimensional diffusion equation:
\begin{equation} \nabla^2\left(\frac{\Delta\rho}{-a\gamma}\right) = \textit{Re}\textit{Sc} . \end{equation}The density flux through the surface of the particle is balanced by the changing external field:
\begin{equation} 4{\rm \pi} a^2 J=V_p\epsilon \textit{Re}\textit{Sc}, \end{equation}
where 
$J$ is the flux normal to the particle surface and 
$V_p$ is the volume of the particle. Based on the integration of (3.9) with the boundary conditions in (3.10) one can calculate the density distribution inside the particle:
\begin{equation} \frac{\Delta\rho(r)}{-a\gamma} = \frac16 \epsilon r^2 \textit{Re}\textit{Sc} + c, \end{equation}
where 
$c$ is the integration constant. Based on (3.11) the averaged density on the surface is 
${\Delta \rho _{\delta }}/{(-a\gamma)}= \frac {1}{6}\epsilon \textit {Re}\textit {Sc} + c$ and the averaged density inside the particle is 
${\Delta \rho _{p}}/{-a\gamma }= \frac {1}{10}\epsilon \textit {Re}\textit {Sc} + c$. The difference of the averaged density of the interstitial pore water and the density in the boundary layer then is 
${\Delta \rho _{p}}/{-a\gamma }=-44Re+{\Delta \rho _{\delta }}/{-a\gamma }$ (for 
$\textit {Sc}=700$), which implies that the scaling relationship 
$\Delta \rho _{\delta }$ also applies for 
$\Delta \rho _{p}$ if shifted by 
$44Re$. Indeed, best fit to our simulations gives close results to the analytical solution (figure 6g–i):
$$\begin{gather} \frac{\Delta\rho_{p,{R2}}}{-a\gamma}=\frac{\rho_p-\rho_B}{-a\gamma}\approx{-}39 (\textit{Re}+\textit{Re}^{0.43}\textit{Fr}^{0.26}), \end{gather}$$
$$\begin{gather}\frac{\Delta\rho_{p,{R3}}}{-a\gamma}= \frac{\rho_p-\rho_B}{-a\gamma}\approx{-}40 (\textit{Re}+\textit{Re}^{0.45}), \end{gather}$$
for the non-dimensional interstitial fluid density contrast with respect to the ambient fluid. In regime 3, the effects of 
$\textit {Fr}$ become vanishingly small which, however, does not imply that effects of stratification are negligible. The delay of mass adaptation depends on the boundary layer thickness determined through the viscous forces.
3.4. Drag enhancement versus mass adaptation
The settling of porous particles in stratification is reduced compared to settling in an unstratified fluid by both the delayed mass adaptation of the interstitial pore fluid and external buoyancy-induced forces represented through the stratification drag coefficient 
$C_D^S$.
Note that we attribute all external forces, the additional Archimedes-like forces of the density boundary layer (
$F_{\rho \rho }$), the inertial force (
$F_{\rho u}$) as well as forces due to vorticity (
$F_{\rho \omega }$) to drag while the term density adaptation refers to the density adjustment of the pore fluid inside the particle only. We consider a scaling law for the normalised stratification drag 
$C_D^N = C_D^S / C_D^H$ depending on 
$\textit {Re}$ and 
$\textit {Fr}$ (figure 7d), with the best fit yielding
\begin{equation} C_{D,R2}^N-1=3.25\textit{Re}^{0.51}\textit{Fr}^{{-}1.1}\end{equation}and
\begin{equation} \mathrm C_{D,R3}^N-1=2.51\textit{Re}^{0.1}\textit{Fr}^{{-}1.45}. \end{equation}
In R2, the drag coefficient is close to that for the solid case 
$C_{D}\sim \textit {Ri}^{0.5}$ confirming the results of Yick et al. (Reference Yick, Stocker, Peakock and Torres2009). In R3 the increase in drag is largely independent of 
$\textit {Re}$ and is mainly controlled by 
$\textit {Fr}$, i.e. the ratio of inertia and stratification strength.
Figure 7. Scaling relationship for the normalised drag coefficient of a porous particle settling in stratification. (a) The normalised drag coefficient 
$C_D^N$ as function of 
$\textit {Fr}$ for 
$\textit {Re}=0.05$ and 
$\textit {Re}=0.5$. The global scaling of 
$C_D^N$ as a function of 
$\textit {Re}$ and 
$\textit {Fr}$ for regime 2 (b) and regime 3 (c). The experimental results are shown as circles with error bars. Blue and orange coloured symbols refer to regime 2 and regime 3, respectively (cf. figure 17).
In the case of porous particles, the gradients in the density boundary layer are steepened compared to the case of solid particles. To investigate the effect of these steepened gradients on the exerted external forces on the porous particles, we applied a force decomposition ((2.9) and figure 8). In regime R2, stronger gradients in the density boundary layer and in the wake are found to translate into an additional augmentation of 
$F_{\rho \rho }$ and 
$F_{\rho \omega }$ as compared to the case of solid particles (figure 8a,b).
Figure 8. The normalised drag coefficient of solid and porous spheres settling in linear stratification 
$C_D^N$ as a function of the Froude number 
$\textit {Fr}$ for (a) solid, 
$Re = 0.05$, (b) porous, 
$Re = 0.05$, (c) solid, 
$Re = 0.5$ and (d) porous, 
$Re = 0.5$. Here 
$F$ refers to stratification-induced forces. Dots indicate at which Fr numbers model runs were performed.
To quantify within the two regimes, we found the vortical force to scale as
\begin{equation} F_{\rho\omega,R2}\sim (\textit{Re}\textit{Fr})^{{-}0.69} \end{equation}and
\begin{equation} F_{\rho\omega,R3} \sim \textit{Re}^{{-}0.62}\textit{Fr}^{{-}1}, \end{equation}
which is very similar to the values of a solid particle which Zhang et al. (Reference Zhang, Mercier and Magnaudet2019) found to scale as 
$F_{\rho \omega,R2}\sim (\textit {Re}\textit {Fr})^{-0.67}$ and 
$F_{\rho \omega,R3} \sim (\textit {Re}\textit {Fr})^{-1}$. Differences in regime 3 are likely to be associated with a smooth transition from the intermediate- to low-Reynolds-number regimes (Zhang et al. (Reference Zhang, Mercier and Magnaudet2019) found 
$F_{\rho \omega,R3} \sim \textit {Re}^{-0.5}\textit {Fr}^{-1}$ for intermediate Reynolds numbers). Overall the external forces associated with the toroidal eddies follow similar scaling relationships between porous and solid particles. However, and with exception for very small Froude number at 
$Re=0.5$, we found an increase in the magnitude of 
$F_{\rho \omega }$ in the case of porous particles. This effect is associated with the strongly increased density contrast in the wake of the particles fed by the interstitial pore water (figure 6b) which results in an increased vortex production in the external field.
In the case of solid particles, the Archimedes-like force 
$F_{\rho \rho }$ was found to play a secondary role for stratification drag enhancement (Zhang et al. Reference Zhang, Mercier and Magnaudet2019). We found 
$F_{\rho \rho }$ to scale as
\begin{equation} F_{\rho\rho,R2}\sim\textit{Re}^{{-}0.47}\textit{Fr}^{{-}1.7} \end{equation}and
\begin{equation} F_{\rho\rho,R3} \sim\textit{Re}^{{-}0.59}\textit{Fr}^{{-}2}, \end{equation}
contrasting with the solid-particle case where, in both regimes, 
$F_{\rho \rho } \sim \textit {Fr}^{-2}$ (Zhang et al. Reference Zhang, Mercier and Magnaudet2019). Differences are mainly associated with the additional 
$\textit {Re}$ dependency that indicates the contribution of enhanced density contrast in the boundary layer fluid (
$\Delta \rho _p \sim \textit {Re}$ for 
$\textit {Sc}=\textrm {const.}$; see also figure 6). The mass force exerted by the density contrast between pore water and ambient fluid is
\begin{equation} F_{m}\sim\frac{4/3 {\rm \pi}a^3 g\Delta\rho_p}{0.5C_D^H\rho W^2{\rm \pi} a^2}\sim \frac{1}{\textit{Fr}^2C_D^H}\frac{\Delta\rho_p}{({-}a\gamma)}. \end{equation}
Mass force 
$F_m$ is reduced due the light fluid travelling in the interior of the porous particle (figures 5c,d and 6). The effect of stratification on the mass force 
$F_m$ is specific for porous particles and does not directly result in a drag enhancement, but decreases the settling velocity through the reduction of the excess density (see also § 4). At low 
$\textit {Re}$ and weak stratification (
$\textit {Fr}>1$) the contribution of 
$F_{m}$ is negligible (figure 9a). However, at larger 
$\textit {Re}$ and stronger stratification (
$\textit {Fr}<1$) the density contrast is substantially increased and 
$F_{m}$ even exceeds the contribution of the external forces 
$C_D^N$ (figure 9b). In direct comparison between solid and porous particles, we found for 
$\textit {Re} = 0.05$ the buoyancy-induced Archimedes-like and mass forces to increase by a factor of 2 to 7 and 8 to 44 for 
$\textit {Re} = 0.5$ (figure 8b). In conclusion, in the case of porous particles the buoyancy forces induced by the water volume travelling with the porous spheres, 
$F_{\rho \rho }$ and 
$F_{m}$, are non-negligible and comprise a range similar to that of 
$F_{\rho \omega }$.
Figure 9. Normalised mass force of the interstitial pore space 
$F_{m}$ for (a) 
$\textit {Re} = 0.05$ and (b) 
$\textit {Re} = 0.5$. The dashed line indicates the separation between the two regimes, R2 and R3. Here 
$F$ refers to stratification-induced forces.
4. Implications for particle settling in stratified marine environments
Under in situ oceanic conditions, it is generally difficult to quantify all variables that determine the sinking velocity of heterogeneously composed marine aggregates. This has impeded direct in situ analysis of the effects of stratification on particle settling. Nevertheless, a limited number of experimental or indirect in situ measurements have suggested an effect of stratification on the sinking velocity of large marine aggregates (
${>}500\,\mathrm {\mu }\textrm {m}$), known as ‘marine snow’ (MacIntyre et al. Reference MacIntyre, Alldredge and Gotschalk1995; Alldredge et al. Reference Alldredge, Cowles, McIntyre, Rins, Donaghay, Greenlaw, Holliday, Dekshenieks, Sullivan and Zaneveld2002; Prairie et al. Reference Prairie, Ziervogel, Camassa, McLaughlin, White, Dewald and Arnosti2015).
To quantify the potential effect of stratification on settling of marine aggregates, we estimate the retention of porous 95 % and impermeable spheres in linear stratifications of varying strength (from 
$N=0.01$ to 
$1\,\textrm {s}^{-1}$). We consider the balance of buoyancy and drag force:
\begin{equation} \frac{4}{3}{\rm \pi} a^3 \rho_a\frac{\textrm{d} W}{\textrm{d} t}=(\rho_a-\rho_B) \frac{4}{3}{\rm \pi} a^3g + \frac{\rho}{2}C_D^{S}{\rm \pi} a^2W^2, \end{equation}
where the Basset history force was neglected. The particle density is defined as 
$\rho _a = (1-\epsilon )\rho _s + \epsilon (\rho _B+\Delta \rho _p)$, with 
$\rho _s$ the hydrated excess density of the solids. The interstitial density and stratification drag enhancement are calculated based on the derived scaling-law equations (3.12), (3.13), (3.14) and (3.15) (see also Appendix C). To determine the effect of the boundary layer fluid in buffering the diffusive exchange of the particle with the ambient stratified fluid, we further applied (3.7) assuming instantaneous homogenisation of the pore water. The parameters are chosen to resemble those of marine aggregates in stratified oceanic environments. Marine aggregates feature hydrated excess densities of the solids (
$\rho _s-\rho _B$) of about 
$\textit {O}(0.01)$ to 
$\textit {O}(100\,\textrm {kg}\,\textrm {m}^{-3})$ (Alldredge & Gotschalk Reference Alldredge and Gotschalk1988), porosities 
${\gg }0.95$ (Alldredge & Gotschalk Reference Alldredge and Gotschalk1988) and sinking velocities of about 4 to 
$416\,\textrm {m d}^{-1}$ (e.g. Iversen & Lampitt Reference Iversen and Lampitt2020), which leads to an intermediate Reynolds regime of 
$0.05\leq \textit {Re}\leq 20$. Considering typical buoyancy frequencies in the ocean and salt lakes (
$N=0.001$ to 
$0.1\,\textrm {s}^{-1}$) (Boehrer & Schultze Reference Boehrer and Schultze2008; Geyer et al. Reference Geyer, Scully and Ralston2008) as well as strong stratification in estuaries and fjords (
$N=0.1$ to 
$0.3\,\textrm {s}^{-1}$) (Geyer et al. Reference Geyer, Scully and Ralston2008), Froude numbers typically fall in the range 
$1\leq \textit {Fr}\leq 100$, and in extreme cases as low as O(0.01). These parameter ranges imply that natural marine aggregates fall into both the R2 and R3 regimes. Further, it is important to notice that our parameterisation represents the flow regime of particles within stratification. This implies, for example, that a porous particle featuring a particular 
$\textit {Re}$ in linear stratification (and hence 
$W$ being constant) will have a higher 
$\textit {Re}$ under homogeneous fluid conditions (assuming 
$\nu$ to be identical). Our scaling laws can thus be applied to a wider 
$\textit {Re}$ range if settling velocities and associated flow regimes are determined in homogeneous fluids (see also Appendix D).
The concurrent processes of delayed mass adaptation and drag enhancement increase the retention time of particles (figure 10a). This increase strongly depends on particle size and buoyancy frequency (figure 10b). Generally, stratification-induced processes start to influence the retention time of particles for 
$N>1\times 10^{-4}\,\textrm {s}^{-1}$. For weak stratification (
$N\leq 0.02\,\textrm {s}^{-1}$) or small particles (
$a<200\,\mathrm {\mu }\textrm {m}$), the increase in retention time is below 10 %, mainly due to drag enhancement (see figure 11). By contrast, mass adaptation dominates the retention enhancement for particles 
${>}0.5\,\textrm {mm}$ and reduces their sinking velocities by 10 % to 100 %. Under strong stratification (
$N>0.2\,\textrm {s}^{-1}$), the sinking velocity of large particles (
$a>1\,\textrm {mm}$) thus can even drop to almost zero, which contrasts with the typical expected increasing sinking velocity with size in non-stratified fluids. For frequently measured intermediate stratification strengths (
$N=0.05\,\textrm {s}^{-1}$) and abundant particle sizes (
$a=1\,\textrm {mm}$) the settling velocity is expected to be reduced by 5 %–15 % due to the combined effect of drag enhancement and delayed mass adaptation.
Figure 10. (a) Time trajectories of a porous sphere (
$a=1\,\textrm {mm}$, 
$\rho _s-\rho _B=10\,\textrm {kg}\,\textrm {m}^{-3}$) settling in a state of instantaneous mass adaptation (dashed curve), delayed mass adaptation (blue) and delayed mass adaptation in the presence of a density boundary layer (orange), in a stratification of 
$N=0.2\,\textrm {s}^{-1}$. For comparison, the trajectory of a solid sphere of equal initial excess density is depicted (black curve). (b) Size dependence of the normalised retention time of porous spheres settling at 
$\textit {Pe}\geq \textit {O}(10)$ for weak (
$N=0.02\,\textrm {s}^{-1}$) and strong (
$N=0.2\,\textrm {s}^{-1}$) stratification. Filled areas indicate the contribution to retention time by delayed mass adaptation (blue) and buffered delayed mass adaptation (orange). The retention time 
$t$ was calculated based on the travelling distance along one particle diameter and scaled by that of a sphere with instantaneous diffusive equilibration 
$t_0$.
Figure 11. Contour lines indicating the increased retention (
$t/t_0-1$) in the phase space of radius 
$a$ and buoyancy frequency 
$N$. Contour line labels represent exponents of basis 10. (a) Solid hydrated density is adjusted to resemble solid particles such as faecal pellets and coccolith aggregates (
$\rho _s-\rho _B=100\,\textrm {kg}\,\textrm {m}^{-3}$). (b) Solid hydrated density is adjusted to resemble light diatom aggregates (
$\rho _s-\rho _B=10\,\textrm {kg}\,\textrm {m}^{-3}$).
In contrast to diffusion-limited retention in two-layer stratification, i.e. sharp pycnoclines (Panah et al. Reference Panah, Blanchette and Khatri2017), the increase in retention time in linear stratification is not driven by a transient diffusive exchange, but instead by the equilibrium of diffusive exchange of pore water and the density change in the surrounding water column, and hence 
$W=\textrm {const.}$ (see also extended discussion in Appendix E). Nevertheless, we found the density boundary layer to buffer diffusive exchange in accordance with studies focusing on two-layer stratification (Camassa et al. Reference Camassa, Khatri, McLaughlin, Prairie, White and Yu2013). Increased retention of marine aggregates through slower settling, however, can occur over several to tens of metres due to the large extension of pycnocline layers under typical oceanic conditions, which are here considered as linear stratification. In addition, the settling and retention are strongly modulated by composition of aggregates (figure 11a,b). Composition of marine aggregates affects 
$\rho _s$, thus 
$\textit {Re}$ and 
$\textit {Fr}$, and hence the location of the regime transition (between R2 and R3) in the phase space (i.e. size–buoyancy frequency). Lighter marine aggregates (e.g. made from diatoms; Alldredge Reference Alldredge1998) settle slowly, which shifts the parameter range to a lower-Reynolds-number regime reducing the effects of delayed mass adaptation and drag enhancement compared with heavier particles (e.g. faecal pellets or coccolith aggregates; Alldredge Reference Alldredge1998). Thus, stratification will lead to a stronger reduction in settling velocity for diatom-composed aggregates than for heavy particles such as coccolith aggregates and faecal pellets.
The increased retention time of marine aggregates and their ongoing microbial degradation may lead to a substantial decrease in carbon export to the deep ocean. To estimate such potential changes in carbon sequestration due to projected increase of ocean stratification (Bopp et al. Reference Bopp2013) requires the representation of stratification-dependent sinking velocity in Earth system models. For this purpose, biogeochemical models need to adequately capture the size spectrum and excess density of marine aggregates (e.g. Maerz et al. Reference Maerz, Six, Stemmler, Ahmerkamp and Ilyina2020). This enables one to apply (3.12), (3.13) and (3.14), (3.15), respectively, to allow for studying the effect of stratification-dependent sinking on biogeochemical cycles in regional- to global-scale models. In summary, the presented results and mechanistic processes in our study demonstrate the potential importance of ocean stratification for the settling velocity of marine aggregates and thus for the strength of the biological carbon pump.
5. Conclusions
Experiments and numerical simulations revealed a pronounced impact of linear stratification on settling velocity of highly porous, impermeable spheres. In the Reynolds and Froude number regime characteristic for marine environments, delayed mass adaptation of the interstitial fluid and stratification drag enhancement concurred in slowing the sinking of porous spheres. The density boundary layer was found to mitigate the diffusive exchange of stratifying agent, serving as a buffer layer around the spheres. The interrelation of the density boundary layer, the diffusive exchange of interstitial and boundary layer fluid and the settling velocity can be rationalised in scaling laws within two stratification regimes. The key parameters are the non-dimensional density contrast (3.12) and (3.13) and the normalised drag coefficient (3.14) and (3.15). Applying the scaling laws in a rather simple model revealed that the retention time strongly depends on particle size and buoyancy frequency. According to our results and for the case of marine particles, settling of smaller particles (
$a < 0.2\,\textrm {mm}$) is mainly affected by enhanced drag that increases the retention by up to 10 %, while for larger particle sizes (
$a > 0.5\,\textrm {mm}$) delayed mass adaptation reduces settling velocity by 10 % and up to 
${\sim }100\,\%$. The results emphasise the potential importance of stratification for vertical carbon fluxes. The relationships derived here can be directly incorporated in biogeochemical models, which perspectively enables the quantitative assessment of the impact of stratification-dependent sinking of marine aggregates on carbon sequestration.
Acknowledgements
We thank R. Naisbit for discussion. We thank J. Magnaudet and J. Zhang for providing the force decomposition results for solid spheres settling in stratification. We further thank M. Mercier and an anonymous reviewer for their constructive comments on the manuscript that helped to improve the presentation of the work.
Funding
The authors acknowledge funding from the Max Planck Society (MPG) for the Multiscale Approach on the Role of Marine Aggregates (MARMA) project. B.L. acknowledges additional funding from the Helmholtz Association (Alfred Wegener Institute Helmholtz Centre for Polar and Marine Research). R.S. acknowledges funding from the Simons Foundation through the Principles of Microbial Ecosystems (PriME) Collaborative (grant 542395).
Declaration of interests
The authors report no conflict of interest.
Author contributions
S. Ahmerkamp, B. Liu and K. Kindler contributed equally to this work.
Appendix A. Reconstructing the density fields
For this investigation, an inverse tomographic algorithm was implemented to reconstruct the three-dimensional density distribution following Liu, Merzkirch & Oberste-Lehn (Reference Liu, Merzkirch and Oberste-Lehn1988), based on the speckle displacement. The reconstruction of finite-density perturbations based on the schlieren method has been used in various applications (Kindler et al. Reference Kindler, Goldhahn, Leopold and Raffel2007; Yick et al. Reference Yick, Stocker and Peakock2007). Assuming a light ray propagating perpendicular to a perturbed density field (
$L_p$), it is deflected around the particle due to changes in the refractive index field. In mathematical terms, the deflection reads
\begin{equation} \tan{\psi_{(x,z)}} = \int\frac{1}{n(x,z)} \frac{\partial n}{\partial (x,z)} {\textrm{d}y}, \end{equation}
where 
$\psi$ is the deflection angle, 
$n$ the index of refraction and 
$y$ the line of sight. Assuming that the fluctuations of the refractive index field are much smaller than the mean of the refractive index (
$n_0$), the deflection angle is related to the Radon transform of the density perturbation (
$\rho '$) by a constant 
$a$:
\begin{equation} \psi_{(x,z)} \approx a \int \frac{\partial \rho'}{\partial (x,z)}{\textrm{d}y}. \end{equation} This deflection can be visualised using the background oriented schlieren technique. Assuming a slice of the density field and rotational symmetry, polar coordinates may be introduced: 
$x' = x \cos \theta + y \sin \theta$, 
$y' = - x \sin \theta + y \cos \theta$. The Fourier transformation of the deflection angles and the density field reads
\begin{equation} \left.\begin{gathered} \mathscr{F} (\psi) \{k_{x'},\theta \} = \int_{-\infty}^{\infty}\psi (x',\theta) \exp({-2{\rm \pi} k_{x'}x'})\,{\textrm{d}x}', \\ \mathscr{F} (\rho') \{k_{x'},\theta \}= a \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \rho' \exp({-2{\rm \pi} k_{x'}x'})\,{\textrm{d}x}'\,{\textrm{d}y}', \end{gathered}\right\} \end{equation}
where 
${k_x'}$ denotes the wavenumber in the Fourier domain. From (A1) it follows that 
$\mathscr {F} (\rho ') = {\mathscr {F}(\psi ) }/{\textrm {i} 2{\rm \pi} k_{x'}}$ and calculating the inverse Fourier transform yields
\begin{equation} \rho' (x,y) = a \int_{0}^{\rm \pi} \int_{-\infty}^{\infty} \underbrace{\frac{|k_{x'}|}{ 2 {\rm \pi}k_{x'} i}}_{{=}Q(k_x')} \cdot \mathscr{F} (\psi) \{k_{x'},\theta \} \,\textrm{d} k_{x'} \,\textrm{d}\theta + \rho_B, \end{equation}
where 
$Q(k_{x'})$ is a ramp filter applied in the Fourier domain. Applying the convolution theorem and introducing the undisturbed background density (
$\rho _B$) yields
\begin{equation} \rho (x,y) = a \int_{0}^{\rm \pi} Q(k_x) * \psi (k_{x'},\theta ) \,\textrm{d}\theta + \rho_{B}. \end{equation}
Figure 12. Parameter spaces of recent literature for two-layer stratification (a) and linear stratification (b). The dashed line indicates the most recent in situ measurements performed by Iversen & Lampitt (Reference Iversen and Lampitt2020). Note that there is a difference in the Reynolds number calculations, i.e. for two-layer stratification the Reynolds number is calculated based on the settling in an unstratified fluid while in the cases of linear stratification the Reynolds number is based on the stationary settling velocity in stratification.
Figure 13. (a) Schmidt numbers 
$\textit {Sc}$ for different glycerol–water mixtures based on empirical relationships. Note that density and viscosity for experiments were measured. (b) Reynolds 
$\textit {Re}$ and Froude 
$\textit {Fr}$ numbers of the schlieren experiments relevant to this study.
The inverse Radon transformation and integrations were performed in Matlab (Mathworks 2017b). The density slices were stacked vertically and azimuthally averaged to a two-dimensional plane through the particle centre section. To reduce high-frequency noise, the ramp filter was convoluted with a Hamming window. Examples of the density perturbations are shown in figure 15.
Figure 14. The effects of permeability on the drag coefficient (a) and the density contrast of the particle (b) for 
$\textit {Re}=1$ and 
$\textit {Fr}=1$. The shaded areas indicate the range of permeabilities of the fibre particles which were used for the experiments. Red line indicates the parameterisation of the model (
$\textit {Da}=10^{-12}$).
Figure 15. Density contrast reconstructed from experimental observations (a) for settling porous spheres compared to the lattice Boltzmann simulations (b) for two contrasting cases: (a) 
$\textit {Re}=0.66$, 
$\textit {Fr}=1.37$ and (b) 
$\textit {Re}=0.1$, 
$\textit {Fr}=0.88$. The black lines indicate isopycnals. The dashed lines in (a) indicate scaled density contrasts of 
$\textit {O}(10^{-3})$. Note that the shown axis limits of the simulation are limited to 
$5a$, while the total model domain is 
$52a$.
Appendix B. Lattice Boltzmann method
Because of the axial symmetry of the problem, a quasi-two-dimensional model can be employed using a lattice Boltzmann method for axial symmetric flows (Guo et al. Reference Guo, Han, Shi and Zheng2009; Zheng et al. Reference Zheng, Guo, Shi and Zheng2010a,Reference Zheng, Shi, Guo and Zhengb). The evolution equations for velocity 
$\boldsymbol {u}$ and density perturbations 
$\rho '$ read
$$\begin{gather} {f}_{i}(\boldsymbol{x}+\boldsymbol{e}_{i}\delta_t,t+\delta_t)-{f}_{i}(\boldsymbol{x},t)={-}\frac{1}{\tau_f} (\,{f}_{i}(\boldsymbol{x},t)- {f}_{i}^{(eq)}(\boldsymbol{x},t))+{F}_{i}(\boldsymbol{x},t), \end{gather}$$
$$\begin{gather}{g}_{i}(\boldsymbol{x}+\boldsymbol{e}_{i}\delta_t,t+\delta_t)-{g}_{i}(\boldsymbol{x},t)={-}\frac{1}{\tau_g}({g}_{i}(\boldsymbol{x},t)- {g}_{i}^{(eq)}(\boldsymbol{x},t))+{G}_{i}(\boldsymbol{x},t), \end{gather}$$
where 
${f}_{i}(x,t)$ and 
${g}_{i}(x,t)$ are the single-particle distribution functions at position 
$\boldsymbol {x}$ and time 
$t$ along the direction represented by the subscript 
$i$ for fluid and density perturbation, respectively.
Using the two-dimensional D2Q9 model (Qian, D'Humières & Lallemand Reference Qian, D'Humières and Lallemand1992), the discrete velocities were defined as 
${\boldsymbol {e}_i=(e_{ir},e_{iz}): i=0,1,\ldots,8}$ with 
$\boldsymbol {e}_0=0$, 
$\boldsymbol {e}_1=-\boldsymbol {e}_3=(1,0)$, 
$\boldsymbol {e}_2=-\boldsymbol {e}_4=(0,1)$, 
$\boldsymbol {e}_5=-\boldsymbol {e}_7=(1,1)$ and 
$\boldsymbol {e}_6=-\boldsymbol {e}_8=(-1,1)$. The equilibrium distribution functions read
$$\begin{gather} {f}_{i}^{(eq)} = r{\omega}_{i} \rho\left[1 + \frac{\boldsymbol{e}_{i}\boldsymbol{\cdot}\boldsymbol{u}}{c_s2} +\frac{1}{2\epsilon}\left(\frac{\boldsymbol{e}_{i}\boldsymbol{\cdot}\boldsymbol{u}}{c_s2}\right)2- \frac{1}{2\epsilon}\frac{u^2}{c_s^2}\right], \end{gather}$$
$$\begin{gather}{g}_{i}^{(eq)} = r\rho'{\omega}_{i} \left[\epsilon+\frac{\boldsymbol{e}_{i}\boldsymbol{\cdot}\boldsymbol{u}}{c_s^2}\right], \end{gather}$$
with the unit tensor 
$\boldsymbol {I}$ and weighting factors 
$\omega _0=4/9$, 
$\omega _{1,\dots,4}=1/9$ and 
$\omega _{5,\dots,8}=1/36$. Equations (B3) and (B4) also hold for the fluid region by setting 
$\epsilon =1$. The force terms are given by
$$\begin{gather} {F}_{i} =\delta t\left(1-\frac{1}{2\tau_f}\right) \frac{(\boldsymbol{e}_{i}-\boldsymbol{u})\boldsymbol{\cdot}\boldsymbol{\tilde\alpha}}{c_s^2} {f}_{i}^{(eq)}, \end{gather}$$
$$\begin{gather}{G}_{i}= \delta t\left(1-\frac{1}{2\tau_g}\right) {\omega}_{i} \boldsymbol{e}_{i}\boldsymbol{\cdot}\boldsymbol{\beta}, \end{gather}$$with
$$\begin{gather} \tilde \alpha_z = \alpha_z , \end{gather}$$
$$\begin{gather}\tilde \alpha_r = \alpha_r +\frac{c_s^2}{r}\left(1-\frac{(2\tau_f-1)u_r}{r}\right), \end{gather}$$
where 
$\boldsymbol {\alpha }=(\alpha _r,\alpha _z)=-(\varepsilon \nu /K)\boldsymbol {u}$ represents the Darcy force in the porous medium that vanishes in the fluid domain. The additional term 
$\boldsymbol{\beta}=(\beta_r,\beta_z)=(\rho', 0)$ is the same for both the fluid and porous domains. The non-dimensional relaxation times were defined by 
$\tau _f=\nu /(c_s^2\delta t)+1/2$ and 
$\tau _g=D/(c_s^2\delta t)+1/2$, where 
$c_s$ is the speed of sound (assuming 
$c_s^2=1/3$) and 
$\delta t$ is the time increment (Succi Reference Succi2001).
The hydrodynamic variables mass density (
$\rho$), momentum (
$\boldsymbol {j}$) and density perturbation (
$\rho '$) were computed by solving
$$\begin{gather} \rho=\frac{1}{r}\sum_{\textrm{i}} {f}_{i}, \end{gather}$$
$$\begin{gather}\boldsymbol{u}=\frac{\displaystyle\sum_{\textrm{i}}\boldsymbol{e}_{i} {f}_{i}}{\dfrac{r^2+\tau_f\delta tc_s^2}{r}+r\epsilon\dfrac{\delta t}{2} \dfrac{\nu}{K}}, \end{gather}$$
$$\begin{gather}\rho^\prime=\frac{1}{r\epsilon}\sum_{\textrm{i}}{g}_{i}{.} \end{gather}$$ A grid refinement technique was applied with a hierarchy of 
$m$ nested grids with successively finer resolution approaching the system centre (figure 17a). Each nested grid consisted of 
$x_n=128 n$ and 
$y_n=257 n$ grid nodes. The resolution of the finest grid was defined through the amount of grid nodes placed into the particle: 
$\delta _x = {1}/{(80 n)}$. For the coarser grids, 
$\delta _x$ and 
$\delta _t$ were then sequentially scaled by a factor of 2. Values for 
$n\in [1,2,3,4]$ and 
$m\in [4,6]$ were selected based on the convergence for the parameter choice of the model runs. Lower limit of domain size was 
$n=1$ with 
$m=6$ resulting in a grid width of 
$x_n \times \delta _x \times 2^{n-1} = 51.2 a$ and a resolution of 
$0.013 a$. To ensure the continuity of pressure, velocity, their derivatives and density at the interface between coarser and finer grids, the interpolation method of Liu & Khalili (Reference Liu and Khalili2008, Reference Liu and Khalili2009) was used. The non-local regularised scheme of Liu & Khalili (Reference Liu and Khalili2008, Reference Liu and Khalili2009) was applied to satisfy the far-field boundary condition (constant velocity 
$W$ far from the aggregate). In this method, the macroscopic strain tensor 
$\boldsymbol {S}$ obtained using a non-local finite-difference scheme was used to evaluate the regularised distribution functions. For validation purposes, we also computed the flow around solid spheres satisfying the no-slip condition at the surface, using a non-equilibrium extrapolation scheme (Guo & Zheng Reference Guo and Zheng2002). The criteria for stationarity are based on both the drag coefficient and the interstitial density contrast. The model was considered to be converged when changes were less than 
$10^{-7}$ between 1000 time steps. Examples are given in figure 17(c–f) for the lower and upper bounds of parameters chosen.
Figure 16. Time trajectories of the aggregate settling experiments in stratification for the two cases corresponding to figure 15: (a) 
$\textit {Re}=0.66$, 
$\textit {Fr}=1.37$ and (b) 
$\textit {Re}=0.1$, 
$\textit {Fr}=0.88$. The black crosses indicate the location of the aggregate during experiments, while orange represents the direct comparison with the model (4.1) (see also figure 10). For comparison, the modelled trajectory of a solid sphere of equal initial excess density settling in the absence of stratification (black curve) and in stratification (dashed line) is depicted.
Figure 17. (a) Schematic of the grid set-up. (b) Drag coefficient based on model runs of a homogeneous fluid for various 
$\textit {Re}$ in comparison with White (Reference White2005) and Stokes’ law. Examples of the model convergence for (c,e) a stratified fluid based on the drag coefficient 
$C^S_D$ and (d,f) the density contrast with respect to an unstratified fluid 
${\Delta \rho _p}/{-a\gamma }$ (see § 2.2 for further information).
Appendix C. Application of derived scaling laws
The scaling laws derived for 
$\Delta \rho _p$ ((3.12), (3.13)), 
$\Delta \rho _\delta$ ((3.3), (3.4)) and 
$C^N_D$ ((3.14), (3.15)) allow for a straightforward implementation of stratification effects into force balance models and direct calculation of the increase in retention time. The potential applications are manifold which is the reason for describing the implementation in more detail here. The balance of buoyancy and drag force is given by
\begin{equation} \frac{4}{3}{\rm \pi} a^3 \rho_a\frac{\textrm{d} W}{\textrm{d} t}=(\rho_a -\rho_B) \frac{4}{3}{\rm \pi} a^3g + \frac{\rho}{2}C_D^{S}{\rm \pi} a^2W^2, \end{equation}
where the Basset history force was neglected. The particle density is defined as 
$\rho _a = (1-\epsilon )\rho _s + \epsilon (\rho _B+\Delta \rho _p)$, with 
$\rho _s$ the hydrated excess density of the solids, while the total excess density is given by
\begin{align} \rho_a -\rho_B &= (1-\epsilon)\rho_s + \epsilon (\rho_B+\Delta\rho_p)-\rho_B \end{align}
\begin{align} &= (1-\epsilon)(\rho_s-\rho_B) + \epsilon\Delta\rho_p. \end{align}Therefore, in regime R2 the excess density of the particle can be calculated as
\begin{equation} \rho_a -\rho_B = (1-\epsilon)\Delta\rho_s + 39a\gamma\epsilon(\textit{Re} +\textit{Re}^{0.43}\textit{Fr}^{0.26}) \\ \end{equation}and in R3 as
\begin{equation} \rho_a -\rho_B = (1-\epsilon)\Delta\rho_s + 40a\gamma\epsilon(\textit{Re} +\textit{Re}^{0.45}\textit{Fr}^{0.06}), \end{equation}
where 
$\Delta \rho _s =\rho _s-\rho _B$. For biogeochemical applications, it might be beneficial to calculate the diffusive exchange using the non-stationary equation (C6) which can be represented through a simple shell model:
\begin{equation} \frac{\partial\rho_p}{\partial t} \sim 3\frac{D}{\epsilon a^2} \left(\rho_p-\rho_B-\alpha\textit{Re}^{n}\textit{Fr}^{m} \right), \end{equation}
with 
$\alpha = 39 a \gamma$, 
$n=0.43$, 
$m=0.26$ in regime R2 and 
$\alpha = 40 a \gamma$, 
$n=0.45$, 
$m\approx 0$ in R3. We solved the ordinary differential equation in a proof-of-principle application to estimate the effect of the density boundary layer in buffering the diffusive exchange. It is important to notice that 
$\rho _p$ is assumed to be constant within the entire particle volume (instantaneous homogenisation). Therefore, the full dynamic and complexity at the surface is not captured; however, this is a powerful approach for obtaining an order-of-magnitude approximation. If higher accuracies are desired, the analytical solutions given by (3.11) and (3.10) may be used.
Appendix D. Settling experiments in linear stratification
To determine settling velocities of aggregates in stratification, we performed experiments with porous fibre aggregates resembling marine snow (see § 2). In total, 11 different density stratifications were generated by adjusting glycerol–water mixtures (figure 13a) yielding buoyancy frequencies ranging from 
$N=0.15$ to 
$0.50\,\textrm {s}^{-1}$. Numbers of 3 to 5 aggregates with a radius of 
$r=6.8\pm 0.2\,\textrm {mm}$ were released into the density-stratified fluid and the settling velocities were determined (see § 2 and Appendix A). The particle settling velocity was found to be highly variable spanning almost two orders of magnitude ranging from 
$W=0.5\times 10^{-3}$ to 
$8.4\times 10^{-3}\,\textrm {m}\,\textrm {s}^{-1}$. The theoretical settling velocity for a solid sphere of equal initial excess density settling in the absence of stratification ranges between 
$W=1.3\times 10^{-3}$ and 
$8.6\times 10^{-3}\,\textrm {m}\,\textrm {s}^{-1}$. The wide range of settling velocities and buoyancy frequencies resulted in a Reynolds number range 
$\textit {Re}=0.04$ to 1.36 and a Froude number range 
$\textit {Fr}=0.16$ to 5.8 (figure 13b and table 1). The measured Reynolds numbers for porous particles in stratification are decreased by on average 37 % and for strong stratification by up to 67 % when compared with theoretical values for a solid sphere with equal initial excess density settling in a homogeneous fluid 
$\textit {Re}_H$.
The time trajectories of the experiments are directly compared with model results ((4.1); cf figure 16). Overall, the time trajectories of the model and experiments match well even though some deviations are visible especially at the lower part of the settling chamber where density gradients became weaker. In both modelling and experimental results, the retention of the aggregate is strongly enhanced compared to a solid sphere of equal initial excess density settling in the absence of stratification. For the examples depicted, the increase in retention time 
$t/t_0-1$ within experiments was 29 % and 66 %, while the model predicted slightly lower values of 20 % and 26 % (figure 16). The increase in retention time for experiments with porous particles in stratification ranged between 10 % and 200 % (table 1).
Appendix E. Extended discussion and comparison with two-layer stratification
Settling of porous particles in two-layer stratification is a topic that has attracted much attention in recent years. Studies have provided new insights into the retention at sharp density interfaces and their potential biological implications (e.g. Prairie et al. Reference Prairie, Ziervogel, Arnosti, Camassa, Falcon, Khatri, McLaughlin, White and Yu2013; Camassa et al. Reference Camassa, Khatri, McLaughlin, Prairie, White and Yu2013). Independent of ‘diffusion-limited’ or ‘entrainment’ regime, the falling of a porous particle through two-layer stratifications is a transient problem, in which the particle decelerates when approaching the sudden increase of the ambient density followed by an acceleration due to the equilibration of the interstitial pore water of the porous particle with the ambient density. In this special case, it is of particular interest to study the change of the density inside the particle, the length scale of the density transition between the layers and the time the particle is retarded at the density interface (Prairie et al. Reference Prairie, Ziervogel, Arnosti, Camassa, Falcon, Khatri, McLaughlin, White and Yu2013; Panah et al. Reference Panah, Blanchette and Khatri2017). In these studies, the buoyancy frequencies were usually above unity and Reynolds numbers were in an intermediate range 
$1< Re<100$ (see figure 12), resembling heavy particles and stratifications that are only found in some aquatic environments, such as estuaries and salt lakes, which potentially limit the application range. By contrast, in typical oceanic conditions the changes in stratification exceed diffusive relaxation times and associated length scales of marine aggregates. Under these conditions and while settling, the changes of density in the external field are equivalent to the changes of density in the density boundary layer and interstitial pore water leading to a stationarity of the settling velocity (
$W=\textrm {const.}$). Centred around this stationarity, we rationalised the diffusive exchange processes in scaling laws which allow for an implementation in larger modelling frameworks such as Earth system models or regional biogeochemical models.
We tested whether the results of Panah et al. (Reference Panah, Blanchette and Khatri2017) can also be extrapolated to the case of weaker linear stratifications and lighter particles when assuming that the two-layer stratification changes over a distance of 100 particle radii (potentially similar to linear stratification). For that, we parameterised the individual retention time scalings from Panah et al. (Reference Panah, Blanchette and Khatri2017) based on the thickness of the two-layer stratification transition zone (
$\gamma =100$; Panah et al. (Reference Panah, Blanchette and Khatri2017), their (20)), the Péclet number (Panah et al. (Reference Panah, Blanchette and Khatri2017), their (21)), the stratification parameter (Panah et al. (Reference Panah, Blanchette and Khatri2017), their (23)) and the Reynolds number (Panah et al. (Reference Panah, Blanchette and Khatri2017), their (24)) and subsequently summed the individual effects to get the overall increase in retention time. When comparing these results with the outcomes of our scaling laws (figure 18), we found substantial deviations spanning up to four orders of magnitude. The deviations reduce for stronger stratifications and larger particles, i.e. the parameter range that was studied by Panah et al. (Reference Panah, Blanchette and Khatri2017). Thus, we conclude that the applicability of their model is somewhat limited when it comes to linear stratification, i.e. it is restricted to stronger stratification 
$N>1\,\textrm {s}^{-1}$.
Figure 18. Deviation between the results of Panah et al. (Reference Panah, Blanchette and Khatri2017) and those of our study. (a) The retention increase 
$\log (t/t_0-1)$ of a relatively heavy particle (
$\rho _s=100\,\textrm {kg}\,\textrm {m}^{-3}$). (b) The retention increase of a light particle (
$\rho _s=10\,\textrm {kg}\,\textrm {m}^{-3}$). The contour lines are 
$\log (t/t_0-1)$ based on Panah et al. (Reference Panah, Blanchette and Khatri2017) (their (20), (21), (23), (24)) and the colours indicate the relative deviation from results of our study. Red colour implies that the model of Panah et al. (Reference Panah, Blanchette and Khatri2017) predicts a higher retention time compared with ours, blue indicates that our retention time estimates are higher and white implies matching results.
Appendix F. Permeability effects
We tested variations of the Darcy number (at 
$\textit {Re}=1$ and 
$\textit {Fr}=1$) to investigate if there is a potential effect on the derived scaling laws (cf. figure 14). We observed that a Darcy number of 
$\textit {Da} = 10^{-5}$ is required for the permeability of the particle to substantially affect the particle density contrast as well as the drag coefficient. At this 
$\textit {Da}$ the effects of the stratification on drag are diminishing. In a dedicated study, Dörgens et al. (Reference Dörgens, Ahmerkamp, Müssig, Stocker, Kuypers, Khalili and Kindler2015) found the permeabilities of the laboratory fibre aggregates to range between 
$k=10^{-11}$ and 
$10^{-12}$ m resulting in 
$\textit {Da}\approx 10^{-6}$ (see figure 14, grey bar). This implies that the experiments with the laboratory particles are only slightly influenced by permeability effects (
${<}5\,\%$ for the drag coefficient).
Open access
