In this manuscript, the authors present spyro, a full waveform inversion finite element solver based on Firedrake. The solver features wavespeed-adapted triangular/tetrahedral meshes, a fully-explicit time stepping scheme based on a mass-lumping technique, and conforming elements of variable orders (up to degree 5 in 2D, and 3 in 3D).

The paper is well-written and fits the scope of this journal. The numerical results also look promising and the proposed solver seems to have several advantages over other FWI methods. I only have a few remarks/questions:

-p3, line 61: spectral triangular elements are known at least up to degree 9 (see, e.g., Mulder 2013: New triangular mass-lumped finite elements of degree six for wave propagation, Cui et al. 2017: High order mass-lumping finite elements on simplexes, and Liu et al. 2017: Higher-order triangular spectral element method with optimized cubature points for seismic wavefield modeling).

-p6, line 161: \sigma_i is not defined here. I suggest to give the definition of \sigma_i and \Psi_i here and explain that p_i and w only need to be computed in the boundary layer.

-Section 2.2: is there a reference for the derivation of the adjoint equations? Especially the boundary conditions given on page 8, line 207 need some explanation. The adjoint problem has 2 boundary conditions while the original problem had only 1.

-p9, line 230: definition of F is missing.

-Section 3.1: is the stiffness matrix assembled and stored or are the computations matrix-free?

-p12, top line: In 3D, alpha is computed taking the cube root?

-p12, definitions A_{n+1}, A_n, A_{n-1}: please double check these definitions. 

--In A_n top-right: should be M_{omega, 1} instead of M_{omega}.

--In A_n bottom-right: should be 0.

--In A_{n-1} second row: second and third term should be swapped.

--In A_{n-1} bottom-right: should have a minus sign.

-p14, equation 36: the definition of G is not really clear. What is its continuous counterpart? Also, the right-hand-side should be a vector. Please give a more precise definition.

-p16, line 396: 32 nodes instead of 50.

-p17, Algorithm 1: step 10 is done via L_BFGS and step 11 via ROL? Or are both used for both steps?

-p17, line 441: equation (36) instead of (19)?

-p18, Figure 5: This figure is rather unclear. Does gradient.py do steps 8+9 of Algorithm 1? This figure also contains several equations, whereas I would expect steps of an algorithm.

-Section 5: it seems that the z-coordinate is always given first. Please explicitly mention this somewhere.

-p22, line 548-550 and Figure 8: for a fixed p, the relation between C and G is actually linear. In Figure 8, straight lines are only expected when using a log-log plot. I would suggest to use a log-scaling for the horizontal axis in Figure 8 or remove the linear fits.

-p26, final paragraph: please double-check the domain sizes.

-p32, Table 3: it would be convenient to also have the final J here as an additional column.

-Section 6: with 2 full pages, the conclusions seem to be overly long. Please try to make it shorter and more concise. A summary of the results and corresponding conclusions should be the main focus. Ideally, this section has one or just a few clear takeaway messages.

-p41, Appendix: The definition of superscripts x_k/x_l are missing. Also, does there need to be a summation over l in the last two equations?