Please see attached supplement

The manuscript “Consistent Point Data Assimilation in Firedrake and Icepack” documents the implementation in Firedrake of a UFL-compatible point-interpolation operator of a finite element function. The key benefit of the proposed implementation is that it enables the computation of gradients and second-order derivatives using the hybrid automatic differentiation capabilities of dolfin-adjoint. The parallel implementation is achieved by use of PETSc data structures for easily computing collisions between points and elements of a finite element mesh. Although from the software point of view the contribution is remarkable, the overall presentation of the manuscript could be improved.

Major issues:

1. The use of unstructured point data in inverse problems and data assimilation is well established. For example, the hIPPYlib software [1,2,3], an inverse problem Python library based on FEniCS (of which I am a developer), already implements point-wise observations while automating the computations of first- and second-order derivatives by use of FEniCS symbolic differentiation of weak forms. Furthermore, the three examples presented in the paper are related to some applications of hIPPYlib, including the inversion for a diffusion coefficient in the Poisson equation [1], poroelasticity [4], and ice sheet flow [5].
2. In the numerical results, the choice of the interpolation methods to map the discrete observation to the grid is naïve. A better approach would have been to compare with Gaussian process interpolation (and its 2D/3D counterparts). Gaussian processes do account for uncertainty in the measurement and provide a spatially varying estimate of the variance of the interpolator, thus allowing for more consistent handling of uncertainty measurement and of the fact that the interpolation is less reliable far away from observation points.
3. The authors claim that the numerical results demonstrate Bayesian consistency of the proposed interpolation method. Why this claim is true, the numerical results do not support such claim as the authors do not solve a Bayesian inverse problem. In particular,
   1. The authors only compute the MAP point, while the solution of the Bayesian inverse problem requires the ability to characterize (either by sampling or using variational inference the posterior distribution.
   2. The authors use Tikhonov regularization on the spatial gradient of the unknown parameter as the logarithmic of the prior distribution. However, since the inverse of the Laplace operator in two or three spatial dimensions is not a trace-class operator, the Bayesian inverse problem is not well-posed with that choice of prior [6].
   3. The L-curve, Morozov, or cross-validation criterion are commonly used in the deterministic formulation of inverse problems. In the Bayesian formulation, hyper-priors are often used to tune the variance and correlation length of the prior distribution [7].
4. The definition of Bayesian consistency is not interpreted correctly. Below formula (31) the authors claim that as the number of observation points goes to infinity then q_est should converge toward q_true. This is not what Bayesian consistency is. In fact, there are several cases in which it is impossible to recover the true parameter exactly no matter how many observation points are used, simply because there is an intrinsic loss of information in the continuous map from the parameter to the state variable.
5. It is not clear what is the solution approach in the second example. Are the authors employing a Kalman filter or are they solving for a MAP point and then invoking the Laplace approximation to the posterior to estimate the posterior PDF?
6. Given that the main contribution of the paper is the implementation in UFL/Firedrake of a point-wise observation operator, it would be good to numerically demonstrate the scalability of the implementation with respect to the number of observation points and the number of MPI processes for both 2D and 3D geometries.
7. The overall quality of the exposition is not suitable for publication.
   1. There are several inconsistencies in the mathematical notation that make the manuscript difficult to follow. For example, sometimes the authors use J^point and sometimes simply J to denote the point-wise misfit, and similarly sometimes use J^field and sometimes J’ to denote the interpolation-based data fidelity.
   2. Section 2 is a mix between the very basic theory of finite elements and more advanced concepts like dual basis. However, the section lacks of mathematical rigor in the notation and there is an imbalance of content. In particular, the description of the foundation of finite elements is too long and well-known (e.g. Fig 1 could be removed), while the most advanced most relevevant to this work (the dual basis and the Charlet triplet) lack of details.
   3. The presentation of the paper is very informal (e.g. the use of “our”, to refer to “our” finite element space, “our” finite element function, “our” inverse problem) and there are several typos in the mathematical formulas (e.g. there appear to be an extra semicolon in Eq (11) and (12); or inconsistency between \tilde{x} in Eq (13) and the text below where \hat{x} is used.

In consideration of the major issues noted above I do not believe that the manuscript is suitable for publication in the present form.

Also, given the fact that in general the use of point-wise observations is well understood by the inverse problem and data assimilation communities, I am wondering if the manuscript should be more focused on the true (and remarkable) contribution that is the implementation of point-wise observation within UFL/Firedrake. As such, a manuscript focused on the implementation details and algorithms as well as the assessment of its numerical performances may be better suited to a software journal such as SoftwareX or TOMs.

References:

[1] Umberto Villa, Noemi Petra, and Omar Ghattas. "HIPPYlib: an extensible software framework for large-scale inverse problems governed by PDEs: part I: deterministic inversion and linearized Bayesian inference." ACM Transactions on Mathematical Software (TOMS) 47, no. 2 (2021): 1-34.

[2] Kim, Ki-Tae, Umberto Villa, Matthew Parno, Youssef Marzouk, Omar Ghattas, and Noemi Petra. "hIPPYlib-MUQ: A Bayesian inference software framework for integration of data with complex predictive models under uncertainty." ACM Transactions on Mathematical Software (2023).

[3] Villa, Umberto, Noemi Petra, and Omar Ghattas. "hIPPYlib: An extensible software framework for large-scale inverse problems." Journal of Open Source Software 3, no. 30 (2018).

[4] Alghamdi, Amal, Marc A. Hesse, Jingyi Chen, and Omar Ghattas. "Bayesian poroelastic aquifer characterization from InSAR surface deformation data. Part I: Maximum a posteriori estimate." Water Resources Research 56, no. 10 (2020): e2020WR027391.

[5] Olalekan Babaniyi, Ruanui Nicholson, Umberto Villa, and Noémi Petra. "Inferring the basal sliding coefficient field for the Stokes ice sheet model under rheological uncertainty." The Cryosphere 15, no. 4 (2021): 1731-1750.

[6] Andrew M  Stuart "Inverse problems: a Bayesian perspective." Acta numerica 19 (2010): 451-559.

[7] Daniela Calvetti and Erkki Somersalo. "Hypermodels in the Bayesian imaging framework." Inverse Problems 24, no. 3 (2008): 034013.