Resolution dependency of sinking Lagrangian particles in ocean general circulation models

Any type of non-buoyant material in the ocean is transported horizontally by currents during its sinking journey. This lateral transport can be far from negligible for small sinking velocities. To estimate its magnitude and direction, the material is often modelled as a set of Lagrangian particles advected by current velocities that are obtained from Ocean General Circulation Models (OGCMs). State-of-the-art OGCMs are strongly eddying, similar to the real ocean, providing results with a spatial resolution on the order of 10 km on a daily frequency. While the importance of eddies in OGCMs is well-appreciated in the physical oceanographic community, other marine research communities may not. Further, many long term climate modelling simulations (e.g. in paleoclimate) rely on lower spatial resolution models that do not capture mesoscale features. To demonstrate how much the absence of mesoscale features in low-resolution models influences the Lagrangian particle transport, we simulate the transport of sinking Lagrangian particles using low- and high-resolution global OGCMs, and assess the lateral transport differences resulting from the difference in spatial and temporal model resolution. We find major differences between the transport in the non-eddying OGCM and in the eddying OGCM. Addition of stochastic noise to the particle trajectories in the non-eddying OGCM parameterises the effect of eddies well in some cases (e.g. in the North Pacific gyre). The effect of a coarser temporal resolution (once every 5 days versus monthly) is smaller compared to a coarser spatial resolution (0.1° versus 1° horizontally). We recommend to use sinking Lagrangian particles, representing e.g. marine snow, microplankton or sinking plastic, only with velocity fields from eddying Eulerian OGCMs, requiring high-resolution models in e.g. paleoceanographic studies. To increase the accessibility of our particle trace simulations, we launch planktondrift.science.uu.nl, an online tool to reconstruct the surface origin of sedimentary particles in a specific location.

We used backward in time calculation to be consistent with the studies in [26,14]. These studies focus on particles which are found at the bottom of the geological ocean.
Author's response Good point.

Reviewer 2
We would like to thank the reviewer for his/her careful reading and his/her constructive comments.
Please find our point by point reply and the changes in the manuscript below.
On behalf of the authors,

Peter Nooteboom
The manuscript describes how statistics of Lagrangian sinking particles differ when advected off-line in an eddy permitting or coarse resolution ocean model. While realizing the need for the work presented for some communities (paleo?), from a physical oceanography point of view, the content is not novel and mostly qualitative. It is well known that ocean models that are not eddy resolving underestimate horizontal dispersion. The literature on the subject is very large and some references should indeed be added. I note also that the manuscript in its current form reads like a physical oceanography work, and the appeal to the broader and more biologically oriented community of PLOS One is unclear to me. I do not have major comments on what is in the manuscript (there is nothing wrong per se), but I would have hoped for a more in-depth justification of the work and its limitations considering the lack of novelty, or a more in-depth investigation of the reliability of the high and low resolution results. The outcome of this work confirms deductions that have been made already in the absence of sinking, and many Lagrangian studies have been performed in this regard comparing runs at the exact resolutions used here (to mention few relevant to the content of this manuscript: Putman and He, 2013; Doos, Rupolo and Brodeau, 2011, and many intercomparison papers focused on the Eulerian fields).

Author's response
We would like to point out that the manuscript is intended for a broad audience, since the results are relevant for work that concerns the application of sinking particles in Eulerian models to biogenic matter that is important for the carbon flux to the ocean bottom [12,13,29], microplankton which is used for climate reconstructions [14,26] and sinking plastic [11]. We have adapted the manuscript accordingly and hope that the intended audience becomes clear.

Author's response
Indeed, a comparison of the back-tracking analysis with Eulerian models of higher horizontal resolution and different vertical resolution is interesting. Especially, it can be interesting for the modelling of sinking particles in the present-day ocean. It is not entirely clear how much of a difference it would make globally. The two model resolutions (1 • and 0.1 • horizontally) in this manuscript are respectively the highest resolutions that are typically used in palaeo and present-day models respectively. Due to long spin-up times, palaeo simulations at 10km resolution do not exist yet, and palaeo simulations at 2km resolution would be unfeasible. However, this issue is now addressed in the discussion, as it is interesting for future work.

