Review of "Islet: Interpolation semi-Lagrangian element-based transport" by Andrew Bradley, Peter Bosler, and Oksana Guba.

The paper presents a new set of basis functions for semi-Lagrangian advection method in spectral element models. The presented method is novel and its performance in atmospheric test cases is encouraging. The paper therefore warrants to be published. The main drawback of the paper is that it is quite tedious to read. This paper essentially presents two things: 1) a novel optmized basis functions for 1D semi-Lagrangian advection schemes, and 2) its implementation in an atmospheric model using three different (sub-element) grids. Both of these are quite complex topics, and their discussion is intertwined (e.g. in section 2) and the reader is easily lost in details. The overall presentation should be improved before the paper can be accepted for publication. Considering the amount of work, I recommend a major review.

Major comments:

For clarity, I suggest that you clearly define the two considered problems, 1D SL advection method, and the implementation in 3-grid atmospheric models, already in the introduction. The introduction is now quite short and actually does not mention many important concepts relevant for the paper.

Section 2 is rather difficult to follow, consider revising. I suggest to start by defining the 1D discretization with N elements, the interpolant functions within each element, and the properties of the interpolant functions (e.g. basis functions, continuity and symmetry). Presently, the interpolant functions first appear only in section 3.1. A figure could clarify the concepts, including the source and target elements/nodes. I would also define a symbol for the basis functions themselves, instead of using L from section 3.1 (L is the interpolant function itself). The discussion of stability becomes comprehensible only after the discretization has been introduced.

As I understand, the 3 axioms of advection methods are: 1) global conservation, 2) preservation of constant tracers (sometimes called local conservation, or tracer consistency), and 3) monotonicity (i.e. no spurious overshoots appear). These properties apply to both the SL tracer advection scheme and the remap operators between grids. In section 1.2 these concepts seem to be mixed and referred to by different names (plausibly due to historical reasons; the so called "property preservation" is just a combination of 1, 2, and 3). Consider revising.

Figures 20 and 21 shows larger diffusion for shorter time steps. Could you elaborate on this? Is the traditional SL property that longer time steps reduce diffusion (while the solution can degrade in other metrics)? How would one choose the right time step in practice? In the case of Eulerian transport it is easy: take the maximum stable one, and you are guaranteed to satisfy all the necessary properties.

Minor comments:

Throughout the manuscript the authors use the terms "mixing ratio" and "tracer" for the advected quantity q_i, seemingly interchangeably. For the sake of clarity I would prefer just to use "tracer".

l35: what is the definition of "local conservation" here?

l155: "This structure arises as follows. Consider a continuous discretization using a nodal n_p -basis, n_p = d+1, with n_p the number

of nodes. The grid has N elements. Each row of the space-time matrix corresponds to a target node."

This description is too brief to be understandable, please elaborate. This is my interpretation of the discretization: The 1D domain is divided into N elements. The solution in each element is approximated by a continuous function, defined by n_p basis functions. Thus a function f in an element e can be written as f_e(x) = sum_i=1^{n_p} f_i psi_i(x) where psi_i and f_i denote the i-th basis function and its corresponding coefficient. Each basis function is associated with a node x_i within the element; The basis is Lagrangian (a.k.a. nodal), i.e. f_i(x_i) = 1. Furthermore, the discretization of the function is continous across element interfaces, implying that the neighboring elements share a (exactly one?) basis function. Furthermore, the basis is assumed to be symmetric about the center point of the element.

l156: "Each row of the space-time matrix corresponds to a target node."

You should define the space-time matrix for this statement to be comprehensible.

l177: Have you defined the basis to be symmetric somewhere?

l195: "L provides a basis for degree-d polynomials." I would say that the basis functions are the \Pi_i functions defined in the equation of L; L itself is the interpolant function defined by the basis and the specific nodal values y(i).

l195: "These are supported by n = d + 1 points, each an element in the n-vector xn." To be consistent with the literature I would use the term "node" instead of "point".

l204: "Given a departure point x" These properties define the interpolant functions, thus there's no need to say that x is a departure point, it can be any point within the element.

l208: I think this constraint is equivalent to saying that the basis must Lagrangian or nodal?

l246: only here you define a d-degree polynomial. This would be useful already in section 2.

l257: what is Runge's phenomenon? what is Lebesque constant? help the reader to understand the rationale behind your work.

Section 3.6: the description of the search algorithm is quite technical and could perhaps be moved to the appendix; it is not necessary to follow the main storyline of the paper.

Section 4: mention TTPR already in the introduction as it seems to be relevant for the entire Islet method.

l439: earlier tracer tendencies were denoted by f_i \Delta t

l442: "immediate element neighbors" Are these neighbors that share an edge or vertex?

l464: "In contrast, ..." Meaning unclear, please revise.

l507: What is "tracer density"? Is it just \rho? Then, for clarity, I'd call it "density"

section 4.3: while computational efficiency is important, I would move this section to the appendix, as it

l570: "is proportional to the number of grid points" Should read: "proportional to the square of number of grid points"

l572: what is a naive h-halo exchange and how does it differ from what is proposed here?

section5, l606: please mention the problem domain (full sphere?) and the equations that are being solved (pure advection on the tracer grid?)

l631: why can you not use the dynamical core to compute the density? the chosen approach seems rather ad-hoc; can you guarantee that density values are realistic?

