This paper presents the applicability of Isogeometric Analysis for problems in geoscience. More particularly, it is herein considered the modelling of the lithosphere as a thin shell subject to gravitational and buoyancy forces exerted by the topographic features above the lithosphere and the upper layers of the athenosphere, respectively. The paper is overall very well written, easy to read, and the findings are accompanied/supported by numerical evidence. Nevertheless, I would like to have the remarks below addressed in an adequate manner before I can suggest publication of the paper:

1. On page 2 the authors mention "Further features and capabilities of isogeometric analysis presented in this paper are the exact representation of curved domains, the coupling of multiple patches". The coupling of multiple patches is not a feature of IGA but rather a challenge of the method, from which standard finite element methods do not suffer. I would not enlist the coupling of multiple patches as a feature, but I would highlight that patch coupling is essential in IGA for problems of practical relevance,
2. On page 5 the auhors provide the parametric description of the continuum using contravariant coordinates. They refer to these as curvilinear coordinates in general. Since the superscript is used, I would recommend calling these as contravariant and not simply curvilinear coordinates,
3. On Section 2.1.1 Shell Configuration the authors expland in very detail in the variational formulation of the Kirchhoff-Love shell. These equations can be found in many resources, see for instance in Kiendl, Josef, et al. "Isogeometric shell analysis with Kirchhoff–Love elements." Computer methods in applied mechanics and engineering 198.49-52 (2009): 3902-3914. I would recommend reducing this part to the absolute minimum and just citing an appropriate resource,
4. On page 7 the authors mention the following: "In the one-dimensional case with w = w(x), the bending term reduces further, which leads to a fourth-order differential equation for an Euler–Bernoulli beam when considering the strong formulation of the problem without boundary conditions". I think this is an oversimplification of what is in fact going on. The mathematical and mechanical effect on reducing a structural element from a 3d continuum to a beam introduces many additional mechanical effects (torsion, bending, twisting, etc.) that can not be understood as a reduction from a plate (or shell). I do not see a reason why this section is added and I would recommend removing it (see Bauer, A. M., et al. "Nonlinear isogeometric spatial Bernoulli beam." Computer Methods in Applied Mechanics and Engineering 303 (2016): 101-127. for more information),
5. On page 7 the authors mention "We model the lithosphere as a thin elastic plate". Is this a thin elastic plate or a thin shell? The equations the authors laid down in Section 2.1.1 refer to a shell,
6. On page 7 the authors refer to a "Kirchhoff-Love plate". To my knowledge "Kirchhoff-Love" is the designation for a shell, the corresponding plate is called just "Kirchhoff",
7. On page 8 the authors write the following "Instead of working with the actual topographic elevation, we use a mass representation r obtained by taking the mass above the mid-surface of the lithosphere and normalizing it by some depth-independent reference density ϱr.We choose the reference density as the mean rock density from the current depth of the mid-surface to its initial depth and assume that it is homogeneous". This approach is not clear to me. I would like to ask the authors to add some reference if it is standard in the literature? I am not sure what this modelling implies,
8. On page 9 the authors mention the following: "For the full Neumann problem, it can be shown using Korn’s inequality that the variational problem is well-posed, provided that A is a bounded domain with piecewise smooth boundary and the coefficient (ϱm −ϱr)g in the buoyancy term is bounded from below by a positive number.". Is there any reference for this statement? When the authors mention that something can be shown it is advisable to either add a reference or a proof, if there is no reference,
9. On pages 9-10 the authors mention the following: "From a modeling point of view, the results may not reflect the physical reality since the Earth consists of different regimes and tectonic plates that interact with each other in a complex manner. Furthermore, due to the large scale of the simulation, the effects of flexural rigidity will not be visible. Nevertheless, we assume that the entire lithosphere can be modeled as a single spherical shell to showcase the capabilities of isogeometric analysis in numerical simulations on curved domains, especially on a spherical domain". Is this to be understood as an academic study without applicability to geoscience? I am not sure how to interpret this statement as it mentions that the results may not reflect the physical reality and that the study aims to showcase the capabilities of isogeometric analysis in numerical simulations on curved domains,
10. On page 10 the authors mention the following: "to both parametrize the domain and construct finite element approximations
    
    of solutions to the equations.". I would recommend rephrasing this to "to both parametrize the domain and construct finite element approximations of solutions to the corresponding partial differential equations."