Changes in manuscript
We added a sentence in the method section to be clear about the choice of the model resolutions: 'These configurations represent the differences between state-of-the-art, global OGCM resolutions of the past (1 • horizontally and monthly model output) and the present-day (0.1 • horizontally and model output on a daily scale). The single effect of model output with lower temporal resolution compared to the state-of-the-art present-day OGCMs is investigated in a separate configuration R 0.1m .' Two paragraphs in the discussion are added now, which also include the references to Stewart el al. 2017 [23] and Bracco et al. 2018 [4]: 'Future work could investigate the sinking Lagrangian particles in other configurations with different OGCM resolutions. The model output of configuration R 0.1 could be coarsened to a 1 • grid before the back-tracking analysis, to separate out the effects on the Lagrangian analysis of (a) a coarsened grid (see also [10,15]) and (b) a lower resolution of the underlying Eulerian model. Moreover, particle trajectories could be sensitive to the vertical resolution of the OGCM (e.g. [4,23]). Some applications require different horizontal model resolutions than the configurations used in this paper. An OGCM with a ∼1km horizontal resolution has a better representation of the spatial/temporal submesoscale that can be important for sinking Lagrangian particles which represent the carbon flux to the ocean bottom [3]. Although the mesoscale flow contains most of the energy that is responsible for the horizontal tracer dispersion [24], submesoscale (1-20km) dynamics have proven to be of importance for the vertical advection of iron in specific regions with strong flow-bathymetric interactions [18].

'
Secondly, what would happen is a reanalysis data-set was used instead? The most recent ones are run at 25km-¿10km near the poles and fields are saved monthly or more frequently. Do the statistics in the 10km run compare well with those from a reanalysis product, where at least surface eddies are constrained by satellite observations and the water column dynamics assimilate the large ARGO dataset? The investigation of the time-dependence of the results (advecting particles using fields saved every month or every 5 days is also superficial and should be expanded).

Author's response
The suggested comparison with sinking particle statistics in reanalysis data is interesting, but outside the scope of this paper. Here, we focus only on the differences in spatial resolution and not the type of Eulerian model. We chose to study the dependence of 5-daily and monthly output fields, because this is the typical timescale difference of model output between present-day and palaeomodels. For an analysis of 3D Lagrangian particles with more temporal resolutions we cite [16].
Thirdly, I understand the argument for justifying the simulation done for small plastic particles, but as far as I know organic particles small enough to have a sinking velocity of 6-20 m/day near the surface do not make it to the ocean bottom as such (vice versa particles that may reach the bottom with these sinking velocities were likely much larger and heavier at formation). Aggregation, grazing and bacterial activities are dominant processes on those sizes. The review by Boyd et al. in Nature in 2019 has a nice image summarizing the fate of sinking particles as function of their size, and shows that sinking is not the dominant process determining the fate of particles as small as those used in the manuscript (they do not reach more than 400-600 m depth in particle form). It is unclear, therefore, what is the paleo relevance that justifies the need for coarse resolution ocean models.
Author's response The review in [3] is mainly about the processes which are responsible for the carbon flux to the ocean bottom. As measured in e.g. [17], the fast sinking particles contribute most to the carbon flux, however slow sinking particles contribute significantly to the number flux. As explained in [3], the mechanisms that are responsible for this carbon flux vary in space and time. As a consequence, the slowly sinking particles can contribute significantly to the carbon flux [27,1].
It could be interesting to make the sinking speed dependent on dynamics such as aggregation, grazing and bacterial activities. As such, one could obtain a more detailed and realistic sinking of particles. However, we chose to keep the sinking constant to get a better idea of the resolution dependency of these results, instead of the sensitivity of the results to parameters which determine the sinking speed.
Finally, the sinking particles that are used for climate reconstructions are not necessarily subject to the mechanisms explained in [3], which increase the sinking speed of small particles, as is explained for dinoflagellate cysts in [14]. Therefore, particles which are sinking at velocities on the order of 6-25 m/day are relevant for these palaeoclimate reconstructions. For an analysis of the effect of sinking velocities on the particle trajectories, we refer to [14].

