Introduction

We thank Professor Steinberger for his commentary and suggestions, which follow below. Before responding to each of his comments in turn we feel it is worthwhile discussing the rationale for this paper in general terms.

We would like to emphasise that our goal with this paper is to establish, in one place, sensitivities of calculated surface deflections predicted by simulations of mantle convection that are parametrised with popular modelling choices. It is not to identify an optimal, e.g. most Earth-like, simulation. We readily acknowledge here, and in the introduction and discussion of the manuscript, that various insights into sensitivity of calculated surface deflections already exist in the literature and within the modelling community. However, the insights are dispersed and have been developed using different numeric and analytic approaches and code. Consequently, from an observationalist's perspective, it is often unclear whether variations in misfit between model predictions and observations of, for example, uplift and subsidence, relate to assumed mantle properties (e.g., density and viscosity) or to modelling choices (e.g. the approach used to solve the governing equations). We suggest that establishing impact of modelling choices benefits from development of self-consistent and systematic testing, such as that presented in this contribution. These tests include exploration of the impact of solving the equations of motion analytically and numerically, gravitation, viscosity and density. We are not aware of another (singular) paper that contains such systematic tests. Therefore, we suggest that it is novel, and will be of interest to those working to relate surface observations to the history of mantle convection, such as ourselves.

In summary, we would like to make it clear here and in the revised manuscript that our goal (in this paper) is not to identify an optimal geodynamic model. Instead it is an attempt to establish a framework to assess the consequences of popular modelling choices on calculated surface deflections in a self-consistent, systematic, way.

Response to individual comments  (reviewer comments are bold text)

I don't think there is anything principally wrong with this paper, but this is also where the problem starts, because the methods used are really decades-old and I frankly don't see much novelty in what is shown here. I mean, you can obviously always find a new combination of modelling assumptions that has not been previously used before, or comparison that has not been previously done, or analyze results in a slightly different way as done before, but that doesn't really advance our knowledge unless it is geared towards addressing a specific question. 

We entirely agree that many of the methodologies we examined are established. Our goal with this paper, as stated in the manuscript and above, is to quantify the impact of widely made methodological choices on predicted surface deflections (e.g. approaches to solving the equations of motion, assumptions about gravitation, viscosity, etc.). As far as we are aware there is no (singular) paper that assesses the impact of the various modelling choices we examine on surface deflections in the systematic way that we do. We think that doing so is useful, especially, for those of us trying to understand how observations of vertical lithospheric motions can be used to assess predictions from whole mantle simulations.

The specific question we address is ‘how do widely used modelling choices impact calculated surface deflections?’, which we will endeavour to make clearer in the revised paper.

What is more, some of the modeling assumptions are really questionable or suboptimal, and I think more realistic models have been run before. For example, radial viscosity structure - overall viscosities are much too high to be compatible with the Haskell average for postglacial rebound. Of course, this is acknowledged by the authors, but others have previously done computations with more realistic viscosities. Some more realistic viscosity profiles are shown as dashed lines in the figures, they are not used in the computations. Also, lateral viscosity variations are much smaller than would be realistic, I think, and also much smaller than what has been previously used, so I am not sure how much is gained compared to ignoring them altogether, as with the spectral method, and whether it is all that meaningful to say lateral viscosity variations lead to <25% changes, if they are much too small.

We largely agree with the reviewer’s comment. However, we feel that it is somewhat beside the point of this paper, which, for clarity, is not to identify the most Earth-like simulation of mantle convection.  Instead, we seek to establish the impact of systematically changing, for instance, how viscosity is parameterised on calculated surface deflections. We stress that we do not advocate for a specific viscosity or aim to identify the most realistic model(s) in this paper. Nonetheless, we suggest that the viscosities tested in this paper are reasonably Earth-like and we include tests of viscosities that are broadly like those indicated by the dashed lines in Figures 2c and 8a. Crucially, for this paper, we compare the impact on surface deflections in a systematic way and are explicit about the viscosities used. In that sense we are of the view that the results provide meaningful assessment of the impact of changing viscosity on calculated surface deflections, and thus we fulfil our goals here.

Thermal expansivity is assumed constant whereas, in reality, it is probably strongly decreasing with depth.

To aid comparisons between numeric and analytic (propagator matrix) solutions of the equations of motion we assume incompressible convection in the models tested. As far as we are aware, constant thermal expansivity is the only choice that is self-consistent with the Boussinesq approximation of incompressible convection (see e.g., Ricard et al., 2022, GJI). It is worth noting that differences between compressible and incompressible mantle simulations are small, even in subduction zones (e.g. Lee and King, 2009, G-Cubed).  Nonetheless, we think that it would be interesting to compare surface deflections predicted by incompressible models and compressible models with depth-dependent thermal expansivity. However, we feel that doing so is beyond the scope of this paper.

