Settling of highly porous and impermeable particles in linear stratification: implications for marine aggregates

Abstract The settling velocity of porous particles in linear stratification is affected by the diffusive exchange between interstitial and ambient water. The extent to which buoyancy and interstitial mass adaptation alters the settling velocity depends on the ratio of the diffusive and viscous time scales. We conducted schlieren experiments and lattice Boltzmann simulations for highly porous (95 %) but impermeable spheres settling in linear stratification. For a parameter range that resembles marine porous particles, ‘marine aggregates’, i.e. low Reynolds numbers ($0.05\leq \textit {Re}\leq 10$), intermediate Froude numbers ($0.1\leq \textit {Fr}\leq 100$) and Schmidt number of salt ($\textit {Sc}=700$), we observe delayed mass adaptation of the interstitial fluid due to lower-density fluid being dragged by a particle that forms a density boundary layer around the particle. The boundary layer buffers the diffusive exchange of stratifying agent with the ambient fluid, leading to an enhanced density contrast of the interstitial pore fluid. Stratification-related drag enhancement by means of additional buoyancy of dragging lighter fluid and buoyancy-induced vorticity resembles earlier findings for solid spheres. However, the exchange between density boundary layer and pore fluid substantially increases stratification drag for small $\textit {Fr}$. To estimate the effect of stratification on marine aggregates settling in the ocean, we derived scaling laws and show that small particles ($\leq$0.5 mm) experience enhanced drag which increases retention times by 10 % while larger porous particle (>0.5 mm) settling is dominated by delayed mass adaptation that diminishes settling velocity by 10 % up to almost 100 %. The derived relationships facilitate the integration of stratification-dependent settling velocities into biogeochemical models.

The settling velocity of porous particles in linear stratification is affected by the diffusive exchange between interstitial and ambient water. The extent to which buoyancy and interstitial mass adaptation alters the settling velocity depends on the ratio of the diffusive and viscous time scales. We conducted schlieren experiments and lattice Boltzmann simulations for highly porous (95 %) but impermeable spheres settling in linear stratification. For a parameter range that resembles marine porous particles, 'marine aggregates', i.e. low Reynolds numbers (0.05 ≤ Re ≤ 10), intermediate Froude numbers (0.1 ≤ Fr ≤ 100) and Schmidt number of salt (Sc = 700), we observe delayed mass adaptation of the interstitial fluid due to lower-density fluid being dragged by a particle that forms a density boundary layer around the particle. The boundary layer buffers the diffusive exchange of stratifying agent with the ambient fluid, leading to an enhanced density contrast of the interstitial pore fluid. Stratification-related drag enhancement by means of additional buoyancy of dragging lighter fluid and buoyancy-induced vorticity resembles earlier findings for solid spheres. However, the exchange between density boundary layer and pore fluid substantially increases stratification drag for small Fr.
To estimate the effect of stratification on marine aggregates settling in the ocean, we derived scaling laws and show that small particles (≤0.5 mm) experience enhanced drag which increases retention times by 10 % while larger porous particle (>0.5 mm)