11. On page 11 the authors mention the following: "In the following, we consider splines on the unit interval [0, 1] with an open knot sequence, i.e., the first and last knot values have multiplicity p+1. Then the knot sequence has the form". I would like to ask the authors to mention the interpolation effect that is enforced when the first and last know values have the prescribed multiplicity. Why are the authors restricting themselves to this choice and what is the implication of doing so?
12. In the caption of Figure 5 the authors write "bivariate quadratic ...", I would recommend replacing "quadratic" with "biquadratic",
13. Regarding Section 3.1.2 NURBS basis functions, I feel that it is spent quite a lot of real estate on detailing the B-Spline and NURBS basis functions, information that can be found easily in many other resources. Moreover, the authors do not cite these other resources in an adequate manner, but only Pielg and Tiller. I would like to ask the authors to add more references as needed,
14. On page 16 the authors mention the following: "We can stack the control points on top of each other". I would recommend rephrasing this to "the control point coordinates are organized in vectors" That sound more appropriate for a scientific publication,
15. On the same page the authors mention "There are various methods that can be employed to achieve this." when referring to NURBS multipatches. There are so many publications and studies on this matter, I would like to ask the authors to add adequate references. The reference list is herein quite poor,
16. On the same page the authors mention the following: "by replacing shape functions at the boundary of each patch that coincides
    
    with the boundary of another patch with interface functions that span over multiple patches". Could the authors clarify whether this is kind of a bending strip? See Kiendl et al. 2010
17. On the same page the authors mention the following: "Another approach is to stitch shape functions at the patch interfaces together by imposing the C1 condition". Could the authors clarifys whether it is herein meant that the C^0 continuity is already assumed (meaning that the patches are conforming along their interface) and that the authors are just constraining the second row of control points in such a way as to achieve C^1 (or G^1) continuity? It is still unclear to me how the C^1 continuity is enforced (is this a Penalty, Lagrange Multipliers, Mortar, or any other approach?),
18. On the same page the authors mention the following: "Aside from that, the computation of the null space
    
    corresponding to the C1 constraint is generally a difficult task numerically.". Do the authors herein mean that the coupled system becomes non-convex?
19. Yet on the same page, the authors mention the following: "In this work, we will mainly consider planar domains that result from joining convex quadrilaterals along the sides. It has been shown that the class of bilinear G1 parametrizations is analysis-suitable, so that optimal convergence can be achieved in this setting". I understood that the authors would model the whole earth's lithosphere as a spherical shell. How can the domains be considered planar in such a case?
20. On page 20 the authors mention the following: "It differs from the single-patch result due to the missing data outside of the simulation domain that are replaced by Neumann boundary conditions.". I think it is natural that the results naturally differ. Could the authors clarify whether the point of this comparison is to show that one can obtain satisfactory results also when using less DOFs? If not, could they please clarify what the exact point of this discussion is?
21. On page 21 the authors mention the following: "When topographic load is the sought parameter, its initial value is set to 1km". Could the authors clarify why is the value of the topographic load measured in Km? Km should measure distances, right?
22. I tried reproducing the results, but I stumbled upon some difficulties with the provided open-source repositories. I am able to run script script_central_java.m provided that all necessary third-party software is installed, but I am unable to execute script script_spherical_earth.m because I receive an error message at line 156 of file op_KL_bending_stress.m from GeoPDEs, namely: "Error using .^ Invalid data type. Argument must be numeric, char, or logical." I tried to debug the issue, but it is probably due to incompatibility of the version of the GeoPDEs library with the one used to produce the results in this study. Could the authors also upload your code also on GitHub and provide clear instructions on the versions of the third-party libraries needed to reproduce the results shown in this study?
23. Otherwise, the results section is really great and I am impressed with the results shown in Figure 12. Great work, congratulations!

 This paper describes a refinement to finite element modeling of thin-plate isostasy that incorporates isogeometric representation of surfaces such as topography and Moho flexure. I am unfamiliar with the existing literature on spherical finite element approaches to modeling isostasy (though presumably there are many such), and without a head-to-head comparison of results from those other methods, it is unclear to me how much of an improvement is afforded by the isogeometric approach to isostatic modeling presented here relative to (say) representing topographic and internal loads or the Moho response as stress fields acting on the plate. The requirement of fewer degrees of freedom sounds promising though (although it is not rigorously explored here), even if it is unclear that a more realistic representation of surface curvature affords other advantages.   Effective elastic thickness (Te) estimation may benefit from improved finite element approaches in a couple of ways: (1) Finite-element methods are commonly used to generate synthetic data with known property variations as a means of testing inverse methods (e.g., Kirby, “Spectral methods for the estimation of the effective elastic thickness of the lithosphere”), and (2) The inverse methods currently rely on the linearity of flexural response afforded by uniform-rigidity spherical harmonic-domain approaches, which of course introduces some error in the elastic plate response when Te is varying. Assessments of inverse methods using synthetic data indicate that the errors are acceptably small for purposes of mapping Te variation using spatially-localized power-spectral methods, but the assumption of locally-uniform Te likely imposes limits on the resolution/bias properties of such methods. The modeling approach described here might be useful for both synthetic and inverse applications, if it were extended to permit both topographic and internal loading of the lithosphere, and to calculate the resulting gravitational perturbations.