Changes in manuscript
Reference [3] is cited in the discussion, to explain the relevance of future work on sinking particles in OGCMs of 1km scale resolution (see 'changes in manuscript' under the first detailed comment; page 6).
Overall, in my opinion this is an exercise that deserves publication if better motivated, with a more in-depth discussion of implications and limitations, and the exploration of a greater range of sinking velocities, saving times, etc.
Author's response Modification in the paper based on the reviewer's comments led to a more in-depth discussion of implications and limitations. The comparison of R 0.1 with R 1md does not depend much on the sinking velocity (approximately the same value of c s for both sinking velocities minimizes the W d between these configurations, and matches their average travel distance and the surface area of the particle distributions optimally). We think that the exploration of more sinking velocities and saving times is out of scope for this manuscript and not necessary to answer the research question posed.

Reviewer 3
We would like to thank the reviewer for his/her careful reading and his/her constructive comments.
Please find our point by point reply and the changes in the manuscript below.
On behalf of the authors,

Peter Nooteboom
This manuscript considers the problem of identifying the source location of sinking particles found at the bottom of the ocean and examines the effect of coarse spatial (and to some extent temporal) resolution. The authors find that ocean general circulation models (OGCMs) at a resolution typical of state-of-the-art paleoclimate applications (1 degree) do not capture the distributions of potential source locations well. This can be remedied to some degree by using a Smagorinski parameterisation for the eddy effects that are not directly modeled. I find the experiment reasonably designed and the manuscript overall well written, although the organization is a bit odd. Some of the discussion could also benefit from a broader perspective, considering alternative explanations, and some of the conclusions are stated too strongly. Therefore, I recommend this submission for publication subject to minor revisions, which should address the points below. Detailed comments: 1. Experiment design: The case with highest spatial and temporal resolution was chosen as the reference case. This is reasonable, since one might anticipate its results to be closest to being accurate. However, it makes it difficult to disentangle effects of spatial and temporal resolutions. You might want to consider using R0.1m as the reference case, which can then be compared to (a) a case with the same spatial and different temporal resolution and (b) cases with different spatial and the same temporal resolution.

