Review of "Intercomparisons of five ocean particle tracking software packages" by Xiong and MacCready

In this manuscript, the authors present a new particle tracking code and compare it to three existing offline codes, one online code, and an online dye simulation. They focus on two regional domains on the west coast of the US. They show that the new code compares well in terms of accuracy and efficiency.

The manuscript is generally well-written. However, I have some significant concerns that prevent me from recommending its publication in GMD. Most importantly:

1. It is unclear what this new Lagrangian tracking code adds. What is the advantage of the new code over the previous codes? What does it add/improve on the other codes? That could be much more explicit. 

2. The title is far too general, and suggests a much wider scope than that the manuscript can deliver. It's therefore not appropriate for this specific manuscript.

3. It's not at all clear why the LiveOcean and Hood Canal are such good testcases for Lagrangian models. I would have expected a much more thorough discussion of why these are specifically suited. Now, it seems as if the authors had these models lying around and decided to do the comparison; instead of selecting an optimal case for the comparison.

4. It's unclear why some of the choices for the tracker code have been made. E.g., why does it employ nearest neighbour interpolation? That is not very customary for Lagrangian codes. And why then also use 4th order Runge-Kutta integration? Why aim for such high accuracy in time, when spatial accuracy is low?

5. The argument for smoothing the AK_s diffusivity field is unclear; what is the advantage of this?

6. The discussion of the Well mixed condition test in section 2.1.1 could be more elaborate. What is the equation that is tested. Why? How is e.g. the non-significant range defined in Figure 2?

7. One of the most difficulty things to do for Lagrangian codes is boundary conditions near land, and avoiding stuck particles. While the strategies of each code is listed in table 1, there is no discussion of how well the tracker code performs near boundaries. This would be important information for potential users, especially in domains like the Hood Canal where there is so complicated topography.

8. Why is there no comparison to online floats or dye in the LiveOcean domain of figure 6 and 7? For a complete picture, that would be useful here too.

Other minor comments:

- line 9: Make explicit that these numbers (200m and 1000m) are the resolution and not the domain sizes)

- line 15: This sentence is very vague; please rephrase in terms of conclusions/outcomes

- line 26: Explain why offline tracking is more frequently applied?

- Line 35: There are some recent articles that compare different tracking codes in e.g. the Agulhas: https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2019JC015753

- line 38: what is meant here with performance? Speed? Memory? I/O? Accuracy? Reproducibility?

- line 40: Is LiveOcean really 'well-established'? What does that even mean, when it comes from the developers of the model?

- line 127: 'studied' instead of 'studies'

- line 145: what was the convenience why the seeding strategy was different for the online particles?

- line 209: The point that dye spreads faster than Lagrangian particles is not new, and could be related to Markovian dynamics (I.e. dye that enters a grid cell from one side can leave it on the other side within a timestep)?

- line 256: 'growing differences in location' is slightly awkward phrasing?

- line 278: is the laptop the same as the Apple M1 Pro?

- line 284: why does LTRANS scale so poorly for large numbers of particles? This is very surprising for a Fortran code?

The authors present a Lagrangian particle tracking tool, Tracker, and compare it to several other existing particle tracking codes, as well as an Eulerian transport model run as an integrated part of ROMS. In general, I find comparisons like this to be very interesting and useful. In my personal opinion, the field of Lagrangian transport modelling would probably benefit from increased attention to details of implementation and comparison between codes. Hence, the topic should (in my opinion) be of interest to the readers of GMD. However, I find that the present manuscript is a bit lacking in some aspects, and my recommendation is that it should be reconsidered after revisions. 

The abstract states that Tracker is introduced, and also compared to other models. However, there are also references to previous papers that use Tracker (Brasseale & MacCready, 2021 and Stone et al., 2022), even though these papers do not seem to explicitly use the name "Tracker" (from a quick search of the documents). If the current manuscript is the definitive introduction of Tracker, it should in my opinion contain additional details on the implementation to properly document the model.