The radially varying gravity, although it was computed for the radial density distribution in the model, does not correspond to the actual radial gravity, which is more close to constant.

We are unsure what the reviewer would like us to do here. Perhaps he is thinking of gravitational attraction according to PREM (etc.), which has gravity relatively constant through most of the mantle. We acknowledge that, in contrast, gravitational acceleration in the numerical models, which were computed from the radial density distributions as the reviewer says, decrease from the surface through the upper mantle. They increase, similarly to PREM, through the lower mantle to the core-mantle boundary (see Figure 4b in our main manuscript).   We note that we also tested simulations in which gravity is parameterised to be constant, and thus test models that likely span the range of actual gravitational acceleration, which we think is satisfactory for our purposes here. 

The only way how the computations are evaluated is by computing surface deflections. Others who have run these kind of models previously, have, for example also compared model geoid with observations, or the predicted density structure with tomography. None of this is done here, so I find it questionable how realistic the model is. Even the predicted topography patterns are not compared to any observations, except in terms of spectral power, and here, just as by looking at the resulting maps, it is clear that the predicted topography amplitudes are much too high to be compatible with observations, and no suitable efforts are being made to resolve this discrepancy - an issue that has, for example, been addressed in previous work of co-author Richards. The introduction talks about the fit between predictions and observations, however, none of this is attempted here, in contrast to previous work by others. In conclusion, I suggest rejection, because of the lack of novelty, and because the models are unrealistic.

As the reviewer points out, many other predictions and observations could be (and have been) used by us and others to assess the realism of mantle convection simulations. However, we politely point out that our goal with this paper is not to identify realistic simulations. Instead, our focus is on assessing the impact of modelling assumptions on predicted surface deflections.

We agree that estimated surface deflections are of higher amplitude than independent estimates of dynamic topography from, for instance, oceanic age-depth residuals, here and in the manuscript. However, as we explain in the manuscript, our ‘approach purposefully avoids isolating passive or plate-driven surface deflections and sub-plate support from numeric simulations or analytic solutions. Since the central focus of this work is merely on quantifying contrasts that arise from choices made when simulating mantle convection, we wish to avoid incorporating additional modifications [e.g. conversion of surface deflections into dynamic topography] where possible’. In other words, we seek to directly assess model output by comparing surface deflections. We note that they are computed directly from the analytic propagator matrix solutions, and by simple (linear) conversion of surface stresses calculated numerically. We think that generating estimates of dynamic topography requires an unavoidable post-hoc conversion of the calculated deflections. Doing so requires, for instance, assessment of lithospheric contributions to surface topography, which are especially poorly understood in the continents. We think that it is instead better for our purposes to focus on assessing impact of modelling choices using the directly computed surface deflections and would like to retain that approach.

With regards to the reviewer’s other comments. First, we think that it is important that this paper is placed within the context of existing observational and theoretical work and would like to retain that material in the introductory parts of the paper. Secondly, we suggest that the novelty of this paper resides in its systematic testing of predicted surface deflections from numeric and analytic approaches to solving the equations of motion of mantle convection. In that sense, whether the models are agreed as being realistic or not is not strictly necessary for this contribution to be a useful one. Nonetheless, we suggest that the model parameterisations embrace a range of purportedly Earth-like values (e.g. for plate motions, viscosity and gravitation) and that the assessment of the impact of different modelling choices would be minimally changed by using other, perhaps more physically accurate, inputs. We will attempt to clarify the rationale for this paper and its goals in the revised manuscript.

Response to “further comments”



line 9: What do you mean with "can vary by ~10%" - do you mean the difference between the numerical (TERRA) and analytical (propagator matrix) approach can be up to 10%?

We will clarify in the revised text that we mean the mean difference in maximum amplitude of predicted surface deflections. 



lines 90 and 270: What do you mean with "radial gravitation" - that the gravity value depends on radius? Because it can also mean that gravity is in radial direction, but that should be an obvious ingredient in any mantle convection model.

We mean that gravitational acceleration can vary as a function of radius (rather than assuming a constant 10 m/s²), which we will clarify in the revised manuscript.

lines 97/98: You switch here from small l to big L, but it appears to be the same.

L = maximum degree, l = degree. What is currently in the manuscript is correct. 

line 108/109: "... amplitudes of calculated deflections will of course not reflect estimated amplitudes of dynamic topography" - why not? I think, quite obviously, amplitudes of calculated deflections should reflect estimated amplitudes of dynamic topography. If they don't and you don't even aim towards it, it shows that your model has a problem and is not realistic.

Dynamic topography need not necessarily equal calculated surface deflections, as we compute them. For instance we have not removed lithospheric isostatic contributions to surface deflections.

Fig. 2c: The range of viscosities plotted is not enough to include all the dashed viscosity profiles. It should go down to 10**20 Pas like in other figures.

We will revise the figure as suggested.  