scaling figures 6-8, 11-16: for easier readability do not use red line twice for np=6 and np=12. x axis label is missing.

l703: "has even more accuracy" Can you quantify this? Is the difference significant?

l711: "For each value ..." unclear sentence, please revise.

l741: 0 in m(0) stands for time t=0?

l780: what is a terminator?

l897: Here you define the basis functions \phi_i. This notation would be useful throughout the manuscript, already in Sect 2 and 3.1.

abstract and intro: you define an abbreviation "dycore". I would omit it as it does not really save space, and it is not used frequently in the paper.

Typos:

l745: "shows the results"

l815: SYPD numbers printed above ...

This manuscript  describes an interpolation-based semi-Lagrangian (SL) method for the transport problem on spectral-element (SE) domains.  The SL  transport  schemes are  widely used  for multi-tracer transport  in atmospheric models due to their  accuracy and computational efficiency.  The classical SL method employs interpolation at the upstream locations of the backward trajectories to estimate the advecting scalar values at the new time level. However, such an approach is not conservative per se, for practical applications an arbitrary procedure  known as the “mass-fixing” usually employed for global conservation — which may have an adverse effect for climate-scale (long term) integration due to the local mass drifting. On the other hand, a finite-volume formulation of  the SL method is conservative by design, where the upstream interpolation over the Lagrangian element is replaced by integration constrained to be locally  (hence globally) mass conservative. The conservative data transfer from regular Eulerian grid to the deformed Lagrangian grid  often referred to as the remapping (re-zoning),  a limiter or shape-preserving scheme is usually employed for physically realizable solutions.  A wide body of literature is available for both conservative and classical SL methods. 

Implementation of  conservative SL method on spherical domains tiled with high-order spectral-elements are very challenging. Authors have proposed an interpolation-based SL method Islet for the SE discretization.  Instead of using the “unstable” native high-order  interpolator (basis function) they have devised a cumbersome numerical procedure   which employs an alternative grid system within each  spectral element, adding another layer of complexity. The Islet method is not conservative, nevertheless, the global conservation is achieved by mass-fixing. The authors argue that the Islet method can handle tracer transport as well as  the remapping between  physics  & dynamics grids, and incorporate shape-preservation filters. 

The manuscript is very long, the Islet interpolation as described by the authors is extremely complex.   Authors failed to explain the core interpolation algorithm with clarity, there are many statements in the manuscript which leads to ambiguity.  The numerical analysis part is very intense  maybe more suitable for a  computational math journal  (e.g., SIAM / JCP)  than  the  GMD. The subject covered could be split into a two-part paper, one describing the basic algorithm and analysis with more details and rigor, and the second part for implementation and validation with standard tests. This would be helpful for better reading. The current  manuscript is written in an awkward manner  and is unacceptable for publication.

Recommendation:      Major revision, possibly resubmit  as a two-part manuscript. Authors should address the following questions. 

(1) The stability associated with the SL method is that the deformational Courant number  (Lipschitz condition) should not exceed unity, in plain  language, the trajectories should not cross intersect (see, Staniforth & Cote’s 1992 MWR paper). Is the  cubic ISL method (lines  115-120) unstable due to this condition? Need some explanation. 

(2) The SL  transport scheme can be stabilized using a limiter,  filter  or with an explicit diffusion (see, Ullrich & Norman, QJRMS, 2014). You can use the native high-order SE interpolation (basis function) for the SL transport combined with the limiter which you are already using for the Islet method. It will be interesting to see how the Islet method compares with this simple SL-SE scheme employing 4x4 GLL grid (I guess that is the SE grid choice made for the operational E3SM). 

(3) It is not convincing to have 3 grid systems (physics: FV, dynamics: GLL, transport: tweaked GLL)  in a SE modeling framework. The Fig.5  shows such a grid configuration, and it appears to be very  challenging.  At  a very high (NH)  resolution the data movement is a major issue for an element-based  Galerkin model (DG/SE).   A typical climate model may have O(100) tracers, an additional   tracer grid with more  DOF than the dynamic grid can exacerbate this problem. This will limit the use of Islet scheme, how do you address it? 

(4) With real data you have velocity information only available at the  GLL (dynamics) grid, the way you find   the 2D trajectory information   using the 3D Cartesian coordinates  leads to additional computational overhead when the method is extended to the 3D application (line 470-475). This needs some justification, why not use the spherical (u,v) components or corresponding contravariant vectors?  It is not clear that the maximum eigenvalue required for the interpolation is the tracer data dependent, in that case you have a serious  computational overhead for the  multi-tracer applications, Please clarify!   What is the computational halo  requirement for an  SE  stencil with  NxN GLL points, when the shape preserving limiter is applied? 

(5) What is the special advantage of using  Islet  method? It seems you have introduced a complex numerical method for a relatively simple linear transport problem.  If mass-fixing is the way to go, one could use the RBF-based (Kriging type)  interpolator which provides very accurate solution, and no need for  the expensive search for max eigenvalue etc. 

(6)  The results are looking good, authors should limit the number of figures and make an effort to compare the results  with  that of other high-order element-based schemes. Why the results from your own previous papers  (Bosler et. al. 2019, SIAM J Sci. Computing; Guba et al. 2014, JCP) discussed? These results  should be compared and  the relative merits should be discussed. 

(7) There are many undefined terms (e.g. CAAS) and notations which I am going to list, this should be fixed.