Introduction
Settling of porous particles in stratified fluids is a ubiquitous process in many environments, for example in the ocean (Asper et al. 1992;Li, Yuan & Wang 2003). 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 2016;Maerz et al. 2020). Since ocean stratification is projected to intensify under ongoing climate change (Bopp et al. 2013), 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 O(0.01) ≤ Re ≤ O(10) showed increased retention times in stratification (Asper 1985;Alldredge & Gotschalk 1988;Riebesell 1992;Alldredge & Crocker 1995;MacIntyre, Alldredge & Gotschalk 1995;Alldredge et al. 2002). We here define the Reynolds number as Re = aW/ν for sphere radius a, stationary sinking velocity W and kinematic molecular viscosity ν. Stratification, expressed in terms of the Brunt-Väisälä frequency, N, typically ranges from N = O(0.001 s −1 )-O(0.1 s −1 ) under oceanic conditions up to N O(1 s −1 ) in estuaries and fjords (Boehrer & Schultze 2008;Geyer, Scully & Ralston 2008). The Brunt-Väisälä frequency is defined as N = √ −(g/ρ 0 )γ , where g is the gravitational acceleration, ρ 0 a reference density and γ = ∂ρ/∂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 1992;Alldredge & Crocker 1995), or the buoyancy of the interstitial fluid carried by the porous particles to cause increased retention during settling in stratification (Kindler, Khalili & Stocker 2010). Such retention is particularly associated with pycnocline layers, which are defined through strong vertical density gradients (MacIntyre et al. 1995;Alldredge et al. 2002). 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. 2010;Camassa et al. 2013;Prairie et al. 2013Prairie et al. , 2015Panah, Blanchette & Khatri 2017). 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 ξ = δρ/ ρ 0 (comparing the density increment over the pycnocline layer, δρ, with a particle's excess density in the lighter phase, ρ 0 ) and the particle permeability to flow k (Kindler et al. 2010;Camassa et al. 2013;Prairie et al. 2013Prairie et al. , 2015Panah et al. 2017;Prairie et al. 2019). In the case of impermeability to flow (which is a fair approximation for most marine 931 A9-2 Porous particles settling in stratification δ ρ p ρ δ ρ z r W a Figure 1. Conceptual representation of stationary, stratified flow around a porous and impermeable sphere of radius a and density ρ 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 δ and density ρ δ . The red dot indicates the stagnation point. aggregates; Moradi et al. 2018) 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 (Re ≤ 1) and large Schmidt number (Sc = ν/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 1999;Higginson, Dalziel & Linden 2003;Yick et al. 2009;Zhang, Mercier & Magnaudet 2019). Stratification effects also cause extra stresses at the sphere's surface by the formation of toroidal vortices due to baroclinic torque (List 1971;Ardekani & Stocker 2010;Zhang et al. 2019) which has been shown to be the dominant contribution to drag enhancement at low to intermediate Re (Zhang et al. 2019). 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 Re, Sc and the Froude number, Fr = W/(aN) (Zhang et al. 2019).
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, ξ < 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. 1999;Camassa et al. 2010). At moderate Reynolds numbers (5 < Re < 15), the wake pinches off rapidly and retreats upwards, while in the Stokes regime (Re 1), the stem of entrained fluid dissolves due to diffusion and, hence, retardation increases (Camassa et al. 2010). If, however, a particle is impermeable to flow and the density increment exceeds the excess density of the particle, hence ξ > 1, settling itself becomes diffusion-limited (Li et al. 2003;Kindler et al. 2010;Panah et al. 2017), 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. 2020), which is much larger than typical sizes of marine aggregates (O(1 μm) ≤ a ≤ O(1 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 = 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. 2019). 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 ( = 0.95, ratio of void to total volume) and impermeable spheres (here defined as k 10 −12 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.

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 cm with height h c = 38 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 1965). The fluids in the buckets were aqueous glycerol mixtures with contents that ranged from 60 to 90 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 kg 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 T = ±0.5 • C. Highly porous, but impermeable particles were produced using fibres as described in Dörgens et al. (2015). Briefly, in a semi-manual process, polyester fibres with 20 μm diameter were agglomerated to nearly perfect spheres. The porosities of 10 replicates were in the range = 92 % to 94 % with mean radius a = 0.68 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 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 cm × 5 cm, at a resolution of 27 μm 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 (Ĩ ref ) was subtracted from each image including the aggregate (I) and used to normalise:Ĩ = ( In the reference imageĨ, 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 2007). 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 mm, f number = 1.4 to 16) installed L c = 20 cm in front of the tank. The field of view was 10 cm × 5 cm, at a resolution of 49 μm 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 px × 24 px with 50 % overlap, resulting in a 12 px × 12 px grid. In this set-up, the maximal traceable shifts were ∼10 μm at a spatial resolution of 588 μ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 O(1000-10 000). Permeabilities of the fibre particles were measured to be k = 3.87 × 10 −12 m resulting in Darcy numbers of O(10 −6 ) (Dörgens et al. 2015); see also figure 13 and table 1 for an overview of the experiment parameters.

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. 2000;Larrazabal, Torres & Castillo 2003;Yick et al. 2009;Liu, Kindler & Khalili 2012). We considered the flow of a linearly stratified fluid at constant velocity W past a stationary sphere: where length is scaled by a, velocity by the undisturbed velocity W (resembling the constant settling velocity), pressure perturbations by ρ 0 W 2 and density perturbations ρ by −aγ . The density perturbation represents the density contrast induced by the particle in comparison to an undisturbed linearly stratified fluid: where ρ B is the density of the undisturbed ambient fluid in linear stratification and −aγ is the vertical density change along one particle radius. Therefore the density contrast ρ/(−aγ ) can be interpreted as the deflection distance of the pycnoclines from their equilibrium position.
Here, u = (u, w) is the fluid velocity vector with radial and vertical components, p is the pressure, j is the vertical unit vector, positive upwards, is the porosity of the particle and Da = k/a 2 is the Darcy number, with k the permeability. The Darcy number was adjusted to Da = 10 −12 , which is comparable to permeabilities of k = 10 −18 to 10 −16 m (for particle radii in the range a = 10 −4 to 10 −3 m) and, therefore, resembling that for typical impermeable marine particles (Kiørboe et al. 2002) (see also figure 14 for tests on permeability effects). The fluid domain is defined by = 1, A = 0, where A is a binary switch, and the particle by = 0.95, A = 1 (Basu & Khalili 1999). 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 ≤ Re ≤ 10 and 0.1 ≤ Fr ≤ 100 at Schmidt number Sc = ν/D = 700, where ν 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 S D was computed as the sum of the pressure and viscous drag coefficients, C S P and C S V , respectively: where n is the unit vector normal to the sphere's surface S and μ the dynamic viscosity. Henceforth, we express the stratification drag force normalised by the drag force in the absence of stratification (C H D ): In order to explore the mechanisms governing the drag enhancement by stratification, a force decomposition scheme was applied following Zhang et al. (2019). The velocity and pressure fields were decomposed into the form u = u H + u ρ and p = p H + p ρ with where (u H , p H ) are the velocity and pressure in the homogeneous fluid and (u ρ , p ρ ) are density-induced contribution in stratified fluid. The buoyancy-induced pressure disturbance in the ambient fluid can further be decomposed into the form p ρ = p ρρ + p ρu + p ρω , where p ρω obeys a Laplace equation, while the remaining two contributions satisfy Fr 2 ∇ρ · j, (2.7) with p ρρ and p ρu vanishing in the far field. These two Poisson equations were solved by a lattice Boltzmann method (Chai & Shi 2008). The external stratification-induced forces can then be expressed as is the vorticity-induced stress tensor. The contribution F ρρ is an additional Archimedes-like force due to the deflection of the isopycnals, F ρu is an inertial force associated with u p and F ρω results from the vorticity-induced baroclinic torque (vortical force). The contributions of F ρu and F ρρ are directly calculated by solving the Poisson equations (2.7) and (2.8); F ρω is subsequently determined as 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 2005) 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. 2019) (figure 3b) showing excellent agreement for the normalised stratification drag C N D , as well as the two force components F ρρ and F ρω at Re = 0.05. Results of additional tests on the convergence criteria for various Re and Fr are shown in figure 17.
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. 2009;Zhang et al. 2019). 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 Re = 0.66 and Fr = 1.37, the boundary layer is narrow, resulting from the dominance of buoyant forces, while for Re = 0.1 and 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 Sc = 700, we subsequently focus on the numerical model and refer to experimental validations whenever possible.

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 1971;Ardekani & Stocker 2010) and the recent analysis of solid-particle settling in density gradients (Zhang et al. 2019).
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. (2019), the flow structure shows a strong dependency on stratification strength and viscosity for large Péclet numbers Pe = aW/D = ScRe ≥ 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  determined by viscous forces, corresponding to regime 3, R3, Fr Re −1 (Zhang et al. 2019).
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 Fr and large Re, while as the stagnation point recedes from the surface (Fr 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 Fr, and large Re.
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 Fr and 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.

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, δ, buffers the diffusive exchange between the external density field and the interstitial liquid.
Here, δ 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 δ by 60 % and an 80 % threshold will decrease δ 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 Fr ≈ Re −1 (figure 6a). The best fit to our results ( figure 6b,c)  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 δ ∼ (ν/N) 1/2 , which implies δ/a ∼ (Fr/Re) 1/2 (Yick et al. 2009). In regime R3 the influence of Fr, i.e. buoyant force, becomes negligible Re -0.66 Fr 0.61 Re -0.71 Fr 0.14 (e) Figure 6. Scaling of the boundary layer thickness δ/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 ρ δ 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 ρ p /(−aγ ) 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 Re and 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). and the boundary layer thickness is largely defined by the interaction of inertial and viscous forces, represented through 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 ρ δ /(−aγ ) was evaluated as the average within the shell with width δ/a surrounding the particle (figure 6d-f ). The best collapse (in terms of least square errors) was found for where ρ δ represents the density of the density boundary layer and ρ δ 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).

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 ρ p was evaluated as the spatial average of the interstitial density contrast: where ρ = ρ p − ρ B . The interstitial pore water density was found to increase as ρ p ∼ 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 ρ δ is reduced with respect to an undisturbed density field ρ. 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: where ρ δ can be exchanged by our scaling relationships (3.3) and (3.4) for the two regimes (figure 6d-f ): with α = 39aγ , n = 0.43, m = 0.26 in regime R2 and α = 40aγ , n = 0.45, m ≈ 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: Under these conditions a direct relationship between ρ p /−aγ and ρ δ /−aγ can be derived. The transport equation (2.2) inside the porous particle reduces to the non-dimensional diffusion equation: The density flux through the surface of the particle is balanced by the changing external field: 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: where c is the integration constant. Based on (3.11) the averaged density on the surface is ρ δ /(−aγ ) = 1 6 ReSc + c and the averaged density inside the particle is ρ p /−aγ = 1 10 ReSc + c. The difference of the averaged density of the interstitial pore water and the density in the boundary layer then is ρ p /−aγ = −44Re + ρ δ /−aγ (for Sc = 700), which implies that the scaling relationship ρ δ also applies for ρ p if shifted by 44Re. Indeed, best fit to our simulations gives close results to the analytical solution (figure 6g-i): for the non-dimensional interstitial fluid density contrast with respect to the ambient fluid. In regime 3, the effects of 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.

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 Note that we attribute all external forces, the additional Archimedes-like forces of the density boundary layer (F ρρ ), the inertial force (F ρu ) as well as forces due to vorticity (F ρω ) 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 In R2, the drag coefficient is close to that for the solid case C D ∼ Ri 0.5 confirming the results of Yick et al. (2009). In R3 the increase in drag is largely independent of Re and is mainly controlled by Fr, i.e. the ratio of inertia and stratification strength. 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 ρρ and F ρω as compared to the case of solid particles ( figure 8a,b).
To quantify within the two regimes, we found the vortical force to scale as which is very similar to the values of a solid particle which Zhang et al. (2019) found to scale as F ρω,R2 ∼ (ReFr) −0.67 and F ρω,R3 ∼ (ReFr) −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. (2019) found F ρω,R3 ∼ Re −0.5 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 ρω 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 ρρ was found to play a secondary role for stratification drag enhancement (Zhang et al. 2019  Re = 0.05, solid . (3.20) 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 Re and weak stratification (Fr > 1) the contribution of F m is negligible ( figure 9a). However, at larger Re and stronger stratification (Fr < 1) the density contrast is substantially increased and F m even exceeds the contribution of the external forces C N D (figure 9b). In direct comparison between solid and porous particles, we found for 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 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 ρρ and F m , are non-negligible and comprise a range similar to that of F ρω .

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 μm), known as 'marine snow ' (MacIntyre et al. 1995;Alldredge et al. 2002;Prairie et al. 2015).
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 s −1 ). We consider the balance of buoyancy and drag force: where the Basset history force was neglected. The particle density is defined as ρ a = (1 − )ρ s + (ρ B + ρ p ), with ρ 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 ( 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 .
represents the flow regime of particles within stratification. This implies, for example, that a porous particle featuring a particular Re in linear stratification (and hence W being constant) will have a higher Re under homogeneous fluid conditions (assuming ν to be identical). Our scaling laws can thus be applied to a wider 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 × 10 −4 s −1 . For weak stratification (N ≤ 0.02 s −1 ) or small particles (a < 200 μ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 mm and reduces their sinking velocities by 10 % to 100 %. Under strong stratification (N > 0.2 s −1 ), the sinking velocity of large particles (a > 1 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 s −1 ) and abundant particle sizes (a = 1 mm) the settling velocity is expected to be reduced by 5 %-15 % due to the combined effect of drag enhancement and delayed mass adaptation.
In contrast to diffusion-limited retention in two-layer stratification, i.e. sharp pycnoclines (Panah et al. 2017), 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 = 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 . 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 ρ s , thus Re and 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 1998) 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 1998). 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. 2013) 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. 2020). 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.

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 931 A9-19 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 mm) is mainly affected by enhanced drag that increases the retention by up to 10 %, while for larger particle sizes (a > 0.5 mm) delayed mass adaptation reduces settling velocity by 10 % and up to ∼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. refractive index (n 0 ), the deflection angle is related to the Radon transform of the density perturbation (ρ ) by a constant a: the deflection angles and the density field reads where k x denotes the wavenumber in the Fourier domain. From (A1) it follows that F (ρ ) = F (ψ)/i2πk x and calculating the inverse Fourier transform yields where Q(k x ) is a ramp filter applied in the Fourier domain. Applying the convolution theorem and introducing the undisturbed background density (ρ B ) yields 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. ρ read where f i (x, t) and g i (x, t) are the single-particle distribution functions at position 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 1992), the discrete velocities were defined as e i = (e ir , e iz ) : i = 0, 1, . . . , 8 with e 0 = 0, e 1 = −e 3 = (1, 0), e 2 = −e 4 = (0, 1), e 5 = −e 7 = (1, 1) and e 6 = −e 8 = (−1, 1). The equilibrium distribution functions read with the unit tensor I and weighting factors ω 0 = 4/9, ω 1,...,4 = 1/9 and ω 5,...,8 = 1/36. Equations (B3) and (B4) also hold for the fluid region by setting = 1. The force terms Δρ p /(-aγ)  Therefore, in regime R2 the excess density of the particle can be calculated as where ρ s = ρ s − ρ 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: with α = 39aγ , n = 0.43, m = 0.26 in regime R2 and α = 40aγ , n = 0.45, m ≈ 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 ρ 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.
ranging from W = 0.5 × 10 −3 to 8.4 × 10 −3 m 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 × 10 −3 and 8.6 × 10 −3 m s −1 . The wide range of settling velocities and buoyancy frequencies resulted in a Reynolds number range Re = 0.04 to 1.36 and a Froude number range 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 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. stratification parameter (Panah et al. (2017), their (23)) and the Reynolds number (Panah et al. (2017), 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. (2017). 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 s −1 .