This is my first review of a manuscript by Oelerich et al. titled “An estimate of the eddy diffusivity tensor from observed and simulated Lagrangian trajectories in the Benguela Upwelling System”. The authors use observed and simulated Lagrangian particle trajectories to estimate the diffusion tensor in the Benguela upwelling region. The authors find that the background mean flow significantly impacts the trajectories and reliable diffusion tensor can only be estimated once the mean flow impact is removed. Likewise, the authors find that motions below 0.3 deg (unresolved by the OSCAR product) contribute significantly to the estimated diffusivities. Although some questions remain as the results are based on relatively short timeseries, I think the manuscript is well written, methods are well established, and the results support the conclusions drawn by the authors. Therefore, I only have a few editorial comments and consider this to be a minor revision.

General

The authors use a drifter dataset that they have contributed in collecting. However, these drifters were all deployed around the same time, and I wonder if the Global Drifter Program has data in the region from other years? I think the manuscript is okay even without additional data, but I would imagine that downloading a few additional GDP trajectories and repeating the computation would be quite straightforward and would make the manuscript more robust and likely more influential.

Also, the POP simulation is for a different year than the trajectories/OSCAR data. It is likely that there is a fair bit of interannual variability (see e.g. https://www.science.org/doi/10.1126/sciadv.aav5014), and, therefore, it doesn’t make sense to make any direct comparison between the model results and the observational results. It would increase the robustness of the results, if the authors would, for example, show timeseries of (E)KE or some other meaningful parameter/index in the region, such that the reader gets an idea how ‘normal’ the years 1996/2016 might have been (altimetry derived eddy trajectories might be a source for interesting data – see also the comment below). The OSCAR product is available for both 1996 and 2016, so one could, for example, show the (E)KE in the upwelling region throughout the OSCAR timeseries.

The authors maybe interested to read a recent paper by Zhang and Wolfe https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023MS004004 in terms of timescales. In the future, it could also be interesting to estimate diffusivity from eddy tracks in the Benguela upwelling region.

L225 This might be a naive question, but is K really symmetric by construction? I have worked before only with diffusivity tensor inversion, and then one needs to still separate the symmetric and antisymmetric part of the tensor. So I am wondering that if the authors actually rotate K by the angle alpha (equation 20), are the off-diagonal components 0? It is not obvious to me that the equation 18 would automatically produce a symmetric tensor.

L475 The last two sentences of this paragraph are a bit puzzling. I would probably rather use ‘entrain’ than ‘pull’ (I guess it is still the wind driven upwelling that brings the cold water to the surface, and these waters are then entrained into westward propagating eddies). It is also somewhat well known that the eddies organize into bands (see e.g. citations below) and I wonder if these results also reflect this fact. i.e. there are mixing barriers (fronts) linked to preferred eddy paths.

https://agupubs.onlinelibrary.wiley.com/doi/10.1002/2016JC012348 (recent paper in upwelling region of Chile)

https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2005GL022728 (original paper)

https://agupubs.onlinelibrary.wiley.com/doi/full/10.1002/2014JC010088 (bands in SST)

Figures:

Figs. 4-5.: The labels are confusing, I would suggest ‘U included’ and ‘U subtracted’

I wonder if instead of using the individual trajectories in Fig. 3, it would be more informative to show heat maps/contours of trajectory concentration (bin the count of the trajectories). Especially for POP, there are so many trajectories that it is hard to say if there is an underlying pattern in the saturated parts of the figure.

I would suggest log scaling for all the plots with diffusivity tensor terms (in python one can use symlog to have log scaling for both positive and negative values).

Style:

There are a few occasions with citations in brackets within brackets, like L107. I would suggest removing the inner brackets around the citations in these cases.

L505 This is a rather long sentence. I would suggest breaking it into two.

Analyzing the drifter and numerical Lagrangian data, this manuscript studies the horizontal 2D diffusivity tensor with a particular focus on the scale effect.

Also, it is interesting to know the departure of the single-and pair-diffusivity relation from the theory.

However, I am confused about issues in the manuscript, such as the definition of average.

I will recommend this manuscript's publication after addressing the following questions and comments.

1) Definition of mean. 

Line 113: " ... since it should be independent from the mean flow"

This statement depends on what "mean" means. If a time average to decompose roughly the internal waves and mesoscale eddies then the mesoscale mean flow also introduces pair diffusion. 

Equation (12): I am confused about the definition of U. A five-day average is used?

In parameterizations, Eulerian spatial average should be considered, this manuscript uses time mean of Lagrangian data. They may not have a direct link.

Line 506: only time average is used, right? Then how can we know the spatial information of the deformation radius? This looks like a pure conjecture. 

2) Inhomogeneity influences the diffusivity a lot, especially for the current Lagrangian data, which captures information of different locations at different times. Is there any way to check the homogeneity of the residual velocity field?

3) Figure 3: It is strange to observe that the total-mean does not equal eddy?

4) Line 417: What about the parameterization effect, such as the unavoidable numerical dissipation, that enters the difference between measurement and numerics?

Oelerich et al. examine a longstanding problem, attempting to estimate a subgrid eddy tracer diffusivity tensor using Lagrangian data, both from observations and an eddy-rich model experiment.    They explore Lagrangian particle statistics and their relation to tensor anisotropy in the presence of a background mean flow, including the presence of flow structures and implied links between flow kinematics, transport and mixing.    Overall, this endeavour is worth pursuing and the basic message - to carefully consider the form of the eddy tracer diffusivity tensor in regions of significant mean flow - is worth continued communication to the modelling community.     However, I have some concerns over the clarity and substantive novelty of this message.    I wonder whether the manuscript could be revised to address my major concerns, as I would like to see these results published, if possible.  

Major comments:

1. Novelty:  I would like to see the aspects which are novel, compared to, e.g. Klocker and Abernathey (2014; https://doi.org/10.1175/JPO-D-13-0159.1), emphasized more in this study.    

2. Clarification of the similarities and differences between Eulerian and Lagrangian representations of the eddy tracer diffusivity(*) tensor.    Typically, climate models use an Eulerian or quasi-Lagrangian  (e.g. if moving vertical coordinate) form, but much of the manuscript refers to either a Lagrangian tensor or is not clear whether Eulerian or Lagrangian is referred to.  Similar goes for references to diffusivity, diffusion and similar in the Introduction, e.g. lines 70-95.

3. As has been stated in the other reviews, I am also concerned about the definition of the "mean" flow and associated implications.   This needs further exploration and justification.     Also, the time analysis window mismatch needs addressing, as suggested in the other review, for example.

4. In several places, e.g. line 199, "submesoscale" motions are mentioned.    These usually have a relatively strong vertical velocity component, compared to mesoscale eddies.    Although a 3x3 tensor is mentioned early in the manuscript, only a 2x2 tensor appears to be diagnosed.   How should this be reconciled?   Also, are there any implications for using 2-particle pair statistics with a 3x3 tensor, e.g. that there might be insufficient statistical information to capture flow deformation gradients in 3D or even in 2D?  Surely 2-particle stats only sample deformation rate gradient in one direction, don't they?

I would be happy to review a revised version of this manuscript.

I broadly agree with the comments already made by the other two reviewers.