Fig. 2g: The (continuous) color bar does not correspond to the discontinous coloring in the map. I think the discontinous one is better, because it makes it somewhat easier to distinguish magnitudes. Actually, I don't really like any of these, supposedly scientific, color bars with two colors only, because they make it hard to see anything except the difference between positive and negative.

Thank you for pointing out the rendering issue in the scale bar for Figure 2g. We will correct it in the revised manuscript. With regards to the general comment about colouring, the colour bars we use are (discretely) graduated between, for instance red and blue, or brown and blue. In that sense, more than two colours are used in each scale bar. We think that drawing the eye to positive and negative amplitudes is quite helpful and would like to retain the colour palettes used. Nonetheless, we note that we include histograms and various other statistical information, which we think helps to contextualise the grids. We also note that we have made the grid files available for the interested reader to plot and interrogate as they wish, which we hope is satisfactory.

Eq. (3). I don't understand why you use this admittance formalism. I think it has been thoroughly shown before (including by Colli et al., 2016), that, because of decorrelation of gravity and topography due to different depth-dependent sensitivity kernels, this formalism doesn't lead to meaningful results for density anomalies in the convecting Earth's mantle. Why don't you refer to directly obtained residual topography estimates?

We will add a statement about decorrelation and sensitivity of admittance to this part of the revised manuscript.  We note that we already refer to Colli et al. (2016) here. As we state in the manuscript, we regard Equation 3 as a rule of thumb that Colli et al., (2016)’s results suggest is valid for upper mantle contributions across most spherical harmonic degrees of interest. It is used to generate a curve for visualisation purposes, giving a simple representation of residual topography (see e.g. Holdt et al., 2022). We would therefore like to retain it, rather than showing residual topography estimates, which will make the already quite busy figures even busier.



Fig. 4 and section 5.1.2: Self-gravitation, as the term is usually used, is the gravitation resulting from the non-spherical equipotential surfaces, and this is not the same as using a self-consistent radial gravity profile, which is what you seem to mean here. Hence also the comparison with Ricard (2015) who seems to be concerned with self-gravitation in the usual sense of the word, seems irrelevant.

We agree that our use of the term ‘self-gravitation’ was unclear and needs adjustment. To be clear, with regards to the tests discussed in Section 5.1.2 and shown in Figure 4, we have not modified the governing equations being solved in the propagator matrix code to, for instance, remove or include additional terms relating to self-gravitation. Instead, we (continue to) solve the equations that account for self-gravitation (e.g. the terms that relate to density and gravitational potential in the momentum equation remain). In fact, what we test is the impact of assuming a constant or ‘self-consistent’ (PREM-like) gravity profile and test the impact of including (or excluding) the geoid perturbations on calculated surface deflections. We will endeavour to clarify our approach in the revised manuscript, where we will, for instance, instead refer to ‘self-consistent radial gravity profiles’ as the reviewer suggests. We suggest that the argument we present concerning the differences between models with and without these features is reasonable because these tests incorporate assessment of the impact of self-gravitation. Nonetheless, for clarity and simplicity, we have removed the reference to Ricard (2015) from this part of the manuscript.

Also, I don't think it is very meaningful to compute changes resulting from the gravity of the deflected surface, without including the effect of the internal density anomalies that cause those surface deflections. For one thing, as I already said, the surface deflections are much too large, and the effect of surface deflection and internal density anomalies on gravity largely compensate each other. So, if all these are taken into account, the effect would be even much smaller.

We agree that the impact of density anomalies and boundary deflections on the geoid almost cancel out. We would like to emphasise here and in the revised manuscript that we do not compute deflections that account for gravity perturbations of the deflected surface in isolation. Instead, we either remove or include the impact of geoid perturbations as a whole, i.e., density anomaly and boundary deflection contributions to geoid, since both are wrapped into the "u₃” term (see Equations 5 & 6 of the main manuscript). In our view, it is useful to establish impact on surface deflections from incorporation of gravity from deflected surfaces even when g = 10 m/s² everywhere, which is often assumed in mantle convection simulations and would like to retain these tests in the manuscript.

Line 544 and following: Here, the figure number is missing.

We refer to the figure number on the preceding lines but will add them here too in the revised manuscript for clarity. 

Fig. S1. I don't think it is in any way useful to plot the surface density. It would be more useful to show them at more depths in the mantle. Regarding the densities at depth 270 km, I am wondering why there seem to be positive density anomalies both beneath ridges and beneath subduction zones.

Unfortunately, we disagree with the reviewer’s first assertion. We think that surface densities should be shown because they feature in conversion of surface stress to surface deflections. With regards to the reviewer’s second point about relatively high densities localised beneath some spreading centres (e.g. in the Pacific at 0 Ma at 270 km depth), we suggest that that is a consequence of the interaction of subducted material and mantle flow in the simulations.