For example, are the u and v components of the current assumed to be on separate grid points, and interpolated independently? How about the vertical current component? Is variable horizontal diffusivity supported? It is stated that vertical diffusion uses reflection at the boundaries, but what about edge cases where the particle is so far outside one boundary that reflection would take it outside the other boundary (can happen with non-zero probability, since the random walk uses Gaussian numbers).

I would also say that a description of the implementation would be useful for a paper introducing a new model, even if one would hope that the actual performance in terms of accuracy is independent of implementation detail. What python libraries are used? How are the particles stored (simple arrays or custom objects?) Are particle positions recorded in lon-lat or in x-y coordinates in meters? If lon-lat, then how are displacements in meters converted to lon-lat? Are any assumptions made (and hard-coded?) about e.g. the radius of the Earth?

Finally, I would be interested to see a bit more discussion of some of the choices that were made in the implementation. For example, why combine 4th-order Runge-Kutta with nearest-neighbour interpolation? 4th-order Runge-Kutta has requirements when it comes to continuous derivatives of the velocity field, which are not met by nearest-neighbour. Based on earlier work I have done, I would guess that you could get a better ratio of accuracy to computational effort with e.g. 2nd-order Runge-Kutta if you use nearest neigbour interpolation (see, e.g, https://doi.org/10.5194/gmd-13-5935-2020). This might not be very important in practice, though, so feel free to ignore.

I would also like to see more discussion on the vertical diffusivity implementation. Why was a 3-point Hanning window chosen? And it says that this is applied in the calculation of the vertical derivative, but does that mean that the smoothing is not applied for the diffusivity values themselves? I think these values should be chosen consistently, otherwise you might have trouble with the well-mixed condition. Did you do any testing of this point? Note also that the well-mixed condition is only a neccessary condition, not a sufficient one. A perhaps more stringent test of the implementation would be to compare to a dedicated 1-D solver of the diffusion equation with variable diffusivity. The comparison to the dye in Fig. 5 is good, but since this is a 3D case with comparison only in the vertical it is a bit hard to reason about the cause of the discrepancy.

The second point of the paper is a comparison of different particle models and a dye study. Here, I would have liked to see a bit more discussion and investigation of details. For example, the paper states that "a very good inter-model match was achieved", with reference to Fig. 3. However, I would say that there is something odd in Fig. 3, as the trajectories of the centers of mass already from the start appear to separate very fast. This should be straightforward to investigate further, for example by running a single timestep with a single particle and no diffusivity, and checking if all the models move the particle the same distance and direction. Certainly those models that use the same interpolation and numerical integrator should have exactly identical behaviour in the case of no diffusivity.

It is further stated that "vertical distributions of particles (Figure 5) exhibit similar evolutions among all tracking codes", whereas I would say that the distributions are quite different. In particular, Tracker shows quite a spikey distribution, and also the dye study has a lot of spikes in the distribution. This last point was particularly surprising to me, as I would have expected an Eulerian diffusion solver to produce a smoother concentration field than a particle based model. Of course, as this is 3D with advection, that makes it a bit harder to reason about, but this could be investigated. In any case, I would say that it is clear from the left panel of Fig. 5 that these models are _not_ equivalent, or at least not run with equivalent setups.

Looking at Figs. 6 and 8, I am a bit surprised about the large vertical fluctuations in position, particularly in the mid-watercolumn release. Panel g in Fig 6 shows that with Tracker, the center of mass moves more than 60 meters downwards in less than 12 hours, which is a much larger displacement than with any of the other models. I'm also curious about what the mechanism behind this transport is. How large is the vertical current component? How stratified is the water column? Discussion on this point would be appreciated.

Finally, I think it would be good in the interest of reproducibility to provide the actual setups used for the different models, perhaps in a separate github repo for the paper. That would make it much easier for others who might want to look into the comparison.