The authors present and apply a new methodology for data assimilation that includes non-Euclidian metric and has a potential for addressing systematic and random errors. Results from applications to linear advection-diffusion equation and to Lorenz-63 model show that the new system is performing better than several referent DA systems. A limited discussion of computational cost and applicability of the new method is also included. The methodology and mathematical formalism are sound, but the choice of experiments and presentation can be improved to highlight the pros and cons of the new method. I recommend the manuscript to be accepted subject to major changes. My comments are listed below.

General comments:

(1) The inclusion of 3D-Var in the suite of experiments is not clearly justified. 3D-Var is a deterministic/variational method, while the new method is ensemble-based. It is already understood how ensemble and variational methods work and differ. In addition 3D-Var is not very much competitive with Ensemble Kalman Filter (EnKF) or Particle Filter (PF) methods anyway, so showing that the new method outperforms 3D-Var is not saying much. If there is a desire to compare the new method with deterministic methods, then using 4D-Var instead of 3D-Var would be more advisable. Please justify using 3D-Var or exclude it from the experiments.

(2) The linear advection-diffusion and Lorenz-63 applications include different sets of data assimilation (DA) systems: 3D-Var is used in advection-diffusion and EnKF and PF are used in Lorenz-63 model. This does not allow clear understanding how the new system performs relative to other DA algorithm and what is its true capability in these different applications (e.g., linear/diffusion and nonlinear/chaotic). Please include EnKF and PF in linear advection-diffusion application, and 3D-Var in Lorenz-63 application (conditional on the comment-1 above) .

(3) There is no clear understanding if the new method is better suited for handling systematic or random errors. Please add the standalone experiments: (i) random error only and (ii) systematic error only, using the same setup as in the random + systematic error experiment.  This will help understand better the true potential of the new system.

Specific comments:

(4) One should be aware that systematic errors can be addressed within variational and ensemble methods. Your comparison includes other DA methodologies without that capability. Could you please briefly comment on that in the manuscript?

(5) Abstract includes “Eulerian penalization of error in Euclidian space”. However this is not a common terminology which requires some explanation (as done in main text). Please adjust Abstract to avoid this terminology until explained.

(6) p.1-3: In general, the manuscript includes new terminology that I would not consider common. I would prefer to see some basic explanation of uncommon terminology (such as Wasserstein metric/space, Riemannian manifold, entropic regularization, Sinkhorn algorithm), to allow the reader to go through the manuscript without the need to immediately check the references.

(7) p.5, L.111-115: I am not sure that 3D-Var the authors refer to is the same as the commonly used 3D-Var. Aside from using a pre-assigned static background error covariance matrix, standard 3D-Var is not including the full matrix inversion that allows direct solution of the linear problem (e.g., in Eq.(26) of Lorenc (1986) the matrix B*H_tr*R_inv*H=0). Please refrain from using 3D-Var to describe the system you use, or include some clarification.  

(8) p.6, L.146: Is the cost with subscript “q” and superscript “q” a typo? Can they differ?

(9) p.8, L178-180: Can this be explained in somewhat more detail? Otherwise, not clear why Cholesky decomposition is relevant here.

(10) p.11, Eq.(13): How is assured that Kw, Kv are not zero, which is required for the element-wise ratio?

(11) p.13, Eq.(14): A vector dot-product is typically used for advection term. Is this a typo or you are using an element-wise product? Please explain.

(12) p.13, L. 290: How difficult would be to obtain the optimal parameters eta in realistic high-dimensional situation?

This paper introduces an optimal transport framework for updating discrete representations of posterior probability density functions during ensemble data assimilation. This work further provides proof-of-concept assimilation experiments comparing the algorithm the authors introduce (and call EnRDA) to standard DA methods (3D-var, particle filters, ensemble Kalman filters). The authors mention two serious issues that will need to be overcome for EnRDA to become a viable strategy for users in the DA community: 1) the high computational expense (that scales super-cubicly with the ensemble size) associated with computing optimal transport maps and 2) the need for the observation operator to be bijective, e.g. all state dimensions must be observable.

• The reader would benefit from more discussion of the Sinkhorn algorithm and how/why gamma is chosen as well as how eta is chosen in practice. On line 290, you mention that these parameters are set by “trial and error”.  Can you offer the reader more guidance on how to make these choices in practice, even if ad-hoc? Do these or should these parameters vary at different assimilation time steps? On line 290 you have gamma=3 and in the caption of  Figure 5 you have gamma=0.003. Are both correct? Are these for the same DA experiment at different times? Or  different experiments? This should be clarified in the text. In Figure 3 you demonstrate that one order of magnitude change in gamma results in quite different joint distributions with the same prior and observation pdfs. What does a 3 order of magnitude change in gamma do?
• The barriers to wide usage of this approach are quite high, yet if overcome the EnRDA could prove very a powerful DA method. As such I believe these barriers and the research advances needed to overcome them warrant a longer discussion than you offer in section 5. 
• (Typographical mistake — you have two periods ending a sentence on line  398.)