Author's response
The case with the highest spatial and temporal resolution was indeed chosen as a reference case because it is assumed to be the closest to reality. While we understand the point of the reviewer, we still like to make the comparison here with R 0.1 as reference case. Since anyone who wants to apply Lagrangian particles may often choose between 0.1 • OGCMs (often written with output of daily order and 1 • OGCMs (with monthly output), we regard the difference between these two configurations as most important.
The differences due to the temporal averaging are still shown by the differences between configurations R 0.1m and R 0.1 . The differences due to the resolution of the Eulerian model can still be inferred from configurations R 0.1m and R 1m .
2. Abstract, 4th sentence: State-of-the-art OGCMS run on regional grids nowadays often provide output at much greater resolution, e.g. 1 km in space and 3-hourly in time is fairly common. The authors may want to specify here (in the abstract) that they are considering global models only. It would also be good to discuss implications of the present findings for using a 10-km model with output every 5 days versus one with even higher spatial and temporal resolution. See also comment #22 below.
Author's response Good point. We changed the manuscript accordingly.

Changes in manuscript
The sentence in the abstract is changed: 'In this study, we simulate the transport of sinking Lagrangian particles using low-and highresolution :::::: global : OGCMs, and assess the lateral transport differences resulting from the difference in spatial and temporal model resolution.' We also added two paragraphs in the discussion to discuss the possibility of increasing the resolution more (see reviewer 2): 'Future work could investigate the sinking Lagrangian particles in other configurations with different OGCM resolutions. The model output of configuration R 0.1 could be coarsened to a 1 • grid before the back-tracking analysis, to separate out the effects on the Lagrangian analysis of (a) a coarsened grid (see also [10,15]) and (b) a lower resolution of the underlying Eulerian model. Moreover, particle trajectories could be sensitive to the vertical resolution of the OGCM (e.g. [4,23]). Some applications require different horizontal model resolutions than the configurations used in this paper. An OGCM of ∼1km resolution has a better representation of the spatial/temporal submesoscale that can be important for sinking Lagrangian particles which represent the carbon flux to the ocean bottom [3]. Although the mesoscale flow contains most of the energy that is responsible for the tracer dispersion [24], submesoscale (1-20km) dynamics have proven to be of importance for the vertical advection of iron in specific regions with strong flowbathymetric interactions [18].' Table 1: The phrase 5-daily is not entirely clear, as it could mean 5 times per day or once every 5 days. I recommend rephrasing it to clarify in these places. Further mentions should then be self-explanatory.

Abstract and
Author's response Good point.

Changes in manuscript
We changed '5-daily' to 'once every 5 days' in the abstract and the

Author's response
Thank you for these suggestions of similar studies that have been done.

Changes in manuscript
We added [29,21,5] to the other relevant literature on the topic in the introduction: 'The spatial resolution of the OGCM determines if the flow is eddying, which played an important role in simulations of sinking particles near the northern Gulf of Mexico [12] and in the Benguela region [13], and for passive tracers near Sellafield [21] and globally [5] (0.25 • versus 1 • resolution).' and 'The Lagrangian techniques are used to model the sinking particle trajectories in the modern ocean [30,19,20,13,28,29], specifically for sinking microplankton [26,14] and microplastic [7].' Paper [10] is cited in the discussion, as they did an experiment with coarsening of the model output before the calculation of Lagrangian trajectories, which relates well to the points raised in the discussion: 'The model output of configuration R 0.1 could be coarsened on 1 • grid before the back-tracking analysis, to separate out the effects on the Lagrangian analysis in (a) a coarsened grid [10] and (b) a lower resolution of the underlying Eulerian model.'

Author's response
It is indeed about 6 years, which is the time it takes until all released particles to reach the surface.

Changes in manuscript
The sentence is changed into: 'This means that we release particles at the bottom of the ocean every three days for more than a year, and compute their trajectories in the changing flow field back in time (similar to [14,25,29,12]) until the particles reached the surface.' 6. Lines 106 109: This description is a bit unclear. First, there seem to be 9 different model configurations, depending on whether the choice for cs is considered part of the configuration. I suggest to use 4 main configurations, whereby R1m should consistently be chosen with the same value of cs (Given the results, I would suggest cs = 2), with 5 additional experiments to investigate the sensitivity of the results to cs. It is not helpful that in the figures R1m does not consistently refer to the same model configuration. Secondly, line 108 back-tracking from a single release location appears to contradict line 61, which specifies a grid of release locations. Reading further, it becomes clearer what is meant, but I recommend editing the description here. Also, is the 130 particles used at each location a typo? For releases every 3 days for 6 years (cf. caption of Fig. 1), there should be 730 releases.

Author's response
We already use four main configurations, as is shown in table 1 and explained in Line 106. In the text it is convenient to refer to a configuration with and one without diffusion added to the dynamics (the 'd' in the subscript stands for 'diffusion'), so we choose to stay with these four configurations. In the revised manuscript, we will consistently use a value of c s = 2 in the plots.
The explanation on the back-tracking can indeed be more clear. We released a particle every three days, but not more than 130 particles per release location. At 6 m/day sinking velocity, the travel time can reach up to 4 years (see for example figure 7 of [14]), which implies that one has to wait until the sixth year before all particles reached the surface.

Changes in manuscript
No changes with respect to the decision of the configurations are made. See comment 5 for changes with respect to the explanation of the analysis.
7. Lines 110 117, Fig. 1, and Results organization: On these lines and in the figure, the three metrics are listed in an order of increasing complexity. The lateral distance is a metric of an individual trajectory; the area metric describes the variability among launches at a single location; and the Wasserstein distance also takes into consideration where in space the distribution sits. This makes sense. But the Results are organized in the opposite order, starting with the most complex metric. I highly recommend changing that and following the order presented on lines 110 117, making it easier for the reader to interpret the additional elements captured in each metric.

Author's response
We understand why the Wasserstein distance is considered by the reviewer as the most complex metric.

Changes in manuscript
We changed the order of the results accordingly, consistent with the order of the metrics in  The travel time also depends greatly on the depth of the release location. If the depth is 6km, the travel time could be respectively 2.8 and 0.7 years for sinking velocities 6m/day and 25m/day. However, the vertical ocean flow induce great variability in the travel time from one release location (see again figure 7 of [14]).

Changes in manuscript
We changed a sentence: 'The travel time of the particles is perhaps too short ::: (at ::::: most :: a ::::: few :::::: years) : for the errors in R 0.1m to grow substantially, and remains lower ::::::: smaller : compared to R 1m .' 9. Line 130: Well, they are not the same; one can clearly distinguish them in Fig. 1d. It might be better to state that they are not significantly different. Also, if you consider these two curves to be the same, then so are the curves for the larger sinking velocity, which is not how they are discussed in the next paragraph.

Author's response
This was a typo. The lines are indeed not the same.

Changes in manuscript
We changed the sentence: The global average W d between R 0.1m and the reference case (W d (R 0.1 , R 0.1m )) is the same as ::::::: slightly :::::: larger :::::::::: compared ::: to the check of R 0.1 with itself (the global average W d (R 0.1 , R 0.1 ); dashed versus dotted in Fig 2d).' 10. Lines 131 134: I recommend explaining this earlier, maybe at the beginning of the Wasserstein distance comparison section, before the reader looks at the plot and starts wondering why that quantity isnt identically 0.

Author's response
We think that the explanation is clear as it is now. If a reader wonders why this value is unequal to 0 in the figure, the explanation is also included in the figure caption.
11. Line 136: There seems to be a strong difference between cs = 0 and cs positive. It is not clear that cs = 2 is significantly different from the other positive cs values. With that in mind, the statement that there is a minimizing value is misleading.

Author's response
We changed the sentence and hope that it is not experienced as misleading anymore.

Author's response
We added a figure in the supporting material with the mean surface eddy kinetic energy.

Changes in manuscript
The caption of the added supporting figure: 'Geographic plot of the time mean eddy kinetic energy at the surface. The eddy kinetic energy is defined as 1 2 u · u , where the bar denotes the time mean and u the deviation from the time mean velocity vector u (so u( x, t) = u( x) + u ( x, t)).' 13. Lines 152 153: This sentence is unclear. More likely than what? Than in the other model? Than when the particles arent near each other? Since the latter is obvious, I assume you meant the former. A clearer wording might be Nearby particles are more likely to follow similar pathways in R0.1m than in R0.1. However, it is also not clear how this observation is relevant to the experiment at hand. The particles within each analysed cluster are not located close to each other, since they are separated in time. The link needs to be made more explicitly.
Author's response Indeed, the particles are more likely to stay close to each other in R 0.1m compared to R 0.1 .
We added the reference to Fig 5 in the caption of Fig 3 and 4. 15. Fig. 3: Why is cs = 5 used here for panel (c)? It would help to pick a consistent value for all figures showing R1md. Also, Fig. 3d suggests that a value of 3.5 would be best to approximate R0.1, or a value of 2 to approximate R0.1m. This is consistent for both sinking velocities. Could you comment on this?
Author's response The standard deviation of the noise is approximately between 6 √ c s and 60 √ c s (see lines 99-105 of the original manuscript). It implies that at exactly these values of c s , the strength of the noise compares best to the scales of the eddies in R 0.1 and R 0.1m .

Changes in manuscript
We consistently used c s = 2 in the revised manuscript figures. We added a sentence to the manuscript: 'Interestingly, the value of c s that approximates configuration R 0.1 (c s ≈ 3.5) and R 0.1m (c s ≈ 2) best, is the same for both sinking velocities 6 and 25 m day −1 . These values of c s must result in a similar scale of the fluctuations (σ x , σ y ) in configurations R 0.1 and R 0.1m . It also indicates that, given c s , the subgrid-scale parameterisation performance is similar for both sinking velocities. Nevertheless, the particle distributions match better with the reference configuration if the sinking velocity is higher (according to the W d in Fig 2), because a lower particle travel time leads to less spread of the particle trajectories and a lower lateral displacement in all configurations.' 16. Lines 177 179: Maybe I misunderstand what is meant by a more smoothed pattern, but I dont think that smoothness explains this. Rather it is the lowering of the peaks in the Southern Ocean.

Author's response
The lowering of the peaks in the Southern Ocean is what we mean with 'less extreme.'

Changes in manuscript
We changed these sentences: 'The average lateral displacement becomes globally more smoothed and less extreme :::::::::: (especially ::: the ::::::::: Southern ::::::: Ocean :::::: peaks :::: are ::::::: lower) if the Smagorinsky diffusion is added to the flow dynamics in R 1md compared to R 1m (Fig 4c). The more smoothed :::: less :::::::: extreme : pattern of the travel distances explains why the globally averaged lateral travel distance is minimal at c s = 0.25 (for w f = 6 :: m ::::: day 1 in Fig 4d), and not at c s = 0.' 17. Lines 189 190: See comment #13 it is not clear that this is relevant. The explanation might be instead that particles dont travel as far and hence dont spread out as much. That would be supported by Fig. 5b. This, in turn, is a result of the velocities generally having lower peaks due to longer temporal averaging.

Author's response
These lines refer to the non-eddying configuration, not to the eddying configuration with monthly averaged model output. Fig 5 does not include any information on the non-eddying configurations, only on the difference between the reference configuration with the eddying configuration with monthly model output.

Changes in manuscript
To make clear that these lines are about the differences between the eddying reference configuration R 0.1 and the non-eddying configuration, we add 'non-eddying' before R 1m .
18. Lines 193 194: The two panels of Fig. S1 should opposite behavior in a very similar location. In one case, the blue particles are closer together than the red ones, in the other they are farther apart. This needs to be addressed in the text.

Changes in manuscript
It is now addressed in the text: 'This could result in notably different back-tracked particle distributions, especially if the shear of the flow field is high (see for instance the location 45.5 • S, 39.5 • E on planktondrift.science.uu.nl or Fig S1 for two similar locations with opposite behaviour).' 19. Lines 197 200: This is not entirely true. Fig. 4 shows that in some regions (mainly the Southern Ocean), individual trajectories disperse more in the absence of eddies.

Author's response
With 'disperse less' we meant here that the trajectories disperse less in comparison with each other. In Fig 4, one can see that the particles travel further from their release location in some areas (such as the Southern Ocean), yet they originate from surface locations which are close to each other (i.e. the particle distribution consists of particles which are located close to each other; as can be seen in Fig 3b).

Changes in manuscript
We changed 'disperse less' into 'spread less.' 20. Fig. 6: It seems that R1m appears mostly in black rather than yellow, because only the outlines of the markers are visible.

Author's response
Agreed.

Changes in manuscript
The dots are made entirely black.
21. Lines 242 247: These conclusions are overstated. I am not convinced that the mean flow field in coarser models is wrong in most areas (line 242). At a minimum, this claim needs a reference. In addition, in Fig. 3 and Fig. 4, it looks like the diffusion is sufficient in most areas! So, I dont understand the justification for the claim in lines 246 247. Lastly, these impacts seem to be highly dependent on the sinking velocity, which needs to be acknowledged.

Author's response
The lack of low resolution models to represent the time-mean flow was discussed in the introduction of the manuscript, and supported with several references (e.g. [9,2,8]). Fig 3 and 4 show that the addition of diffusion in the low resolution model can increase the surface area of the particle distributions, up to a point where they match well with the reference high-resolution configuration. However, the Wasserstein distance of Fig 2 shows that the particle distributions are still distant from the reference configuration. This result highlights that the coarse models do show a different surface origin location compared to the high resolution models. The differences of the particle distributions between the configurations are indeed influenced by the sinking velocities. However, this does not imply that the effects of eddies is easier to capture at a higher sinking velocity. In contrast, a higher sinking speed leads to a lower particle travel time in all configurations. As a consequence, the particle distributions will be less spread out and the lateral travel distance will be lower, which is easier to represent.
22. Lines 248 253: Will eddying (i.e., a resolution of 10 km) suffice? Why wouldnt we expect this effect to continue as the model resolution is increased and eddies of smaller and smaller size are being resolved?
Author's response Indeed, a comparison of the back-tracking analysis with Eulerian models of higher horizontal resolution is interesting. Especially, it can be interesting for the modelling of sinking particles in the present-day. The two model resolutions (1 • and 0.1 • horizontally) in this manuscript are the highest resolutions that are typically used globally in palaeo models and present-day, respectively. Due to long spin-up times, palaeo simulations at 10km resolution do not exist yet, and palaeo simulations at 2km resolution require even more computational time.

Changes in manuscript
The paragraph which is added to the discussion addresses this point (see comment 2 for the specific text).
23. Lines 256 257: What does it mean that these models do not produce a flow field that is representative? I would imagine that paleo modelers would consider their models to represent the time period being modeled.

Author's response
These sentences should indeed be rephrased. These models are to some extend representative to the time period, but they do not have the eddying characteristics and do not solve the mean flow correctly as is discussed.