With that said though, the application of the method to data from Europe (and elsewhere) has the potential to confuse unsophisticated readers (those who are not deeply-steeped in isostatic inversion methods). I can anticipate some readers believing that (e.g.) the Te or reference density estimates presented in Figure 11 have some bearing on reality, when really they are just parametric artefacts of imposing a single-input/single-output assumption on an inherently multiple-input/multiple-output problem (as first described in Forsyth, JGR 1985). The paper will be much more valuable to readers if that is clearly communicated in the abstract and conclusions, and if results given in §5 are caveated as a demonstration of application of the numerical technique, recognizing that additional enhancements to the modeling approach must be made (e.g., calculating gravity, incorporating additional seismic/density fields) before meaningful estimates of Te may result.

The additional comments below expand a bit on this key point, and I hope will be useful/constructive in helping the authors to clarify that the results given in the paper are useful for demonstrating the method rather than bearing on physical properties.

Line 18: “Globally C1 finite element spaces…” The audience for EGUsphere is not mostly made up of mathematicians, to my knowledge, so it might be worthwhile to define terms like C1 upon introduction…

Line 40: There has been significant change to European Moho maps since Grad et al. (2009), with the addition of hundreds (or thousands) of seismic stations and technical advances like ambient-noise Rayleigh wave imaging. Crust 1.0 (Laske et al. 2013) would be better, and even this has been improved on in more recent global models (e.g. SPiRaL, Simmons et al., GJI 2021). Given the limited utility of the simplified version of isostasy in the modeling performed in this paper, it has little impact: i.e., the parameter estimates derived in section 5 have little bearing on reality owing to oversimplification in the assumptions of load-forcing and response, for reasons described in the comment below regarding lines 163-165, so it doesn’t matter whether the data used to constrain the model is the most accurate available. In future efforts though, assuming the method is extended to accommodate loading at the Moho, internal density variations and gravity prediction, it may be worth incorporating more recent estimates of data fields.

Lines 53-56: It’s worth being a bit careful in the wording here, as the notion of an “elastic lithosphere” is a red flag that can raises suspicions in the lithosphere community. The bending-stress moments that are supported by the nebulously-defined “lithosphere” are actually dynamically-maintained flow- and frictionally-supported stresses (hence the term “effective” in “effective elastic thickness”, see e.g. McNutt & Menard JGR 1978, McNutt JGR 1984). The thin elastic plate approximation is useful solely because the bending moments supported in a thin elastic plate can be physically related to the bending moments supported in the real (dynamically-supported) lithosphere, and the two observationally appear the same for identical bending moments (e.g., Willett et al. Nature 1985; Burov et al. JGR 1995; Brown & Phillips JGR 2000). This enables valuable constraints to be placed on unknown rheological properties.

Lines 163-165: This treatment of the Moho as a surface of flexural displacement is problematic in this context. The Moho is a compositional boundary separating silica-rich (“crustal”) from ultramafic (“mantle”) rocks, and the thickness of the crust is determined by a combination of surface processes (erosion, deposition, volcanic extrusion, faulting), intrusive/magmatic processes (which are of order 10-50X greater volume than extrusions), tectonic strain, and lower-crustal flow. Each of these will have different transfer functions relating surface elevation to Moho deflection at flexural isostatic wavelengths, leading to Forsyth’s (JGR 1985) conceptualization of “surface” and internal or “subsurface” loads acting on the lithosphere. (The Moho can approximate a flexural deflection in response to treatment of topography as a load at very long wavelengths where the transfer functions approach those of Airy isostasy… But only if other buoyancy phenomena such as crustal density variations, thermal variations, or variable normal stress forcing by asthenospheric buoyancy are negligible contributors to vertical stress.) The practical consequence of this is that any estimate of Te arrived at by assuming the Moho is a flexural response to topographic loading will be meaningless for purposes of understanding long-term rheological support of stress by the lithosphere (i.e., the bending moment implied by the resulting estimate of Te will have no knowable relationship to the true bending moment). This is illustrated for example by comparing the map of Te in Europe (Figure 11a) with any well-tested inverse model incorporating a solution for surface and internal loading (e.g., Pérez-Gussinyé & Watts, Nature 2005), and by recognizing that the low Te’s in Figure 11 are incompatible with a reasonable expectation of rheological controls on long-term stress-support in the region (consistent with the earlier finding by Forsyth, JGR 1985, that inversion using an assumption of topographic loading only yields extreme underestimates of Te).

Figures (generally): Matlab has some (albeit limited) geographical mapping capabilities in its figure-making, and for those figures in which geographical regions are represented it would be very helpful to include (at the very least) the coastlines to help orient the reader. Flipping back and forth between the topography shown in Figure 10 and the maps in Figure 11 to orient oneself, for example, is cumbersome.